Functions | Variables
PhotonSim_mp.C File Reference
#include <cmath>
#include <iostream>
#include <iomanip>
#include <vector>
#include <unistd.h>
#include <TF1.h>
#include <TH1D.h>
#include <TH2D.h>
#include <TFile.h>
#include <TProfile.h>
#include <TCanvas.h>
#include <TRandom2.h>
#include <TVector3.h>
#include <TString.h>
#include <TLatex.h>
#include <TLine.h>
#include <TROOT.h>
#include <TMath.h>
#include <TThread.h>
#include "include/RootUtils.hh"
#include "include/PhotonTracer.hh"

Go to the source code of this file.

Functions

vector< Double_t > MakeBins (bool isAbs)
 
void * RunOneTrace (void *arg)
 
void PhotonSim_mp ()
 

Variables

const int nThreads = 10
 
int threadIdx [nThreads]
 
Long64_t iterations [nThreads]
 
TH1D * hCollectionZ [nThreads]
 
TH1D * hCollectionT [nThreads]
 
TH2D * hCollectionRate [nThreads]
 
TH2D * hCollectionZ_NBounces [nThreads]
 
TH1D * hFiberAngle [nThreads]
 
TH1D * hNBounces [nThreads]
 

Function Documentation

vector<Double_t> MakeBins ( bool  isAbs)

Definition at line 27 of file PhotonSim_mp.C.

References bins, and e.

Referenced by RunOneTrace().

28 {
29  vector<Double_t> bins;
30 
31  Double_t highEdge, binWidth;
32 
33  if (isAbs) {
34  Double_t currEdge = 0.0;
35  bins.push_back(currEdge);
36  highEdge = 50.0; binWidth = 2.5;
37  while( currEdge < highEdge - 1e-3 ) {currEdge += binWidth; bins.push_back(currEdge);}
38  highEdge = 100.0; binWidth = 5.0;
39  while( currEdge < highEdge - 1e-3 ) {currEdge += binWidth; bins.push_back(currEdge);}
40  highEdge = 200.0; binWidth = 10.0;
41  while( currEdge < highEdge - 1e-3 ) {currEdge += binWidth; bins.push_back(currEdge);}
42  highEdge = 500.0; binWidth = 20.0;
43  while( currEdge < highEdge - 1e-3 ) {currEdge += binWidth; bins.push_back(currEdge);}
44  }
45  else {
46  Double_t currEdge = -500.0;
47  bins.push_back(currEdge);
48  highEdge = -200.0; binWidth = 20.0;
49  while( currEdge < highEdge - 1e-3 ) {currEdge += binWidth; bins.push_back(currEdge);}
50  highEdge = -100.0; binWidth = 10.0;
51  while( currEdge < highEdge - 1e-3 ) {currEdge += binWidth; bins.push_back(currEdge);}
52  highEdge = -50.0; binWidth = 5.0;
53  while( currEdge < highEdge - 1e-3 ) {currEdge += binWidth; bins.push_back(currEdge);}
54  highEdge = 0.0; binWidth = 2.5;
55  while( currEdge < highEdge - 1e-3 ) {currEdge += binWidth; bins.push_back(currEdge);}
56  highEdge = 50.0; binWidth = 2.5;
57  while( currEdge < highEdge - 1e-3 ) {currEdge += binWidth; bins.push_back(currEdge);}
58  highEdge = 100.0; binWidth = 5.0;
59  while( currEdge < highEdge - 1e-3 ) {currEdge += binWidth; bins.push_back(currEdge);}
60  highEdge = 200.0; binWidth = 10.0;
61  while( currEdge < highEdge - 1e-3 ) {currEdge += binWidth; bins.push_back(currEdge);}
62  highEdge = 500.0; binWidth = 20.0;
63  while( currEdge < highEdge - 1e-3 ) {currEdge += binWidth; bins.push_back(currEdge);}
64  }
65 
66  return bins;
67 }
bins
const Binning bins
Definition: nus17_fd_arrays.h:114
e
Float_t e
Definition: plot.C:35
void PhotonSim_mp ( )

Definition at line 202 of file PhotonSim_mp.C.

References demo5::c1, demo5::c2, chisquared::c3, chisquared::c4, plot_validation_datamc::Clone(), om::cout, allTimeWatchdog::endl, file, hCollectionRate, hCollectionT, hCollectionZ, hCollectionZ_NBounces, hFiberAngle, hNBounces, iterations, nThreads, RunOneTrace(), threadIdx, and generate_CCQE_events::threads.

203 {
204  TH1::AddDirectory(false);
205 
206  TThread* threads[nThreads];
207  for (int iThread = 0; iThread < nThreads; ++iThread)
208  {
209  TString threadName = "thread_";
210  threadName += iThread;
211  threadIdx[iThread] = iThread;
212  threads[iThread] = new TThread(threadName, RunOneTrace, (void*)&threadIdx[iThread]);
213  threads[iThread]->Run();
214  }
215  TThread::Ps();
216  for (int iThread = 0; iThread < nThreads; ++iThread)
217  {
218  threads[iThread]->Join();
219  }
220 
221  TH1D* _hCollectionZ = (TH1D*)hCollectionZ[0]->Clone("dZ_CollectionRate");
222  _hCollectionZ->Reset();
223 
224  TH1D* _hCollectionT = (TH1D*)hCollectionT[0]->Clone("dT_CollectionRate");
225  _hCollectionT->Reset();
226 
227  TH2D* _hCollectionRate = (TH2D*)hCollectionRate[0]->Clone("dT_dZ_CollectionRate");
228  _hCollectionRate->Reset();
229 
230  TH2D* _hCollectionZ_NBounces = (TH2D*)hCollectionZ_NBounces[0]->Clone("CollectionZ_NBounces");
231  _hCollectionZ_NBounces->Reset();
232 
233  TH1D* _hNBounces = (TH1D*)hNBounces[0]->Clone("CollectionZ");
234  _hNBounces->Reset();
235 
236  TH1D* _hFiberAngle = (TH1D*)hFiberAngle[0]->Clone("FiberAngle");
237  _hFiberAngle->Reset();
238 
239  Long64_t totIter = 0;
240  for (int iThread = 0; iThread < nThreads; ++iThread)
241  {
242  totIter += iterations[iThread];
243  _hCollectionZ->Add(hCollectionZ[iThread], 1);
244  _hCollectionT->Add(hCollectionT[iThread], 1);
245  _hCollectionRate->Add(hCollectionRate[iThread], 1);
246  _hCollectionZ_NBounces->Add(hCollectionZ_NBounces[iThread], 1);
247  _hNBounces->Add(hNBounces[iThread], 1);
248  _hFiberAngle->Add(hFiberAngle[iThread], 1);
249  cout << "Integral of hCollectionZ[" << iThread << "] = " << hCollectionZ[iThread]->Integral() << endl;
250  }
251  cout << "Integral of _hCollectionZ = " << _hCollectionZ->Integral() << endl;
252  cout << "Total Collection Rate = " << 100.0*(_hCollectionZ->Integral()/totIter) << "%" << endl;
253 
254  _hCollectionZ->Scale(1.0/totIter);
255  _hCollectionT->Scale(1.0/totIter);
256  _hCollectionRate->Scale(1.0/totIter);
257  _hCollectionZ_NBounces->Scale(1.0/totIter);
258 
259  TFile* file = new TFile("dT_dZ_CollectionRate.root","recreate");
260 
261  _hCollectionRate->SetDirectory(file);
262  file->cd();
263  _hCollectionRate->Write();
264 
265  TCanvas* c1 = new TCanvas("CollectionRate", "CollectionRate");
266  c1->cd();
267  _hCollectionRate->Draw("colz");
268 
269  TCanvas* c2 = new TCanvas("CollectionZ", "CollectionZ");
270  c2->cd();
271  _hCollectionZ->Draw();
272 
273 
274  TCanvas* c3 = new TCanvas("CollectionT", "CollectionT");
275  c3->cd();
276  _hCollectionT->Draw();
277 
278  TCanvas* c4 = new TCanvas("CollectionZ_NBounces","CollectionZ_NBounces");
279  c4->cd();
280  _hCollectionZ_NBounces->Draw("colz");
281 
282  TCanvas* c5 = new TCanvas("FiberAngle", "FiberAngle");
283  c5->cd();
284  _hFiberAngle->Draw();
285 
286  TCanvas* c6 = new TCanvas("NBounces","NBounces");
287  c6->cd();
288  _hNBounces->Draw();
289 }
hCollectionRate
TH2D * hCollectionRate[nThreads]
Definition: PhotonSim_mp.C:74
hCollectionZ
TH1D * hCollectionZ[nThreads]
Definition: PhotonSim_mp.C:72
plot_validation_datamc.Clone
def Clone(obj)
Definition: plot_validation_datamc.py:15
RunOneTrace
void * RunOneTrace(void *arg)
Definition: PhotonSim_mp.C:79
demo5.c2
c2
Definition: demo5.py:33
nThreads
const int nThreads
Definition: PhotonSim_mp.C:69
hNBounces
TH1D * hNBounces[nThreads]
Definition: PhotonSim_mp.C:77
om::cout
OStream cout
Definition: OStream.cxx:6
threadIdx
int threadIdx[nThreads]
Definition: PhotonSim_mp.C:70
hCollectionZ_NBounces
TH2D * hCollectionZ_NBounces[nThreads]
Definition: PhotonSim_mp.C:75
hCollectionT
TH1D * hCollectionT[nThreads]
Definition: PhotonSim_mp.C:73
file
TFile * file
Definition: cellShifts.C:17
demo5.c1
c1
Definition: demo5.py:24
hFiberAngle
TH1D * hFiberAngle[nThreads]
Definition: PhotonSim_mp.C:76
generate_CCQE_events.threads
threads
Definition: generate_CCQE_events.py:66
iterations
Long64_t iterations[nThreads]
Definition: PhotonSim_mp.C:71
void* RunOneTrace ( void *  arg)

Definition at line 79 of file PhotonSim_mp.C.

References om::cout, stan::math::fabs(), msf_helper::generator, hCollectionRate, hCollectionT, hCollectionZ, hCollectionZ_NBounces, hFiberAngle, hNBounces, MECModelEnuComparisons::i, iterations, art::left(), MakeBins(), ana::weight, and test::z.

Referenced by PhotonSim_mp().

80 {
81  TThread::Lock();
82 
83  sleep(1);
84  double z = 750.0;
85  TRandom2 generator(0);
86 
87  //double fiberSeparation(6.15);
88  double fiberSeparation(5.0);
89  PhotonTracer tracer(fiberSeparation);
90 
91  int tID = *(int*)arg;
92 
93  //cout << "tID = " << tID << endl;
94 
95  vector<Double_t> Zbins = MakeBins(false);
96  vector<Double_t> ZbinsAbs = MakeBins(true);
97 
98  TString hName = "dZ_CollectionRate";
99  hName += tID;
100  hCollectionZ[tID] = new TH1D(hName, ";#Delta Z;Collection Rate", 500, -500, 500);
101  hCollectionZ[tID]->SetDirectory(0);
102 
103  hName = "dT_CollectionRate";
104  hName += tID;
105  hCollectionT[tID] = new TH1D(hName, ";#Delta T;Collection Rate", 100, 0, 100);
106  hCollectionT[tID]->SetDirectory(0);
107 
108  hName = "dT_dZ_CollectionRate";
109  hName += tID;
110  //hCollectionRate[tID] = new TH2D(hName,";#Delta Z;#Delta T;Collection Rate", ZbinsAbs.size()-1, &ZbinsAbs[0], 30, 0, 60);
111  hCollectionRate[tID] = new TH2D(hName,";#Delta Z;#Delta T;Collection Rate", 500, 0, 500, 30, 0, 60);
112  hCollectionRate[tID]->SetDirectory(0);
113 
114  hName = "CollectionZ_NBounces";
115  hName += tID;
116  hCollectionZ_NBounces[tID] = new TH2D(hName, ";N Bounces;#DeltaZ", 100, 0, 100, 500, 0, 500);
117  hCollectionZ_NBounces[tID]->SetDirectory(0);
118 
119  hName = "NBounces";
120  hName += tID;
121  hNBounces[tID] = new TH1D(hName, ";N Bounces", 100, 0, 100);
122  hNBounces[tID]->SetDirectory(0);
123 
124  hName = "FiberAngle";
125  hName += tID;
126  hFiberAngle[tID] = new TH1D(hName, ";Angle to Fiber", 500, -0.1, TMath::Pi()+0.1);
127  hFiberAngle[tID]->SetDirectory(0);
128 
129  TThread::UnLock();
130 
131  double fracCollected(0.0);
132  double frac100(0.0);
133  double frac150(0.0);
134  double frac200(0.0);
135  double frac250(0.0);
136  double frac300(0.0);
137  double frac350(0.0);
138  double frac400(0.0);
139  double frac450(0.0);
140  double frac500(0.0);
141 
142  double trapEff = 0.054; //From Kuraray documents
143 
144  iterations[tID] = 1e6;
145  //double weight = 1.0/(double)iterations[tID];
146  double weight = 1.0;
147 
148  cout.setf(ios::left);
149  for (Long64_t i = 0; i < iterations[tID]; ++i)
150  {
151  if (i%2000 == 0 && tID == 0)
152  {
153  double percent_done = 100.0*((double)i/(double)iterations[tID]);
154  cout << "percent done: " << setprecision(2) << fixed << percent_done << "%\r" << flush;
155  }
156  //if ( tracer.TraceOnePhoton( z, -999, -999, false ) )
157  if ( tracer.TraceOnePhoton( z, 3.6/2, 5.64/2, false ) )
158  {
159  //tracer.GetPhoton().Print();
160  double phoZ = z - tracer.GetPhoton().GetR0().Z();
161  double phoT = tracer.GetPhoton().GetTime();
162  //double phoW = tracer.GetPhoton().GetWeight();
163  //smear phoT by emission time
164  //phoT += generator.Exp(9.0);
165  hCollectionZ[tID]->Fill( phoZ, weight );
166  hCollectionT[tID]->Fill( phoT, weight );
167  hCollectionRate[tID]->Fill( fabs(phoZ), phoT, weight);
168  hCollectionZ_NBounces[tID]->Fill( tracer.GetPhoton().GetNReflections(), fabs(phoZ), weight );
169  hNBounces[tID]->Fill( tracer.GetPhoton().GetNReflections(), weight );
170  hFiberAngle[tID]->Fill( tracer.GetPhoton().GetFiberAngle() );
171  fracCollected += weight;
172  if (fabs(phoZ) > 100) frac100 += weight;
173  if (fabs(phoZ) > 150) frac150 += weight;
174  if (fabs(phoZ) > 200) frac200 += weight;
175  if (fabs(phoZ) > 250) frac250 += weight;
176  if (fabs(phoZ) > 300) frac300 += weight;
177  if (fabs(phoZ) > 350) frac350 += weight;
178  if (fabs(phoZ) > 400) frac400 += weight;
179  if (fabs(phoZ) > 450) frac450 += weight;
180  if (fabs(phoZ) > 500) frac500 += weight;
181  }
182  }
183 
184  /*
185  if (tID == 0)
186  {
187  cout << "Percent collected = " << setprecision(4) << 100.0*fracCollected << "%" << endl;
188  cout << "Percent past 100 = " << setprecision(4) << 100.0*frac100/fracCollected << "%" << endl;
189  cout << "Percent past 150 = " << setprecision(4) << 100.0*frac150/fracCollected << "%" << endl;
190  cout << "Percent past 200 = " << setprecision(4) << 100.0*frac200/fracCollected << "%" << endl;
191  cout << "Percent past 250 = " << setprecision(4) << 100.0*frac250/fracCollected << "%" << endl;
192  cout << "Percent past 300 = " << setprecision(4) << 100.0*frac300/fracCollected << "%" << endl;
193  cout << "Percent past 350 = " << setprecision(4) << 100.0*frac350/fracCollected << "%" << endl;
194  cout << "Percent past 400 = " << setprecision(4) << 100.0*frac400/fracCollected << "%" << endl;
195  cout << "Percent past 450 = " << setprecision(4) << 100.0*frac450/fracCollected << "%" << endl;
196  cout << "Percent past 500 = " << setprecision(4) << 100.0*frac500/fracCollected << "%" << endl;
197  }
198  */
199  return 0;
200 }
hCollectionRate
TH2D * hCollectionRate[nThreads]
Definition: PhotonSim_mp.C:74
hCollectionZ
TH1D * hCollectionZ[nThreads]
Definition: PhotonSim_mp.C:72
stan::math::fabs
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
msf_helper.generator
generator
Definition: msf_helper.py:26
ana::weight
const Var weight
Definition: PIDCutTuningDefs.h:96
MakeBins
vector< Double_t > MakeBins(bool isAbs)
Definition: PhotonSim_mp.C:27
test.z
z
Definition: test.py:28
hNBounces
TH1D * hNBounces[nThreads]
Definition: PhotonSim_mp.C:77
om::cout
OStream cout
Definition: OStream.cxx:6
art::left
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:104
hCollectionZ_NBounces
TH2D * hCollectionZ_NBounces[nThreads]
Definition: PhotonSim_mp.C:75
hCollectionT
TH1D * hCollectionT[nThreads]
Definition: PhotonSim_mp.C:73
hFiberAngle
TH1D * hFiberAngle[nThreads]
Definition: PhotonSim_mp.C:76
iterations
Long64_t iterations[nThreads]
Definition: PhotonSim_mp.C:71

Variable Documentation

TH2D* hCollectionRate[nThreads]

Definition at line 74 of file PhotonSim_mp.C.

Referenced by PhotonSim(), PhotonSim_mp(), and RunOneTrace().

TH1D* hCollectionT[nThreads]

Definition at line 73 of file PhotonSim_mp.C.

Referenced by PhotonSim(), PhotonSim_mp(), and RunOneTrace().

TH1D* hCollectionZ[nThreads]

Definition at line 72 of file PhotonSim_mp.C.

Referenced by PhotonSim(), PhotonSim_mp(), and RunOneTrace().

TH2D* hCollectionZ_NBounces[nThreads]

Definition at line 75 of file PhotonSim_mp.C.

Referenced by PhotonSim_mp(), and RunOneTrace().

TH1D* hFiberAngle[nThreads]

Definition at line 76 of file PhotonSim_mp.C.

Referenced by PhotonSim_mp(), and RunOneTrace().

TH1D* hNBounces[nThreads]

Definition at line 77 of file PhotonSim_mp.C.

Referenced by PhotonSim_mp(), and RunOneTrace().

Long64_t iterations[nThreads]

Definition at line 71 of file PhotonSim_mp.C.

Referenced by PhotonSim(), PhotonSim_mp(), RunOneTrace(), TestSpacing(), TransverseEff(), and ana::nueccinc::NueCCIncCrossSectionAnalysis::~NueCCIncCrossSectionAnalysis().

const int nThreads = 10

Definition at line 69 of file PhotonSim_mp.C.

Referenced by PhotonSim_mp(), and lem::dec::Forest::Train().

int threadIdx[nThreads]

Definition at line 70 of file PhotonSim_mp.C.

Referenced by PhotonSim_mp().

Generated on Thu Jan 21 2021 14:16:58 for NOvA by   doxygen 1.8.11