#include "lumifast.h"
// constructor
LumiFast::LumiFast() {

  //  fBdl       = 0.275; // changed on oct 15th after edge strips were removed
  fBdl       = 0.296;
  fysep      = 8.4;
  fyoffset   = 0.;
  fxoffset   = 0.;
  fnxstrips  = 14;
  fnystrips_up = 9;
  fnystrips_dn = 9;
  fzmagnet   = 300.;
  fzspectro  = 1200;
  fEnergy_res= 0.18;
  fPos_res   = 0.07;
  fEscale_up = 1.;
  fEscale_dn = 1.;
  fEwindowMin= 17.;
  fEwindowMax= 23.;

  fIsScattering=1;
  fIsEnergyLoss=1;


  // set distribution functions
  fThetaDist = new TF1("fThetaDist","x*exp(-4200.*x)",0.,0.002);
 

  // set gamma generator
  fGenInternal = new GenGamma();
  fGenerator = fGenInternal;

}

LumiFast::~LumiFast() {
  delete fGenInternal;
  delete fThetaDist;
}




void LumiFast::NewEvent() {
  fGenerator->MakeEvent();
  ProcessEvent();
}

void LumiFast::ProcessEvent() {

  // this processes the event stored in the generator
  // after consecutive calls the event is not regenerated so 
  // the photons with the same energy and z are reconstructed
  //


  float Ehit = fGenerator->GetEnergy();
  float z = fGenerator->GetZ();
  
  float theta1 =0.,theta2=0.;
  float phi1=0.,phi2=0.;

  if (fIsScattering) {
    theta1=fThetaDist->GetRandom();
    theta2=fThetaDist->GetRandom();
    phi1 = 2*3.1415*gRandom->Rndm();
    phi2 = 2*3.1415*gRandom->Rndm();
  }


  fEnergy_e = z*Ehit;
  fEnergy_p = (1-z)*Ehit;

  // do radiation
  if (fIsEnergyLoss) {
    fEnergy_e*=(1-fGenerator->GetX1());
    fEnergy_p*=(1-fGenerator->GetX2());
  }
  
  if (fEnergy_e==0.) fEnergy_e=0.001; // this is to avoid division by zero
  if (fEnergy_p==0.) fEnergy_p=0.001;

  


  fxpos_gam = fGenerator->GetXPos();
  fypos_gam = fGenerator->GetYPos();

  // smear with resolution
  fEnergy_up = fEscale_up * gRandom->Gaus(fEnergy_e,Eres(fEnergy_e));
  fEnergy_dn = fEscale_dn * gRandom->Gaus(fEnergy_p,Eres(fEnergy_p)); 
 
  // get positions
  fypos_up = fypos_gam + fzspectro*theta1*cos(phi1)+GetBeta()/fEnergy_e;
  fypos_dn = fypos_gam - fzspectro*theta2*cos(phi2)-GetBeta()/fEnergy_p; 
  fxpos_up = fxpos_gam + fzspectro*theta1*sin(phi1);
  fxpos_dn = fxpos_gam + fzspectro*theta2*sin(phi2);

  // smear positions
  fypos_up = gRandom->Gaus(fypos_up,fPos_res);
  fypos_dn = gRandom->Gaus(fypos_dn,fPos_res);

  fxpos_up = gRandom->Gaus(fxpos_up,fPos_res);
  fxpos_dn = gRandom->Gaus(fxpos_dn,fPos_res);

  return;
}


int LumiFast::IsCoincidence() {
  
  return IsHitUp()&&IsHitDn()&&GetEnergyGam()<GetEWindowMax()&&GetEnergyGam()>GetEWindowMin();
}

int LumiFast::IsHitUp() {
  if (GetXPosUp()<GetXMax()&&GetXPosUp()>GetXMin()&&GetYPosUp()<GetYMaxUp()&&GetYPosUp()>GetYMinUp()) {
    return 1;
  }
  else 
    return 0;
}

int LumiFast::IsHitDn() {
  if (GetXPosDn()<GetXMax()&&GetXPosDn()>GetXMin()&&GetYPosDn()<GetYMaxDn()&&GetYPosDn()>GetYMinDn())
  {
    return 1;
  }
  else { 
    return 0;
  }
}

void LumiFast::Show() {
  // display event - debugging tool
  cout << "      True       Reconstructed" <<endl;

  cout << "Eup  "<<GetEnergyTrueUp()<< "  " << GetEnergyUp()<<endl;
  cout << "Edn  "<<GetEnergyTrueDn()<< "  " << GetEnergyDn()<<endl;
  cout << "xgam "<<GetXPosGamTrue()<<  "  " << GetXPosGam()<<endl;
  cout << "ygam "<<GetYPosGamTrue()<<  "  " << GetYPosGam()<<endl;
  cout << "z    "<<GetZTrue()      <<  "  " << GetZ()      <<endl<<endl;
  cout << "xdn="<<GetXPosDn() << " ydn="<< GetYPosDn() <<endl;
  cout << "xup="<<GetXPosUp() << " yup="<< GetYPosUp() <<endl;
  cout << "iscoincidence="<<IsCoincidence()<<" ishitup=" << IsHitUp()<<" ishitdn="<<IsHitDn()<<endl;
  return;
}