#include "TRandom.h"
#include "TRandom2.h"
#include "TBaseClass.h"
#include <vector>
#include <string>
#include <iostream>
#include <cmath>
#include "TMath.h"
#include "TApplication.h"
#include "TMatrix.h"
#include "TVector3.h"
#include "Vec.hh"
#include "ana/search/SearchEvent.hh"
#include "TVector3.h"
#include "TLegend.h"
#include "ana/search/DetResponse.hh"

using namespace std;

TApplication app("app",NULL,NULL);



double gaus2dlog(double *x, double *par){
  double xxx=pow(10,x[0]);
  double_t arg = (xxx-x[1])/par[2];
  double res=arg*arg*exp(-0.5*arg*arg);
  //return res/par[0];
  return par[0]*res;

}

TSpline3 to_spline( TH1D& h )
{
  vector<double> xs;
  vector<double> ys;
  
  for (int b = 1; b <= h.GetNbinsX() ; b++ )
    {
      double bc = h.GetBinCenter(b);
      if ( bc > log10(180) ) break;
      xs.push_back(bc);
      ys.push_back( h.GetBinContent( b));
    }
  
  xs.push_back( log10(180) );
  ys.push_back( 0 );
  
  TSpline3 spline( "s", &(xs[0]), &ys[0], xs.size() , "e1b1"  );
  spline.SetLineColor(4);
  spline.SetLineWidth(2);
  return spline;
}

double d_omega_d_loga ( double loga ) 
{
  double a = pow( 10, loga) * pi / 180;
  return sin(a)  * 2 * pi * log(10) * a;
}

double d_loga_d_omega ( double loga ) 
{
  return 1.0 / d_omega_d_loga ( loga );
}	

int main(){

  

  int ngen=int(1e6);
  TRandom2 rnd;
  rnd.SetSeed(0);
 
  double sigma_src= 1; //degree
  double sigma_det= 1; //degree


  // Detector resolution
  int nbinloga=100;
  double xmin=-3;
  double xmax=3; 
  TH1D hdPdloga("hdPdloga","dP/dloga",nbinloga,xmin,xmax);
  TH1D hdPdOmegaa("hdPdOmegaa","dP/dOmegaa",nbinloga,xmin,xmax);
 

  //Extension of the source
  int nbinlogPsi=nbinloga;
   TH1D hdPdlogPsi("hdPdlogPsi","dP/dlogPsi",nbinloga,xmin,xmax);
  TH1D hdPdOmegaPsi("hdPdOmegaPsi","dP/dOmegaPsi",nbinloga,xmin,xmax);


 //Detector resolution + Extension of the source


  // ----------------------------------------------------
  // RANDOM generation using histo
  // -------------------------------------------------
  Timer timer_random("t_random");
  timer_random.start();

  TH1D hdPdlogCsi_rnd("hdPdlogCsi_rnd","dP/dlogCsi rnd",nbinloga,xmin,xmax);
  TH1D hdPdOmegaCsi_rnd("hdPdOmegaCsi_rnd","dP/dOmegaCsi rnd",nbinloga,xmin,xmax);
  
  Double_t a,loga,Psi,logPsi,Csi,logCsi;
  Double_t Omegaa,OmegaPsi,OmegaCsi;
  Vec Vrec,Vgen,Vsrc;
  
  for(int i=0;i<ngen;i++){
    
    a=sigma_det/deg * sqrt( -2 * log(rnd.Rndm()) ); //rad
    Psi = sigma_src/deg * sqrt( -2 * log(rnd.Rndm()) ); //rad
    Vrec.set_angles( a, rnd.Rndm() * 2 * pi );
    Vsrc.set_angles( Psi, rnd.Rndm() * 2 * pi );
    Csi=angle_between(Vrec,Vsrc);
    hdPdloga.Fill(log10(a*deg));
    hdPdlogPsi.Fill(log10(Psi*deg));
    hdPdlogCsi_rnd.Fill(log10(Csi*deg)); 
  }
  
  hdPdloga.Scale(1/hdPdloga.Integral());
  hdPdlogPsi.Scale(1/hdPdlogPsi.Integral());
  hdPdlogCsi_rnd.Scale(1/hdPdlogCsi_rnd.Integral());
  
   for(int ia=1;ia<=nbinloga;ia++){
    
    loga=hdPdloga.GetBinCenter(ia);
    hdPdOmegaa.SetBinContent(ia,hdPdloga.GetBinContent(ia)*d_loga_d_omega( loga ));
    
    logPsi=hdPdlogPsi.GetBinCenter(ia);
    hdPdOmegaPsi.SetBinContent(ia,hdPdlogPsi.GetBinContent(ia)*d_loga_d_omega( logPsi ));
    
    logCsi=hdPdlogCsi_rnd.GetBinCenter(ia);
    hdPdOmegaCsi_rnd.SetBinContent(ia,hdPdlogCsi_rnd.GetBinContent(ia)*d_loga_d_omega( logCsi ));
  }
  
  hdPdOmegaa.Scale(1/hdPdOmegaa.Integral());
  hdPdOmegaPsi.Scale(1/hdPdOmegaPsi.Integral());
  hdPdOmegaCsi_rnd.Scale(1/hdPdOmegaCsi_rnd.Integral());

  timer_random.stop();
  timer_random.report();
  // ----------------------------------------------------
  // Convolution 
  // ----------------------------------------------------

  // ----------------------------------------------------
  //Convolution using histo
  // ----------------------------------------------------

  TH1D hdPdlogCsi_conv_fromhisto("hdPdlogCsi_conv_fromhisto","dP/dlogCsi conv from histo",nbinloga,xmin,xmax);
  TH1D hdPdOmegaCsi_conv_fromhisto("hdPdOmegaCsi_conv_fromhisto","dP/dOmegaCsi conv from histo",nbinloga,xmin,xmax);
  int iCsi,ilogCsi;
  int npts=100; 
  int npts_phi=100;

  for(int ia=0;ia<=nbinloga;ia++){
    loga=hdPdloga.GetBinCenter(ia); //log(dec)
    a=pow(10,loga); //deg
    Vrec.set_angles(a/deg,0);
    for(int ipsi=1;ipsi<=nbinlogPsi;ipsi++){
      logPsi=hdPdlogPsi.GetBinCenter(ipsi); //log(deg)
      Psi=pow(10,logPsi); //deg
      for (double Phi=0; Phi<=360.; Phi+=360./static_cast<double>(npts_phi)) { //let's add running through the angle between two vectors
      	Vsrc.set_angles(Psi/deg,Phi/deg);
	Csi=angle_between(Vrec,Vsrc);//rad
	logCsi=log10(Csi*deg);
	ilogCsi= hdPdlogCsi_conv_fromhisto.GetXaxis()->FindBin(logCsi);
	hdPdlogCsi_conv_fromhisto.AddBinContent(ilogCsi,hdPdloga.GetBinContent(ia)*hdPdlogPsi.GetBinContent(ipsi));
      }
    }
    
  }
  
  hdPdlogCsi_conv_fromhisto.Scale(1/hdPdlogCsi_conv_fromhisto.Integral());
  
  
  for(int ia=1;ia<=nbinloga;ia++){
    
    logCsi=hdPdlogCsi_conv_fromhisto.GetBinCenter(ia);
    hdPdOmegaCsi_conv_fromhisto.SetBinContent(ia,hdPdlogCsi_conv_fromhisto.GetBinContent(ia)*d_loga_d_omega( logCsi ));
    
    
    
  }
  
  hdPdOmegaCsi_conv_fromhisto.Scale(1/hdPdOmegaCsi_conv_fromhisto.Integral());


  // ----------------------------------------------------
  //Convolution using spline
  // ----------------------------------------------------
  
  TSpline3 sdPdloga=to_spline(hdPdloga);
  TSpline3 sdPdOmegaa=to_spline(hdPdOmegaa);
  TSpline3 sdPdlogPsi=to_spline(hdPdlogPsi);
  TSpline3 sdPdOmegaPsi=to_spline(hdPdOmegaPsi);
  

  TH1D hdPdlogCsi_conv_fromspline("hdPdlogCsi_conv_fromspline","dP/dlogCsi conv from spline",nbinloga,xmin,xmax);
  TH1D hdPdOmegaCsi_conv_fromspline("hdPdOmegaCsi_conv_fromspline","dPdOmegaCsi conv from spline",nbinloga,xmin,xmax);

  
  double logphimin=xmin;
  double logphimax=log(360)/log(10); //since the phi distribution is uniform, we need to cut it manually
  double Phi, logPhi;

  
  for(int ia=1;ia<=npts;ia++){
    loga=xmin+(xmax-xmin)/(double)npts*ia;
    a=pow(10,loga); //deg
    Vrec.set_angles(a/deg,0.);
    for(int ipsi=1;ipsi<=npts;ipsi++){
      logPsi=xmin+(xmax-xmin)/(double)npts*ipsi;
      Psi=pow(10,logPsi); //deg
      //for (double Phi=0; Phi<=360.; Phi+=360./static_cast<double>(npts_phi)) { //let's add running through the angle between two vectors
      for (double iphi=1; iphi<=npts;iphi++){
	logPhi = logphimin+(logphimax-logphimin)/(double)npts*iphi;
	Phi = pow(10,logPhi);
	Vsrc.set_angles(Psi/deg,Phi/deg);
	//Vsrc.set_angles(Psi/deg,0);
	Csi=angle_between(Vrec,Vsrc);//rad
      	logCsi=log10(Csi*deg);
        ilogCsi= hdPdlogCsi_conv_fromspline.GetXaxis()->FindBin(logCsi);
        hdPdlogCsi_conv_fromspline.AddBinContent(ilogCsi,sdPdloga.Eval(loga)*sdPdlogPsi.Eval(logPsi)*Phi);
	}
    }
  }

  hdPdlogCsi_conv_fromspline.Scale(1/hdPdlogCsi_conv_fromspline.Integral());


 for(int ia=1;ia<=nbinloga;ia++){
   
   logCsi=hdPdlogCsi_conv_fromspline.GetBinCenter(ia);
   hdPdOmegaCsi_conv_fromspline.SetBinContent(ia,hdPdlogCsi_conv_fromspline.GetBinContent(ia)*d_loga_d_omega( logCsi ));
 }
  
  hdPdOmegaCsi_conv_fromspline.Scale(1/hdPdOmegaCsi_conv_fromspline.Integral());


  // ----------------------------------------------------
  //Convolution using TF
  // ----------------------------------------------------
  Timer timer_conv("t_conv");
  timer_conv.start();
  
  TH1D hdPdlogCsi_conv_fromfunc("hdPdlogCsi_conv_fromfunc","dP/dlogCsi conv from function",nbinloga,xmin,xmax);
  TH1D hdPdOmegaCsi_conv_fromfunc("hdPdOmegaCsi_conv_fromfunc","dP/dOmegaCsi conv from function",nbinloga,xmin,xmax);


 TF1 fdPdloga("fdPdloga",gaus2dlog,xmin,xmax,3);
 fdPdloga.SetParameter(0,0.135); // to check
 fdPdloga.SetParameter(1,0.);
 fdPdloga.SetParameter(2,sigma_det);
 cout<<fdPdloga.Integral(xmin,xmax)<<endl;
 TF1 fdPdlogPsi("hdPdlogpsi",gaus2dlog,xmin,xmax,3);
 fdPdlogPsi.SetParameter(0,0.135);// to check
 fdPdlogPsi.SetParameter(1,0.);
 fdPdlogPsi.SetParameter(2,sigma_src);
 cout<<fdPdlogPsi.Integral(xmin,xmax)<<endl;



 int npoints=100;
 int npoints2=100;

 for(int ia=0;ia<=npoints;ia++){
   loga=xmin+(xmax-xmin)/(double)npoints*ia;
   a=pow(10,loga); //deg
   Vrec.set_angles(a/deg,0);
   for(int ipsi=1;ipsi<=npoints2;ipsi++){
     logPsi=xmin+(xmax-xmin)/(double)npoints2*ipsi;
     Psi=pow(10,logPsi); //deg
     for (double Phi=0; Phi<=360.; Phi+=360./static_cast<double>(npts_phi)) { //let's add running through the angle between two vectors
       Vsrc.set_angles(Psi/deg,Phi/deg);
       //Vsrc.set_angles(Psi/deg,0.);
       Csi=angle_between(Vrec,Vsrc);//rad
       if(Csi*deg>=1e-3 && Csi*deg<=1e3){
	 logCsi=log10(Csi*deg);
	 ilogCsi= hdPdlogCsi_conv_fromfunc.GetXaxis()->FindBin(logCsi);
	 hdPdlogCsi_conv_fromfunc.AddBinContent(ilogCsi,fdPdloga.Eval(loga)*fdPdlogPsi.Eval(logPsi));
	 
       }
     }
   }
   
 }
  
  hdPdlogCsi_conv_fromfunc.Scale(1/hdPdlogCsi_conv_fromfunc.Integral());


 for(int ia=1;ia<=nbinloga;ia++){
   
   logCsi=hdPdlogCsi_conv_fromfunc.GetBinCenter(ia);
   hdPdOmegaCsi_conv_fromfunc.SetBinContent(ia,hdPdlogCsi_conv_fromfunc.GetBinContent(ia)*d_loga_d_omega( logCsi ));
 }
  
  hdPdOmegaCsi_conv_fromfunc.Scale(1/hdPdOmegaCsi_conv_fromfunc.Integral());
     
   timer_conv.stop();
   timer_conv.report();
  // ------------------------------------------------------------------
  //   PLOT
  // -----------------------------------------------------------------


  TCanvas c1;
  c1.Divide(3,2);

  c1.cd(1);
  gPad->SetLogy();
  hdPdloga.GetXaxis()->SetTitle("loga");
  hdPdloga.Draw("HIST"); 
  sdPdloga.Draw("SAME");
  fdPdloga.SetLineColor(6);
  fdPdloga.Draw("SAME");
  auto legend1 = new TLegend(0.55,0.7,0.9,0.9);
  legend1->AddEntry(&hdPdloga,"dP/dloga rnd","l");
  legend1->AddEntry(&sdPdloga,"dP/dloga spline","l"); 
  legend1->AddEntry(&fdPdloga,"dP/dloga function","l");
  legend1->Draw();


  c1.cd(2);
  gPad->SetLogy();
  hdPdlogPsi.GetXaxis()->SetTitle("log#psi");
  hdPdlogPsi.Draw("HIST"); 
  sdPdlogPsi.Draw("SAME");
  fdPdlogPsi.SetLineColor(6);
  fdPdlogPsi.Draw("SAME");
  auto legend2 = new TLegend(0.55,0.7,0.9,0.9);
  legend2->AddEntry(&hdPdlogPsi,"dP/dlog#psi rnd","l");
  legend2->AddEntry(&sdPdlogPsi,"dP/dlog#psi spline","l"); 
  legend2->AddEntry(&fdPdlogPsi,"dP/dlog#psi function","l");
  legend2->Draw();

  c1.cd(3); 
  gPad->SetLogy();
  hdPdlogCsi_rnd.GetXaxis()->SetTitle("log#xi");
  hdPdlogCsi_rnd.Draw("HIST");

  hdPdlogCsi_conv_fromhisto.SetLineColor(1);
  hdPdlogCsi_conv_fromhisto.Draw("HISTSAME");
  
  hdPdlogCsi_conv_fromspline.SetLineColor(2);
  hdPdlogCsi_conv_fromspline.Draw("HISTSAME");

  hdPdlogCsi_conv_fromfunc.SetLineColor(4);
  hdPdlogCsi_conv_fromfunc.Draw("HISTSAME");
  
  auto legend3 = new TLegend(0.55,0.7,0.9,0.9);
  legend3->AddEntry(&hdPdlogCsi_rnd,"dP/dlog#xi rnd","l");
  legend3->AddEntry(&hdPdlogCsi_conv_fromhisto,"dP/dlog#xi conv from histo","l"); 
  legend3->AddEntry(&hdPdlogCsi_conv_fromspline,"dPd/log#xi conv from spline","l");
  legend3->AddEntry(&hdPdlogCsi_conv_fromfunc,"dP/dlog#xi conv from function","l");
  legend3->Draw();



  c1.cd(4); 
  gPad->SetLogy();
  hdPdOmegaa.GetXaxis()->SetTitle("loga");
  hdPdOmegaa.Draw("HIST"); 
  c1.cd(5);
  gPad->SetLogy();
  hdPdOmegaPsi.GetXaxis()->SetTitle("log#psi");
  hdPdOmegaPsi.Draw("HIST");
  c1.cd(6);
  gPad->SetLogy();
  hdPdOmegaCsi_rnd.GetXaxis()->SetTitle("log#xi");
  hdPdOmegaCsi_rnd.Draw("HIST");
  
  hdPdOmegaCsi_conv_fromhisto.SetLineColor(1);
  hdPdOmegaCsi_conv_fromhisto.Draw("HISTSAME");

  hdPdOmegaCsi_conv_fromspline.SetLineColor(2);
  hdPdOmegaCsi_conv_fromspline.Draw("HISTSAME");  

  hdPdOmegaCsi_conv_fromfunc.SetLineColor(4);
  hdPdOmegaCsi_conv_fromfunc.Draw("HISTSAME");

 
  auto legend4 = new TLegend(0.55,0.7,0.9,0.9);
  legend4->AddEntry(&hdPdOmegaCsi_rnd,"dP/d#Omega_{#xi} rnd","l");
  legend4->AddEntry(&hdPdOmegaCsi_conv_fromhisto,"dP/ conv from histo)","l");
  legend4->AddEntry(&hdPdOmegaCsi_conv_fromspline,"dP/d#Omega_{#xi} conv from spline","l");
  legend4->AddEntry(&hdPdOmegaCsi_conv_fromfunc,"dP/ conv from function","l");
  legend4->Draw();



  TCanvas c2;
  c2.Divide(2,1);
  c2.cd(1);
  gPad->SetLogy();
  hdPdloga.SetLineColor(1);
  hdPdloga.Draw("HIST");
  hdPdlogPsi.SetLineColor(2);
  hdPdlogPsi.Draw("HISTSAME");
  hdPdlogCsi_rnd.SetLineColor(3);
  hdPdlogCsi_rnd.Draw("HISTSAME");

  auto legend5 = new TLegend(0.55,0.7,0.9,0.9);
  legend5->AddEntry(&hdPdloga,"dP/dloga","l");
  legend5->AddEntry(&hdPdlogPsi,"dP/dlog#Psi)","l");
  legend5->AddEntry(&hdPdlogCsi_rnd,"dPdlog#xi","l");
  legend5->Draw();


  c2.cd(2);
  gPad->SetLogy();
  hdPdOmegaa.SetLineColor(1);
  hdPdOmegaa.Draw("HIST");
  hdPdOmegaPsi.SetLineColor(2);
  hdPdOmegaPsi.Draw("HISTSAME");
  hdPdOmegaCsi_rnd.SetLineColor(3);
  hdPdOmegaCsi_rnd.Draw("HISTSAME");

  auto legend6 = new TLegend(0.55,0.7,0.9,0.9);
  legend6->AddEntry(&hdPdOmegaa,"dP/d#Omega_{a}","l");
  legend6->AddEntry(&hdPdOmegaPsi,"dP/d#Omega_{#Psi})","l");
  legend6->AddEntry(&hdPdOmegaCsi_rnd,"dPd#Omega_{#xi}","l");
  legend6->Draw();

  TCanvas c3;
  c3.Divide(2,1);
  c3.cd(1);
  TH1D h_ratio_loga=hdPdlogCsi_rnd/hdPdlogCsi_conv_fromhisto;
  TH1D h_ratio_loga_func=hdPdlogCsi_rnd/hdPdlogCsi_conv_fromfunc;
  h_ratio_loga.Draw("HIST"); 
  h_ratio_loga_func.SetLineColor(2); 
  h_ratio_loga_func.Draw("HISTSAME");
  c3.cd(2);
  TH1D h_ratio_Omega=hdPdOmegaCsi_rnd/hdPdOmegaCsi_conv_fromhisto;
  TH1D h_ratio_Omega_func=hdPdOmegaCsi_rnd/hdPdOmegaCsi_conv_fromfunc;
  h_ratio_Omega.Draw("HIST"); 
  h_ratio_Omega_func.SetLineColor(2); 
  h_ratio_Omega_func.Draw("HISTSAME");
  app.Run(true);
  c1.Close();
  c2.Close();
}