| Geant4 Cross Reference |
1 #include <TROOT.h>
2 #include <TTree.h>
3 #include <TCanvas.h>
4 #include <TApplication.h>
5 #include <TString.h>
6 #include <TBranch.h>
7 #include <TStyle.h>
8 #include <TPie.h>
9 #include <TFile.h>
10 #include <TLegend.h>
11 #include <TH1D.h>
12 #include <TH2D.h>
13 #include <iostream>
14 #include <fstream>
15 #include <string>
16 #include <sstream>
17 #include <vector>
18 #include <stdio.h>
19 #include <stdlib.h>
20 #include <iomanip>
21 #include <algorithm>
22
23 using namespace std;
24 TApplication myapp("app",NULL,NULL);
25
26 //////////////////////////////////////////////////////
27
28 int main(int argc, char *argv[]){
29 string ifilename = "";
30 int no_events = 10000000; // number of primaries
31
32 if(argc == 1)
33 {
34 cout << "No input file selected. Please use: ./read_tree_spectrum [input_file.root]" << endl;
35 return 0;
36 }
37 else
38 {
39 ifilename = argv[1];
40 if (argc == 3) no_events = atoi(argv[2]);
41 }
42
43 TString ifilenameforroot = ifilename;
44 TFile *f = new TFile(ifilenameforroot, "READ");
45
46 int no_pixel = 317;
47
48 // Histogram definition
49 // Detector deposits
50 TH1D *pixel_dep = new TH1D ("pixel", "Avg. energy deposit per pixel [MeV]", no_pixel, 0, no_pixel); //(nbinx, xdown, xup, nbiny, ylow, yup) for drawing the hexagon
51 TH1D *part_per_pixel = new TH1D ("Part per pixel", "No. particles per pixel", no_pixel, 0, no_pixel);
52
53 // Spectra for energies (initial, incident and deposited)
54 int bin_number = 500;
55 int max_val = 10000;
56 int min_val = 0;
57 double bin_width = (max_val-min_val)/(double)bin_number;
58 TH1D *initial_energies = new TH1D ("energies4", "GCR Protons impacting spectra on the detector", bin_number, min_val, max_val);
59 TH1D *incident_energies = new TH1D ("energies2", " ", bin_number, min_val, max_val);
60 TH1D *dep_energies = new TH1D ("energies3", " ", bin_number, min_val, max_val);
61
62 // Spectra for energies (initial, incident and deposited) - for primaries (protons)
63 TH1D *initial_energies_P = new TH1D ("energies4_P", " ", bin_number, min_val, max_val);
64 TH1D *incident_energies_P = new TH1D ("energies2_P", " ", bin_number, min_val, max_val);
65 TH1D *dep_energies_P = new TH1D ("energies3_P", "Primary GCR Protons impacting spectra on the detector", bin_number, min_val, max_val);
66
67 // Incident particles
68 TH1D *inc_particles = new TH1D ("energies6", "GCR Protons particle fluxes on the detector composition", 5, 0, 5);
69
70 // Variables to store the tuples' values
71 int eventID = 0;
72 char vol_name[500];
73 int trackID = 0;
74 double x = 0;
75 double y = 0;
76 double z = 0;
77 double theta = 0;
78 double phi = 0;
79 int parentID = 0;
80 int pixel_number = 0;
81 double step_energy_dep = 0;
82 int step_number = 0;
83 double init_kinetic_energy = 0;
84 double kinetic_energy = 0;
85 char particle_name[500];
86 char pre_step_name[500];
87 char post_step_name[500];
88
89 // Define tuple elements
90 TTree *mytree = (TTree*)f->Get("TES_Tuple");
91 Long64_t nentries = mytree->GetEntries();
92 mytree->SetBranchAddress("eventID", &eventID);
93 mytree->SetBranchAddress("vol_name", &vol_name);
94 mytree->SetBranchAddress("trackID", &trackID);
95 mytree->SetBranchAddress("x", &x);
96 mytree->SetBranchAddress("y", &y);
97 mytree->SetBranchAddress("z", &z);
98 mytree->SetBranchAddress("theta", &theta);
99 mytree->SetBranchAddress("phi", &phi);
100 mytree->SetBranchAddress("parentID", &parentID);
101 mytree->SetBranchAddress("pixel_number", &pixel_number);
102 mytree->SetBranchAddress("step_energy_dep", &step_energy_dep);
103 mytree->SetBranchAddress("step_number", &step_number);
104 mytree->SetBranchAddress("init_kinetic_energy", &init_kinetic_energy);
105 mytree->SetBranchAddress("kinetic_energy", &kinetic_energy);
106 mytree->SetBranchAddress("particle_name", &particle_name);
107 mytree->SetBranchAddress("pre_step_name", &pre_step_name);
108 mytree->SetBranchAddress("post_step_name", &post_step_name);
109
110 cout << "Reading " << nentries << " from the TTree." << endl;
111 string vol_name_old = "";
112
113 int eventID_old = 0;
114 int trackID_old = 0;
115 double total_step_dep = 0;
116 double total_step_dep_P = 0;
117 double kine = 0; // kinetic energy
118 double kine_P = 0; // kinetic energy for protons
119 double in_kine = 0; // initial kinetic energy
120 double in_kine_P = 0; // initial kinetic energy for protons
121 bool entered = false;
122 bool entered_P = false;
123 char arr[50];
124 string str = "";
125 int pixel_number_old = 0;
126
127 for (Long64_t i=0; i<nentries; i++)
128 {
129 mytree->GetEntry(i);
130
131 if (entered || entered_P)
132 {
133 if (eventID != eventID_old || (eventID == eventID_old && trackID != trackID_old))
134 {
135 if (entered && total_step_dep > 0)
136 {
137 initial_energies->Fill(in_kine);
138 incident_energies->Fill(kine);
139 dep_energies->Fill(total_step_dep);
140 inc_particles->Fill(arr, 1);
141 total_step_dep = 0;
142 }
143 if (entered_P && total_step_dep_P > 0)
144 {
145 initial_energies_P->Fill(in_kine_P);
146 incident_energies_P->Fill(kine_P);
147 dep_energies_P->Fill(total_step_dep_P);
148 total_step_dep_P = 0;
149 }
150 entered = false;
151 entered_P = false;
152 }
153 }
154 str = (string)particle_name;
155
156 if((string)vol_name == "Bipxl")
157 {
158 if (vol_name_old != "Bipxl")
159 {
160 if (!entered)
161 {
162 entered = true;
163 kine = kinetic_energy;
164 in_kine = init_kinetic_energy;
165
166 if (str != "proton" && str != "gamma" && str != "e+" && str != "e-")
167 {
168 str = "Other particles";
169 }
170 strcpy(arr, str.c_str());
171
172 // Same but onyl for primaries
173 if (!entered_P)
174 {
175 if (str == "proton" && parentID == 0)
176 {
177 entered_P = true;
178 kine_P = kinetic_energy;
179 in_kine_P = init_kinetic_energy;
180 }
181 }
182 }
183 }
184
185 // Sum the total energy deposit from all particles
186 total_step_dep += step_energy_dep;
187
188 // Sum the total energy deposit from the primaries
189 if (str == "proton" && parentID == 0) total_step_dep_P += step_energy_dep;
190 // hist_det_count->Fill(x, y, step_energy_dep);
191
192 if (pixel_number != pixel_number_old)
193 {
194 part_per_pixel->Fill(-pixel_number);
195 }
196
197 pixel_dep->Fill(-pixel_number, step_energy_dep);
198 }
199 eventID_old = eventID;
200 trackID_old = trackID;
201 vol_name_old = (string)vol_name;
202 pixel_number_old = pixel_number;
203 }
204
205 cout << "Starts plotting..." << endl;
206
207 // Draw histograms
208 TCanvas *c1 = new TCanvas ("c1", "Energy dep", 0, 0, 1000, 900);
209 double radius = 27.0;
210 double pi = 3.1415;
211 double T = no_events/(0.407*4*pi*pi*radius*radius); //equivalent time
212 double S = 2.1; //cm2
213 dep_energies->Scale(1.0/(T*S*bin_width)); //scale based on the equivalent time and detector surface
214 initial_energies->Scale(1.0/(T*S*bin_width));
215 incident_energies->Scale(1.0/(T*S*bin_width));
216 gStyle->SetOptStat(0);
217 c1->SetLogx();
218 c1->SetLogy();
219 dep_energies->SetLineColor(kRed);
220 initial_energies->Draw("hist");
221 initial_energies->SetLineColor(kBlack);
222 incident_energies->Draw("histsame");
223 dep_energies->Draw("histsame");
224 incident_energies->SetLineColor(kBlue);
225 initial_energies->GetYaxis()->SetRangeUser(1.0e-6, 1.0);
226 initial_energies->GetXaxis()->SetTitle("Energy [MeV]");
227 initial_energies->GetYaxis()->SetTitle("Counts/cm2/s/MeV");
228 auto legend1 = new TLegend(0.75,0.8,0.9,0.9);
229 legend1->AddEntry(dep_energies,"Edep");
230 legend1->AddEntry(incident_energies, "Einc");
231 legend1->AddEntry(initial_energies, "Ei");
232 legend1->Draw("same");
233
234 TCanvas *c2 = new TCanvas ("c2", "Proton energy dep", 50, 0, 1000, 900);
235 dep_energies_P->Scale(1.0/(T*S*bin_width));
236 initial_energies_P->Scale(1.0/(T*S*bin_width));
237 incident_energies_P->Scale(1.0/(T*S*bin_width));
238 gStyle->SetOptStat(0);
239 c2->SetLogx();
240 c2->SetLogy();
241 dep_energies_P->Draw("hist");
242 dep_energies_P->SetLineColor(kRed);
243 initial_energies_P->Draw("histsame");
244 initial_energies_P->SetLineColor(kBlack);
245 incident_energies_P->Draw("histsame");
246 incident_energies_P->SetLineColor(kBlue);
247 dep_energies_P->GetYaxis()->SetRangeUser(1.0e-6, 1.0);
248 dep_energies_P->GetXaxis()->SetTitle("Energy [MeV]");
249 dep_energies_P->GetYaxis()->SetTitle("Counts/cm2/s/MeV");
250 auto legend2 = new TLegend(0.75,0.8,0.9,0.9);
251 legend2->AddEntry(dep_energies,"Edep");
252 legend2->AddEntry(incident_energies, "Einc");
253 legend2->AddEntry(initial_energies, "Ei");
254 legend2->Draw("same");
255
256 TCanvas *cpie = new TCanvas("Particless distribution", "Particles distribution", 100, 0, 1000, 1000);
257 string label1 = inc_particles->GetXaxis()->GetBinLabel(1);
258 string label2 = inc_particles->GetXaxis()->GetBinLabel(2);
259 string label3 = inc_particles->GetXaxis()->GetBinLabel(3);
260 string label4 = inc_particles->GetXaxis()->GetBinLabel(4);
261 string label5 = inc_particles->GetXaxis()->GetBinLabel(5);
262 Float_t val1 = inc_particles->GetBinContent(1);
263 Float_t val2 = inc_particles->GetBinContent(2);
264 Float_t val3 = inc_particles->GetBinContent(3);
265 Float_t val4 = inc_particles->GetBinContent(4);
266 Float_t val5 = inc_particles->GetBinContent(5);
267 Float_t vals[] = {val1,val2,val3,val4,val5};
268 Int_t colors[] = {1,2,3,4,5};
269 Int_t nvals = sizeof(vals)/sizeof(vals[0]);
270 TPie *pie1 = new TPie("pie1", "Particles distribution",nvals,vals,colors);
271 pie1->SetLabelsOffset(.01);
272 pie1->SetRadius(.2);
273 pie1->SetEntryLabel(0, label1.c_str());
274 pie1->SetEntryLabel(1, label2.c_str());
275 pie1->SetEntryLabel(2, label3.c_str());
276 pie1->SetEntryLabel(3, label4.c_str());
277 pie1->SetEntryLabel(4, label5.c_str());
278 pie1->SetLabelFormat("%txt (%perc)");
279 pie1->Draw();
280
281 TCanvas *c4_num = new TCanvas ("c4_num", "Particle count per pixel", 150, 0, 1000, 900);
282 part_per_pixel->GetXaxis()->SetTitle("Pixel number");
283 part_per_pixel->GetYaxis()->SetTitle("Counts");
284 part_per_pixel->SetFillColor(kBlack);
285 part_per_pixel->SetMarkerColor(kBlack);
286 part_per_pixel->Draw();
287
288 // Draw the histogram of the single pixel deposition
289 TCanvas *c5 = new TCanvas ("c5", "Avg. energy deposit per pixel", 250, 0, 1000, 900);
290 pixel_dep->Divide(part_per_pixel);
291 pixel_dep->GetXaxis()->SetTitle("Pixel number");
292 pixel_dep->GetYaxis()->SetTitle("Avg. energy deposit [MeV]");
293 pixel_dep->SetFillColor(kBlack);
294 pixel_dep->SetMarkerColor(kBlack);
295 pixel_dep->Draw("");
296
297 myapp.Run(true);
298 f->Close();
299 return 0;
300 }