JMarkovPathGenerator.cc

Go to the documentation of this file.

#define MPG_VERSION 2
 
/*******************************************************************************
 * \file Markov Chain Monte Carlo to generate photon paths
 
  The result is saved to a file.
*******************************************************************************/
#include<iostream>
#include<sstream>
#include<iomanip>
#include<vector>
#include<cmath>
#include<cstdlib>
#include<string>
 
#include "TRandom3.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TFile.h"
#include "TTree.h"
#include "TBranch.h"
#include "TGraph.h"
 
// JPP includes
#include "Jeep/JParser.hh"
#include "JMarkov/JMarkovPathGenerator.hh"
#include "JMarkov/JPhotonPathWriter.hh"
#include "JPhysics/KM3NeT.hh"
//#include "JMarkov/JScatteringModel.hh"
 
// this little line removes the make_field macro defined in JParser.hh
// so that we can call the
// JPARSER::make_field(T,string) function ourselves
#undef make_field
using namespace JPP;
 
using namespace JMARKOV ;
using namespace std ;
 
//ClassImp(ScatteringModel)
 
 
int main( int argc, char** argv ) {
  cout << "JMarkovPathGenerator version " << MPG_VERSION << endl 
       << "Written by Martijn Jongen" << endl 
       << endl ; 
  cout << "Type '" << argv[0] << " -h!' to display the command-line options." << endl ;
  cout << endl ;
 
  string outfile = "out.paths"  ; // output file
  string smoutfile = ""        ; // output file to save scattering model to
  double lHG = 60.241          ; // scattering length for Henyey-Greenstein model
  double g   = 0.924           ; // parameter g in Heyney-Greenstein model
  double lR  = 294.118         ; // scattering length for Rayleigh scattering
  double a   = 0.853           ; // a parameter for Rayleigh scatteing
  double lA  = 50              ; // absorption length
  double d = 37                ; // distance to target
  int npaths = 100             ; // number of paths to simulate
  int nscat = 2                ; // number of scatterings
  int interval = 1000          ; // number of MCMC steps between saved paths
  int burnIn = 50000           ; // number of MCMC steps for initialization
  double stepsize = 10         ; // size of the MCMC steps
  bool sourceNB                ;
  double target_zenith         ;
 
  try {
    JParser<string> zap ; // this argument parser can handle strings    
    zap["o"] = make_field(outfile,"output file name") ;
    zap["O"] = make_field(smoutfile,"OPTIONAL: output file name for scattering model") ;
    zap["-lH"] = make_field(lHG,"scattering length in m for Henyey-Greenstein scattering") ;
    zap["g"] = make_field(g,"parameter g for Henyey-Greenstein function") ;
    zap["R"] = make_field(lR,"scattering length in m for Rayleigh scattering") ;
    zap["a"] = make_field(a,"parameter a for Rayleigh scattering") ;
    zap["A"] = make_field(lA,"absorption length in m") ;
    zap["d"] = make_field(d,"distance between source and target in m") ;
    zap["n"] = make_field(npaths,"number of paths to generate") ;
    zap["N"] = make_field(nscat,"number of scatterings") ;
    zap["i"] = make_field(interval,"number of MCMC steps to take between saving paths") ;
    zap["b"] = make_field(burnIn,"number of burn-in steps") ;
    zap["s"] = make_field(stepsize,"step size for the MCMC steps") ;
    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") ;
    zap["target_zenith"] = make_field(target_zenith,"[degrees] OPTIONAL: set to  use a realistic PMT acceptance") = -1 ;
 
    if (zap.read(argc, argv) != 0) {
      return 1 ;
    }
  }
  catch(const exception &error) {
    // ignore exceptions
  }
 
  // inform user of the settings
  cout << "output file name   = '" << outfile << "'." << endl ;
  cout << "absorption length  = " << lA << " m" << endl ;
  cout << "distance to target = " << d  << " m" << endl ;
  cout << "npaths             = " << npaths << endl ;
  cout << "nscat              = " << nscat << endl ;
  cout << "interval           = " << interval << " steps" << endl ;
  cout << "burn-in            = " << burnIn << " steps" << endl ;
  cout << "step size          = " << stepsize << " m" << endl ;
  cout << endl ;
 
  // Henyey-Greenstein scattering model
  JScatteringModel smHG ;
  cout << "Henyey-Greenstein scattering length = " << lHG << " m" << endl ;
  cout << "g = " << g << endl ;
  smHG.setScatteringLength(lHG) ;
  smHG.setScatteringProfileHG(g) ;
  cout << endl ;
 
  // Rayleigh scattering model
  JScatteringModel smR ;
  cout << "Rayleigh scattering length = " << lR << " m" << endl ;
  cout << "a = " << a << endl ;
  smR.setScatteringLength(lR) ;
  smR.setScatteringProfileRayleigh(a) ;
  cout << endl ;
 
  // combined scattering model
  cout << "Combining HG and Rayleigh scattering into a single effective model." << endl ;
  JScatteringModel sm ;
  setEffectiveScatteringModel( smHG, smR, sm ) ;
  cout << "Effective scattering length = " << sm.getScatteringLength() << " m" << endl ;
  sm.setAbsorptionLength(lA) ;
 
  cout << "Integral over scattering profile = " << sm.hscat->Integral("width")*2*M_PI << " (should be 1)" << endl ;
  cout << endl ;
 
  // set source distribution
  if( sourceNB ) {
    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 ;
    // we oversample each bin
    const int n = 100 ;
    for( Int_t xbin=1; xbin<=sm.hsource->GetNbinsX() ; ++xbin ) {
      double val = 0 ;
      double xmin = sm.hsource->GetXaxis()->GetBinLowEdge(bin) ;
      double xmax = sm.hsource->GetXaxis()->GetBinUpEdge(bin) ;
      for( int i=0 ; i<n ; ++i ) {
        double ct = xmin + (xmax-xmin)*(i+0.5)/n ;
        if( ct >= sqrt(0.5) ) val += 1 ;
      }
      val /= n ;
      for( Int_t ybin=1 ; ybin<=sm.hsource->GetNbinsY() ; ++ybin ) {
        Int_t bin = sm.hsource->GetBin(xbin,ybin) ;
        sm.hsource->SetBinContent(bin,val) ;
      }
    }
    sm.hsource->Scale( sqrt(2)/(2*M_PI*(sqrt(2)-1)) ) ;
    cout << endl ;
  }
 
  cout << "Integral over source profile = " << sm.hsource->Integral("width") << " (should be 1)" << endl ;
  cout << endl ;
 
  // set PMT efficiency
  if( target_zenith >= 0 ) {
    cout << "Setting target to a KM3NeT PMT." << endl 
         << "Its orientation is rotated " << target_zenith << " degrees w.r.t. the negative z-axis" << endl
         << "(the rotation is in the yz-plane)" << endl ;
    target_zenith *= M_PI/180 ; // convert to radians
    JGEOMETRY3D::JVersor3D pmtdir( sin(target_zenith), 0, cos(target_zenith) ) ;
    for( Int_t xbin=1; xbin<=sm.htarget->GetXaxis()->GetNbins() ; ++xbin ) {
      double ct = sm.htarget->GetXaxis()->GetBinCenter(xbin) ;
      double theta = acos(ct) ;
      for( Int_t ybin=1; ybin<=sm.htarget->GetYaxis()->GetNbins() ; ++ybin ) {
        Int_t bin = sm.htarget->GetBin(xbin,ybin) ;
        double phi = sm.htarget->GetYaxis()->GetBinCenter(ybin) ;
        
        JGEOMETRY3D::JVersor3D testdir( cos(phi)*sin(theta), 
                                        sin(phi)*sin(theta), 
                                        cos(theta) ) ;
        double effct = -testdir.getDot(pmtdir) ;
        double val = KM3NET::getAngularAcceptance(effct) ;
        // safety catch: I never want my probability to be exactly 0
        // so I just set it to a really low number!
        if( val == 0 ) val = 1e-10 ;
        sm.htarget->SetBinContent(bin,val) ;
      }
    }
    double maxval = sm.htarget->GetBinContent( sm.htarget->GetMaximumBin() ) ;
    sm.htarget->Scale(1.0/maxval) ;
    cout << endl ;
  }
 
  // create output root file containing the scattering model
  if( smoutfile != "" ) {
    cout << "Saving scattering model to '" << smoutfile << "'." << endl ;
    TFile* fout = new TFile(smoutfile.c_str(),"recreate") ;
    sm.Write() ;
    // we also save the source, scattering and target histograms to a separate
    // folder, so that we can view them immediately from a TBrowser
    fout->mkdir("ScatteringModel_ingredients")->cd() ;
    sm.hsource->Write() ;
    sm.hscat->Write() ;
    sm.htarget->Write() ;
    fout->cd() ;
    fout->Close() ;
    cout << endl ;
  }
 
  // generate an ensemble of paths
  cout << "Generating ensemble" << endl ;
  JMarkovPathGenerator mpg ;
  vector<JPhotonPath> ensemble = mpg.generateEnsemble( npaths, nscat, sm, burnIn, interval, d, stepsize ) ;
  cout << "Done generating ensemble." << endl ;
  double acceptance = mpg.getFractionAccepted() ;
  cout << 100*acceptance << "% of steps were accepted" << endl
       << "(as a rule of thumb, ~23% is optimal for high-dimensional spaces)" << endl 
       << endl ;
 
  // write the generated photon paths to a file 
  cout << "Writing generated ensemble to '" << outfile << "'." << endl ;
  JMARKOV::JPhotonPathWriter writer ;
  writer.open(outfile.c_str()) ;
  for( vector<JPhotonPath>::iterator it=ensemble.begin() ; it!=ensemble.end() ; ++it ) {
    writer.put( *it ) ;
  }
  writer.close() ;
  cout << endl ;
 
  cout << "Done!" << endl ;
 
  return 0 ;
}