PhotonSim.C
Go to the documentation of this file.
1 #include <cmath>
2 #include <iostream>
3 #include <iomanip>
4 #include <vector>
5 
6 #include <TF1.h>
7 #include <TH1D.h>
8 #include <TH2D.h>
9 #include <TFile.h>
10 #include <TProfile.h>
11 #include <TCanvas.h>
12 #include <TRandom2.h>
13 #include <TVector3.h>
14 #include <TString.h>
15 #include <TLatex.h>
16 #include <TLine.h>
17 #include <TROOT.h>
18 #include <TMath.h>
19 
20 #include "include/RootUtils.hh"
21 #include "include/PhotonTracer.hh"
22 
23 using namespace std;
24 using namespace TMath;
25 
26 vector<Double_t> MakeBins(bool isAbs)
27 {
28  vector<Double_t> bins;
29 
30  Double_t highEdge, binWidth;
31 
32  if (isAbs) {
33  Double_t currEdge = 0.0;
34  bins.push_back(currEdge);
35  highEdge = 50.0; binWidth = 2.5;
36  while( currEdge < highEdge - 1e-3 ) {currEdge += binWidth; bins.push_back(currEdge);}
37  highEdge = 100.0; binWidth = 5.0;
38  while( currEdge < highEdge - 1e-3 ) {currEdge += binWidth; bins.push_back(currEdge);}
39  highEdge = 200.0; binWidth = 10.0;
40  while( currEdge < highEdge - 1e-3 ) {currEdge += binWidth; bins.push_back(currEdge);}
41  highEdge = 500.0; binWidth = 20.0;
42  while( currEdge < highEdge - 1e-3 ) {currEdge += binWidth; bins.push_back(currEdge);}
43  }
44  else {
45  Double_t currEdge = -500.0;
46  bins.push_back(currEdge);
47  highEdge = -200.0; binWidth = 20.0;
48  while( currEdge < highEdge - 1e-3 ) {currEdge += binWidth; bins.push_back(currEdge);}
49  highEdge = -100.0; binWidth = 10.0;
50  while( currEdge < highEdge - 1e-3 ) {currEdge += binWidth; bins.push_back(currEdge);}
51  highEdge = -50.0; binWidth = 5.0;
52  while( currEdge < highEdge - 1e-3 ) {currEdge += binWidth; bins.push_back(currEdge);}
53  highEdge = 0.0; binWidth = 2.5;
54  while( currEdge < highEdge - 1e-3 ) {currEdge += binWidth; bins.push_back(currEdge);}
55  highEdge = 50.0; binWidth = 2.5;
56  while( currEdge < highEdge - 1e-3 ) {currEdge += binWidth; bins.push_back(currEdge);}
57  highEdge = 100.0; binWidth = 5.0;
58  while( currEdge < highEdge - 1e-3 ) {currEdge += binWidth; bins.push_back(currEdge);}
59  highEdge = 200.0; binWidth = 10.0;
60  while( currEdge < highEdge - 1e-3 ) {currEdge += binWidth; bins.push_back(currEdge);}
61  highEdge = 500.0; binWidth = 20.0;
62  while( currEdge < highEdge - 1e-3 ) {currEdge += binWidth; bins.push_back(currEdge);}
63  }
64 
65  return bins;
66 }
67 
68 
69 void PhotonSim()
70 {
71  double z = 750.0;
72  TRandom2 generator(0);
73 
74  //double fiberSeparation = 6.15;
75  double fiberSeparation = 5.0;
76  PhotonTracer tracer(fiberSeparation);
77 
78  vector<Double_t> Zbins = MakeBins(false);
79  vector<Double_t> ZbinsAbs = MakeBins(true);
80 
81  TFile* file = new TFile("dT_dZ_CollectionRate.root","recreate");
82  TH1D* hCollectionZ = new TH1D("dZ_CollectionRate", ";#Delta Z;Collection Rate", Zbins.size()-1, &Zbins[0]);
83  TH1D* hCollectionT = new TH1D("dT_CollectionRate", ";#Delta T;Collection Rate", 100, 0, 100);
84  TH2D* hCollectionRate = new TH2D("dT_dZ_CollectionRate",";#Delta Z;#Delta T;Collection Rate", ZbinsAbs.size()-1, &ZbinsAbs[0], 20, 0, 100);
85 
86  double fracCollected(0.0);
87  double frac100(0.0);
88  double frac150(0.0);
89  double frac200(0.0);
90  double frac250(0.0);
91  double frac300(0.0);
92  double frac350(0.0);
93  double frac400(0.0);
94  double frac450(0.0);
95  double frac500(0.0);
96 
97  double trapEff = 0.054; //From Kuraray documents
98 
99  Long64_t iterations = 1e4;
100  double weight = 1.0/(double)iterations;
101 
102  cout.setf(ios::left);
103  for (Long64_t i = 0; i < iterations; ++i)
104  {
105  if (i%2000 == 0)
106  {
107  double percent_done = 100.0*((double)i/(double)iterations);
108  cout << "percent done: " << setprecision(2) << fixed << percent_done << "%\r" << flush;
109  }
110  if ( tracer.TraceOnePhoton( z, 1.8, 2.82 ) )
111  {
112  double phoZ = z - tracer.GetPhoton().GetR0().Z();
113  double phoT = tracer.GetPhoton().GetTime();
114  //double phoW = tracer.GetPhoton().GetWeight();
115  //smear phoT by emission time
116  phoT += generator.Exp(9.0);
117  hCollectionZ->Fill( phoZ, weight );
118  hCollectionT->Fill( phoT, weight );
119  hCollectionRate->Fill( fabs(phoZ), phoT, weight);
120 
121  fracCollected += weight;
122  if (fabs(phoZ) > 100) frac100 += weight;
123  if (fabs(phoZ) > 150) frac150 += weight;
124  if (fabs(phoZ) > 200) frac200 += weight;
125  if (fabs(phoZ) > 250) frac250 += weight;
126  if (fabs(phoZ) > 300) frac300 += weight;
127  if (fabs(phoZ) > 350) frac350 += weight;
128  if (fabs(phoZ) > 400) frac400 += weight;
129  if (fabs(phoZ) > 450) frac450 += weight;
130  if (fabs(phoZ) > 500) frac500 += weight;
131  }
132  }
133 
134  cout << "Percent collected = " << setprecision(4) << 100.0*fracCollected << "%" << endl;
135  cout << "Percent past 100 = " << setprecision(4) << 100.0*frac100/fracCollected << "%" << endl;
136  cout << "Percent past 150 = " << setprecision(4) << 100.0*frac150/fracCollected << "%" << endl;
137  cout << "Percent past 200 = " << setprecision(4) << 100.0*frac200/fracCollected << "%" << endl;
138  cout << "Percent past 250 = " << setprecision(4) << 100.0*frac250/fracCollected << "%" << endl;
139  cout << "Percent past 300 = " << setprecision(4) << 100.0*frac300/fracCollected << "%" << endl;
140  cout << "Percent past 350 = " << setprecision(4) << 100.0*frac350/fracCollected << "%" << endl;
141  cout << "Percent past 400 = " << setprecision(4) << 100.0*frac400/fracCollected << "%" << endl;
142  cout << "Percent past 450 = " << setprecision(4) << 100.0*frac450/fracCollected << "%" << endl;
143  cout << "Percent past 500 = " << setprecision(4) << 100.0*frac500/fracCollected << "%" << endl;
144 
145  hCollectionRate->SetDirectory(file);
146  file->cd();
147  hCollectionRate->Write();
148 
149  TCanvas* c1 = new TCanvas("CollectionRate", "CollectionRate");
150  c1->cd();
151  hCollectionRate->Draw("colz");
152 
153  TCanvas* c2 = new TCanvas("CollectionZ", "CollectionZ");
154  c2->cd();
155  hCollectionZ->Draw();
156 
157 
158  TCanvas* c3 = new TCanvas("CollectionT", "CollectionT");
159  c3->cd();
160  hCollectionT->Draw();
161 
162  TCanvas* c4 = new TCanvas("cRefXY", "cRefXY");
163  c4->cd();
164  tracer.GetRefXY()->Draw("colz");
165 
166  TCanvas* c5 = new TCanvas("cAttenXY", "cAttenXY");
167  c5->cd();
168  tracer.GetAttenXY()->Draw("colz");
169 
170  TCanvas* c6 = new TCanvas("cAbsXY", "cAbsXY");
171  c6->cd();
172  tracer.GetAbsXY()->Draw("colz");
173  return;
174 }
175 
TH2D * hCollectionRate[nThreads]
Definition: PhotonSim_mp.C:74
TH1D * hCollectionZ[nThreads]
Definition: PhotonSim_mp.C:72
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
const Var weight
void PhotonSim()
Definition: PhotonSim.C:69
vector< Double_t > MakeBins(bool isAbs)
Definition: PhotonSim.C:26
c2
Definition: demo5.py:33
const Binning bins
z
Definition: test.py:28
OStream cout
Definition: OStream.cxx:6
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:104
TH1D * hCollectionT[nThreads]
Definition: PhotonSim_mp.C:73
TFile * file
Definition: cellShifts.C:17
c1
Definition: demo5.py:24
Float_t e
Definition: plot.C:35
Long64_t iterations[nThreads]
Definition: PhotonSim_mp.C:71