44 cout <<
"JMarkovPathGenerator version " <<
MPG_VERSION << endl
45 <<
"Written by Martijn Jongen" << endl
47 cout <<
"Type '" << argv[0] <<
" -h!' to display the command-line options." << endl ;
50 string outfile =
"out.paths" ;
51 string smoutfile =
"" ;
62 double stepsize = 10 ;
64 double target_zenith ;
68 zap[
"o"] =
make_field(outfile,
"output file name") ;
69 zap[
"O"] =
make_field(smoutfile,
"OPTIONAL: output file name for scattering model") ;
70 zap[
"-lH"] =
make_field(lHG,
"scattering length in m for Henyey-Greenstein scattering") ;
71 zap[
"g"] =
make_field(g,
"parameter g for Henyey-Greenstein function") ;
72 zap[
"R"] =
make_field(lR,
"scattering length in m for Rayleigh scattering") ;
73 zap[
"a"] =
make_field(a,
"parameter a for Rayleigh scattering") ;
74 zap[
"A"] =
make_field(lA,
"absorption length in m") ;
75 zap[
"d"] =
make_field(d,
"distance between source and target in m") ;
76 zap[
"n"] =
make_field(npaths,
"number of paths to generate") ;
77 zap[
"N"] =
make_field(nscat,
"number of scatterings") ;
78 zap[
"i"] =
make_field(interval,
"number of MCMC steps to take between saving paths") ;
79 zap[
"b"] =
make_field(burnIn,
"number of burn-in steps") ;
80 zap[
"s"] =
make_field(stepsize,
"step size for the MCMC steps") ;
81 zap[
"sourceNB"] =
make_field(sourceNB,
"Use a nanobeacon profile as source (currently a uniform distribution in a 45 degree cone around the positive z-direction") ;
82 zap[
"target_zenith"] =
make_field(target_zenith,
"[degrees] OPTIONAL: set to use a realistic PMT acceptance") = -1 ;
84 if (zap.
read(argc, argv) != 0) {
88 catch(
const exception &error) {
93 cout <<
"output file name = '" << outfile <<
"'." << endl ;
94 cout <<
"absorption length = " << lA <<
" m" << endl ;
95 cout <<
"distance to target = " << d <<
" m" << endl ;
96 cout <<
"npaths = " << npaths << endl ;
97 cout <<
"nscat = " << nscat << endl ;
98 cout <<
"interval = " << interval <<
" steps" << endl ;
99 cout <<
"burn-in = " << burnIn <<
" steps" << endl ;
100 cout <<
"step size = " << stepsize <<
" m" << endl ;
105 cout <<
"Henyey-Greenstein scattering length = " << lHG <<
" m" << endl ;
106 cout <<
"g = " << g << endl ;
107 smHG.setScatteringLength(lHG) ;
108 smHG.setScatteringProfileHG(g) ;
113 cout <<
"Rayleigh scattering length = " << lR <<
" m" << endl ;
114 cout <<
"a = " << a << endl ;
115 smR.setScatteringLength(lR) ;
116 smR.setScatteringProfileRayleigh(a) ;
120 cout <<
"Combining HG and Rayleigh scattering into a single effective model." << endl ;
122 setEffectiveScatteringModel( smHG, smR, sm ) ;
124 sm.setAbsorptionLength(lA) ;
126 cout <<
"Integral over scattering profile = " << sm.hscat->Integral(
"width")*2*M_PI <<
" (should be 1)" << endl ;
131 cout <<
"Setting source distribution to a preliminary approximation of a nanobeacon profile. Light is emitted uniformly within a cone around the positive z-axis." << endl ;
134 for( Int_t xbin=1; xbin<=sm.hsource->GetNbinsX() ; ++xbin ) {
136 double xmin = sm.hsource->GetXaxis()->GetBinLowEdge(bin) ;
137 double xmax = sm.hsource->GetXaxis()->GetBinUpEdge(bin) ;
138 for(
int i=0 ; i<
n ; ++i ) {
139 double ct = xmin + (xmax-
xmin)*(i+0.5)/
n ;
140 if( ct >= sqrt(0.5) ) val += 1 ;
143 for( Int_t ybin=1 ; ybin<=sm.hsource->GetNbinsY() ; ++ybin ) {
144 Int_t bin = sm.hsource->GetBin(xbin,ybin) ;
145 sm.hsource->SetBinContent(bin,val) ;
148 sm.hsource->Scale( sqrt(2)/(2*M_PI*(sqrt(2)-1)) ) ;
152 cout <<
"Integral over source profile = " << sm.hsource->Integral(
"width") <<
" (should be 1)" << endl ;
156 if( target_zenith >= 0 ) {
157 cout <<
"Setting target to a KM3NeT PMT." << endl
158 <<
"Its orientation is rotated " << target_zenith <<
" degrees w.r.t. the negative z-axis" << endl
159 <<
"(the rotation is in the yz-plane)" << endl ;
160 target_zenith *= M_PI/180 ;
162 for( Int_t xbin=1; xbin<=sm.htarget->GetXaxis()->GetNbins() ; ++xbin ) {
163 double ct = sm.htarget->GetXaxis()->GetBinCenter(xbin) ;
164 double theta = acos(ct) ;
165 for( Int_t ybin=1; ybin<=sm.htarget->GetYaxis()->GetNbins() ; ++ybin ) {
166 Int_t bin = sm.htarget->GetBin(xbin,ybin) ;
167 double phi = sm.htarget->GetYaxis()->GetBinCenter(ybin) ;
172 double effct = -testdir.getDot(pmtdir) ;
176 if( val == 0 ) val = 1e-10 ;
177 sm.htarget->SetBinContent(bin,val) ;
180 double maxval = sm.htarget->GetBinContent( sm.htarget->GetMaximumBin() ) ;
181 sm.htarget->Scale(1.0/maxval) ;
186 if( smoutfile !=
"" ) {
187 cout <<
"Saving scattering model to '" << smoutfile <<
"'." << endl ;
188 TFile* fout =
new TFile(smoutfile.c_str(),
"recreate") ;
192 fout->mkdir(
"ScatteringModel_ingredients")->cd() ;
193 sm.hsource->Write() ;
195 sm.htarget->Write() ;
202 cout <<
"Generating ensemble" << endl ;
205 cout <<
"Done generating ensemble." << endl ;
207 cout << 100*acceptance <<
"% of steps were accepted" << endl
208 <<
"(as a rule of thumb, ~23% is optimal for high-dimensional spaces)" << endl
212 cout <<
"Writing generated ensemble to '" << outfile <<
"'." << endl ;
214 writer.
open(outfile.c_str()) ;
221 cout <<
"Done!" << endl ;
Utility class to parse command line options.
The JMarkovPathGenerator generates ensembles of photon paths using a Markov Chain Monte Carlo (MCMC)...
virtual void open(const char *file_name) override
Open file.
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...
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
double getFractionAccepted()
get the fraction of accepted steps during the last call to generateEnsemble (after burn-in) ...
then JMuonMCEvt f $INPUT_FILE o $INTERMEDIATE_FILE d
virtual bool put(const T &object) override
Object output.
int read(const int argc, const char *const argv[])
Parse the program's command line options.
double getAngularAcceptance(const double x)
Get angular acceptance of PMT.
Virtual base class for a scattering model.
virtual void close()
Close file.
virtual double getScatteringLength()
Data structure for normalised vector in three dimensions.