#include "GenGamma.h"

// code Bethe-Heitler distribution

GenGamma::GenGamma() {

  // reset the event
  fEnergy=20.; // to prevent zero's in the analysis
  fXPos=0.;
  fYPos=0.;
  fZ =0.5;
 
  // set up parameters

  // beam parameters
  fEproton=920.;
  fEelectron=27.5;

  // beam profile parameters
  fXTilt=0.;
  fXWidth=1.5;
  fYTilt=0.;
  fYWidth=0.7;

  fEnergyMin=5.;

  // the radiation parameters
  fRadXmin =0.005;
  fRadXmean = 0.0585;

  // set up Bethe-Heitler
  fBethe = new TF1("fBethe","4*[2]*([0]-x)/([0])/x *([0]/([0]-x)+([0]-x)/[0]-2/3)*(log(4*[1]*[0]*([0]-x)/[3]/x)-1/2)",fEnergyMin,fEelectron);
  fBethe->SetParameter(0,fEelectron); //electron energy
  fBethe->SetParameter(1,fEproton);  //proton energy
  fBethe->SetParameter(2,1/137.*((2.8*1e-15)*(2.8*1e-15)/(1e-28))); //alpha QED*re^2/1barn
  fBethe->SetParameter(3,1*0.0005); //Mp*Me

  
  SetEdist(fBethe); // make it to be the default energy distribution

  // code yhit distribution
  fydistgaus = new TF1("fydistgaus","1/sqrt(2*3.1415)/[1]*exp((-(x-[0])**2)/2/[1]**2)",-5.,5.);
  fydistgaus->SetParameter(0,fYTilt);
  fydistgaus->SetParameter(1,fYWidth);

  SetYDist(fydistgaus); // make it a default distribution

  // code xdist distribution
  fxdistgaus = new TF1("fxdistgaus","1/sqrt(2*3.1415)/[1]*exp((-(x-[0])**2)/2/[1]**2)",-5.,5.);
  fxdistgaus->SetParameter(0,fXTilt);
  fxdistgaus->SetParameter(1,fXWidth);

  SetXDist(fxdistgaus); // make it a default distribution

  // code z=Ee/Egamma distribution
  SetZDist(new TF1("zdist","(1/0.995169*9/7*(x^2+(1-x)^2)+1/1.035*2/3*x*(1-x))",0.,1.));


}

GenGamma::~GenGamma() {

  delete fBethe; 
  delete fydistgaus; 
  delete fxdistgaus;
  delete fzdist;

}

float GenGamma::EnergyFracLoss() {
  //fractional energy loss of electrons in the window
  // the distribution function is given by GEANT studies
  // 31% events radiate more than 0.5% of energy 
  // with distribution function f(x)=1/x
  // EnergyLoss = x*EnergyGenerated
  //
  // to check the sensitivity on the parameters the
  // two numbers frac=#of events radiating more than
  // some fraction x-min are related through the 
  // conversion rate and the <x> - the mean energy loss in the window
  // 
  // the formula is frac = <x>/(1-xmin) *(-ln(xmin));

  if (fRadXmin==0.) {
    fRadXmin=0.005;
    cout << "parameter RadXmin=0; changed to 0.5%"<<endl;
  }
  float frac =-fRadXmean/(1-fRadXmin)*log(fRadXmin);
  if (frac>1.||frac<0.) cout << "Warning: GenGamma:EnergyFracLoss() fraction= "<<frac<<endl;
  float ran1,ran2; //random numbers
  if ((ran1=gRandom->Rndm())>frac) {
    return 0.; // no radiation
  }
  else {
    ran2=gRandom->Rndm();
    float x=exp(ran2 * log(fRadXmin));
    return x;
  }
  return 0.; // should never get here
}


void GenGamma::SetEnergyMin(float Emin) {
  fEnergyMin=Emin;
  fBethe->SetRange(Emin,fBethe->GetXmax());
  fBethe->Update();
}

void GenGamma::SetXTilt(float xtilt) {
  fXTilt=xtilt;
}

void GenGamma::SetYTilt(float ytilt) {
  fYTilt=ytilt;
}

void GenGamma::SetXWidth(float xwidth) {
  fXWidth=xwidth;
  fxdistgaus->SetParameter(1,xwidth);
  fxdistgaus->Update();
}

void GenGamma::SetYWidth(float ywidth) {
  fYWidth=ywidth;
  fydistgaus->SetParameter(1,ywidth);
  fydistgaus->Update();
}

void GenGamma::SetEProton(float Eproton) {
  fEproton=Eproton;
  fBethe->SetParameter(1,fEproton);
  fBethe->Update();
}

void GenGamma::SetEElectron(float Eelectron) {
  fEelectron=Eelectron;
  fBethe->SetParameter(0,fEelectron);
  fBethe->Update();
}

float GenGamma::GetXsec() {
  // return the total cross-section corresponding Energy distribution function
  float xsec;
  xsec = fEdist->Integral(fEdist->GetXmin(),fEdist->GetXmax());
  return xsec;
}

void GenGamma::MakeEvent() {
  fEnergy = fEdist->GetRandom();
  fXPos   = fXTilt+fxdist->GetRandom();
  fYPos   = fYTilt+fydist->GetRandom();
  fZ      = fzdist->GetRandom();
  fx1     = EnergyFracLoss();
  fx2     = EnergyFracLoss();
}

void GenGamma::Show() {
  cout << "Energy="<<fEnergy<<" x="<<fXPos<<" y="<<fYPos<<" z="<<fZ<<endl;
}