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