TEvtGen/EvtGenBase/EvtHepMCEvent.cpp

   1 //--------------------------------------------------------------------------
   2 //
   3 // Environment:
   4 //      This software is part of the EvtGen package. If you use all or part
   5 //      of it, please give an appropriate acknowledgement.
   6 //
   7 // Copyright Information: See EvtGen/COPYRIGHT
   8 //      Copyright (C) 2011      University of Warwick, UK
   9 //
  10 // Module: EvtHepMCEvent
  11 //
  12 // Description: Create an HepMC::GenEvent for the complete EvtParticle
  13 //              decay tree.
  14 //
  15 // Modification history:
  16 //
  17 //    John Back       June 2011            Module created
  18 //
  19 //------------------------------------------------------------------------
  20
  21 #include "EvtGenBase/EvtPatches.hh"
  22 #include "EvtGenBase/EvtHepMCEvent.hh"
  23 #include "EvtGenBase/EvtParticle.hh"
  24 #include "EvtGenBase/EvtPDL.hh"
  25
  26 #include "HepMC/Units.h"
  27
  28 EvtHepMCEvent::EvtHepMCEvent() :
  29   _theEvent(0),
  30   _translation(0.0, 0.0, 0.0, 0.0)
  31 {
  32 }
  33
  34 EvtHepMCEvent::~EvtHepMCEvent() {
  35   this->deleteEvent();
  36 }
  37
  38 void EvtHepMCEvent::deleteEvent() {
  39
  40   if (_theEvent != 0) {
  41     _theEvent->clear();
  42     delete _theEvent; _theEvent = 0;
  43   }
  44
  45 }
  46
  47 void EvtHepMCEvent::constructEvent(EvtParticle* baseParticle) {
  48
  49   EvtVector4R origin(0.0, 0.0, 0.0, 0.0);
  50   this->constructEvent(baseParticle, origin);
  51
  52 }
  53
  54 void EvtHepMCEvent::constructEvent(EvtParticle* baseParticle, EvtVector4R& translation) {
  55
  56   // This class does not take ownership of the base particle pointer.
  57   // Rather, it uses the base particle to construct the event.
  58
  59   this->deleteEvent();
  60   if (baseParticle == 0) {return;}
  61
  62   _theEvent = new HepMC::GenEvent(HepMC::Units::GEV, HepMC::Units::MM);
  63   _translation = translation;
  64
  65   // Use the recursive function addVertex to add a vertex with incoming/outgoing
  66   // particles. Adds a new vertex for any EvtParticles with decay daughters.
  67   // All particles are in the rest frame of the base particle ("lab frame").
  68
  69   HepMC::GenParticle* hepMCGenParticle = this->createGenParticle(baseParticle, EvtHepMCEvent::LAB);
  70
  71   this->addVertex(baseParticle, hepMCGenParticle);
  72
  73 }
  74
  75 HepMC::GenParticle* EvtHepMCEvent::createGenParticle(EvtParticle* theParticle, int frameType) {
  76
  77   // Create an HepMC GenParticle, with the 4-momenta in the frame given by the frameType integer
  78   HepMC::GenParticle* genParticle = 0;
  79
  80   if (theParticle != 0) {
  81
  82     // Set the particle status integer to either stable or decayed
  83     int status(EvtHepMCEvent::STABLE);
  84     int nDaug = theParticle->getNDaug();
  85     if (nDaug > 0) {status = EvtHepMCEvent::DECAYED;}
  86
  87     // Get the 4-momentum (E, px, py, pz) for the EvtParticle.
  88     EvtVector4R p4(0.0, 0.0, 0.0, 0.0);
  89
  90     if (frameType == EvtHepMCEvent::RESTFRAME) {
  91       p4 = theParticle->getP4Restframe();
  92     } else if (frameType == EvtHepMCEvent::LAB) {
  93       p4 = theParticle->getP4Lab();
  94     } else {
  95       p4 = theParticle->getP4();
  96     }
  97
  98     // Convert this to the HepMC 4-momentum
  99     double E = p4.get(0);
 100     double px = p4.get(1);
 101     double py = p4.get(2);
 102     double pz = p4.get(3);
 103
 104     HepMC::FourVector hepMC_p4(px, py, pz, E);
 105
 106     // Get the particle PDG integer id
 107     int PDGInt = EvtPDL::getStdHep(theParticle->getId());
 108
 109     genParticle = new HepMC::GenParticle(hepMC_p4, PDGInt, status);
 110
 111   }
 112
 113   return genParticle;
 114
 115 }
 116
 117 void EvtHepMCEvent::addVertex(EvtParticle* inEvtParticle, HepMC::GenParticle* inGenParticle) {
 118
 119   // This is a recursive function that adds GenVertices to the GenEvent for
 120   // the incoming EvtParticle and its daughters. We use two separate
 121   // pointers for the EvtParticle and GenParticle information: the former
 122   // to obtain the PDGId, 4-momenta, daughter and vertex positions, the latter to
 123   // set the incoming particle to the vertex. Note that the outgoing particle for
 124   // one vertex might be the incoming particle for another vertex - this needs to
 125   // be the same GenParticle pointer, hence the reason for using it as a 2nd argument
 126   // in this function.
 127
 128   if (_theEvent == 0 || inEvtParticle == 0 || inGenParticle == 0) {return;}
 129
 130   // Create the decay vertex
 131   HepMC::FourVector vtxCoord = this->getVertexCoord(inEvtParticle);
 132   HepMC::GenVertex* theVertex = new HepMC::GenVertex(vtxCoord);
 133
 134   // Add the vertex to the event
 135   _theEvent->add_vertex(theVertex);
 136
 137   // Set the incoming particle
 138   theVertex->add_particle_in(inGenParticle);
 139
 140   // Set the outgoing particles (decay products)
 141   int nDaug = inEvtParticle->getNDaug();
 142   int iDaug(0);
 143   // Loop over the daughters
 144   for (iDaug = 0; iDaug < nDaug; iDaug++) {
 145
 146     EvtParticle* evtDaughter = inEvtParticle->getDaug(iDaug);
 147     HepMC::GenParticle* genDaughter = this->createGenParticle(evtDaughter, EvtHepMCEvent::LAB);
 148
 149     if (genDaughter != 0) {
 150
 151       // Add a new GenParticle (outgoing) particle daughter to the vertex
 152       theVertex->add_particle_out(genDaughter);
 153
 154       // Find out if the daughter also has decay products.
 155       // If so, recursively run this function again.
 156       int nDaugProducts = evtDaughter->getNDaug();
 157
 158       if (nDaugProducts > 0) {
 159
 160         // Recursively process daughter particles and add their vertices to the event
 161         this->addVertex(evtDaughter, genDaughter);
 162
 163       } // Have daughter products
 164
 165     } // hepMCDaughter != 0
 166
 167   } // Loop over daughters
 168
 169 }
 170
 171 HepMC::FourVector EvtHepMCEvent::getVertexCoord(EvtParticle* theParticle) {
 172
 173   HepMC::FourVector vertexCoord(0.0, 0.0, 0.0, 0.0);
 174
 175   if (theParticle != 0 && theParticle->getNDaug() != 0) {
 176
 177     // Get the position (t,x,y,z) of the EvtParticle, offset by the translation vector.
 178     // This position will be the point where the particle decays. So we ask
 179     // the position of the (1st) daughter particle.
 180     EvtParticle* daugParticle = theParticle->getDaug(0);
 181
 182     if (daugParticle != 0) {
 183
 184       EvtVector4R vtxPosition = daugParticle->get4Pos() + _translation;
 185
 186       // Create the HepMC 4 vector of the position (x,y,z,t)
 187       vertexCoord.setX(vtxPosition.get(1));
 188       vertexCoord.setY(vtxPosition.get(2));
 189       vertexCoord.setZ(vtxPosition.get(3));
 190       vertexCoord.setT(vtxPosition.get(0));
 191
 192     }
 193
 194   }
 195
 196   return vertexCoord;
 197
 198 }