Jpp  19.1.0-rc.1
the software that should make you happy
JBeaconSimulator.cc
Go to the documentation of this file.
1 /**
2  * \author mjongen
3  */
4 
5 /*******************************************************************************
6  JBeaconSimulator
7 
8  Note that this program is still under development.
9 
10  It uses the JMarkov library to simulate the response of a given PMT to an
11  instantaneous flash of a nanobeacon.
12 
13  -s "source", a string like "2 1", indicating the string and floor of the
14  DOM on which the beacon is activated
15 
16  -t "target", a string like "1 10 3", indicating the string, floor and DAQ
17  channel of the target PMT
18 
19 *******************************************************************************/
20 const double JBSversion = 3.21 ;
21 
22 // c++ standard library
23 #include <string>
24 
25 // JPP
26 #include "Jeep/JParser.hh"
29 #include "JDetector/JModuleLocation.hh"
30 #include "JGeometry3D/JAxis3D.hh"
36 #include "JPhysics/JPDFToolkit.hh"
37 #include "JPhysics/JDispersion.hh"
38 #include "JPhysics/KM3NeT.hh"
39 
40 // root
41 #include "TH1D.h"
42 #include "TRandom3.h"
43 #include "TFile.h"
44 // namespaces
45 using namespace std ;
46 using namespace JLANG ;
47 using namespace JDETECTOR ;
48 using namespace JMARKOV ;
49 using namespace JPHYSICS ;
50 using namespace KM3NET ;
51 
52 /// find module with a given string and floor number
53 const JModule& getModule( const JDetector& detector, const JModuleLocation& location ) ;
54 
55 /// get mean of vector content
56 double getMean( vector<double>& v ) {
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 }
62 
63 /// get standard deviation of vector content
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 }
71 
72 /// scale vector content
73 void scale( vector<double>& v, double c ) {
74  for( vector<double>::iterator it=v.begin() ; it!=v.end() ; ++it )
75  *it = c * (*it) ;
76 }
77 
79  double mean = getMean(v) ;
80  double sigma = getSigma(v) ;
81  cout << "Value = " << mean << " +- " << sigma/sqrt(1.0*v.size()) << endl ;
82 }
83 
84 // write a vector of generated photon paths to a file
85 void writePaths( string ofname, vector<JPhotonPath> paths ) {
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 }
94 
95 /**
96  Returns a good guess of where the actual center of the DOM is based on the PMT positions
97  and directions.
98  **/
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 }
116 
117 
118 /**
119  Get the position and orientation of the nanobeacon on a DOM.
120 
121  The local coordinate system of the DOM is deduced from the PMT
122  positions.
123 
124  r = DOM radius [m]
125  beacon_theta = theta angle of the beacon in the DOM's coordinate system [rad]
126  phi_theta = phi angle of the beacon in the DOM's coordinate system [rad]
127  **/
128 JAxis3D getBeaconAxis( const JModule& dom, double r, double beacon_theta, double beacon_phi ) {
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 }
151 
152 /**
153  Finds photon paths from a nanobeacon source to a target PMT that
154  have a non-zero probability. This means that
155  - the emission angle has to be within the nanobeacon opening angle
156  - the photon has to hit the PMT at an angle less than 90 degrees
157  - the photon may not intersect the target
158 
159  It works by maximizing a score function. If we find a score of 0,
160  the path meets all criteria.
161  **/
163  public :
164 
166 
167  JPhotonPath solve(int nscat, const JDirectedSource* src, const JPMTTarget* trg ) {
168  //cout << "Running FreePathSolver::solve with nscat = " << nscat << endl ;
169 
170  JPhotonPath p(min(2,nscat)) ;
171 
172  // start at source
173  p.front() = src->getPosition() ;
174  // end on PMT surface
175  p.back() = trg->getPosition() ;
176  if( trg->getRadius()>0 ) p.back() = trg->getPosition() + trg->getRadius() * JVector3D(trg->getDirection()) ;
177 
178  if( nscat == 0 ) return p ;
179 
180  // create intermediate points
181  for( int nv=1 ; nv<(int)p.size()-1 ; ++nv ) {
182  double factor = 1.0 * nv / (p.size()-1) ;
183  p[nv] = (1-factor)*p.front() + factor*p.back() ;
184  }
185 
186  int maxnsteps = -1 ;
187 
188  if( nscat == 1 ) {
189  double score = getScore(p,src,trg) ;
190  double dl = 1 ; // m
191  const double dlmin = 0.0001 ;
192  JPosition3D dx(1,0,0) ;
193  JPosition3D dy(0,1,0) ;
194  JPosition3D dz(0,0,1) ;
195  int istep = 0 ;
196  while( score<0 ) {
197  if( ++istep == maxnsteps ) break ;
198  double initscore = score ;
199 
200  // try step in x-direction
201  p[1] = p[1] + dl*dx ;
202  if( getScore(p,src,trg) < score ) {
203  // try moving in the opposite direction
204  p[1] = p[1] - 2*dl*dx ;
205  if( getScore(p,src,trg) < score ) p[1] = p[1] + dl*dx ;
206  }
207 
208  // try step in y-direction
209  p[1] = p[1] + dl*dy ;
210  if( getScore(p,src,trg) < score ) {
211  // try moving in the opposite direction
212  p[1] = p[1] - 2*dl*dy ;
213  if( getScore(p,src,trg) < score ) p[1] = p[1] + dl*dy ;
214  }
215 
216  // try step in z-direction
217  p[1] = p[1] + dl*dz ;
218  if( getScore(p,src,trg) < score ) {
219  // try moving in the opposite direction
220  p[1] = p[1] - 2*dl*dz ;
221  if( getScore(p,src,trg) < score ) p[1] = p[1] + dl*dz ;
222  }
223 
224  // see if it worked
225  score = getScore(p,src,trg) ;
226  // check if we are stuck
227  if( score == initscore ) {
228  dl = 0.5 * dl ;
229  //cout << "dl = " << dl << endl ;
230  if( dl<dlmin ) break ;
231  }
232  if( p.getLength() > 10000 ) break ;
233  }
234  }
235 
236  if( nscat>1) {
237  //cout << "nscat more than one" << endl ;
238  //cout << "starting with " << endl ;
239  //for( vector<JPosition3D>::iterator it=p.begin() ; it!= p.end() ; ++it ) cout << *it << endl ;
240  //cout << endl ;
241 
242  double score = getScore(p,src,trg) ;
243  double dl = 1 ; // m
244  const double dlmin = 1e-10 ;
245 
246  const int nsm = 6 ;
247  JPosition3D single_moves[nsm] = { JPosition3D(1,0,0), JPosition3D(0,1,0), JPosition3D(0,0,1),
248  JPosition3D(-1,0,0), JPosition3D(0,-1,0), JPosition3D(0,0,-1) } ;
249  /*ctor< pair<JPosition3D,JPosition3D> > moves ;
250  for( int i=0 ; i<nsm ; ++i ) {
251  for( int j=0 ; j<nsm ; ++j ) {
252  moves.push_back( pair<JPosition3D,JPosition3D>(single_moves[i],single_moves[j]) ) ;
253  }
254  }*/
255 
256  int i = 0 ;
257 
258  while( score<0 ) {
259  if( ++i>100 ) break ;
260  //cout << "@ " << JPosition3D(p[1]-src->getPosition()) << ", " << JPosition3D(p[2]-src->getPosition()) << ", score = " << score << ", dl = " << dl << endl ;
261  double initscore = score ;
262 
263  for( int i=0; i<nsm; ++i ) {
264  // try move
265  p[1] = p[1] + dl * single_moves[i] ;
266  double newscore = getScore(p,src,trg) ;
267  if( newscore < score ) {
268  // if the score is reduced, undo the move
269  p[1] = p[1] - dl * single_moves[i] ;
270  } else {
271  //cout << "First vertex adds " << single_moves[i]
272  // << ", score difference = " << newscore - score << endl ;
273  score = newscore ;
274  }
275  }
276 
277  for( int i=0; i<nsm; ++i ) {
278  // try move
279  p[2] = p[2] + dl * single_moves[i] ;
280  double newscore = getScore(p,src,trg) ;
281  if( newscore < score ) {
282  // if the score is reduced, undo the move
283  p[2] = p[2] - dl * single_moves[i] ;
284  } else {
285  //cout << "Second vertex adds " << single_moves[i]
286  // << ", score difference = " << newscore - score << endl ;
287  score = newscore ;
288  }
289  }
290 
291  // see if it worked
292  score = getScore(p,src,trg) ;
293  // check if we are stuck
294  if( score == initscore ) {
295  //cout << "Stuck" << endl ;
296  dl = 0.5*dl ;
297  if( dl<dlmin ) break ;
298  }
299  if( p.getLength() > 10000 ) {
300  //cout << "Too long" << endl ;
301  break ;
302  }
303  }
304 
305  // fill in the gaps
306  JPhotonPath p2(nscat) ;
307  p2[0] = p[0] ;
308  p2[1] = p[1] ;
309  p2[2] = p[2] ;
310  p2.back() = p.back() ;
311  for( int nv=3 ; nv<nscat+1 ; ++nv ) {
312  double factor = 1.0 * (nv-2) / (nscat-1) ;
313  p2[nv] = (1-factor)*p2[2] + factor*p2.back() ;
314  }
315  return p2 ;
316  }
317 
318  return p ;
319 
320  }
321 
322  protected :
323 
324  double sourceScore( const JPhotonPath& p, const JDirectedSource* src ) {
325  double ct = JVersor3D(p[1]-src->getPosition()).getDot(src->getDirection()) ;
326  const double ctmin = cos(0.5*src->getOpeningAngle()) ;
327  double val = max(ctmin-ct,0.0) / (1-ctmin) ;
328  return -val*val ;
329  }
330 
331  double targetScore( const JPhotonPath& p, const JPMTTarget* trg ) {
332  double ct = JVersor3D(p[p.size()-2]-p[p.size()-1]).getDot(trg->getDirection()) ;
333  const double ctmin = cos(0.5*trg->getOpeningAngle()) ;
334  double val = max(ctmin-ct,0.0) / (1-ctmin) ;
335  return -val*val ;
336  }
337 
338  double intersectScore( const JPhotonPath& p, const JPMTTarget* trg ) {
339  double score = 0 ;
340 
341  // check whether the photon path goes through the target sphere
342  // we do not need to consider the last vertex, because if that intersects the sphere
343  // the target score would be 0.
344  if( p.n == 2 ) return score ;
345 
346  for( int nv=1 ; nv<(int)p.size()-2 ; ++nv ) {
347  JPhotonPath seg(0) ;
348  seg[0] = p[nv] ;
349  seg[1] = p[nv+1] ;
350  if( seg.hitsSphere(trg->getPosition(),trg->getRadius()) ) {
351  double d = JAxis3D( seg[0], JVersor3D(seg[1]-seg[0]) ).getDistance(trg->getPosition()) ;
352  double R = trg->getRadius() ;
353  double val = (R-d)/R ;
354  score -= val*val ;
355  //cout << "Segment from " << nv << " to " << nv+1 << "intersects sphere"
356  // << ", distance = " << d << ", val = " << val << endl ;
357  }
358  }
359 
360  /*JPhotonPath pshort(p) ;
361  pshort.resize( pshort.size()-1 ) ;
362  --pshort.n ;
363  if( pshort.hitsSphere(trg->getPosition(),trg->getRadius()) ) {
364  JVersor3D hit_dir( pshort.getSphereHitPosition(trg->getPosition(),trg->getRadius()) - trg->getPosition() ) ;
365  double val = hit_dir.getDot(trg->getDirection()) ;
366  score -= val * val ;
367  }
368 
369  // extra penalty when a vertex is inside the target
370  for( int nv=1 ; nv<p.n-1 ; ++nv ) {
371  double r = p[nv].getDistance(trg->getPosition()) ;
372  if( r < trg->getRadius() ) {
373  double val = ( trg->getRadius() - r ) / trg->getRadius() ;
374  score -= val * val ;
375  }
376  }*/
377 
378  return score ;
379  }
380 
381  // total score
382  double getScore( const JPhotonPath& p, const JDirectedSource* src, const JPMTTarget* trg ) {
383  double score = sourceScore(p,src) + targetScore(p,trg) + intersectScore(p,trg) ;
384  return score ;
385  }
386 
387 } ;
388 
389 int main( int argc, char** argv ) {
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 }
926 
927 
928 const JModule& getModule( const JDetector& detector, const JModuleLocation& location ) {
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 }
JAxis3D getBeaconAxis(const JModule &dom, double r, double beacon_theta, double beacon_phi)
Get the position and orientation of the nanobeacon on a DOM.
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...
int main(int argc, char **argv)
double getMean(vector< double > &v)
get mean of vector content
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
double getSigma(vector< double > &v)
get standard deviation of vector content
void scale(vector< double > &v, double c)
scale vector content
void writePaths(string ofname, vector< JPhotonPath > paths)
Direct access to module in detector data structure.
Auxiliary methods for PDF calculations.
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2158
Properties of KM3NeT PMT and deep-sea water.
Finds photon paths from a nanobeacon source to a target PMT that have a non-zero probability.
double sourceScore(const JPhotonPath &p, const JDirectedSource *src)
JPhotonPath solve(int nscat, const JDirectedSource *src, const JPMTTarget *trg)
double targetScore(const JPhotonPath &p, const JPMTTarget *trg)
double getScore(const JPhotonPath &p, const JDirectedSource *src, const JPMTTarget *trg)
double intersectScore(const JPhotonPath &p, 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
const JPMT & getPMT(const int index) const
Get PMT.
Definition: JModule.hh:172
Axis object.
Definition: JAxis3D.hh:41
double getDistance(const JVector3D &pos) const
Get distance.
Definition: JAxis3D.hh:213
const JDirection3D & getDirection() const
Get direction.
double getDot(const JAngle3D &angle) const
Get dot product.
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
JVector3D & add(const JVector3D &vector)
Add vector.
Definition: JVector3D.hh:142
double getY() const
Get y position.
Definition: JVector3D.hh:104
double getDistance(const JVector3D &pos) const
Get distance to point.
Definition: JVector3D.hh:270
double getZ() const
Get z position.
Definition: JVector3D.hh:115
JVector3D & div(const double factor)
Scale vector.
Definition: JVector3D.hh:190
double getX() const
Get x position.
Definition: JVector3D.hh:94
Data structure for normalised vector in three dimensions.
Definition: JVersor3D.hh:28
double getDY() const
Get y direction.
Definition: JVersor3D.hh:106
double getDX() const
Get x direction.
Definition: JVersor3D.hh:95
double getDot(const JVersor3D &versor) const
Get dot product.
Definition: JVersor3D.hh:156
double getDZ() const
Get z direction.
Definition: JVersor3D.hh:117
virtual void open(const char *file_name) override
Open file.
virtual void close()=0
Close device.
General exception.
Definition: JException.hh:24
virtual bool put(const T &object)=0
Object output.
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...
JVersor3D getDirection() const
return the source direction
double getOpeningAngle() const
return the opening angle in radians
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.
double getOpeningAngle() const
return the PMT opening angle in radians
A photon path.
Definition: JPhotonPath.hh:38
bool hitsSphere(const JPosition3D &pos, double r)
Returns whether the photon path intersects a sphere of radius r at position pos.
Definition: JPhotonPath.hh:176
string getString()
get a string with the path vertex positions
Definition: JPhotonPath.hh:136
double getLength()
get the total path length
Definition: JPhotonPath.hh:104
Implementation of the JScatteringModel interface with Rayleigh scattering.
Virtual base class for a scattering model.
void setPosition(JPosition3D &_pos)
const JPosition3D & getPosition() const
Virtual base class for a light detector ("photon target").
const JVersor3D & getDirection() const
const JPosition3D & getPosition() const
void setPosition(JPosition3D &_pos)
void setRadius(double _r)
Utility class to parse command line options.
Definition: JParser.hh:1714
int read(const int argc, const char *const argv[])
Parse the program's command line options.
Definition: JParser.hh:2008
virtual double getIndexOfRefractionGroup(const double lambda) const
Index of refraction for group velocity.
Implementation of dispersion for water in deep sea.
Definition: JDispersion.hh:28
const double sigma[]
Definition: JQuadrature.cc:74
const double xmax
Definition: JQuadrature.cc:24
const double xmin
Definition: JQuadrature.cc:23
file Auxiliary data structures and methods for detector calibration.
Definition: JAnchor.hh:12
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Auxiliary classes and methods for language specific functionality.
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.
Auxiliary methods for light properties of deep-sea water.
const double getSpeedOfLight()
Get speed of light.
data_type w[N+1][M+1]
Definition: JPolint.hh:867
data_type r[M+1]
Definition: JPolint.hh:868
data_type v[N+1][M+1]
Definition: JPolint.hh:866
Name space for KM3NeT.
Definition: Jpp.hh:16
Definition: JSTDTypes.hh:14