280 #include "TauolaHepMCEvent.h"
287 #include "CLHEP/Random/RandomEngine.h"
289 #include "HepMC/GenEvent.h"
290 #include "HepMC/IO_HEPEVT.h"
291 #include "HepMC/HEPEVT_Wrapper.h"
303 for(
int i = 0;
i < *lenv;
i++)
305 *rvec++ = instance->
flat();
324 : fPolarization(
false), fPSet(0), fIsInitialized(
false), fMDTAU(-1), fSelectDecayByEvent(
false)
330 <<
"The RandomNumberProducer module requires the RandomNumberGeneratorService\n"
331 "which appears to be absent. Please add that service to your configuration\n"
332 "or remove the modules that require it." << std::endl;
393 <<
"Attempt to override Tauola an existing ParameterSet\n"
412 <<
"Attempt to initialize Tauola with an empty ParameterSet\n"
420 Tauola::setDecayingParticle(15);
432 fMDTAU = cards.getParameter<
int >(
"mdtau" );
434 if ( fMDTAU == 0 || fMDTAU == 1 )
436 Tauola::setSameParticleDecayMode( cards.getParameter<
int >(
"pjak1" ) ) ;
437 Tauola::setOppositeParticleDecayMode( cards.getParameter<
int >(
"pjak2" ) ) ;
440 Tauola::setTauLifetime(0.0);
455 fPDGs.push_back( Tauola::getDecayingParticle() );
458 Tauola::initialize();
466 if ( fMDTAU != 0 && fMDTAU != 1 )
471 Log::LogWarning(
false);
483 <<
"Attempt to run random number generator of un-initialized Tauola\n"
491 <<
"Attempt to run random number generator of un-initialized Tauola\n"
504 int NPartBefore = evt->particles_size();
505 int NVtxBefore = evt->vertices_size();
522 TauolaHepMCEvent * t_event =
new TauolaHepMCEvent(evt);
530 t_event->decayTaus();
542 for (
int iv=NVtxBefore+1; iv<=evt->vertices_size(); iv++ )
544 HepMC::GenVertex* GenVtx = evt->barcode_to_vertex(-iv);
546 HepMC::GenVertex* ProdVtx = GenPart->production_vertex();
547 HepMC::FourVector PMom = GenPart->momentum();
548 double mass = GenPart->generated_mass();
551 double lifetime = PData->lifetime().value();
555 double ct = -lifetime *
std::log(prob);
556 double VxDec = GenVtx->position().x();
557 VxDec += ct * (PMom.px()/mass);
558 VxDec += ProdVtx->position().x();
559 double VyDec = GenVtx->position().y();
560 VyDec += ct * (PMom.py()/mass);
561 VyDec += ProdVtx->position().y();
562 double VzDec = GenVtx->position().z();
563 VzDec += ct * (PMom.pz()/mass);
564 VzDec += ProdVtx->position().z();
565 double VtDec = GenVtx->position().t();
566 VtDec += ct * (PMom.e()/mass);
567 VtDec += ProdVtx->position().t();
568 GenVtx->set_position( HepMC::FourVector(VxDec,VyDec,VzDec,VtDec) );
574 std::vector<int> BCodes;
576 for (HepMC::GenVertex::particle_iterator pitr= GenVtx->particles_begin(HepMC::children);
577 pitr != GenVtx->particles_end(HepMC::children); ++pitr)
579 if ( (*pitr)->barcode() > 10000 )
581 BCodes.push_back( (*pitr)->barcode() );
584 if ( BCodes.size() > 0 )
586 for (
size_t ibc=0; ibc<BCodes.size(); ibc++ )
589 int nbc = p1->barcode() - 10000 + NPartBefore;
590 p1->suggest_barcode( nbc );
625 if ( mdtau == 101 || mdtau == 201 )
631 Tauola::setSameParticleDecayMode( 1 ) ;
632 Tauola::setOppositeParticleDecayMode( 1 ) ;
636 if ( mdtau == 102 || mdtau == 202 )
642 Tauola::setSameParticleDecayMode( 2 ) ;
643 Tauola::setOppositeParticleDecayMode( 2 ) ;
647 if ( mdtau == 111 || mdtau == 211 )
654 Tauola::setSameParticleDecayMode( 1 ) ;
655 Tauola::setOppositeParticleDecayMode( 0 ) ;
659 if ( mdtau == 112 || mdtau == 212 )
666 Tauola::setSameParticleDecayMode( 2 ) ;
667 Tauola::setOppositeParticleDecayMode( 0 ) ;
671 if ( mdtau == 121 || mdtau == 221 )
678 Tauola::setSameParticleDecayMode( 0 ) ;
679 Tauola::setOppositeParticleDecayMode( 1 ) ;
683 if ( mdtau == 122 || mdtau == 222 )
690 Tauola::setSameParticleDecayMode( 0 ) ;
691 Tauola::setOppositeParticleDecayMode( 2 ) ;
695 if ( mdtau == 140 || mdtau == 240 )
701 Tauola::setSameParticleDecayMode( 3 ) ;
702 Tauola::setOppositeParticleDecayMode( 3 ) ;
706 if ( mdtau == 141 || mdtau == 241 )
713 Tauola::setSameParticleDecayMode( 3 ) ;
714 Tauola::setOppositeParticleDecayMode( 0 ) ;
718 if ( mdtau == 142 || mdtau == 242 )
725 Tauola::setSameParticleDecayMode( 0 ) ;
726 Tauola::setOppositeParticleDecayMode( 3 ) ;
743 for (
int i=0;
i<22;
i++ )
747 if ( sumBra == 0. )
return ;
748 for (
int i=0;
i<22;
i++ )
750 double newBra =
taubra_.gamprt[
i] / sumBra;
751 Tauola::setTauBr(
i+1, newBra );
756 double sumHadronBra = sumBra - sumLeptonBra;
758 for (
int i=0;
i<2;
i++ )
763 for (
int i=2;
i<22;
i++ )
782 Tauola::setSameParticleDecayMode( mode );
785 Tauola::setOppositeParticleDecayMode( mode );
796 Tauola::setSameParticleDecayMode( modeL );
797 Tauola::setOppositeParticleDecayMode( 0 );
805 Tauola::setSameParticleDecayMode( 0 );
806 Tauola::setOppositeParticleDecayMode( modeL );
814 Tauola::setSameParticleDecayMode( modeL );
815 Tauola::setOppositeParticleDecayMode( modeH );
823 Tauola::setSameParticleDecayMode( modeH );
824 Tauola::setOppositeParticleDecayMode( modeL );
832 Tauola::setSameParticleDecayMode( 1 );
833 Tauola::setOppositeParticleDecayMode( modeH );
841 Tauola::setSameParticleDecayMode( modeH );
842 Tauola::setOppositeParticleDecayMode( 1 );
850 Tauola::setSameParticleDecayMode( 2 );
851 Tauola::setOppositeParticleDecayMode( modeH );
859 Tauola::setSameParticleDecayMode( modeH );
860 Tauola::setOppositeParticleDecayMode( 2 );
868 Tauola::setSameParticleDecayMode( modeH );
869 Tauola::setOppositeParticleDecayMode( jaki_.jak2 );
877 Tauola::setSameParticleDecayMode( modeH );
878 Tauola::setOppositeParticleDecayMode( 0 );
886 Tauola::setSameParticleDecayMode( 0 );
887 Tauola::setOppositeParticleDecayMode( modeH );
896 Tauola::setSameParticleDecayMode( 0 );
897 Tauola::setOppositeParticleDecayMode( 0 );
929 if ( prob > 0. && prob <= sumBra )
936 for (
int i=1;
i<NN;
i++ )
952 return (
double)instance->
flat();
T getParameter(std::string const &) const
void selectDecayByMDTAU()
HepMC::GenEvent * decay(HepMC::GenEvent *)
static PFTauRenderPlugin instance
void setPSet(const edm::ParameterSet &)
void rmarin_(int *, int *, int *)
std::vector< double > fScaledLeptonBrRatios
void getData(T &iHolder) const
std::vector< int > fLeptonModes
static TauolaInterface * fInstance
double TauolappInterface_RandGetter()
edm::ESHandle< HepPDT::ParticleDataTable > fPDGTable
virtual CLHEP::HepRandomEngine & getEngine() const =0
Use this to get the random number engine, this is the only function most users should call...
void ranmar_(float *rvec, int *lenv)
HepPDT::ParticleData ParticleData
CLHEP::HepRandomEngine * fRandomEngine
edm::ParameterSet * fPSet
static TauolaInterface * getInstance()
std::vector< double > fScaledHadronBrRatios
void init(const edm::EventSetup &)
std::vector< int > fHadronModes