Jpp test-rotations-old
the software that should make you happy
Loading...
Searching...
No Matches
Functions
JMarkovPathIntegrator.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 "Jeep/JParser.hh"
#include "JMarkov/JPhotonPath.hh"
#include "JMarkov/JPhotonPathReader.hh"
#include "JMarkov/JPhotonPathWriter.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 47 of file JMarkovPathIntegrator.cc.

47 {
48 cout << "JMarkovPathIntegrator" << endl
49 << "Written by Martijn Jongen" << endl
50 << endl ;
51
52 gRandom->SetSeed(0) ;
53
54 string ifname_paths ;
55 string ifname_model ;
56 string ofname ;
57 int nsamples ;
58 int nscat ;
59 int nbins ;
60
61 try {
62 JParser<string> zap ; // this argument parser can handle strings
63 zap["f"] = make_field(ifname_paths,"input file name (binary file containing JPhotonPaths)") ;
64 zap["m"] = make_field(ifname_model,"input file name (root file containing JScatteringModel)") ;
65 zap["o"] = make_field(ofname,"OPTIONAL file name for root output") = "" ;
66 zap["n"] = make_field(nscat,"number of scatterings") ;
67 zap["N"] = make_field(nsamples,"number of samples") = 10000 ;
68 zap["nbins"] = make_field(nbins,"nbins of vertex distribution histograms") = 1000 ;
69
70 if (zap.read(argc, argv) != 0) {
71 return 1 ;
72 }
73 }
74 catch(const exception &error) {
75 cout << error.what() << endl ;
76 cout << "Type '" << argv[0] << " -h!' to display the command-line options." << endl ;
77 exit(1) ;
78 }
79
80 cout << "SETTINGS: " << endl ;
81 cout << " Input file (paths) = '" << ifname_paths << "'" << endl ;
82 cout << " Input file (model) = '" << ifname_model << "'" << endl ;
83 cout << " Output file = '" << ofname << "'" << endl ;
84 cout << " nsamples = " << nsamples << endl ;
85 cout << " nscat = " << nscat << endl ;
86 cout << " nbins = " << nbins << endl ;
87 cout << endl ;
88
89 int nvert = nscat + 2 ; // the number of vertices in each path
90
91 // read in the scatteringmodel
92 cout << "Loading scattering model" << endl ;
93 JMARKOV::JScatteringModel* sm ;
94 // open the root file
95 TFile* f = new TFile(ifname_model.c_str(),"read") ;
96 if( f->IsZombie() ) {
97 cerr << "Could not open file '" << ifname_model << "'." << endl ;
98 exit(1) ;
99 }
100 sm = (JMARKOV::JScatteringModel*) f->Get("JMARKOV::JScatteringModel") ;
101 if( sm == NULL ) {
102 cerr << "Could not read JScatteringModel from file!" << endl ;
103 exit(1) ;
104 }
105 cout << "JScatteringModel loaded" << endl
106 << " absorption length = " << sm->getAbsorptionLength() << " m" << endl
107 << " scattering length = " << sm->getScatteringLength() << " m" << endl
108 << endl ;
109
110
111 // read the paths from the file
112 JMARKOV::JPhotonPathReader reader ;
113 reader.open(ifname_paths.c_str()) ;
114 if( !reader.is_open() ) {
115 cerr << "FATAL ERROR: unable to open input file '" << ifname_paths << "'." << endl ;
116 exit(1) ;
117 }
118 int nread = 0 ;
119 JMARKOV::JPhotonPath* p = NULL ;
120 cout << "Reading file" << endl ;
121 vector<JMARKOV::JPhotonPath> paths ;
122 while( reader.hasNext() ) {
123 p = reader.next() ;
124 if( p->size() == (unsigned int)nvert )
125 paths.push_back(*p) ;
126 else
127 cout << p->size() << endl ;
128 ++nread ;
129 }
130 reader.close() ;
131
132 int npaths = paths.size() ;
133 if( npaths == 0 ) {
134 cerr << "FATAL ERROR: No JPhotonPaths with nscat=" << nscat << " in '" << ifname_paths << "'." << endl ;
135 exit(1) ;
136 }
137 cout << "Done reading file. Selected " << npaths << " / " << nread << " JPhotonPaths." << endl ;
138 cout << endl ;
139
140 // find the minimal and maximal values of each coordinate of each vertex
141 JMARKOV::JPhotonPath minvals(nscat) ;
142 JMARKOV::JPhotonPath maxvals(nscat) ;
143 double INF = 1.0/0.0 ;
144 double MININF = -1.0/0.0 ;
145 // initialize to +- infinity
146 for( int i=0 ; i<nvert ; ++i ) {
147 minvals[i] = JGEOMETRY3D::JPosition3D( INF, INF, INF ) ;
148 maxvals[i] = JGEOMETRY3D::JPosition3D( MININF, MININF, MININF ) ;
149 }
150
151 // loop over the paths once, to get the minimal and maximal values of x, y and z for each vertex
152 for( vector<JMARKOV::JPhotonPath>::iterator p=paths.begin() ; p!=paths.end() ; ++p ) {
153 // loop over the vertices
154 for( int i=0; i<nvert; ++i ) {
155 minvals[i] = JGEOMETRY3D::JPosition3D( min(minvals[i].getX(),p->at(i).getX()),
156 min(minvals[i].getY(),p->at(i).getY()),
157 min(minvals[i].getZ(),p->at(i).getZ()) ) ;
158 maxvals[i] = JGEOMETRY3D::JPosition3D( max(maxvals[i].getX(),p->at(i).getX()),
159 max(maxvals[i].getY(),p->at(i).getY()),
160 max(maxvals[i].getZ(),p->at(i).getZ()) ) ;
161 }
162 }
163 cout << endl ;
164
165 cout << "Allocating histograms." << endl ;
166 // allocate histograms for the vertex positions
167 vector<TH1F*> hX(nvert) ; // x positions of vertices
168 vector<TH1F*> hY(nvert) ; // y positions of vertices
169 vector<TH1F*> hZ(nvert) ; // z positions of vertices
170 for( int n=0; n<nvert; ++n ) {
171 char hname[200] ;
172 char htitle[200] ;
173
174 // we multiply the min/max values by this factor to ensure that no
175 // entries end up in the overflow bin
176 double fac = 1.001 ;
177
178 sprintf(hname, "hX_nscat_%i_vertex_%i", nscat, n ) ;
179 sprintf(htitle, "X of vertex %i (nscat = %i)", n, nscat ) ;
180 hX[n] = new TH1F(hname,htitle,nbins,fac*minvals[n].getX(),fac*maxvals[n].getX()) ;
181
182 sprintf(hname, "hY_nscat_%i_vertex_%i", nscat, n ) ;
183 sprintf(htitle, "Y of vertex %i (nscat = %i)", n, nscat ) ;
184 hY[n] = new TH1F(hname,htitle,nbins,fac*minvals[n].getY(),fac*maxvals[n].getY()) ;
185
186 sprintf(hname, "hZ_nscat_%i_vertex_%i", nscat, n ) ;
187 sprintf(htitle, "Z of vertex %i (nscat = %i)", n, nscat ) ;
188 hZ[n] = new TH1F(hname,htitle,nbins,fac*minvals[n].getZ(),fac*maxvals[n].getZ()) ;
189 }
190 cout << endl ;
191
192 // fill the histograms
193 cout << "Filling histograms" << endl ;
194 for( vector<JMARKOV::JPhotonPath>::iterator p=paths.begin() ; p!=paths.end() ; ++p ) {
195 // loop over the vertices
196 for( int n=0; n<nvert; ++n ) {
197 hX[n]->Fill( p->at(n).getX() ) ;
198 hY[n]->Fill( p->at(n).getY() ) ;
199 hZ[n]->Fill( p->at(n).getZ() ) ;
200 }
201 }
202 cout << endl ;
203
204 // normalize the histograms to 1
205 cout << "Normalizing histograms." << endl ;
206 for( int n=0 ; n<nvert ; ++n ) {
207 hX[n]->Scale( 1.0 / hX[n]->Integral() ) ;
208 hY[n]->Scale( 1.0 / hY[n]->Integral() ) ;
209 hZ[n]->Scale( 1.0 / hZ[n]->Integral() ) ;
210 }
211 cout << endl ;
212
213 // #tmp temporarily save the paths with extremely high contributions
214 // to a file
215 JMARKOV::JPhotonPathWriter writer ;
216 writer.open("high_contribution_paths.paths") ;
217
218 // MC-integrate the path probability, using the vertex distributions as
219 // our sampling function
220 cout << "Computing scattering probability" << endl ;
221 double Pscat = 0 ;
222 double Pscat_err ;
223 vector<double> integralConts(nsamples) ;
224 vector<double> rhos(nsamples) ;
225 vector<double> ws(nsamples) ;
226 TH1F* hRhos ;
227 TH1F* hWs ;
228 TH1F* hIntegralConts ;
229 TH1F* hIntegralConts_zoom ;
230
231 // define test path
232 JMARKOV::JPhotonPath testpath(nscat) ;
233 // starting vertex
234 testpath[0] = JGEOMETRY3D::JPosition3D( hX[0]->GetMean(1),
235 hY[0]->GetMean(1),
236 hZ[0]->GetMean(1) ) ;
237 // ending vertex
238 testpath[nvert-1] = JGEOMETRY3D::JPosition3D( hX[nvert-1]->GetMean(1),
239 hY[nvert-1]->GetMean(1),
240 hZ[nvert-1]->GetMean(1) ) ;
241
242 // integrate over the free vertices
243 for( int i=0 ; i<nsamples; ++i ) {
244 double w = 1 ; // contribution to the integral
245 // loop over all vertices except the first and last one
246 for( int n=1 ; n<nvert-1 ; ++n ) {
247 double x = hX[n]->GetRandom() ;
248 double y = hY[n]->GetRandom() ;
249 double z = hZ[n]->GetRandom() ;
250 Int_t xbin = hX[n]->GetXaxis()->FindBin(x) ;
251 Int_t ybin = hY[n]->GetXaxis()->FindBin(y) ;
252 Int_t zbin = hZ[n]->GetXaxis()->FindBin(z) ;
253 double wx = hX[n]->GetBinContent(xbin) / hX[n]->GetBinWidth(xbin) ;
254 double wy = hY[n]->GetBinContent(ybin) / hY[n]->GetBinWidth(ybin) ;
255 double wz = hZ[n]->GetBinContent(zbin) / hZ[n]->GetBinWidth(zbin) ;
256 w *= wx * wy * wz ;
257 testpath[n] = JGEOMETRY3D::JPosition3D(x,y,z) ;
258 }
259 double rho = sm->getRho(testpath) ;
260 rhos[i] = rho ;
261 ws[i] = w ;
262 integralConts[i] = rho/w ;
263 Pscat += rho / w ;
264
265 if( rho/w > 0.001 ) {
266 /*cout << "Found high rho/w = " << rho/w << endl ;
267 cout << "rho = " << rho << ", w = " << w << endl ;
268 cout << "Rerunning calculation of rho in verbose mode" << endl ;
269 sm->getRho(testpath,true) ;
270 cout << "Path: " << endl ;
271 for( JMARKOV::JPhotonPath::iterator it=testpath.begin() ; it!=testpath.end() ; ++it ) {
272 cout << "( " << it->getX() << ", " << it->getY() << ", " << it->getZ() << " ) " ;
273 }
274 cout << endl ;
275 cout << endl ;*/
276 writer.put( testpath ) ;
277 }
278 }
279 Pscat /= nsamples ;
280 writer.close() ;
281
282 // statistical characteristics of the integral contributions
283 double sigma = 0 ; // standard deviation
284 {
285 double max ;
286 char hname[200] ;
287 char htitle[300] ;
288
289 // histogram of path probability densities (whole range)
290 max = *(max_element( rhos.begin(), rhos.end() )) ;
291 sprintf( hname, "hRhos_nscat%i", nscat ) ;
292 sprintf( htitle, "Path probability density for nscat = %i", nscat ) ;
293 hRhos = new TH1F(hname,htitle,100,0,1.01*max) ;
294
295 // histogram of path probability densities (whole range)
296 max = *(max_element( ws.begin(), ws.end() )) ;
297 sprintf( hname, "hWs_nscat%i", nscat ) ;
298 sprintf( htitle, "Path weights for nscat = %i", nscat ) ;
299 hWs = new TH1F(hname,htitle,100,0,1.01*max) ;
300
301 // histogram of integral contributions (whole range)
302 max = *(max_element( integralConts.begin(), integralConts.end() )) ;
303 sprintf( hname, "hIntegralConts_nscat%i", nscat ) ;
304 sprintf( htitle, "Contributions to the integral for nscat = %i", nscat ) ;
305 hIntegralConts = new TH1F(hname,htitle,100,0,1.01*max) ;
306
307 // histogram of integral contributions (up to 10*mean)
308 max = 10*Pscat ;
309 sprintf( hname, "hIntegralConts_nscat%i_zoom", nscat ) ;
310 sprintf( htitle, "Contributions to the integral for nscat = %i (zoom)", nscat ) ;
311 hIntegralConts_zoom = new TH1F(hname,htitle,100,0,1.01*max) ;
312
313 for( int i=0; i<nsamples; ++i ) {
314 hRhos->Fill( rhos[i] ) ;
315 hWs->Fill( ws[i] ) ;
316 hIntegralConts->Fill( integralConts[i] ) ;
317 hIntegralConts_zoom->Fill( integralConts[i] ) ;
318 }
319
320 // get sample variance
321 double var = 0 ;
322 for( int i=0; i<nsamples; ++i ) {
323 var += (integralConts[i]-Pscat) * (integralConts[i]-Pscat) ;
324 }
325 var /= nsamples ;
326 // this is now the estimated width of the distribution
327 sigma = sqrt(var) ;
328 cout << "sigma = " << sigma << " (standard deviation of integral contributions)" << endl ;
329 }
330 Pscat_err = sigma / sqrt(nsamples) ;
331 cout << "Pscat(" << nscat << " scatterings) = " << Pscat << endl ;
332 cout << "Error estimate 1: sigma/sqrt(N) = "
333 << sigma << " / sqrt(" << nsamples << ") = "
334 << Pscat_err << endl ;
335 // write some nice grep-able output
336 cout << "MPI_PSCAT " << Pscat << endl
337 << "MPI_ERROR " << Pscat_err << endl
338 << endl ;
339
340 // optionally write the histograms to a root file
341 if( ofname != "" ) {
342 TFile* fout = new TFile(ofname.c_str(),"recreate") ;
343 char dname[200] ;
344 sprintf( dname, "nscat%i", nscat ) ;
345 gDirectory->mkdir(dname)->cd() ;
346
347 hRhos->Write() ;
348 hWs->Write() ;
349 hIntegralConts->Write() ;
350 hIntegralConts_zoom->Write() ;
351
352 for( int n=0; n<nvert; ++n ) {
353 hX[n]->Write() ;
354 hY[n]->Write() ;
355 hZ[n]->Write() ;
356 }
357 fout->cd() ;
358 fout->Close() ;
359 delete fout ;
360 cout << "Output written to '" << ofname << "'." << endl ;
361 cout << endl ;
362 }
363
364 cout << "Done!" << endl ;
365 return 0 ;
366}
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
JLANG::JAccessibleBinaryOutputStream::open
virtual void open(const char *file_name) override
Open file.
Definition JAccessibleBinaryStream.hh:92
JLANG::JAccessibleOutputStream::close
virtual void close()
Close file.
Definition JAccessibleStream.hh:142
JLANG::JObjectOutput::put
virtual bool put(const T &object)=0
Object output.
JMARKOV::JPhotonPath
A photon path.
Definition JPhotonPath.hh:38
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
GAUSS_HERMITE::sigma
const double sigma[]
Definition JQuadrature.cc:74
JTOOLS::n
const int n
Definition JPolint.hh:791
Generated by doxygen 1.12.0