| Geant4 Cross Reference |
1 #include <TROOT.h> 1
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 prima
31
32 if(argc == 1)
33 {
34 cout << "No input file selected. Please
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. e
51 TH1D *part_per_pixel = new TH1D ("Part per p
52
53 // Spectra for energies (initial, incident a
54 int bin_number = 500;
55 int max_val = 10000;
56 int min_val = 0;
57 double bin_width = (max_val-min_val)/(double
58 TH1D *initial_energies = new TH1D ("energies
59 TH1D *incident_energies = new TH1D ("energie
60 TH1D *dep_energies = new TH1D ("energies3",
61
62 // Spectra for energies (initial, incident a
63 TH1D *initial_energies_P = new TH1D ("energi
64 TH1D *incident_energies_P = new TH1D ("energ
65 TH1D *dep_energies_P = new TH1D ("energies3_
66
67 // Incident particles
68 TH1D *inc_particles = new TH1D ("energies6",
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_na
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", &parent
101 mytree->SetBranchAddress("pixel_number", &pi
102 mytree->SetBranchAddress("step_energy_dep",
103 mytree->SetBranchAddress("step_number", &ste
104 mytree->SetBranchAddress("init_kinetic_energ
105 mytree->SetBranchAddress("kinetic_energy", &
106 mytree->SetBranchAddress("particle_name", &p
107 mytree->SetBranchAddress("pre_step_name", &p
108 mytree->SetBranchAddress("post_step_name", &
109
110 cout << "Reading " << nentries << " from the
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 fo
119 double in_kine = 0; // initial kinetic e
120 double in_kine_P = 0; // initial kinetic e
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 =
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_
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
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 ==
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
186 total_step_dep += step_energy_dep;
187
188 // Sum the total energy deposit from the
189 if (str == "proton" && parentID == 0) to
190 // hist_det_count->Fill(x, y, step_energ
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_ener
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
209 double radius = 27.0;
210 double pi = 3.1415;
211 double T = no_events/(0.407*4*pi*pi*radius*r
212 double S = 2.1; //cm2
213 dep_energies->Scale(1.0/(T*S*bin_width));
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
226 initial_energies->GetXaxis()->SetTitle("Ener
227 initial_energies->GetYaxis()->SetTitle("Coun
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 ene
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_widt
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.0
248 dep_energies_P->GetXaxis()->SetTitle("Energy
249 dep_energies_P->GetYaxis()->SetTitle("Counts
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 dist
257 string label1 = inc_particles->GetXaxis()->G
258 string label2 = inc_particles->GetXaxis()->G
259 string label3 = inc_particles->GetXaxis()->G
260 string label4 = inc_particles->GetXaxis()->G
261 string label5 = inc_particles->GetXaxis()->G
262 Float_t val1 = inc_particles->GetBinContent(
263 Float_t val2 = inc_particles->GetBinContent(
264 Float_t val3 = inc_particles->GetBinContent(
265 Float_t val4 = inc_particles->GetBinContent(
266 Float_t val5 = inc_particles->GetBinContent(
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 dis
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", "Pa
282 part_per_pixel->GetXaxis()->SetTitle("Pixel
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 de
289 TCanvas *c5 = new TCanvas ("c5", "Avg. energ
290 pixel_dep->Divide(part_per_pixel);
291 pixel_dep->GetXaxis()->SetTitle("Pixel numbe
292 pixel_dep->GetYaxis()->SetTitle("Avg. energy
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 }