PPS/src/PPSUtilities.cc

0001 #include "Utilities/PPS/interface/PPSUtilities.h"
0002 #include "HepMC/GenParticle.h"
0003 #include "H_BeamParticle.h"
0004 #include "TLorentzVector.h"
0005
0006 TLorentzVector PPSTools::HectorParticle2LorentzVector(H_BeamParticle hp, int direction) {
0007   double partP = sqrt(pow(hp.getE(), 2) - ProtonMassSQ);
0008   double theta = sqrt(pow(hp.getTX(), 2) + pow(hp.getTY(), 2)) * urad;
0009   double pz = partP * cos(theta);
0010   double px = tan((double)hp.getTX() * urad) * pz;  //it is equivalente to PartP*sin(theta)*cos(phi);
0011   double py = tan((double)hp.getTY() * urad) * pz;  //it is equivalente to partP*sin(theta)*sin(phi);
0012   pz *= direction;
0013   return TLorentzVector(px, py, pz, hp.getE());
0014 }
0015
0016 H_BeamParticle PPSTools::LorentzVector2HectorParticle(TLorentzVector p) {
0017   H_BeamParticle h_p;
0018   h_p.set4Momentum(p.Px(), p.Py(), abs(p.Pz()), p.E());
0019   return h_p;
0020 }
0021 void PPSTools::LorentzBoost(HepMC::GenParticle& p, const string& frame, FullBeamInfo const& bi) {
0022   TLorentzVector p_out;
0023   p_out.SetPx(p.momentum().px());
0024   p_out.SetPy(p.momentum().py());
0025   p_out.SetPz(p.momentum().pz());
0026   LorentzBoost(p_out, frame, bi);
0027   p.set_momentum(HepMC::FourVector(p_out.Px(), p_out.Py(), p_out.Pz(), p_out.E()));
0028 }
0029 void PPSTools::LorentzBoost(H_BeamParticle& h_p, int dir, const string& frame, FullBeamInfo const& bi) {
0030   TLorentzVector p_out = HectorParticle2LorentzVector(h_p, dir);
0031   LorentzBoost(p_out, frame, bi);
0032   h_p = LorentzVector2HectorParticle(p_out);
0033 }
0034 void PPSTools::LorentzBoost(TLorentzVector& p_out, const string& frame, FullBeamInfo const& bi) {
0035   const long double microrad = 1.e-6;
0036   //
0037   double px_P, py_P, pz_P;
0038   double px_N, py_N, pz_N;
0039   double fBoostAngle1 = 0.;
0040   double fBoostAngle2 = 0.;
0041   if (frame == "LAB") {
0042     fBoostAngle1 = bi.fCrossingAngleBeam1;
0043     fBoostAngle2 = bi.fCrossingAngleBeam2;
0044   }
0045   if (frame == "MC") {
0046     fBoostAngle1 = -bi.fCrossingAngleBeam1;
0047     fBoostAngle2 = -bi.fCrossingAngleBeam2;
0048   }
0049   px_P = bi.fBeamMomentum * sin(fBoostAngle2 * microrad);
0050   px_N = bi.fBeamMomentum * sin(fBoostAngle1 * microrad);
0051   pz_P = bi.fBeamMomentum * cos(fBoostAngle2 * microrad);
0052   pz_N = bi.fBeamMomentum * cos(fBoostAngle1 * microrad);
0053   py_P = 0.;
0054   py_N = 0.;
0055
0056   TLorentzVector BeamP, BeamN, projVect;
0057   BeamP.SetPx(px_P);
0058   BeamP.SetPy(py_P);
0059   BeamP.SetPz(pz_P);
0060   BeamP.SetE(bi.fBeamEnergy);
0061   BeamN.SetPx(px_N);
0062   BeamN.SetPy(py_N);
0063   BeamN.SetPz(-pz_N);
0064   BeamN.SetE(bi.fBeamEnergy);
0065   projVect = BeamP + BeamN;
0066   TVector3 beta;
0067   TLorentzVector boosted = p_out;
0068   beta = projVect.BoostVector();
0069   boosted.Boost(beta);
0070   p_out = boosted;
0071 }
0072 void PPSTools::Get_t_and_xi(const TLorentzVector* proton, double& t, double& xi, LimitedBeamInfo const& bi) {
0073   t = 0.;
0074   xi = -1.;
0075   if (!proton)
0076     return;
0077   double mom =
0078       sqrt((proton->Px()) * (proton->Px()) + (proton->Py()) * (proton->Py()) + (proton->Pz()) * (proton->Pz()));
0079   if (mom > bi.fBeamMomentum)
0080     mom = bi.fBeamMomentum;
0081   double energy = proton->E();
0082   double theta = (proton->Pz() > 0) ? proton->Theta() : TMath::Pi() - proton->Theta();
0083   t = -2. * (ProtonMassSQ - bi.fBeamEnergy * energy + bi.fBeamMomentum * mom * cos(theta));
0084   xi = (1.0 - energy / bi.fBeamEnergy);
0085 }