| Geant4 Cross Reference |
1
1
2 #include "TFile.h"
2 #include "TFile.h"
3 #include "TTree.h"
3 #include "TTree.h"
4 #include "TH1D.h"
4 #include "TH1D.h"
5 #include "TH2D.h"
5 #include "TH2D.h"
6 #include "TH1I.h"
6 #include "TH1I.h"
7 #include "TF1.h"
7 #include "TF1.h"
8 #include "TProfile2D.h"
8 #include "TProfile2D.h"
9 #include "TGraph.h"
9 #include "TGraph.h"
10 #include "TGraph2D.h"
10 #include "TGraph2D.h"
11 #include "TGraphErrors.h"
11 #include "TGraphErrors.h"
12 #include "TGaxis.h"
12 #include "TGaxis.h"
13 #include "TAxis.h"
13 #include "TAxis.h"
14 #include "TMath.h"
14 #include "TMath.h"
15 #include "TRandom3.h"
15 #include "TRandom3.h"
16 #include "TCanvas.h"
16 #include "TCanvas.h"
17 #include "TLegend.h"
17 #include "TLegend.h"
18 #include <iostream>
18 #include <iostream>
19 #include <fstream>
19 #include <fstream>
20 #include <vector>
20 #include <vector>
21 #include "string.h"
21 #include "string.h"
22 #include <sstream>
22 #include <sstream>
23
23
24 #define pi 3.1415926535
24 #define pi 3.1415926535
25
25
26 using namespace std;
26 using namespace std;
27
27
28 void scattAnalysis()
28 void scattAnalysis()
29 {
29 {
30 gROOT->Reset();
30 gROOT->Reset();
31
31
32 //----------------------------------input--- 32 //----------------------------------input--------------------------------
33 //open a simulation result file
33 //open a simulation result file
34 TFile* file = TFile::Open("output.root");
34 TFile* file = TFile::Open("output.root");
35
35
36 //scattering histogram parameters
36 //scattering histogram parameters
37 Double_t thetalimit = 0.5; //deg
37 Double_t thetalimit = 0.5; //deg
38 Double_t Nbin = 120.;
38 Double_t Nbin = 120.;
39 Double_t thetaMin = 0.; //deg
39 Double_t thetaMin = 0.; //deg
40 Double_t thetaMax = 30.; //deg
40 Double_t thetaMax = 30.; //deg
41
41
42 //options
42 //options
43 Bool_t IwantTotalScatterig = false;
43 Bool_t IwantTotalScatterig = false;
44 Bool_t IWantPlotComptonScattering = false; 44 Bool_t IWantPlotComptonScattering = false;
45 Bool_t IWantPlotPlotComparison = false;
45 Bool_t IWantPlotPlotComparison = false;
46 Bool_t IWantPlotProcessDistrib = false;
46 Bool_t IWantPlotProcessDistrib = false;
47
47
48 //export
48 //export
49 Bool_t IwantExportScattDistrib = false;
49 Bool_t IwantExportScattDistrib = false;
50 Char_t ExportFileName[128];
50 Char_t ExportFileName[128];
51 sprintf(ExportFileName,"SimResults.txt");
51 sprintf(ExportFileName,"SimResults.txt");
52 //---------------------------------------- 52 //-----------------------------------------------------------------------
53
53
54 //definitions
54 //definitions
55 Char_t cutproc[256];
55 Char_t cutproc[256];
56 Double_t val, angle;
56 Double_t val, angle;
57
57
58 //cut definition
58 //cut definition
59 if (IwantTotalScatterig) {
59 if (IwantTotalScatterig) {
60 sprintf(cutproc,"(processID == 1 || proces 60 sprintf(cutproc,"(processID == 1 || processID == 2)");
61 } else {
61 } else {
62 sprintf(cutproc,"processID == 1");
62 sprintf(cutproc,"processID == 1");
63 }
63 }
64
64
65 //define a tree from the ntuple contained in 65 //define a tree from the ntuple contained in the input file
66 TTree* scatt = (TTree*)file->Get("scatt");
66 TTree* scatt = (TTree*)file->Get("scatt");
67 Int_t N = scatt->GetEntries();
67 Int_t N = scatt->GetEntries();
68
68
69 //Rayleigh/Total scattering
69 //Rayleigh/Total scattering
70 TH1D* thetaDistr = new TH1D("thetaDistr", 70 TH1D* thetaDistr = new TH1D("thetaDistr", "", Nbin, thetaMin, thetaMax);
71 scatt->Project("thetaDistr", "theta", cutp 71 scatt->Project("thetaDistr", "theta", cutproc);
72
72
73 Double_t Norig = thetaDistr->Integral();
73 Double_t Norig = thetaDistr->Integral();
74
74
75 for (Int_t i=1; i<=Nbin; i++) {
75 for (Int_t i=1; i<=Nbin; i++) {
76 angle = thetaDistr->GetBinCenter(i);
76 angle = thetaDistr->GetBinCenter(i);
77 if (angle >= thetalimit) {
77 if (angle >= thetalimit) {
78 val = thetaDistr->GetBinContent(i);
78 val = thetaDistr->GetBinContent(i);
79 thetaDistr->SetBinContent(i,val/TMath::S 79 thetaDistr->SetBinContent(i,val/TMath::Sin(angle*pi/180));
80 }
80 }
81 }
81 }
82
82
83 Double_t Ndiv = thetaDistr->Integral();
83 Double_t Ndiv = thetaDistr->Integral();
84 thetaDistr->Scale(Norig/Ndiv);
84 thetaDistr->Scale(Norig/Ndiv);
85
85
86 cout << endl << "Scattering events: " << the 86 cout << endl << "Scattering events: " << thetaDistr->Integral() << endl << endl;
87
87
88 thetaDistr->SetLineColor(kRed);
88 thetaDistr->SetLineColor(kRed);
89 thetaDistr->SetLineWidth(2);
89 thetaDistr->SetLineWidth(2);
90 thetaDistr->SetXTitle("#theta (deg)");
90 thetaDistr->SetXTitle("#theta (deg)");
91 thetaDistr->GetXaxis()->CenterTitle();
91 thetaDistr->GetXaxis()->CenterTitle();
92 thetaDistr->GetXaxis()->SetTitleOffset(1.1 92 thetaDistr->GetXaxis()->SetTitleOffset(1.1);
93 thetaDistr->GetXaxis()->SetRangeUser(0., t 93 thetaDistr->GetXaxis()->SetRangeUser(0., thetaMax);
94 thetaDistr->SetYTitle("Counts (a.u.)");
94 thetaDistr->SetYTitle("Counts (a.u.)");
95 thetaDistr->GetYaxis()->CenterTitle();
95 thetaDistr->GetYaxis()->CenterTitle();
96 thetaDistr->GetYaxis()->SetTitleOffset(1.5 96 thetaDistr->GetYaxis()->SetTitleOffset(1.5);
97 thetaDistr->SetTitle("");
97 thetaDistr->SetTitle("");
98 TCanvas* c1 = new TCanvas("c1","",1200,800 98 TCanvas* c1 = new TCanvas("c1","",1200,800);
99 c1->cd();
99 c1->cd();
100 gStyle->SetOptStat(kFALSE); 100 gStyle->SetOptStat(kFALSE);
101 thetaDistr->Draw("HIST");
101 thetaDistr->Draw("HIST");
102
102
103
103
104 //Compton scattering
104 //Compton scattering
105 if (IWantPlotComptonScattering) {
105 if (IWantPlotComptonScattering) {
106 TH1D* thetaDistrCompt = new TH1D("thetaD 106 TH1D* thetaDistrCompt = new TH1D("thetaDistrCompt", "", Nbin, thetaMin, thetaMax);
107 scatt->Project("thetaDistrCompt", "theta 107 scatt->Project("thetaDistrCompt", "theta", "processID == 2");
108
108
109 Double_t NorigC = thetaDistrCompt->Integ 109 Double_t NorigC = thetaDistrCompt->Integral();
110
110
111 for (Int_t i=1; i<=Nbin; i++) {
111 for (Int_t i=1; i<=Nbin; i++) {
112 angle = thetaDistrCompt->GetBinCenter( 112 angle = thetaDistrCompt->GetBinCenter(i);
113 if (angle >= thetalimit) {
113 if (angle >= thetalimit) {
114 val = thetaDistrCompt->GetBinContent(i 114 val = thetaDistrCompt->GetBinContent(i);
115 thetaDistrCompt->SetBinContent(i,val/T 115 thetaDistrCompt->SetBinContent(i,val/TMath::Sin(angle*pi/180));
116 }
116 }
117 }
117 }
118
118
119 Double_t NdivC = thetaDistrCompt->Integr 119 Double_t NdivC = thetaDistrCompt->Integral();
120 thetaDistrCompt->Scale(NorigC/NdivC);
120 thetaDistrCompt->Scale(NorigC/NdivC);
121
121
122 thetaDistrCompt->SetLineColor(kBlue);
122 thetaDistrCompt->SetLineColor(kBlue);
123 thetaDistrCompt->SetLineWidth(2);
123 thetaDistrCompt->SetLineWidth(2);
124
124
125 if (IWantPlotPlotComparison) {
125 if (IWantPlotPlotComparison) {
126 thetaDistrCompt->Draw("HIST SAME");
126 thetaDistrCompt->Draw("HIST SAME");
127 TLegend* leg = new TLegend(0.7,0.75,0. 127 TLegend* leg = new TLegend(0.7,0.75,0.85,0.85);
128 leg->SetBorderSize(0);
128 leg->SetBorderSize(0);
129 leg->AddEntry(thetaDistr,"Rayleigh","lp" 129 leg->AddEntry(thetaDistr,"Rayleigh","lp");
130 leg->AddEntry(thetaDistrCompt,"Compton 130 leg->AddEntry(thetaDistrCompt,"Compton","lp");
131 leg->Draw();
131 leg->Draw();
132 } else {
132 } else {
133 thetaDistrCompt->SetXTitle("#theta (deg) 133 thetaDistrCompt->SetXTitle("#theta (deg)");
134 thetaDistrCompt->GetXaxis()->CenterTit 134 thetaDistrCompt->GetXaxis()->CenterTitle();
135 thetaDistrCompt->GetXaxis()->SetTitleO 135 thetaDistrCompt->GetXaxis()->SetTitleOffset(1.1);
136 thetaDistrCompt->GetXaxis()->SetRangeU 136 thetaDistrCompt->GetXaxis()->SetRangeUser(0., thetaMax);
137 thetaDistr->SetYTitle("Counts (a.u.)") 137 thetaDistr->SetYTitle("Counts (a.u.)");
138 thetaDistrCompt->GetYaxis()->CenterTit 138 thetaDistrCompt->GetYaxis()->CenterTitle();
139 thetaDistrCompt->GetYaxis()->SetTitleO 139 thetaDistrCompt->GetYaxis()->SetTitleOffset(1.5);
140 thetaDistrCompt->SetTitle("");
140 thetaDistrCompt->SetTitle("");
141 TCanvas* c2 = new TCanvas("c2","",1200 141 TCanvas* c2 = new TCanvas("c2","",1200,800);
142 c2->cd();
142 c2->cd();
143 thetaDistrCompt->Draw("HIST SAME");
143 thetaDistrCompt->Draw("HIST SAME");
144 }
144 }
145
145
146 }
146 }
147
147
148
148
149 //process distribution (0->transport, 1->Ray 149 //process distribution (0->transport, 1->Rayleigh, 2->Compton, 3->Photoel,
150 // 4->Pair prod, 5->Nuc 150 // 4->Pair prod, 5->Nuclear, 6->Diffraction)
151 if (IWantPlotProcessDistrib) {
151 if (IWantPlotProcessDistrib) {
152 TH1D* procDistr = new TH1D("procDistr", "p 152 TH1D* procDistr = new TH1D("procDistr", "process distribution", 70, 0, 7);
153 scatt->Project("procDistr", "processID") 153 scatt->Project("procDistr", "processID");
154 procDistr->SetLineColor(kBlack);
154 procDistr->SetLineColor(kBlack);
155 procDistr->SetLineWidth(2);
155 procDistr->SetLineWidth(2);
156 procDistr->SetXTitle("process");
156 procDistr->SetXTitle("process");
157 procDistr->GetXaxis()->CenterTitle();
157 procDistr->GetXaxis()->CenterTitle();
158 procDistr->GetXaxis()->SetTitleOffset(1. 158 procDistr->GetXaxis()->SetTitleOffset(1.1);
159 procDistr->SetYTitle("Counts");
159 procDistr->SetYTitle("Counts");
160 procDistr->GetYaxis()->CenterTitle();
160 procDistr->GetYaxis()->CenterTitle();
161 procDistr->GetYaxis()->SetTitleOffset(1. 161 procDistr->GetYaxis()->SetTitleOffset(1.2);
162 TCanvas* c3 = new TCanvas("c3","",1200,800 162 TCanvas* c3 = new TCanvas("c3","",1200,800);
163 c3->cd();
163 c3->cd();
164 procDistr->Draw();
164 procDistr->Draw();
165 }
165 }
166
166
167
167
168 //export scattering data
168 //export scattering data
169 if (IwantExportScattDistrib) {
169 if (IwantExportScattDistrib) {
170 ofstream f(ExportFileName);
170 ofstream f(ExportFileName);
171 if (!f) {
171 if (!f) {
172 cout << "Error opening the file!";
172 cout << "Error opening the file!";
173 return;
173 return;
174 }
174 }
175 for (Int_t i=1; i<=Nbin; i++) {
175 for (Int_t i=1; i<=Nbin; i++) {
176 angle = thetaDistr->GetBinCenter(i);
176 angle = thetaDistr->GetBinCenter(i);
177 if (angle >= thetalimit) {
177 if (angle >= thetalimit) {
178 val = thetaDistr->GetBinContent(i);
178 val = thetaDistr->GetBinContent(i);
179 f << angle << " " << val << endl;
179 f << angle << " " << val << endl;
180 }
180 }
181 }
181 }
182 f.close();
182 f.close();
183 cout << "writing successful!" << endl << e 183 cout << "writing successful!" << endl << endl;
184 }
184 }
185
185
186 }
186 }
187
187
188 188