/*
 *  SolidIBDKinematicsNew.C
 *  ANRA
 *
 *  Created by Gregory Lehaut on 26/11/14.
 *  Copyright 2014 LPC Caen, CNRS/IN2P3. All rights reserved.
 *
 */

#include "SolidIBDKinematics.hh"
#include <iostream>
#include <fstream>
#include <math.h>

#include "Randomize.hh"

#include "TStyle.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TMultiGraph.h"
#include "TCanvas.h"
#include "TFrame.h"
#include "TPad.h"
#include "TLegend.h"
#include "TH2D.h"

#include "globals.hh"


/////////////////////////////////////////////////////////////////////////////
/* Class implementation */

protonibd::protonibd(){
	_mnu = 0.000;
	_mp = 938.272;  //MeV
	_me = 0.51099;  //MeV
	_mn = 939.565;  //MeV
	_Tnu = 0.;      //MeV
	_Tp = 0.;       //
}

protonibd::~protonibd(){}

protonibd_Vogel0::protonibd_Vogel0()
:protonibd(){
	
	double Vud = 0.97425;
	double deltainner = 0.024;
	double Gf2 = 1.1663788e-11; //MeV^(-2)
	_sigma0 = Gf2*Vud*Vud/M_PI*(1.+deltainner);
	
	_delta = _mn - _mp;
	_f = 1.;
	_g = 1.26;
}

protonibd_Vogel0::~protonibd_Vogel0(){}


void protonibd::computeInputKinematic( ){
	_s = (_Tp+_Tnu+_mnu+_mp)*(_Tp+_Tnu+_mnu+_mp)
	-(sqrt((_Tnu+_mnu)*(_Tnu+_mnu)-_mnu*_mnu)+sqrt((_Tp+_mp)*(_Tp+_mp)-_mp*_mp))
	*(sqrt((_Tnu+_mnu)*(_Tnu+_mnu)-_mnu*_mnu)+sqrt((_Tp+_mp)*(_Tp+_mp)-_mp*_mp));
	if (sqrt(_s)-_me-_mn < 0) {
		_gammaCM = 0.;
		_betaCM = 0.;
		_pf = 0.;
		return;
	}
	_betaCM = (sqrt((_Tnu+_mnu)*(_Tnu+_mnu)-_mnu*_mnu)
			   +sqrt((_Tp+_mp)*(_Tp+_mp)-_mp*_mp))
	/(_Tp+_Tnu+_mnu+_mp);
	_gammaCM = 1./sqrt(1.-_betaCM*_betaCM);
	_pf = sqrt((sqrt(_s)-_me-_mn)*
			   (sqrt(_s)-_me+_mn)*
			   (sqrt(_s)+_me+_mn)*
			   (sqrt(_s)+_me-_mn))/(sqrt(_s)*2.);
}

double protonibd::getEmaxPositron( double enu){
	_Tnu = enu;
	computeInputKinematic();
	if (_betaCM!= 0.) {
		return( _gammaCM*sqrt(_pf*_pf+_me*_me)+_gammaCM*_betaCM*_pf - _me );	
	}
	else {
		return( 0. );
	}
}

double protonibd::getEminPositron( double enu){
	_Tnu = enu;
	computeInputKinematic();
	if (_betaCM!= 0.) {
		return( _gammaCM*sqrt(_pf*_pf+_me*_me)-_gammaCM*_betaCM*_pf - _me );	
	}
	else {
		return( 0. );
	}
	
}

double protonibd::getEmaxNeutron( double enu){
	_Tnu = enu;
	computeInputKinematic();
	if (_betaCM!= 0.) {
		return( _gammaCM*sqrt(_pf*_pf+_mn*_mn)+_gammaCM*_betaCM*_pf - _mn );	
	}
	else {
		return( 0. );
	}
}

double protonibd::getEminNeutron( double enu){
	_Tnu = enu;
	computeInputKinematic();
	if (_betaCM!= 0.) {
		return( _gammaCM*sqrt(_pf*_pf+_mn*_mn)-_gammaCM*_betaCM*_pf - _mn );	
	}
	else {
		return( 0. );
	}
	
}

double protonibd::getThreshold(){
    return( ( (_mn + _me)*(_mn + _me) - _mp*_mp ) / (2.*_mp) ); // in the Lab frame
}


double protonibd::getCX( double enu ){
	_Tnu = enu;
	computeInputKinematic();
	if (_betaCM != 0.) {
		return(1.);
	}
	else {
		return(0.);
	}
}

double protonibd_Vogel0::getCX( double enu ){
	_Tnu = enu;
	computeInputKinematic();
	if (_betaCM != 0.) {
		_Te = enu-_delta;
		return( _sigma0*(_f*_f+3.*_g*_g)*_Te*sqrt(_Te*_Te-_me*_me) );
	}
	else {
		return(0.);
	}
}

void protonibd::genereEvent( double enu ){
	_Tnu = enu;
	computeInputKinematic();
	if (_betaCM != 0.) {
		_costhetastar = drand48()*2.-1.;
		_Te = (sqrt(_pf*_pf+_me*_me)*_gammaCM+_gammaCM*_betaCM*_pf*_costhetastar)-_me ;
		_Tn = (sqrt(_pf*_pf+_mn*_mn)*_gammaCM-_gammaCM*_betaCM*_pf*_costhetastar)-_mn ;
		_thetae = acos((sqrt(_pf*_pf+_me*_me)*_gammaCM*_betaCM+_gammaCM*_pf*_costhetastar)/
					   sqrt((sqrt(_pf*_pf+_me*_me)*_gammaCM*_betaCM+_gammaCM*_pf*_costhetastar)*
							(sqrt(_pf*_pf+_me*_me)*_gammaCM*_betaCM+_gammaCM*_pf*_costhetastar)+
							(_pf*_pf*(1.-_costhetastar*_costhetastar))));
		
		_thetan = acos((sqrt(_pf*_pf+_mn*_mn)*_gammaCM*_betaCM-_gammaCM*_pf*_costhetastar)/
					   sqrt((sqrt(_pf*_pf+_mn*_mn)*_gammaCM*_betaCM-_gammaCM*_pf*_costhetastar)*
							(sqrt(_pf*_pf+_mn*_mn)*_gammaCM*_betaCM-_gammaCM*_pf*_costhetastar)+
							(_pf*_pf*(1.-_costhetastar*_costhetastar))));
	}
	else {
		_Tn = 0.;
		_Te = 0.;
		_thetae = 0.;
		_thetan = 0.;
	}

	
}

void protonibd_Vogel0::genereEvent( double enu ){
	_Tnu = enu;
	_Te = _Tnu-_delta;
	_ve = sqrt(_Te*_Te-_me*_me)/(_Te);
	
	computeInputKinematic();
	if (_betaCM != 0.) {
		
		_Aparameter = _sigma0*0.5*(_f*_f+3.*_g*_g)*sqrt(_Te*_Te-_me*_me)*(_Te);
		_Bparameter = _sigma0*0.5*(_f*_f-_g*_g)*_ve*sqrt(_Te*_Te-_me*_me)*(_Te);
		_costhetae = (sqrt(_Bparameter*_Bparameter+
						   (4.*drand48()-2.)*_Aparameter*_Bparameter+
						   _Aparameter*_Aparameter)
					  -_Aparameter)/_Bparameter;
		
		
		double atri = _pf*_pf*(_gammaCM*_gammaCM-_gammaCM*_gammaCM*_costhetae*_costhetae+_costhetae*_costhetae);
		double btri = 2.*_gammaCM*_gammaCM*_betaCM*_pf*sqrt((_pf*_pf+_me*_me))*(1.-_costhetae*_costhetae);
		double ctri = ( _gammaCM*_gammaCM*_betaCM*_betaCM
					   *(_pf*_pf+_me*_me)
					   *(1.-_costhetae*_costhetae)-_pf*_pf*_costhetae*_costhetae );
		double deltatri = btri*btri - 4.*atri*ctri;
		
		if (deltatri < 0.) {
			_Tn = 0.;
			_Te = 0.;
			_thetae = 0.;
			_thetan = 0.;
			return;
		}
		else if (deltatri == 0. ){
			_costhetastar = -btri/(2.*atri);
		}
		else {
			if (_costhetae < 0. ){
				_costhetastar = (-btri-sqrt(deltatri))/(2.*atri);
			}
			else {
				_costhetastar = (-btri+sqrt(deltatri))/(2.*atri);
			}
		}
		
		_Te = (sqrt(_pf*_pf+_me*_me)*_gammaCM+_gammaCM*_betaCM*_pf*_costhetastar)-_me ;
		_Tn = (sqrt(_pf*_pf+_mn*_mn)*_gammaCM-_gammaCM*_betaCM*_pf*_costhetastar)-_mn ;
		_thetae = acos((sqrt(_pf*_pf+_me*_me)*_gammaCM*_betaCM+_gammaCM*_pf*_costhetastar)/
					   sqrt((sqrt(_pf*_pf+_me*_me)*_gammaCM*_betaCM+_gammaCM*_pf*_costhetastar)*
							(sqrt(_pf*_pf+_me*_me)*_gammaCM*_betaCM+_gammaCM*_pf*_costhetastar)+
							(_pf*_pf*(1.-_costhetastar*_costhetastar))));
		_thetan = acos((sqrt(_pf*_pf+_mn*_mn)*_gammaCM*_betaCM-_gammaCM*_pf*_costhetastar)/
					   sqrt((sqrt(_pf*_pf+_mn*_mn)*_gammaCM*_betaCM-_gammaCM*_pf*_costhetastar)*
							(sqrt(_pf*_pf+_mn*_mn)*_gammaCM*_betaCM-_gammaCM*_pf*_costhetastar)+
							(_pf*_pf*(1.-_costhetastar*_costhetastar))));
	}
	else {
		_Tn = 0.;
		_Te = 0.;
		_thetae = 0.;
		_thetan = 0.;
	}
}


/////////////////////////////////////////////////////////////////////////////
/* Test function */

void isotropeKinematic(){
	protonibd * process = new protonibd;
	
	TGraph * geminpositron = new TGraph();
	TGraph * gemaxpositron = new TGraph();
	
	TMultiGraph * gpositron = new TMultiGraph();
	gpositron->SetTitle("; E_{nu} (MeV); E_{e}^{lab} (MeV)");
	gemaxpositron->SetLineColor( 2 );
	geminpositron->SetLineColor( 2 );
	gemaxpositron->SetLineWidth( 2 );
	geminpositron->SetLineWidth( 2 );
	gpositron->Add(gemaxpositron, "L");
	gpositron->Add(geminpositron, "L");

	TGraph * geminneutron = new TGraph();
	TGraph * gemaxneutron = new TGraph();
	
	TMultiGraph * gneutron = new TMultiGraph();
	gneutron->SetTitle("; E_{nu} (MeV); E_{e}^{lab} (MeV)");
	gemaxneutron->SetLineColor( 2 );
	geminneutron->SetLineColor( 2 );
	gemaxneutron->SetLineWidth( 2 );
	geminneutron->SetLineWidth( 2 );
	gneutron->Add(gemaxneutron, "L");
	gneutron->Add(geminneutron, "L");
	
	
	TH2D * h2EnuEe = new TH2D("h2EnuEe","; E_{#nu} (MeV); E_{e} (MeV)",500,0., 10., 500,0., 10. );
	TH2D * h2EnuEn = new TH2D("h2EnuEn","; E_{#nu} (MeV); E_{n} (MeV)",500,0., 10., 500,0., .5 );
	
	TH2D * h2Theta_n_Theta_e = new TH2D("h2Theta_n_Theta_e", "#frac{d#sigma}{d cos #theta^{*}=c} ; #theta_{e}; #theta_{n}",
										360., 0., 180., 
										180., 0., 90.);
	
    double Tth = process->getThreshold();
    for (double E = Tth; E < 10; E+= 0.0005) {
        gemaxpositron->SetPoint( gemaxpositron->GetN(), E, process->getEmaxPositron(E) );
		geminpositron->SetPoint( geminpositron->GetN(), E, process->getEminPositron(E) );
		
		
		gemaxneutron->SetPoint( gemaxneutron->GetN(), E, process->getEmaxNeutron(E) );
		geminneutron->SetPoint( geminneutron->GetN(), E, process->getEminNeutron(E) );
	}
    double enu;
	for (int i = 0; i < 100000000; i++) {
        enu = Tth + drand48()*(10. - Tth);
		process->genereEvent(enu);
		h2EnuEe->Fill( enu, process->getEPositron() );
		h2EnuEn->Fill( enu, process->getENeutron() );
		h2Theta_n_Theta_e->Fill( process->getThetaPositron()*180./M_PI,
								process->getThetaNeutron()*180./M_PI );
	}
	gStyle->SetPalette(1);
	TCanvas * can = new TCanvas();
	can->Divide( 2, 1, 1e-20, 1e-20 );
	can->GetPad(1)->cd();
	h2EnuEe->DrawCopy("colz");
	gpositron->DrawClone("L");
	
	can->GetPad(2)->cd();
	h2EnuEn->DrawCopy("colz");
	gneutron->DrawClone("L");
	
	TCanvas * can2 = new TCanvas();
	can2->GetPad(0)->cd();
	h2Theta_n_Theta_e->DrawCopy("colz");
	
}

void isotropeKinematic_Vogel(){
	protonibd * process = new protonibd_Vogel0;
	
	TGraph * geminpositron = new TGraph();
	TGraph * gemaxpositron = new TGraph();
	
	TMultiGraph * gpositron = new TMultiGraph();
	gpositron->SetTitle("; E_{nu} (MeV); E_{e}^{lab} (MeV)");
	gemaxpositron->SetLineColor( 2 );
	geminpositron->SetLineColor( 2 );
	gemaxpositron->SetLineWidth( 2 );
	geminpositron->SetLineWidth( 2 );
	gpositron->Add(gemaxpositron, "L");
	gpositron->Add(geminpositron, "L");
	
	TGraph * geminneutron = new TGraph();
	TGraph * gemaxneutron = new TGraph();
	
	TMultiGraph * gneutron = new TMultiGraph();
	gneutron->SetTitle("; E_{nu} (MeV); E_{e}^{lab} (MeV)");
	gemaxneutron->SetLineColor( 2 );
	geminneutron->SetLineColor( 2 );
	gemaxneutron->SetLineWidth( 2 );
	geminneutron->SetLineWidth( 2 );
	gneutron->Add(gemaxneutron, "L");
	gneutron->Add(geminneutron, "L");
	
	
	TH2D * h2EnuEe = new TH2D("h2EnuEe","; E_{#nu} (MeV); E_{e} (MeV)",500,0., 10., 500,0., 10. );
	TH2D * h2EnuEn = new TH2D("h2EnuEn","; E_{#nu} (MeV); E_{n} (MeV)",500,0., 10., 500,0., .5 );
	
	TH2D * h2Theta_n_Theta_e = new TH2D("h2Theta_n_Theta_e", "#frac{d#sigma}{d cos #theta_{e}} =a+b cos#theta_{e} ; #theta_{e}; #theta_{n}",
										360., 0., 180., 
										180., 0., 90.);
	
	TGraph * gCX = new TGraph();
	gCX->SetTitle("Vogel 0th; E_{#nu} (MeV); #sigma");
	
	
    double Tth = process->getThreshold();
    for (double E = Tth; E < 10; E+= 0.0005) {
        gCX->SetPoint( gCX->GetN(), E, process->getCX(E) );
		gemaxpositron->SetPoint( gemaxpositron->GetN(), E, process->getEmaxPositron(E) );
		geminpositron->SetPoint( geminpositron->GetN(), E, process->getEminPositron(E) );
		
		
		gemaxneutron->SetPoint( gemaxneutron->GetN(), E, process->getEmaxNeutron(E) );
		geminneutron->SetPoint( geminneutron->GetN(), E, process->getEminNeutron(E) );
	}
    double enu;
	for (int i = 0; i < 10000000; i++) {
        enu = Tth + drand48()*(10. - Tth);
        process->genereEvent(enu);
		h2EnuEe->Fill( enu, process->getEPositron() );
		h2EnuEn->Fill( enu, process->getENeutron() );
		h2Theta_n_Theta_e->Fill( process->getThetaPositron()*180./M_PI,
								process->getThetaNeutron()*180./M_PI );
	}
	gStyle->SetPalette(1);
	
	TCanvas * canCX = new TCanvas();
	canCX->cd();
	gCX->DrawClone("AL");
	
	TCanvas * can = new TCanvas();
	can->Divide( 2, 1, 1e-20, 1e-20 );
	can->GetPad(1)->cd();
	h2EnuEe->DrawCopy("colz");
	gpositron->DrawClone("L");
	
	can->GetPad(2)->cd();
	h2EnuEn->DrawCopy("colz");
	gneutron->DrawClone("L");
	
	TCanvas * can2 = new TCanvas();
	can2->GetPad(0)->cd();
	h2Theta_n_Theta_e->DrawCopy("colz");
	
}

void testEnergyCosTheta( double enu ){
	protonibd * processV = new protonibd_Vogel0;
	protonibd * process = new protonibd;

	TH1D * h1cosThetaE = new TH1D("h1cosThetaE", "; cos #theta_{e}; ", 200, -1., 1.);
	TH1D * h1cosThetaEV = new TH1D("h1cosThetaEv", "; cos #theta_{e}; ", 200, -1., 1.);
	
	TH1D * h1cosThetaN = new TH1D("h1cosThetaN", "; cos #theta_{n}; ", 200, -1., 1.);
	TH1D * h1cosThetaNV = new TH1D("h1cosThetaNv", "; cos #theta_{n}; ", 200, -1., 1.);
	h1cosThetaEV->SetLineColor( 2 );
	h1cosThetaNV->SetLineColor( 2 );
	
	
	for (int i = 0; i < 10000000; i++) {
		process->genereEvent(enu);
		processV->genereEvent(enu);
		h1cosThetaE->Fill( cos( process->getThetaPositron()) );
		h1cosThetaEV->Fill( cos( processV->getThetaPositron()) );
		h1cosThetaN->Fill( cos( process->getThetaNeutron()) );
		h1cosThetaNV->Fill( cos( processV->getThetaNeutron()) );
	}
	TCanvas * can = new TCanvas();
	can->Divide( 2, 1, 1e-20, 1e-20 );
	can->GetPad(1)->cd();
	h1cosThetaEV->DrawCopy();
	h1cosThetaE->DrawCopy("same");
	can->GetPad(2)->cd();
	h1cosThetaNV->DrawCopy();
	h1cosThetaN->DrawCopy("same");
}

void meanCosThetaEN( int stat ){
	protonibd * processV = new protonibd_Vogel0;
	protonibd * process = new protonibd;
	
	double thetaE;
	double thetaEv;
	double thetaN;
	double thetaNv;
	
	TGraph * gthetaE = new TGraph();
	TGraph * gthetaEv = new TGraph();
	TGraph * gthetaN = new TGraph();
	TGraph * gthetaNv = new TGraph();
	
	TMultiGraph * mgthetaE = new TMultiGraph();
	TMultiGraph * mgthetaN = new TMultiGraph();
	
	
	mgthetaE->SetTitle( ";E_{#nu} (MeV); <cos #theta_{e}>" );
	mgthetaN->SetTitle( ";E_{#nu} (MeV); <cos #theta_{n}>" );
	gthetaE->SetLineColor( 4 );
	gthetaEv->SetLineColor( 2 );
	gthetaN->SetLineColor( 4 );
	gthetaNv->SetLineColor( 2 );
	
	mgthetaE->Add( gthetaE, "L" );
	mgthetaE->Add( gthetaEv, "L" );
	
	mgthetaN->Add( gthetaN, "L" );
	mgthetaN->Add( gthetaNv, "L" );
	
	double dE = 0.0005;
    double Tth = process->getThreshold();
    for (double E = Tth; E < 10.; E+= dE) {
        if (E > 3.) {
			dE = 0.2;
		}
		thetaE = 0.;
		thetaEv = 0.;
		thetaN = 0.;
		thetaNv = 0.;
		for (int i = 0; i < stat; i++) {
			process->genereEvent(E);
			processV->genereEvent(E);
			if (process->getENeutron() != 0.) {
				thetaE += cos(process->getThetaPositron());
				thetaEv += cos(processV->getThetaPositron());
				thetaN += cos(process->getThetaNeutron());
				thetaNv += cos(processV->getThetaNeutron());
			}
			
		}
		if (process->getENeutron() != 0.){
			gthetaE->SetPoint( gthetaE->GetN(), E, thetaE/double(stat) );
			gthetaEv->SetPoint( gthetaEv->GetN(), E, thetaEv/double(stat) );
			gthetaN->SetPoint( gthetaN->GetN(), E, thetaN/double(stat) );
			gthetaNv->SetPoint( gthetaNv->GetN(), E, thetaNv/double(stat) );
		}
	}
	TCanvas * can = new TCanvas("can", "can");
	can->Divide( 2, 1, 1e-20, 1e-20 );
	can->GetPad(1)->cd();
	mgthetaE->DrawClone("ALP");
	can->GetPad(2)->cd();
	mgthetaN->DrawClone("ALP");
	
	
}