Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JMarkovPathReader.cc
Go to the documentation of this file.
1 /*******************************************************************************
2  Markov path reader
3  * \file Reads and analyzes photon path ensembles
4 
5 *******************************************************************************/
6 // C++ standard library
7 #include<iostream>
8 #include<sstream>
9 #include<iomanip>
10 #include<vector>
11 #include<cmath>
12 #include<cstdlib>
13 
14 // ROOT
15 #include "TRandom3.h"
16 #include "TFile.h"
17 #include "TPolyLine3D.h"
18 #include "TPolyMarker3D.h"
19 #include "TAxis3D.h"
20 #include "TView.h"
21 #include "TView3D.h"
22 #include "TCanvas.h"
23 #include "TPad.h"
24 #include "TH1F.h"
25 #include "TH2F.h"
26 
27 // JPP
28 #include "Jeep/JParser.hh"
29 #include "JMarkov/JPhotonPath.hh"
30 #include "JMarkov/JPhotonPathReader.hh"
31 #include "JMarkov/JScatteringModel.hh"
32 
33 // namespaces
34 using namespace std ;
35 
36 // this little line removes the make_field macro defined in JParser.hh
37 // so that we can call the
38 // JPARSER::make_field(T,string) function ourselves
39 #undef make_field
40 using namespace JPP;
41 
42 int main( int argc, char** argv ) {
43  cout << "MarkovPathReader" << endl
44  << "Written by Martijn Jongen" << endl
45  << endl ;
46 
47  vector<string> ifnames ;
48  string ifname_model ;
49  string ofname ;
50  vector<double> Pscat ;
51 
52  try {
53  JParser<string> zap ; // this argument parser can handle strings
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  // NOTE: multiple files can be passed by calling e.g.
59  // ./JMarkovPathreader [other arguments] -f "file1.paths
60  // file2.paths
61  // file3.paths"
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  // counter for number of paths of each length
75  //map<int,int> npaths ;
76 
77  // these will contain the minimal and maximal values of the vertex
78  // coordinates
79  // they are accessed as e.g. hX[nscat][vertex_index]
80  /*vector< vector<double> > Xmin ; // minimum x position of vertices
81  vector< vector<double> > Xmax ; // maximum x position of vertices
82  vector< vector<double> > Ymin ; // minimum y position of vertices
83  vector< vector<double> > Ymax ; // maximum y position of vertices
84  vector< vector<double> > Zmin ; // minimum z position of vertices
85  vector< vector<double> > Zmax ; // maximum z position of vertices
86  vector<double> Lmax ; // maximum path length
87  Xmin.resize(nmax) ;
88  Xmax.resize(nmax) ;
89  Ymin.resize(nmax) ;
90  Ymax.resize(nmax) ;
91  Zmin.resize(nmax) ;
92  Zmax.resize(nmax) ;
93  Lmax.resize(nmax) ;
94  for( int n=0 ; n<nmax ; ++n ) {
95  Xmin[n].resize(n+2) ;
96  Xmax[n].resize(n+2) ;
97  Ymin[n].resize(n+2) ;
98  Ymax[n].resize(n+2) ;
99  Zmin[n].resize(n+2) ;
100  Zmax[n].resize(n+2) ;
101  Lmax[n] = -1.0/0.0 ;
102  for( int i=0 ; i<n+2 ; ++i ) {
103  Xmin[n][i] = 1.0/0.0 ;
104  Ymin[n][i] = 1.0/0.0 ;
105  Zmin[n][i] = 1.0/0.0 ;
106  Xmax[n][i] = -1.0/0.0 ;
107  Ymax[n][i] = -1.0/0.0 ;
108  Zmax[n][i] = -1.0/0.0 ;
109  }
110  }*/
111 
112  // read in the scattering model
113  JMARKOV::JScatteringModel* sm = NULL ;
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  }
119  sm = (JMARKOV::JScatteringModel*) f->Get("JMARKOV::JScatteringModel") ;
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
126  << " scattering length = " << sm->getScatteringLength() << " m" << endl
127  << endl ;
128 
129  // read the photon paths into memory
130  // we keep track of them in a vector of vectors (they are ordered by nscat)
131  vector< vector<JMARKOV::JPhotonPath> > paths ;
132 
133  // read the file once to determine the minimal and maximal values
134  // of the vertex coordinates
135  cout << "Reading files." << endl ;
136  for( vector<string>::iterator it=ifnames.begin() ; it!=ifnames.end() ; ++it ) {
137  cout << "Reading '" << *it << "'." << endl ;
138  // read the paths from the file
139  JMARKOV::JPhotonPathReader reader ;
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 ;
146  JMARKOV::JPhotonPath* p = NULL ;
147  while( reader.hasNext() ) {
148  p = reader.next() ;
149  int nscat = p->n-2 ; //p->size()-1 ;
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() ; // maximal number of scatterings
162 
163 
164  Pscat.resize( nscatmax ) ;
165  vector<double> weight(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 ;
185  vector<JMARKOV::JPhotonPath> posmin(nscatmax) ;
186  vector<JMARKOV::JPhotonPath> posmax(nscatmax) ;
187  vector<double> Lmin(nscatmax) ; // minimum path length
188  vector<double> Lmax(nscatmax) ; // maximum path length
189  // initialize to +- infinity
190  double INF = 1.0/0.0 ;
191  double MININF = -1.0/0.0 ;
192  for( int nscat=0 ; nscat<nscatmax ; ++nscat ) {
193  posmin[nscat] = JMARKOV::JPhotonPath(nscat) ;
194  posmax[nscat] = JMARKOV::JPhotonPath(nscat) ;
195  Lmax[nscat] = 0 ;
196  Lmin[nscat] = INF ;
197  int nvert = nscat + 2 ;
198  for( int j=0 ; j<nvert ; ++j ) {
199  posmin[nscat][j] = JGEOMETRY3D::JPosition3D( INF, INF, INF ) ;
200  posmax[nscat][j] = JGEOMETRY3D::JPosition3D( MININF, MININF, MININF ) ;
201  }
202  }
203 
204  // loop over the paths to find the histogram ranges
205  for( int nscat=0 ; nscat<(int)paths.size() ; ++nscat ) {
206  // loop over paths with nscat scatterings
207  for( vector<JMARKOV::JPhotonPath>::iterator p=paths[nscat].begin() ; p!=paths[nscat].end() ; ++p ) {
208  double length = p->getLength() ;
209  if( length > Lmax[nscat] ) Lmax[nscat] = length ;
210  if( length < Lmin[nscat] ) Lmin[nscat] = length ;
211  // loop over the vertices
212  int nvert = nscat + 2 ;
213  for( int j=0 ; j<nvert ; ++j ) {
214  posmin[nscat][j] = JGEOMETRY3D::JPosition3D( min(posmin[nscat][j].getX(),p->at(j).getX()),
215  min(posmin[nscat][j].getY(),p->at(j).getY()),
216  min(posmin[nscat][j].getZ(),p->at(j).getZ()) ) ;
217  posmax[nscat][j] = JGEOMETRY3D::JPosition3D( max(posmax[nscat][j].getX(),p->at(j).getX()),
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  // create pointers to the histograms
227  vector< vector<TH1F*> > hX(nscatmax) ; // x position of vertices
228  vector< vector<TH1F*> > hY(nscatmax) ; // y position of vertices
229  vector< vector<TH1F*> > hZ(nscatmax) ; // z position of vertices
230  vector< vector<TH2F*> > hXY(nscatmax) ; // xy position of vertices
231  vector< vector<TH2F*> > hXZ(nscatmax) ; // xz position of vertices
232  vector< vector<TH2F*> > hYZ(nscatmax) ; // yz position of vertices
233  vector<TH1F*> hL(nscatmax) ; // path length
234  vector<TH1F*> hL_zoom(nscatmax) ; // path length
235  vector<TH1F*> hCosThetaSource(nscatmax) ; // source angle
236  vector<TH1F*> hCosThetaTarget(nscatmax) ; // target angle
237  vector<TH2F*> hCosThetaTarget_L(nscatmax) ; // target angle versus path length
238  // number of bins
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  // allocate the histograms
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  // loop over the paths again, filling the histograms
318  cout << "Filling histograms." << endl ;
319  for( int nscat=0 ; nscat<(int)paths.size() ; ++nscat ) {
320  int nvert = nscat + 2 ;
321  // loop over paths with nscat scatterings
322  for( vector<JMARKOV::JPhotonPath>::iterator p=paths[nscat].begin() ; p!=paths[nscat].end() ; ++p ) {
323  double L = p->getLength() ;
324  hL[nscat]->Fill(L) ;
325  hL_zoom[nscat]->Fill(L) ;
326  // source angle
327  {
328  JGEOMETRY3D::JVersor3D seg(p->at(1)-p->at(0)) ;
329  hCosThetaSource[nscat]->Fill(seg.getDZ()) ;
330  }
331  // target angle
332  {
333  JGEOMETRY3D::JVersor3D seg(p->at(nvert-1)-p->at(nvert-2)) ;
334  hCosThetaTarget[nscat]->Fill(seg.getDZ()) ;
335  hCosThetaTarget_L[nscat]->Fill(seg.getDZ(),L) ;
336  }
337  // loop over vertices
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  // normalize the histograms to 1
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  // write results to a file
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 }
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1410
JLANG::JAbstractObjectIterator::next
virtual const pointer_type & next()
Get next element.
Definition: JAbstractObjectIterator.hh:71
JMARKOV::JPhotonPath
A photon path.
Definition: JPhotonPath.hh:38
JLANG::JAccessibleInputStream::close
virtual void close()
Close file.
Definition: JAccessibleStream.hh:77
JLANG::JAccessibleBinaryInputStream::open
virtual void open(const char *file_name)
Open file.
Definition: JAccessibleBinaryStream.hh:52
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1836
JLANG::JAccessibleInputStream::is_open
virtual bool is_open() const
Check is file is open.
Definition: JAccessibleStream.hh:53
JPARSER::JParser::read
int read(const int argc, const char *const argv[])
Parse the program&#39;s command line options.
Definition: JParser.hh:1673
JParser.hh
Utility class to parse command line options.
JMARKOV::JScatteringModel
Virtual base class for a scattering model.
Definition: JScatteringModel.hh:45
JLANG::JAbstractObjectIterator::hasNext
virtual bool hasNext()
Check availability of next element.
Definition: JAbstractObjectIterator.hh:56
JGEOMETRY3D::JPosition3D
Data structure for position in three dimensions.
Definition: JPosition3D.hh:35
JMARKOV::JScatteringModel::getScatteringLength
virtual double getScatteringLength()
Definition: JScatteringModel.hh:72
JGEOMETRY3D::JVersor3D
Data structure for normalised vector in three dimensions.
Definition: JVersor3D.hh:23
main
int main(int argc, char *argv[])
Definition: Main.cpp:15
Generated by   doxygen 1.8.5