Jpp 19.3.0-rc.1
the software that should make you happy
Loading...
Searching...
No Matches
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*******************************************************************************/
20const 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"
38#include "JPhysics/KM3NeT.hh"
39
40// root
41#include "TH1D.h"
42#include "TRandom3.h"
43#include "TFile.h"
44// namespaces
45using namespace std ;
46using namespace JLANG ;
47using namespace JDETECTOR ;
48using namespace JMARKOV ;
49using namespace JPHYSICS ;
50using namespace KM3NET ;
51
52/// find module with a given string and floor number
53const JModule& getModule( const JDetector& detector, const JModuleLocation& location ) ;
54
55/// get mean of vector content
56double 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
73void 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
85void 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 **/
128JAxis3D 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
389int 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) ;
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
928const 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 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:2142
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.
const JPosition3D & getPosition() const
Get position.
Data structure for vector in three dimensions.
Definition JVector3D.hh:36
double getY() const
Get y position.
Definition JVector3D.hh:104
JVector3D & add(const JVector3D &vector)
Add vector.
Definition JVector3D.hh:142
double getLength() const
Get length.
Definition JVector3D.hh:246
JVector3D & div(const double factor)
Scale vector.
Definition JVector3D.hh:190
double getDistance(const JVector3D &pos) const
Get distance to point.
Definition JVector3D.hh:270
double getZ() const
Get z position.
Definition JVector3D.hh:115
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()
Close file.
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.
bool hitsSphere(const JPosition3D &pos, double r)
Returns whether the photon path intersects a sphere of radius r at position pos.
string getString()
get a string with the path vertex positions
double getLength()
get the total path length
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
void setPosition(JPosition3D &_pos)
const JPosition3D & getPosition() const
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
virtual double getIndexOfRefractionGroup(const double lambda) const
Index of refraction for group velocity.
Implementation of dispersion for water in deep sea.
file Auxiliary data structures and methods for detector calibration.
Definition JAnchor.hh:12
const JModule & getModule(const JType< JDetector_t > type, const int id, const JLocation &location=JLocation())
Get module according given detector type.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
JPosition3D getPosition(const JModule &first, const JModule &second)
Get position to go from first to second module.
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.
Name space for KM3NeT.
Definition Jpp.hh:16