Jpp  18.0.0-rc.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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
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) ;
681  printResult(result) ;
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 }
const double xmax
Definition: JQuadrature.cc:24
Utility class to parse command line options.
Definition: JParser.hh:1514
General exception.
Definition: JException.hh:23
data_type w[N+1][M+1]
Definition: JPolint.hh:778
The JMarkovPhotonTracker generates ensembles of photon paths by tracking photons from a source to a t...
The JMarkovPathGenerator generates ensembles of photon paths using a Markov Chain Monte Carlo (MCMC)...
int main(int argc, char *argv[])
Definition: Main.cc:15
virtual void open(const char *file_name) override
Open file.
Auxiliary methods for PDF calculations.
Data structure for a composite optical module.
Definition: JModule.hh:68
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 setRadius(double _r)
double getDot(const JNeutrinoDirection &first, const JNeutrinoDirection &second)
Dot product.
Definition: JAstronomy.hh:409
void writePaths(string ofname, vector< JPhotonPath > paths)
double getRadius() const
const JPosition3D & getPosition() const
Detector data structure.
Definition: JDetector.hh:89
Implementation of dispersion for water in deep sea.
Definition: JDispersion.hh:26
const JDirection3D & getDirection() const
Get direction.
double getSigma(vector< double > &v)
get standard deviation of vector content
A photon path.
Definition: JPhotonPath.hh:38
Router for direct addressing of module data in detector data structure.
Implementation of the JTargetModel class that represents a spherically symmetric target.
int getNsteps()
get the number of Markov steps taken in the last call to generateEnsemble (after burn-in) ...
o $QUALITY_ROOT d $DEBUG!CHECK_EXIT_CODE JPlot1D f
Definition: JDataQuality.sh:76
std::vector< JPhotonPath > trackPhotons(unsigned long int n, JSourceModel *src, JScatteringModel *sm, JTargetModel *trg, double lambda_abs)
Track n photons.
exit
Definition: JPizza.sh:36
data_type r[M+1]
Definition: JPolint.hh:779
Implementation of JMarkovEnsembleIntegrator interface with 1D histograms.
double getDistance(const JVector3D &pos) const
Get distance to point.
Definition: JVector3D.hh:270
#define KM3NET
void setPosition(JPosition3D &_pos)
double getDot(const JAngle3D &angle) const
Get dot product.
double getMean(vector< double > &v)
get mean of vector content
void setCoordinateRemapping(bool val=true)
activate or deactive coordinate remapping
int getNacceptedSteps()
get the number of accepted steps taken during the last call to generateEnsemble (after burn-in) ...
const double JBSversion
skip continue
Definition: JTuneHV.sh:92
Implementation of the JScatteringModel interface with scattering according to the Henyey-Greenstein f...
const double sigma[]
Definition: JQuadrature.cc:74
Axis object.
Definition: JAxis3D.hh:38
then usage set_variable ACOUSTICS_WORKDIR $WORKDIR set_variable FORMULA sin([0]+2 *$PI *([1]+[2]*x)*x)" set_variable DY 1.0e-8 mkdir $WORKDIR for DETECTOR in $DETECTORS[*]
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
void setTangentialStepSize_deg(double val)
Set the average step size theta in degrees for steps in the tangential direction for the scattering v...
Implementation of the JScatteringModel interface with Rayleigh scattering.
Finds photon paths from a nanobeacon source to a target PMT that have a non-zero probability.
const JPosition3D & getPosition() const
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...
string getString()
get a string with the path vertex positions
Definition: JPhotonPath.hh:136
Properties of KM3NeT PMT and deep-sea water.
JVersor3D getDirection() const
return the source direction
Implementation of the JSourceModel class that represents a directed source with a flat intensity dist...
Data structure for vector in three dimensions.
Definition: JVector3D.hh:34
double getOpeningAngle() const
return the opening angle in radians
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1989
const JVersor3D & getDirection() const
double getFractionAccepted_radial()
get the fraction of steps that were accepted when a radial step was performed during the last call to...
double getOpeningAngle() const
return the PMT opening angle in radians
double getFractionAccepted()
get the fraction of accepted steps during the last call to generateEnsemble (after burn-in) ...
return result
Definition: JPolint.hh:764
double intersectScore(const JPhotonPath &p, const JPMTTarget *trg)
void setPosition(JPosition3D &_pos)
double getY() const
Get y position.
Definition: JVector3D.hh:104
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:130
const JPMT & getPMT(const int index) const
Get PMT.
Definition: JModule.hh:173
Implementation of the JTargetModel class that represents a single PMT on a DOM.
double getFractionAccepted_tangential()
get the fraction of steps that were accepted when a tangential step was performed during the last cal...
Direct access to module in detector data structure.
$WORKDIR ev_configure_dqsimulator txt echo process $DQ_SIMULATOR $i $SOURCE_HOST[$index] csh c(setenv ROOTSYS $ROOTSYS &&source $JPP_DIR/setenv.csh $JPP_DIR &&($DQ_SIMULATOR\-u\$NAME\$\-H\$SERVER\$\-M\$LOGGER\$\-d $DEBUG</dev/null > &/dev/null &))'
p2
Definition: module-Z:fit.sh:74
then JCookie sh JDataQuality D $DETECTOR_ID R
Definition: JDataQuality.sh:41
const double xmin
Definition: JQuadrature.cc:23
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...
double getDistance(const JVector3D &pos) const
Get distance.
Definition: JAxis3D.hh:213
virtual double getIndexOfRefractionGroup(const double lambda) const
Index of refraction for group velocity.
Implementation of the JScatteringModel interface with that combines two scattering models into one ef...
then JMuonMCEvt f $INPUT_FILE o $INTERMEDIATE_FILE d
Definition: JMuonPath.sh:47
virtual bool put(const T &object) override
Object output.
const double getSpeedOfLight()
Get speed of light.
int read(const int argc, const char *const argv[])
Parse the program&#39;s command line options.
Definition: JParser.hh:1811
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Utility class to parse command line options.
double getX() const
Get x position.
Definition: JVector3D.hh:94
Virtual base class for a scattering model.
virtual void close()
Close file.
Virtual base class for a light detector (&quot;photon target&quot;).
void printResult(vector< double > &v)
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
double targetScore(const JPhotonPath &p, const JPMTTarget *trg)
JPhotonPath solve(int nscat, const JDirectedSource *src, const JPMTTarget *trg)
void setRadialStepSize_m(double val)
set the average step size in [m] in the radial direction for the scattering vertices ...
JVector3D & div(const double factor)
Scale vector.
Definition: JVector3D.hh:190
do set_variable DETECTOR_TXT $WORKDIR detector
data_type v[N+1][M+1]
Definition: JPolint.hh:777
Data structure for normalised vector in three dimensions.
Definition: JVersor3D.hh:26
double sourceScore(const JPhotonPath &p, const JDirectedSource *src)
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.
double getZ() const
Get z position.
Definition: JVector3D.hh:115
double getLength()
get the total path length
Definition: JPhotonPath.hh:104
JVector3D & add(const JVector3D &vector)
Add vector.
Definition: JVector3D.hh:142
double getScore(const JPhotonPath &p, const JDirectedSource *src, const JPMTTarget *trg)
const JModule & getModule(const JDetector &detector, const JModuleLocation &location)
find module with a given string and floor number
vector< double > integrate(int N, int nscat, JSourceModel *src, JScatteringModel *sm, JTargetModel *trg, double lambda_abs)
Integrate with N samples.
void setTargetStepSize_deg(double val)
set the average step size in degrees for the impact point on the target
void scale(vector< double > &v, double c)
scale vector content