43 cout <<
"MarkovPathReader" << endl
44 <<
"Written by Martijn Jongen" << endl
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") ;
63 if (zap.
read(argc, argv) != 0) {
67 catch(
const exception &error) {
68 cout << error.what() << endl ;
69 cout <<
"Type '" << argv[0] <<
" -h!' to display the command-line options." << endl ;
114 TFile*
f =
new TFile(ifname_model.c_str(),
"read") ;
115 if( f->IsZombie() ) {
116 cerr <<
"Could not open file '" << ifname_model <<
"'." << endl ;
121 cerr <<
"Could not read JScatteringModel from file!" << endl ;
124 cout <<
"JScatteringModel loaded" << endl
125 <<
" absorption length = " << sm->getAbsorptionLength() <<
" m" << endl
135 cout <<
"Reading files." << endl ;
137 cout <<
"Reading '" << *it <<
"'." << endl ;
140 reader.
open(it->c_str()) ;
142 cerr <<
"FATAL ERROR: unable to open input file '" << *it <<
"'." << endl ;
150 if( nscat >= (
int)paths.size() ) paths.resize(nscat+1) ;
151 paths[nscat].push_back(*p) ;
154 cout <<
"--> read " << nread <<
" paths." << endl ;
158 cout <<
"Done reading paths." << endl ;
161 const int nscatmax = paths.size() ;
164 Pscat.resize( nscatmax ) ;
166 for(
int nscat=0 ; nscat<nscatmax ; ++nscat ) {
168 if( Pscat[nscat]>0 && paths[nscat].size()>0 ) {
169 weight[nscat] = Pscat[nscat] / paths[nscat].size() ;
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]
184 cout <<
"Determining histogram ranges." << endl ;
190 double INF = 1.0/0.0 ;
191 double MININF = -1.0/0.0 ;
192 for(
int nscat=0 ; nscat<nscatmax ; ++nscat ) {
197 int nvert = nscat + 2 ;
198 for(
int j=0 ;
j<nvert ; ++
j ) {
205 for(
int nscat=0 ; nscat<(int)paths.size() ; ++nscat ) {
209 if( length > Lmax[nscat] ) Lmax[nscat] = length ;
210 if( length < Lmin[nscat] ) Lmin[nscat] = length ;
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()) ) ;
225 cout <<
"Allocating histograms." << endl ;
239 int nbins_L = 10000 ;
240 int nbins_ct = 10000 ;
245 for(
int nscat=0 ; nscat<nscatmax ; ++nscat ) {
246 if( paths[nscat].size() == 0 )
continue ;
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]) ;
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) ;
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) ;
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) ;
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") ;
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()) ;
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()) ;
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()) ;
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") ;
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") ;
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") ;
318 cout <<
"Filling histograms." << endl ;
319 for(
int nscat=0 ; nscat<(int)paths.size() ; ++nscat ) {
320 int nvert = nscat + 2 ;
325 hL_zoom[nscat]->Fill(L) ;
329 hCosThetaSource[nscat]->Fill(seg.getDZ()) ;
334 hCosThetaTarget[nscat]->Fill(seg.getDZ()) ;
335 hCosThetaTarget_L[nscat]->Fill(seg.getDZ(),L) ;
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) ;
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] ) ;
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 ;
380 sprintf( dname,
"nscat%i", nscat ) ;
381 gDirectory->mkdir(dname)->cd() ;
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() ;
400 cout <<
"Output written to '" << ofname <<
"'." << endl ;
402 cout <<
"Done!" << endl ;
Utility class to parse command line options.
virtual bool hasNext() override
Check availability of next element.
virtual const pointer_type & next() override
Get next element.
o $QUALITY_ROOT d $DEBUG!JPlot1D f
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
int read(const int argc, const char *const argv[])
Parse the program's command line options.
Virtual base class for a scattering model.
Data structure for position in three dimensions.
virtual double getScatteringLength()
Data structure for normalised vector in three dimensions.
double getLength()
get the total path length
std::vector< double > weight