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