#include "JppTrackPdf.hh"

#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>

#include "util.hh"
#include "Jeep/JMessage.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JPhysics/JPDFToolkit.hh"
#include "Vec.hh"
#include "Mat.hh"

#include "JGeometry3D/JPosition3D.hh"
#include "JGeometry3D/JRotation3D.hh"
#include "JGeometry3D/JDirection3D.hh"
#include "JGeometry3D/JAxis3D.hh"


#include "rootutil.hh"

using namespace JPP;


template<typename T>
void load_jpp_pdf( T& pdf , string fn )
{
  string f = fn;
  pdf.load(f.c_str() );
  pdf.compile();
  pdf.setExceptionHandler(new Function1D_t::JDefaultResult(JMATH::zero)); //??
  cout << "loaded " << f << endl;
}


JppTrackPdf::JppTrackPdf( string muon_direct_pdf,
			  string muon_scattered_pdf,
			  string em_direct_pdf,
			  string em_scattered_pdf,
			  string delta_direct_pdf,
			  string delta_scattered_pdf,
			  const double blur,
			  const double time_window_start,
			  const double time_window_end,
			  const double k40_rate,
			  bool   v ):
  _k40_rate( k40_rate ), verbose( v );
{ 
  pdfs.resize(6);

  load_jpp_track_pdf( pdfs[0], muon_direct_pdf     );
  load_jpp_track_pdf( pdfs[1], muon_scattered_pdf  );
  load_jpp_track_pdf( pdfs[2], em_direct_pdf       );
  load_jpp_track_pdf( pdfs[3], em_scattered_pdf    );
  load_jpp_track_pdf( pdfs[4], delta_direct_pdf    );
  load_jpp_track_pdf( pdfs[5], delta_scattered_pdf );
  
}


double JppTrackPdf::npe(double E,
			double D,
			double cos_pmt)
{
  const uint steps = 100;

  Vec dir = Vec(-1, 0, 0);

  const double pmt_a = acos(cos_pmt);

  std::vector<Vec> dirs;

  for(uint i = 0; i < steps; i++){

    const double phi = 2 * pi * i / steps;
    const Vec v = Vec(0, 0, 1).rotate_x(pmt_a).rotate_z(phi);
    dirs.push_back( v );
  }

  print(dirs);

  Mat Rot = Mat( dir );

  std::vector<Vec> rot_dirs;

  for(uint idir = 0; idir < dirs.size(); idir ++){
    rot_dirs[idir] = Rot * dirs[idir];
  }

  print(rot_dirs);

  return 1.0;
  
}


double JppTrackPdf::npe(double E,
			double R,
			double theta,
			double phi,
			double dt )
{
  const double deltas = getDeltaRaysFromMuon(E);
  
  terms = {
    pdfs[0](R, theta, phi, dt).V          ,
    pdfs[1](R, theta, phi, dt).V          ,
    pdfs[2](R, theta, phi, dt).V * E      ,
    pdfs[3](R, theta, phi, dt).V * E      ,
    pdfs[4](R, theta, phi, dt).V * deltas ,
    pdfs[5](R, theta, phi, dt).V * deltas};

  return sum(terms);
}

double JppTrackPdf::dnpe_dt(double E,
			    double R,
			    double theta,
			    double phi,
			    double dt)
{
  const double deltas = getDeltaRaysFromMuon(E);

  terms = {
    pdfs[0](R, theta, phi, dt).f          ,
    pdfs[1](R, theta, phi, dt).f          ,
    pdfs[2](R, theta, phi, dt).f * E      ,
    pdfs[3](R, theta, phi, dt).f * E      ,
    pdfs[4](R, theta, phi, dt).f * deltas ,
    pdfs[5](R, theta, phi, dt).f * deltas };

  return sum(terms);
}

double JppTrackPdf::d2npe_dt2(double E,
			      double R,
			      double theta,
			      double phi,
			      double dt)
{

  const double deltas = getDeltaRaysFromMuon(E);
  
  terms = {
    pdfs[0](R, theta, phi, dt).fp         ,
    pdfs[1](R, theta, phi, dt).fp         ,
    pdfs[2](R, theta, phi, dt).fp * E     ,
    pdfs[3](R, theta, phi, dt).fp * E     ,
    pdfs[4](R, theta, phi, dt).fp * deltas,
    pdfs[5](R, theta, phi, dt).fp * deltas};

  return sum(terms);

}



// double JppTrackPdf::npe( const Trk& trk,
// 			 const Hit& hit )
// {
//   HitParams p( trk, hit);
//   return npe( trk.E, p.R, p.theta, p.phi, p.dt );
// }


// double JppTrackPdf::dnpe_dt( const Trk& trk,
// 			     const Hit& hit )
// {
//   HitParams p( trk, hit);
//   return dnpe_dt( trk.E, p.R, p.theta, p.phi, p.dt );
// }


// double JppTrackPdf::d2npe_dt2( const Trk& trk,
// 			       const Hit& hit )
// {
//   HitParams p( trk, hit);
//   return d2npe_dt2( trk.E, p.R, p.theta, p.phi, p.dt );
// }



// double JppTrackPdf::eval_L(Trk& trk){
  
//   {
//     if (trk.E > 1e14 ) trk.E = 1e14;
    
//     //const double alpha = trk.E * 1e-6; // energy scale factor ( H3 is for 1 PeV ).
    
//     double L=0;
//     double L0=0;
//     double L1=0;
//     int nhit =0;
    
//     foreach( dom, domsums )
//       {
// 	{
// 	  TIMER( hit-pmt-loop )
	    
// 	    vector<Vec>::const_iterator ipmtdir = dom.hit_pmts.begin();
// 	  vector< vector<double> >::const_iterator ihit_times = dom.hit_times.begin();
	  
// 	  for(;
// 	      ipmtdir != dom.hit_pmts.end();
// 	      ipmtdir++, ihit_times++){
// 	    nhit ++;
// 	    Hit_params this_hit(trk, dom.pos, *ipmtdir);//much of this can be pulled out of loop
// 	    double mu1 = npe(trk.E, this_hit.D, this_hit.cd, this_hit.theta, this_hit.phi);
	    
// 	    //	    foreach(hit_time, *ihit_times){ //Need to take into account same pmt being hit several times
// 	    //Hit_params this_hit(trk, dom.pos, *ipmtdir, hit_time);//much of this can be pulled out of loop
// 	    //double mu1 = dnpe_dt(trk.E, this_hit.D, this_hit.cd, this_hit.theta, this_hit.phi, this_hit.dt);
	    
// 	    const double mu  = mu1 + k40prob;
// 	    double logp1 = log( 1 - exp( - mu ) );
	    
// 	    //if (do_clamp) clamp( logp1, minlogp, maxlogp );
// 	    L -= logp1;
// 	    L1 -= logp1;
// 	  }
// 	}
	
// 	{
// 	  TIMER( unhit-pmt-loop )
	    
// 	    for(vector<Vec>::const_iterator ipmtdir = dom.hit_pmts.begin();
// 		ipmtdir != dom.hit_pmts.end(); ipmtdir++){
// 	      Hit_params this_hit(trk, dom.pos, *ipmtdir);//much of this can be pulled out of loop
// 	      double mu1 = npe(trk.E, this_hit.D, this_hit.cd, this_hit.theta, this_hit.phi);

// 	      const double mu  = mu1 + k40prob;
// 	      double logp0 = -mu;
	      
// 	      //if (do_clamp) clamp( logp1, minlogp, maxlogp );
// 	      L -= logp0;
// 	      L0 -= logp0;
// 	    }
// 	}
	
// 	if ( L != L ) 
// 	  {
// 	    print("likelihood is nan @ r,z,a =");
// 	    print("trk Energy: ", trk.E);
// 	    print("Likelihood    hits:", L1);
// 	    print("Likelihood no hits:", L0);
// 	  }
	
//       } // dom
    
//     //test_params(trk); //A-ok
//     return L;
//   } 
// };


// void JppTrackPdf::test_params(Trk& trk){
  
//   {
//     foreach( dom, domsums )
//       {
// 	{
// 	  TIMER( hit-pmt-loop )
	    
// 	    vector<Vec>::const_iterator ipmtdir = dom.hit_pmts.begin();
// 	  vector< vector<double> >::const_iterator ihit_times = dom.hit_times.begin();
	  
// 	  for(;
// 	      ipmtdir != dom.hit_pmts.end();
// 	      ipmtdir++, ihit_times++){
// 	    foreach(hit_time, *ihit_times){ //Need to take into account same pmt being hit several times
	      
// 	      Hit_params this_hit(trk, dom.pos, *ipmtdir, hit_time);//much of this can be pulled out of loop
	      
// 	      if(round(this_hit.D * 1000) != round(calc_D(trk.pos, dom.pos) * 1000))
// 		print("test_params: bad D, ", this_hit.D , " vs. ", calc_D(trk.pos, dom.pos));
// 	      if(round(this_hit.R * 1000) != round(calc_R(trk.pos, dom.pos, trk.dir) * 1000))
// 		print("test_params: bad R, ", this_hit.R , " vs. ", calc_R(trk.pos, dom.pos, trk.dir));
// 	      if(round(this_hit.cd * 1000) != round(calc_cd(trk.pos, dom.pos, trk.dir) * 1000))
// 		print("test_params: bad cd, ", this_hit.cd , " vs. ", calc_cd(trk.pos, dom.pos, trk.dir));
// 	      if(round(this_hit.dt * 1000) != round(calc_dt(trk.pos, dom.pos, trk.t, hit_time) * 1000))
// 		print("test_params: bad dt, ", this_hit.dt , " vs. ", calc_dt(trk.pos, dom.pos, trk.t, hit_time));
// 	    }
// 	  }
// 	}
	
// 	{
// 	  TIMER( unhit-pmt-loop )
	    
// 	    for(vector<Vec>::const_iterator ipmtdir = dom.hit_pmts.begin();
// 		ipmtdir != dom.hit_pmts.end(); ipmtdir++){
// 	      Hit_params this_hit(trk, dom.pos, *ipmtdir);//much of this can be pulled out of loop
	      
// 	      if(round(this_hit.D * 1000) != round(calc_D(trk.pos, dom.pos) * 1000))
// 		print("test_params: bad D, ", this_hit.D , " vs. ", calc_D(trk.pos, dom.pos));
// 	      if(round(this_hit.R * 1000) != round(calc_R(trk.pos, dom.pos, trk.dir) * 1000))
// 		print("test_params: bad R, ", this_hit.R , " vs. ", calc_R(trk.pos, dom.pos, trk.dir));
// 	      if(round(this_hit.cd * 1000) != round(calc_cd(trk.pos, dom.pos, trk.dir) * 1000))
// 		print("test_params: bad cd, ", this_hit.cd , " vs. ", calc_cd(trk.pos, dom.pos, trk.dir));
// 	    }
// 	}
//       } // dom
//   } 
// };