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"
7 using namespace std;
8 
9 namespace Photospp
10 {
11 
12 vector<PhotosParticle*> HEPEVT_struct::m_particle_list;
13 int HEPEVT_struct::ME_channel=0;
14 int HEPEVT_struct::decay_idx=0;
15 
16 void HEPEVT_struct::clear(){
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 
45 void HEPEVT_struct::add_particle(int i,PhotosParticle * particle,
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
84  if(Photos::forceMassList)
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 
110 int HEPEVT_struct::set(PhotosBranch *branch)
111 {
112  HEPEVT_struct::clear();
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 
161 void HEPEVT_struct::get(){
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  {
238  PhotosParticle *p = m_particle_list.at(i);
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 
407 void HEPEVT_struct::prepare()
408 {
409  check_ME_channel();
410 }
411 
412 void HEPEVT_struct::complete()
413 {
414 
415 }
416 
417 void HEPEVT_struct::check_ME_channel()
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
vector< PhotosParticle * > getMothers()
Definition: PhotosBranch.h:33
virtual double getPz()=0
vector< PhotosParticle * > getDaughters()
Definition: PhotosBranch.h:36
virtual void setPz(double pz)=0
virtual std::vector< PhotosParticle * > getDaughters()=0
virtual double getE()=0
virtual double getPy()=0
virtual int getPdgID()=0
virtual void setE(double e)=0
virtual double getVirtuality()
virtual void createHistoryEntry()=0
void boostDaughtersToRestFrame(PhotosParticle *boost)
virtual int getStatus()=0
void boostDaughtersFromRestFrame(PhotosParticle *boost)
void boostToRestFrame(PhotosParticle *boost)
virtual void setPx(double px)=0
virtual void setPy(double py)=0
virtual double getMass()=0
virtual PhotosParticle * createNewParticle(int pdg_id, int status, double mass, double px, double py, double pz, double e)=0
virtual void addDaughter(PhotosParticle *daughter)=0
void boostFromRestFrame(PhotosParticle *boost)
virtual double getPx()=0
PhotosParticle * getDecayingParticle()
Definition: PhotosBranch.h:30