eicsmear/erhic/ParticleMC.h

0001 /**
0002  \file
0003  Declaration of class erhic::ParticleMC.
0004
0005  \author    Thomas Burton
0006  \date      2011-10-10
0007  \copyright 2011 Brookhaven National Lab
0008  */
0009
0010 #ifndef INCLUDE_EICSMEAR_ERHIC_PARTICLEMC_H_
0011 #define INCLUDE_EICSMEAR_ERHIC_PARTICLEMC_H_
0012
0013 #include <string>
0014
0015 #include <TLorentzVector.h>
0016 #include <TRef.h>
0017 #include <TVector3.h>
0018
0019 #include "eicsmear/erhic/Pid.h"
0020 #include "eicsmear/erhic/VirtualParticle.h"
0021
0022 namespace erhic {
0023
0024   class ParticleMCeA: public TObject {
0025   public:
0026     explicit ParticleMCeA() {};
0027     virtual ~ParticleMCeA() {};
0028
0029     // Let them all be public;
0030     //added by liang to include add on particle data structure
0031     Int_t charge;
0032     Int_t massNum;
0033     Int_t NoBam; ///< 0, 2 create in hard collisions not affected by INC: 0 inside nucleus, 2 outside nucleus;
0034                              ///< 3 create in evaporation; 4 heaviest remnant in evaporation;
0035                              ///< 11-34 created in INC, 10+n, n=IDCH cascade generated in INC
0036     ClassDef(ParticleMCeA, 1)
0037 };
0038
0039 class EventMC;
0040 /**
0041  A particle produced by a Monte Carlo generator.
0042  */
0043 class ParticleMCbase : public VirtualParticle {
0044  public:
0045   explicit ParticleMCbase();
0046   /**
0047    Destructor
0048    */
0049   virtual ~ParticleMCbase() {};
0050
0051   /**
0052    Print the contents of Particle to standard output.
0053    The format is that of the input Monte Carlo i.e.
0054    I KS id orig daughter ldaughter px py pz E m xv yv zv.
0055    Inherited from TObject. The argument is unused.
0056    */
0057   virtual void Print(Option_t* = "") const;
0058
0059   /**
0060    Returns the particle index in an event, in the range [1, N].
0061    */
0062   virtual UInt_t GetIndex() const;
0063
0064   /**
0065    Returns the status of the particle.
0066    The meaning of the status code depends on the generator.
0067    For PYTHIA, see the description of variable K(I,1) in the manual.
0068    */
0069   virtual UShort_t GetStatus() const;
0070
0071   /**
0072    Returns the index of this particle's parent in an event.
0073    */
0074   virtual UShort_t GetParentIndex() const;
0075   virtual UShort_t GetParentIndex1() const;
0076
0077   /**
0078    Returns the index of this particle's first child particle.
0079    Returns 0 if this particle has no children.
0080    */
0081   virtual UShort_t GetChild1Index() const;
0082
0083   /**
0084    Returns the index of this particle's last child particle.
0085    Returns 0 if this particle has zero or one children.
0086    */
0087   virtual UShort_t GetChildNIndex() const;
0088
0089   /**
0090    Returns the number of children of this particle.
0091    Returns 0 if the particle did not decay.
0092    */
0093   virtual UInt_t GetNChildren() const;
0094
0095   /**
0096    Returns the x component of 3-momentum.
0097    */
0098   virtual Double_t GetPx() const;
0099
0100   /**
0101    Returns the y component of 3-momentum.
0102    */
0103   virtual Double_t GetPy() const;
0104
0105   /**
0106    Returns the z component of 3-momentum.
0107    */
0108   virtual Double_t GetPz() const;
0109
0110   /**
0111    Returns invariant mass (GeV/c<sup>2</sup>).
0112    */
0113   virtual Double_t GetM() const;
0114
0115   /**
0116    Returns momentum transverse to the beam direction.
0117    */
0118   virtual Double_t GetPt() const;
0119
0120   /**
0121    Returns the origin point of the particle (cm).
0122    (0,0,0) indicates a particle originating in the collision.
0123    */
0124   virtual TVector3 GetVertex() const;
0125
0126   /**
0127    Returns the identity information of this particle's parent.
0128
0129    */
0130   virtual Pid GetParentId() const;
0131
0132   /**
0133    Returns the total momentum (GeV).
0134    */
0135   virtual Double_t GetP() const;
0136
0137   /**
0138    Returns the polar angle in the range [0, pi] radians.
0139    */
0140   virtual Double_t GetTheta() const;
0141
0142   /**
0143    Returns the polar angle in the range [0, 2pi] radians.
0144    */
0145   virtual Double_t GetPhi() const;
0146
0147   /**
0148    Returns the rapidity.
0149    */
0150   virtual Double_t GetRapidity() const;
0151
0152   /**
0153    Returns the pseudorapidity.
0154    */
0155   virtual Double_t GetEta() const;
0156
0157   /**
0158    Returns the variable z.
0159    z = (P.p_h)/(P.q).
0160    */
0161   virtual Double_t GetZ() const;
0162
0163   /**
0164    Returns Feynman-x.
0165    x<sub>F</sub> = 2*p<sub>z</sub>/sqrt(s).
0166    */
0167   virtual Double_t GetXFeynman() const;
0168
0169   /**
0170    Returns the angle with respect to the exchange boson.
0171    Defined in the beam hadron's rest frame.
0172    Given in the range [0,pi] radians.
0173    */
0174   virtual Double_t GetThetaVsGamma() const;
0175
0176   /**
0177    Returns the p<sub>T</sub> with respect to the exchange boson.
0178    Defined in the beam hadron's rest frame.
0179    */
0180   virtual Double_t GetPtVsGamma() const;
0181
0182   /**
0183    Returns a pointer to the event containing this particle.
0184    Returns NULL if this particle has not been associated
0185    with an event.
0186    */
0187   const EventMC* GetEvent() const;
0188
0189   /**
0190    Set the event with which to associate this particle.
0191    */
0192   void SetEvent(EventMC* event);
0193
0194   /**
0195    Returns the (E,p) 4-vector in the lab frame.
0196    */
0197   virtual TLorentzVector Get4Vector() const;
0198
0199   /**
0200    Returns the (E,p) 4-vector in the lab frame.
0201    */
0202   virtual TLorentzVector PxPyPzE() const { return Get4Vector(); }
0203
0204   /**
0205    Returns the (E,p) 4-vector in the hadron-boson frame.
0206    This frame is defined such that
0207    <UL>
0208    <LI> the beam hadron is at rest </LI>
0209    <LI> the z direction is the exchange boson momentum vector </LI>
0210    <LI> the y direction is defined as q x e`, where q is the boson and
0211    e` is the scattered lepton momentum </LI>
0212    <LI> x is defined to complete the right-handed coordinate system </LI>
0213    </UL>
0214
0215    \note Due to details of the implementation and how ROOT handles reading
0216          events in a TTree, this function will not work for the exchange
0217          boson when using a TTree::Draw (or similar) statement. It does
0218          work for the exchange boson if reading events manually via
0219          TTree::GetEntry.
0220    */
0221   virtual TLorentzVector Get4VectorInHadronBosonFrame() const;
0222
0223   /**
0224    Returns the energy of the particle in the lab frame.
0225    */
0226   virtual Double_t GetE() const;
0227
0228   virtual void SetE(Double_t);
0229
0230   virtual void SetM(Double_t);
0231
0232   virtual void SetP(Double_t);
0233
0234   virtual void SetPt(Double_t);
0235
0236   virtual void SetPz(Double_t);
0237
0238   virtual void SetPhi(Double_t);
0239
0240   virtual void SetTheta(Double_t);
0241
0242   virtual void SetStatus(UShort_t);
0243
0244   /**
0245    Returns the ID of the particle.
0246    */
0247   virtual Pid Id() const;
0248
0249   /**
0250    Returns the ID of the particle.
0251    */
0252   virtual Pid GetPdgCode() const { return Id(); }
0253
0254   /**
0255    Sets quantities derived from the four-momentum (E, px, py, pz), namely
0256    <UL>
0257    <LI> total momentum </LI>
0258    <LI> transverse momentum </LI>
0259    <LI> rapidity </LI>
0260    <LI> pseudorapidity </LI>
0261    <LI> theta </LI>
0262    <LI> phi </LI>
0263    </UL>
0264    This should be called if (E, px, py, pz) are manually altered
0265    in order to propagate the changes to these other quantities.
0266    */
0267   virtual void ComputeDerivedQuantities();
0268
0269   /**
0270    Sets quantities that depend on the properties of the event or
0271    associations of one particle with another, namely
0272    <UL>
0273    <LI> z </LI>
0274    <LI> Feynman x </LI>
0275    <LI> angle and pt with respect to the exchanged boson </LI>
0276    <LI> azimuthal angle around the exchanged boson </LI>
0277    <LI> parent pdg code </LI>
0278    </UL>
0279    Important: this particle is assumed to be in the same frame of
0280    reference as those contained in the event that is passed as an
0281    argument.
0282    */
0283   virtual void ComputeEventDependentQuantities(EventMC&);
0284
0285   /**
0286    Sets the index of the particle i.e. its position in the track list
0287    (in principle this can be any
0288    integer you require to associated with the particle).
0289    */
0290   virtual void SetIndex(int i) { I = i; }
0291
0292   /** Sets the status code of the particle (generally final state
0293    particles are given status == 1 */
0294   virtual void SetStatus(int i) { KS = i; }
0295
0296   /** Sets the ID of the particle. In order to make use of class Pid this
0297    should be the PDG code of the particle, but in principle can be any
0298    value you wish to use to identify it. */
0299   virtual void SetId(int i) { id = i; }
0300
0301   /** Sets the index of this particle's parent if it has one.
0302    By default this is zero, indicating no parent. */
0303   virtual void SetParentIndex(int i) { orig = i; }
0304
0305   /** Sets the index of this particle's other parent if it has one.
0306    By default this is zero, indicating no parent. */
0307   virtual void SetParentIndex1(int i) { orig1 = i; }
0308
0309   /** Sets the index of this particle's first child. By default this
0310    is zero, indicating no children. */
0311   virtual void SetChild1Index(int i) { daughter = i; }
0312
0313   /** Sets the index of this particle's last child.
0314    By default this is zero, indication zero or one children. */
0315   virtual void SetChildNIndex(int i) { ldaughter = i; }
0316
0317   /**
0318    Sets the four-momentum of the particle.
0319    Changes are propagated to derived quantities.
0320    */
0321   virtual void Set4Vector(const TLorentzVector&);
0322
0323   /**
0324    Sets the origin coordinates
0325    */
0326   virtual void SetVertex(const TVector3&);
0327
0328   /**
0329    Sets the ID of this particle's parent. See comments in SetId()
0330    */
0331   virtual void SetParentId(int i) { parentId = i; }
0332
0333 // protected:
0334
0335   UShort_t    I;             ///< Particle index in event
0336   //UShort_t    KS;            ///< Particle status code: see PYTHIA manual
0337   //modified by liang to include negative KS values
0338   Int_t       KS;            ///< Particle status code: see PYTHIA manual
0339   Int_t       id;            ///< PDG particle code
0340   // Looks strange (see also access methods); keep like this for sort of
0341   // backward compatibility;
0342   UShort_t    orig;          ///< I of parent particle
0343   UShort_t    orig1;         ///< I of parent particle1
0344   UShort_t    daughter;      ///< I of first child particle
0345   UShort_t    ldaughter;     ///< I of last child particle
0346
0347   // AYK, 2017/04/02: changed {px,py,pz} & {xv,yv,zv} variable types
0348   // from Double32_t (float) to Double_t (double);
0349   Double_t px;           ///< x component of particle momentum
0350   Double_t py;           ///< y component of particle momentum
0351   Double_t pz;           ///< z component of particle momentum
0352
0353   // Reducing file size.
0354   // The following variables could be removed, at the cost of increased
0355   // running time to calculate them on-the-fly:
0356   // p - compute from px/y/z
0357   // m - accessible via PID with TDatabasePDG
0358   // E - compute from p and m
0359   // parentId - accessible via event record and parent index IF a
0360   //          persistent pointer to the containing event can be stored
0361   // daughter - IF daughter particles always start immediately after this one
0362   // pt - from px, py
0363   // theta - from pt, pz
0364   // phi - from px, py
0365   // rapidity - from E, pz
0366   // eta - from theta
0367   // xFeynman - if event is accessible
0368   // thetaGamma - duplicated in pHadronBoson
0369   // ptVsGamma - duplicated in pHadronBoson
0370   // phiPrf - duplicated in pHadronBoson
0371
0372   Double32_t E;            ///< Energy of particle
0373   Double32_t m;            ///< Invariant mass of particle
0374   Double32_t pt;           ///< Transverse momentum of particle
0375   Double_t xv;           ///< x coordinate of particle production vertex
0376   Double_t yv;           ///< y coordinate of particle production vertex
0377   Double_t zv;           ///< z coordinate of particle production vertex
0378
0379   // Can be deprecated if parent particle itself is available
0380   Int_t parentId;          ///< PDG code of this particle's parent
0381
0382   Double32_t p;            ///< Total momentum of particle
0383   Double32_t theta;        ///< Polar angle
0384   Double32_t phi;          ///< Azimuthal angle
0385   Double32_t rapidity;     ///< Rapidity of particle
0386   Double32_t eta;          ///< Pseudorapidity of particle
0387   Double32_t z;            ///< Fraction of virtual photon energy
0388                            ///< carried by particle
0389   Double32_t xFeynman;     ///< Feynman x = p<sub>z</sub>/(2sqrt(s))
0390   Double32_t thetaGamma;   ///< Angle between particle and the exchange
0391                            ///< boson in the hadron beam rest frame
0392   Double32_t ptVsGamma;    ///< pt w.r.t. the virtual photon in the
0393                            ///< hadron beam rest frame
0394   Double32_t thetaGammaHCM;   ///< Angle between particle and the exchange
0395                               ///< boson in the hadron-boson rest frame
0396   Double32_t ptVsGammaHCM;    ///< pt w.r.t. the virtual photon in the
0397                               ///< hadron boson rest frame
0398   Double32_t phiPrf;       ///< Azimuthal angle around virtual
0399                            ///< photon in hadron beam rest frame
0400
0401   TRef event;  ///< Persistent reference to the event containing
0402                ///< this particle.
0403
0404   ClassDef(ParticleMCbase, 5)
0405 };
0406
0407 // The unfortunate consequence of adding 'eA' pointer to the ParticleMC
0408 // class is that default copy ctor does not work; the easiest workaround
0409 // (unless somebody would like to copy all the other variables by hand)
0410 // is to offload these variables to an intermediate base class;
0411 class ParticleMC : public ParticleMCbase {
0412  public:
0413   /**
0414    Default constructor.
0415    Optionally pass a string with particle information in a HEPEVT
0416    format, namely:
0417      "index status id parent firstChild lastChild px py pz E m xv yv zv"
0418    */
0419   explicit ParticleMC(): eA(0) {};
0420   explicit ParticleMC(const std::string&, bool eAflag);
0421   // So this is the evil copy ctor, which is damn simple now;
0422  ParticleMC(const ParticleMC &src): ParticleMCbase(src) {
0423     eA = src.eA ? new ParticleMCeA(*src.eA) : 0; //orig1 = src.orig1;
0424   };
0425
0426   ~ParticleMC();
0427
0428   /**
0429    Returns a pointer to the parent of this particle.
0430    This is the particle with index GetParentIndex() in the
0431    event containing this particle (obtainable via GetEvent()).
0432    Returns NULL if this particle has no parent or if it cannot
0433    be accessed via GetEvent().
0434    */
0435   virtual const ParticleMC* GetParent() const;
0436
0437   /**
0438    Returns a pointer to the nth child particle of this particle,
0439    where n is in the range [0, GetNChildren()).
0440    GetChild(0) returns the child whose index in this particle's
0441    event is GetChild1Index().
0442    GetChild(GetNChildren()-1) returns the child whose index
0443    is GetChildNIndex().
0444    Returns NULL if there is no such child particle or it
0445    cannot be accessed via the event for some reason (see GetEvent()).
0446    */
0447   virtual const ParticleMC* GetChild(UShort_t) const;
0448
0449   /**
0450    Returns true if n in the range [0, N), where N is the number
0451    of children of this particle.
0452    Returns false otherwise.
0453    Equivalent to GetChild(n) != NULL.
0454    */
0455   virtual Bool_t HasChild(Int_t) const;
0456
0457   // FIXME: let it be public?; NB: TRef (like for ParticleMCbase.event) is
0458   // not needed here, since it is just a POD instance with a single reference;
0459   ParticleMCeA *eA;
0460   //UShort_t    orig1;          ///< I of parent particle1
0461
0462   ClassDef(ParticleMC, 2)
0463 };
0464
0465 inline UInt_t ParticleMCbase::GetIndex() const {
0466   return I;
0467 }
0468
0469 inline UShort_t ParticleMCbase::GetStatus() const {
0470   return KS;
0471 }
0472
0473 inline UShort_t ParticleMCbase::GetParentIndex() const {
0474   return orig;
0475 }
0476 inline UShort_t ParticleMCbase::GetParentIndex1() const {
0477   return orig1;
0478 }
0479
0480 inline UShort_t ParticleMCbase::GetChild1Index() const {
0481   return daughter;
0482 }
0483
0484 inline UShort_t ParticleMCbase::GetChildNIndex() const {
0485   return ldaughter;
0486 }
0487
0488 inline Double_t ParticleMCbase::GetPx() const {
0489   return px;
0490 }
0491
0492 inline Double_t ParticleMCbase::GetPy() const {
0493   return py;
0494 }
0495
0496 inline Double_t ParticleMCbase::GetPz() const {
0497   return pz;
0498 }
0499
0500 inline Double_t ParticleMCbase::GetM() const {
0501   return m;
0502 }
0503
0504 inline Double_t ParticleMCbase::GetPt() const {
0505   return pt;
0506 }
0507
0508 inline TVector3 ParticleMCbase::GetVertex() const {
0509   return TVector3(xv, yv, zv);
0510 }
0511
0512 inline erhic::Pid ParticleMCbase::GetParentId() const {
0513   return erhic::Pid(parentId);
0514 }
0515
0516 inline Double_t ParticleMCbase::GetP() const {
0517   return p;
0518 }
0519
0520 inline Double_t ParticleMCbase::GetTheta() const {
0521   return theta;
0522 }
0523
0524 inline Double_t ParticleMCbase::GetPhi() const {
0525   return phi;
0526 }
0527
0528 inline Double_t ParticleMCbase::GetRapidity() const {
0529   return rapidity;
0530 }
0531
0532 inline Double_t ParticleMCbase::GetEta() const {
0533   return eta;
0534 }
0535
0536 inline Double_t ParticleMCbase::GetZ() const {
0537   return z;
0538 }
0539
0540 inline Double_t ParticleMCbase::GetXFeynman() const {
0541   return xFeynman;
0542 }
0543
0544 inline Double_t ParticleMCbase::GetThetaVsGamma() const {
0545   return thetaGamma;
0546 }
0547
0548 inline Double_t ParticleMCbase::GetPtVsGamma() const {
0549   return ptVsGamma;
0550 }
0551
0552 inline Pid ParticleMCbase::Id() const {
0553   return Pid(id);
0554 }
0555
0556 inline UInt_t ParticleMCbase::GetNChildren() const {
0557   if (0 == daughter) return 0;
0558   if (0 == ldaughter) return 1;
0559   return ldaughter - daughter + 1;
0560 }
0561
0562 inline Double_t ParticleMCbase::GetE() const {
0563   return E;
0564 }
0565
0566 inline void ParticleMCbase::SetE(Double_t e) {
0567   E = e;
0568 }
0569
0570 inline void ParticleMCbase::SetM(Double_t mass) {
0571   m = mass;
0572 }
0573
0574 inline void ParticleMCbase::SetP(Double_t momentum) {
0575   p = momentum;
0576 }
0577
0578 inline void ParticleMCbase::SetPt(Double_t momentum) {
0579   pt = momentum;
0580 }
0581
0582 inline void ParticleMCbase::SetPz(Double_t momentum) {
0583   pz = momentum;
0584 }
0585
0586 inline void ParticleMCbase::SetPhi(Double_t value) {
0587   phi = value;
0588 }
0589
0590 inline void ParticleMCbase::SetTheta(Double_t value) {
0591   theta = value;
0592 }
0593
0594 inline void ParticleMCbase::SetStatus(UShort_t status) {
0595   KS = status;
0596 }
0597
0598 }  // namespace erhic
0599
0600 #endif  // INCLUDE_EICSMEAR_ERHIC_PARTICLEMC_H_