PhotosHEPEVTParticle.cxx
1#include "PhotosHEPEVTParticle.h"
2#include "Log.h"
3
4namespace Photospp
5{
6
8{
9 // Cleanup particles that do not belong to event
10 for(unsigned int i=0;i<cache.size();i++)
11 if(cache[i]->m_barcode<0)
12 delete cache[i];
13}
14
15PhotosHEPEVTParticle::PhotosHEPEVTParticle(int pdgid, int status, double px, double py, double pz, double e, double m, int ms, int me, int ds, int de){
16 m_px = px;
17 m_py = py;
18 m_pz = pz;
19 m_e = e;
21
22 m_pdgid = pdgid;
23 m_status = status;
24
25 m_first_mother = ms;
26 m_second_mother = me;
28 m_daughter_end = de;
29
30 m_barcode = -1;
31 m_event = NULL;
32}
33
34/** Add a new daughter to this particle */
36{
37 if(!m_event) Log::Fatal("PhotosHEPEVTParticle::addDaughter - particle not in event record");
38
39 std::vector<PhotosParticle*> mothers = daughter->getMothers();
40
41 mothers.push_back( (PhotosParticle*)this );
42
43 daughter->setMothers(mothers);
44
45 int bc = daughter->getBarcode();
46
47 if(m_daughter_end < 0)
48 {
50 m_daughter_end = bc;
51 }
52 // if it's in the middle of the event record
53 else if(m_daughter_end != bc-1)
54 {
55 PhotosHEPEVTParticle *newPart = m_event->getParticle(bc);
56
57 // Move all particles one spot down the list, to make place for new particle
58 for(int i=bc-1;i>m_daughter_end;i--)
59 {
60 PhotosHEPEVTParticle *move = m_event->getParticle(i);
61 move->setBarcode(i+1);
62 m_event->setParticle(i+1,move);
63 }
64
65 m_daughter_end++;
66 newPart->setBarcode(m_daughter_end);
67 m_event->setParticle(m_daughter_end,newPart);
68
69 // Now: correct all pointers before new particle
70 for(int i=0;i<m_daughter_end;i++)
71 {
72 PhotosHEPEVTParticle *check = m_event->getParticle(i);
73 int m = check->getDaughterRangeEnd();
74 if(m!=-1 && m>m_daughter_end)
75 {
76 check->setDaughterRangeEnd(m+1);
78 }
79 }
80 }
81 else m_daughter_end = bc;
82}
83
85
86 // If this particle has not yet been added to the event record
87 // then add it to the mothers' event record
88 if(m_barcode<0 && mothers.size()>0)
89 {
90 PhotosHEPEVTEvent *evt = ((PhotosHEPEVTParticle*)mothers[0])->m_event;
91 evt->addParticle(this);
92 }
93
94 if(mothers.size()>2) Log::Fatal("PhotosHEPEVTParticle::setMothers: HEPEVT does not allow more than two mothers!");
95
96 if(mothers.size()>0) m_first_mother = mothers[0]->getBarcode();
97 if(mothers.size()>1) m_second_mother = mothers[1]->getBarcode();
98}
99
101
102 // This particle must be inside some event record to be able to add daughters
103 if(m_event==NULL) Log::Fatal("PhotosHEPEVTParticle::setDaughters: particle not inside event record.");
104
105 int beg = 65535, end = -1;
106
107 for(unsigned int i=0;i<daughters.size();i++)
108 {
109 int bc = daughters[i]->getBarcode();
110 if(bc<0) Log::Fatal("PhotosHEPEVTParticle::setDaughters: all daughters has to be in event record first");
111 if(bc<beg) beg = bc;
112 if(bc>end) end = bc;
113 }
114 if(end == -1) beg = -1;
115
116 m_daughter_start = beg;
117 m_daughter_end = end;
118}
119
120std::vector<PhotosParticle*> PhotosHEPEVTParticle::getMothers(){
121
122 std::vector<PhotosParticle*> mothers;
123
124 PhotosParticle *p1 = NULL;
125 PhotosParticle *p2 = NULL;
126
127 // 21.XI.2013: Some generators set same mother ID in both indices if there is only one mother
128 if(m_first_mother == m_second_mother) m_second_mother = -1;
129
130 if(m_first_mother>=0) p1 = m_event->getParticle(m_first_mother);
131 if(m_second_mother>=0) p2 = m_event->getParticle(m_second_mother);
132
133 if(p1) mothers.push_back(p1);
134 if(p2) mothers.push_back(p2);
135
136 return mothers;
137}
138
139// WARNING: this method also corrects daughter indices
140// if such were not defined
141std::vector<PhotosParticle*> PhotosHEPEVTParticle::getDaughters(){
142
143 std::vector<PhotosParticle*> daughters;
144
145 if(!m_event) return daughters;
146
147 // Check if m_daughter_start and m_daughter_end are set
148 // If not - try to get list of daughters from event
149 if(m_daughter_end<0)
150 {
151 int min_d=65535, max_d=-1;
152 for(int i=0;i<m_event->getParticleCount();i++)
153 {
154 if(m_event->getParticle(i)->isDaughterOf(this))
155 {
156 if(i<min_d) min_d = i;
157 if(i>max_d) max_d = i;
158 }
159 }
160 if(max_d>=0)
161 {
162 m_daughter_start = min_d;
163 m_daughter_end = max_d;
164 m_status = 2;
165 }
166 }
167
168 // If m_daughter_end is still not set - there are no daughters
169 // Otherwsie - get daughters
170 if(m_daughter_end>=0)
171 {
172 for(int i=m_daughter_start;i<=m_daughter_end;i++)
173 {
174 PhotosParticle *p = m_event->getParticle(i);
175 if(p==NULL)
176 {
177 Log::Warning()<<"PhotosHEPEVTParticle::getDaughters(): No particle with index "<<i<<endl;
178 return daughters;
179 }
180
181 daughters.push_back(p);
182 }
183 }
184
185 return daughters;
186}
187
188std::vector<PhotosParticle*> PhotosHEPEVTParticle::getAllDecayProducts()
189{
190 std::vector<PhotosParticle*> list;
191
192 if(!hasDaughters()) // if this particle has no daughters
193 return list;
194
195 std::vector<PhotosParticle*> daughters = getDaughters();
196
197 // copy daughters to list of all decay products
198 list.insert(list.end(),daughters.begin(),daughters.end());
199
200 // Now, get all daughters recursively, without duplicates.
201 // That is, for each daughter:
202 // 1) get list of her daughters
203 // 2) for each particle on this list:
204 // a) check if it is already on the list
205 // b) if it's not, add her to the end of the list
206 for(unsigned int i=0;i<list.size();i++)
207 {
208 std::vector<PhotosParticle*> daughters2 = list[i]->getDaughters();
209
210 if(!list[i]->hasDaughters()) continue;
211 for(unsigned int j=0;j<daughters2.size();j++)
212 {
213 bool add=true;
214 for(unsigned int k=0;k<list.size();k++)
215 if( daughters2[j]->getBarcode() == list[k]->getBarcode() )
216 {
217 add=false;
218 break;
219 }
220
221 if(add) list.push_back(daughters2[j]);
222 }
223 }
224
225 return list;
226}
227
229
230 if(!m_event) return true;
231 if(m_daughter_end < 0) return true;
232
233 PhotosHEPEVTParticle *buf = m_event->getParticle(m_daughter_start);
234
235 int first_mother_idx = buf->getFirstMotherIndex();
236 int second_mother_idx = buf->getSecondMotherIndex();
237
238 double px =0.0, py =0.0, pz =0.0, e =0.0;
239 double px2=0.0, py2=0.0, pz2=0.0, e2=0.0;
240
241 for(int i=m_daughter_start;i<=m_daughter_end;i++)
242 {
243 buf = m_event->getParticle(i);
244 px += buf->getPx();
245 py += buf->getPy();
246 pz += buf->getPz();
247 e += buf->getE ();
248 }
249
250 if(first_mother_idx>=0)
251 {
252 buf = m_event->getParticle(first_mother_idx);
253 px2 += buf->getPx();
254 py2 += buf->getPy();
255 pz2 += buf->getPz();
256 e2 += buf->getE();
257 }
258
259 if(second_mother_idx>=0)
260 {
261 buf = m_event->getParticle(second_mother_idx);
262 px2 += buf->getPx();
263 py2 += buf->getPy();
264 pz2 += buf->getPz();
265 e2 += buf->getE();
266 }
267 // 3-momentum // test HepMC style
268 double dp = sqrt( (px-px2)*(px-px2) + (py-py2)*(py-py2) + (pz-pz2)*(pz-pz2) );
269 // virtuality test as well.
270 double m1 = sqrt( fabs( e*e - px*px - py*py - pz*pz ) );
271 double m2 = sqrt( fabs( e2*e2 - px2*px2 - py2*py2 - pz2*pz2 ) );
272
273 if( fabs(m1-m2) > 0.0001 || dp > 0.0001*(e+e2))
274 {
275 Log::RedirectOutput( Log::Warning()<<"Momentum not conserved in vertex: " );
276 if(first_mother_idx >=0) m_event->getParticle(first_mother_idx) ->print();
277 if(second_mother_idx>=0) m_event->getParticle(second_mother_idx)->print();
278 cout<<" TO: "<<endl;
279 for(int i=m_daughter_start;i<=m_daughter_end;i++) m_event->getParticle(i)->print();
281 return false;
282 }
283
284 return true;
285}
286
288 int pdg_id, int status, double mass,
289 double px, double py, double pz, double e){
290
291 // New particles created using this method are added to cache
292 // They will be deleted when this particle will be deleted
293
294 cache.push_back(new PhotosHEPEVTParticle(pdg_id,status,px,py,pz,e,mass,-1,-1,-1,-1));
295 return cache.back();
296}
297
299{
300 Log::Warning()<<"PhotosParticle::createHistoryEntry() not implemented for HEPEVT."<<endl;
301}
302
304{
305 Log::Warning()<<"PhotosHEPEVTParticle::createSelfDecayVertex() not implemented for HEPEVT."<<endl;
306}
307
309{
310 int bc = p->getBarcode();
311 if(bc==m_first_mother || bc==m_second_mother) return true;
312
313 return false;
314}
315
317{
318 int bc = p->getBarcode();
319 if(bc>=m_daughter_start && bc<=m_daughter_end) return true;
320
321 return false;
322}
323
325 char buf[256];
326 sprintf(buf,"P: (%2i) %6i %2i | %11.4e %11.4e %11.4e %11.4e | %11.4e | M: %2i %2i | D: %2i %2i\n",
327 m_barcode, m_pdgid, m_status, m_px, m_py, m_pz, m_e, m_generated_mass,
328 m_first_mother, m_second_mother, m_daughter_start, m_daughter_end);
329
330 cout<<buf;
331}
332
333/******** Getter and Setter methods: ***********************/
334
336 m_pdgid = pdg_id;
337}
338
340 m_status = status;
341}
342
344 m_generated_mass = mass;
345}
346
350
354
358
360 return m_px;
361}
362
364 return m_py;
365}
366
368 return m_pz;
369}
370
372 return m_e;
373}
374
376 m_px = px;
377}
378
380 m_py = py;
381}
382
383
385 m_pz = pz;
386}
387
389 m_e = e;
390}
391
395
397 m_barcode = barcode;
398}
399
403
407
409 return m_second_mother;
410}
411
415
417 return m_daughter_end;
418}
419
420} // 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 addParticle(PhotosHEPEVTParticle *p)
void setMothers(std::vector< PhotosParticle * > mothers)
void setDaughters(std::vector< PhotosParticle * > daughters)
bool isMotherOf(PhotosHEPEVTParticle *p)
PhotosHEPEVTParticle(int pdgid, int status, double px, double py, double pz, double e, double m, int ms, int me, int ds, int de)
std::vector< PhotosParticle * > getAllDecayProducts()
void setEvent(PhotosHEPEVTEvent *event)
std::vector< PhotosParticle * > getDaughters()
void createSelfDecayVertex(PhotosParticle *out)
void addDaughter(PhotosParticle *daughter)
bool isDaughterOf(PhotosHEPEVTParticle *p)
PhotosHEPEVTParticle * createNewParticle(int pdg_id, int status, double mass, double px, double py, double pz, double e)
std::vector< PhotosParticle * > getMothers()
virtual void setMothers(std::vector< PhotosParticle * > mothers)=0
virtual std::vector< PhotosParticle * > getMothers()=0
virtual int getBarcode()=0
STL class.
STL class.