| Geant4 Cross Reference |
1 #include "TROOT.h" 1
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 us
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 branc
49 TFile *f = new TFile(ifilenameforroot, "READ
50 TTree *mytree = (TTree*)f->Get("XraySPO");
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_na
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_
74 mytree->SetBranchAddress("parentID", &parent
75 mytree->SetBranchAddress("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_
86 double bin_width_phi = (max_val_phi-min_val_
87
88 TH1D *angle_theta = new TH1D ("theta", "Effi
89 TH1D *angle_phi = new TH1D ("phi", "Efficien
90
91 TH1I *reflections = new TH1I ("reflections",
92
93 // Read all entries and fill the histograms
94 Long64_t nentries = mytree->GetEntries();
95 cout << "Reading " << nentries << " from the
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" &&
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_reflect
145 reflections->Fill(num_reflections);
146
147 // Check if the eventID has not alre
148 if (eventID == check) cout << "Event
149
150 theta_exit = pi-theta_exit;
151 if (phi_exit <= 0) phi_exit = -pi-ph
152 else phi_exit = pi-phi_exit;
153
154 if (theta_exit >= min_val_theta && t
155 {
156 if (phi_exit >= min_val_phi && phi
157 {
158 count_exit+=1;
159 angle_theta->Fill(to_deg(theta_e
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)
172 Float_t val2 = reflections->GetBinContent(3)
173 Float_t val3 = reflections->GetBinContent(4)
174 Float_t val4 = reflections->GetBinContent(5)
175 Float_t val5 = reflections->GetBinContent(6)
176 cout << "Values are: " << val0 << " " << val
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",
184 TPie *pie1 = new TPie("pie1", "100 keV, 1deg
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: " << coun
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,
204 int c1_id = c1->GetCanvasID();
205 cout << "The first canvas ID is: " << c1_id
206 c1->SetLogy();
207 angle_theta->SetFillColor(kBlack);
208 angle_theta->SetMarkerColor(kBlack);
209 angle_theta->Draw("");
210 angle_theta->GetXaxis()->SetTitle("Theta [de
211 angle_theta->GetYaxis()->SetTitle("Normalize
212 angle_theta->GetXaxis()->SetTitleSize(0.03);
213 angle_theta->GetYaxis()->SetTitleSize(0.03);
214 cout << "integral of angle_theta from histo:
215
216 TCanvas *c2 = new TCanvas("c2", "c2", 0, 100
217 int c2_id = c2->GetCanvasID();
218 cout << "The second canvas ID is: " << c2_id
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
225 angle_phi->GetXaxis()->SetTitleSize(0.03);
226 angle_phi->GetYaxis()->SetTitleSize(0.03);
227 cout << "integral of angle_phi is:" << angle
228
229 myapp.Run(true);
230
231 f->Close();
232 return 0;
233 }
234