// -*- C++ -*-
#include "FinalStateBreitCurrent.hh"
#include "Rivet/Cmp.hh"
#include "Rivet/Projections/DISKinematics.hh"
#include "Rivet/Rivet.hh"
#include "Rivet/Math/Constants.hh"
#include "Rivet/Tools/ParticleIdUtils.hh"
#include "math.h"

int hemispherea=0;
double elima=0;
 
namespace Rivet {


  int FinalStateBreitCurrent::compare(const Projection& p) const {
    return mkNamedPCmp(p, "Kinematics");
  }

  FinalStateBreitCurrent::FinalStateBreitCurrent(const DISKinematics& kinematicsp){//definition for just boost and not cuts on hemispheres
    setName("FinalStateBreitCurrent");
    addProjection(kinematicsp, "Kinematics");
    addProjection(Beam(), "Beam");
    addProjection(DISLepton(), "Lepton");
  }

  FinalStateBreitCurrent::FinalStateBreitCurrent(const DISKinematics& kinematicsp,const int& hemisphere, const double& elim){//definition for version with cuts on hemispheres
    setName("FinalStateBreitCurrent");
    addProjection(kinematicsp, "Kinematics");
    addProjection(Beam(), "Beam");
    addProjection(DISLepton(), "Lepton");
    hemispherea=hemisphere;
    elima=elim;
  }

  void FinalStateBreitCurrent::project(const Event& e) {
    const DISKinematics& diskin = applyProjection<DISKinematics>(e, "Kinematics");
    const LorentzTransform breitboost = diskin.boostBreit();
    const DISLepton& dislep = diskin.applyProjection<DISLepton>(e, "Lepton");
    const GenParticle& dislepGP = dislep.out().genParticle();
    const FinalState& fs = dislep.applyProjection<FinalState>(e, "FS");
    Particle _inHadron;
    LorentzTransform _hcm;
    
    //This is mostly to find the breit shifted PGamma to make sure that the 
    //final state particles I store are in the chosen hemisphere:

    // Identify beam hadron
    const ParticlePair& inc = applyProjection<Beam>(e, "Beam").beams();
    bool firstIsHadron  = PID::isHadron(inc.first.pdgId());
    bool secondIsHadron = PID::isHadron(inc.second.pdgId());
    if (firstIsHadron && !secondIsHadron) {
      _inHadron = inc.first;
    } else if (!firstIsHadron && secondIsHadron) {
      _inHadron = inc.second;
    } else {
       //help!
      throw Error("DISKinematics projector could not find the correct beam hadron");
    }
    
    // several key DIS kinematics needed to work out the boost:
    FourMomentum pLepIn = dislep.in().momentum();
    FourMomentum pLepOut = dislep.out().momentum();
    FourMomentum pGamma = pLepIn - pLepOut;
    FourMomentum pHad = _inHadron.momentum();
    FourMomentum tothad = pGamma + pHad;
    double Q2 = -pGamma.mass2();
    double dx = diskin.x();

    // Calculate boost vector for boost into HCM-system:
    LorentzTransform tmp;
    tmp.setBoost(-tothad.boostVector());
    
    // Rotate so the photon is in x-z plane in HCM rest frame:
    FourMomentum pGammaHCM = tmp.transform(pGamma);
    tmp.preMult(Matrix3(Vector3::mkZ(), -pGammaHCM.azimuthalAngle()));
    pGammaHCM = tmp.transform(pGamma);
    
     // Rotate so the photon is along the positive z-axis:
    const double rot_angle = pGammaHCM.polarAngle() * (pGammaHCM.px() >= 0 ? -1 : 1);
    tmp.preMult(Matrix3(Vector3::mkY(), rot_angle));
    pGammaHCM = tmp.transform(pGamma);
 
    // Finally rotate so outgoing lepton at phi = 0:
    FourMomentum pLepOutHCM = tmp.transform(pLepOut);
    tmp.preMult(Matrix3(Vector3::mkZ(), -pLepOutHCM.azimuthalAngle()));
    _hcm = tmp;
    
    // Boost to Breit frame (use opposite convention for photon --- along *minus* z):
    tmp.preMult(Matrix3(Vector3::mkX(), PI));
    const double bz = 1 - 2*dx;
    LorentzTransform _breit = LorentzTransform(Vector3::mkZ() * bz).combine(tmp);
    
    //calculate the photons 3 momentum, aswell as the pre&post boost energys for cuts in the next stage:
    FourMomentum pGammaB = _breit.transform(pGamma);
    double bq= pGammaB.mod();
    double q= sqrt(Q2);
    Vector3 pgam = pGammaB.vector3();
    
    
    //And now I loop over and actually shift the particles:
    
    // Fill the particle list with all particles _other_ than the DIS scattered
    // lepton, with momenta boosted into the Breit frame.
    _theParticles.clear();
    _theParticles.reserve(fs.particles().size());
    for (ParticleVector::const_iterator p = fs.particles().begin(); p != fs.particles().end(); ++p) {
      const GenParticle& loopGP = p->genParticle();
      if (&loopGP != &dislepGP) { //< Ensure that we skip the DIS lepton
        Particle temp = *p;
        const FourMomentum breitMom = breitboost.transform(temp.momentum());
	Vector3 bp =breitMom.vector3();
	if((hemispherea==0)||(hemispherea==1&&(dot(pgam,bp)>0))||(hemispherea==2&&(dot(pgam,bp)<0))){//Make sure in the current hemisphere
	  if((elima==0.)||bq>((elima)*q)){//Make sure over the energy limit >.25Q in the breit frame
	    temp.setMomentum(breitMom);
	    _theParticles.push_back(temp);
	  }
	}
      }
    }
  }
  
  
  
  
  
}