/******************************************************************************* Markov path reader * \file Reads and analyzes photon path ensembles *******************************************************************************/ // C++ standard library #include #include #include #include #include #include // ROOT #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" // JPP #include "Jeep/JParser.hh" #include "JMarkov/JPhotonPath.hh" #include "JMarkov/JPhotonPathReader.hh" #include "JMarkov/JScatteringModel.hh" // namespaces using namespace std ; // this little line removes the make_field macro defined in JParser.hh // so that we can call the // JPARSER::make_field(T,string) function ourselves #undef make_field using namespace JPP; int main( int argc, char** argv ) { cout << "MarkovPathReader" << endl << "Written by Martijn Jongen" << endl << endl ; vector ifnames ; string ifname_model ; string ofname ; vector Pscat ; try { JParser zap ; // this argument parser can handle strings zap["f"] = make_field(ifnames,"input file names (binary files containing JPhotonPaths)") ; zap["m"] = make_field(ifname_model,"input file name (root file containing JScatteringModel)") ; zap["o"] = make_field(ofname,"file name for root output") ; zap["P"] = make_field(Pscat,"list of Pscat for each nscat") ; // NOTE: multiple files can be passed by calling e.g. // ./JMarkovPathreader [other arguments] -f "file1.paths // file2.paths // file3.paths" if (zap.read(argc, argv) != 0) { return 1 ; } } catch(const exception &error) { cout << error.what() << endl ; cout << "Type '" << argv[0] << " -h!' to display the command-line options." << endl ; exit(1) ; } // counter for number of paths of each length //map npaths ; // these will contain the minimal and maximal values of the vertex // coordinates // they are accessed as e.g. hX[nscat][vertex_index] /*vector< vector > Xmin ; // minimum x position of vertices vector< vector > Xmax ; // maximum x position of vertices vector< vector > Ymin ; // minimum y position of vertices vector< vector > Ymax ; // maximum y position of vertices vector< vector > Zmin ; // minimum z position of vertices vector< vector > Zmax ; // maximum z position of vertices vector Lmax ; // maximum path length Xmin.resize(nmax) ; Xmax.resize(nmax) ; Ymin.resize(nmax) ; Ymax.resize(nmax) ; Zmin.resize(nmax) ; Zmax.resize(nmax) ; Lmax.resize(nmax) ; for( int n=0 ; nIsZombie() ) { cerr << "Could not open file '" << ifname_model << "'." << endl ; exit(1) ; } sm = (JMARKOV::JScatteringModel*) f->Get("JMARKOV::JScatteringModel") ; if( sm == NULL ) { cerr << "Could not read JScatteringModel from file!" << endl ; exit(1) ; } cout << "JScatteringModel loaded" << endl << " absorption length = " << sm->getAbsorptionLength() << " m" << endl << " scattering length = " << sm->getScatteringLength() << " m" << endl << endl ; // read the photon paths into memory // we keep track of them in a vector of vectors (they are ordered by nscat) vector< vector > paths ; // read the file once to determine the minimal and maximal values // of the vertex coordinates cout << "Reading files." << endl ; for( vector::iterator it=ifnames.begin() ; it!=ifnames.end() ; ++it ) { cout << "Reading '" << *it << "'." << endl ; // read the paths from the file JMARKOV::JPhotonPathReader reader ; reader.open(it->c_str()) ; if( !reader.is_open() ) { cerr << "FATAL ERROR: unable to open input file '" << *it << "'." << endl ; exit(1) ; } int nread = 0 ; JMARKOV::JPhotonPath* p = NULL ; while( reader.hasNext() ) { p = reader.next() ; int nscat = p->n-2 ; //p->size()-1 ; if( nscat >= (int)paths.size() ) paths.resize(nscat+1) ; paths[nscat].push_back(*p) ; ++nread ; } cout << "--> read " << nread << " paths." << endl ; reader.close() ; } cout << endl ; cout << "Done reading paths." << endl ; cout << endl ; const int nscatmax = paths.size() ; // maximal number of scatterings Pscat.resize( nscatmax ) ; vector weight(nscatmax) ; for( int nscat=0 ; nscat0 && paths[nscat].size()>0 ) { weight[nscat] = Pscat[nscat] / paths[nscat].size() ; } } cout << setw(5) << "nscat" << setw(15) << "npaths" << setw(15) << "Pscat" << endl ; for( int nscat=0; nscat 0 ) { cout << setw(5) << nscat << setw(15) << paths[nscat].size() << setw(15) << Pscat[nscat] << endl ; } } cout << "Determining histogram ranges." << endl ; vector posmin(nscatmax) ; vector posmax(nscatmax) ; vector Lmin(nscatmax) ; // minimum path length vector Lmax(nscatmax) ; // maximum path length // initialize to +- infinity double INF = 1.0/0.0 ; double MININF = -1.0/0.0 ; for( int nscat=0 ; nscat::iterator p=paths[nscat].begin() ; p!=paths[nscat].end() ; ++p ) { double length = p->getLength() ; if( length > Lmax[nscat] ) Lmax[nscat] = length ; if( length < Lmin[nscat] ) Lmin[nscat] = length ; // loop over the vertices int nvert = nscat + 2 ; for( int j=0 ; jat(j).getX()), min(posmin[nscat][j].getY(),p->at(j).getY()), min(posmin[nscat][j].getZ(),p->at(j).getZ()) ) ; posmax[nscat][j] = JGEOMETRY3D::JPosition3D( max(posmax[nscat][j].getX(),p->at(j).getX()), max(posmax[nscat][j].getY(),p->at(j).getY()), max(posmax[nscat][j].getZ(),p->at(j).getZ()) ) ; } } } cout << endl ; cout << "Allocating histograms." << endl ; // create pointers to the histograms vector< vector > hX(nscatmax) ; // x position of vertices vector< vector > hY(nscatmax) ; // y position of vertices vector< vector > hZ(nscatmax) ; // z position of vertices vector< vector > hXY(nscatmax) ; // xy position of vertices vector< vector > hXZ(nscatmax) ; // xz position of vertices vector< vector > hYZ(nscatmax) ; // yz position of vertices vector hL(nscatmax) ; // path length vector hL_zoom(nscatmax) ; // path length vector hCosThetaSource(nscatmax) ; // source angle vector hCosThetaTarget(nscatmax) ; // target angle vector hCosThetaTarget_L(nscatmax) ; // target angle versus path length // number of bins int nbins_L = 10000 ; int nbins_ct = 10000 ; int nbins_X = 1000 ; int nbins_Y = 1000 ; int nbins_Z = 1000 ; // allocate the histograms for( int nscat=0 ; nscatSetOption("colz") ; int nvert = nscat + 2 ; hX[nscat].resize(nvert) ; hY[nscat].resize(nvert) ; hZ[nscat].resize(nvert) ; hXY[nscat].resize(nvert) ; hXZ[nscat].resize(nvert) ; hYZ[nscat].resize(nvert) ; for( int j=0 ; jSetOption("colz") ; sprintf(hname, "hXZ_nscat_%i_vertex_%i", nscat, j ) ; sprintf(htitle, "XZ of vertex %i in paths with %i scatterings.", j, nscat ) ; hXZ[nscat][j] = new TH2F(hname,htitle, nbins_X,1.01*posmin[nscat][j].getX(),1.01*posmax[nscat][j].getX(), nbins_Z,1.01*posmin[nscat][j].getZ(),1.01*posmax[nscat][j].getZ() ) ; hXZ[nscat][j]->SetOption("colz") ; sprintf(hname, "hYZ_nscat_%i_vertex_%i", nscat, j ) ; sprintf(htitle, "YZ of vertex %i in paths with %i scatterings.", j, nscat ) ; hYZ[nscat][j] = new TH2F(hname,htitle, nbins_Y,1.01*posmin[nscat][j].getY(),1.01*posmax[nscat][j].getY(), nbins_Z,1.01*posmin[nscat][j].getZ(),1.01*posmax[nscat][j].getZ() ) ; hYZ[nscat][j]->SetOption("colz") ; } } cout << endl ; // loop over the paths again, filling the histograms cout << "Filling histograms." << endl ; for( int nscat=0 ; nscat<(int)paths.size() ; ++nscat ) { int nvert = nscat + 2 ; // loop over paths with nscat scatterings for( vector::iterator p=paths[nscat].begin() ; p!=paths[nscat].end() ; ++p ) { double L = p->getLength() ; hL[nscat]->Fill(L) ; hL_zoom[nscat]->Fill(L) ; // source angle { JGEOMETRY3D::JVersor3D seg(p->at(1)-p->at(0)) ; hCosThetaSource[nscat]->Fill(seg.getDZ()) ; } // target angle { JGEOMETRY3D::JVersor3D seg(p->at(nvert-1)-p->at(nvert-2)) ; hCosThetaTarget[nscat]->Fill(seg.getDZ()) ; hCosThetaTarget_L[nscat]->Fill(seg.getDZ(),L) ; } // loop over vertices for( int j=0; jat(j).getX() ; double y = p->at(j).getY() ; double z = p->at(j).getZ() ; hX[nscat][j]->Fill(x) ; hY[nscat][j]->Fill(y) ; hZ[nscat][j]->Fill(z) ; hXY[nscat][j]->Fill(x,y) ; hXZ[nscat][j]->Fill(x,z) ; hYZ[nscat][j]->Fill(y,z) ; } } } cout << endl ; // normalize the histograms to 1 cout << "Normalizing histograms." << endl ; for( int nscat=0 ; nscatScale( weight[nscat] ) ; hL_zoom[nscat]->Scale( weight[nscat] ) ; hCosThetaSource[nscat]->Scale( weight[nscat] ) ; hCosThetaTarget[nscat]->Scale( weight[nscat] ) ; for( int j=0 ; jScale( weight[nscat] ) ; hY[nscat][j]->Scale( weight[nscat] ) ; hZ[nscat][j]->Scale( weight[nscat] ) ; hXY[nscat][j]->Scale( weight[nscat] ) ; hXZ[nscat][j]->Scale( weight[nscat] ) ; hYZ[nscat][j]->Scale( weight[nscat] ) ; } } // write results to a file cout << "Writing results to '" << ofname << "'." << endl ; TFile* fout = new TFile(ofname.c_str(),"recreate") ; for( int nscat=0; nscatmkdir(dname)->cd() ; hL[nscat]->Write() ; hL_zoom[nscat]->Write() ; hCosThetaSource[nscat]->Write() ; hCosThetaTarget[nscat]->Write() ; hCosThetaTarget_L[nscat]->Write() ; for( int j=0; jWrite() ; hY[nscat][j]->Write() ; hZ[nscat][j]->Write() ; hXY[nscat][j]->Write() ; hXZ[nscat][j]->Write() ; hYZ[nscat][j]->Write() ; } fout->cd() ; } fout->Close() ; delete fout ; cout << "Output written to '" << ofname << "'." << endl ; cout << "Done!" << endl ; return 0 ; }