29 #include "JDetector/JModuleLocation.hh"
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 ) {
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)
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 ;
396 JModuleLocation source_location ;
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) ;
475 load(detectorFile, detector);
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
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 ;
529 JModuleLocation target_location(target_pmt[0],target_pmt[1]) ;
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 ;
JAxis3D getBeaconAxis(const JModule &dom, double r, double beacon_theta, double beacon_phi)
Get the position and orientation of the nanobeacon on a DOM.
JPosition3D getDOMcenter(const JModule &dom)
Returns a good guess of where the actual center of the DOM is based on the PMT positions and directio...
int main(int argc, char **argv)
double getMean(vector< double > &v)
get mean of vector content
void printResult(vector< double > &v)
const JModule & getModule(const JDetector &detector, const JModuleLocation &location)
find module with a given string and floor number
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.
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
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)
Router for direct addressing of module data in detector data structure.
Data structure for a composite optical module.
const JPMT & getPMT(const int index) const
Get PMT.
double getDistance(const JVector3D &pos) const
Get distance.
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.
JVector3D & add(const JVector3D &vector)
Add vector.
double getY() const
Get y position.
double getDistance(const JVector3D &pos) const
Get distance to point.
double getZ() const
Get z position.
JVector3D & div(const double factor)
Scale vector.
double getX() const
Get x position.
Data structure for normalised vector in three dimensions.
double getDY() const
Get y direction.
double getDX() const
Get x direction.
double getDot(const JVersor3D &versor) const
Get dot product.
double getDZ() const
Get z direction.
virtual void open(const char *file_name) override
Open file.
virtual void close()=0
Close device.
virtual bool put(const T &object)=0
Object output.
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
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
const JPosition3D & getPosition() const
void setPosition(JPosition3D &_pos)
void setRadius(double _r)
Utility class to parse command line options.
int read(const int argc, const char *const argv[])
Parse the program's command line options.
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.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Auxiliary classes and methods for language specific functionality.
double getPhotonPathProbabilityDensity(JPhotonPath &p, JSourceModel *src, JScatteringModel *sm, JTargetModel *trg, double lambda_abs)
Return the probability density for a photon path with the given ingredients.
Auxiliary methods for light properties of deep-sea water.
const double getSpeedOfLight()
Get speed of light.