JBeaconSimulator.cc File Reference

#include <string>
#include "Jeep/JParser.hh"
#include "JDetector/JModuleRouter.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleLocation.hh"
#include "JGeometry3D/JAxis3D.hh"
#include "JMarkov/JScatteringModel.hh"
#include "JMarkov/JMarkovPhotonTracker.hh"
#include "JMarkov/JMarkovPathGenerator.hh"
#include "JMarkov/JMarkovIntegrator.hh"
#include "JMarkov/JPhotonPathWriter.hh"
#include "JPhysics/JPDFToolkit.hh"
#include "JPhysics/JDispersion.hh"
#include "JPhysics/KM3NeT.hh"
#include "TH1D.h"
#include "TRandom3.h"
#include "TFile.h"

Go to the source code of this file.

Classes
class	FreePathSolver
	Finds photon paths from a nanobeacon source to a target PMT that have a non-zero probability. More...

Functions
const JModule &	getModule (const JDetector &detector, const JModuleLocation &location)
	find module with a given string and floor number More...
double	getMean (vector< double > &v)
	get mean of vector content More...
double	getSigma (vector< double > &v)
	get standard deviation of vector content More...
void	scale (vector< double > &v, double c)
	scale vector content More...
void	printResult (vector< double > &v)
void	writePaths (string ofname, vector< JPhotonPath > paths)
JPosition3D	getDOMcenter (const JModule &dom)
	Returns a good guess of where the actual center of the DOM is based on the PMT positions and directions. More...
JAxis3D	getBeaconAxis (const JModule &dom, double r, double beacon_theta, double beacon_phi)
	Get the position and orientation of the nanobeacon on a DOM. More...
int	main (int argc, char **argv)

Variables
const double	JBSversion = 3.21

Function Documentation

getModule()

const JModule & getModule	(	const JDetector &	detector,
		const JModuleLocation &	location
	)

find module with a given string and floor number

Definition at line 928 of file JBeaconSimulator.cc.

                                                                                        {
  for( unsigned int i=0; i<detector.size(); ++i ) {
    if( detector[i].getLocation() == location ) {
      return detector[i] ;
    }
  }
  cerr << "FATAL ERROR: no DOM found with location " << location << endl ;
  exit(1) ;
}

getMean()

double getMean

(

vector< double > &

v

)

get mean of vector content

Definition at line 56 of file JBeaconSimulator.cc.

                                    {
  double sum = 0 ;
  for( vector<double>::iterator it=v.begin() ; it!=v.end() ; ++it ) 
    sum += *it ;
  return sum / v.size() ;
}

getSigma()

double getSigma

(

vector< double > &

v

)

get standard deviation of vector content

Definition at line 64 of file JBeaconSimulator.cc.

                                     {
  double mean = getMean(v) ;
  double varsum = 0 ;
  for( vector<double>::iterator it=v.begin() ; it!=v.end() ; ++it ) 
    varsum += (*it-mean)*(*it-mean) ;
  return sqrt( varsum / v.size() ) ;
}

scale()

void scale	(	vector< double > &	v,
		double	c
	)

scale vector content

Definition at line 73 of file JBeaconSimulator.cc.

                                          {
  for( vector<double>::iterator it=v.begin() ; it!=v.end() ; ++it )
    *it = c * (*it) ;
}

printResult()

void printResult

(

vector< double > &

v

)

Definition at line 78 of file JBeaconSimulator.cc.

                                      {
  double mean = getMean(v) ;
  double sigma = getSigma(v) ;
  cout << "Value = " << mean << " +- " << sigma/sqrt(1.0*v.size()) << endl ;
}

writePaths()

void writePaths	(	string	ofname,
		vector< JPhotonPath >	paths
	)

Definition at line 85 of file JBeaconSimulator.cc.

                                                            {
  JMARKOV::JPhotonPathWriter writer ;
  writer.open(ofname.c_str()) ;
  for( vector<JPhotonPath>::iterator it=paths.begin() ; it!=paths.end() ; ++it ) {
    writer.put( *it ) ;
  }
  writer.close() ;
  cout << "paths written to '" << ofname << "'." << endl ;
}

getDOMcenter()

JPosition3D getDOMcenter

(

const JModule &

dom

)

Returns a good guess of where the actual center of the DOM is based on the PMT positions and directions.

Definition at line 99 of file JBeaconSimulator.cc.

                                               {
  double r = 0 ;
  for( int pmt1=0 ; pmt1<31 ; ++pmt1 ) {
    for( int pmt2=pmt1+1; pmt2<31 ; ++pmt2 ) { 
      double num = (dom.getPMT(pmt1).getPosition()-dom.getPMT(pmt2).getPosition()).getLength() ;
      double den = (JVector3D(dom.getPMT(pmt1).getDirection())-JVector3D(dom.getPMT(pmt2).getDirection())).getLength() ;
      r += num/den ;
    }
  }
  r /= 0.5*31*30 ;
  JPosition3D center(0,0,0) ;
  for( int pmt=0 ; pmt<31 ; ++pmt ) {
    center.add( dom.getPMT(pmt).getPosition() - r*JVector3D(dom.getPMT(pmt).getDirection()) ) ;
  }
  center.div(31) ;
  return center ;
}

getBeaconAxis()

JAxis3D getBeaconAxis	(	const JModule &	dom,
		double	r,
		double	beacon_theta,
		double	beacon_phi
	)

Get the position and orientation of the nanobeacon on a DOM.

The local coordinate system of the DOM is deduced from the PMT positions.

r = DOM radius [m] beacon_theta = theta angle of the beacon in the DOM's coordinate system [rad] phi_theta = phi angle of the beacon in the DOM's coordinate system [rad]

Definition at line 128 of file JBeaconSimulator.cc.

                                                                                              {
  JPosition3D dom_center = getDOMcenter(dom) ; 
 
  // PMT 22 defines the negative z-direction of the DOM coordinate system
  const JVersor3D neg_zdir( dom.getPMT(22).getDirection() ) ;
  const JVersor3D zdir( -neg_zdir.getDX(), -neg_zdir.getDY(), -neg_zdir.getDZ() ) ;
  // PMT 0 points in the negative y-direction of the DOM coordinate system (+some z-component)
  JVector3D tmp ;
  const JVersor3D xdir(  tmp.cross( zdir, dom.getPMT(0).getDirection() ) ) ;
  const JVersor3D ydir(  tmp.cross( zdir, xdir )  ) ;
 
  // direction of beacon's position w.r.t. the DOM center
  JVector3D _bdir( cos(beacon_theta)*JVector3D(zdir) 
                   + sin(beacon_theta)*cos(beacon_phi)*JVector3D(xdir)
                   + sin(beacon_theta)*sin(beacon_phi)*JVector3D(ydir)
                   ) ;
  JVersor3D bdir(_bdir) ; // start out in DOM z-direction
 
  JPosition3D bpos( dom_center + r*JVector3D(bdir) ) ;
  JAxis3D ax(bpos,zdir) ;
 
  return ax ;
}

main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 389 of file JBeaconSimulator.cc.

                                   {
  cout << "JBeaconSimulator version " << JBSversion << endl
       << endl ;
 
  string ofname ;
  string detectorFile ;
  
  JModuleLocation source_location ;
  vector<int> target_pmt ;
  bool useTracker ;
  unsigned int nscatmax ;
  unsigned long int nphotons ;
  double Labs ;
  double opening_angle_deg ;
  double beacon_theta_deg ;
  double beacon_phi_deg ;
  double lHG ;
  double gHG ;
  double lR  ;
  double aR  ;
  vector<double> stepsizes ;
  double target_step_size_deg ;
  double tangential_step_size_deg ;
  const vector<double> empty_vector ;
  int nsamples ;
  bool write_paths ;
  bool no_remapping ;
 
  try { 
    JParser<string> zap;
    
    zap["o"] = make_field(ofname) ; // output file name
    zap["a"] = make_field(detectorFile) ;
    zap["s"] = make_field(source_location) ;
    zap["t"] = make_field(target_pmt) ;
    zap["m"] = make_field(nscatmax) = 5 ;
    zap["N"] = make_field(nphotons) ;
    zap["N-integration-samples"] = make_field(nsamples) = 1e6 ;
    zap["tracker"] = make_field(useTracker) ;
    zap["write-paths"] = make_field(write_paths) ;
    zap["Labs"] = make_field(Labs) = 50 ;
    zap["opening-angle"] = make_field(opening_angle_deg) = 45 ;
    zap["beacon-theta-deg"] = make_field(beacon_theta_deg) = 56.2893-15 ; // 15 degrees above theta of upper row PMTs
    zap["beacon-phi-deg"] = make_field(beacon_phi_deg) = -59.9998 ; // phi of PMT 3
    zap["LscatHG"] = make_field(lHG) = 60.241 ;
    zap["gHG"] = make_field(gHG) = 0.924 ;
    zap["LscatR"] = make_field(lR) = 294.118 ;
    zap["aR"] = make_field(aR) = 0.853 ;
    zap["stepsizes"] = make_field(stepsizes) = empty_vector ;
    zap["target-step-size-deg"] = make_field(target_step_size_deg) = 1.0 ;
    zap["tangential-step-size-deg"] = make_field(tangential_step_size_deg) = 1.0 ;
    zap["no-remapping"] = make_field(no_remapping) ;
 
    if (zap.read(argc, argv) != 0) exit(1) ;
  }
  catch(const exception &error) {
    cerr << error.what() << endl ;
    exit(1) ;
  }
 
  // remove extension from output file name
  {
    vector<string> extensions ;
    extensions.push_back(".root")  ;
    extensions.push_back(".paths") ;
    for( vector<string>::iterator it=extensions.begin(); it!=extensions.end(); ++it ) {
      size_t pos = ofname.find(*it) ;
      if( pos != string::npos ) {
        ofname = ofname.substr(0,pos) ;
      }
    }
  }
 
  // use default step sizes when no new ones (or too few are specified)
  if( !useTracker ) {
    const int ndefvals = 6 ;
    double defvals[ndefvals] = { 8, 5, 4, 3, 2, 1 } ;
    while( stepsizes.size()<nscatmax ) {
      double val = defvals[min(ndefvals-1,(int)stepsizes.size())] ;
      stepsizes.push_back(val) ;
    }
  }
  
  // load detector
  JDetector detector;
  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    cout << "FATAL ERROR: could not open detector file '" << detectorFile << "'." ;
    exit(1) ;
  }
  JModuleRouter moduleRouter(detector) ;
 
 
  // initialize gRandom
  gRandom = new TRandom3(0) ;
 
  // radius of a DOM (used for source and target model)
  const double DOM_radius = 0.2159 ; // 0.2159 [m] (=17 inch diameter glass sphere)
 
 
  // find source position and direction (stored as a JAxis3D)
  // the source is a NB, whose position is deduced from the PMT positions in the detector file
  // note that the distance from the DOM center is DOM_radius + 0.01 mm, to make sure that the 
  // initial vertex does not intersect the DOM (which is a problem when the source is located
  // on the target DOM)
  JModule source_module( getModule(detector,source_location) ) ;
  JAxis3D source_axis( getBeaconAxis(source_module,DOM_radius+0.00001,beacon_theta_deg*M_PI/180,beacon_phi_deg*M_PI/180) ) ;
  JVersor3D surface_normal( source_axis.getPosition() - getDOMcenter(source_module) ) ;
 
  cout << "SOURCE (NANOBEACON)" << endl 
       << "Positioned on string " << source_location.getString() << " floor " << source_location.getFloor()
       << " (DOM center @ " << getDOMcenter(source_module) << ")" << endl
       << "Position in local DOM coordinates: theta = " 
       << beacon_theta_deg << " deg, "
       << "phi = " << beacon_phi_deg << " deg" << endl 
       << "distance w.r.t. DOM center = " << (source_axis.getPosition()-getDOMcenter(source_module)).getLength() << " m" << endl 
       << "unit vector from DOM center to beacon position = " << JVersor3D(source_axis.getPosition()-getDOMcenter(source_module)) << endl
       << "beam direction = " << source_axis.getDirection() << endl
       << "beam opening angle = " << opening_angle_deg << " deg" << endl 
       << "surface normal     = " << surface_normal << endl
       << endl ;
 
  // find target position and direction (stored as a JAxis3D)
  if( target_pmt.size() != 3 ) {
    cerr << "FATAL ERROR: provide a target PMT in the format "
         << "\"string floor daq_channel\"." << endl ;
    exit(1) ;
  }
 
  // require valid PMT index for the MCMC
  if( !useTracker ) {
    if( target_pmt[2]<0 || target_pmt[2]>30 ) {
      cerr << "FATAL ERROR: invalid target pmt index (" 
           << target_pmt[2] << ")" << endl ;
      exit(2) ;
    }
  }
 
  JModuleLocation target_location(target_pmt[0],target_pmt[1]) ;
  JModule target_module( getModule(detector,target_location) ) ;
  JAxis3D target_axis( getDOMcenter(target_module), target_module.getPMT(target_pmt[2]).getDirection() ) ;
 
  // source model
  JDirectedSource* src = new JDirectedSource( source_axis.getDirection(), 
                                              opening_angle_deg*M_PI/180,
                                              true,
                                              surface_normal ) ;
  src->setPosition( source_axis.getPosition() ) ;
 
  // initialize scattering model
  cout << "PHYSICS" << endl ;
  JHenyeyGreensteinScattering sm1( lHG, gHG ) ; 
  cout << "Henyey-Greenstein scattering length = " << lHG << " [m], g = " << gHG << endl ;
  JRayleighScattering sm2( lR, aR ) ;
  cout << "Rayleigh scattering length = " << lR << " [m], a = " << aR << endl ;
  JScatteringModel* sm = new JCombinedScattering( &sm1, &sm2 ) ;
  cout << "Absorption length = " << Labs << " [m]" << endl ;
 
  // water properties
  // for now this is not customizable
  const double water_depth = 3000 ; // m
  const double pressure    = 0.1*water_depth ; // atm
  JDispersion jd(pressure) ;
  const double nbwavelength = 470 ; // nm
  const double index_of_refraction = jd.getIndexOfRefractionGroup(nbwavelength) ;
  const double cw = JTOOLS::getSpeedOfLight() / index_of_refraction ;
  cout << "Assuming speed of light in water = " << cw << " [m/ns]" << endl ;
  cout << endl ;
 
 
  // storage container for generated photon paths
  vector<JPhotonPath> paths ;
 
  // target model (for now this is not customizable)
  // we use the approximation that a PMT is a spherical cap on a DOM, with a uniform
  // acceptance, i.e. the photon always causes a hit when it impacts the cap
  cout << "TARGET" << endl ;
  const double pmt_opening_angle_deg = 30 ;
  const double pmt_opening_angle_rad = pmt_opening_angle_deg*M_PI/180.0 ;
  const double pmt_ctmin = cos(0.5*pmt_opening_angle_rad) ;
  cout << "DOM: string " << target_pmt[0] 
       << ", floor " << target_pmt[1] 
       << ", @" << getDOMcenter(target_module) << endl
       << "Target radius = " << DOM_radius << " [m]" << endl 
       << "PMT: " ;
  if( useTracker ) cout << "all 31 PMTs" ;
  else             cout << "PMT " << target_pmt[2] << ", direction = " << target_axis.getDirection() ;
  cout << endl
       << "PMTs are simulated as spherical caps" << endl
       << "PMT opening angle = " << pmt_opening_angle_deg << " deg." << endl
       << endl ;
 
  JTargetModel* trg ;
  if( useTracker ) {
    // for the tracker we use an isotropic sphere that represents an entire DOM
    // later on, we will divide the hits over the different PMTs
    trg = new JIsotropicTarget() ; 
    trg->setRadius(DOM_radius) ;
    trg->setPosition(target_axis.getPosition()) ;
  } else {
    // for the MCMC we use a single PMT target
    trg = new JPMTTarget(target_axis,pmt_opening_angle_rad,DOM_radius) ;
  }
 
  unsigned long int n_generated_photons = nphotons ;
  if( useTracker ) { // generate paths by photon tracking
    //n_generated_photons = 1000e6 ;
    cout << "Generating photon paths by photon tracking" << endl ;
    JMarkovPhotonTracker mpt ;
    paths = mpt.trackPhotons(
         n_generated_photons,
         src, sm, trg,
         Labs ) ;
    cout << "Photon hits: " << paths.size() << " / " << n_generated_photons << endl ;
    cout << endl ;
  }
 
  vector<double> Pscat(nscatmax+1,1) ;
  vector<double> Perr(nscatmax+1,1) ;
  vector<JMarkovEnsembleIntegrator1D*> integrators ;
  if( !useTracker ) { // generate paths by MCMC
    JMarkovPathGenerator mpg ;
    if( no_remapping ) {
      // deactivate coordinate remapping.
      // That will most likely lead to trouble with the 1/r^2 singularities 
      // in the path probability density, meaning the space will be undersampled
      // for cases where vertices are very close together.
      mpg.setCoordinateRemapping(false) ;
      cout << "WARNING: coordinate remapping is deactivated!" << endl ;
    }
    mpg.setTargetStepSize_deg(target_step_size_deg) ;
    mpg.setTangentialStepSize_deg(tangential_step_size_deg) ;
    int nsteps_save = 250 ;
    int nsteps_burn_in = 10*nsteps_save ;
    for( int nscat=0 ; nscat<=(int)nscatmax ; ++nscat ) {
      cout << "nscat = " << nscat << endl ;
 
      double stepsize = 1 ; // default step size
      if( nscat>0 && (int)stepsizes.size()>=nscat ) stepsize = stepsizes[nscat-1] ;
      mpg.setRadialStepSize_m(stepsize) ;
 
      if( nscat != 0 ) cout << "step size = " << stepsize << " [m]" << endl ;
 
       // create start_path
      FreePathSolver fps ;
      JPhotonPath start_path = fps.solve(nscat,src,(JPMTTarget*)trg) ;
      cout << "path = " << endl ;
      cout << start_path.getString() << endl ;
      double start_rho = getPhotonPathProbabilityDensity(start_path,src,sm,trg,Labs) ;
      cout << "start_rho = " << start_rho << endl ;
      if( start_rho == 0 ) {
        cout << "No possible photon path found with nscat = " << nscat <<". Skipping nscat = " << nscat << "." << endl 
             << "WARNING: note that this does not necessarily mean there are no such paths, simply that the search algorithm could not find one." << endl 
             << endl;
        continue ;
      }
 
      /*cout << "Found the following start path: " << endl ;
      for( vector<JPosition3D>::iterator it=start_path.begin() ; it!= start_path.end() ; ++it ) {
        cout << *it << endl ;
      }
      cout << endl ;*/
 
      vector<JPhotonPath> these_paths = mpg.generateEnsemble( 
                   n_generated_photons,
                   start_path,
                   src, sm, trg,
                   Labs,
                   nsteps_burn_in, nsteps_save ) ;
      int nsteps = mpg.getNsteps() ;
      int naccepted = mpg.getNacceptedSteps() ;
      if( nscat != 0 ) {
        cout << "Accepted " << naccepted << " / " << nsteps << " steps (" << mpg.getFractionAccepted()*100.0 << "%)" << endl 
             << "Accepted " << 100.0*mpg.getFractionAccepted_radial() << "% of the radial steps." << endl
             << "Accepted " << 100.0*mpg.getFractionAccepted_tangential() << "% of the tangential steps." << endl ;
        cout << "NOTE: as a rule of thumb, the optimal value is around ~23% (this can be achieved by modifying the various step sizes)" << endl ;
      }
      paths.insert(paths.end(), these_paths.begin(), these_paths.end()) ;
      cout << endl ;
 
      if( these_paths.size() == 0 ) {
        cerr << "Something went wrong. The MCMC generator returned no paths." << endl ;
        Pscat[nscat] = 0.0 ;
        Perr[nscat]  = 0.0 ;
      } else {
        cout << "Creating JMarkovEnsembleIntegrator" << endl ;
        JMarkovEnsembleIntegrator1D* jmi = new JMarkovEnsembleIntegrator1D(these_paths,trg,100,100,200) ;
 
        cout << "Integrating" << endl ;
        vector<double> result = jmi->integrate(nsamples,nscat,src,sm,trg,Labs) ;
        printResult(result) ;
        Pscat[nscat] = getMean(result) ;
        Perr[nscat]  = getSigma(result) / sqrt(nsamples) ;
        //delete jmi ;
        integrators.push_back(jmi) ;
      }
      cout << endl ;
    }
  }
 
  if( paths.size()==0 ) {
    cerr << "FATAL ERROR: no paths were generated!" << endl ;
    exit(1) ;
  }
 
 
  // get path lengths and directions into an array
  // also find what the maximal number of scatterings is
  unsigned int npaths = paths.size() ;
  vector<double> path_lengths(npaths) ;
  vector<double> xs(npaths) ;
  vector<double> ys(npaths) ;
  vector<double> zs(npaths) ;
  for( unsigned int i=0; i<npaths; ++i ) {
    path_lengths[i] = paths[i].getLength() ;
    xs[i] = paths[i].back().getX() ;
    ys[i] = paths[i].back().getY() ;
    zs[i] = paths[i].back().getZ() ;
    int nscat = paths[i].size()-2 ;
    if( nscat>(int)nscatmax ) nscatmax = nscat ;
  }
 
  // find min and max x, y and z
  double xmin = trg->getPosition().getX() - 1.1*DOM_radius ; // *min_element(xs.begin(),xs.end() ) ;
  double xmax = trg->getPosition().getX() + 1.1*DOM_radius ; //*max_element(xs.begin(),xs.end() ) ;
  double ymin = trg->getPosition().getY() - 1.1*DOM_radius ; //*min_element(ys.begin(),ys.end() ) ;
  double ymax = trg->getPosition().getY() + 1.1*DOM_radius ; //*max_element(ys.begin(),ys.end() ) ;
  double zmin = trg->getPosition().getZ() - 1.1*DOM_radius ; //*min_element(zs.begin(),zs.end() ) ;
  double zmax = trg->getPosition().getZ() + 1.1*DOM_radius ; //*max_element(zs.begin(),zs.end() ) ;
 
  // find out which PMT is hit
  vector<int> PMT_is_hit ;
 
  if( !useTracker ) {
    // in the MCMC we know which PMT is hit
    PMT_is_hit = vector<int>(npaths,target_pmt[2]) ;
  }
 
  if( useTracker ) {
    PMT_is_hit = vector<int>(npaths,-1) ;
    vector<unsigned int> nHitsPerPMT(31,0) ;
    // in the tracker we have to figure out for every hit whether the PMT was hit or not
    for( unsigned int i=0; i<npaths; ++i ) {
 
      // find out which PMTs were hit (should be only one)
      vector<int> hitPMTs ;
      JVersor3D hit_direction = JVersor3D( paths[i].back() - target_axis.getPosition() ) ;
      for( unsigned int pmt=0 ; pmt<31 ; ++pmt ) {
        double ct = target_module.getPMT(pmt).getDirection().getDot(hit_direction) ;
        if( ct>pmt_ctmin ) hitPMTs.push_back(pmt) ;
      }
      if( hitPMTs.size()==1 ) {
        PMT_is_hit[i] = hitPMTs[0] ;
        ++nHitsPerPMT[ hitPMTs[0] ] ;
      }
      if( hitPMTs.size()>1 ) {
        cerr << "FATAL ERROR: hit " << hitPMTs.size() << " PMTs simultaneously" << endl ;
        exit(1) ;
      }
    }
    // print number of hits separated by PMT
    for( int pmt=0 ; pmt<31 ; ++pmt ) {
      if( nHitsPerPMT[pmt]>0 ) {
        cout << "PMT" << setw(2) << pmt << ": " << setw(15) << nHitsPerPMT[pmt] << " hits" << endl ;
      }
    }
    cout << endl ;
  }
 
  // range and binning for pulse profile(s)
  const double dt = 0.1 ; // bin width in ns 
  const unsigned int margin = 5 ; // ns
  double Lmin = *min_element(path_lengths.begin(),path_lengths.end() ) ;
  double Lmax = *max_element(path_lengths.begin(),path_lengths.end() ) ;
  double tmin = floor( Lmin / cw - margin ) ;
  double tmax = ceil( Lmax / cw + margin )  ;
  int nbinst = (int) round( (tmax-tmin)/dt ) ;
 
  int pmtfront = 0  ;
  int pmtback  = 31 ;
  if( !useTracker ) {
    // in the MCMC there is only one PMT
    pmtfront = target_pmt[2]   ;
    pmtback  = target_pmt[2]+1 ;
  }
 
  for( int pmt=pmtfront ; pmt<pmtback ; ++pmt ) {
    // open file
    char fname[500] ;
    sprintf( fname, "%s_PMT%i.root", ofname.c_str(), pmt ) ;
    TFile* f = new TFile(fname,"recreate") ;
 
    // create histograms
    TH1D hpulse("hpulse","Instantaneous flash pulse profile;time since flash [ns];hit probability / ns",
                nbinst,tmin,tmax) ;
    vector<TH1D*> hpulse_per_nscat(nscatmax+1) ;
    for( int nscat=0 ; nscat<=(int)nscatmax ; ++nscat ) {
      char hname[200] ;
      sprintf(hname,"hpulse_nscat%i",nscat) ;
      hpulse_per_nscat[nscat] = new TH1D(hname,"Instantaneous flash pulse profile;time since flash [ns];hit probability / ns",
                                         nbinst,tmin,tmax ) ;
    }
    TH2D hXY("hXY","hit position;X;Y",
             100,xmin,xmax,
             100,ymin,ymax ) ;
    hXY.SetOption("colz") ;
    TH2D hXZ("hXZ","hit position;X;Z",
             100,xmin,xmax,
             100,zmin,zmax ) ;
    hXZ.SetOption("colz") ;
    TH2D hYZ("hYZ","hit position;Y;Z",
             100,ymin,ymax,
             100,zmin,zmax ) ;
    hYZ.SetOption("colz") ;
    TH1D hnscat("hnscat",";N_{scat};probability",nscatmax+2,-0.5,nscatmax+1.5) ;
 
    TH3D hregion("hregion","Photon paths;X;Y;Z",
                 100,
                 trg->getPosition().getX() - 2*DOM_radius,
                 trg->getPosition().getX() + 2*DOM_radius,
                 100,
                 trg->getPosition().getY() - 2*DOM_radius,
                 trg->getPosition().getY() + 2*DOM_radius,
                 100,
                 trg->getPosition().getZ() - 2*DOM_radius,
                 trg->getPosition().getZ() + 2*DOM_radius) ;
 
    // fill histograms
    vector<JPhotonPath> paths_to_write ;
    for( unsigned int i=0 ; i<npaths ; ++i ) {
      if( PMT_is_hit[i] != pmt ) continue ;
      
      paths_to_write.push_back(paths[i]) ;
 
      const double path_length = paths[i].getLength() ;
      double t = path_length/cw ;
      int nscat = paths[i].size()-2 ;
 
      double w = 1 ;
      if( useTracker ) w = 1.0 / n_generated_photons ;
      else             w = Pscat[nscat] / n_generated_photons ; 
 
      hpulse.Fill(t,w) ;
      hnscat.Fill(nscat,w) ;
      hpulse_per_nscat[nscat]->Fill(t,w) ;
      
      hXY.Fill(xs[i],ys[i],w) ;
      hXZ.Fill(xs[i],zs[i],w) ;
      hYZ.Fill(ys[i],zs[i],w) ;
 
      if( paths.size()<10000 ) {
        const double dl = 0.001 ; // 1 mm
        const double dw = w*dl/path_length ;
        double l_along_seg = 0.5*dl ;
        int nfilled = 0 ;
        for( int nv=0; nv<(int)paths[i].size()-1; ++nv ) {
          double seg_length = paths[i][nv+1].getDistance(paths[i][nv]) ;
          JVersor3D seg_dir( paths[i][nv+1] - paths[i][nv] ) ;
          while( l_along_seg<seg_length ) {
            JPosition3D pos( paths[i][nv] + l_along_seg*JVector3D(seg_dir) ) ;
            hregion.Fill( pos.getX(), pos.getY(), pos.getZ(), dw ) ;
            if( pos.getDistance( trg->getPosition() ) < trg->getRadius() ) {
              cerr << "Woohoooo" << endl ;
              exit(1) ;
            }
            ++nfilled ;
            l_along_seg += dl ;
          }
          l_along_seg -= seg_length ;
        }
      }
    }
    // scale hregion histogram, so that SetShowProjection3D will no show an
    // empty histogram
    hregion.Scale(1e20) ;
 
    // write paths 
    if( write_paths ) {
      char path_fname[500] ;
      sprintf( path_fname, "%s_PMT%i.paths", ofname.c_str(), pmt ) ;
      cout << "Writing paths." << endl ;
      writePaths( path_fname, paths_to_write ) ;
      cout << endl ;
    }
 
    // for the tracker, our error is just the Poisson error
    if( useTracker ) {
      for( Int_t bin=1 ; bin<=hnscat.GetNbinsX() ; ++bin ) {
        double val = hnscat.GetBinContent(bin) ;
        double nentries = max(1.0, val*n_generated_photons ) ;
        hnscat.SetBinError(bin, sqrt(nentries)/n_generated_photons) ;
      }
    }
 
    // for the MCMC, our error is sort of the error on the integral calculation
    if( !useTracker ) {
      for( int nscat=0 ; nscat<=(int)nscatmax ; ++nscat ) {
        Int_t bin = hnscat.GetXaxis()->FindBin(nscat) ;
        double val = hnscat.GetBinContent(bin) ;
        double err = Perr[nscat]/Pscat[nscat]*val ;
        hnscat.SetBinError(bin,err) ;
      }
    }
    
    // scale the pulse profile histograms by the bin width
    hpulse.Scale(1.0/dt) ;
    for( int nscat=0 ; nscat<=(int)nscatmax ; ++nscat ) {
      hpulse_per_nscat[nscat]->Scale(1.0/dt) ;
    }
    
    // make cumulative pulse profile (so you can more easily see what
    // the total hit probability is)
    TH1D* hpulse_cum = (TH1D*) hpulse.Clone("hpulse_cum") ;
    hpulse_cum->GetYaxis()->SetTitle("cumulative hit probability") ;
    for( Int_t bin=1 ; bin<=hpulse_cum->GetNbinsX() ; ++bin ) {
      double val = hpulse_cum->GetBinContent(bin-1) ; 
      val += hpulse_cum->GetBinContent(bin) * hpulse_cum->GetXaxis()->GetBinWidth(bin) ;
      hpulse_cum->SetBinContent(bin,val) ;
    }
    
    f->Write() ;
 
    f->mkdir("EnsembleIntegratorHistograms")->cd() ;
    for( vector<JMarkovEnsembleIntegrator1D*>::iterator it=integrators.begin(); it!=integrators.end(); ++it ) {
      (*it)->writeHistograms() ;
    }
    f->cd() ;
 
    f->Close() ;
    //delete f ;
    cout << "Output in '" << fname << "'." << endl ;
  }
  cout << "Done!" << endl ;
  cout << endl ;
}

Variable Documentation

JBSversion

const double JBSversion = 3.21

Author

mjongen

Definition at line 20 of file JBeaconSimulator.cc.