#include "TObject.h"
#include "TH2D.h"
#include <cmath>
#include <iostream>
#include "TRandom.h"
#include "TMath.h"
#include <vector> 

#include "Vec.hh"
#include "Mat.hh"
#include "Trk.hh"
#include "foreach.hh"
#include "rootutil.hh"

#include "Evt.hh"
#include "EventFile.hh"


using std::vector;
inline TH1D convolute( TH1D& h1, TH1D& h2 )
{
  TH1D r( "conv","conv",
	  h1.GetNbinsX(),
	  h1.GetBinLowEdge(1)+ h2.GetBinLowEdge(1),
	  h1.GetBinLowEdge(h1.GetNbinsX()+1 )+ h2.GetBinLowEdge(h2.GetNbinsX()+1) );


  for ( int i =1 ; i <= h1.GetNbinsX(); i++ )
    {
      double a1 = h1.GetBinLowEdge(i);
      double b1 = h1.GetBinLowEdge(i+1);
      double w1 = h1.GetBinContent(i);
      
      for ( int j =1 ; j <= h2.GetNbinsX(); j++ )
	{
	  double a2 = h2.GetBinLowEdge(j);
	  double b2 = h2.GetBinLowEdge(j+1);
	  double w2 = h2.GetBinContent(j);

	  int bin0 = r.FindBin( a1+a2 );
	  int bin1 = r.FindBin( b1+b2 );

	  for ( int bin = bin0; bin<=bin1; bin++)
	    {
	      double f = 1.0;
	      if ( r.GetBinLowEdge(bin)   < a1+a2 ) f = 1-( a1+a2 - r.GetBinLowEdge(bin) ) / r.GetBinWidth(bin) ;
	      if ( r.GetBinLowEdge(bin+1) > b1+b2 ) f = ( b1+b1 - r.GetBinLowEdge(bin) ) / r.GetBinWidth(bin) ;
                
           //cout << a1+a2 << " to " << b1+b2 << " : " << bin << " " << r.GetBinLowEdge(bin) << " " <<  r.GetBinLowEdge(bin+1) << " " << f << endl;
           r.AddBinContent( bin, w1*w2 * f ) ;   
	    }
        //cout << "----" << endl;
	  
	}
    }
  return r;
}

inline double addlogs( double a, double b )
{
  return log10(pow(10,a)+pow(10,b));
}


inline TH2D convolute_log2d(TH1D& h1, TH1D& h2 )
{
  TH2D r( "conv2d","conv2d",
	  h1.GetNbinsX(),
	  h1.GetBinLowEdge(1),
	  h1.GetBinLowEdge(h1.GetNbinsX()+1 ),

	  h2.GetNbinsX(),
	  h2.GetBinLowEdge(1),
	  h2.GetBinLowEdge(h1.GetNbinsX()+1 )
	  );
  
  for ( int i =1 ; i <= h1.GetNbinsX(); i++ )
    {
      double a1 = h1.GetBinLowEdge(i);
      double b1 = h1.GetBinLowEdge(i+1);
      double w1 = h1.GetBinContent(i);
      
      for ( int j =1 ; j <= h2.GetNbinsX(); j++ )
	{
	  double a2 = h2.GetBinLowEdge(j);
	  double b2 = h2.GetBinLowEdge(j+1);
	  double w2 = h2.GetBinContent(j);

	  double A = addlogs(a1,a2);
	  double B = addlogs(b1,b2);
	  
	  TAxis& xax = *r.GetXaxis();
	  
	  int bin0 = xax.FindBin( A );
	  int bin1 = xax.FindBin( B );

	  for ( int xbin = bin0; xbin<=bin1; xbin++)
	    {
	      int bin = r.GetBin( xbin, i );
	      
	      double f = 1.0;
	      if ( xax.GetBinLowEdge(xbin)   < A ) f = 1-( A - xax.GetBinLowEdge(xbin) ) / xax.GetBinWidth(xbin) ;
	      if ( xax.GetBinLowEdge(xbin+1) > B ) f = ( B - xax.GetBinLowEdge(xbin) ) / xax.GetBinWidth(xbin) ;
	      r.AddBinContent( bin, w1*w2 * f ) ;   
	    }	  
	}
    }

  return r;
}
		

		
inline TH1D convolute_log( TH1D& h1, TH1D& h2 )
{
  TH1D r( "conv","conv",
	  h1.GetNbinsX(),
	  h1.GetBinLowEdge(1),
	  h1.GetBinLowEdge(h1.GetNbinsX()+1 ) );

  for ( int i =1 ; i <= h1.GetNbinsX(); i++ )
    {
      double a1 = h1.GetBinLowEdge(i);
      double b1 = h1.GetBinLowEdge(i+1);
      double w1 = h1.GetBinContent(i);
      
      for ( int j =1 ; j <= h2.GetNbinsX(); j++ )
	{
	  double a2 = h2.GetBinLowEdge(j);
	  double b2 = h2.GetBinLowEdge(j+1);
	  double w2 = h2.GetBinContent(j);

	  double A = addlogs(a1,a2);
	  double B = addlogs(b1,b2);
	  
	  int bin0 = r.FindBin( A );
	  int bin1 = r.FindBin( B );

	  for ( int bin = bin0; bin<=bin1; bin++)
	    {
	      double f = 1.0;
	      if ( r.GetBinLowEdge(bin)   < A ) f = 1-( A - r.GetBinLowEdge(bin) ) / r.GetBinWidth(bin) ;
	      if ( r.GetBinLowEdge(bin+1) > B ) f = ( B - r.GetBinLowEdge(bin) ) / r.GetBinWidth(bin) ;

	      r.AddBinContent( bin, w1*w2 * f ) ;   
	    }	  
	}
    }
  return r;
}





struct Cylinder
{
  double height;
  double radius;
  Vec center;
  
  double top()    { return center.z + height/2 ; }
  double bottom() { return center.z - height/2 ; }

  Cylinder() {}
  
  Cylinder( double height, double radius, Vec center ) :
    height(height), radius(radius), center(center) {}

  Cylinder( double height, double radius,
	    double center_x=0, double center_y=0, double center_z=0 ) :
    height(height), radius(radius), center( Vec( center_x, center_y, center_z) ) {}

  const char* __str__() const
  {
    return Form("Cylinder height:%4.3f radius:%4.3f center: %s",
		height, radius, center.__str__() );
  }

  
  /*! returns true if pos is in the cylinder */
  bool contains (const Vec& pos )
  {
    Vec v = pos - center;
    if ( abs(v.z) > height / 2 ) return false;
    return ( v.lenxy() <= radius );
  }

  /*! 
    dz = -1 means downgoing, dz=+1 means upgoing, dz=0 horizontal
  */

  double area( double dz ) // is this correct?
  {
    const double s = sqrt( 1 - dz*dz);
    
    return
      pi*radius*radius    * abs(dz) +
      2 * height * radius * s;
  }


  /*! generate a random point on the surface, where a 
      track with direction d will
      intersect with the surface */

  Vec random_on_surface( const Vec& d )
  {
    
    const double s = sqrt( 1 - d.z*d.z);
    
    const double atop  = pi*radius*radius * abs(d.z) ;
    const double aside = 2 * height * radius * s;
    
    Vec v;
    
    if ( gRandom->Rndm() < atop / (atop+aside ) ) // generate on the top
      {
	const double r = radius * sqrt( gRandom->Rndm() );
	const double z   = d.z < 0 ? top() : bottom();
	const double phi = 2 * pi * gRandom->Rndm();
	v.set( r*sin(phi), r*cos(phi), z );
      }
    else // on the side tube
      {
	const double z = height * (gRandom->Rndm()-0.5);
	const double y = radius * (gRandom->Rndm()-0.5);
	const double x = -sqrt( radius*radius - y*y );
	const double phidir = atan2( d.y, d.x );
	v.set(x,y,z);
	v.rotate_z( phidir );
      }
    
    v.x += center.x;
    v.y += center.y;
    return v;
  }
  

  vector<Vec> intersect( const Trk& trk )
  {
    return intersect( trk.pos, trk.dir );
  }
  
  vector<Vec> intersect( const Vec& pos, const Vec& dir )
  {
    // --- check the top-view : does the line intersect the circle? ---

    Vec l = center-pos;
    l.z = 0;
    
    const double lxy = l.lenxy();
    
    Vec dirxy( dir.x, dir.y, 0 );
    dirxy.normalize();
    const double dxy = dirxy.dot(l);
    const double kxy = sqrt( lxy*lxy - dxy*dxy ) ;

    if ( kxy > radius) return {};
    
    const double a   = sqrt( radius*radius - kxy*kxy );
    
    // intersection in the xy plane with the tube:
    // d1 < d2, which means the track reaches vd1 before it reaches vd2
    
    const double sin_theta = sqrt( 1 - dir.z*dir.z );
    const double d1 = ( dxy - a ) / sin_theta; // after d1 m, it hits the first side
    const double d2 = ( dxy + a ) / sin_theta; // after d2 m, it hits the second side

    const Vec vd1 = pos + d1 * dir;
    const Vec vd2 = pos + d2 * dir;
    
    // does it pass over or under the can?
    if ( vd1.z > top()    && vd2.z > top()    ) return {}; // (1)
    if ( vd1.z < bottom() && vd2.z < bottom() ) return {};

    vector<Vec> r;
    
    if ( vd1.z >= bottom() &&  vd1.z <= top() ) // it hits the tube
      {
	r = {vd1};
      }
    else if ( vd1.z > top() ) // it must be downgoing because (1) is false, and it must hit the top circle
      {
	const double h1 = ( pos.z - top() ) / dir.z;
	const Vec vh1 = pos - h1 * dir;
	r = {vh1};
	
      }
    else if ( vd1.z < bottom() ) // it must be upgoing and hit the bottom circle
      {
	const double h1 = (pos.z - bottom()) / dir.z;
	const Vec vh1 = pos - h1 * dir;
	r = { vh1 };
      }
    
    if ( vd2.z >= bottom() &&  vd2.z <= top() ) // it hits the tube
      {
	r.push_back( vd2 );
      }
    else if ( vd2.z > top() ) // it must be up because (1) is false, and it must hit the top circle
      {
	const double h1 = (pos.z - top()) / dir.z;
	const Vec vh1 = pos - h1 * dir;
	r.push_back(vh1);
      }
    else if ( vd2.z < bottom() ) // it must be downgoing and hit the bottom circle
      {
	const double h1 = (pos.z - bottom()) / dir.z;
	const Vec vh1 = pos - h1 * dir;
	r.push_back(vh1);
      }

    return r;
  }

  ClassDefNV(Cylinder, 1)
};
      
inline std::ostream& operator<<( std::ostream& out, const Cylinder& c )
{
  return out << c.__str__();
}


template <typename F>
inline double generate_random(const F& cdf_func,
		       double xmin, double xmax,
		       double precission = 1e-6)
{
  double delta = ( xmax-xmin ) / 4;
  double x     = xmin + 2 * delta;
  double r     = gRandom->Rndm();

  for ( ; delta > precission ; delta /= 2 )
    {
      const int sign = 2 * ( r > cdf_func( x ) ) -1;
      x += sign * delta;
    }

  return x;
}

inline int random_bin( TH1& h )
{
  double r = gRandom->Rndm();

  return 1+TMath::BinarySearch(h.GetNbinsX(),
			       h.GetIntegral(),
			       r );
}


inline double random_in_bin( TH1& h, int bin )
{
  return h.GetBinLowEdge( bin ) + gRandom->Rndm() * h.GetBinWidth( bin );
}


struct RandomFromHist
{
  TH1D h1;
  vector<TH1D> h2;
  
  RandomFromHist() {}
  
  RandomFromHist( TH2D& H )
  {
    h1 = *H.ProjectionX();

    h2.resize(1); // add one empty histogram so that we start at 1
    
    for (int b = 1; b <= h1.GetNbinsX(); b++ )
      {
	// slice at bin b
	TH1D* h = H.ProjectionY( Form("%s%d", H.GetName(), b ) , b, b );
	h2.push_back( *h );
      }
  }

  // this is what we need for ct vs multiplity
  void get_random( double& x, int& biny )
  {
    int b1 = random_bin( h1 );
    x      = random_in_bin( h1, b1 );
    biny   = random_bin( h2[b1] );
  }
 
};



struct MuonFlux : public TObject
{  
  // depth- and zenith angle dependence of flux
  double  K0a =  0.0072;
  double  K0b = -1.927;
  double  K1a = -0.581;
  double  K1b =  0.034;
  double ni0a = -0.0771;
  double ni0b =  0.524;
  double ni0c =  2.068;
  double ni1a =  0.03;
  double ni1b =  0.47;
 

  double dN_dt_domega(double depth, double zenith, int multiplicity ) const
  {
    const double ct = cos(zenith);
    
    const double k0 = pow(depth,K0b) * K0a;
    const double k1 = (K1a * depth + K1b );
    const double k  = k0 * ct * exp( k1/ct );
    
    const double ni0 = ni0a*depth*depth + ni0b*depth + ni0c;
    const double ni1 = ni1a*exp(ni1b*depth);
    
    const double ni  = ni0*exp(ni1/ct);
    
    return k / pow(multiplicity,ni);
  }

 
  //======================== single muon energy distribution ================

  double  g0 = -0.232;
  double  g1 =  3.961;
  double e0a =  0.0304;
  double e0b =  0.359;
  double e1a = -0.0077;
  double e1b =  0.659;
  double BETA = 0.42;
  
  double BetaChi(double depth, double zenith ) const
  {
    return BETA* depth / cos( zenith );
  }
 
  double Gamma_singlemu( double depth ) const
  {
    return g0*log(depth) + g1;
  }
  
  double Epsilon_singlemu(double depth, double zenith) const
  {
    const double e0 = e0a * exp(e0b*depth);
    const double e1 = e1a*depth + e1b;
    return e0/cos(zenith) + e1;
  }

  
  /*! Function of the form N * (A+ E )^{-gamma}, including CDF, and function 
      to get random numbers from it. It is unsed to describe the energy dependence
      for both single- and multi-muon events.
  */

  struct EDist
  {
    double A, gamma, N ;
    double bias;

    EDist( double A_, double gamma_) :
      A(A_)                                 ,
      gamma(gamma_)                         ,
      N( pow ( A, gamma_-1 ) * (gamma_-1) ) {}

    
    double PDF ( double E ) const
    {
      return N * pow( A + E, -gamma );
    }
    
    double CDF( double E ) const
    {
      return  1 + N * pow( A + E, 1-gamma) / ( 1-gamma );
    }

    double CDF_inv( double y ) const
    {
      return pow ( (1-gamma) * (y-1) / N , 1/(1-gamma) ) - A ;
    }

    double random( double rmin = 0, double rmax = 1) const
    {
      return CDF_inv( gRandom->Uniform(rmin,rmax) );
    }

    
    
  };


  /*! Function of the form N * R (R0+ R )^{-gamma}, including CDF, and function 
      to get random numbers from it. It is used to describe the distance of tracks
      to the bundle axis.
  */

  struct RDist
  {
    double R0, alpha, N ;

    RDist( double R0_, double alpha_) :
      R0(R0_)                               ,
      alpha(alpha_)
    {
      N = pow ( R0, alpha-1) / R0;
    }
      
    double PDF ( double R ) const
    {
      return N * R * pow( R0 + R, -alpha );
    }
    
    double CDF( double R ) const
    {
      return  1 - N * pow ( R0+R, 1-alpha) * ( (alpha-1) * R + R0);
    }

    double random() const
    {
      return generate_random( [&](double r) { return CDF(r);},
			      0, 10000 );
    }
  
  };


  
  EDist energy_distribution( double depth, double zenith,
			     int multiplicity, double radius = -1 )
  {
    double gamma, epsilon;

    if ( multiplicity == 1 )
      {
	gamma   = Gamma_singlemu( depth );
	epsilon = Epsilon_singlemu( depth, zenith );
      }
    else
      {
	if ( radius < 0 ) std::cout << " error " << std::endl;
	gamma   = Gamma_multimu(radius, depth, multiplicity );
	epsilon = Epsilon_multimu(radius, depth, zenith );
      }

    const double bx = BetaChi( depth, zenith );
    const double A  = ( 1 - exp(-bx) )  * epsilon;
    return EDist(A, gamma);
  }
  
  TH1D energy_distribution_histogram_log( double depth, double zenith,
					  int multiplicity )
  {
    TH1D r("h","h",100,0,5);

    for (int i=1; i <= r.GetNbinsX(); i++)
      {
	double e = pow(10, r.GetBinCenter(i) );
	double dp_de = energy_distribution( depth, zenith, multiplicity , 1 ).PDF( e/1000 );
	r.SetBinContent(i, e* dp_de );
      }
 
    return r;
  }
  TH1D energy_distribution_histogram( double depth, double zenith,
				      int multiplicity )
  {
    TH1D r("h","h",1000,0,10000);

    for (int i=1; i <= r.GetNbinsX(); i++)
      {
	double e = r.GetBinCenter(i);
	double dp_de = energy_distribution( depth, zenith, multiplicity , 1 ).PDF( e/1000 );
	r.SetBinContent(i, e* dp_de );
      }

    return r;
  }

  
  // ================== multi muon radial distribtion =================
  
  // latteral spread 
  double     RHO0a = - 1.786;
  double     RHO0b =  28.26;
  double      RHO1 = - 1.06;
  double    THETA0 =   1.3;
  double  Fpiccolo =  10.4;
  double   ALPHA0a = - 0.448;
  double   ALPHA0b =   4.969;
  double   ALPHA1a =   0.0194;
  double   ALPHA1b =   0.276;
  
  double Alpha(double depth, int multiplicity) const
  {
    int mult = ( multiplicity < 4 ) ? multiplicity : 4;
    double alpha0 = ALPHA0a*mult + ALPHA0b;
    double alpha1 = ALPHA1a*mult + ALPHA1b;
    return alpha0 * exp(alpha1*depth);
  }
  
  double AverageRadius(double depth, double zenith, int multiplicity) const
  {
    int mult = ( multiplicity < 4 ) ? multiplicity : 4;
    double rho0 = RHO0a * mult + RHO0b;
    double espo = exp((zenith - THETA0) * Fpiccolo) + 1;
    return rho0 * pow(depth,RHO1)  / espo ; 
  }
  
  double Rpeak(double depth, double zenith, int multiplicity ) const
  {
    double alpha = Alpha(depth, multiplicity );
    double R0 = (alpha - 3) * AverageRadius(zenith, depth, multiplicity ) * 0.5;
    return R0 / (alpha - 1);
  }

  double dNdR(double r,
	      double  depth, double zenith, int multiplicity) const
  {
    double alpha = Alpha(depth, multiplicity);
    double R0 = (alpha - 3) * AverageRadius(depth,zenith,multiplicity) * 0.5;
    double AlphaLessTwo = alpha - 2;
    double C = (alpha - 1) * AlphaLessTwo * pow(R0, AlphaLessTwo);
    double denom = pow((r + R0), alpha);
    return C * r / denom;
  }

  RDist radius_distribution( double  depth, double zenith, int multiplicity )
  {
    const double alpha = Alpha(depth, multiplicity);
    const double R0 = (alpha - 3) * AverageRadius(depth,zenith,multiplicity) * 0.5;

    return RDist( R0, alpha );
  }

  
   /*! integral of dNdR - thanks to wolfram */
  double Ncum_vs_R( double R,
		    double  depth, double zenith, int multiplicity) const
  {
    double alpha = Alpha(depth, multiplicity);
    double R0    = (alpha-3) * AverageRadius(depth,zenith,multiplicity) * 0.5;
    double C     = (alpha-1) * (alpha-2) * pow(R0, (alpha-2) );

    const double N = (alpha-2) * (alpha-1);
    return 1 - C * pow( R0+R, 1-alpha) * (R0+ (alpha-1)*R ) /N;
  }

  
  //========================= multi-muon energy distribution =================
  
  // multi-muon parametrisation
  double  A0 =  0.0033;
  double  A1 =  0.0079;
  double B0a =  0.0407;
  double B0b =  0.0283;
  double B1a = -0.312;
  double B1b =  6.124;
  double  Q0 =  0.0543;
  double  Q1 = -0.365;
  double C0a = -0.069;
  double C0b =  0.488;
  double  C1 = -0.117;
  double D0a = -0.398;
  double D0b =  3.955;
  double D1a =  0.012;
  double D1b = -0.35;

  
  double Gamma_multimu(double r,
		       double depth, int multiplicity)
  {
    int mult = ( multiplicity < 4 ) ? multiplicity : 4;
    double a  = A0 * depth + A1;
    double b0 = B0a * mult + B0b;
    double b1 = B1a * mult + B1b;
    double b  = b0 * depth + b1;
    double q  = Q0 * depth + Q1;
    double func_gamma = 1 - (exp(q*r)*.5);
    return a*r + b*func_gamma;
  }
  
  double Epsilon_multimu(double r,
			 double depth, double zenith )
  {
    double c0 = C0a*depth + C0b;
    double c  = c0 * exp(C1*r);
    
    double d0 = D0a*depth + D0b;
    double d1 = D1a*depth + D1b;
    double d  = d0 * pow(r, d1);
    return c*zenith + d;
  }
  
  void dum();  
  ClassDef(MuonFlux,1);
};


struct EnergySampler
{
  double depth;
  double zenith;
  int    multiplicity;
  
  TH2D h_r_vs_e1; // y:radius, x:energy of 1 muon.

  vector<TH1D> h_en; // energy of bundle;
  vector<TH2D> h_en_vs_e1; // x: En, y: E(n-1)

  MuonFlux F;

  EnergySampler() {}
  
  EnergySampler( double depth_,
		 double zenith_,
		 double multiplicity_ ):
    
    depth(depth_), zenith(zenith_), multiplicity( multiplicity_)
  {
    h_en.resize(100);
    h_en_vs_e1.resize(100);
		
    
    int ebins = 440;
    double emin  = -1;
    double emax  = 10;

    int rbins = 1000;
    double rmin = 0;
    double rmax = 300;

    h_r_vs_e1 = TH2D(" h_r_vs_e1"," h_r_vs_e1",
		     ebins, emin, emax,
		     rbins, rmin, rmax );
  
    for (int i = 1; i <= h_r_vs_e1.GetXaxis()->GetNbins(); i++ )
      for (int j = 1; j <= h_r_vs_e1.GetYaxis()->GetNbins(); j++ )
	{
	  double loge = h_r_vs_e1.GetXaxis()->GetBinCenter( i );
	  double e    = pow(10,loge);
	  double r    = h_r_vs_e1.GetYaxis()->GetBinCenter( j );
	  double w = F.energy_distribution( depth, zenith, multiplicity , r ).PDF( e/1000 );
	  w *= F.radius_distribution( depth, zenith, multiplicity ).PDF(r);
	  
	  int b = h_r_vs_e1.GetBin(i,j);
	  h_r_vs_e1.SetBinContent( b, w * e  );
	}
    
    //----------------------------------------
    h_en[1] = * h_r_vs_e1.ProjectionX();
    h_en[1].Scale( 1.0 / h_en[1].Integral() );
    
    int M = multiplicity;
    //if ( M == 4 ) M = 100; 
    
    for (int i=2; i<= M ; i++ )
      {
	// the following must produce a 2d histrogram
	h_en_vs_e1[i] = convolute_log2d( h_en[i-1], h_en[1] );
	h_en[i] = *h_en_vs_e1[i].ProjectionX();
	cout << " convoluted " << i << endl;
      } 
  }

  vector<double> generate()
  {
    vector<double> r;
    
    double log_e = h_en[multiplicity].GetRandom();

    int M = multiplicity;

    for (int i=M; i>0 ; i--)
      {
	// get slice of h_en_vs_e1[i] at log_etot
	// ..and then throw random log_enext

	// r.push_back( pow(10, log_enext) - pow(10, log_e ) );
	
	// log_e = log_enext;
      }
      


  }


  
  void dus();
  
  ClassDefNV( EnergySampler,1 );
};


extern "C" {
  void initialize_music_sw_(int& i);
  
  void muon_transport_sw_(const double& x,   const double& y,      const double& z,
			  const double& dx,  const double& dy,     const double& dz,
			  const double& e,   const double& depth,  const double& time,
			  const double& rho,
			  const int& multiple_scattering_flag,
			  const int& deflection_flag,
			  double& fx,  double& fy,     double& fz,
			  double& fdx, double& fdy,    double& fdz,
			  double& fe,  double& fdepth, double& ftime );
  
  // x,y,z are in cm
  // energy is gev
  // time is in ns
  
  void rluxgo_(int&, int&, int&, int& );
  void rmarin_(int&, int&, int& );
}

#include "Trk.hh"



/*! using the Music muon propagator, let the track 
    fly from 'length' meters through water, 
    simulating energy loss and multiple scattering. */

inline Trk propagate_music( const Trk& trk , double length )
{
  if (!trk.is_muon())
    {
      print("warning: propagate_music called with non-muon track.");
    }

  
  static bool init = false;
  if (!init)
    {
      int een(1);
      initialize_music_sw_( een );
      init = true;
    }
  
  double rho = 1.03;
  int multiple_scattering = 1;
  int deflection          = 1;
  
  Trk trk_out ( trk );
  double outlen;
  
  muon_transport_sw_ ( trk.pos.x * 100 , trk.pos.y * 100, trk.pos.z * 100,
		       trk.dir.x, trk.dir.y, trk.dir.z,
		       trk.E,
		       length * 100 ,
		       trk.t,
		       rho, multiple_scattering, deflection,
		       trk_out.pos.x , trk_out.pos.y, trk_out.pos.z ,
		       trk_out.dir.x , trk_out.dir.y, trk_out.dir.z,
		       trk_out.E,
		       outlen,
		       trk_out.t );

  trk_out.pos /= 100;
  return trk_out;
}

inline void propagate_music( vector<Trk>& trks , double length )
{ 
  foreach( t, trks )
    {
      t = propagate_music( t, length );
    }
}


/*! holds all generation information, measuring generation efficiency, 
    debuggin and reweighting */

struct MupageEventInfo: public TObject
{
  double depth;
  int multiplicity;
  Vec axis_dir;
  Vec axis_pos;
  
  double depth_ratio;
  
  vector<double> radii;
  
  ClassDef( MupageEventInfo,1 );
};


struct Mupage : public TObject
{
  double depth            = 2.4; // in km
  double Emin             = 1e-3; // in TeV
  double Emax             = 1e3;
  double theta_min        = 0;
  double theta_max        = 85;
  double multiplicity_min = 1;
  double multiplicity_max = 100;

  Cylinder can;

  int _nbinsct = 100;
  double _ctmin;
  double _ctmax;

  double total_rate; // integral over multiplicity and zenith angle within the generation boundaries
  MuonFlux flux;
  TH2D hrate;
  RandomFromHist ct_mult_dist;

  MupageEventInfo info;


  bool bias_m1;
  double _fraction_above_emin_mult1;

  TH2D h_trials;
  TH2D h_success;
  
  
  Mupage() {}
  
  Mupage( double depth_, double can_height = 0, double can_radius = 1 ):
    depth(depth_)
  {
    can.height = can_height;
    can.radius = can_radius;
    can.center = Vec(0,0,-can_height/2);
    init();
  }
  
  void init()
  {
    _ctmin = cos( theta_max /180 * pi );
    _ctmax = cos( theta_min /180 * pi );
    
    hrate = TH2D ("hrate", ("hrate at depth "+str(depth)).c_str(),
		  _nbinsct, _ctmin, _ctmax,
		  multiplicity_max - multiplicity_min +1,
		  multiplicity_min - 0.5,
		  multiplicity_max + 0.5 );

    h_trials  = * (TH2D*) hrate.Clone("h_trials");
    h_success = * (TH2D*) hrate.Clone("h_success");
      
    
    total_rate = 0;
    for (int bx=1 ; bx<=hrate.GetNbinsX(); bx++ )
      for (int by=1 ; by<=hrate.GetNbinsY(); by++ )
	{
	  const int b   = hrate.GetBin(bx,by);
	  const double ct = hrate.GetXaxis()->GetBinCenter( bx );
	  const double zenith = acos( ct );

	  const double can_area = can.area( -ct );
	  
	  const double mult   = by; // in this case.	  
	  double f      = can_area * flux.dN_dt_domega( depth, zenith, mult );

	  // for the multiplicity == 1 bin, we immmediately implement the energy cut.
	  // So when we generate, we can simply pick E > Emin.
	  if ( mult == 1 && bias_m1 )
	    {
	      _fraction_above_emin_mult1 = 1 - flux.energy_distribution( depth, zenith, 1, 0.0 ).CDF( Emin );
	      
	      f *= _fraction_above_emin_mult1;
	    }
	  
	  
	  hrate.SetBinContent( b, f );

	  const double omega_bin = hrate.GetXaxis()->GetBinWidth( bx ) * 2 *pi;
	  total_rate += omega_bin * f;
	}

    // make the object to efficiently draw random numbers
    // from the ct_vs_mult distribution.
    
    ct_mult_dist = RandomFromHist( hrate );
  }

  /*! event rate per second */
  double event_rate()
  {
    return total_rate;
  }

  bool dummy_call()
  {
    return theta_max >0;
  }
    
  void generate( std::vector<Trk>& trks ) 
  {
    Timer _1("generate()");
        
    double ct, zenith;
    int multiplicity;

    ct_mult_dist.get_random( ct , multiplicity );
    zenith = acos( ct );

    
    // -------------------------------------------------------------
    //  generate impact point on the can
    // -------------------------------------------------------------
   
    info.axis_dir.set_angles( pi-acos( ct ), 2 * pi *  gRandom->Rndm() );
    info.axis_pos = can.random_on_surface( info.axis_dir );
    info.depth = depth + ( can.top() - info.axis_pos.z ) / 1000; // depth of actual impact on can

    
    // -------------------------------------------------------------
    //  generate the energetics
    // -------------------------------------------------------------

    double E[200];
    double r[200];
    
    if ( multiplicity == 1 && bias_m1 )
      {
	double f = flux.energy_distribution( info.depth, zenith, 1, 0.0 ).CDF( Emin );
	E[0] = 1000 * flux.energy_distribution( info.depth, zenith, 1 , 0 ).random( f );
      }
    else
      {
	auto rdist = flux.radius_distribution( info.depth, zenith, multiplicity );

	
	for (int i=0; i< multiplicity; i++)
	  {
	    r[i]= rdist.random();
	  }

	int ntrials = 10;
	double Etot;
	for (int n = 0; n< ntrials; n++ )
	  {
	    Etot = 0;
	    for (int i=0; i< multiplicity; i++)
	      {
		Etot += E[i]= flux.energy_distribution( info.depth, zenith, multiplicity, r[i] ).random();
	      }

	    if ( Etot > Emin ) break;
	  }
	
	h_trials.Fill( ct, multiplicity );
	if ( Etot < Emin )
	  {
	    trks.resize(0);
	    return;
	  }
	h_success.Fill( ct, multiplicity );
      }



    // ---------------------------------------------------------
    // Is the event still alive at the depth of impact ?
    // ---------------------------------------------------------
    
    if (info.depth == depth )
      {
	info.depth_ratio = 1;
      }
    else // we hit the side
      {  
	const double f0 = flux.dN_dt_domega( depth       , zenith , multiplicity );
	const double f1 = flux.dN_dt_domega( info.depth  , zenith , multiplicity );
    
	info.depth_ratio = f1/f0;
    
	if ( f1 > f0 )
	  {
	    cout << " programming error " << endl;
	    exit(1);
	  }

	if ( gRandom->Rndm() > info.depth_ratio )
	  {
	    trks.resize(0);
	    return;
	  }
      }


    
    // --------------------------------------------------------------------------
    //  Construct the bundle
    // --------------------------------------------------------------------------

    // the matrix R will rotate (0,0,1) to info.axis_dir
    Mat R( info.axis_dir );

    static Trk track;
    
    for (int i =0; i< multiplicity ; i++ )
      {
	track.dir = R * Vec(0,0,1);
	double phi = 2*pi* gRandom->Rndm();
	track.pos.set( r[i]*sin(phi), r[i]*cos(phi), 0 );
	track.E = E[i];
	
	// the shower axis is now exactly on the can, but others may miss the
	// can or reach it earlier or later.

	vector<Vec> v = can.intersect( track );
	if ( v.size() == 0 )  continue;

	double d = (v[0]-track.pos).len();
	track.pos += track.dir * d;
	track.t   += c_light   * d;

	track.id = trks.size();
	trks.push_back( track );
      }
  }


  
  void generate_event( Evt& evt )
  {
    static int id_counter = 0;
    evt.id = id_counter++;

    generate( evt.mc_trks );
  }

  int generate_event_file( string filename,
			    int nevents =10000 ) // before cuts
  { 
    EventFile f;
    f.set_output( filename );

    f.evt.mc_trks.reserve( 200 ) ;
    
    for (int i=0; i< nevents; i++)
      {
	generate_event( f.evt );
	if (f.evt.mc_trks.size() > 0 )
	  {
	    f.write();
	  }
	
      }
    
    return f.output_tree->GetEntries();
  }


  
  ClassDef(Mupage, 1);
};