22 #include "TPolyLine3D.h"
23 #include "TPolyMarker3D.h"
47 int main(
int argc,
char** argv ) {
48 cout <<
"JMarkovPathIntegrator" << endl
49 <<
"Written by Martijn Jongen" << endl
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 ;
70 if (zap.
read(argc, argv) != 0) {
74 catch(
const exception &error) {
75 cout << error.what() << endl ;
76 cout <<
"Type '" << argv[0] <<
" -h!' to display the command-line options." << endl ;
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 ;
89 int nvert = nscat + 2 ;
92 cout <<
"Loading scattering model" << endl ;
95 TFile*
f =
new TFile(ifname_model.c_str(),
"read") ;
97 cerr <<
"Could not open file '" << ifname_model <<
"'." << endl ;
102 cerr <<
"Could not read JScatteringModel from file!" << endl ;
105 cout <<
"JScatteringModel loaded" << endl
106 <<
" absorption length = " << sm->getAbsorptionLength() <<
" m" << endl
113 reader.
open(ifname_paths.c_str()) ;
115 cerr <<
"FATAL ERROR: unable to open input file '" << ifname_paths <<
"'." << endl ;
120 cout <<
"Reading file" << endl ;
124 if( p->size() == (
unsigned int)nvert )
125 paths.push_back(*p) ;
127 cout << p->size() << endl ;
132 int npaths = paths.size() ;
134 cerr <<
"FATAL ERROR: No JPhotonPaths with nscat=" << nscat <<
" in '" << ifname_paths <<
"'." << endl ;
137 cout <<
"Done reading file. Selected " << npaths <<
" / " << nread <<
" JPhotonPaths." << endl ;
143 double INF = 1.0/0.0 ;
144 double MININF = -1.0/0.0 ;
146 for(
int i=0 ; i<nvert ; ++i ) {
154 for(
int i=0; i<nvert; ++i ) {
156 min(minvals[i].getY(),p->at(i).getY()),
157 min(minvals[i].getZ(),p->at(i).getZ()) ) ;
159 max(maxvals[i].getY(),p->at(i).getY()),
160 max(maxvals[i].getZ(),p->at(i).getZ()) ) ;
165 cout <<
"Allocating histograms." << endl ;
170 for(
int n=0;
n<nvert; ++
n ) {
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()) ;
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()) ;
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()) ;
193 cout <<
"Filling histograms" << endl ;
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() ) ;
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() ) ;
216 writer.
open(
"high_contribution_paths.paths") ;
220 cout <<
"Computing scattering probability" << endl ;
228 TH1F* hIntegralConts ;
229 TH1F* hIntegralConts_zoom ;
236 hZ[0]->GetMean(1) ) ;
239 hY[nvert-1]->GetMean(1),
240 hZ[nvert-1]->GetMean(1) ) ;
243 for(
int i=0 ; i<nsamples; ++i ) {
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) ;
259 double rho = sm->getRho(testpath) ;
262 integralConts[i] = rho/
w ;
265 if( rho/w > 0.001 ) {
276 writer.
put( testpath ) ;
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) ;
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) ;
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) ;
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) ;
313 for(
int i=0; i<nsamples; ++i ) {
314 hRhos->Fill( rhos[i] ) ;
316 hIntegralConts->Fill( integralConts[i] ) ;
317 hIntegralConts_zoom->Fill( integralConts[i] ) ;
322 for(
int i=0; i<nsamples; ++i ) {
323 var += (integralConts[i]-Pscat) * (integralConts[i]-Pscat) ;
328 cout <<
"sigma = " << sigma <<
" (standard deviation of integral contributions)" << endl ;
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 ;
336 cout <<
"MPI_PSCAT " << Pscat << endl
337 <<
"MPI_ERROR " << Pscat_err << endl
342 TFile* fout =
new TFile(ofname.c_str(),
"recreate") ;
344 sprintf( dname,
"nscat%i", nscat ) ;
345 gDirectory->mkdir(dname)->cd() ;
349 hIntegralConts->Write() ;
350 hIntegralConts_zoom->Write() ;
352 for(
int n=0;
n<nvert; ++
n ) {
360 cout <<
"Output written to '" << ofname <<
"'." << endl ;
364 cout <<
"Done!" << endl ;
Utility class to parse command line options.
int main(int argc, char *argv[])
virtual void open(const char *file_name) override
Open file.
virtual bool hasNext() override
Check availability of next element.
virtual const pointer_type & next() override
Get next element.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
virtual bool put(const T &object) override
Object output.
int read(const int argc, const char *const argv[])
Parse the program's command line options.
Utility class to parse command line options.
Virtual base class for a scattering model.
virtual void close()
Close file.
Data structure for position in three dimensions.
virtual double getScatteringLength()