46 using namespace JLANG ;
60 return sum /
v.size() ;
68 varsum += (*it-mean)*(*it-mean) ;
69 return sqrt( varsum /
v.size() ) ;
81 cout <<
"Value = " << mean <<
" +- " << sigma/sqrt(1.0*
v.size()) << endl ;
87 writer.
open(ofname.c_str()) ;
92 cout <<
"paths written to '" << ofname <<
"'." << endl ;
101 for(
int pmt1=0 ; pmt1<31 ; ++pmt1 ) {
102 for(
int pmt2=pmt1+1; pmt2<31 ; ++pmt2 ) {
110 for(
int pmt=0 ; pmt<31 ; ++pmt ) {
133 const JVersor3D zdir( -neg_zdir.getDX(), -neg_zdir.getDY(), -neg_zdir.getDZ() ) ;
137 const JVersor3D ydir( tmp.cross( zdir, xdir ) ) ;
141 + sin(beacon_theta)*cos(beacon_phi)*
JVector3D(xdir)
142 + sin(beacon_theta)*sin(beacon_phi)*
JVector3D(ydir)
147 JAxis3D ax(bpos,zdir) ;
178 if( nscat == 0 )
return p ;
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() ;
189 double score = getScore(p,src,trg) ;
191 const double dlmin = 0.0001 ;
197 if( ++istep == maxnsteps ) break ;
198 double initscore = score ;
201 p[1] = p[1] + dl*dx ;
202 if( getScore(p,src,trg) < score ) {
204 p[1] = p[1] - 2*dl*dx ;
205 if( getScore(p,src,trg) < score ) p[1] = p[1] + dl*dx ;
209 p[1] = p[1] + dl*dy ;
210 if( getScore(p,src,trg) < score ) {
212 p[1] = p[1] - 2*dl*dy ;
213 if( getScore(p,src,trg) < score ) p[1] = p[1] + dl*dy ;
217 p[1] = p[1] + dl*dz ;
218 if( getScore(p,src,trg) < score ) {
220 p[1] = p[1] - 2*dl*dz ;
221 if( getScore(p,src,trg) < score ) p[1] = p[1] + dl*dz ;
225 score = getScore(p,src,trg) ;
227 if( score == initscore ) {
230 if( dl<dlmin ) break ;
242 double score = getScore(p,src,trg) ;
244 const double dlmin = 1e-10 ;
259 if( ++i>100 ) break ;
261 double initscore = score ;
263 for(
int i=0; i<nsm; ++i ) {
265 p[1] = p[1] + dl * single_moves[i] ;
266 double newscore = getScore(p,src,trg) ;
267 if( newscore < score ) {
269 p[1] = p[1] - dl * single_moves[i] ;
277 for(
int i=0; i<nsm; ++i ) {
279 p[2] = p[2] + dl * single_moves[i] ;
280 double newscore = getScore(p,src,trg) ;
281 if( newscore < score ) {
283 p[2] = p[2] - dl * single_moves[i] ;
292 score = getScore(p,src,trg) ;
294 if( score == initscore ) {
297 if( dl<dlmin ) break ;
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() ;
327 double val = max(ctmin-ct,0.0) / (1-ctmin) ;
334 double val = max(ctmin-ct,0.0) / (1-ctmin) ;
344 if( p.
n == 2 )
return score ;
346 for(
int nv=1 ; nv<(int)p.size()-2 ; ++nv ) {
353 double val = (R-d)/R ;
383 double score = sourceScore(p,src) + targetScore(p,trg) + intersectScore(p,trg) ;
389 int main(
int argc,
char** argv ) {
390 cout <<
"JBeaconSimulator version " <<
JBSversion << endl
394 string detectorFile ;
399 unsigned int nscatmax ;
400 unsigned long int nphotons ;
402 double opening_angle_deg ;
403 double beacon_theta_deg ;
404 double beacon_phi_deg ;
410 double target_step_size_deg ;
411 double tangential_step_size_deg ;
426 zap[
"N-integration-samples"] =
make_field(nsamples) = 1e6 ;
428 zap[
"write-paths"] =
make_field(write_paths) ;
430 zap[
"opening-angle"] =
make_field(opening_angle_deg) = 45 ;
431 zap[
"beacon-theta-deg"] =
make_field(beacon_theta_deg) = 56.2893-15 ;
432 zap[
"beacon-phi-deg"] =
make_field(beacon_phi_deg) = -59.9998 ;
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) ;
442 if (zap.
read(argc, argv) != 0) exit(1) ;
444 catch(
const exception &error) {
445 cerr << error.what() << endl ;
452 extensions.push_back(
".root") ;
453 extensions.push_back(
".paths") ;
455 size_t pos = ofname.find(*it) ;
456 if( pos != string::npos ) {
457 ofname = ofname.substr(0,pos) ;
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) ;
478 cout <<
"FATAL ERROR: could not open detector file '" << detectorFile <<
"'." ;
485 gRandom =
new TRandom3(0) ;
488 const double DOM_radius = 0.2159 ;
497 JAxis3D source_axis(
getBeaconAxis(source_module,DOM_radius+0.00001,beacon_theta_deg*M_PI/180,beacon_phi_deg*M_PI/180) ) ;
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
514 if( target_pmt.size() != 3 ) {
515 cerr <<
"FATAL ERROR: provide a target PMT in the format "
516 <<
"\"string floor daq_channel\"." << endl ;
522 if( target_pmt[2]<0 || target_pmt[2]>30 ) {
523 cerr <<
"FATAL ERROR: invalid target pmt index ("
524 << target_pmt[2] <<
")" << endl ;
535 opening_angle_deg*M_PI/180,
541 cout <<
"PHYSICS" << endl ;
543 cout <<
"Henyey-Greenstein scattering length = " << lHG <<
" [m], g = " << gHG << endl ;
545 cout <<
"Rayleigh scattering length = " << lR <<
" [m], a = " << aR << endl ;
547 cout <<
"Absorption length = " << Labs <<
" [m]" << endl ;
551 const double water_depth = 3000 ;
552 const double pressure = 0.1*water_depth ;
554 const double nbwavelength = 470 ;
557 cout <<
"Assuming speed of light in water = " << cw <<
" [m/ns]" << endl ;
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]
574 <<
"Target radius = " << DOM_radius <<
" [m]" << endl
576 if( useTracker ) cout <<
"all 31 PMTs" ;
577 else cout <<
"PMT " << target_pmt[2] <<
", direction = " << target_axis.getDirection() ;
579 <<
"PMTs are simulated as spherical caps" << endl
580 <<
"PMT opening angle = " << pmt_opening_angle_deg <<
" deg." << endl
592 trg =
new JPMTTarget(target_axis,pmt_opening_angle_rad,DOM_radius) ;
595 unsigned long int n_generated_photons = nphotons ;
598 cout <<
"Generating photon paths by photon tracking" << endl ;
604 cout <<
"Photon hits: " << paths.size() <<
" / " << n_generated_photons << endl ;
619 cout <<
"WARNING: coordinate remapping is deactivated!" << endl ;
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 ;
628 double stepsize = 1 ;
629 if( nscat>0 && (
int)stepsizes.size()>=nscat ) stepsize = stepsizes[nscat-1] ;
632 if( nscat != 0 ) cout <<
"step size = " << stepsize <<
" [m]" << endl ;
637 cout <<
"path = " << endl ;
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
659 nsteps_burn_in, nsteps_save ) ;
663 cout <<
"Accepted " << naccepted <<
" / " << nsteps <<
" steps (" << mpg.
getFractionAccepted()*100.0 <<
"%)" << 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 ;
668 paths.insert(paths.end(), these_paths.begin(), these_paths.end()) ;
671 if( these_paths.size() == 0 ) {
672 cerr <<
"Something went wrong. The MCMC generator returned no paths." << endl ;
676 cout <<
"Creating JMarkovEnsembleIntegrator" << endl ;
679 cout <<
"Integrating" << endl ;
685 integrators.push_back(jmi) ;
691 if( paths.size()==0 ) {
692 cerr <<
"FATAL ERROR: no paths were generated!" << endl ;
699 unsigned int npaths = paths.size() ;
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 ;
733 for(
unsigned int i=0; i<npaths; ++i ) {
738 for(
unsigned int pmt=0 ; pmt<31 ; ++pmt ) {
740 if( ct>pmt_ctmin ) hitPMTs.push_back(pmt) ;
742 if( hitPMTs.size()==1 ) {
743 PMT_is_hit[i] = hitPMTs[0] ;
744 ++nHitsPerPMT[ hitPMTs[0] ] ;
746 if( hitPMTs.size()>1 ) {
747 cerr <<
"FATAL ERROR: hit " << hitPMTs.size() <<
" PMTs simultaneously" << endl ;
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 ;
761 const double dt = 0.1 ;
762 const unsigned int margin = 5 ;
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 ) ;
773 pmtfront = target_pmt[2] ;
774 pmtback = target_pmt[2]+1 ;
777 for(
int pmt=pmtfront ; pmt<pmtback ; ++pmt ) {
780 sprintf( fname,
"%s_PMT%i.root", ofname.c_str(), pmt ) ;
781 TFile* f =
new TFile(fname,
"recreate") ;
784 TH1D hpulse(
"hpulse",
"Instantaneous flash pulse profile;time since flash [ns];hit probability / ns",
787 for(
int nscat=0 ; nscat<=(int)nscatmax ; ++nscat ) {
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",
793 TH2D hXY(
"hXY",
"hit position;X;Y",
796 hXY.SetOption(
"colz") ;
797 TH2D hXZ(
"hXZ",
"hit position;X;Z",
800 hXZ.SetOption(
"colz") ;
801 TH2D hYZ(
"hYZ",
"hit position;Y;Z",
804 hYZ.SetOption(
"colz") ;
805 TH1D hnscat(
"hnscat",
";N_{scat};probability",nscatmax+2,-0.5,nscatmax+1.5) ;
807 TH3D hregion(
"hregion",
"Photon paths;X;Y;Z",
820 for(
unsigned int i=0 ; i<npaths ; ++i ) {
821 if( PMT_is_hit[i] != pmt ) continue ;
823 paths_to_write.push_back(paths[i]) ;
825 const double path_length = paths[i].getLength() ;
826 double t = path_length/cw ;
827 int nscat = paths[i].size()-2 ;
830 if( useTracker )
w = 1.0 / n_generated_photons ;
831 else w = Pscat[nscat] / n_generated_photons ;
834 hnscat.Fill(nscat,
w) ;
835 hpulse_per_nscat[nscat]->Fill(t,
w) ;
837 hXY.Fill(xs[i],ys[i],
w) ;
838 hXZ.Fill(xs[i],zs[i],
w) ;
839 hYZ.Fill(ys[i],zs[i],
w) ;
841 if( paths.size()<10000 ) {
842 const double dl = 0.001 ;
843 const double dw =
w*dl/path_length ;
844 double l_along_seg = 0.5*dl ;
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 ) {
851 hregion.Fill( pos.
getX(), pos.
getY(), pos.
getZ(), dw ) ;
853 cerr <<
"Woohoooo" << endl ;
859 l_along_seg -= seg_length ;
865 hregion.Scale(1e20) ;
869 char path_fname[500] ;
870 sprintf( path_fname,
"%s_PMT%i.paths", ofname.c_str(), pmt ) ;
871 cout <<
"Writing paths." << endl ;
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) ;
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) ;
896 hpulse.Scale(1.0/dt) ;
897 for(
int nscat=0 ; nscat<=(int)nscatmax ; ++nscat ) {
898 hpulse_per_nscat[nscat]->Scale(1.0/dt) ;
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) ;
913 f->mkdir(
"EnsembleIntegratorHistograms")->cd() ;
915 (*it)->writeHistograms() ;
921 cout <<
"Output in '" << fname <<
"'." << endl ;
923 cout <<
"Done!" << endl ;
929 for(
unsigned int i=0; i<
detector.size(); ++i ) {
930 if(
detector[i].getLocation() == location ) {
934 cerr <<
"FATAL ERROR: no DOM found with location " << location << endl ;