PPSUtilities.cc

Go to the documentation of this file.

#include "Utilities/PPS/interface/PPSUtilities.h"
#include "HepMC/GenParticle.h"
#include "H_BeamParticle.h"
#include "TLorentzVector.h"

TLorentzVector PPSTools::HectorParticle2LorentzVector(H_BeamParticle hp, int direction) {
  double partP = sqrt(pow(hp.getE(), 2) - ProtonMassSQ);
  double theta = sqrt(pow(hp.getTX(), 2) + pow(hp.getTY(), 2)) * urad;
  double pz = partP * cos(theta);
  double px = tan((double)hp.getTX() * urad) * pz;  //it is equivalente to PartP*sin(theta)*cos(phi);
  double py = tan((double)hp.getTY() * urad) * pz;  //it is equivalente to partP*sin(theta)*sin(phi);
  pz *= direction;
  return TLorentzVector(px, py, pz, hp.getE());
}

H_BeamParticle PPSTools::LorentzVector2HectorParticle(TLorentzVector p) {
  H_BeamParticle h_p;
  h_p.set4Momentum(p.Px(), p.Py(), abs(p.Pz()), p.E());
  return h_p;
}
void PPSTools::LorentzBoost(HepMC::GenParticle& p, const string& frame, FullBeamInfo const& bi) {
  TLorentzVector p_out;
  p_out.SetPx(p.momentum().px());
  p_out.SetPy(p.momentum().py());
  p_out.SetPz(p.momentum().pz());
  LorentzBoost(p_out, frame, bi);
  p.set_momentum(HepMC::FourVector(p_out.Px(), p_out.Py(), p_out.Pz(), p_out.E()));
}
void PPSTools::LorentzBoost(H_BeamParticle& h_p, int dir, const string& frame, FullBeamInfo const& bi) {
  TLorentzVector p_out = HectorParticle2LorentzVector(h_p, dir);
  LorentzBoost(p_out, frame, bi);
  h_p = LorentzVector2HectorParticle(p_out);
}
void PPSTools::LorentzBoost(TLorentzVector& p_out, const string& frame, FullBeamInfo const& bi) {
  const long double microrad = 1.e-6;
  //
  double px_P, py_P, pz_P;
  double px_N, py_N, pz_N;
  double fBoostAngle1 = 0.;
  double fBoostAngle2 = 0.;
  if (frame == "LAB") {
    fBoostAngle1 = bi.fCrossingAngleBeam1;
    fBoostAngle2 = bi.fCrossingAngleBeam2;
  }
  if (frame == "MC") {
    fBoostAngle1 = -bi.fCrossingAngleBeam1;
    fBoostAngle2 = -bi.fCrossingAngleBeam2;
  }
  px_P = bi.fBeamMomentum * sin(fBoostAngle2 * microrad);
  px_N = bi.fBeamMomentum * sin(fBoostAngle1 * microrad);
  pz_P = bi.fBeamMomentum * cos(fBoostAngle2 * microrad);
  pz_N = bi.fBeamMomentum * cos(fBoostAngle1 * microrad);
  py_P = 0.;
  py_N = 0.;

  TLorentzVector BeamP, BeamN, projVect;
  BeamP.SetPx(px_P);
  BeamP.SetPy(py_P);
  BeamP.SetPz(pz_P);
  BeamP.SetE(bi.fBeamEnergy);
  BeamN.SetPx(px_N);
  BeamN.SetPy(py_N);
  BeamN.SetPz(-pz_N);
  BeamN.SetE(bi.fBeamEnergy);
  projVect = BeamP + BeamN;
  TVector3 beta;
  TLorentzVector boosted = p_out;
  beta = projVect.BoostVector();
  boosted.Boost(beta);
  p_out = boosted;
}
void PPSTools::Get_t_and_xi(const TLorentzVector* proton, double& t, double& xi, LimitedBeamInfo const& bi) {
  t = 0.;
  xi = -1.;
  if (!proton)
    return;
  double mom =
      sqrt((proton->Px()) * (proton->Px()) + (proton->Py()) * (proton->Py()) + (proton->Pz()) * (proton->Pz()));
  if (mom > bi.fBeamMomentum)
    mom = bi.fBeamMomentum;
  double energy = proton->E();
  double theta = (proton->Pz() > 0) ? proton->Theta() : TMath::Pi() - proton->Theta();
  t = -2. * (ProtonMassSQ - bi.fBeamEnergy * energy + bi.fBeamMomentum * mom * cos(theta));
  xi = (1.0 - energy / bi.fBeamEnergy);
}