#include "SolidGenerateIBD.hh"
#include "globals.hh"
#include "Randomize.hh"

#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TLorentzVector.h"

#include "G4ThreeVector.hh"
#include "G4DynamicParticle.hh"
#include "G4SystemOfUnits.hh"
#include "G4PhysicalConstants.hh"

SolidGenerateIBD::SolidGenerateIBD(const char* RootFile){
	
	
	_ratefile = new TFile(RootFile);
	_anuerate = (TH1D*)_ratefile->Get("anuerate");

    _ibdKinematics = new protonibd_Vogel0();



	
}

SolidGenerateIBD::~SolidGenerateIBD(){
	
}

double SolidGenerateIBD::randomWithinInterval(double lowerLimit, double upperLimit){
	double randomValue;
	randomValue=lowerLimit+G4UniformRand()*(upperLimit-lowerLimit); 
	return randomValue;
}

double SolidGenerateIBD::randomNeutrinoEnergy(){
	bool valid = false;
	int counter = 0.;
    double enu;
    double threshold = _ibdKinematics->getThreshold();

	while (valid == false){
        enu = randomWithinInterval(threshold, 15.00);
		double prob = _anuerate->GetBinContent(_anuerate->FindBin(enu));
		//cout<<" E_a1, prob,  "<<E_a1<<" "<<prob<<std::endl;
		if ( prob  >  randomWithinInterval(0.,1.)) valid = true;
		++counter;
	}
	return enu;
}

void SolidGenerateIBD::NewInteraction(G4DynamicParticle* positron, G4DynamicParticle* neutron,
                                      G4ThreeVector antineutrinoMomentumDirection, G4double energy)
{

    if( energy < 0. ) // default IBD
        energy = randomNeutrinoEnergy();

    _ibdKinematics->genereEvent(energy);
	
    G4double ePositron      = _ibdKinematics->getEPositron();//MeV
    G4double eNeutron       = _ibdKinematics->getENeutron(); //MeV
    G4double thetaPositron  = _ibdKinematics->getThetaPositron();   //theta in pi
    G4double thetaNeutron   = _ibdKinematics->getThetaNeutron();//theta in pi
    G4double phiPositron    = randomWithinInterval(0, 2.0*pi);   //phi in pi
    G4double phiNeutron     = phiPositron - pi;
	
	positron->SetKineticEnergy(ePositron);
	neutron->SetKineticEnergy(eNeutron);

	G4ThreeVector* positronMomentumDirection= new G4ThreeVector(antineutrinoMomentumDirection);
	//G4cout << "Positron original momentum direction x=" << positronMomentumDirection->getX()
	//	   <<" y= "	<< positronMomentumDirection->getY()
	//	   <<" z= "	<< positronMomentumDirection->getZ()
	//	   <<G4endl;

	positronMomentumDirection->rotate(antineutrinoMomentumDirection.orthogonal(),thetaPositron);
	positronMomentumDirection->rotate(antineutrinoMomentumDirection,phiPositron);
	//G4cout << "Positron rotated1 momentum direction x=" << positronMomentumDirection->getX()
	//	   <<" y= "	<< positronMomentumDirection->getY()
	//	   <<" z= "	<< positronMomentumDirection->getZ()
	//	   <<G4endl;

	G4ThreeVector* neutronMomentumDirection= new G4ThreeVector(antineutrinoMomentumDirection);
	//G4cout << "Neutron original momentum direction x=" << neutronMomentumDirection->getX()
	//	   <<" y= "	<< neutronMomentumDirection->getY()
	//	   <<" z= "	<< neutronMomentumDirection->getZ()
	//	   <<G4endl;
	neutronMomentumDirection->rotate(antineutrinoMomentumDirection.orthogonal(),thetaNeutron);
	neutronMomentumDirection->rotate(antineutrinoMomentumDirection,phiNeutron);
	//G4cout << "Neutron rotated1 momentum direction x=" << neutronMomentumDirection->getX()
	//	   <<" y= "	<< neutronMomentumDirection->getY()
	//	   <<" z= "	<< neutronMomentumDirection->getZ()
	//	   <<G4endl;

	//G4cout << "Antineutrino momentum direction x=" << antineutrinoMomentumDirection.getX()
	//	   <<" y= "	<< antineutrinoMomentumDirection.getY()
	//	   <<" z= "	<< antineutrinoMomentumDirection.getZ()
	//	   <<G4endl;

	positron->SetMomentumDirection(*positronMomentumDirection);
	neutron->SetMomentumDirection(*neutronMomentumDirection);


	// Print out positron and neutron data
	/**
	G4cout << "e+ " << "Kinetic energy: "
			<<qb1_lab.E()-m_e
			<<" Direction: "
			<<qb1_lab.Px()/qb1_lab.P()<<" "
			<<qb1_lab.Py()/qb1_lab.P()<<" "
			<<qb1_lab.Pz()/qb1_lab.P()<<G4endl;

	G4cout << "n " << "Kinetic energy: "
			<<E_b2-m_n
			<<" Direction: "
			<<qb2_lab.Px()/Pmag_b2<<" "
			<<qb2_lab.Py()/Pmag_b2<<" "
			<<qb2_lab.Pz()/Pmag_b2<<G4endl;
	**/

}