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