172 int daughters_start = hep.jmohep[hep.nhep-1][0];
174 bool isParticleCreated = (new_particles>0);
176 std::vector<PhotosParticle*> new_particles_list;
187 if(hep.jmohep[hep.nhep-1][1]>0)
188 daughters_start = hep.jmohep[hep.nhep-1][1];
190 index = particle_count;
193 for(;new_particles>0; new_particles--, index++){
202 new_particle =
m_particle_list.at(0)->createNewParticle(hep.idhep[index],
213 mother->addDaughter(new_particle);
216 new_particles_list.push_back(new_particle);
226 if( isParticleCreated )
228 std::vector<PhotosParticle*> daughters;
236 for(
int i=daughters_start;i<particle_count;i++)
240 daughters.push_back(p);
246 if(daughters.size()==0) special =
false;
250 for(
unsigned int i=0;i<daughters.size();i++)
252 if(daughters[i]->getStatus()==1)
261 std::vector<PhotosParticle*> daughters2 = daughters[i]->getDaughters();
263 if(daughters2.size()!=1 ||
264 daughters2[0]->getPdgID() != daughters[i]->getPdgID() )
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;
278 for(
int i=0; i < particle_count-daughters_start; i++)
280 px1+=daughters[i]->getPx();
281 py1+=daughters[i]->getPy();
282 pz1+=daughters[i]->getPz();
283 e1 +=daughters[i]->getE();
288 for(
int i=0; i < particle_count-daughters_start; i++)
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();
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);
305 for(
unsigned int i=0;i<new_particles_list.size();i++)
309 new_particles_list[i]->getPx(),
310 new_particles_list[i]->getPy(),
311 new_particles_list[i]->getPz(),
312 new_particles_list[i]->getE() );
314 boosted->boostToRestFrame(p1);
315 boosted->boostFromRestFrame(p2);
317 new_particles_list[i]->createSelfDecayVertex(boosted);
322 Log::Warning()<<
"Hidden interaction, all daughters self decay."
323 <<
"Potentially over simplified solution applied."<<endl;
331 for(index=daughters_start; index < particle_count && index < (int)
m_particle_list.size(); index++){
335 if(hep.idhep[index]!=particle->
getPdgID())
336 Log::Fatal(
"HEPEVT_struct::get(): Something is wrong with the HEPEVT struct",5);
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;
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]);
379 particled->boostDaughtersToRestFrame(particled);
384 particled->setPx( particle->
getPx() );
385 particled->setPy( particle->
getPy() );
386 particled->setPz( particle->
getPz() );
387 particled->setE ( particle->
getE() );
392 particled->boostToRestFrame(p1);
393 particled->boostFromRestFrame(p2);
396 particled->boostDaughtersFromRestFrame(particled);
424 if(hep.idhep[0]*hep.idhep[1]>0)
return;
426 Log::Debug(900)<<
"ME_channel: Mothers PDG: "<<hep.idhep[0]<<
" "<<hep.idhep[1]<<endl;
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;
445 for(
int i=firstDaughter; i<hep.nhep;i++)
447 int pdg = abs(hep.idhep[i]);
448 if(pdg==11 || pdg==13 || pdg==15)
450 if(firstPDG==0) firstPDG=hep.idhep[i];
453 secondPDG=hep.idhep[i];
455 if(firstPDG*secondPDG>0) secondPDG=0;
461 if(
ME_channel==0 && firstPDG!=0 && secondPDG!=0 &&
468 for(
int i=firstDaughter; i<hep.nhep;i++)
470 int pdg = abs(hep.idhep[i]);
471 if(pdg>=11 && pdg<=16)
473 if(firstPDG==0) firstPDG=hep.idhep[i];
476 secondPDG=hep.idhep[i];
478 if(firstPDG*secondPDG>0) secondPDG=0;
484 firstPDG =abs(firstPDG);
485 secondPDG=abs(secondPDG);
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)
507 Log::Debug(901,
false)<<
" but set to 0: wrong intermediate particle: "<<pdg<<endl;
517 default: Log::Error()<<
"Unimplemented ME channel: "<<
ME_channel<<endl;
break;
virtual void setE(double e)=0
virtual std::vector< PhotosParticle * > getDaughters()=0
virtual void setPy(double py)=0
virtual void setPx(double px)=0
virtual void setPz(double pz)=0
void boostDaughtersFromRestFrame(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