42 {
43 cout << "MarkovPathReader" << endl
44 << "Written by Martijn Jongen" << endl
45 << endl ;
46
48 string ifname_model ;
49 string ofname ;
51
52 try {
54 zap[
"f"] =
make_field(ifnames,
"input file names (binary files containing JPhotonPaths)") ;
55 zap[
"m"] =
make_field(ifname_model,
"input file name (root file containing JScatteringModel)") ;
56 zap[
"o"] =
make_field(ofname,
"file name for root output") ;
57 zap[
"P"] =
make_field(Pscat,
"list of Pscat for each nscat") ;
58
59
60
61
62
63 if (zap.
read(argc, argv) != 0) {
64 return 1 ;
65 }
66 }
67 catch(const exception &error) {
68 cout << error.what() << endl ;
69 cout << "Type '" << argv[0] << " -h!' to display the command-line options." << endl ;
70 exit(1) ;
71 }
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
114 TFile* f = new TFile(ifname_model.c_str(),"read") ;
115 if( f->IsZombie() ) {
116 cerr << "Could not open file '" << ifname_model << "'." << endl ;
117 exit(1) ;
118 }
120 if( sm == NULL ) {
121 cerr << "Could not read JScatteringModel from file!" << endl ;
122 exit(1) ;
123 }
124 cout << "JScatteringModel loaded" << endl
125 << " absorption length = " << sm->getAbsorptionLength() << " m" << endl
127 << endl ;
128
129
130
132
133
134
135 cout << "Reading files." << endl ;
137 cout << "Reading '" << *it << "'." << endl ;
138
140 reader.open(it->c_str()) ;
141 if( !reader.is_open() ) {
142 cerr << "FATAL ERROR: unable to open input file '" << *it << "'." << endl ;
143 exit(1) ;
144 }
145 int nread = 0 ;
147 while( reader.hasNext() ) {
148 p = reader.next() ;
150 if( nscat >= (int)paths.size() ) paths.resize(nscat+1) ;
151 paths[nscat].push_back(*p) ;
152 ++nread ;
153 }
154 cout << "--> read " << nread << " paths." << endl ;
155 reader.close() ;
156 }
157 cout << endl ;
158 cout << "Done reading paths." << endl ;
159 cout << endl ;
160
161 const int nscatmax = paths.size() ;
162
163
164 Pscat.resize( nscatmax ) ;
166 for( int nscat=0 ; nscat<nscatmax ; ++nscat ) {
167 weight[nscat] = 0 ;
168 if( Pscat[nscat]>0 && paths[nscat].size()>0 ) {
169 weight[nscat] = Pscat[nscat] / paths[nscat].size() ;
170 }
171 }
172
173
174 cout << setw(5) << "nscat" << setw(15) << "npaths" << setw(15) << "Pscat" << endl ;
175 for( int nscat=0; nscat<nscatmax; ++nscat ) {
176 if( paths[nscat].size() > 0 ) {
177 cout << setw(5) << nscat
178 << setw(15) << paths[nscat].size()
179 << setw(15) << Pscat[nscat]
180 << endl ;
181 }
182 }
183
184 cout << "Determining histogram ranges." << endl ;
189
190 double INF = 1.0/0.0 ;
191 double MININF = -1.0/0.0 ;
192 for( int nscat=0 ; nscat<nscatmax ; ++nscat ) {
195 Lmax[nscat] = 0 ;
196 Lmin[nscat] = INF ;
197 int nvert = nscat + 2 ;
198 for(
int j=0 ;
j<nvert ; ++
j ) {
201 }
202 }
203
204
205 for( int nscat=0 ; nscat<(int)paths.size() ; ++nscat ) {
206
209 if( length > Lmax[nscat] ) Lmax[nscat] = length ;
210 if( length < Lmin[nscat] ) Lmin[nscat] = length ;
211
212 int nvert = nscat + 2 ;
213 for(
int j=0 ;
j<nvert ; ++
j ) {
215 min(posmin[nscat][
j].getY(),p->at(
j).getY()),
216 min(posmin[nscat][
j].getZ(),p->at(
j).getZ()) ) ;
218 max(posmax[nscat][
j].getY(),p->at(
j).getY()),
219 max(posmax[nscat][
j].getZ(),p->at(
j).getZ()) ) ;
220 }
221 }
222 }
223 cout << endl ;
224
225 cout << "Allocating histograms." << endl ;
226
238
239 int nbins_L = 10000 ;
240 int nbins_ct = 10000 ;
241 int nbins_X = 1000 ;
242 int nbins_Y = 1000 ;
243 int nbins_Z = 1000 ;
244
245 for( int nscat=0 ; nscat<nscatmax ; ++nscat ) {
246 if( paths[nscat].size() == 0 ) continue ;
247
248 char hname[200] ;
249 char htitle[200] ;
250
251 sprintf( hname, "hL_nscat_%i", nscat ) ;
252 sprintf( htitle, "path length distribution of paths with %i scatterings.", nscat ) ;
253 hL[nscat] = new TH1F(hname,htitle,nbins_L,Lmin[nscat],1.01*Lmax[nscat]) ;
254
255 sprintf( hname, "hL_zoom_nscat_%i", nscat ) ;
256 sprintf( htitle, "path length distribution of paths with %i scatterings.", nscat ) ;
257 hL_zoom[nscat] = new TH1F(hname,htitle,nbins_L,Lmin[nscat],Lmin[nscat]+10) ;
258
259 double ctmax = 1+1e-8 ;
260 sprintf( hname, "hCosThetaSource_nscat_%i", nscat ) ;
261 sprintf( htitle, "Source angle with %i scatterings;cos(#theta_{S});a.u.", nscat ) ;
262 hCosThetaSource[nscat] = new TH1F(hname,htitle,nbins_ct,-ctmax,ctmax) ;
263
264 sprintf( hname, "hCosThetaTarget_nscat_%i", nscat ) ;
265 sprintf( htitle, "Target angle with %i scatterings;cos(#theta_{T});a.u.", nscat ) ;
266 hCosThetaTarget[nscat] = new TH1F(hname,htitle,nbins_ct,-ctmax,ctmax) ;
267
268 sprintf( hname, "hCosThetaTarget_L_nscat_%i", nscat ) ;
269 sprintf( htitle, "%i scatterings;cos(#theta_{T});path length L", nscat ) ;
270 hCosThetaTarget_L[nscat] = new TH2F(hname,htitle,20,-ctmax,ctmax,200,Lmin[nscat],1.01*Lmax[nscat]) ;
271 hCosThetaTarget_L[nscat]->SetOption("colz") ;
272
273 int nvert = nscat + 2 ;
274 hX[nscat].resize(nvert) ;
275 hY[nscat].resize(nvert) ;
276 hZ[nscat].resize(nvert) ;
277 hXY[nscat].resize(nvert) ;
278 hXZ[nscat].resize(nvert) ;
279 hYZ[nscat].resize(nvert) ;
280 for(
int j=0 ;
j<nvert; ++
j ) {
281 sprintf(hname,
"hX_nscat_%i_vertex_%i", nscat,
j ) ;
282 sprintf(htitle,
"X of vertex %i in paths with %i scatterings.",
j, nscat ) ;
283 hX[nscat][
j] =
new TH1F(hname,htitle,nbins_X,1.01*posmin[nscat][
j].getX(),1.01*posmax[nscat][
j].getX()) ;
284
285 sprintf(hname,
"hY_nscat_%i_vertex_%i", nscat,
j ) ;
286 sprintf(htitle,
"Y of vertex %i in paths with %i scatterings.",
j, nscat ) ;
287 hY[nscat][
j] =
new TH1F(hname,htitle,nbins_Y,1.01*posmin[nscat][
j].getY(),1.01*posmax[nscat][
j].getY()) ;
288
289 sprintf(hname,
"hZ_nscat_%i_vertex_%i", nscat,
j ) ;
290 sprintf(htitle,
"Z of vertex %i in paths with %i scatterings.",
j, nscat ) ;
291 hZ[nscat][
j] =
new TH1F(hname,htitle,nbins_Z,1.01*posmin[nscat][
j].getZ(),1.01*posmax[nscat][
j].getZ()) ;
292
293 sprintf(hname,
"hXY_nscat_%i_vertex_%i", nscat,
j ) ;
294 sprintf(htitle,
"XY of vertex %i in paths with %i scatterings.",
j, nscat ) ;
295 hXY[nscat][
j] =
new TH2F(hname,htitle,
296 nbins_X,1.01*posmin[nscat][
j].getX(),1.01*posmax[nscat][
j].getX(),
297 nbins_Y,1.01*posmin[nscat][
j].getY(),1.01*posmax[nscat][
j].getY() ) ;
298 hXY[nscat][
j]->SetOption(
"colz") ;
299
300 sprintf(hname,
"hXZ_nscat_%i_vertex_%i", nscat,
j ) ;
301 sprintf(htitle,
"XZ of vertex %i in paths with %i scatterings.",
j, nscat ) ;
302 hXZ[nscat][
j] =
new TH2F(hname,htitle,
303 nbins_X,1.01*posmin[nscat][
j].getX(),1.01*posmax[nscat][
j].getX(),
304 nbins_Z,1.01*posmin[nscat][
j].getZ(),1.01*posmax[nscat][
j].getZ() ) ;
305 hXZ[nscat][
j]->SetOption(
"colz") ;
306
307 sprintf(hname,
"hYZ_nscat_%i_vertex_%i", nscat,
j ) ;
308 sprintf(htitle,
"YZ of vertex %i in paths with %i scatterings.",
j, nscat ) ;
309 hYZ[nscat][
j] =
new TH2F(hname,htitle,
310 nbins_Y,1.01*posmin[nscat][
j].getY(),1.01*posmax[nscat][
j].getY(),
311 nbins_Z,1.01*posmin[nscat][
j].getZ(),1.01*posmax[nscat][
j].getZ() ) ;
312 hYZ[nscat][
j]->SetOption(
"colz") ;
313 }
314 }
315 cout << endl ;
316
317
318 cout << "Filling histograms." << endl ;
319 for( int nscat=0 ; nscat<(int)paths.size() ; ++nscat ) {
320 int nvert = nscat + 2 ;
321
324 hL[nscat]->Fill(L) ;
325 hL_zoom[nscat]->Fill(L) ;
326
327 {
329 hCosThetaSource[nscat]->Fill(seg.getDZ()) ;
330 }
331
332 {
334 hCosThetaTarget[nscat]->Fill(seg.getDZ()) ;
335 hCosThetaTarget_L[nscat]->Fill(seg.getDZ(),L) ;
336 }
337
338 for(
int j=0;
j<nvert; ++
j ) {
339 double x = p->at(
j).getX() ;
340 double y = p->at(
j).getY() ;
341 double z = p->at(
j).getZ() ;
342 hX[nscat][
j]->Fill(x) ;
343 hY[nscat][
j]->Fill(y) ;
344 hZ[nscat][
j]->Fill(z) ;
345 hXY[nscat][
j]->Fill(x,y) ;
346 hXZ[nscat][
j]->Fill(x,z) ;
347 hYZ[nscat][
j]->Fill(y,z) ;
348 }
349 }
350 }
351 cout << endl ;
352
353
354 cout << "Normalizing histograms." << endl ;
355 for( int nscat=0 ; nscat<nscatmax ; ++nscat ) {
356 if( paths[nscat].size() == 0 ) continue ;
357 int nvert = nscat + 2 ;
358 hL[nscat]->Scale( weight[nscat] ) ;
359 hL_zoom[nscat]->Scale( weight[nscat] ) ;
360 hCosThetaSource[nscat]->Scale( weight[nscat] ) ;
361 hCosThetaTarget[nscat]->Scale( weight[nscat] ) ;
362 for(
int j=0 ;
j<nvert; ++
j ) {
363 hX[nscat][
j]->Scale( weight[nscat] ) ;
364 hY[nscat][
j]->Scale( weight[nscat] ) ;
365 hZ[nscat][
j]->Scale( weight[nscat] ) ;
366 hXY[nscat][
j]->Scale( weight[nscat] ) ;
367 hXZ[nscat][
j]->Scale( weight[nscat] ) ;
368 hYZ[nscat][
j]->Scale( weight[nscat] ) ;
369 }
370 }
371
372
373 cout << "Writing results to '" << ofname << "'." << endl ;
374 TFile* fout = new TFile(ofname.c_str(),"recreate") ;
375 for( int nscat=0; nscat<nscatmax; ++nscat ) {
376 if( paths[nscat].size() == 0 ) continue ;
377 int nvert = nscat + 2 ;
378
379 char dname[200] ;
380 sprintf( dname, "nscat%i", nscat ) ;
381 gDirectory->mkdir(dname)->cd() ;
382 hL[nscat]->Write() ;
383 hL_zoom[nscat]->Write() ;
384 hCosThetaSource[nscat]->Write() ;
385 hCosThetaTarget[nscat]->Write() ;
386 hCosThetaTarget_L[nscat]->Write() ;
387 for(
int j=0;
j<nvert ; ++
j ) {
388 hX[nscat][
j]->Write() ;
389 hY[nscat][
j]->Write() ;
390 hZ[nscat][
j]->Write() ;
391 hXY[nscat][
j]->Write() ;
392 hXZ[nscat][
j]->Write() ;
393 hYZ[nscat][
j]->Write() ;
394 }
395 fout->cd() ;
396 }
397 fout->Close() ;
398 delete fout ;
399
400 cout << "Output written to '" << ofname << "'." << endl ;
401
402 cout << "Done!" << endl ;
403
404 return 0 ;
405}
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Data structure for position in three dimensions.
Data structure for normalised vector in three dimensions.
double getLength()
get the total path length
Virtual base class for a scattering model.
virtual double getScatteringLength()
Utility class to parse command line options.
int read(const int argc, const char *const argv[])
Parse the program's command line options.