src/erhic/ParticleMC.cxx

0001 /**
0002  \file
0003  Implementation of class erhic::ParticleMC.
0004
0005  \author    Thomas Burton
0006  \date      2011-10-10
0007  \copyright 2011 Brookhaven National Lab
0008  */
0009
0010 #include "eicsmear/erhic/ParticleMC.h"
0011
0012 #include <iostream>
0013 #include <limits>
0014 #include <sstream>
0015 #include <stdexcept>
0016 #include <string>
0017
0018 #include <TLorentzRotation.h>
0019 #include <TParticlePDG.h>
0020 #include <TRotation.h>
0021
0022 #include "eicsmear/erhic/EventBase.h"
0023 #include "eicsmear/functions.h"
0024
0025 namespace {
0026
0027 /*
0028  Returns the boost to transform to the rest frame of the vector "rest".
0029  If z is non-NULL, rotate the frame so that z AFTER BOOSTING
0030  defines the positive z direction of that frame.
0031  e.g. To boost a gamma-p interaction from the lab frame to the
0032  proton rest frame with the virtual photon defining z:
0033  computeBoost(protonLab, photonLab);
0034  */
0035 TLorentzRotation computeBoost(const TLorentzVector& rest,
0036                               const TLorentzVector* z) {
0037   TLorentzRotation toRest(-(rest.BoostVector()));
0038   if (z) {
0039     TRotation rotate;
0040     TLorentzVector boostedZ(*z);
0041     boostedZ *= toRest;
0042     rotate.SetZAxis(boostedZ.Vect());
0043     // We need the rotation of the frame, so take the inverse.
0044     // See the TRotation documentation for more explanation.
0045     rotate = rotate.Inverse();
0046     toRest.Transform(rotate);
0047   }  // if
0048   return toRest;
0049 }
0050
0051 }  // anonymous namespace
0052
0053 namespace erhic {
0054
0055   ParticleMCbase::ParticleMCbase()
0056 : I(0)
0057 , KS(0)
0058 , id(0)
0059 , orig(0)
0060 , orig1(0)
0061 , daughter(0)
0062 , ldaughter(0)
0063 , px(0.)
0064 , py(0.)
0065 , pz(0.)
0066 , E(0.)
0067 , m(0.)
0068 , pt(0.)
0069 , xv(0.)
0070 , yv(0.)
0071 , zv(0.)
0072 , parentId(std::numeric_limits<Int_t>::min())
0073 , p(0.)
0074 , theta(0.)
0075 , phi(0.)
0076 , rapidity(0.)
0077 , eta(0.)
0078 , z(0.)
0079 , xFeynman(0.)
0080 , thetaGamma(0.)
0081 , ptVsGamma(0.)
0082 , thetaGammaHCM(0.)
0083 , ptVsGammaHCM(0.)
0084 , phiPrf(0.)
0085   {
0086   }
0087
0088   ParticleMC::ParticleMC(const std::string &line, bool eAflag): ParticleMCbase(), eA(0)
0089 {
0090   // Initialise to nonsense values to make input errors easy to spot
0091   if (!line.empty()) {
0092     static std::stringstream ss;
0093     ss.str("");
0094     ss.clear();
0095     ss << line;
0096     if (eAflag) {
0097       eA = new ParticleMCeA();
0098
0099       ss >>
0100     //changed by liang to add particle mother1
0101     I >> KS >> id >> orig1 >> orig >> daughter >> ldaughter >>
0102     px >> py >> pz >> E >> m >> xv >> yv >> zv
0103     //added by liang to include add particle data structure
0104      >> eA->massNum >> eA->charge >> eA->NoBam;
0105
0106       //eA->pz = pz; orig1 = eA->orig1;
0107     } else {
0108       ss >>
0109     I >> KS >> id >> orig >> daughter >> ldaughter >>
0110     px >> py >> pz >> E >> m >> xv >> yv >> zv;
0111       orig1 = 0;
0112     } //if
0113       // We should have no stream errors and should have exhausted
0114       // the whole of the stream filling the particle.
0115     if (ss.fail() || !ss.eof()) {
0116       throw std::runtime_error("Bad particle input: " + line);
0117     }  // if
0118     ComputeDerivedQuantities();
0119   }  // if
0120 }
0121
0122 ParticleMC::~ParticleMC()
0123 {
0124   if (eA) delete eA;
0125 }
0126
0127   // FIXME: may also want to print out ParticleMCeA variables?;
0128 void ParticleMCbase::Print(Option_t* /* option */) const {
0129   std::cout << I << '\t' << KS << '\t' << id << '\t' << orig << '\t' <<
0130   daughter << '\t' << ldaughter << '\t' << px << '\t' << py << '\t' << pz
0131   << '\t' << E << '\t' << m << '\t' << xv << '\t' << yv << '\t' << zv <<
0132   std::endl;
0133 }
0134
0135 void ParticleMCbase::ComputeDerivedQuantities() {
0136   // Calculate quantities that depend only on the properties already read.
0137   pt = sqrt(pow(px, 2.) + pow(py, 2.));
0138   p = sqrt(pow(pt, 2.) + pow(pz, 2.));
0139   // Rapidity and pseudorapidity
0140   Double_t Epluspz = E + pz;
0141   Double_t Eminuspz = E - pz;
0142   Double_t Ppluspz = p + pz;
0143   Double_t Pminuspz = p - pz;
0144   if (Eminuspz <= 0.0 || Pminuspz == 0.0 ||
0145      Ppluspz == 0.0 || Epluspz <= 0.0) {
0146     // Dummy values to avoid zero or infinite arguments in calculations
0147     rapidity = -19.;
0148     eta = -19.;
0149   } else {
0150     rapidity = 0.5 * log(Epluspz / Eminuspz);
0151     eta = 0.5 * log(Ppluspz / Pminuspz);
0152   }  // if
0153   theta = atan2(pt, pz);
0154   phi = TVector2::Phi_0_2pi(atan2(py, px));
0155 }
0156
0157 void ParticleMCbase::ComputeEventDependentQuantities(EventMC& event) {
0158   try {
0159     // Get the beam hadon, beam lepton and exchange boson.
0160     const TLorentzVector& hadron = event.BeamHadron()->Get4Vector();
0161     const TLorentzVector& lepton = event.ScatteredLepton()->Get4Vector();
0162     // Determine the exchange boson 4-vector from the scattered lepton,
0163     // since we're not always guaranteed to have one (weak neutral current e.g.)
0164     // const TLorentzVector& boson = event.ExchangeBoson()->Get4Vector();
0165     const TLorentzVector boson = event.BeamLepton()->Get4Vector() - lepton;
0166
0167     // Calculate z using the 4-vector definition,
0168     // so we don't care about frame of reference.
0169     z = hadron.Dot(Get4Vector()) / hadron.Dot(boson);
0170     // Calculate properties in the proton rest frame.
0171     // We want pT and angle with respect to the virtual photon,
0172     // so use that to define the z axis.
0173     TLorentzRotation toHadronRest = computeBoost(hadron, &boson);
0174     // Boost this particle to the proton rest frame and calculate its
0175     // pT and angle with respect to the virtual photon:
0176     TLorentzVector p = (Get4Vector() *= toHadronRest);
0177     thetaGamma = p.Theta();
0178     ptVsGamma =  p.Pt();
0179     // Calculate phi angle around virtual photon according
0180     // to the HERMES convention.
0181     TLorentzVector bosonPrf = (TLorentzVector(boson) *= toHadronRest);
0182     TLorentzVector leptonPrf = (TLorentzVector(lepton) *= toHadronRest);
0183     phiPrf = computeHermesPhiH(p, leptonPrf, bosonPrf);
0184     // Feynman x with xF = 2 * pz / W in the boson-hadron CM frame.
0185     // First boost to boson-hadron centre-of-mass frame.
0186     // Use the photon to define the z direction.
0187     TLorentzRotation boost = computeBoost(boson + hadron, &boson);
0188     xFeynman = 2. * (Get4Vector() *= boost).Pz() / sqrt(event.GetW2());
0189
0190     thetaGammaHCM = (Get4Vector() *= boost).Theta();
0191     ptVsGammaHCM =  (Get4Vector() *= boost).Pt();
0192
0193     // Determine the PDG code of the parent particle, if the particle
0194     // has a parent and the parent is present in the particle array.
0195     // The index of the particles from the Monte Carlo runs from [1,N]
0196     // while the index in the array runs from [0,N-1], so subtract 1
0197     // from the parent index to find its position.
0198     if (event.GetNTracks() > unsigned(orig - 1)) {
0199       parentId = event.GetTrack(orig - 1)->Id();
0200     }  // if
0201   }  // try
0202   catch(std::exception& e) {
0203     std::cerr <<
0204     "Exception in Particle::ComputeEventDependentQuantities: " <<
0205     e.what() << std::endl;
0206   }  // catch
0207 }
0208
0209 TLorentzVector ParticleMCbase::Get4Vector() const {
0210   return TLorentzVector(px, py, pz, E);
0211 }
0212
0213 const EventMC* ParticleMCbase::GetEvent() const {
0214   return static_cast<EventMC*>(event.GetObject());
0215 }
0216
0217 const ParticleMC* ParticleMC::GetChild(UShort_t u) const {
0218   // Look up this particle's child via the event
0219   // containing it and the index of the child in that event.
0220   if (!GetEvent()) {
0221     return NULL;
0222   }  // if
0223      // index is in the range [1,N]
0224   unsigned idx = daughter + u;
0225   if (daughter < 1 ||  // If first daughter index = 0, it has no children
0226       u >= GetNChildren()) {  // Insufficient children
0227     return NULL;
0228   }  // if
0229   --idx;  // Convert [1,N] --> [0,N)
0230   const ParticleMC* p(NULL);
0231   // Check this index is within the # of particles in the event
0232   if (idx < GetEvent()->GetNTracks()) {
0233     p = GetEvent()->GetTrack(idx);
0234   }  // if
0235   return p;
0236 }
0237
0238 const ParticleMC* ParticleMC::GetParent() const {
0239   // Look up this particle's parent via the event
0240   // containing it and the index of the parent in that event.
0241   const ParticleMC* p(NULL);
0242   if (GetEvent()) {
0243     if (GetEvent()->GetNTracks() >= GetParentIndex()) {
0244       p = GetEvent()->GetTrack(GetParentIndex() - 1);
0245     }  // if
0246   }  // if
0247   return p;
0248 }
0249
0250 Bool_t ParticleMC::HasChild(Int_t pdg) const {
0251   for (UInt_t i(0); i < GetNChildren(); ++i) {
0252     if (!GetChild(i)) {
0253       continue;
0254     }  // if
0255     if (pdg == GetChild(i)->Id()) {
0256       return true;
0257     }  // if
0258   }  // for
0259   return false;
0260 }
0261
0262 TLorentzVector ParticleMCbase::Get4VectorInHadronBosonFrame() const {
0263   double p_(0.), e_(0.), px_(ptVsGammaHCM), py_(0.), pz_(0.);
0264   // Calculate mangitude of momentum from pT and polar angle in
0265   // hadron-boson frame. If theta is ~parallel to the beam just set
0266   // p to whatever pT is (not that it makes all that much sense).
0267   if (thetaGammaHCM > 1.e-6) {
0268     p_ = ptVsGammaHCM / sin(thetaGammaHCM);
0269   }  // if
0270   // Deal with virtual particles later, so check if particle is off mass-shell
0271   if (!(m < 0.)) {
0272     e_ = sqrt(pow(p_, 2.) + pow(m, 2.));
0273   } else {
0274     e_ = sqrt(pow(p_, 2.) - pow(m, 2.));
0275     if (TMath::IsNaN(e_)) {
0276       e_ = 0.;
0277     }  // if
0278   }  // if
0279   // Calculate pZ from pT and theta, unless it's very close to the beam,
0280   // in which case set it to p.
0281   if (thetaGammaHCM > 1.e-6) {
0282     pz_ = ptVsGammaHCM / tan(thetaGammaHCM);
0283   }  // if
0284   // Calculate px and py, unless a dummy phi value is present.
0285   if (phiPrf > -100.) {
0286     px_ = ptVsGammaHCM * cos(phiPrf);
0287     py_ = ptVsGammaHCM * sin(phiPrf);
0288   }  // if
0289   // If we ended up with no energy, it's likely ths is the exchange boson,
0290   // as nothing will have happened above. If so, try to reference the event
0291   // record to get the necessary information, as we don't have enough
0292   // in the particle itself.
0293   // Note that this appears not to work in a TTree as ROOT doesn't read
0294   // the event when it reads the particle, though I'm not certain of the
0295   // exact cause.
0296   if (m < 0. && GetEvent()) {
0297     if (GetEvent()->ExchangeBoson()) {
0298       // Don't check if GetEvent()->ExchangeBoson() == this, in case this
0299       // particle is a copy of the event track that isn't part of a copied
0300       // event.
0301       if (GetEvent()->ExchangeBoson()->GetIndex() == GetIndex()) {
0302         // If it's the exchange boson just set E == nu.
0303         e_ = GetEvent()->GetNu();
0304         // Calculate p, careful about negative mass.
0305         p_ = sqrt(pow(e_, 2.) + pow(m, 2.));
0306         pz_ = p_ * cos(thetaGammaHCM);
0307       }  // if
0308     }  // if
0309   }  // if
0310   return TLorentzVector(px_, py_, pz_, e_);
0311 }
0312
0313 void ParticleMCbase::SetEvent(EventMC* e) {
0314   event = e;
0315 }
0316
0317 void ParticleMCbase::Set4Vector(const TLorentzVector& v) {
0318   E = v.Energy();
0319   px = v.Px();
0320   py = v.Py();
0321   pz = v.Pz();
0322   m = v.M();
0323   ComputeDerivedQuantities();  // Rapidity etc
0324   // If an event reference is set, recalculate event-dependent quantities
0325   // like z, xF
0326   EventMC* ev = static_cast<EventMC*>(event.GetObject());
0327   if (ev) {
0328     ComputeEventDependentQuantities(*ev);
0329   }  // if
0330 }
0331
0332 void ParticleMCbase::SetVertex(const TVector3& v) {
0333   xv = v.X();
0334   yv = v.Y();
0335   zv = v.Z();
0336 }
0337 }  // namespace erhic