| Geant4 Cross Reference |
1
2 #include "TFile.h"
3 #include "TTree.h"
4 #include "TH1D.h"
5 #include "TH2D.h"
6 #include "TH1I.h"
7 #include "TF1.h"
8 #include "TProfile2D.h"
9 #include "TGraph.h"
10 #include "TGraph2D.h"
11 #include "TGraphErrors.h"
12 #include "TGaxis.h"
13 #include "TAxis.h"
14 #include "TMath.h"
15 #include "TRandom3.h"
16 #include "TCanvas.h"
17 #include "TLegend.h"
18 #include <iostream>
19 #include <fstream>
20 #include <vector>
21 #include "string.h"
22 #include <sstream>
23
24 #define pi 3.1415926535
25
26 using namespace std;
27
28 void ADXRD()
29 {
30 gROOT->Reset();
31
32 //----------------------------------input--------------------------------
33 //open a simulation result file
34 TFile* file = TFile::Open("output.root");
35
36 //histogram parameters
37 Double_t binXsize = 0.5; //mm
38 Double_t binXStart = -200.; //mm
39 Double_t binXEnd = -binXStart;
40 Int_t nbinsX = binXEnd*2./binXsize;
41
42 Double_t binYsize = 0.5; //mm
43 Double_t binYStart = -200.; //mm
44 Double_t binYEnd = binXEnd;
45 Int_t nbinsY = binYEnd*2./binYsize;
46
47 Double_t Nbin = 92;
48 Double_t thetaMin = 1; //deg
49 Double_t thetaMax = 27; //deg
50
51 //cuts
52 Double_t Emin = 0.; //keV
53 Double_t Emax = 140.; //keV
54 Double_t angCut = 0.; //deg
55 Bool_t IWantOnlyScatt = true;
56 Bool_t IWantOnlyRayleighScatt = false;
57 Bool_t IWantOnlyComptonScatt = false;
58 Bool_t IWantOnlyDiffraction = false;
59
60 //plot options
61 Bool_t IWantThetaDistribution = true;
62 Bool_t IWantEdistrib = true;
63
64 Bool_t IWantPlot2D = true;
65 Bool_t IWantProfiles = false;
66
67 Bool_t IWantBoxAnalysis = false;
68 Bool_t IWantPlotEtheta = true;
69
70 Int_t nbp = 2;
71 gStyle->SetOptStat(kFALSE);
72 //gStyle->SetPalette(52); //52->gray, 53->hot
73
74 //Energy Distribution Analysis
75 Double_t NbinE = 140;
76 Double_t EMin = 0; //keV
77 Double_t EMax = 140; //keV
78 Double_t Angle = 2.58; //deg
79 Double_t DeltaAngle = 2.; //deg
80 Bool_t ApplyThetaCut = false;
81
82 //export
83 Bool_t IwantExportScattering = false;
84 Char_t ExportScattFileName[128];
85 sprintf(ExportScattFileName,"scatt.txt");
86
87 Bool_t IwantExportVariables = false;
88 Char_t ExportVarFileName[128];
89 sprintf(ExportVarFileName,"var.txt");
90
91 Bool_t IwantExportImage = false;
92 Char_t ExportImageFileName[128];
93 sprintf(ExportImageFileName,"image.txt");
94 //-----------------------------------------------------------------------
95
96 //tree definition
97 Double_t x,y,energy;
98 Int_t kind,ID,nri,nci,ndi;
99 TTree* t1 = (TTree*)file->Get("part");
100 t1->SetBranchAddress("e",&energy);
101 t1->SetBranchAddress("posx",&x);
102 t1->SetBranchAddress("posy",&y);
103 t1->SetBranchAddress("type",&kind);
104 t1->SetBranchAddress("trackID",&ID);
105 t1->SetBranchAddress("NRi",&nri);
106 t1->SetBranchAddress("NCi",&nci);
107 t1->SetBranchAddress("NDi",&ndi);
108 Int_t N = t1->GetEntries();
109 const Int_t Nval = N;
110
111 //scatt cut
112 Char_t cutScatt[128];
113 if (IWantOnlyScatt) {
114 if (IWantOnlyRayleighScatt) {
115 sprintf(cutScatt,"type==0 & trackID==1 & e>=%f & e<=%f & NRi>=1 & NCi==0 & NDi==0",Emin,Emax);
116 } else if (IWantOnlyComptonScatt) {
117 sprintf(cutScatt,"type==0 & trackID==1 & e>=%f & e<=%f & NRi==0 & NCi>=1 & NDi==0",Emin,Emax);
118 } else if (IWantOnlyDiffraction) {
119 sprintf(cutScatt,"type==0 & trackID==1 & e>=%f & e<=%f & NRi==0 & NCi==0 & NDi>1",Emin,Emax);
120 } else {
121 sprintf(cutScatt,"type==0 & trackID==1 & e>=%f & e<=%f & (NRi+NCi+NDi)>=1",Emin,Emax);
122 }
123 } else {
124 sprintf(cutScatt,"type==0 & trackID==1 & e>=%f & e<=%f",Emin,Emax);
125 }
126
127
128 //spatial distribution
129 TH2D* SpatialDistribution = new TH2D("SpatialDistribution", "Spatial Distribution", nbinsX, binXStart, binXEnd, nbinsY, binYStart, binYEnd);
130 t1->Project("SpatialDistribution","posy:posx",cutScatt);
131 if (IWantPlot2D) {
132 SpatialDistribution->SetXTitle("X (mm)");
133 SpatialDistribution->GetXaxis()->CenterTitle();
134 SpatialDistribution->GetXaxis()->SetTitleOffset(1.2);
135 SpatialDistribution->SetYTitle("Y (mm)");
136 SpatialDistribution->GetYaxis()->CenterTitle();
137 SpatialDistribution->GetYaxis()->SetTitleOffset(1.3);
138 TCanvas* c1 = new TCanvas("c1","",1000,1100);
139 c1->cd();
140 SpatialDistribution->SetContour(50);
141 SpatialDistribution->Draw("colz");
142 }
143
144 //profiles (projections transformed in profiles by myself)
145 if (IWantProfiles) {
146 TCanvas* c2 = new TCanvas("c2","",1200,900);
147 c2->Divide(2,1);
148
149 Int_t bcx = Int_t(nbinsX/2.);
150 Int_t bcy = Int_t(nbinsY/2.);
151 Double_t sf = 1./(2.*nbp+1.);
152
153 TH1D* profileX = SpatialDistribution->ProjectionX("px", bcx-nbp, bcx+nbp);
154 profileX->Scale(sf);
155 profileX->SetTitle("X profile");
156 c2->cd(1);
157 //profileX->Fit("gaus");
158 profileX->Draw();
159
160 TH1D* profileY = SpatialDistribution->ProjectionY("py", bcy-nbp, bcy+nbp);
161 profileY->Scale(sf);
162 profileY->SetTitle("Y profile");
163 c2->cd(2);
164 //profileY->Fit("gaus");
165 profileY->Draw();
166
167 cout << endl;
168 cout << "profileX FWHM: " << profileX->GetRMS()*2.35 << " mm" << endl;
169 cout << "profileY FWHM: " << profileY->GetRMS()*2.35 << " mm" << endl;
170 cout << endl;
171 }
172
173
174 //theta distribution
175 TH1D* thetaDistr = new TH1D("thetaDistr", "", Nbin, thetaMin, thetaMax);
176 t1->Project("thetaDistr","acos(momz)*180/acos(-1)",cutScatt);
177
178 Double_t Norig = thetaDistr->Integral();
179
180 Double_t val, angle;
181 for (Int_t i=1; i<=Nbin; i++) {
182 angle = thetaDistr->GetBinCenter(i)*acos(-1)/180;
183 val = thetaDistr->GetBinContent(i);
184 thetaDistr->SetBinContent(i,val/TMath::Sin(angle));
185 }
186
187 Double_t Ndiv = thetaDistr->Integral();
188 thetaDistr->Scale(Norig/Ndiv);
189
190 cout << endl << "total counts of thetaDistr: " << thetaDistr->Integral() << endl;
191
192 if (IWantThetaDistribution) {
193 thetaDistr->SetLineColor(kBlack);
194 thetaDistr->SetLineWidth(2);
195 thetaDistr->SetXTitle("#theta (degree)");
196 thetaDistr->GetXaxis()->CenterTitle();
197 thetaDistr->GetXaxis()->SetTitleOffset(1.1);
198 //thetaDistr->GetXaxis()->SetRangeUser(2., 12.);
199 thetaDistr->SetYTitle("Counts");
200 thetaDistr->GetYaxis()->CenterTitle();
201 thetaDistr->GetYaxis()->SetTitleOffset(1.2);
202 TCanvas* c3 = new TCanvas("c3","",1100,1100);
203 c3->cd();
204 thetaDistr->Draw("HIST");
205 }
206
207
208 //Energy Distribution Analysis
209 if (IWantEdistrib) {
210 //Theta cut
211 Char_t cutTheta[256];
212
213 if (ApplyThetaCut) {
214 sprintf(cutTheta,"type==0 & TMath::Abs(acos(momz)-%f) <= %f",Angle*0.017453293,DeltaAngle*0.017453293);
215 } else {
216 sprintf(cutTheta,"");
217 }
218
219 //Energy distribution of scattered photons
220 TH1D* edistrib = new TH1D("edistrib", "", NbinE, EMin, EMax);
221 t1->Project("edistrib", "e", cutTheta);
222 edistrib->SetLineColor(kRed);
223 edistrib->SetLineWidth(2);
224 edistrib->SetXTitle("E (keV)");
225 edistrib->GetXaxis()->CenterTitle();
226 edistrib->GetXaxis()->SetTitleOffset(1.1);
227 edistrib->SetYTitle("Counts (a. u.)");
228 edistrib->GetYaxis()->CenterTitle();
229 edistrib->GetYaxis()->SetTitleOffset(1.3);
230 edistrib->SetTitle("Energy spectrum of the particles impinging on the detector");
231 TCanvas* c4 = new TCanvas("c4","",1200,800);
232 c4->cd();
233 gStyle->SetOptStat(kFALSE);
234 edistrib->Draw();
235
236 Int_t Ntot = edistrib->Integral();
237 cout << "total counts of edistrib: " << Ntot << endl;
238 }
239
240
241 //Box Score Analysis
242 if (IWantBoxAnalysis) {
243 Int_t xBins = 40;
244 Int_t yBins = 40;
245 Int_t zBins = 20;
246
247 Double_t xlow = -25; //mm
248 Double_t xup = 25;
249 Double_t ylow = -25;
250 Double_t yup = 25;
251 Double_t zlow = -10;
252 Double_t zup = 10;
253
254 Double_t rngadd = 1.; //Box plot range margin (mm)
255
256 TH3D* XYZ = new TH3D("XYZ","Hot spots", xBins, xlow-rngadd, xup+rngadd, zBins, zlow-rngadd, zup+rngadd, yBins, ylow-rngadd, yup+rngadd);
257 t1->Project("XYZ", "posy:posz:posx",cutScatt);
258 XYZ->SetFillColor(2);
259 XYZ->SetXTitle("x (mm)");
260 XYZ->GetXaxis()->CenterTitle();
261 XYZ->GetXaxis()->SetTitleOffset(1.8);
262 XYZ->SetYTitle("z (mm)");
263 XYZ->GetYaxis()->CenterTitle();
264 XYZ->GetYaxis()->SetTitleOffset(2.4);
265 XYZ->SetZTitle("y (mm)");
266 XYZ->GetZaxis()->CenterTitle();
267 XYZ->GetZaxis()->SetTitleOffset(1.8);
268 gStyle->SetCanvasPreferGL(kTRUE);
269 TCanvas* c5 = new TCanvas("c5","Box Score Analysis",1200,800);
270 c5->cd();
271 XYZ->Draw("glbox");
272 }
273
274
275 //energy-angle correlation of impinging photons
276 if (IWantPlotEtheta) {
277 TH2D* Etheta = new TH2D("Etheta", "", NbinE, EMin, EMax, Nbin, thetaMin, thetaMax);
278 t1->Project("Etheta","acos(momz)*180/acos(-1):e");
279 Etheta->SetXTitle("E (keV)");
280 Etheta->GetXaxis()->CenterTitle();
281 Etheta->GetXaxis()->SetTitleOffset(1.2);
282 Etheta->SetYTitle("#theta (degree)");
283 Etheta->GetYaxis()->CenterTitle();
284 Etheta->GetYaxis()->SetTitleOffset(1.3);
285 TCanvas* c6 = new TCanvas("c6","",1000,1100);
286 c6->cd();
287 Etheta->SetContour(50);
288 Etheta->Draw("colz");
289 }
290
291
292 //export scattering data
293 if (IwantExportScattering) {
294 //open a txt file
295 ofstream f(ExportScattFileName);
296 if(!f) {
297 cout << "Error opening the file!";
298 return;
299 }
300 //variables extraction from tree
301 for (Int_t i=1; i<=Nbin; i++) {
302 Double_t ang = thetaDistr->GetBinCenter(i);
303 Double_t scatt = thetaDistr->GetBinContent(i);
304 if (ang >= angCut) {
305 f << ang << " " << scatt << endl;
306 }
307 }
308 //close the txt file
309 f.close();
310 cout << "writing the file with the scattering data successful!" << endl << endl;
311 }
312
313
314 //export (x,y) variables
315 if (IwantExportVariables) {
316 //open a txt file
317 ofstream f(ExportVarFileName);
318 if(!f) {
319 cout << "Error opening the file!";
320 return;
321 }
322 for (Int_t i=0; i<Nval; i++) {
323 t1->GetEntry(i);
324 if (IWantOnlyScatt) {
325 if (kind==0 && ID==1 && energy>=Emin && energy<=Emax && nri+nci+ndi>=1) {
326 f << x << " " << y << endl;
327 }
328 } else {
329 f << x << " " << y << endl;
330 }
331 }
332 //close the txt file
333 f.close();
334 cout << "writing the file with the (x,y) variables successful!" << endl << endl;
335 }
336
337
338 //export Image
339 if (IwantExportImage) {
340 //open a txt file
341 ofstream fout(ExportImageFileName);
342 if (!fout) {
343 cout << "Error opening the file!";
344 return;
345 }
346 //variables extraction from histogram
347 Double_t counts = 0.;
348 for (Int_t i=1; i<=nbinsY; i++) {
349 for (Int_t j=1; j<=nbinsX; j++) {
350 counts = SpatialDistribution->GetBinContent(j,i);
351 fout << counts << " ";
352 }
353 fout << endl;
354 }
355 //close the txt file
356 fout.close();
357 cout << "writing text image file successful!" << endl << endl;
358 }
359
360 }
361
362