/**
 * \author mjongen
 */

/*******************************************************************************
  JBeaconSimulator

  Note that this program is still under development.

  It uses the JMarkov library to simulate the response of a given PMT to an
  instantaneous flash of a nanobeacon.

  -s   "source", a string like "2 1", indicating the string and floor of the 
       DOM on which the beacon is activated
 
  -t   "target", a string like "1 10 3", indicating the string, floor and DAQ
       channel of the target PMT

*******************************************************************************/
const double JBSversion = 3.21 ;

// c++ standard library
#include <string>

// JPP
#include "Jeep/JParser.hh"
#include "JDetector/JModuleRouter.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleLocation.hh"
#include "JGeometry3D/JAxis3D.hh"
#include "JMarkov/JScatteringModel.hh"
#include "JMarkov/JMarkovPhotonTracker.hh"
#include "JMarkov/JMarkovPathGenerator.hh"
#include "JMarkov/JMarkovIntegrator.hh"
#include "JMarkov/JPhotonPathWriter.hh"
#include "JPhysics/JPDFToolkit.hh"
#include "JPhysics/JDispersion.hh"
#include "JPhysics/KM3NeT.hh"

// root
#include "TH1D.h"
#include "TRandom3.h"
#include "TFile.h"
// namespaces
using namespace std ;
using namespace JLANG ;   
using namespace JDETECTOR ;
using namespace JMARKOV ;
using namespace JPHYSICS ;
using namespace KM3NET ;

/// find module with a given string and floor number
const JModule& getModule( const JDetector& detector, const JModuleLocation& location ) ;

/// get mean of vector content
double getMean( vector<double>& v ) {
  double sum = 0 ;
  for( vector<double>::iterator it=v.begin() ; it!=v.end() ; ++it ) 
    sum += *it ;
  return sum / v.size() ;
}

/// get standard deviation of vector content
double getSigma( vector<double>& v ) {
  double mean = getMean(v) ;
  double varsum = 0 ;
  for( vector<double>::iterator it=v.begin() ; it!=v.end() ; ++it ) 
    varsum += (*it-mean)*(*it-mean) ;
  return sqrt( varsum / v.size() ) ;
}

/// scale vector content
void scale( vector<double>& v, double c ) {
  for( vector<double>::iterator it=v.begin() ; it!=v.end() ; ++it )
    *it = c * (*it) ;
}

void printResult( vector<double>& v ) {
  double mean = getMean(v) ;
  double sigma = getSigma(v) ;
  cout << "Value = " << mean << " +- " << sigma/sqrt(1.0*v.size()) << endl ;
}

// write a vector of generated photon paths to a file 
void writePaths( string ofname, vector<JPhotonPath> paths ) {
  JMARKOV::JPhotonPathWriter writer ;
  writer.open(ofname.c_str()) ;
  for( vector<JPhotonPath>::iterator it=paths.begin() ; it!=paths.end() ; ++it ) {
    writer.put( *it ) ;
  }
  writer.close() ;
  cout << "paths written to '" << ofname << "'." << endl ;
}

/**
   Returns a good guess of where the actual center of the DOM is based on the PMT positions
   and directions.
 **/
JPosition3D getDOMcenter( const JModule& dom ) {
  double r = 0 ;
  for( int pmt1=0 ; pmt1<31 ; ++pmt1 ) {
    for( int pmt2=pmt1+1; pmt2<31 ; ++pmt2 ) { 
      double num = (dom.getPMT(pmt1).getPosition()-dom.getPMT(pmt2).getPosition()).getLength() ;
      double den = (JVector3D(dom.getPMT(pmt1).getDirection())-JVector3D(dom.getPMT(pmt2).getDirection())).getLength() ;
      r += num/den ;
    }
  }
  r /= 0.5*31*30 ;
  JPosition3D center(0,0,0) ;
  for( int pmt=0 ; pmt<31 ; ++pmt ) {
    center.add( dom.getPMT(pmt).getPosition() - r*JVector3D(dom.getPMT(pmt).getDirection()) ) ;
  }
  center.div(31) ;
  return center ;
}


/**
   Get the position and orientation of the nanobeacon on a DOM.

   The local coordinate system of the DOM is deduced from the PMT
   positions.

   r            = DOM radius [m]
   beacon_theta = theta angle of the beacon in the DOM's coordinate system [rad]
   phi_theta    = phi angle of the beacon in the DOM's coordinate system [rad]
 **/
JAxis3D getBeaconAxis( const JModule& dom, double r, double beacon_theta, double beacon_phi ) {
  JPosition3D dom_center = getDOMcenter(dom) ; 

  // PMT 22 defines the negative z-direction of the DOM coordinate system
  const JVersor3D neg_zdir( dom.getPMT(22).getDirection() ) ;
  const JVersor3D zdir( -neg_zdir.getDX(), -neg_zdir.getDY(), -neg_zdir.getDZ() ) ;
  // PMT 0 points in the negative y-direction of the DOM coordinate system (+some z-component)
  JVector3D tmp ;
  const JVersor3D xdir(  tmp.cross( zdir, dom.getPMT(0).getDirection() ) ) ;
  const JVersor3D ydir(  tmp.cross( zdir, xdir )  ) ;

  // direction of beacon's position w.r.t. the DOM center
  JVector3D _bdir( cos(beacon_theta)*JVector3D(zdir) 
		   + sin(beacon_theta)*cos(beacon_phi)*JVector3D(xdir)
		   + sin(beacon_theta)*sin(beacon_phi)*JVector3D(ydir)
		   ) ;
  JVersor3D bdir(_bdir) ; // start out in DOM z-direction

  JPosition3D bpos( dom_center + r*JVector3D(bdir) ) ;
  JAxis3D ax(bpos,zdir) ;

  return ax ;
}

/**
   Finds photon paths from a nanobeacon source to a target PMT that 
   have a non-zero probability. This means that 
   - the emission angle has to be within the nanobeacon opening angle
   - the photon has to hit the PMT at an angle less than 90 degrees
   - the photon may not intersect the target

   It works by maximizing a score function. If we find a score of 0, 
   the path meets all criteria.
 **/
class FreePathSolver {
  public :
  
  FreePathSolver() {}

  JPhotonPath solve(int nscat, const JDirectedSource* src, const JPMTTarget* trg ) {
    //cout << "Running FreePathSolver::solve with nscat = " << nscat << endl ;

    JPhotonPath p(min(2,nscat)) ;

    // start at source
    p.front() = src->getPosition() ; 
    // end on PMT surface
    p.back()  = trg->getPosition() ;
    if( trg->getRadius()>0 ) p.back() = trg->getPosition() + trg->getRadius() * JVector3D(trg->getDirection()) ;

    if( nscat == 0 ) return p ;

    // create intermediate points
    for( int nv=1 ; nv<(int)p.size()-1 ; ++nv ) {
      double factor = 1.0 * nv / (p.size()-1) ;
      p[nv] = (1-factor)*p.front() + factor*p.back() ;
    }

    int maxnsteps = -1 ;

    if( nscat == 1 ) {
      double score = getScore(p,src,trg) ;
      double dl = 1 ; // m
      const double dlmin = 0.0001 ;
      JPosition3D dx(1,0,0) ;
      JPosition3D dy(0,1,0) ;
      JPosition3D dz(0,0,1) ;
      int istep = 0 ;
      while( score<0 ) {
	if( ++istep == maxnsteps ) break ;
	double initscore = score ;

	// try step in x-direction
	p[1] = p[1] + dl*dx ;
	if( getScore(p,src,trg) < score ) {
	  // try moving in the opposite direction
	  p[1] = p[1] - 2*dl*dx ;
	  if( getScore(p,src,trg) < score ) p[1] = p[1] + dl*dx ;
	}

	// try step in y-direction
	p[1] = p[1] + dl*dy ;
	if( getScore(p,src,trg) < score ) {
	  // try moving in the opposite direction
	  p[1] = p[1] - 2*dl*dy ;
	  if( getScore(p,src,trg) < score ) p[1] = p[1] + dl*dy ; 
	}

	// try step in z-direction
	p[1] = p[1] + dl*dz ;
	if( getScore(p,src,trg) < score ) {
	  // try moving in the opposite direction
	  p[1] = p[1] - 2*dl*dz ;
	  if( getScore(p,src,trg) < score ) p[1] = p[1] + dl*dz ;
	}

	// see if it worked
	score = getScore(p,src,trg) ;
	// check if we are stuck
	if( score == initscore ) {
	  dl = 0.5 * dl ;
	  //cout << "dl = " << dl << endl ;
	  if( dl<dlmin ) break ;
	}
	if( p.getLength() > 10000 ) break ;
      }
    }

    if( nscat>1) {
      //cout << "nscat more than one" << endl ;
      //cout << "starting with " << endl ;
      //for( vector<JPosition3D>::iterator it=p.begin() ; it!= p.end() ; ++it ) cout << *it << endl ;
      //cout << endl ;

      double score = getScore(p,src,trg) ;
      double dl = 1 ; // m
      const double dlmin = 1e-10 ;

      const int nsm = 6 ;
      JPosition3D single_moves[nsm] = { JPosition3D(1,0,0), JPosition3D(0,1,0), JPosition3D(0,0,1),
      					JPosition3D(-1,0,0), JPosition3D(0,-1,0), JPosition3D(0,0,-1) } ;
      /*ctor< pair<JPosition3D,JPosition3D> > moves ;
	for( int i=0 ; i<nsm ; ++i ) {
    	for( int j=0 ; j<nsm ; ++j ) {
	moves.push_back( pair<JPosition3D,JPosition3D>(single_moves[i],single_moves[j]) ) ;
	}
	}*/

      int i = 0 ;
      
      while( score<0 ) {
	if( ++i>100 ) break ;
	//cout << "@ " << JPosition3D(p[1]-src->getPosition()) << ", " << JPosition3D(p[2]-src->getPosition()) << ", score = " << score << ", dl = " << dl << endl ;
	double initscore = score ;

	for( int i=0; i<nsm; ++i ) {
	  // try move
	  p[1] = p[1] + dl * single_moves[i]  ;
	  double newscore = getScore(p,src,trg) ;
	  if( newscore < score ) {
	    // if the score is reduced, undo the move
	    p[1] = p[1] - dl * single_moves[i]  ;
	  } else {
	    //cout << "First vertex adds " << single_moves[i] 
	    //	 << ", score difference = " << newscore - score << endl ;
	    score = newscore ;
	  }
	}

	for( int i=0; i<nsm; ++i ) {
	  // try move
	  p[2] = p[2] + dl * single_moves[i]  ;
	  double newscore = getScore(p,src,trg) ;
	  if( newscore < score ) {
	    // if the score is reduced, undo the move
	    p[2] = p[2] - dl * single_moves[i]  ;
	  } else {
	    //cout << "Second vertex adds " << single_moves[i] 
	    //	 << ", score difference = " << newscore - score << endl ;
	    score = newscore ;
	  }
	}

	// see if it worked
	score = getScore(p,src,trg) ;
	// check if we are stuck
	if( score == initscore ) {
	  //cout << "Stuck" << endl ;
	  dl = 0.5*dl ;
	  if( dl<dlmin ) break ;
	}
	if( p.getLength() > 10000 ) {
	  //cout << "Too long" << endl ;
	  break ;
	}
      }

      // fill in the gaps
      JPhotonPath p2(nscat) ;
      p2[0] = p[0] ;
      p2[1] = p[1] ;
      p2[2] = p[2] ;
      p2.back() = p.back() ;
      for( int nv=3 ; nv<nscat+1 ; ++nv ) {
	double factor = 1.0 * (nv-2) / (nscat-1) ;
	p2[nv] = (1-factor)*p2[2] + factor*p2.back() ;
      }
      return p2 ;
    }
  
    return p ;
  
  }

  protected :

  double sourceScore( const JPhotonPath& p, const JDirectedSource* src ) {
    double ct = JVersor3D(p[1]-src->getPosition()).getDot(src->getDirection()) ;
    const double ctmin = cos(0.5*src->getOpeningAngle()) ;
    double val =  max(ctmin-ct,0.0) / (1-ctmin) ;
    return -val*val ;
  }

  double targetScore( const JPhotonPath& p, const JPMTTarget* trg ) {
    double ct = JVersor3D(p[p.size()-2]-p[p.size()-1]).getDot(trg->getDirection()) ;
    const double ctmin = cos(0.5*trg->getOpeningAngle()) ;
    double val =  max(ctmin-ct,0.0) / (1-ctmin) ;
    return -val*val ;
  }

  double intersectScore( const JPhotonPath& p, const JPMTTarget* trg ) {
    double score = 0 ;

    // check whether the photon path goes through the target sphere
    // we do not need to consider the last vertex, because if that intersects the sphere
    // the target score would be 0.
    if( p.n == 2 ) return score ;

    for( int nv=1 ; nv<(int)p.size()-2 ; ++nv ) {
      JPhotonPath seg(0) ;
      seg[0] = p[nv] ;
      seg[1] = p[nv+1] ;
      if( seg.hitsSphere(trg->getPosition(),trg->getRadius()) ) {
	double d = JAxis3D( seg[0], JVersor3D(seg[1]-seg[0]) ).getDistance(trg->getPosition()) ;
	double R = trg->getRadius() ;
	double val = (R-d)/R ;
	score -= val*val ;
	//cout << "Segment from " << nv << " to " << nv+1 << "intersects sphere"
	//     << ", distance = " << d << ", val = " << val << endl ;
      }
    }

    /*JPhotonPath pshort(p) ;
    pshort.resize( pshort.size()-1 ) ;
    --pshort.n ;
    if( pshort.hitsSphere(trg->getPosition(),trg->getRadius()) ) {
      JVersor3D hit_dir( pshort.getSphereHitPosition(trg->getPosition(),trg->getRadius()) - trg->getPosition() ) ;
      double val = hit_dir.getDot(trg->getDirection()) ;
      score -= val * val ; 
    }
    
    // extra penalty when a vertex is inside the target
    for( int nv=1 ; nv<p.n-1 ; ++nv ) {
      double r =  p[nv].getDistance(trg->getPosition()) ;
      if( r < trg->getRadius() ) {
	double val = ( trg->getRadius() - r ) / trg->getRadius() ;
	score -= val * val ;
      } 
      }*/

    return score ;
  }

  // total score 
  double getScore( const JPhotonPath& p, const JDirectedSource* src, const JPMTTarget* trg  ) {
    double score = sourceScore(p,src) + targetScore(p,trg) + intersectScore(p,trg) ;
    return score ;
  }

} ;

int main(  int argc, char** argv ) {
  cout << "JBeaconSimulator version " << JBSversion << endl
       << endl ;

  string ofname ;
  string detectorFile ;
  
  JModuleLocation source_location ;
  vector<int> target_pmt ;
  bool useTracker ;
  unsigned int nscatmax ;
  unsigned long int nphotons ;
  double Labs ;
  double opening_angle_deg ;
  double beacon_theta_deg ;
  double beacon_phi_deg ;
  double lHG ;
  double gHG ;
  double lR  ;
  double aR  ;
  vector<double> stepsizes ;
  double target_step_size_deg ;
  double tangential_step_size_deg ;
  const vector<double> empty_vector ;
  int nsamples ;
  bool write_paths ;
  bool no_remapping ;

  try { 
    JParser<string> zap;
    
    zap["o"] = make_field(ofname) ; // output file name
    zap["a"] = make_field(detectorFile) ;
    zap["s"] = make_field(source_location) ;
    zap["t"] = make_field(target_pmt) ;
    zap["m"] = make_field(nscatmax) = 5 ;
    zap["N"] = make_field(nphotons) ;
    zap["N-integration-samples"] = make_field(nsamples) = 1e6 ;
    zap["tracker"] = make_field(useTracker) ;
    zap["write-paths"] = make_field(write_paths) ;
    zap["Labs"] = make_field(Labs) = 50 ;
    zap["opening-angle"] = make_field(opening_angle_deg) = 45 ;
    zap["beacon-theta-deg"] = make_field(beacon_theta_deg) = 56.2893-15 ; // 15 degrees above theta of upper row PMTs
    zap["beacon-phi-deg"] = make_field(beacon_phi_deg) = -59.9998 ; // phi of PMT 3
    zap["LscatHG"] = make_field(lHG) = 60.241 ;
    zap["gHG"] = make_field(gHG) = 0.924 ;
    zap["LscatR"] = make_field(lR) = 294.118 ;
    zap["aR"] = make_field(aR) = 0.853 ;
    zap["stepsizes"] = make_field(stepsizes) = empty_vector ;
    zap["target-step-size-deg"] = make_field(target_step_size_deg) = 1.0 ;
    zap["tangential-step-size-deg"] = make_field(tangential_step_size_deg) = 1.0 ;
    zap["no-remapping"] = make_field(no_remapping) ;

    if (zap.read(argc, argv) != 0) exit(1) ;
  }
  catch(const exception &error) {
    cerr << error.what() << endl ;
    exit(1) ;
  }

  // remove extension from output file name
  {
    vector<string> extensions ;
    extensions.push_back(".root")  ;
    extensions.push_back(".paths") ;
    for( vector<string>::iterator it=extensions.begin(); it!=extensions.end(); ++it ) {
      size_t pos = ofname.find(*it) ;
      if( pos != string::npos ) {
	ofname = ofname.substr(0,pos) ;
      }
    }
  }

  // use default step sizes when no new ones (or too few are specified)
  if( !useTracker ) {
    const int ndefvals = 6 ;
    double defvals[ndefvals] = { 8, 5, 4, 3, 2, 1 } ;
    while( stepsizes.size()<nscatmax ) {
      double val = defvals[min(ndefvals-1,(int)stepsizes.size())] ;
      stepsizes.push_back(val) ;
    }
  }
  
  // load detector
  JDetector detector;
  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    cout << "FATAL ERROR: could not open detector file '" << detectorFile << "'." ;
    exit(1) ;
  }
  JModuleRouter moduleRouter(detector) ;


  // initialize gRandom
  gRandom = new TRandom3(0) ;

  // radius of a DOM (used for source and target model)
  const double DOM_radius = 0.2159 ; // 0.2159 [m] (=17 inch diameter glass sphere)


  // find source position and direction (stored as a JAxis3D)
  // the source is a NB, whose position is deduced from the PMT positions in the detector file
  // note that the distance from the DOM center is DOM_radius + 0.01 mm, to make sure that the 
  // initial vertex does not intersect the DOM (which is a problem when the source is located
  // on the target DOM)
  JModule source_module( getModule(detector,source_location) ) ;
  JAxis3D source_axis( getBeaconAxis(source_module,DOM_radius+0.00001,beacon_theta_deg*M_PI/180,beacon_phi_deg*M_PI/180) ) ;
  JVersor3D surface_normal( source_axis.getPosition() - getDOMcenter(source_module) ) ;

  cout << "SOURCE (NANOBEACON)" << endl 
       << "Positioned on string " << source_location.getString() << " floor " << source_location.getFloor()
       << " (DOM center @ " << getDOMcenter(source_module) << ")" << endl
       << "Position in local DOM coordinates: theta = " 
       << beacon_theta_deg << " deg, "
       << "phi = " << beacon_phi_deg << " deg" << endl 
       << "distance w.r.t. DOM center = " << (source_axis.getPosition()-getDOMcenter(source_module)).getLength() << " m" << endl 
       << "unit vector from DOM center to beacon position = " << JVersor3D(source_axis.getPosition()-getDOMcenter(source_module)) << endl
       << "beam direction = " << source_axis.getDirection() << endl
       << "beam opening angle = " << opening_angle_deg << " deg" << endl 
       << "surface normal     = " << surface_normal << endl
       << endl ;

  // find target position and direction (stored as a JAxis3D)
  if( target_pmt.size() != 3 ) {
    cerr << "FATAL ERROR: provide a target PMT in the format "
	 << "\"string floor daq_channel\"." << endl ;
    exit(1) ;
  }

  // require valid PMT index for the MCMC
  if( !useTracker ) {
    if( target_pmt[2]<0 || target_pmt[2]>30 ) {
      cerr << "FATAL ERROR: invalid target pmt index (" 
	   << target_pmt[2] << ")" << endl ;
      exit(2) ;
    }
  }

  JModuleLocation target_location(target_pmt[0],target_pmt[1]) ;
  JModule target_module( getModule(detector,target_location) ) ;
  JAxis3D target_axis( getDOMcenter(target_module), target_module.getPMT(target_pmt[2]).getDirection() ) ;

  // source model
  JDirectedSource* src = new JDirectedSource( source_axis.getDirection(), 
					      opening_angle_deg*M_PI/180,
					      true,
					      surface_normal ) ;
  src->setPosition( source_axis.getPosition() ) ;

  // initialize scattering model
  cout << "PHYSICS" << endl ;
  JHenyeyGreensteinScattering sm1( lHG, gHG ) ; 
  cout << "Henyey-Greenstein scattering length = " << lHG << " [m], g = " << gHG << endl ;
  JRayleighScattering sm2( lR, aR ) ;
  cout << "Rayleigh scattering length = " << lR << " [m], a = " << aR << endl ;
  JScatteringModel* sm = new JCombinedScattering( &sm1, &sm2 ) ;
  cout << "Absorption length = " << Labs << " [m]" << endl ;

  // water properties
  // for now this is not customizable
  const double water_depth = 3000 ; // m
  const double pressure    = 0.1*water_depth ; // atm
  JDispersion jd(pressure) ;
  const double nbwavelength = 470 ; // nm
  const double index_of_refraction = jd.getIndexOfRefractionGroup(nbwavelength) ;
  const double cw = JTOOLS::getSpeedOfLight() / index_of_refraction ;
  cout << "Assuming speed of light in water = " << cw << " [m/ns]" << endl ;
  cout << endl ;


  // storage container for generated photon paths
  vector<JPhotonPath> paths ;

  // target model (for now this is not customizable)
  // we use the approximation that a PMT is a spherical cap on a DOM, with a uniform
  // acceptance, i.e. the photon always causes a hit when it impacts the cap
  cout << "TARGET" << endl ;
  const double pmt_opening_angle_deg = 30 ;
  const double pmt_opening_angle_rad = pmt_opening_angle_deg*M_PI/180.0 ;
  const double pmt_ctmin = cos(0.5*pmt_opening_angle_rad) ;
  cout << "DOM: string " << target_pmt[0] 
       << ", floor " << target_pmt[1] 
       << ", @" << getDOMcenter(target_module) << endl
       << "Target radius = " << DOM_radius << " [m]" << endl 
       << "PMT: " ;
  if( useTracker ) cout << "all 31 PMTs" ;
  else             cout << "PMT " << target_pmt[2] << ", direction = " << target_axis.getDirection() ;
  cout << endl
       << "PMTs are simulated as spherical caps" << endl
       << "PMT opening angle = " << pmt_opening_angle_deg << " deg." << endl
       << endl ;

  JTargetModel* trg ;
  if( useTracker ) {
    // for the tracker we use an isotropic sphere that represents an entire DOM
    // later on, we will divide the hits over the different PMTs
    trg = new JIsotropicTarget() ; 
    trg->setRadius(DOM_radius) ;
    trg->setPosition(target_axis.getPosition()) ;
  } else {
    // for the MCMC we use a single PMT target
    trg = new JPMTTarget(target_axis,pmt_opening_angle_rad,DOM_radius) ;
  }

  unsigned long int n_generated_photons = nphotons ;
  if( useTracker ) { // generate paths by photon tracking
    //n_generated_photons = 1000e6 ;
    cout << "Generating photon paths by photon tracking" << endl ;
    JMarkovPhotonTracker mpt ;
    paths = mpt.trackPhotons(
	 n_generated_photons,
	 src, sm, trg,
	 Labs ) ;
    cout << "Photon hits: " << paths.size() << " / " << n_generated_photons << endl ;
    cout << endl ;
  }

  vector<double> Pscat(nscatmax+1,1) ;
  vector<double> Perr(nscatmax+1,1) ;
  vector<JMarkovEnsembleIntegrator1D*> integrators ;
  if( !useTracker ) { // generate paths by MCMC
    JMarkovPathGenerator mpg ;
    if( no_remapping ) {
      // deactivate coordinate remapping.
      // That will most likely lead to trouble with the 1/r^2 singularities 
      // in the path probability density, meaning the space will be undersampled
      // for cases where vertices are very close together.
      mpg.setCoordinateRemapping(false) ;
      cout << "WARNING: coordinate remapping is deactivated!" << endl ;
    }
    mpg.setTargetStepSize_deg(target_step_size_deg) ;
    mpg.setTangentialStepSize_deg(tangential_step_size_deg) ;
    int nsteps_save = 250 ;
    int nsteps_burn_in = 10*nsteps_save ;
    for( int nscat=0 ; nscat<=(int)nscatmax ; ++nscat ) {
      cout << "nscat = " << nscat << endl ;

      double stepsize = 1 ; // default step size
      if( nscat>0 && (int)stepsizes.size()>=nscat ) stepsize = stepsizes[nscat-1] ;
      mpg.setRadialStepSize_m(stepsize) ;

      if( nscat != 0 ) cout << "step size = " << stepsize << " [m]" << endl ;

       // create start_path
      FreePathSolver fps ;
      JPhotonPath start_path = fps.solve(nscat,src,(JPMTTarget*)trg) ;
      cout << "path = " << endl ;
      cout << start_path.getString() << endl ;
      double start_rho = getPhotonPathProbabilityDensity(start_path,src,sm,trg,Labs) ;
      cout << "start_rho = " << start_rho << endl ;
      if( start_rho == 0 ) {
	cout << "No possible photon path found with nscat = " << nscat <<". Skipping nscat = " << nscat << "." << endl 
	     << "WARNING: note that this does not necessarily mean there are no such paths, simply that the search algorithm could not find one." << endl 
	     << endl;
	continue ;
      }

      /*cout << "Found the following start path: " << endl ;
      for( vector<JPosition3D>::iterator it=start_path.begin() ; it!= start_path.end() ; ++it ) {
	cout << *it << endl ;
      }
      cout << endl ;*/

      vector<JPhotonPath> these_paths = mpg.generateEnsemble( 
		   n_generated_photons,
		   start_path,
		   src, sm, trg,
		   Labs,
		   nsteps_burn_in, nsteps_save ) ;
      int nsteps = mpg.getNsteps() ;
      int naccepted = mpg.getNacceptedSteps() ;
      if( nscat != 0 ) {
	cout << "Accepted " << naccepted << " / " << nsteps << " steps (" << mpg.getFractionAccepted()*100.0 << "%)" << endl 
	     << "Accepted " << 100.0*mpg.getFractionAccepted_radial() << "% of the radial steps." << endl
	     << "Accepted " << 100.0*mpg.getFractionAccepted_tangential() << "% of the tangential steps." << endl ;
	cout << "NOTE: as a rule of thumb, the optimal value is around ~23% (this can be achieved by modifying the various step sizes)" << endl ;
      }
      paths.insert(paths.end(), these_paths.begin(), these_paths.end()) ;
      cout << endl ;

      if( these_paths.size() == 0 ) {
	cerr << "Something went wrong. The MCMC generator returned no paths." << endl ;
	Pscat[nscat] = 0.0 ;
	Perr[nscat]  = 0.0 ;
      } else {
	cout << "Creating JMarkovEnsembleIntegrator" << endl ;
	JMarkovEnsembleIntegrator1D* jmi = new JMarkovEnsembleIntegrator1D(these_paths,trg,100,100,200) ;

	cout << "Integrating" << endl ;
	vector<double> result = jmi->integrate(nsamples,nscat,src,sm,trg,Labs) ;
	printResult(result) ;
	Pscat[nscat] = getMean(result) ;
	Perr[nscat]  = getSigma(result) / sqrt(nsamples) ;
	//delete jmi ;
	integrators.push_back(jmi) ;
      }
      cout << endl ;
    }
  }

  if( paths.size()==0 ) {
    cerr << "FATAL ERROR: no paths were generated!" << endl ;
    exit(1) ;
  }


  // get path lengths and directions into an array
  // also find what the maximal number of scatterings is
  unsigned int npaths = paths.size() ;
  vector<double> path_lengths(npaths) ;
  vector<double> xs(npaths) ;
  vector<double> ys(npaths) ;
  vector<double> zs(npaths) ;
  for( unsigned int i=0; i<npaths; ++i ) {
    path_lengths[i] = paths[i].getLength() ;
    xs[i] = paths[i].back().getX() ;
    ys[i] = paths[i].back().getY() ;
    zs[i] = paths[i].back().getZ() ;
    int nscat = paths[i].size()-2 ;
    if( nscat>(int)nscatmax ) nscatmax = nscat ;
  }

  // find min and max x, y and z
  double xmin = trg->getPosition().getX() - 1.1*DOM_radius ; // *min_element(xs.begin(),xs.end() ) ;
  double xmax = trg->getPosition().getX() + 1.1*DOM_radius ; //*max_element(xs.begin(),xs.end() ) ;
  double ymin = trg->getPosition().getY() - 1.1*DOM_radius ; //*min_element(ys.begin(),ys.end() ) ;
  double ymax = trg->getPosition().getY() + 1.1*DOM_radius ; //*max_element(ys.begin(),ys.end() ) ;
  double zmin = trg->getPosition().getZ() - 1.1*DOM_radius ; //*min_element(zs.begin(),zs.end() ) ;
  double zmax = trg->getPosition().getZ() + 1.1*DOM_radius ; //*max_element(zs.begin(),zs.end() ) ;

  // find out which PMT is hit
  vector<int> PMT_is_hit ;

  if( !useTracker ) {
    // in the MCMC we know which PMT is hit
    PMT_is_hit = vector<int>(npaths,target_pmt[2]) ;
  }

  if( useTracker ) {
    PMT_is_hit = vector<int>(npaths,-1) ;
    vector<unsigned int> nHitsPerPMT(31,0) ;
    // in the tracker we have to figure out for every hit whether the PMT was hit or not
    for( unsigned int i=0; i<npaths; ++i ) {

      // find out which PMTs were hit (should be only one)
      vector<int> hitPMTs ;
      JVersor3D hit_direction = JVersor3D( paths[i].back() - target_axis.getPosition() ) ;
      for( unsigned int pmt=0 ; pmt<31 ; ++pmt ) {
	double ct = target_module.getPMT(pmt).getDirection().getDot(hit_direction) ;
	if( ct>pmt_ctmin ) hitPMTs.push_back(pmt) ;
      }
      if( hitPMTs.size()==1 ) {
	PMT_is_hit[i] = hitPMTs[0] ;
	++nHitsPerPMT[ hitPMTs[0] ] ;
      }
      if( hitPMTs.size()>1 ) {
	cerr << "FATAL ERROR: hit " << hitPMTs.size() << " PMTs simultaneously" << endl ;
	exit(1) ;
      }
    }
    // print number of hits separated by PMT
    for( int pmt=0 ; pmt<31 ; ++pmt ) {
      if( nHitsPerPMT[pmt]>0 ) {
	cout << "PMT" << setw(2) << pmt << ": " << setw(15) << nHitsPerPMT[pmt] << " hits" << endl ;
      }
    }
    cout << endl ;
  }

  // range and binning for pulse profile(s)
  const double dt = 0.1 ; // bin width in ns 
  const unsigned int margin = 5 ; // ns
  double Lmin = *min_element(path_lengths.begin(),path_lengths.end() ) ;
  double Lmax = *max_element(path_lengths.begin(),path_lengths.end() ) ;
  double tmin = floor( Lmin / cw - margin ) ;
  double tmax = ceil( Lmax / cw + margin )  ;
  int nbinst = (int) round( (tmax-tmin)/dt ) ;

  int pmtfront = 0  ;
  int pmtback  = 31 ;
  if( !useTracker ) {
    // in the MCMC there is only one PMT
    pmtfront = target_pmt[2]   ;
    pmtback  = target_pmt[2]+1 ;
  }

  for( int pmt=pmtfront ; pmt<pmtback ; ++pmt ) {
    // open file
    char fname[500] ;
    sprintf( fname, "%s_PMT%i.root", ofname.c_str(), pmt ) ;
    TFile* f = new TFile(fname,"recreate") ;

    // create histograms
    TH1D hpulse("hpulse","Instantaneous flash pulse profile;time since flash [ns];hit probability / ns",
		nbinst,tmin,tmax) ;
    vector<TH1D*> hpulse_per_nscat(nscatmax+1) ;
    for( int nscat=0 ; nscat<=(int)nscatmax ; ++nscat ) {
      char hname[200] ;
      sprintf(hname,"hpulse_nscat%i",nscat) ;
      hpulse_per_nscat[nscat] = new TH1D(hname,"Instantaneous flash pulse profile;time since flash [ns];hit probability / ns",
					 nbinst,tmin,tmax ) ;
    }
    TH2D hXY("hXY","hit position;X;Y",
	     100,xmin,xmax,
	     100,ymin,ymax ) ;
    hXY.SetOption("colz") ;
    TH2D hXZ("hXZ","hit position;X;Z",
	     100,xmin,xmax,
	     100,zmin,zmax ) ;
    hXZ.SetOption("colz") ;
    TH2D hYZ("hYZ","hit position;Y;Z",
	     100,ymin,ymax,
	     100,zmin,zmax ) ;
    hYZ.SetOption("colz") ;
    TH1D hnscat("hnscat",";N_{scat};probability",nscatmax+2,-0.5,nscatmax+1.5) ;

    TH3D hregion("hregion","Photon paths;X;Y;Z",
		 100,
		 trg->getPosition().getX() - 2*DOM_radius,
		 trg->getPosition().getX() + 2*DOM_radius,
		 100,
		 trg->getPosition().getY() - 2*DOM_radius,
		 trg->getPosition().getY() + 2*DOM_radius,
		 100,
		 trg->getPosition().getZ() - 2*DOM_radius,
		 trg->getPosition().getZ() + 2*DOM_radius) ;

    // fill histograms
    vector<JPhotonPath> paths_to_write ;
    for( unsigned int i=0 ; i<npaths ; ++i ) {
      if( PMT_is_hit[i] != pmt ) continue ;
      
      paths_to_write.push_back(paths[i]) ;

      const double path_length = paths[i].getLength() ;
      double t = path_length/cw ;
      int nscat = paths[i].size()-2 ;

      double w = 1 ;
      if( useTracker ) w = 1.0 / n_generated_photons ;
      else             w = Pscat[nscat] / n_generated_photons ; 

      hpulse.Fill(t,w) ;
      hnscat.Fill(nscat,w) ;
      hpulse_per_nscat[nscat]->Fill(t,w) ;
      
      hXY.Fill(xs[i],ys[i],w) ;
      hXZ.Fill(xs[i],zs[i],w) ;
      hYZ.Fill(ys[i],zs[i],w) ;

      if( paths.size()<10000 ) {
	const double dl = 0.001 ; // 1 mm
	const double dw = w*dl/path_length ;
	double l_along_seg = 0.5*dl ;
	int nfilled = 0 ;
	for( int nv=0; nv<(int)paths[i].size()-1; ++nv ) {
	  double seg_length = paths[i][nv+1].getDistance(paths[i][nv]) ;
	  JVersor3D seg_dir( paths[i][nv+1] - paths[i][nv] ) ;
	  while( l_along_seg<seg_length ) {
	    JPosition3D pos( paths[i][nv] + l_along_seg*JVector3D(seg_dir) ) ;
	    hregion.Fill( pos.getX(), pos.getY(), pos.getZ(), dw ) ;
	    if( pos.getDistance( trg->getPosition() ) < trg->getRadius() ) {
	      cerr << "Woohoooo" << endl ;
	      exit(1) ;
	    }
	    ++nfilled ;
	    l_along_seg += dl ;
	  }
	  l_along_seg -= seg_length ;
	}
      }
    }
    // scale hregion histogram, so that SetShowProjection3D will no show an
    // empty histogram
    hregion.Scale(1e20) ;

    // write paths 
    if( write_paths ) {
      char path_fname[500] ;
      sprintf( path_fname, "%s_PMT%i.paths", ofname.c_str(), pmt ) ;
      cout << "Writing paths." << endl ;
      writePaths( path_fname, paths_to_write ) ;
      cout << endl ;
    }

    // for the tracker, our error is just the Poisson error
    if( useTracker ) {
      for( Int_t bin=1 ; bin<=hnscat.GetNbinsX() ; ++bin ) {
	double val = hnscat.GetBinContent(bin) ;
	double nentries = max(1.0, val*n_generated_photons ) ;
	hnscat.SetBinError(bin, sqrt(nentries)/n_generated_photons) ;
      }
    }

    // for the MCMC, our error is sort of the error on the integral calculation
    if( !useTracker ) {
      for( int nscat=0 ; nscat<=(int)nscatmax ; ++nscat ) {
	Int_t bin = hnscat.GetXaxis()->FindBin(nscat) ;
	double val = hnscat.GetBinContent(bin) ;
	double err = Perr[nscat]/Pscat[nscat]*val ;
	hnscat.SetBinError(bin,err) ;
      }
    }
    
    // scale the pulse profile histograms by the bin width
    hpulse.Scale(1.0/dt) ;
    for( int nscat=0 ; nscat<=(int)nscatmax ; ++nscat ) {
      hpulse_per_nscat[nscat]->Scale(1.0/dt) ;
    }
    
    // make cumulative pulse profile (so you can more easily see what
    // the total hit probability is)
    TH1D* hpulse_cum = (TH1D*) hpulse.Clone("hpulse_cum") ;
    hpulse_cum->GetYaxis()->SetTitle("cumulative hit probability") ;
    for( Int_t bin=1 ; bin<=hpulse_cum->GetNbinsX() ; ++bin ) {
      double val = hpulse_cum->GetBinContent(bin-1) ; 
      val += hpulse_cum->GetBinContent(bin) * hpulse_cum->GetXaxis()->GetBinWidth(bin) ;
      hpulse_cum->SetBinContent(bin,val) ;
    }
    
    f->Write() ;

    f->mkdir("EnsembleIntegratorHistograms")->cd() ;
    for( vector<JMarkovEnsembleIntegrator1D*>::iterator it=integrators.begin(); it!=integrators.end(); ++it ) {
      (*it)->writeHistograms() ;
    }
    f->cd() ;

    f->Close() ;
    //delete f ;
    cout << "Output in '" << fname << "'." << endl ;
  }
  cout << "Done!" << endl ;
  cout << endl ;
}


const JModule& getModule( const JDetector& detector,  const JModuleLocation& location ) {
  for( unsigned int i=0; i<detector.size(); ++i ) {
    if( detector[i].getLocation() == location ) {
      return detector[i] ;
    }
  }
  cerr << "FATAL ERROR: no DOM found with location " << location << endl ;
  exit(1) ;
}