#include <iostream>
using namespace std;

#include "TROOT.h"
#include "TFoam.h" 
#include "TH2D.h"
#include "TMath.h"
#include "TRandom3.h"
#include "TCanvas.h"
#include "TFoamIntegrand.h"
#include "TStyle.h"
#include "TFile.h"

double pairCS( double ePlus,
	       double theMinus, double phiMinus,
	       double thePlus , double phiPlus  ) ;

double testCS( double ePlus,
	       double theMinus, double phiMinus,
	       double thePlus , double phiPlus  ) ;

class TFDISTR: public TFoamIntegrand {
  //private:
public:
  static double f_thetaMin   , f_thetaMax;
  
  static double f_ePlusMin   , f_ePlusMax   ;
  static double f_theMinusMin, f_theMinusMax;
  static double f_phiMinusMin, f_phiMinusMax;
  static double f_thePlusMin , f_thePlusMax ;
  static double f_phiPlusMin , f_phiPlusMax ;

  TFDISTR() {};
  Double_t Density( int nDim, Double_t *Xarg ) {
    // Integrand for mFOAM
    // 5-dimensional distribution for FOAM
    
    Double_t ePlus     = ConvEPlus   ( Xarg[0] );
    Double_t theMinus  = ConvTheMinus( Xarg[1] );
    Double_t phiMinus  = ConvPhiMinus( Xarg[2] );
    Double_t thePlus   = ConvThePlus ( Xarg[3] );
    Double_t phiPlus   = ConvPhiPlus ( Xarg[4] );
    
    // cout << "requested " << ePlus << " " << theMinus << " " << phiMinus 
    // 	 << " " << thePlus << " " << phiPlus << endl;

    double val = pairCS( ePlus, theMinus, phiMinus, thePlus, phiPlus );
    //    double val = testCS( ePlus, theMinus, phiMinus, thePlus, phiPlus );
    
    // cout << "CS is " << val << endl;
    return val ;
  }

  static inline Double_t ConvEPlus   ( Double_t inVal ) 
  { return  f_ePlusMin    + inVal * ( f_ePlusMax    - f_ePlusMin   ) ; }
  static inline Double_t ConvTheMinus( Double_t inVal ) 
  { return f_theMinusMin + inVal * ( f_theMinusMax - f_theMinusMin ) ; }
  static inline Double_t ConvPhiMinus( Double_t inVal )
  { return f_phiMinusMin + inVal * ( f_phiMinusMax - f_phiMinusMin ) ; }
  static inline Double_t ConvThePlus( Double_t inVal ) 
  { return f_thePlusMin  + inVal * ( f_thePlusMax  - f_thePlusMin  ) ; }
  static inline Double_t ConvPhiPlus( Double_t inVal ) 
  { return f_phiPlusMin  + inVal * ( f_phiPlusMax  - f_phiPlusMin  ) ; }

  ClassDef(TFDISTR,1) //Class of testing functions for FOAM
};

ClassImp(TFDISTR)

double TFDISTR::f_thetaMin    = 1.0e-05;
double TFDISTR::f_thetaMax    = 5.0e-02;
double TFDISTR::f_ePlusMin    = 0.050 ;
double TFDISTR::f_ePlusMax    = 0.130;
double TFDISTR::f_theMinusMin = TFDISTR::f_thetaMin;
double TFDISTR::f_theMinusMax = TFDISTR::f_thetaMax;
double TFDISTR::f_phiMinusMin = 0.0;
double TFDISTR::f_phiMinusMax = TMath::TwoPi();
double TFDISTR::f_thePlusMin  = TFDISTR::f_thetaMin;
double TFDISTR::f_thePlusMax  = TFDISTR::f_thetaMax;
double TFDISTR::f_phiPlusMin  = 0.0;
double TFDISTR::f_phiPlusMax  = TMath::TwoPi();

double eGamma = 0;
double zMedia = 0;




int pairSim( double E = 0.110, double Z = 12.0, long nEvt = -1 ) {

  TFile rootFile( "foamOut.root", "recreate" );
  rootFile.cd();
  eGamma = E;
  zMedia = Z;

  TH2D  *hThPhiM  = new TH2D("hThPhiM" ,  "Theta-Phi plot M", 
			     180, 0., 0.05, 180, 0, TMath::TwoPi() );
  TH2D  *hThPhiP  = new TH2D("hThPhiP" ,  "Theta-Phi plot P",
			     180, 0., 0.05, 180, 0, TMath::TwoPi() );
  TH2D  *hPhiPsi = new TH2D("hPhiPsi",  "Phi-Psi plot", 
			    360, -TMath::TwoPi(), TMath::TwoPi(),
			    180, 0, TMath::TwoPi() );
  TH2D  *hThP    = new TH2D("hThPi" ,  "Theta-P plot",
			    180, 0., 0.05, 100, 0, 0.150 );
  
  Double_t *MCvect =new Double_t[5]; // 5-dim vector generated in the MC run

  Int_t  kDim     =       5;   // total dimension
  Int_t  nCells   = 5000000;   // Number of Cells
  Int_t  nSampl   =    5000;   // Number of MC events per cell in build-up
  Int_t  nBin     =     500;   // Number of bins in build-up
  Int_t  OptRej   =       1;   // Wted events for OptRej=0; wt=1 for OptRej=1 (default)
  Int_t  OptDrive =       2;   // (D=2) Option, type of Drive =0,1,2 for TrueVol,Sigma,WtMax
  Int_t  EvPerBin =     100;   // Maximum events (equiv.) per bin in buid-up
  Int_t  Chat     =       1;   // Chat level
//=========================================================
  TRandom *PseRan   = new TRandom3();       // Create random number generator
  TFoam   *FoamX    = new TFoam("FoamX");   // Create Simulator
  TFoamIntegrand    *rho= new TFDISTR(); 
  PseRan->SetSeed(4357);
//=========================================================
  cout <<"*****   Demonstration Program for Foam version " <<FoamX->GetVersion() 
       << "    *****" << endl;
  FoamX->SetkDim(        kDim);      // Mandatory!!!
  FoamX->SetnCells(      nCells);    // optional
  FoamX->SetnSampl(      nSampl);    // optional
  FoamX->SetnBin(        nBin);      // optional
  // FoamX->SetOptRej(      OptRej);    // optional
  // FoamX->SetOptDrive(    OptDrive);  // optional
  FoamX->SetEvPerBin(    EvPerBin);  // optional
  // FoamX->SetChat(        Chat);      // optiona

  FoamX->SetRho( rho );   // Set 5-dim distribution
  FoamX->SetPseRan(PseRan);   // Set random number generator
  cout << "Initializing " << endl;
  FoamX->Initialize();        // Initialize simulator, takes a few seconds...
  // From now on FoamX is ready to generate events according to Camel2(x,y)
  long nCalls=FoamX->GetnCalls();
  

  long nEvtReq = 50000;
  if( nEvt > 0 ) nEvtReq = nEvt;
  for( long iEvt = 0; iEvt < nEvtReq; iEvt++ ) {
    if( iEvt % 10000 == 0 && iEvt > 0 ) {
      cout << "Generating event " << iEvt << endl;
    }
    FoamX->MakeEvent();          // generate MC event
    FoamX->GetMCvect( MCvect);   // get generated vector (x,y)
    Double_t pPlus    = TFDISTR::ConvEPlus   ( MCvect[0] ) ;
    Double_t theMinus = TFDISTR::ConvTheMinus( MCvect[1] ) ;
    Double_t phiMinus = TFDISTR::ConvPhiMinus( MCvect[2] ) ;
    Double_t thePlus  = TFDISTR::ConvThePlus ( MCvect[3] ) ;
    Double_t phiPlus  = TFDISTR::ConvPhiPlus ( MCvect[4] ) ;
    //    cout << "Generated theta = " << theMinus << " phi " << phiMinus << endl;
    double phi = phiMinus - phiPlus;
    double psi = phiPlus;

    hThP->Fill   ( thePlus , pPlus ) ;
    hThPhiP->Fill( thePlus , phiPlus ) ;
    hThPhiM->Fill( theMinus, phiMinus ) ;
    hPhiPsi->Fill( phi, psi );
  }
  Double_t mcResult, mcError;
  Double_t eps = 0.0005;
  Double_t Effic, WtMax, AveWt, Sigma;
  Double_t IntNorm, Errel;
  FoamX->Finalize( IntNorm, Errel);     // final printout
  FoamX->GetIntegMC( mcResult, mcError);  // get MC intnegral
  FoamX->GetWtParams(eps, AveWt, WtMax, Sigma); // get MC wt parameters
  Effic=0; if(WtMax>0) Effic=AveWt/WtMax;
  cout << "================================================================" << endl;
  cout << " mcResult= " << mcResult << " +- " << mcError << " RelErr= "<< mcError/mcResult << endl;
  cout << " Dispersion/<wt>= " << Sigma/AveWt << endl;
  cout << "      <wt>/WtMax= " << Effic <<",    for epsilon = "<<eps << endl;
  cout << " nCalls (initialization only) =   " << nCalls << endl;
  cout << "================================================================" << endl;

  TCanvas *cPairs = new TCanvas("Canvas","Canvas for plotting",1300,1000);
  gStyle->SetPalette(1);
  cPairs->Divide( 2, 2 ) ;
  cPairs->cd(1);
  hThPhiM->DrawCopy("colz");
  cPairs->cd(2);
  hThPhiP->DrawCopy("colz");
  cPairs->cd(3);
  hPhiPsi->DrawCopy("colz");
  cPairs->cd(4);
  hThP->DrawCopy("colz");


  FoamX->Write( "FoamX" );

  hThPhiM->Write();
  hThPhiP->Write();
  hPhiPsi->Write();
  hThP->Write();

  rootFile.Write();
  rootFile.Close();

 return 0;
}




const double alpha     = 1.0/137.035;       // Fine structure constant
const double r0_fm     = 2.8179 ;           // electron classical radius in fm
const double fm2GeV_1  = 1.0/0.1973269;     // fermi to GeV^-1 conversion constant
const double r0        = r0_fm * fm2GeV_1;  // electron classical radius in GeV^-1
const double mEl       = 0.510999e-03 ;     // electron mass in GeV
const double pi        = 2.0 * acos(0.);    // number pi


// A function which calculates 5-D pair production diff. cross section
// d3 sigma / dE+ dTheta+ dTheta-
// See Depaola, NIM A452 (2000) 298; 
// Depaola, Astroparticle Physics 10 (1999) 175
double pairCS( double ePlus,
	       double theMinus, double phiMinus,
	       double thePlus , double phiPlus  ) {
 
  if( ePlus >= eGamma ) return 0;
  
  double phi = phiMinus - phiPlus;
  double psi = phiPlus;

  double q2 = ( -2.0 *
		( ePlus * ( eGamma - ePlus ) *
		  ( 1.0 - sin(thePlus) * sin(theMinus) * cos(phi) -
		    cos(thePlus) * cos(theMinus) ) +
		  eGamma * ePlus * ( cos(thePlus) - 1.0 ) +
		  eGamma * ( eGamma - ePlus ) * ( cos(theMinus) - 1.0 )
		  + mEl*mEl ) );

  if( q2 <= 0 ) {
    // cout << "Negative Q2 " << q2 << " "  ;
    // cout  << ePlus << " at " << theMinus << " " << phiMinus
    // 	  << " " << thePlus << " " << phiPlus << endl;
    return 0.0;
  }

  double Mult1 = -2.0 * alpha * zMedia*zMedia * r0 * mEl*mEl *
    ePlus * ( eGamma - ePlus ) /
    ( pow( 2*pi, 2 ) * pow( eGamma, 3 ) * ( q2*q2 ) );

  double Mult2 =  ePlus * ( eGamma - ePlus ) / ( q2*q2 ) ;

  double Mult3 = ( eGamma * sin(theMinus) / ( 1.0 - cos(theMinus) ) ) *
    ( eGamma * sin(thePlus) / ( 1.0 - cos(thePlus) ) ) ;

  double Brack1 = ePlus * sin(theMinus) * cos(psi) / ( 1.0 - cos(theMinus) ) + 
    ( eGamma - ePlus ) * sin(thePlus) * cos(psi+phi) / ( 1.0 - cos(thePlus) );

  double Brack2 = sin(theMinus) * cos(psi) / ( 1.0 - cos(theMinus ) ) - 
    sin(thePlus) * cos(psi+phi) / ( 1.0 - cos(thePlus) ) ;

  double Brack3 = ePlus * sin(thePlus) / ( ( eGamma - ePlus ) * sin(theMinus) ) 
    + ( eGamma - ePlus ) * sin(theMinus ) / ( ePlus * sin(thePlus) ) + 
    2.0 * cos(phi) ; 

  // Mult1 = 5000;
  // Mult3 = 0;

  // Brack1 = pow( 0.5*cos(psi) + 0.8*cos(psi+phi), 2 );
  // Brack2 = 0;
  // Brack3 = 10 + 2.0 * cos(phi) ; 

  double diffCS = sin(thePlus) *  sin(theMinus) * Mult1 * 
    ( 4.0 * Brack1*Brack1 - q2 * Brack2*Brack2 - Mult3 * Brack3 ) ;

  // cout << "q2 is " << q2 << endl;
  // cout << "Mult1 is " << Mult1 << " Mult2 is " << Mult2 << " Mult 3 is " << Mult3 << endl;
  // cout << "Brack1  is " << Brack1 << " Brack2 is " << Brack2 << " Brack3 is " << Brack3 << endl;
  // cout << "q2 is " << q2 << " CS is " << diffCS << endl;
  // cout << endl;

  if ( diffCS < 0 ) {
    cout << endl;
    cout << "Negative CS " << diffCS << " "  ;
    cout  << ePlus << " at " << theMinus << " " << phiMinus 
	  << " " << thePlus << " " << phiPlus << endl;
    return 0.0;
  }

  return diffCS ;
}


double testCS( double ePlus,
	       double theMinus, double phiMinus,
	       double thePlus , double phiPlus  ) {
  return ( 1.0 + cos( 2 * phiMinus ) );

}