Line data Source code
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 0 : _theEvent(0),
30 0 : _translation(0.0, 0.0, 0.0, 0.0)
31 0 : {
32 0 : }
33 :
34 0 : EvtHepMCEvent::~EvtHepMCEvent() {
35 0 : this->deleteEvent();
36 0 : }
37 :
38 : void EvtHepMCEvent::deleteEvent() {
39 :
40 0 : if (_theEvent != 0) {
41 0 : _theEvent->clear();
42 0 : delete _theEvent; _theEvent = 0;
43 0 : }
44 :
45 0 : }
46 :
47 : void EvtHepMCEvent::constructEvent(EvtParticle* baseParticle) {
48 :
49 0 : EvtVector4R origin(0.0, 0.0, 0.0, 0.0);
50 0 : this->constructEvent(baseParticle, origin);
51 :
52 0 : }
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 0 : this->deleteEvent();
60 0 : if (baseParticle == 0) {return;}
61 :
62 0 : _theEvent = new HepMC::GenEvent(HepMC::Units::GEV, HepMC::Units::MM);
63 0 : _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 0 : HepMC::GenParticle* hepMCGenParticle = this->createGenParticle(baseParticle, EvtHepMCEvent::LAB);
70 :
71 0 : this->addVertex(baseParticle, hepMCGenParticle);
72 :
73 0 : }
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 0 : if (theParticle != 0) {
81 :
82 : // Set the particle status integer to either stable or decayed
83 : int status(EvtHepMCEvent::STABLE);
84 0 : int nDaug = theParticle->getNDaug();
85 0 : if (nDaug > 0) {status = EvtHepMCEvent::DECAYED;}
86 :
87 : // Get the 4-momentum (E, px, py, pz) for the EvtParticle.
88 0 : EvtVector4R p4(0.0, 0.0, 0.0, 0.0);
89 :
90 0 : if (frameType == EvtHepMCEvent::RESTFRAME) {
91 0 : p4 = theParticle->getP4Restframe();
92 0 : } else if (frameType == EvtHepMCEvent::LAB) {
93 0 : p4 = theParticle->getP4Lab();
94 0 : } else {
95 0 : p4 = theParticle->getP4();
96 : }
97 :
98 : // Convert this to the HepMC 4-momentum
99 0 : double E = p4.get(0);
100 0 : double px = p4.get(1);
101 0 : double py = p4.get(2);
102 0 : double pz = p4.get(3);
103 :
104 0 : HepMC::FourVector hepMC_p4(px, py, pz, E);
105 :
106 : // Get the particle PDG integer id
107 0 : int PDGInt = EvtPDL::getStdHep(theParticle->getId());
108 :
109 0 : genParticle = new HepMC::GenParticle(hepMC_p4, PDGInt, status);
110 :
111 0 : }
112 :
113 0 : return genParticle;
114 :
115 0 : }
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 0 : if (_theEvent == 0 || inEvtParticle == 0 || inGenParticle == 0) {return;}
129 :
130 : // Create the decay vertex
131 0 : HepMC::FourVector vtxCoord = this->getVertexCoord(inEvtParticle);
132 0 : HepMC::GenVertex* theVertex = new HepMC::GenVertex(vtxCoord);
133 :
134 : // Add the vertex to the event
135 0 : _theEvent->add_vertex(theVertex);
136 :
137 : // Set the incoming particle
138 0 : theVertex->add_particle_in(inGenParticle);
139 :
140 : // Set the outgoing particles (decay products)
141 0 : int nDaug = inEvtParticle->getNDaug();
142 : int iDaug(0);
143 : // Loop over the daughters
144 0 : for (iDaug = 0; iDaug < nDaug; iDaug++) {
145 :
146 0 : EvtParticle* evtDaughter = inEvtParticle->getDaug(iDaug);
147 0 : HepMC::GenParticle* genDaughter = this->createGenParticle(evtDaughter, EvtHepMCEvent::LAB);
148 :
149 0 : if (genDaughter != 0) {
150 :
151 : // Add a new GenParticle (outgoing) particle daughter to the vertex
152 0 : 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 0 : int nDaugProducts = evtDaughter->getNDaug();
157 :
158 0 : if (nDaugProducts > 0) {
159 :
160 : // Recursively process daughter particles and add their vertices to the event
161 0 : this->addVertex(evtDaughter, genDaughter);
162 :
163 0 : } // Have daughter products
164 :
165 0 : } // hepMCDaughter != 0
166 :
167 : } // Loop over daughters
168 :
169 0 : }
170 :
171 : HepMC::FourVector EvtHepMCEvent::getVertexCoord(EvtParticle* theParticle) {
172 :
173 0 : HepMC::FourVector vertexCoord(0.0, 0.0, 0.0, 0.0);
174 :
175 0 : 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 0 : EvtParticle* daugParticle = theParticle->getDaug(0);
181 :
182 0 : if (daugParticle != 0) {
183 :
184 0 : EvtVector4R vtxPosition = daugParticle->get4Pos() + _translation;
185 :
186 : // Create the HepMC 4 vector of the position (x,y,z,t)
187 0 : vertexCoord.setX(vtxPosition.get(1));
188 0 : vertexCoord.setY(vtxPosition.get(2));
189 0 : vertexCoord.setZ(vtxPosition.get(3));
190 0 : vertexCoord.setT(vtxPosition.get(0));
191 :
192 0 : }
193 :
194 0 : }
195 :
196 0 : return vertexCoord;
197 :
198 : }
|
