PhotosHepMC3Particle.cxx
1#include "HepMC3/GenEvent.h"
2#include "HepMC3/GenVertex.h"
3#include "HepMC3/GenParticle.h"
4#include "HepMC3/Print.h"
5#include "PhotosHepMC3Particle.h"
6#include "Log.h"
7#include "Photos.h"
8
9using namespace HepMC3;
10
11namespace Photospp
12{
13
15 m_particle = make_shared<GenParticle>();
16}
17
18PhotosHepMC3Particle::PhotosHepMC3Particle(int pdg_id, int status, double mass){
19 m_particle = make_shared<GenParticle>();
20 m_particle->set_pid(pdg_id);
21 m_particle->set_status(status);
22 m_particle->set_generated_mass(mass);
23}
24
26 m_particle = particle;
27}
28
32 // clear(m_created_particles);
33}
34
35
36//delete the TauolaHepMC3Particle objects
37void PhotosHepMC3Particle::clear(std::vector<PhotosParticle*> v){
38 while(v.size()!=0){
39 PhotosParticle * temp = v.back();
40 v.pop_back();
41 delete temp;
42 }
43}
44
46 return m_particle;
47}
48
49void PhotosHepMC3Particle::setMothers(vector<PhotosParticle*> mothers){
50
51 /******** Deal with mothers ***********/
52
54
55 //If there are mothers
56 if(mothers.size()>0){
57
58 GenParticlePtr part;
59 part=dynamic_cast<PhotosHepMC3Particle*>(mothers.at(0))->getHepMC3();
60
61 //Use end vertex of first mother as production vertex for particle
62 GenVertexPtr production_vertex = part->end_vertex();
63 GenVertexPtr orig_production_vertex = production_vertex;
64
65 if(!production_vertex){ //if it does not exist create it
66 production_vertex = make_shared<GenVertex>();
67 part->parent_event()->add_vertex(production_vertex);
68 }
69
70 //Loop over all mothers to check that the end points to the right place
71 vector<PhotosParticle*>::iterator mother_itr;
72 for(mother_itr = mothers.begin(); mother_itr != mothers.end();
73 mother_itr++){
74
75 GenParticlePtr moth;
76 moth = dynamic_cast<PhotosHepMC3Particle*>(*mother_itr)->getHepMC3();
77
78 if(moth->end_vertex()!=orig_production_vertex)
79 Log::Fatal("PhotosHepMC3Particle::setMothers(): Mother production_vertices point to difference places. Can not override. Please delete vertices first.",1);
80 else
81 production_vertex->add_particle_in(moth);
82
83 //update status info
84 if(moth->status()==PhotosParticle::STABLE)
85 {
86 //AV moth->set_status(PhotosParticle::DECAYED);
87 moth->set_status(2);
88 }
89 }
90 production_vertex->add_particle_out(m_particle);
91 }
92}
93
94
95
97
98 //add to this classes internal list as well.
99 m_daughters.push_back(daughter);
100
101 //this assumes there is already an end vertex for the particle
102
103 if(!m_particle->end_vertex())
104 Log::Fatal("PhotosHepMC3Particle::addDaughter(): This method assumes an end_vertex exists. Maybe you really want to use setDaughters.",2);
105
106 GenParticlePtr daugh = (dynamic_cast<PhotosHepMC3Particle*>(daughter))->getHepMC3();
107 m_particle->end_vertex()->add_particle_out(daugh);
108
109}
110
111void PhotosHepMC3Particle::setDaughters(vector<PhotosParticle*> daughters){
112
113 if(!m_particle->parent_event())
114 Log::Fatal("PhotosHepMC3Particle::setDaughters(): New particle needs the event set before it's daughters can be added",3);
115
117
118 //If there are daughters
119 if(daughters.size()>0){
120
121 //Use production vertex of first daughter as end vertex for particle
122 GenParticlePtr first_daughter;
123 first_daughter = (dynamic_cast<PhotosHepMC3Particle*>(daughters.at(0)))->getHepMC3();
124
125 GenVertexPtr end_vertex;
126 end_vertex=first_daughter->production_vertex();
127 GenVertexPtr orig_end_vertex = end_vertex;
128
129 if(!end_vertex){ //if it does not exist create it
130 end_vertex = make_shared<GenVertex>();
131 m_particle->parent_event()->add_vertex(end_vertex);
132 }
133
134 //Loop over all daughters to check that the end points to the right place
135 vector<PhotosParticle*>::iterator daughter_itr;
136 for(daughter_itr = daughters.begin(); daughter_itr != daughters.end();
137 daughter_itr++){
138
139 GenParticlePtr daug;
140 daug = dynamic_cast<PhotosHepMC3Particle*>(*daughter_itr)->getHepMC3();
141
142
143 if(daug->production_vertex()!=orig_end_vertex)
144 Log::Fatal("PhotosHepMC3Particle::setDaughters(): Daughter production_vertices point to difference places. Can not override. Please delete vertices first.",4);
145 else
146 end_vertex->add_particle_out(daug);
147 }
148 end_vertex->add_particle_in(m_particle);
149 }
150
151}
152
153std::vector<PhotosParticle*> PhotosHepMC3Particle::getMothers(){
154
155 if(m_mothers.size()==0&&m_particle->production_vertex()){
156
157 for(auto p: m_particle->production_vertex()->particles_in() ) {
158 m_mothers.push_back(new PhotosHepMC3Particle(p));
159 }
160 }
161 return m_mothers;
162}
163
164std::vector<PhotosParticle*> PhotosHepMC3Particle::getDaughters(){
165
166 if(m_daughters.size()==0&&m_particle->end_vertex()){
167
168 for(auto p: m_particle->end_vertex()->particles_out() ) {
169
170 // ommit particles if their status code is ignored by Photos
171 if( Photos::isStatusCodeIgnored( p->status() ) ) continue;
172
173 m_daughters.push_back(new PhotosHepMC3Particle(p));
174 }
175 }
176 return m_daughters;
177
178}
179
180std::vector<PhotosParticle*> PhotosHepMC3Particle::getAllDecayProducts(){
181
182 m_decay_products.clear();
183
184 if(!hasDaughters()) // if this particle has no daughters
185 return m_decay_products;
186
187 std::vector<PhotosParticle*> daughters = getDaughters();
188
189 // copy daughters to list of all decay products
190 m_decay_products.insert(m_decay_products.end(),daughters.begin(),daughters.end());
191
192 // Now, get all daughters recursively, without duplicates.
193 // That is, for each daughter:
194 // 1) get list of her daughters
195 // 2) for each particle on this list:
196 // a) check if it is already on the list
197 // b) if it's not, add her to the end of the list
198 for(unsigned int i=0;i<m_decay_products.size();i++)
199 {
200 std::vector<PhotosParticle*> daughters2 = m_decay_products[i]->getDaughters();
201
202 if(!m_decay_products[i]->hasDaughters()) continue;
203 for(unsigned int j=0;j<daughters2.size();j++)
204 {
205 bool add=true;
206 for(unsigned int k=0;k<m_decay_products.size();k++)
207 if( daughters2[j]->getBarcode() == m_decay_products[k]->getBarcode() )
208 {
209 add=false;
210 break;
211 }
212
213 if(add) m_decay_products.push_back(daughters2[j]);
214 }
215 }
216 return m_decay_products;
217}
218
220
221 if(!m_particle->end_vertex()) return true;
222
223 // HepMC version of check_momentum_conservation
224 // with added energy check
225 FourVector sum;
226
227 for(ConstGenParticlePtr p: m_particle->end_vertex()->particles_in() ) {
228 if( Photos::isStatusCodeIgnored(p->status()) ) continue;
229
230 sum += p->momentum();
231 }
232
233 for(ConstGenParticlePtr p: m_particle->end_vertex()->particles_out() ) {
234 if( Photos::isStatusCodeIgnored(p->status()) ) continue;
235
236 sum -= p->momentum();
237 }
238
239 if( sum.length() > Photos::momentum_conservation_threshold ) {
240 Log::Warning()<<"Momentum not conserved in the vertex:"<<endl;
241 Log::RedirectOutput(Log::Warning(false));
242 Print::line(m_particle->end_vertex());
244 return false;
245 }
246
247 return true;
248}
249
251 m_particle->set_pid(pdg_id);
252}
253
255 m_particle->set_generated_mass(mass);
256}
257
259 m_particle->set_status(status);
260}
261
263 return m_particle->pid();
264}
265
267 return m_particle->status();
268}
269
271 return m_particle->id();
272}
273
274
276 int pdg_id, int status, double mass,
277 double px, double py, double pz, double e){
278
279 PhotosHepMC3Particle * new_particle = new PhotosHepMC3Particle();
280 new_particle->getHepMC3()->set_pid(pdg_id);
281 new_particle->getHepMC3()->set_status(status);
282 new_particle->getHepMC3()->set_generated_mass(mass);
283
284 FourVector momentum(px,py,pz,e);
285 new_particle->getHepMC3()->set_momentum(momentum);
286
287 m_created_particles.push_back(new_particle);
288 return new_particle;
289}
290
292
293 if(!m_particle->production_vertex())
294 {
295 Log::Warning()<<"PhotosHepMC3Particle::createHistoryEntry(): particle without production vertex."<<endl;
296 return;
297 }
298
299 GenParticlePtr part = make_shared<GenParticle>(*m_particle);
300 part->set_status(Photos::historyEntriesStatus);
301 m_particle->production_vertex()->add_particle_out(part);
302}
303
305{
306 if(m_particle->end_vertex())
307 {
308 Log::Error()<<"PhotosHepMC3Particle::createSelfDecayVertex: particle already has end vertex!"<<endl;
309 return;
310 }
311
312 if(getHepMC3()->parent_event()==NULL)
313 {
314 Log::Error()<<"PhotosHepMC3Particle::createSelfDecayVertex: particle not in the HepMC event!"<<endl;
315 return;
316 }
317
318 // Add new vertex and new particle to HepMC
319 GenParticlePtr outgoing = make_shared<GenParticle>( *(dynamic_cast<PhotosHepMC3Particle*>(out)->m_particle) );
320 GenVertexPtr v = make_shared<GenVertex>();
321
322 // Copy vertex position from parent vertex
323 v->set_position( m_particle->production_vertex()->position() );
324
325 v->add_particle_in (m_particle);
326 v->add_particle_out(outgoing);
327
328 getHepMC3()->parent_event()->add_vertex(v);
329
330 // If this particle was stable, set its status to 2
331 if(getStatus()==1) setStatus(2);
332}
333
335 Print::line(m_particle);
336}
337
338
339/******** Getter and Setter methods: ***********************/
340
342 return m_particle->momentum().px();
343}
344
346 return m_particle->momentum().py();
347}
348
350 return m_particle->momentum().pz();
351}
352
354 return m_particle->momentum().e();
355}
356
358 //make new momentum as something is wrong with
359 //the HepMC momentum setters
360
361 FourVector momentum(m_particle->momentum());
362 momentum.setPx(px);
363 m_particle->set_momentum(momentum);
364}
365
367 FourVector momentum(m_particle->momentum());
368 momentum.setPy(py);
369 m_particle->set_momentum(momentum);
370}
371
372
374 FourVector momentum(m_particle->momentum());
375 momentum.setPz(pz);
376 m_particle->set_momentum(momentum);
377}
378
380 FourVector momentum(m_particle->momentum());
381 momentum.setE(e);
382 m_particle->set_momentum(momentum);
383}
384
386{
387 return m_particle->generated_mass();
388}
389
390} // namespace Photospp
static void Fatal(string text, unsigned short int code=0)
static void RevertOutput()
static void RedirectOutput(void(*func)(), ostream &where=*out)
Definition: Log.cxx:93
void setMothers(std::vector< PhotosParticle * > mothers)
void setDaughters(std::vector< PhotosParticle * > daughters)
std::vector< PhotosParticle * > getAllDecayProducts()
void clear(std::vector< PhotosParticle * > v)
std::vector< PhotosParticle * > getDaughters()
void createSelfDecayVertex(PhotosParticle *out)
void addDaughter(PhotosParticle *daughter)
std::vector< PhotosParticle * > m_created_particles
PhotosHepMC3Particle * createNewParticle(int pdg_id, int status, double mass, double px, double py, double pz, double e)
std::vector< PhotosParticle * > m_decay_products
std::vector< PhotosParticle * > getMothers()
static bool isStatusCodeIgnored(int status)
Definition: Photos.cxx:387
static double momentum_conservation_threshold