10 #include "HepMC/GenEvent.h"
11 #include "HepMC/IO_HEPEVT.h"
12 #include "HepMC/HEPEVT_Wrapper.h"
91 for (
int ip=0; ip<4000; ip++ )
100 for (
int iv=1; iv<=evt->vertices_size(); iv++ )
103 bool legalVtx =
false;
107 HepMC::GenVertex* vtx = evt->barcode_to_vertex( -iv ) ;
109 if ( vtx->particles_in_size() != 1 )
continue;
110 if ( vtx->particles_out_size() <= 1 )
continue;
112 if ( (*(vtx->particles_in_const_begin()))->pdg_id() == 111 )
continue;
114 if (
fOnlyPDG != 1 && (*(vtx->particles_in_const_begin()))->pdg_id() !=
fOnlyPDG )
125 for ( HepMC::GenVertex::particle_iterator pitr=vtx->particles_begin(HepMC::children);
126 pitr != vtx->particles_end(HepMC::children); ++pitr)
128 if ( (*pitr)->pdg_id() ==
fOnlyPDG )
134 if ( same )
continue;
155 for ( HepMC::GenVertex::particle_iterator pitr=vtx->particles_begin(HepMC::children);
156 pitr != vtx->particles_end(HepMC::children); ++pitr)
160 if (
abs((*pitr)->pdg_id()) >=1 &&
abs((*pitr)->pdg_id()) <=8 )
break;
161 if (
abs((*pitr)->pdg_id()) == 21 )
break;
163 if ( (*pitr)->status() == 1 || (*pitr)->end_vertex() )
171 if ( !legalVtx )
continue;
178 HepMC::HEPEVT_Wrapper::set_event_number( evt->event_number() );
187 HepMC::GenVertex* vtx = evt->barcode_to_vertex( vtxbcode );
197 HepMC::HEPEVT_Wrapper::zero_everything();
203 HepMC::HEPEVT_Wrapper::set_id( index, (*(vtx->particles_in_const_begin()))->pdg_id() );
204 HepMC::FourVector
vec4;
205 vec4 = (*(vtx->particles_in_const_begin()))->momentum();
206 HepMC::HEPEVT_Wrapper::set_momentum( index, vec4.x(), vec4.y(), vec4.z(), vec4.e() );
207 HepMC::HEPEVT_Wrapper::set_mass( index, (*(vtx->particles_in_const_begin()))->generated_mass() );
208 HepMC::HEPEVT_Wrapper::set_position( index, vtx->position().x(), vtx->position().y(),
209 vtx->position().z(), vtx->position().t() );
210 HepMC::HEPEVT_Wrapper::set_status( index, (*(vtx->particles_in_const_begin()))->status() );
211 HepMC::HEPEVT_Wrapper::set_parents( index, 0, 0 );
212 fBarcodes.push_back( (*(vtx->particles_in_const_begin()))->barcode() );
217 for ( HepMC::GenVertex::particle_iterator pitr=vtx->particles_begin(HepMC::children);
218 pitr != vtx->particles_end(HepMC::children); ++pitr)
221 if ( (*pitr)->status() == 1 || (*pitr)->end_vertex() )
224 vec4 = (*pitr)->momentum();
225 HepMC::HEPEVT_Wrapper::set_id( index, (*pitr)->pdg_id() );
226 HepMC::HEPEVT_Wrapper::set_momentum( index, vec4.x(), vec4.y(), vec4.z(), vec4.e() );
227 HepMC::HEPEVT_Wrapper::set_mass( index, (*pitr)->generated_mass() );
228 vec4 = (*pitr)->production_vertex()->position();
229 HepMC::HEPEVT_Wrapper::set_position( index, vec4.x(), vec4.y(), vec4.z(), vec4.t() );
230 HepMC::HEPEVT_Wrapper::set_status( index, (*pitr)->status() );
231 HepMC::HEPEVT_Wrapper::set_parents( index, 1, 1 );
232 fBarcodes.push_back( (*pitr)->barcode() );
235 if ( (*pitr)->end_vertex() )
237 fSecVtxStore.push_back( (*pitr)->end_vertex()->barcode() );
243 int nentries =
index;
247 HepMC::HEPEVT_Wrapper::set_children ( index, 2, lastDau );
252 HepMC::HEPEVT_Wrapper::set_number_entries( nentries );
290 unsigned int vcounter = 0;
307 if ( HepMC::HEPEVT_Wrapper::number_entries() > nentries )
318 int largestBarcode = -1;
321 for (
int ip=1; ip<Nbcodes; ip++ )
326 if ( bcode > largestBarcode ) largestBarcode = bcode;
327 double px = HepMC::HEPEVT_Wrapper::px(ip+1);
328 double py = HepMC::HEPEVT_Wrapper::py(ip+1);
329 double pz = HepMC::HEPEVT_Wrapper::pz(ip+1);
333 if ( prt->end_vertex() )
336 HepMC::GenVertex* endVtx = prt->end_vertex();
338 std::vector<int> secVtxStorage;
339 secVtxStorage.clear();
341 secVtxStorage.push_back( endVtx->barcode() );
343 HepMC::FourVector mom4 = prt->momentum();
346 double bet1[3], bet2[3], gam1, gam2, pb;
347 double mass = mom4.m();
348 bet1[0] = -(mom4.px()/mass);
349 bet1[1] = -(mom4.py()/mass);
350 bet1[2] = -(mom4.pz()/mass);
354 gam1 = mom4.e()/mass;
357 unsigned int vcounter = 0;
359 while ( vcounter < secVtxStorage.size() )
362 HepMC::GenVertex* theVtx = evt->barcode_to_vertex( secVtxStorage[vcounter] );
364 for ( HepMC::GenVertex::particle_iterator pitr=theVtx->particles_begin(HepMC::children);
365 pitr != theVtx->particles_end(HepMC::children); ++pitr)
368 if ( (*pitr)->end_vertex() )
370 secVtxStorage.push_back( (*pitr)->end_vertex()->barcode() );
373 if ( theVtx->particles_out_size() == 1 && (*pitr)->pdg_id() == prt->pdg_id() )
376 (*pitr)->set_momentum( HepMC::FourVector(px,py,pz,e) );
380 HepMC::FourVector dmom4 = (*pitr)->momentum();
383 pb = bet1[0]*dmom4.px() + bet1[1]*dmom4.py() + bet1[2]*dmom4.pz();
384 double dpx = dmom4.px() + bet1[0] * (dmom4.e() + pb/(gam1+1.) );
385 double dpy = dmom4.py() + bet1[1] * (dmom4.e() + pb/(gam1+1.) );
386 double dpz = dmom4.pz() + bet1[2] * (dmom4.e() + pb/(gam1+1.) );
387 double de = gam1*dmom4.e() + pb;
389 pb = bet2[0]*dpx + bet2[1]*dpy + bet2[2]*dpz;
390 dpx += bet2[0] * ( de + pb/(gam2+1.) );
391 dpy += bet2[1] * ( de + pb/(gam2+1.) );
392 dpz += bet2[2] * ( de + pb/(gam2+1.) );
396 (*pitr)->set_momentum( HepMC::FourVector(dpx,dpy,dpz,de) );
402 secVtxStorage.clear();
405 prt->set_momentum( HepMC::FourVector(px,py,pz,e) );
409 int newlyGen = HepMC::HEPEVT_Wrapper::number_entries() - nentries;
411 if ( largestBarcode < evt->particles_size() )
417 for (
int ipp=evt->particles_size(); ipp>largestBarcode; ipp-- )
419 (evt->barcode_to_particle(ipp))->suggest_barcode( ipp+newlyGen );
425 for (
int ipnw=1; ipnw<=newlyGen; ipnw++ )
427 int nbcode = largestBarcode+ipnw;
428 int pdg_id = HepMC::HEPEVT_Wrapper::id( nentries+ipnw );
430 double px = HepMC::HEPEVT_Wrapper::px( nentries+ipnw );
431 double py = HepMC::HEPEVT_Wrapper::py( nentries+ipnw );
432 double pz = HepMC::HEPEVT_Wrapper::pz( nentries+ipnw );
438 NewPart->set_generated_mass( m );
439 NewPart->suggest_barcode( nbcode );
440 vtx->add_particle_out( NewPart ) ;
454 if (
abs((*(vtx->particles_in_const_begin()))->pdg_id()) != 15 )
return false;
456 for ( HepMC::GenVertex::particle_iterator pitr=vtx->particles_begin(HepMC::children);
457 pitr != vtx->particles_end(HepMC::children); ++pitr)
459 if (
abs((*pitr)->pdg_id()) == 11 ||
abs((*pitr)->pdg_id()) == 13 )
void attachParticles(HepMC::GenEvent *, HepMC::GenVertex *, int)
double phoran_(int *idummy)
std::vector< int > fBarcodes
bool fAvoidTauLeptonicDecays
Abs< T >::type abs(const T &t)
void applyToBranch(HepMC::GenEvent *, int)
std::vector< std::string > fSpecialSettings
std::vector< int > fSecVtxStore
void applyToVertex(HepMC::GenEvent *, int)
CLHEP::HepRandomEngine * decayRandomEngine
bool isTauLeptonicDecay(HepMC::GenVertex *)
HepMC::GenEvent * apply(HepMC::GenEvent *)
void configureOnlyFor(int)