HEPEVT_struct.cxx
1#include <vector>
2#include <cmath>
3#include "PhotosBranch.h"
4#include "PhotosParticle.h"
5#include "HEPEVT_struct.h"
6#include "Log.h"
7using namespace std;
8
9namespace Photospp
10{
11
12vector<PhotosParticle*> HEPEVT_struct::m_particle_list;
15
17
18 m_particle_list.clear();
19
20 hep.nevhep=0;
21 hep.nhep=0;
22
23
24 /** for(int i=0; i < NMXHEP; i++){
25
26 hep.isthep[i]=0;
27 hep.idhep[i]=0;
28
29 for(int j=0; j<2; j++){
30 hep.jmohep[i][j]=0;
31 hep.jdahep[i][j]=0;
32 }
33
34 for(int j=0; j<5; j++)
35 hep.phep[i][j]=0;
36
37 for(int j=0; j<4; j++)
38 hep.vhep[i][j]=0;
39
40 ph_phoqed_.qedrad[i]=0;
41
42 }**/
43}
44
46 int first_mother, int last_mother,
47 int first_daughter, int last_daughter){
48
49 if(i>0)
50 i--; //account for fortran indicies begining at 1
51 else
52 Log::Warning()<<"Index given to HEPEVT_struct::add_particle "
53 <<"is too low (it must be > 0)."<<endl;
54
55 //add to our internal list of pointer/index pairs
56 m_particle_list.push_back(particle);
57
58 //now set the element of HEPEVT struct
59 hep.nevhep=0; //dummy
60 hep.nhep = hep.nhep + 1; // 1.II.2014: fix for gcc 4.8.1. Previously it was:
61 // hep.nhep = hep.nhep++; which technically is an undefined operation
62 hep.isthep[i]=particle->getStatus();
63 hep.idhep[i]=particle->getPdgID();
64
65 hep.jmohep[i][0]=first_mother;
66 hep.jmohep[i][1]=last_mother;
67
68 hep.jdahep[i][0]=first_daughter;
69 hep.jdahep[i][1]=last_daughter;
70
71 hep.phep[i][0]=particle->getPx();
72 hep.phep[i][1]=particle->getPy();
73 hep.phep[i][2]=particle->getPz();
74 hep.phep[i][3]=particle->getE();
75
76 // if massFrom4Vector=true (default) - get sqrt(e^2-p^2)
77 // otherwise - get mass from event record
78 if(!Photos::massFrom4Vector) hep.phep[i][4]=particle->getMass();
79 else hep.phep[i][4]=particle->getVirtuality();
80
81 int pdgid = abs(particle->getPdgID());
82
83 // if 'forceMass' for this PDGID was used - overwrite mass
85 {
86 for(unsigned int j=0;j<Photos::forceMassList->size();j++)
87 {
88 if(pdgid == abs(Photos::forceMassList->at(j)->first))
89 {
90 double mass = Photos::forceMassList->at(j)->second;
91
92 // when 'forceMass' is used the mass provided is larger than 0.0
93 // when 'forceMassFromEventRecord' is used mass is -1.0
94 // in this case - get mass from event record
95 if(mass<0.0) mass = particle->getMass();
96 hep.phep[i][4] = mass;
97 }
98 }
99 }
100
101 hep.vhep[i][0]=0;
102 hep.vhep[i][1]=0;
103 hep.vhep[i][2]=0;
104 hep.vhep[i][3]=0;
105
106 hep.qedrad[i]=1;
107
108}
109
111{
113 int idx=1;
114
115 //get mothers
116 vector<PhotosParticle *> mothers = branch->getMothers();
117 int nmothers=mothers.size();
118
119 //check if mid-particle exist
120 decay_idx=0;
121 PhotosParticle *decay_particle = branch->getDecayingParticle();
122 if(decay_particle) decay_idx=nmothers+1;
123
124 //get daughters
125 vector<PhotosParticle *> daughters = branch->getDaughters();
126 int ndaughters=daughters.size();
127
128 for(int i=0;i<nmothers;i++)
129 {
130 if(decay_idx)
131 add_particle(idx++,mothers.at(i),
132 0,0, //mothers
133 decay_idx,decay_idx); //daughters
134 else
135 add_particle(idx++,mothers.at(i),
136 0,0, //mothers
137 nmothers+1,nmothers+ndaughters); //daughters
138 }
139
140 if(decay_particle)
141 add_particle(idx++,decay_particle,
142 1,nmothers, //mothers
143 nmothers+2,nmothers+1+ndaughters); //daughters
144
145 for(int i=0;i<ndaughters;i++)
146 {
147 if(decay_idx)
148 add_particle(idx++,daughters.at(i),
149 decay_idx,decay_idx, //mothers
150 0,0); //daughters
151 else
152 add_particle(idx++,daughters.at(i),
153 1,nmothers, //mothers
154 0,0); //daughters
155 }
156 //Log::RedirectOutput( phodmp_ , Log::Debug(1000) );
157 Log::Debug(1000,false)<<"HEPEVT_struct returning: "<<( (decay_idx) ? decay_idx : 1 )<<" from "<<idx-1<<" particles."<<endl;
158 return (decay_idx) ? decay_idx : 1;
159}
160
162
163 int index = 0;
164
165 //if no new particles have been added to the event record, do nothing.
166 if(hep.nhep == (int) m_particle_list.size())
167 return;
168
169 //phodmp_();
170
171 int particle_count = m_particle_list.size();
172 int daughters_start = hep.jmohep[hep.nhep-1][0];
173 int new_particles = hep.nhep - m_particle_list.size();
174 bool isParticleCreated = (new_particles>0);
175
176 std::vector<PhotosParticle*> new_particles_list; // list of added particles
177 // which need kinematical treatment
178 // in special case
179
180 // we decipher daughters_start from last entry
181 // that is last daughter in ph_hepevt_
182 // another option of this functionality may be
183 // hep.jdahep[ hep.jmohep[hep.nhep-1][0]-1][0];
184 // Update daughters_start if there are two mothers
185 // NOTE: daughters_start is index for C++ arrays, while hep.jmohep
186 // contains indices for Fortran arrays.
187 if(hep.jmohep[hep.nhep-1][1]>0)
188 daughters_start = hep.jmohep[hep.nhep-1][1];
189
190 index = particle_count;
191
192 // Add new particles
193 for(;new_particles>0; new_particles--, index++){
194
195 // 27.11.2014: This sanity check is no longer useful (or needed)
196 // We now allow photos to produce particles other than gamma
197 //if(hep.idhep[index]!=PhotosParticle::GAMMA)
198 // Log::Fatal("HEPEVT_struct::get(): Extra particle added to the HEPEVT struct in not a photon!",6);
199
200 //create a new particle
201 PhotosParticle * new_particle;
202 new_particle = m_particle_list.at(0)->createNewParticle(hep.idhep[index],
203 hep.isthep[index],
204 hep.phep[index][4],
205 hep.phep[index][0],
206 hep.phep[index][1],
207 hep.phep[index][2],
208 hep.phep[index][3]);
209
210 //add into the event record
211 //get mother particle of new particle
212 PhotosParticle * mother = m_particle_list.at(hep.jmohep[index][0]-1);
213 mother->addDaughter(new_particle);
214
215 //add to list of new_particles
216 new_particles_list.push_back(new_particle);
217 }
218
219 // Before we update particles, we check for special cases
220 // At this step, particles are yet unmodified
221 // but new particles are already in the event record
222 bool special=false;
223 PhotosParticle *p1 = NULL;
224 PhotosParticle *p2 = NULL;
225
226 if( isParticleCreated )
227 {
228 std::vector<PhotosParticle*> daughters;
229
230 // in the following we create list of daughters,
231 // later we calculate bool special which is true only if all
232 // daughters self-decay
233 // at peresent warning for mixed self-decay and not self decay
234 // daughters is not printed.
235
236 for(int i=daughters_start;i<particle_count;i++)
237 {
239
240 daughters.push_back(p);
241 }
242
243 // Check if this is a special case
244 special = true;
245
246 if(daughters.size()==0) special = false;
247
248 // special = false if there is a stable particle on the list
249 // or there is a particle without self-decay
250 for(unsigned int i=0;i<daughters.size();i++)
251 {
252 if(daughters[i]->getStatus()==1)
253 {
254 special = false;
255 break;
256 }
257
258 // NOTE: We can use 'getDaughters' here, because vertices
259 // of daughters are not being modified by Photos right now
260 // (so there will be no caching)
261 std::vector<PhotosParticle*> daughters2 = daughters[i]->getDaughters();
262
263 if(daughters2.size()!=1 ||
264 daughters2[0]->getPdgID() != daughters[i]->getPdgID() )
265 {
266 special = false;
267 break;
268 }
269 }
270
271 if( special )
272 {
273 double px1=0.0, py1=0.0, pz1=0.0, e1=0.0;
274 double px2=0.0, py2=0.0, pz2=0.0, e2=0.0;
275
276 // get sum of 4-momenta of unmodified particles
277 // 27.11.2014: range changed. Now it does not depend on added particles
278 for(int i=0; i < particle_count-daughters_start; i++)
279 {
280 px1+=daughters[i]->getPx();
281 py1+=daughters[i]->getPy();
282 pz1+=daughters[i]->getPz();
283 e1 +=daughters[i]->getE();
284 }
285
286 // get sum of 4-momenta of particles in self-decay vertices
287 // 27.11.2014: range changed. Now it does not depend on added particles
288 for(int i=0; i < particle_count-daughters_start; i++)
289 {
290 // since 'allDaughtersSelfDecay()' is true
291 // each of these particles has exactly one daughter
292 px2 += daughters[i]->getDaughters().at(0)->getPx();
293 py2 += daughters[i]->getDaughters().at(0)->getPy();
294 pz2 += daughters[i]->getDaughters().at(0)->getPz();
295 e2 += daughters[i]->getDaughters().at(0)->getE();
296 }
297
298 //cout<<"ORIG: "<<px1<<" "<<py1<<" "<<pz1<<" "<<e1<<endl;
299 //cout<<"SELF: "<<px2<<" "<<py2<<" "<<pz2<<" "<<e2<<endl;
300
301 p1 = m_particle_list.at(0)->createNewParticle(0,-1,0.0,px1,py1,pz1,e1);
302 p2 = m_particle_list.at(0)->createNewParticle(0,-2,0.0,px2,py2,pz2,e2);
303
304 // Finally, boost new particles to appropriate frame
305 for(unsigned int i=0;i<new_particles_list.size();i++)
306 {
307 PhotosParticle *boosted = new_particles_list[i]->createNewParticle( 22, 1,
308 0.0,
309 new_particles_list[i]->getPx(),
310 new_particles_list[i]->getPy(),
311 new_particles_list[i]->getPz(),
312 new_particles_list[i]->getE() );
313
314 boosted->boostToRestFrame(p1);
315 boosted->boostFromRestFrame(p2);
316
317 new_particles_list[i]->createSelfDecayVertex(boosted);
318
319 delete boosted;
320 }
321
322 Log::Warning()<<"Hidden interaction, all daughters self decay."
323 <<"Potentially over simplified solution applied."<<endl;
324 }
325 }
326
327 //otherwise loop over particles which are already in the
328 //event record and modify their 4 momentum
329 //4.03.2012: Fix to prevent kinematical trap in vertex of simultaneous:
330 // z-collinear and non-conservation pf E,p for dauthters of grandmothers
331 for(index=daughters_start; index < particle_count && index < (int) m_particle_list.size(); index++){
332
333 PhotosParticle * particle = m_particle_list.at(index);
334
335 if(hep.idhep[index]!=particle->getPdgID())
336 Log::Fatal("HEPEVT_struct::get(): Something is wrong with the HEPEVT struct",5);
337
338 // If new particles were added - for each daughter create a history entry
339 if(isParticleCreated && Photos::isCreateHistoryEntries)
340 {
341 particle->createHistoryEntry();
342 }
343
344 //check to see if this particle's 4-momentum has been modified
345 bool update=false;
346
347 // don't update particle if difference lower than THRESHOLD * particle energy (default threshold = 10e-8)
348 double threshold = NO_BOOST_THRESHOLD*hep.phep[index][3];
349 if( fabs(hep.phep[index][0]-particle->getPx()) > threshold ||
350 fabs(hep.phep[index][1]-particle->getPy()) > threshold ||
351 fabs(hep.phep[index][2]-particle->getPz()) > threshold ||
352 fabs(hep.phep[index][3]-particle->getE()) > threshold ) update=true;
353
354 if(update)
355 {
356
357 //modify this particle's momentum and it's daughters momentum
358 //Steps 1., 2. and 3. must be executed in order.
359
360 //1. boost the particles daughters into it's (old) rest frame
361 particle->boostDaughtersToRestFrame(particle);
362
363 //2. change this particles 4 momentum
364 particle->setPx(hep.phep[index][0]);
365 particle->setPy(hep.phep[index][1]);
366 particle->setPz(hep.phep[index][2]);
367 particle->setE(hep.phep[index][3]);
368
369 //3. boost the particles daughters back into the lab frame
370 particle->boostDaughtersFromRestFrame(particle);
371
372 if(special && particle->getDaughters().size()>0){
373
374 // Algorithm for special case:
375 // a. get self-daughter of 'particle'
376 PhotosParticle *particled = particle->getDaughters().at(0);
377
378 // b. boost 'particled' daughters to rest frame
379 particled->boostDaughtersToRestFrame(particled);
380
381 // c. copy four momentum of 'particle' into four momentum of
382 // its self-daughter 'particled'
383
384 particled->setPx( particle->getPx() );
385 particled->setPy( particle->getPy() );
386 particled->setPz( particle->getPz() );
387 particled->setE ( particle->getE() );
388
389 // d. boost self daughter to rest-frame of <e1>
390 // boost self daughter from rest-frame of <e2>
391
392 particled->boostToRestFrame(p1);
393 particled->boostFromRestFrame(p2);
394
395 // e. boost the 'particled' daughters back into the lab frame
396 particled->boostDaughtersFromRestFrame(particled);
397 }
398
399 }
400 }
401
402 // cleanup
403 if(p1) delete p1;
404 if(p2) delete p2;
405}
406
408{
410}
411
413{
414
415}
416
418{
419 ME_channel=0;
420
421// Check mothers:
422
423 if(decay_idx==2) return; // Only one mother present
424 if(hep.idhep[0]*hep.idhep[1]>0) return; // Mothers have same sign
425
426 Log::Debug(900)<<"ME_channel: Mothers PDG: "<<hep.idhep[0]<<" "<<hep.idhep[1]<<endl;
427 if(decay_idx)
428 Log::Debug(900,false)<<" Intermediate: "<<hep.idhep[decay_idx-1]<<endl;
429
430 int firstDaughter=3;
431 if(decay_idx==0) firstDaughter=2; // if no intermediate particle - daughters start at idx 2
432
433 // Are mothers in range +/- 1-6; +/- 11-16?
434 int mother1 = abs(hep.idhep[0]);
435 int mother2 = abs(hep.idhep[1]);
436 if( mother1<1 || (mother1>6 && mother1<11) || mother1>16 ) return;
437 if( mother2<1 || (mother2>6 && mother2<11) || mother2>16 ) return;
438
439//Check daughters
440
441 // Z: check for pairs 11 -11 ; 13 -13 ; 15 -15
442 // -------------------------------------------
443 int firstPDG =0;
444 int secondPDG=0;
445 for(int i=firstDaughter; i<hep.nhep;i++)
446 {
447 int pdg = abs(hep.idhep[i]);
448 if(pdg==11 || pdg==13 || pdg==15)
449 {
450 if(firstPDG==0) firstPDG=hep.idhep[i];
451 else
452 {
453 secondPDG=hep.idhep[i];
454 // Just in case two pairs are genereted - verify that we have a pair with oposite signs
455 if(firstPDG*secondPDG>0) secondPDG=0;
456 break;
457 }
458 }
459 }
460
461 if( ME_channel==0 && firstPDG!=0 && secondPDG!=0 &&
462 firstPDG==-secondPDG ) ME_channel=1;
463
464 // W: check for pairs 11 -12; -11 12; 13 -14; -13 14; 15 -16; -15 16
465 // -----------------------------------------------------------------
466 firstPDG =0;
467 secondPDG=0;
468 for(int i=firstDaughter; i<hep.nhep;i++)
469 {
470 int pdg = abs(hep.idhep[i]);
471 if(pdg>=11 && pdg<=16)
472 {
473 if(firstPDG==0) firstPDG=hep.idhep[i];
474 else
475 {
476 secondPDG=hep.idhep[i];
477 // Just in case two pairs are genereted - verify that we have a pair with oposite signs
478 if(firstPDG*secondPDG>0) secondPDG=0;
479 break;
480 }
481 }
482 }
483
484 firstPDG =abs(firstPDG);
485 secondPDG=abs(secondPDG);
486
487 if( ME_channel==0 && firstPDG!=0 && secondPDG!=0 &&
488 ( ( firstPDG==11 && secondPDG==12 ) || (firstPDG == 12 && secondPDG == 11) ||
489 ( firstPDG==13 && secondPDG==14 ) || (firstPDG == 14 && secondPDG == 13) ||
490 ( firstPDG==15 && secondPDG==16 ) || (firstPDG == 16 && secondPDG == 15)
491 )
492 ) ME_channel=2;
493
494 Log::Debug(901)<<"ME_channel: Found ME_channel: "<<ME_channel<<endl;
495
496// Check intermediate particle (if exists):
497
498 // Verify that intermediate particle PDG matches ME_channel found
499 if(ME_channel>0 && decay_idx)
500 {
501 int pdg=hep.idhep[decay_idx-1];
502
503 if(ME_channel==1 && !(pdg==22 || pdg==23) ) ME_channel=0; //gamma/Z
504 if(ME_channel==2 && !(pdg==24 || pdg==-24)) ME_channel=0; //W+/W-
505
506 if(ME_channel==0)
507 Log::Debug(901,false)<<" but set to 0: wrong intermediate particle: "<<pdg<<endl;
508 }
509
510// Check flags
511
512 switch(ME_channel)
513 {
514 case 0: break; // Ok - no channel found
515 case 1: if(!Photos::meCorrectionWtForZ) ME_channel=0; break;
516 case 2: if(!Photos::meCorrectionWtForW) ME_channel=0; break;
517 default: Log::Error()<<"Unimplemented ME channel: "<<ME_channel<<endl; break;
518 }
519 Log::Debug(902)<<"ME_channel: Finally, after flag check, ME_channel is: "<<ME_channel<<endl;
520}
521
522
523} // namespace Photospp
static std::vector< PhotosParticle * > m_particle_list
static void add_particle(int i, PhotosParticle *particle, int first_mother, int last_mother, int first_daughter, int last_daughter)
static int set(PhotosBranch *branch)
static ostream & Debug(unsigned short int code=0, bool count=true)
Definition Log.cxx:33
static void Fatal(string text, unsigned short int code=0)
vector< PhotosParticle * > getDaughters()
vector< PhotosParticle * > getMothers()
virtual void setE(double e)=0
virtual std::vector< PhotosParticle * > getDaughters()=0
void boostToRestFrame(PhotosParticle *boost)
virtual void setPy(double py)=0
virtual int getPdgID()=0
virtual double getPz()=0
virtual double getE()=0
virtual double getPy()=0
virtual double getPx()=0
virtual void createHistoryEntry()=0
virtual double getVirtuality()
virtual void setPx(double px)=0
virtual int getStatus()=0
virtual void setPz(double pz)=0
virtual void addDaughter(PhotosParticle *daughter)=0
virtual double getMass()=0
void boostDaughtersFromRestFrame(PhotosParticle *boost)
void boostFromRestFrame(PhotosParticle *boost)
void boostDaughtersToRestFrame(PhotosParticle *boost)
virtual PhotosParticle * createNewParticle(int pdg_id, int status, double mass, double px, double py, double pz, double e)=0
static bool isCreateHistoryEntries
static vector< pair< int, double > * > * forceMassList