Jpp  master_rocky-37-gf0c5bc59d
the software that should make you happy
Classes | Functions | Variables
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 JModulegetModule (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.

928  {
929  for( unsigned int i=0; i<detector.size(); ++i ) {
930  if( detector[i].getLocation() == location ) {
931  return detector[i] ;
932  }
933  }
934  cerr << "FATAL ERROR: no DOM found with location " << location << endl ;
935  exit(1) ;
936 }

◆ getMean()

double getMean ( vector< double > &  v)

get mean of vector content

Definition at line 56 of file JBeaconSimulator.cc.

56  {
57  double sum = 0 ;
58  for( vector<double>::iterator it=v.begin() ; it!=v.end() ; ++it )
59  sum += *it ;
60  return sum / v.size() ;
61 }
data_type v[N+1][M+1]
Definition: JPolint.hh:866

◆ getSigma()

double getSigma ( vector< double > &  v)

get standard deviation of vector content

Definition at line 64 of file JBeaconSimulator.cc.

64  {
65  double mean = getMean(v) ;
66  double varsum = 0 ;
67  for( vector<double>::iterator it=v.begin() ; it!=v.end() ; ++it )
68  varsum += (*it-mean)*(*it-mean) ;
69  return sqrt( varsum / v.size() ) ;
70 }
double getMean(vector< double > &v)
get mean of vector content

◆ scale()

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

scale vector content

Definition at line 73 of file JBeaconSimulator.cc.

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

◆ printResult()

void printResult ( vector< double > &  v)

Definition at line 78 of file JBeaconSimulator.cc.

78  {
79  double mean = getMean(v) ;
80  double sigma = getSigma(v) ;
81  cout << "Value = " << mean << " +- " << sigma/sqrt(1.0*v.size()) << endl ;
82 }
double getSigma(vector< double > &v)
get standard deviation of vector content
const double sigma[]
Definition: JQuadrature.cc:74

◆ writePaths()

void writePaths ( string  ofname,
vector< JPhotonPath paths 
)

Definition at line 85 of file JBeaconSimulator.cc.

85  {
87  writer.open(ofname.c_str()) ;
88  for( vector<JPhotonPath>::iterator it=paths.begin() ; it!=paths.end() ; ++it ) {
89  writer.put( *it ) ;
90  }
91  writer.close() ;
92  cout << "paths written to '" << ofname << "'." << endl ;
93 }
virtual void open(const char *file_name) override
Open file.
virtual void close()=0
Close device.
virtual bool put(const T &object)=0
Object output.

◆ 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.

99  {
100  double r = 0 ;
101  for( int pmt1=0 ; pmt1<31 ; ++pmt1 ) {
102  for( int pmt2=pmt1+1; pmt2<31 ; ++pmt2 ) {
103  double num = (dom.getPMT(pmt1).getPosition()-dom.getPMT(pmt2).getPosition()).getLength() ;
104  double den = (JVector3D(dom.getPMT(pmt1).getDirection())-JVector3D(dom.getPMT(pmt2).getDirection())).getLength() ;
105  r += num/den ;
106  }
107  }
108  r /= 0.5*31*30 ;
109  JPosition3D center(0,0,0) ;
110  for( int pmt=0 ; pmt<31 ; ++pmt ) {
111  center.add( dom.getPMT(pmt).getPosition() - r*JVector3D(dom.getPMT(pmt).getDirection()) ) ;
112  }
113  center.div(31) ;
114  return center ;
115 }
const JPMT & getPMT(const int index) const
Get PMT.
Definition: JModule.hh:172
const JDirection3D & getDirection() const
Get direction.
Data structure for position in three dimensions.
Definition: JPosition3D.hh:38
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:130
Data structure for vector in three dimensions.
Definition: JVector3D.hh:36
data_type r[M+1]
Definition: JPolint.hh:868

◆ 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.

128  {
129  JPosition3D dom_center = getDOMcenter(dom) ;
130 
131  // PMT 22 defines the negative z-direction of the DOM coordinate system
132  const JVersor3D neg_zdir( dom.getPMT(22).getDirection() ) ;
133  const JVersor3D zdir( -neg_zdir.getDX(), -neg_zdir.getDY(), -neg_zdir.getDZ() ) ;
134  // PMT 0 points in the negative y-direction of the DOM coordinate system (+some z-component)
135  JVector3D tmp ;
136  const JVersor3D xdir( tmp.cross( zdir, dom.getPMT(0).getDirection() ) ) ;
137  const JVersor3D ydir( tmp.cross( zdir, xdir ) ) ;
138 
139  // direction of beacon's position w.r.t. the DOM center
140  JVector3D _bdir( cos(beacon_theta)*JVector3D(zdir)
141  + sin(beacon_theta)*cos(beacon_phi)*JVector3D(xdir)
142  + sin(beacon_theta)*sin(beacon_phi)*JVector3D(ydir)
143  ) ;
144  JVersor3D bdir(_bdir) ; // start out in DOM z-direction
145 
146  JPosition3D bpos( dom_center + r*JVector3D(bdir) ) ;
147  JAxis3D ax(bpos,zdir) ;
148 
149  return ax ;
150 }
JPosition3D getDOMcenter(const JModule &dom)
Returns a good guess of where the actual center of the DOM is based on the PMT positions and directio...
Axis object.
Definition: JAxis3D.hh:41
Data structure for normalised vector in three dimensions.
Definition: JVersor3D.hh:28

◆ main()

int main ( int  argc,
char **  argv 
)

Definition at line 389 of file JBeaconSimulator.cc.

389  {
390  cout << "JBeaconSimulator version " << JBSversion << endl
391  << endl ;
392 
393  string ofname ;
394  string detectorFile ;
395 
396  JModuleLocation source_location ;
397  vector<int> target_pmt ;
398  bool useTracker ;
399  unsigned int nscatmax ;
400  unsigned long int nphotons ;
401  double Labs ;
402  double opening_angle_deg ;
403  double beacon_theta_deg ;
404  double beacon_phi_deg ;
405  double lHG ;
406  double gHG ;
407  double lR ;
408  double aR ;
409  vector<double> stepsizes ;
410  double target_step_size_deg ;
411  double tangential_step_size_deg ;
412  const vector<double> empty_vector ;
413  int nsamples ;
414  bool write_paths ;
415  bool no_remapping ;
416 
417  try {
418  JParser<string> zap;
419 
420  zap["o"] = make_field(ofname) ; // output file name
421  zap["a"] = make_field(detectorFile) ;
422  zap["s"] = make_field(source_location) ;
423  zap["t"] = make_field(target_pmt) ;
424  zap["m"] = make_field(nscatmax) = 5 ;
425  zap["N"] = make_field(nphotons) ;
426  zap["N-integration-samples"] = make_field(nsamples) = 1e6 ;
427  zap["tracker"] = make_field(useTracker) ;
428  zap["write-paths"] = make_field(write_paths) ;
429  zap["Labs"] = make_field(Labs) = 50 ;
430  zap["opening-angle"] = make_field(opening_angle_deg) = 45 ;
431  zap["beacon-theta-deg"] = make_field(beacon_theta_deg) = 56.2893-15 ; // 15 degrees above theta of upper row PMTs
432  zap["beacon-phi-deg"] = make_field(beacon_phi_deg) = -59.9998 ; // phi of PMT 3
433  zap["LscatHG"] = make_field(lHG) = 60.241 ;
434  zap["gHG"] = make_field(gHG) = 0.924 ;
435  zap["LscatR"] = make_field(lR) = 294.118 ;
436  zap["aR"] = make_field(aR) = 0.853 ;
437  zap["stepsizes"] = make_field(stepsizes) = empty_vector ;
438  zap["target-step-size-deg"] = make_field(target_step_size_deg) = 1.0 ;
439  zap["tangential-step-size-deg"] = make_field(tangential_step_size_deg) = 1.0 ;
440  zap["no-remapping"] = make_field(no_remapping) ;
441 
442  if (zap.read(argc, argv) != 0) exit(1) ;
443  }
444  catch(const exception &error) {
445  cerr << error.what() << endl ;
446  exit(1) ;
447  }
448 
449  // remove extension from output file name
450  {
451  vector<string> extensions ;
452  extensions.push_back(".root") ;
453  extensions.push_back(".paths") ;
454  for( vector<string>::iterator it=extensions.begin(); it!=extensions.end(); ++it ) {
455  size_t pos = ofname.find(*it) ;
456  if( pos != string::npos ) {
457  ofname = ofname.substr(0,pos) ;
458  }
459  }
460  }
461 
462  // use default step sizes when no new ones (or too few are specified)
463  if( !useTracker ) {
464  const int ndefvals = 6 ;
465  double defvals[ndefvals] = { 8, 5, 4, 3, 2, 1 } ;
466  while( stepsizes.size()<nscatmax ) {
467  double val = defvals[min(ndefvals-1,(int)stepsizes.size())] ;
468  stepsizes.push_back(val) ;
469  }
470  }
471 
472  // load detector
473  JDetector detector;
474  try {
475  load(detectorFile, detector);
476  }
477  catch(const JException& error) {
478  cout << "FATAL ERROR: could not open detector file '" << detectorFile << "'." ;
479  exit(1) ;
480  }
481  JModuleRouter moduleRouter(detector) ;
482 
483 
484  // initialize gRandom
485  gRandom = new TRandom3(0) ;
486 
487  // radius of a DOM (used for source and target model)
488  const double DOM_radius = 0.2159 ; // 0.2159 [m] (=17 inch diameter glass sphere)
489 
490 
491  // find source position and direction (stored as a JAxis3D)
492  // the source is a NB, whose position is deduced from the PMT positions in the detector file
493  // note that the distance from the DOM center is DOM_radius + 0.01 mm, to make sure that the
494  // initial vertex does not intersect the DOM (which is a problem when the source is located
495  // on the target DOM)
496  JModule source_module( getModule(detector,source_location) ) ;
497  JAxis3D source_axis( getBeaconAxis(source_module,DOM_radius+0.00001,beacon_theta_deg*M_PI/180,beacon_phi_deg*M_PI/180) ) ;
498  JVersor3D surface_normal( source_axis.getPosition() - getDOMcenter(source_module) ) ;
499 
500  cout << "SOURCE (NANOBEACON)" << endl
501  << "Positioned on string " << source_location.getString() << " floor " << source_location.getFloor()
502  << " (DOM center @ " << getDOMcenter(source_module) << ")" << endl
503  << "Position in local DOM coordinates: theta = "
504  << beacon_theta_deg << " deg, "
505  << "phi = " << beacon_phi_deg << " deg" << endl
506  << "distance w.r.t. DOM center = " << (source_axis.getPosition()-getDOMcenter(source_module)).getLength() << " m" << endl
507  << "unit vector from DOM center to beacon position = " << JVersor3D(source_axis.getPosition()-getDOMcenter(source_module)) << endl
508  << "beam direction = " << source_axis.getDirection() << endl
509  << "beam opening angle = " << opening_angle_deg << " deg" << endl
510  << "surface normal = " << surface_normal << endl
511  << endl ;
512 
513  // find target position and direction (stored as a JAxis3D)
514  if( target_pmt.size() != 3 ) {
515  cerr << "FATAL ERROR: provide a target PMT in the format "
516  << "\"string floor daq_channel\"." << endl ;
517  exit(1) ;
518  }
519 
520  // require valid PMT index for the MCMC
521  if( !useTracker ) {
522  if( target_pmt[2]<0 || target_pmt[2]>30 ) {
523  cerr << "FATAL ERROR: invalid target pmt index ("
524  << target_pmt[2] << ")" << endl ;
525  exit(2) ;
526  }
527  }
528 
529  JModuleLocation target_location(target_pmt[0],target_pmt[1]) ;
530  JModule target_module( getModule(detector,target_location) ) ;
531  JAxis3D target_axis( getDOMcenter(target_module), target_module.getPMT(target_pmt[2]).getDirection() ) ;
532 
533  // source model
534  JDirectedSource* src = new JDirectedSource( source_axis.getDirection(),
535  opening_angle_deg*M_PI/180,
536  true,
537  surface_normal ) ;
538  src->setPosition( source_axis.getPosition() ) ;
539 
540  // initialize scattering model
541  cout << "PHYSICS" << endl ;
542  JHenyeyGreensteinScattering sm1( lHG, gHG ) ;
543  cout << "Henyey-Greenstein scattering length = " << lHG << " [m], g = " << gHG << endl ;
544  JRayleighScattering sm2( lR, aR ) ;
545  cout << "Rayleigh scattering length = " << lR << " [m], a = " << aR << endl ;
546  JScatteringModel* sm = new JCombinedScattering( &sm1, &sm2 ) ;
547  cout << "Absorption length = " << Labs << " [m]" << endl ;
548 
549  // water properties
550  // for now this is not customizable
551  const double water_depth = 3000 ; // m
552  const double pressure = 0.1*water_depth ; // atm
553  JDispersion jd(pressure) ;
554  const double nbwavelength = 470 ; // nm
555  const double index_of_refraction = jd.getIndexOfRefractionGroup(nbwavelength) ;
556  const double cw = JTOOLS::getSpeedOfLight() / index_of_refraction ;
557  cout << "Assuming speed of light in water = " << cw << " [m/ns]" << endl ;
558  cout << endl ;
559 
560 
561  // storage container for generated photon paths
562  vector<JPhotonPath> paths ;
563 
564  // target model (for now this is not customizable)
565  // we use the approximation that a PMT is a spherical cap on a DOM, with a uniform
566  // acceptance, i.e. the photon always causes a hit when it impacts the cap
567  cout << "TARGET" << endl ;
568  const double pmt_opening_angle_deg = 30 ;
569  const double pmt_opening_angle_rad = pmt_opening_angle_deg*M_PI/180.0 ;
570  const double pmt_ctmin = cos(0.5*pmt_opening_angle_rad) ;
571  cout << "DOM: string " << target_pmt[0]
572  << ", floor " << target_pmt[1]
573  << ", @" << getDOMcenter(target_module) << endl
574  << "Target radius = " << DOM_radius << " [m]" << endl
575  << "PMT: " ;
576  if( useTracker ) cout << "all 31 PMTs" ;
577  else cout << "PMT " << target_pmt[2] << ", direction = " << target_axis.getDirection() ;
578  cout << endl
579  << "PMTs are simulated as spherical caps" << endl
580  << "PMT opening angle = " << pmt_opening_angle_deg << " deg." << endl
581  << endl ;
582 
583  JTargetModel* trg ;
584  if( useTracker ) {
585  // for the tracker we use an isotropic sphere that represents an entire DOM
586  // later on, we will divide the hits over the different PMTs
587  trg = new JIsotropicTarget() ;
588  trg->setRadius(DOM_radius) ;
589  trg->setPosition(target_axis.getPosition()) ;
590  } else {
591  // for the MCMC we use a single PMT target
592  trg = new JPMTTarget(target_axis,pmt_opening_angle_rad,DOM_radius) ;
593  }
594 
595  unsigned long int n_generated_photons = nphotons ;
596  if( useTracker ) { // generate paths by photon tracking
597  //n_generated_photons = 1000e6 ;
598  cout << "Generating photon paths by photon tracking" << endl ;
600  paths = mpt.trackPhotons(
601  n_generated_photons,
602  src, sm, trg,
603  Labs ) ;
604  cout << "Photon hits: " << paths.size() << " / " << n_generated_photons << endl ;
605  cout << endl ;
606  }
607 
608  vector<double> Pscat(nscatmax+1,1) ;
609  vector<double> Perr(nscatmax+1,1) ;
611  if( !useTracker ) { // generate paths by MCMC
613  if( no_remapping ) {
614  // deactivate coordinate remapping.
615  // That will most likely lead to trouble with the 1/r^2 singularities
616  // in the path probability density, meaning the space will be undersampled
617  // for cases where vertices are very close together.
618  mpg.setCoordinateRemapping(false) ;
619  cout << "WARNING: coordinate remapping is deactivated!" << endl ;
620  }
621  mpg.setTargetStepSize_deg(target_step_size_deg) ;
622  mpg.setTangentialStepSize_deg(tangential_step_size_deg) ;
623  int nsteps_save = 250 ;
624  int nsteps_burn_in = 10*nsteps_save ;
625  for( int nscat=0 ; nscat<=(int)nscatmax ; ++nscat ) {
626  cout << "nscat = " << nscat << endl ;
627 
628  double stepsize = 1 ; // default step size
629  if( nscat>0 && (int)stepsizes.size()>=nscat ) stepsize = stepsizes[nscat-1] ;
630  mpg.setRadialStepSize_m(stepsize) ;
631 
632  if( nscat != 0 ) cout << "step size = " << stepsize << " [m]" << endl ;
633 
634  // create start_path
635  FreePathSolver fps ;
636  JPhotonPath start_path = fps.solve(nscat,src,(JPMTTarget*)trg) ;
637  cout << "path = " << endl ;
638  cout << start_path.getString() << endl ;
639  double start_rho = getPhotonPathProbabilityDensity(start_path,src,sm,trg,Labs) ;
640  cout << "start_rho = " << start_rho << endl ;
641  if( start_rho == 0 ) {
642  cout << "No possible photon path found with nscat = " << nscat <<". Skipping nscat = " << nscat << "." << endl
643  << "WARNING: note that this does not necessarily mean there are no such paths, simply that the search algorithm could not find one." << endl
644  << endl;
645  continue ;
646  }
647 
648  /*cout << "Found the following start path: " << endl ;
649  for( vector<JPosition3D>::iterator it=start_path.begin() ; it!= start_path.end() ; ++it ) {
650  cout << *it << endl ;
651  }
652  cout << endl ;*/
653 
654  vector<JPhotonPath> these_paths = mpg.generateEnsemble(
655  n_generated_photons,
656  start_path,
657  src, sm, trg,
658  Labs,
659  nsteps_burn_in, nsteps_save ) ;
660  int nsteps = mpg.getNsteps() ;
661  int naccepted = mpg.getNacceptedSteps() ;
662  if( nscat != 0 ) {
663  cout << "Accepted " << naccepted << " / " << nsteps << " steps (" << mpg.getFractionAccepted()*100.0 << "%)" << endl
664  << "Accepted " << 100.0*mpg.getFractionAccepted_radial() << "% of the radial steps." << endl
665  << "Accepted " << 100.0*mpg.getFractionAccepted_tangential() << "% of the tangential steps." << endl ;
666  cout << "NOTE: as a rule of thumb, the optimal value is around ~23% (this can be achieved by modifying the various step sizes)" << endl ;
667  }
668  paths.insert(paths.end(), these_paths.begin(), these_paths.end()) ;
669  cout << endl ;
670 
671  if( these_paths.size() == 0 ) {
672  cerr << "Something went wrong. The MCMC generator returned no paths." << endl ;
673  Pscat[nscat] = 0.0 ;
674  Perr[nscat] = 0.0 ;
675  } else {
676  cout << "Creating JMarkovEnsembleIntegrator" << endl ;
677  JMarkovEnsembleIntegrator1D* jmi = new JMarkovEnsembleIntegrator1D(these_paths,trg,100,100,200) ;
678 
679  cout << "Integrating" << endl ;
680  vector<double> result = jmi->integrate(nsamples,nscat,src,sm,trg,Labs) ;
682  Pscat[nscat] = getMean(result) ;
683  Perr[nscat] = getSigma(result) / sqrt(nsamples) ;
684  //delete jmi ;
685  integrators.push_back(jmi) ;
686  }
687  cout << endl ;
688  }
689  }
690 
691  if( paths.size()==0 ) {
692  cerr << "FATAL ERROR: no paths were generated!" << endl ;
693  exit(1) ;
694  }
695 
696 
697  // get path lengths and directions into an array
698  // also find what the maximal number of scatterings is
699  unsigned int npaths = paths.size() ;
700  vector<double> path_lengths(npaths) ;
701  vector<double> xs(npaths) ;
702  vector<double> ys(npaths) ;
703  vector<double> zs(npaths) ;
704  for( unsigned int i=0; i<npaths; ++i ) {
705  path_lengths[i] = paths[i].getLength() ;
706  xs[i] = paths[i].back().getX() ;
707  ys[i] = paths[i].back().getY() ;
708  zs[i] = paths[i].back().getZ() ;
709  int nscat = paths[i].size()-2 ;
710  if( nscat>(int)nscatmax ) nscatmax = nscat ;
711  }
712 
713  // find min and max x, y and z
714  double xmin = trg->getPosition().getX() - 1.1*DOM_radius ; // *min_element(xs.begin(),xs.end() ) ;
715  double xmax = trg->getPosition().getX() + 1.1*DOM_radius ; //*max_element(xs.begin(),xs.end() ) ;
716  double ymin = trg->getPosition().getY() - 1.1*DOM_radius ; //*min_element(ys.begin(),ys.end() ) ;
717  double ymax = trg->getPosition().getY() + 1.1*DOM_radius ; //*max_element(ys.begin(),ys.end() ) ;
718  double zmin = trg->getPosition().getZ() - 1.1*DOM_radius ; //*min_element(zs.begin(),zs.end() ) ;
719  double zmax = trg->getPosition().getZ() + 1.1*DOM_radius ; //*max_element(zs.begin(),zs.end() ) ;
720 
721  // find out which PMT is hit
722  vector<int> PMT_is_hit ;
723 
724  if( !useTracker ) {
725  // in the MCMC we know which PMT is hit
726  PMT_is_hit = vector<int>(npaths,target_pmt[2]) ;
727  }
728 
729  if( useTracker ) {
730  PMT_is_hit = vector<int>(npaths,-1) ;
731  vector<unsigned int> nHitsPerPMT(31,0) ;
732  // in the tracker we have to figure out for every hit whether the PMT was hit or not
733  for( unsigned int i=0; i<npaths; ++i ) {
734 
735  // find out which PMTs were hit (should be only one)
736  vector<int> hitPMTs ;
737  JVersor3D hit_direction = JVersor3D( paths[i].back() - target_axis.getPosition() ) ;
738  for( unsigned int pmt=0 ; pmt<31 ; ++pmt ) {
739  double ct = target_module.getPMT(pmt).getDirection().getDot(hit_direction) ;
740  if( ct>pmt_ctmin ) hitPMTs.push_back(pmt) ;
741  }
742  if( hitPMTs.size()==1 ) {
743  PMT_is_hit[i] = hitPMTs[0] ;
744  ++nHitsPerPMT[ hitPMTs[0] ] ;
745  }
746  if( hitPMTs.size()>1 ) {
747  cerr << "FATAL ERROR: hit " << hitPMTs.size() << " PMTs simultaneously" << endl ;
748  exit(1) ;
749  }
750  }
751  // print number of hits separated by PMT
752  for( int pmt=0 ; pmt<31 ; ++pmt ) {
753  if( nHitsPerPMT[pmt]>0 ) {
754  cout << "PMT" << setw(2) << pmt << ": " << setw(15) << nHitsPerPMT[pmt] << " hits" << endl ;
755  }
756  }
757  cout << endl ;
758  }
759 
760  // range and binning for pulse profile(s)
761  const double dt = 0.1 ; // bin width in ns
762  const unsigned int margin = 5 ; // ns
763  double Lmin = *min_element(path_lengths.begin(),path_lengths.end() ) ;
764  double Lmax = *max_element(path_lengths.begin(),path_lengths.end() ) ;
765  double tmin = floor( Lmin / cw - margin ) ;
766  double tmax = ceil( Lmax / cw + margin ) ;
767  int nbinst = (int) round( (tmax-tmin)/dt ) ;
768 
769  int pmtfront = 0 ;
770  int pmtback = 31 ;
771  if( !useTracker ) {
772  // in the MCMC there is only one PMT
773  pmtfront = target_pmt[2] ;
774  pmtback = target_pmt[2]+1 ;
775  }
776 
777  for( int pmt=pmtfront ; pmt<pmtback ; ++pmt ) {
778  // open file
779  char fname[500] ;
780  sprintf( fname, "%s_PMT%i.root", ofname.c_str(), pmt ) ;
781  TFile* f = new TFile(fname,"recreate") ;
782 
783  // create histograms
784  TH1D hpulse("hpulse","Instantaneous flash pulse profile;time since flash [ns];hit probability / ns",
785  nbinst,tmin,tmax) ;
786  vector<TH1D*> hpulse_per_nscat(nscatmax+1) ;
787  for( int nscat=0 ; nscat<=(int)nscatmax ; ++nscat ) {
788  char hname[200] ;
789  sprintf(hname,"hpulse_nscat%i",nscat) ;
790  hpulse_per_nscat[nscat] = new TH1D(hname,"Instantaneous flash pulse profile;time since flash [ns];hit probability / ns",
791  nbinst,tmin,tmax ) ;
792  }
793  TH2D hXY("hXY","hit position;X;Y",
794  100,xmin,xmax,
795  100,ymin,ymax ) ;
796  hXY.SetOption("colz") ;
797  TH2D hXZ("hXZ","hit position;X;Z",
798  100,xmin,xmax,
799  100,zmin,zmax ) ;
800  hXZ.SetOption("colz") ;
801  TH2D hYZ("hYZ","hit position;Y;Z",
802  100,ymin,ymax,
803  100,zmin,zmax ) ;
804  hYZ.SetOption("colz") ;
805  TH1D hnscat("hnscat",";N_{scat};probability",nscatmax+2,-0.5,nscatmax+1.5) ;
806 
807  TH3D hregion("hregion","Photon paths;X;Y;Z",
808  100,
809  trg->getPosition().getX() - 2*DOM_radius,
810  trg->getPosition().getX() + 2*DOM_radius,
811  100,
812  trg->getPosition().getY() - 2*DOM_radius,
813  trg->getPosition().getY() + 2*DOM_radius,
814  100,
815  trg->getPosition().getZ() - 2*DOM_radius,
816  trg->getPosition().getZ() + 2*DOM_radius) ;
817 
818  // fill histograms
819  vector<JPhotonPath> paths_to_write ;
820  for( unsigned int i=0 ; i<npaths ; ++i ) {
821  if( PMT_is_hit[i] != pmt ) continue ;
822 
823  paths_to_write.push_back(paths[i]) ;
824 
825  const double path_length = paths[i].getLength() ;
826  double t = path_length/cw ;
827  int nscat = paths[i].size()-2 ;
828 
829  double w = 1 ;
830  if( useTracker ) w = 1.0 / n_generated_photons ;
831  else w = Pscat[nscat] / n_generated_photons ;
832 
833  hpulse.Fill(t,w) ;
834  hnscat.Fill(nscat,w) ;
835  hpulse_per_nscat[nscat]->Fill(t,w) ;
836 
837  hXY.Fill(xs[i],ys[i],w) ;
838  hXZ.Fill(xs[i],zs[i],w) ;
839  hYZ.Fill(ys[i],zs[i],w) ;
840 
841  if( paths.size()<10000 ) {
842  const double dl = 0.001 ; // 1 mm
843  const double dw = w*dl/path_length ;
844  double l_along_seg = 0.5*dl ;
845  int nfilled = 0 ;
846  for( int nv=0; nv<(int)paths[i].size()-1; ++nv ) {
847  double seg_length = paths[i][nv+1].getDistance(paths[i][nv]) ;
848  JVersor3D seg_dir( paths[i][nv+1] - paths[i][nv] ) ;
849  while( l_along_seg<seg_length ) {
850  JPosition3D pos( paths[i][nv] + l_along_seg*JVector3D(seg_dir) ) ;
851  hregion.Fill( pos.getX(), pos.getY(), pos.getZ(), dw ) ;
852  if( pos.getDistance( trg->getPosition() ) < trg->getRadius() ) {
853  cerr << "Woohoooo" << endl ;
854  exit(1) ;
855  }
856  ++nfilled ;
857  l_along_seg += dl ;
858  }
859  l_along_seg -= seg_length ;
860  }
861  }
862  }
863  // scale hregion histogram, so that SetShowProjection3D will no show an
864  // empty histogram
865  hregion.Scale(1e20) ;
866 
867  // write paths
868  if( write_paths ) {
869  char path_fname[500] ;
870  sprintf( path_fname, "%s_PMT%i.paths", ofname.c_str(), pmt ) ;
871  cout << "Writing paths." << endl ;
872  writePaths( path_fname, paths_to_write ) ;
873  cout << endl ;
874  }
875 
876  // for the tracker, our error is just the Poisson error
877  if( useTracker ) {
878  for( Int_t bin=1 ; bin<=hnscat.GetNbinsX() ; ++bin ) {
879  double val = hnscat.GetBinContent(bin) ;
880  double nentries = max(1.0, val*n_generated_photons ) ;
881  hnscat.SetBinError(bin, sqrt(nentries)/n_generated_photons) ;
882  }
883  }
884 
885  // for the MCMC, our error is sort of the error on the integral calculation
886  if( !useTracker ) {
887  for( int nscat=0 ; nscat<=(int)nscatmax ; ++nscat ) {
888  Int_t bin = hnscat.GetXaxis()->FindBin(nscat) ;
889  double val = hnscat.GetBinContent(bin) ;
890  double err = Perr[nscat]/Pscat[nscat]*val ;
891  hnscat.SetBinError(bin,err) ;
892  }
893  }
894 
895  // scale the pulse profile histograms by the bin width
896  hpulse.Scale(1.0/dt) ;
897  for( int nscat=0 ; nscat<=(int)nscatmax ; ++nscat ) {
898  hpulse_per_nscat[nscat]->Scale(1.0/dt) ;
899  }
900 
901  // make cumulative pulse profile (so you can more easily see what
902  // the total hit probability is)
903  TH1D* hpulse_cum = (TH1D*) hpulse.Clone("hpulse_cum") ;
904  hpulse_cum->GetYaxis()->SetTitle("cumulative hit probability") ;
905  for( Int_t bin=1 ; bin<=hpulse_cum->GetNbinsX() ; ++bin ) {
906  double val = hpulse_cum->GetBinContent(bin-1) ;
907  val += hpulse_cum->GetBinContent(bin) * hpulse_cum->GetXaxis()->GetBinWidth(bin) ;
908  hpulse_cum->SetBinContent(bin,val) ;
909  }
910 
911  f->Write() ;
912 
913  f->mkdir("EnsembleIntegratorHistograms")->cd() ;
914  for( vector<JMarkovEnsembleIntegrator1D*>::iterator it=integrators.begin(); it!=integrators.end(); ++it ) {
915  (*it)->writeHistograms() ;
916  }
917  f->cd() ;
918 
919  f->Close() ;
920  //delete f ;
921  cout << "Output in '" << fname << "'." << endl ;
922  }
923  cout << "Done!" << endl ;
924  cout << endl ;
925 }
JAxis3D getBeaconAxis(const JModule &dom, double r, double beacon_theta, double beacon_phi)
Get the position and orientation of the nanobeacon on a DOM.
void printResult(vector< double > &v)
const JModule & getModule(const JDetector &detector, const JModuleLocation &location)
find module with a given string and floor number
const double JBSversion
void writePaths(string ofname, vector< JPhotonPath > paths)
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2142
Finds photon paths from a nanobeacon source to a target PMT that have a non-zero probability.
JPhotonPath solve(int nscat, const JDirectedSource *src, const JPMTTarget *trg)
Detector data structure.
Definition: JDetector.hh:96
Router for direct addressing of module data in detector data structure.
Data structure for a composite optical module.
Definition: JModule.hh:75
double getY() const
Get y position.
Definition: JVector3D.hh:104
double getZ() const
Get z position.
Definition: JVector3D.hh:115
double getX() const
Get x position.
Definition: JVector3D.hh:94
General exception.
Definition: JException.hh:24
Implementation of the JScatteringModel interface with that combines two scattering models into one ef...
Implementation of the JSourceModel class that represents a directed source with a flat intensity dist...
Implementation of the JScatteringModel interface with scattering according to the Henyey-Greenstein f...
Implementation of the JTargetModel class that represents a spherically symmetric target.
Implementation of JMarkovEnsembleIntegrator interface with 1D histograms.
vector< double > integrate(int N, int nscat, JSourceModel *src, JScatteringModel *sm, JTargetModel *trg, double lambda_abs)
Integrate with N samples.
The JMarkovPathGenerator generates ensembles of photon paths using a Markov Chain Monte Carlo (MCMC).
int getNsteps()
get the number of Markov steps taken in the last call to generateEnsemble (after burn-in)
void setRadialStepSize_m(double val)
set the average step size in [m] in the radial direction for the scattering vertices
void setCoordinateRemapping(bool val=true)
activate or deactive coordinate remapping
void setTargetStepSize_deg(double val)
set the average step size in degrees for the impact point on the target
void setTangentialStepSize_deg(double val)
Set the average step size theta in degrees for steps in the tangential direction for the scattering v...
double getFractionAccepted_radial()
get the fraction of steps that were accepted when a radial step was performed during the last call to...
int getNacceptedSteps()
get the number of accepted steps taken during the last call to generateEnsemble (after burn-in)
std::vector< JPhotonPath > generateEnsemble(int n, const JPhotonPath &start_path, JSourceModel *src, JScatteringModel *sm, JTargetModel *trg, double lambda_abs, int nsteps_burn_in, int nsteps_save)
Generate an ensemble of n paths with a fixed number of scatterings by MCMC-sampling the given scatter...
double getFractionAccepted_tangential()
get the fraction of steps that were accepted when a tangential step was performed during the last cal...
double getFractionAccepted()
get the fraction of accepted steps during the last call to generateEnsemble (after burn-in)
The JMarkovPhotonTracker generates ensembles of photon paths by tracking photons from a source to a t...
std::vector< JPhotonPath > trackPhotons(unsigned long int n, JSourceModel *src, JScatteringModel *sm, JTargetModel *trg, double lambda_abs)
Track n photons.
Implementation of the JTargetModel class that represents a single PMT on a DOM.
A photon path.
Definition: JPhotonPath.hh:38
string getString()
get a string with the path vertex positions
Definition: JPhotonPath.hh:136
Implementation of the JScatteringModel interface with Rayleigh scattering.
Virtual base class for a scattering model.
void setPosition(JPosition3D &_pos)
Virtual base class for a light detector ("photon target").
const JPosition3D & getPosition() const
void setPosition(JPosition3D &_pos)
void setRadius(double _r)
Utility class to parse command line options.
Definition: JParser.hh:1698
int read(const int argc, const char *const argv[])
Parse the program's command line options.
Definition: JParser.hh:1992
Implementation of dispersion for water in deep sea.
Definition: JDispersion.hh:28
const double xmax
Definition: JQuadrature.cc:24
const double xmin
Definition: JQuadrature.cc:23
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
double getPhotonPathProbabilityDensity(JPhotonPath &p, JSourceModel *src, JScatteringModel *sm, JTargetModel *trg, double lambda_abs)
Return the probability density for a photon path with the given ingredients.
const double getSpeedOfLight()
Get speed of light.
data_type w[N+1][M+1]
Definition: JPolint.hh:867

Variable Documentation

◆ JBSversion

const double JBSversion = 3.21
Author
mjongen

Definition at line 20 of file JBeaconSimulator.cc.