/*
Born approximation (Tsai (RMP 1974) Eq. 2.7) using Mathmore. 

*/

#include "TMath.h"
#include "TH1.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TObject.h"
#include "TROOT.h"
#include "TLegend.h"
#include "TStopwatch.h"

//#include "TLabel.h"
#include "Math/Functor.h"
#include "Math/WrappedFunction.h"
#include "Math/IFunction.h"
#include "Math/Integrator.h"
#include "Math/GSLMCIntegrator.h"
#include "Math/WrappedParamFunction.h"
#include <iostream>
	
const double pi = TMath::Pi();
const double zee = 4.;
const double ay = 9.;
const double mnuc = ay*0.938;
	
const double dee = (0.162)*pow(ay, -2./3.); //GeV2
const double mu = 0.105;
const double me = 0.511E-3; //electron mass in GeV
const double alpha = 1./(137.036);
const double zeered = pow(zee*alpha, 2.);

// integration order: cosplus, costheta, pee
// parameter order: kay, pee, costheta


double crossdiff(const double *x, const double *params){

//	1st integral is over costheta+  = x	
// 	2nd integral is over 2pi*costheta (emission solid angle)
//	3rd integral is over p, muon momentum

// only do coherent scattering, so no change in final nuclear mass
	
	double cosplus = x[0]; //costheta+
	double phi = x[1]; //costheta
	double xi = x[2]; //pee

//  parameters are: k, p, costheta (can define lab frame angles so no phi-dep?)

	double kay = params[0];
	
	double kthresh = (2.*mu*mu+2.*mu*mnuc)/mnuc;
	if (kay<kthresh) return 0;

//-------pee

	double pfunc1 = mnuc*kay - mnuc*mu;
	double pfunc2 = mnuc*mnuc + 2.*kay*mnuc;
		
	double pmax = (kay*pfunc1 + (kay+mnuc)*pow(pfunc1*pfunc1 - mu*mu*pfunc2, 0.5))/pfunc2;
	double pmin = (kay*pfunc1 - (kay+mnuc)*pow(pfunc1*pfunc1 - mu*mu*pfunc2, 0.5))/pfunc2;
	
	if(pmin<0.) pmin = 0;
	
	double pee = xi*(pmax - pmin) + pmin;

	double mueng = pow(pee*pee+mu*mu,0.5);

//----------costheta

	double costhetamax = (1./(kay*pee))*(mu*mnuc - mnuc*(kay - mueng) + kay*mueng) ;

	if(costhetamax<= -1.) costhetamax = -1.;
	if(costhetamax>1.) costhetamax = 1.;

	double costheta = phi*(1.-costhetamax) +costhetamax;
	
	
	double kmin = mnuc*(mu+mueng)/(mnuc - mueng + pee*costheta);
	if (kmin<0.) return 0;
	if (kay<kmin) return 0;

//-----------function

	double dotprod = kay*mueng - kay*pee*costheta;

	double yoo =  mu*mu +mnuc*mnuc + 2*mnuc*(kay-mueng) - 2*dotprod;
		   yoo =  pow(yoo, 0.5);
		   
	double kay_s = (kay*mnuc - dotprod)/yoo;
	double ee_plus_s = (yoo*yoo + mu*mu - mnuc*mnuc)/(2.*yoo);
	double pee_plus_s = pow(ee_plus_s*ee_plus_s - mu*mu , 0.5);
	double pee_init_s = (mnuc/yoo)*pow(kay*kay + pee*pee - 2*pee*kay*costheta , 0.5);
	double ee_s = (dotprod - mu*mu + mueng*mnuc)/yoo;
	double costhetak = (kay_s - ee_s)/(pee_init_s) + (dotprod)/(kay_s*pee_init_s);
	double sinsqthetak = 1. - costhetak*costhetak;
	
	double qsquared = 2*mu*mu - 2*dotprod - 2*ee_plus_s*(kay_s - ee_s) + 2.*pee_plus_s*pee_init_s*cosplus;
	
	double doubleyoo = ee_plus_s - pee_plus_s*cosplus*costhetak;	
	double wai = mu*mu*sinsqthetak + pow(pee_plus_s*cosplus - ee_plus_s*costhetak , 2.);
	       wai = pow(wai, 0.5);
	
	//delta = 0 for coherent scattering
	
	double haitch = -1.*mu*mu*( 0.5*qsquared*(1. - 2.*mueng/mnuc) + 2*mueng*mueng );

	double bee1 = (-2./dotprod)*(  (mu*mu - 0.5*qsquared)*(2.*mueng*(mueng-kay) + 0.5*qsquared*( (kay-2.*mueng)/mnuc + 1.) ) - 0.5*qsquared*kay*kay );
	double bee2 = (qsquared/mnuc)*(mnuc + mueng - kay - 0.5*qsquared/mnuc) + dotprod;
	double bee = bee1 + bee2;
	
	double cee1 = -1.*(mu*mu/(dotprod*dotprod)) *(  2.*(kay- mueng + 0.5*qsquared/mnuc)*(kay-mueng) + 0.5*qsquared );
	double cee2 = qsquared*(1.-mueng/mnuc)/dotprod;
	double cee = cee1 + cee2;

	double bigdee = 1./dotprod;
	
	double form2 = 2*mnuc*zee*zee/pow(1. - qsquared/dee, 2.);

	
	double integrand1 = (-1./(2*pi))*pow(alpha, 3.)*pow(0.19732, 2.)*1.E-26*pee*pee*pee_plus_s/(yoo*kay*mueng*qsquared*qsquared);
	double integrand2 = form2*( (haitch*doubleyoo/(wai*wai*wai*kay_s*kay_s)) + (bee/(wai*kay_s)) + cee + bigdee*kay_s*doubleyoo);

	return integrand1*integrand2*(1.-costhetamax)*2.*pi*(pmax - pmin);  
}

double crossdiffint(double* x){

	double kay = x[0];
	
	double p[1] = {kay};

	double lo[3] = {-1., 0., 0.};
	double hi[3] = {1., 1., 1.};
	
//	k>~ 200 GeV, calls = 120000*kay ?
//  k = 100 GeV, calls = 10000*kay 
//  k = 40 GeV, calls = 10000*kay 
	
// answers are different running script first time	?????
	
	unsigned int calls = 10000*kay;
	
	ROOT::Math::WrappedParamFunction<> funptr(&crossdiff, 3, 1, p);
  	ROOT::Math::GSLMCIntegrator ig1(ROOT::Math::MCIntegration::kVEGAS, 1E-6, 1E-4, calls);
  	ig1.SetFunction(funptr);
	double tsai = ig1.Integral(lo, hi);

	return tsai;
}

void tsaimc2(){
	
 	double kay = 100.;
 
 	double parameters[1] = {kay};
 
 	double answer = crossdiffint(parameters);
 	
 	cout <<answer <<endl;	
}