Jpp 20.0.0-rc.3
the software that should make you happy
Loading...
Searching...
No Matches
Functions
JMarkovPathReader.cc File Reference
#include <iostream>
#include <sstream>
#include <iomanip>
#include <vector>
#include <cmath>
#include <cstdlib>
#include "TRandom3.h"
#include "TFile.h"
#include "TPolyLine3D.h"
#include "TPolyMarker3D.h"
#include "TAxis3D.h"
#include "TView.h"
#include "TView3D.h"
#include "TCanvas.h"
#include "TPad.h"
#include "TH1F.h"
#include "TH2F.h"
#include "Jeep/JParser.hh"
#include "JMarkov/JPhotonPath.hh"
#include "JMarkov/JPhotonPathReader.hh"
#include "JMarkov/JScatteringModel.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 42 of file JMarkovPathReader.cc.

42 {
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}
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
JGEOMETRY3D::JPosition3D
Data structure for position in three dimensions.
Definition JPosition3D.hh:38
JGEOMETRY3D::JVersor3D
Data structure for normalised vector in three dimensions.
Definition JVersor3D.hh:28
JMARKOV::JPhotonPath
A photon path.
Definition JPhotonPath.hh:38
JMARKOV::JPhotonPath::getLength
double getLength()
get the total path length
Definition JPhotonPath.hh:104
JMARKOV::JScatteringModel
Virtual base class for a scattering model.
Definition JScatteringModel.hh:45
JMARKOV::JScatteringModel::getScatteringLength
virtual double getScatteringLength()
Definition JScatteringModel.hh:72
JPARSER::JParser
Utility class to parse command line options.
Definition JParser.hh:1698
JPARSER::JParser::read
int read(const int argc, const char *const argv[])
Parse the program's command line options.
Definition JParser.hh:1992
JTOOLS::j
int j
Definition JPolint.hh:801
Generated by doxygen 1.12.0