| Geant4 Cross Reference |
1 #include "TROOT.h"
2 #include "TTree.h"
3 #include "TCanvas.h"
4 #include "TBranch.h"
5 #include <TFile.h>
6 #include <TH1D.h>
7 #include <iostream>
8 #include <fstream>
9 #include <string>
10 #include <sstream>
11 #include <vector>
12 #include <stdio.h>
13 #include <stdlib.h>
14 #include <iomanip>
15 #include <math.h>
16 #include "TLegend.h"
17 #include "TApplication.h"
18 #include "TString.h"
19 #include <algorithm>
20 #include "TPie.h"
21
22 using namespace std;
23 TApplication myapp("app",NULL,NULL);
24
25 const double pi = 3.14159;
26 double to_deg(double angle)
27 {
28 return angle * 180 / pi;
29 }
30
31 //////////////////////////////////////////////////////
32
33 int main(int argc, char *argv[]){
34 // input = ASCII data filename
35 string ifilename = "";
36
37 if(argc == 1)
38 {
39 cout << "No input file selected. Please use: ./read_tree_spectrum [input_file.root]" << endl;
40 return 0;
41 }
42 else
43 {
44 ifilename = argv[1];
45 }
46 TString ifilenameforroot = ifilename;
47
48 // create the file, the Tree and a few branches
49 TFile *f = new TFile(ifilenameforroot, "READ");
50 TTree *mytree = (TTree*)f->Get("XraySPO"); // name of the ntuple in the rootfile
51
52 // create variables to store column values
53 int eventID = 0;
54 int trackID = 0;
55 double x = 0;
56 double y = 0;
57 double z = 0;
58 double theta = 0;
59 double phi = 0;
60 int parentID = 0;
61 char vol_name[500];
62 char proc_name[500];
63 int num_reflections = 0;
64
65 mytree->SetBranchAddress("eventID", &eventID);
66 mytree->SetBranchAddress("vol_name", &vol_name);
67 mytree->SetBranchAddress("trackID", &trackID);
68 mytree->SetBranchAddress("x", &x);
69 mytree->SetBranchAddress("y", &y);
70 mytree->SetBranchAddress("z", &z);
71 mytree->SetBranchAddress("theta", &theta);
72 mytree->SetBranchAddress("phi", &phi);
73 mytree->SetBranchAddress("proc_name", &proc_name);
74 mytree->SetBranchAddress("parentID", &parentID);
75 mytree->SetBranchAddress("num_reflections", &num_reflections);
76
77 int bin_theta = 50;
78 int bin_phi = 100;
79
80 int min_val_theta = 0;
81 int max_val_theta = 5;
82 int min_val_phi = -100;
83 int max_val_phi = 100;
84
85 double bin_width_theta = (max_val_theta-min_val_theta)/(double)bin_theta;
86 double bin_width_phi = (max_val_phi-min_val_phi)/(double)bin_phi;
87
88 TH1D *angle_theta = new TH1D ("theta", "Efficiency - theta [deg]", bin_theta, min_val_theta, max_val_theta);
89 TH1D *angle_phi = new TH1D ("phi", "Efficiency - phi [deg]", bin_phi, min_val_phi, max_val_phi);
90
91 TH1I *reflections = new TH1I ("reflections", "Number of reflections", 8, -1.5, 6.5);
92
93 // Read all entries and fill the histograms
94 Long64_t nentries = mytree->GetEntries();
95 cout << "Reading " << nentries << " from the TTree." << endl;
96
97 int trackID_old = 0;
98
99 bool passed_entrance = false;
100 bool passed_exit = false;
101 int eventID_old = 0;
102
103 int check = -999;
104
105 int count_enter = 0;
106 int count_exit = 0;
107
108 double theta_exit = 0 ;
109 double phi_exit = 0;
110
111 // std::string new_particle_event = "";
112 for (Long64_t i=0; i<nentries; i++)
113 {
114 mytree->GetEntry(i);
115 if (eventID != eventID_old || eventID == 0)
116 {
117 passed_entrance = false;
118 passed_exit = false;
119
120 if ((string)vol_name == "pDummyEntrance")
121 {
122 passed_entrance = true;
123 passed_exit = false;
124 count_enter += 1;
125 }
126 }
127 else
128 {
129 if ((string)vol_name == "pDummyExit" && passed_entrance == true)
130 {
131 // Plot the angle at the
132 theta_exit = theta;
133 phi_exit = phi;
134 passed_exit = true;
135 }
136
137 if (passed_exit)
138 {
139 if ((string)vol_name == "pDummySphere")
140 {
141 passed_entrance = false;
142 passed_exit = false;
143
144 if (num_reflections < 1) num_reflections = -1;
145 reflections->Fill(num_reflections);
146
147 // Check if the eventID has not already been counted before
148 if (eventID == check) cout << "Evento uguale! " << endl;
149
150 theta_exit = pi-theta_exit;
151 if (phi_exit <= 0) phi_exit = -pi-phi_exit;
152 else phi_exit = pi-phi_exit;
153
154 if (theta_exit >= min_val_theta && theta_exit <= max_val_theta)
155 {
156 if (phi_exit >= min_val_phi && phi_exit <= max_val_phi)
157 {
158 count_exit+=1;
159 angle_theta->Fill(to_deg(theta_exit));
160 angle_phi->Fill(to_deg(phi_exit));
161 }
162 }
163 check = eventID;
164 }
165 }
166 }
167 eventID_old = eventID;
168 }
169
170 // pie chart
171 Float_t val0 = reflections->GetBinContent(1); // number of entries with 0 reflections
172 Float_t val2 = reflections->GetBinContent(3); // number of entries with 1 reflection
173 Float_t val3 = reflections->GetBinContent(4);
174 Float_t val4 = reflections->GetBinContent(5);
175 Float_t val5 = reflections->GetBinContent(6) + reflections->GetBinContent(7); // value for 4+ reflections
176 cout << "Values are: " << val0 << " " << val2 << " " << val3 << " " << val4 << " " << val5 << " " << endl;
177
178
179 Float_t vals[] = {val0,val2,val3,val4,val5};
180 Int_t colors[] = {0,1,2,3,4};
181 Int_t nvals = sizeof(vals)/sizeof(vals[0]);
182
183 TCanvas *cpie = new TCanvas("cpie", "cpie", 0, 0, 900, 900);
184 TPie *pie1 = new TPie("pie1", "100 keV, 1deg., SS model",nvals,vals,colors);
185 pie1->SetLabelsOffset(.01);
186 pie1->SetRadius(.2);
187
188 pie1->SetEntryLabel(0, "0");
189 pie1->SetEntryLabel(1, "1");
190 pie1->SetEntryLabel(2, "2");
191 pie1->SetEntryLabel(3, "3");
192 pie1->SetEntryLabel(4, ">3");
193 pie1->SetLabelFormat("%txt (%perc)");
194 pie1->Draw();
195
196 // Normalization
197 cout << "Count particle: <entered: " << count_enter << "> <exited: " << count_exit << ">" << endl;
198 angle_theta->Scale(1.0/((double)count_enter));
199 angle_phi->Scale(1.0/((double)count_enter));
200 angle_theta->Scale(1.0/(bin_width_theta));
201 angle_phi->Scale(1.0/(bin_width_phi));
202
203 TCanvas *c1 = new TCanvas("c1", "c1", 0, 50, 900, 900);
204 int c1_id = c1->GetCanvasID();
205 cout << "The first canvas ID is: " << c1_id << endl;
206 c1->SetLogy();
207 angle_theta->SetFillColor(kBlack);
208 angle_theta->SetMarkerColor(kBlack);
209 angle_theta->Draw("");
210 angle_theta->GetXaxis()->SetTitle("Theta [deg]");
211 angle_theta->GetYaxis()->SetTitle("Normalized efficiency [deg^-1]");
212 angle_theta->GetXaxis()->SetTitleSize(0.03);
213 angle_theta->GetYaxis()->SetTitleSize(0.03);
214 cout << "integral of angle_theta from histo:" << angle_theta->Integral() << " and from fraction: " << count_exit/count_enter << endl;
215
216 TCanvas *c2 = new TCanvas("c2", "c2", 0, 100, 900, 900);
217 int c2_id = c2->GetCanvasID();
218 cout << "The second canvas ID is: " << c2_id << endl;
219 c2->SetLogy();
220 angle_phi->SetFillColor(kBlack);
221 angle_phi->SetMarkerColor(kBlack);
222 angle_phi->Draw("");
223 angle_phi->GetXaxis()->SetTitle("Phi [deg]");
224 angle_phi->GetYaxis()->SetTitle("Normalized efficiency [deg^-1]");
225 angle_phi->GetXaxis()->SetTitleSize(0.03);
226 angle_phi->GetYaxis()->SetTitleSize(0.03);
227 cout << "integral of angle_phi is:" << angle_theta->Integral() << " and from fraction: " << count_exit/count_enter << endl;
228
229 myapp.Run(true);
230
231 f->Close();
232 return 0;
233 }
234