7 #include "Tauola/Tauola.h"
8 #include "Tauola/TauolaHepMCEvent.h"
9 #include "Tauola/Log.h"
10 #include "Tauola/TauolaHepMCParticle.h"
11 #include "Tauola/TauolaParticle.h"
16 #include "CLHEP/Random/RandomEngine.h"
18 #include "HepMC/GenEvent.h"
19 #include "HepMC/IO_HEPEVT.h"
20 #include "HepMC/HEPEVT_Wrapper.h"
42 for(
int i = 0;
i < *lenv;
i++)
57 fIsInitialized(
false),
59 fSelectDecayByEvent(
false),
63 dolheBosonCorr(
false),
66 if (
fPSet != 0 )
throw cms::Exception(
"TauolappInterfaceError") <<
"Attempt to override Tauola an existing ParameterSet\n" << std::endl;
74 if (
fPSet == 0 )
throw cms::Exception(
"TauolappInterfaceError") <<
"Attempt to initialize Tauola with an empty ParameterSet\n" << std::endl;
80 Tauolapp::Tauola::setDecayingParticle(15);
96 fMDTAU = cards.getParameter<
int >(
"mdtau" );
98 if ( fMDTAU == 0 || fMDTAU == 1 )
100 Tauolapp::Tauola::setSameParticleDecayMode( cards.getParameter<
int >(
"pjak1" ) ) ;
101 Tauolapp::Tauola::setOppositeParticleDecayMode( cards.getParameter<
int >(
"pjak2" ) ) ;
104 Tauolapp::Tauola::spin_correlation.setAll(fPolarization);
111 double lifetime = PData->lifetime().value();
112 Tauolapp::Tauola::setTauLifetime( lifetime );
114 fPDGs.push_back( Tauolapp::Tauola::getDecayingParticle() );
119 std::vector<std::string> par =
fPSet->
getParameter< std::vector<std::string> >(
"parameterSets");
120 for (
unsigned int ip=0; ip<par.size(); ++ip ){
122 if(curSet==
"setNewCurrents") Tauolapp::Tauola::setNewCurrents(
fPSet->
getParameter<
int>(curSet));
128 Tauolapp::Tauola::spin_correlation.setAll(fPolarization);
131 std::vector<std::string> par =
fPSet->
getParameter< std::vector<std::string> >(
"parameterSets");
132 for (
unsigned int ip=0; ip<par.size(); ++ip ){
134 if(curSet==
"spinCorrelationSetAll") Tauolapp::Tauola::spin_correlation.setAll(
fPSet->
getParameter<
bool>(curSet));
135 if(curSet==
"spinCorrelationGAMMA") Tauolapp::Tauola::spin_correlation.GAMMA=
fPSet->
getParameter<
bool>(curSet);
136 if(curSet==
"spinCorrelationZ0") Tauolapp::Tauola::spin_correlation.Z0=
fPSet->
getParameter<
bool>(curSet);
137 if(curSet==
"spinCorrelationHIGGS") Tauolapp::Tauola::spin_correlation.HIGGS=
fPSet->
getParameter<
bool>(curSet);
138 if(curSet==
"spinCorrelationHIGGSH") Tauolapp::Tauola::spin_correlation.HIGGS_H=
fPSet->
getParameter<
bool>(curSet);
139 if(curSet==
"spinCorrelationHIGGSA") Tauolapp::Tauola::spin_correlation.HIGGS_A=
fPSet->
getParameter<
bool>(curSet);
140 if(curSet==
"spinCorrelationHIGGSPLUS") Tauolapp::Tauola::spin_correlation.HIGGS_PLUS=
fPSet->
getParameter<
bool>(curSet);
141 if(curSet==
"spinCorrelationHIGGSMINUS") Tauolapp::Tauola::spin_correlation.HIGGS_MINUS=
fPSet->
getParameter<
bool>(curSet);
142 if(curSet==
"spinCorrelationWPLUS") Tauolapp::Tauola::spin_correlation.W_PLUS=
fPSet->
getParameter<
bool>(curSet);
143 if(curSet==
"spinCorrelationWMINUS") Tauolapp::Tauola::spin_correlation.W_MINUS=
fPSet->
getParameter<
bool>(curSet);
145 if(curSet==
"setHiggsScalarPseudoscalarPDG") Tauolapp::Tauola::setHiggsScalarPseudoscalarPDG(
fPSet->
getParameter<
int>(curSet));
146 if(curSet==
"setHiggsScalarPseudoscalarMixingAngle") Tauolapp::Tauola::setHiggsScalarPseudoscalarMixingAngle(
fPSet->
getParameter<
double>(curSet));
148 if(curSet==
"setRadiation") Tauolapp::Tauola::setRadiation(
fPSet->
getParameter<
bool>(curSet));
149 if(curSet==
"setRadiationCutOff") Tauolapp::Tauola::setRadiationCutOff(
fPSet->
getParameter<
double>(curSet));
151 if(curSet==
"setEtaK0sPi"){
153 if(vpar.size()==3) Tauolapp::Tauola::setEtaK0sPi(vpar[0],vpar[1],vpar[2]);
154 else {
std::cout <<
"WARNING invalid size for setEtaK0sPi: " << vpar.size() <<
" Require 3 elements " << std::endl;}
157 if(curSet==
"setTaukle"){
159 if(vpar.size()==4) Tauolapp::Tauola::setTaukle(vpar[0], vpar[1], vpar[2], vpar[3]);
160 else {
std::cout <<
"WARNING invalid size for setTaukle: " << vpar.size() <<
" Require 4 elements " << std::endl;}
163 if(curSet==
"setTauBr"){
165 std::vector<int> vJAK = cfg.
getParameter<std::vector<int> >(
"JAK");
166 std::vector<double> vBR = cfg.
getParameter<std::vector<double> >(
"BR");
167 if(vJAK.size() == vBR.size()){
168 for(
unsigned int i=0;
i<vJAK.size();
i++) Tauolapp::Tauola::setTauBr(vJAK[
i],vBR[i]);
170 else {
std::cout <<
"WARNING invalid size for setTauBr - JAK: " << vJAK.size() <<
" BR: " << vBR.size() << std::endl;}
179 if ( fMDTAU != 0 && fMDTAU != 1 )
184 Tauolapp::Log::LogWarning(
false);
193 <<
"TauolaInterface::flat: Attempt to generate random number when engine pointer is null\n"
194 <<
"This might mean that the code was modified to generate a random number outside the\n"
195 <<
"event and beginLuminosityBlock methods, which is not allowed.\n";
203 int NPartBefore = evt->particles_size();
204 int NVtxBefore = evt->vertices_size();
215 std::vector<HepMC::GenParticle> particles;
216 std::vector<int> m_idx;
219 for(
unsigned int i=0;
i<spinup.size();
i++){
222 particles.at(particles.size()-1).set_generated_mass(
lhe->
getHEPEUP()->
PUP.at(
i)[4]);
223 particles.at(particles.size()-1).set_status(status > 0 ? (status == 2 ? 3 : status) : 3);
228 std::vector<HepMC::GenParticle*>
match;
229 for(HepMC::GenEvent::particle_const_iterator
iter = evt->particles_begin();
iter != evt->particles_end();
iter++) {
230 if(
abs((*iter)->pdg_id())==15){
234 for(HepMC::GenVertex::particle_iterator mother=(*iter)->production_vertex()->particles_begin(
HepMC::parents); mother!=(*iter)->production_vertex()->particles_end(
HepMC::parents); mother++) {
235 mother_pid = (*mother)->pdg_id();
236 if(mother_pid!=(*iter)->pdg_id()){
238 if(
abs(mother_pid) == 24 ||
239 abs(mother_pid) == 37 ||
240 abs(mother_pid) == 23 ||
241 abs(mother_pid) == 22 ||
242 abs(mother_pid) == 25 ||
243 abs(mother_pid) == 35 ||
244 abs(mother_pid) == 36
247 for(
unsigned int k=0;
k<match.size();
k++){
248 if((*mother)==match.at(
k))isfound=
true;
250 if(!isfound) match.push_back(*mother);
261 auto * t_event =
new Tauolapp::TauolaHepMCEvent(evt);
262 t_event->undecayTaus();
263 t_event->decayTaus();
265 for(
unsigned int j=0;
j<spinup.size();
j++){
266 if(
abs(pdg.at(
j))==15){
267 double diffhelminus=(-1.0*(double)Tauolapp::Tauola::getHelMinus()-spinup.at(
j));
268 double diffhelplus=((double)Tauolapp::Tauola::getHelPlus()-spinup.at(
j));
269 if(pdg.at(
j)==15 && diffhelminus>0.5) ismatch=
false;
270 if(pdg.at(
j)==-15 && diffhelplus>0.5) ismatch=
false;
280 auto * t_event =
new Tauolapp::TauolaHepMCEvent(evt);
281 t_event->undecayTaus();
284 for(HepMC::GenEvent::particle_const_iterator
iter = evt->particles_begin();
iter != evt->particles_end();
iter++) {
286 TLorentzVector ltau((*iter)->momentum().px(),(*iter)->momentum().py(),(*iter)->momentum().pz(),(*iter)->momentum().e());
288 TLorentzVector mother(m->momentum().px(),m->momentum().py(),m->momentum().pz(),m->momentum().e());
289 TVector3 boost=-1.0*mother.BoostVector();
290 TLorentzVector ltau_lab=ltau;
295 Tauolapp::TauolaParticle *tp =
new Tauolapp::TauolaHepMCParticle(p);
297 if((*iter)->pdg_id()==15) helicity*=-1.0;
300 Tauolapp::Tauola::decayOne(tp,
true,0,0,((
double)helicity)*0.999999);
314 auto * t_event =
new Tauolapp::TauolaHepMCEvent(evt);
316 t_event->decayTaus();
320 for (
int iv=NVtxBefore+1; iv<=evt->vertices_size(); iv++ ){
321 HepMC::GenVertex* GenVtx = evt->barcode_to_vertex(-iv);
327 std::vector<int> BCodes;
329 for (HepMC::GenVertex::particle_iterator pitr= GenVtx->particles_begin(HepMC::children);
330 pitr != GenVtx->particles_end(HepMC::children); ++pitr){
331 if ( (*pitr)->barcode() > 10000 ){
332 BCodes.push_back( (*pitr)->barcode() );
335 if ( BCodes.size() > 0 ){
336 for (
size_t ibc=0; ibc<BCodes.size(); ibc++ ){
338 int nbc = p1->barcode() - 10000 + NPartBefore;
339 p1->suggest_barcode( nbc );
372 if ( mdtau == 101 || mdtau == 201 )
376 Tauolapp::jaki_.jak1 = 1;
377 Tauolapp::jaki_.jak2 = 1;
378 Tauolapp::Tauola::setSameParticleDecayMode( 1 ) ;
379 Tauolapp::Tauola::setOppositeParticleDecayMode( 1 ) ;
383 if ( mdtau == 102 || mdtau == 202 )
387 Tauolapp::jaki_.jak1 = 2;
388 Tauolapp::jaki_.jak2 = 2;
389 Tauolapp::Tauola::setSameParticleDecayMode( 2 ) ;
390 Tauolapp::Tauola::setOppositeParticleDecayMode( 2 ) ;
394 if ( mdtau == 111 || mdtau == 211 )
399 Tauolapp::jaki_.jak1 = 1;
400 Tauolapp::jaki_.jak2 = 0;
401 Tauolapp::Tauola::setSameParticleDecayMode( 1 ) ;
402 Tauolapp::Tauola::setOppositeParticleDecayMode( 0 ) ;
406 if ( mdtau == 112 || mdtau == 212 )
411 Tauolapp::jaki_.jak1 = 2;
412 Tauolapp::jaki_.jak2 = 0;
413 Tauolapp::Tauola::setSameParticleDecayMode( 2 ) ;
414 Tauolapp::Tauola::setOppositeParticleDecayMode( 0 ) ;
418 if ( mdtau == 121 || mdtau == 221 )
423 Tauolapp::jaki_.jak1 = 0;
424 Tauolapp::jaki_.jak2 = 1;
425 Tauolapp::Tauola::setSameParticleDecayMode( 0 ) ;
426 Tauolapp::Tauola::setOppositeParticleDecayMode( 1 ) ;
430 if ( mdtau == 122 || mdtau == 222 )
435 Tauolapp::jaki_.jak1 = 0;
436 Tauolapp::jaki_.jak2 = 2;
437 Tauolapp::Tauola::setSameParticleDecayMode( 0 ) ;
438 Tauolapp::Tauola::setOppositeParticleDecayMode( 2 ) ;
442 if ( mdtau == 140 || mdtau == 240 )
446 Tauolapp::jaki_.jak1 = 3;
447 Tauolapp::jaki_.jak2 = 3;
448 Tauolapp::Tauola::setSameParticleDecayMode( 3 ) ;
449 Tauolapp::Tauola::setOppositeParticleDecayMode( 3 ) ;
453 if ( mdtau == 141 || mdtau == 241 )
458 Tauolapp::jaki_.jak1 = 3;
459 Tauolapp::jaki_.jak2 = 0;
460 Tauolapp::Tauola::setSameParticleDecayMode( 3 ) ;
461 Tauolapp::Tauola::setOppositeParticleDecayMode( 0 ) ;
465 if ( mdtau == 142 || mdtau == 242 )
470 Tauolapp::jaki_.jak1 = 0;
471 Tauolapp::jaki_.jak2 = 3;
472 Tauolapp::Tauola::setSameParticleDecayMode( 0 ) ;
473 Tauolapp::Tauola::setOppositeParticleDecayMode( 3 ) ;
490 for (
int i=0;
i<22;
i++ )
494 if ( sumBra == 0. )
return ;
495 for (
int i=0;
i<22;
i++ )
498 Tauolapp::Tauola::setTauBr(
i+1, newBra );
503 double sumHadronBra = sumBra - sumLeptonBra;
505 for (
int i=0;
i<2;
i++ )
510 for (
int i=2;
i<22;
i++ )
528 Tauolapp::jaki_.jak1 =
mode;
529 Tauolapp::Tauola::setSameParticleDecayMode( mode );
531 Tauolapp::jaki_.jak2 =
mode;
532 Tauolapp::Tauola::setOppositeParticleDecayMode( mode );
541 Tauolapp::jaki_.jak1 = modeL;
542 Tauolapp::jaki_.jak2 = 0;
543 Tauolapp::Tauola::setSameParticleDecayMode( modeL );
544 Tauolapp::Tauola::setOppositeParticleDecayMode( 0 );
550 Tauolapp::jaki_.jak1 = 0;
551 Tauolapp::jaki_.jak2 = modeL;
552 Tauolapp::Tauola::setSameParticleDecayMode( 0 );
553 Tauolapp::Tauola::setOppositeParticleDecayMode( modeL );
559 Tauolapp::jaki_.jak1 = modeL;
560 Tauolapp::jaki_.jak2 = modeH;
561 Tauolapp::Tauola::setSameParticleDecayMode( modeL );
562 Tauolapp::Tauola::setOppositeParticleDecayMode( modeH );
568 Tauolapp::jaki_.jak1 = modeH;
569 Tauolapp::jaki_.jak2 = modeL;
570 Tauolapp::Tauola::setSameParticleDecayMode( modeH );
571 Tauolapp::Tauola::setOppositeParticleDecayMode( modeL );
577 Tauolapp::jaki_.jak1 = 1;
578 Tauolapp::jaki_.jak2 = modeH;
579 Tauolapp::Tauola::setSameParticleDecayMode( 1 );
580 Tauolapp::Tauola::setOppositeParticleDecayMode( modeH );
586 Tauolapp::jaki_.jak1 = modeH;
587 Tauolapp::jaki_.jak2 = 1;
588 Tauolapp::Tauola::setSameParticleDecayMode( modeH );
589 Tauolapp::Tauola::setOppositeParticleDecayMode( 1 );
595 Tauolapp::jaki_.jak1 = 2;
596 Tauolapp::jaki_.jak2 = modeH;
597 Tauolapp::Tauola::setSameParticleDecayMode( 2 );
598 Tauolapp::Tauola::setOppositeParticleDecayMode( modeH );
604 Tauolapp::jaki_.jak1 = modeH;
605 Tauolapp::jaki_.jak2 = 2;
606 Tauolapp::Tauola::setSameParticleDecayMode( modeH );
607 Tauolapp::Tauola::setOppositeParticleDecayMode( 2 );
613 Tauolapp::jaki_.jak1 = modeH;
615 Tauolapp::Tauola::setSameParticleDecayMode( modeH );
616 Tauolapp::Tauola::setOppositeParticleDecayMode( Tauolapp::jaki_.
jak2 );
622 Tauolapp::jaki_.jak1 = modeH;
623 Tauolapp::jaki_.jak2 = 0;
624 Tauolapp::Tauola::setSameParticleDecayMode( modeH );
625 Tauolapp::Tauola::setOppositeParticleDecayMode( 0 );
631 Tauolapp::jaki_.jak1 = 0;
632 Tauolapp::jaki_.jak2 = modeH;
633 Tauolapp::Tauola::setSameParticleDecayMode( 0 );
634 Tauolapp::Tauola::setOppositeParticleDecayMode( modeH );
643 Tauolapp::Tauola::setSameParticleDecayMode( 0 );
644 Tauolapp::Tauola::setOppositeParticleDecayMode( 0 );
676 if ( prob > 0. && prob <= sumBra )
683 for (
int i=1;
i<NN;
i++ )
698 HepMC::GenEvent *
event =
new HepMC::GenEvent();
700 HepMC::FourVector momentum_tau1(l.Px(),l.Py(),l.Pz(),l.E());
704 HepMC::GenVertex *vertex =
new HepMC::GenVertex();
705 vertex->add_particle_out(tau1);
706 event->add_vertex(vertex);
711 partHep->set_status(p->status());
713 if(!partHep->end_vertex()){
714 HepMC::GenVertex* vtx =
new HepMC::GenVertex(p->end_vertex()->position());
715 theEvent->add_vertex(vtx);
716 vtx->add_particle_in(partHep);
718 if(p->end_vertex()->particles_out_size()!=0){
719 for(HepMC::GenVertex::particles_out_const_iterator d=p->end_vertex()->particles_out_const_begin(); d!=p->end_vertex()->particles_out_const_end();d++){
721 TLorentzVector
l((*d)->momentum().px(),(*d)->momentum().py(),(*d)->momentum().pz(),(*d)->momentum().e());
723 HepMC::FourVector momentum(
l.Px(),
l.Py(),
l.Pz(),
l.E());
725 daughter->suggest_barcode(theEvent->particles_size()+1);
726 partHep->end_vertex()->add_particle_out(daughter);
734 if ( tau->end_vertex() ) {
735 HepMC::GenVertex::particle_iterator dau;
736 for (dau = tau->end_vertex()->particles_begin(HepMC::children); dau!= tau->end_vertex()->particles_end(HepMC::children); dau++ ) {
737 int dau_pid = (*dau)->pdg_id();
738 if(dau_pid == tau->pdg_id())
return false;
745 std::vector<int> &m_idx){
748 TLorentzVector
t(tau->momentum().px(),tau->momentum().py(),tau->momentum().pz(),tau->momentum().e());
749 TLorentzVector
m(mother->momentum().px(),mother->momentum().py(),mother->momentum().pz(),mother->momentum().e());
750 for(
unsigned int i=0;
i<p.size();
i++){
751 if(tau->pdg_id()==p.at(
i).pdg_id()){
752 if(mother->pdg_id()==p.at(m_idx.at(
i)).pdg_id()){
753 TLorentzVector pm(p.at(m_idx.at(
i)).momentum().px(),p.at(m_idx.at(
i)).momentum().py(),p.at(m_idx.at(
i)).momentum().pz(),p.at(m_idx.at(
i)).momentum().e());
754 if(fabs(
m.M()-pm.M())<
dmMatch)
return spinup.at(
i);
762 if ( tau->production_vertex() ) {
763 HepMC::GenVertex::particle_iterator mother;
764 for (mother = tau->production_vertex()->particles_begin(
HepMC::parents); mother!= tau->production_vertex()->particles_end(
HepMC::parents); mother++ ) {
765 if((*mother)->pdg_id()==tau->pdg_id())
return FirstTauInChain(*mother);
772 if ( tau->production_vertex() ) {
773 HepMC::GenVertex::particle_iterator mother;
774 for (mother = tau->production_vertex()->particles_begin(
HepMC::parents); mother!= tau->production_vertex()->particles_end(
HepMC::parents); mother++ ) {
775 if((*mother)->pdg_id() == tau->pdg_id())
return GetMother(*mother);
783 if(p->end_vertex() && p->production_vertex()){
784 HepMC::GenVertex* PGenVtx=p->production_vertex();
785 HepMC::GenVertex* EGenVtx=p->end_vertex();
786 double VxDec = PGenVtx->position().x()+lab.Px()/prod.Px()*(EGenVtx->position().x()-PGenVtx->position().x());
787 double VyDec = PGenVtx->position().y()+lab.Py()/prod.Py()*(EGenVtx->position().y()-PGenVtx->position().y());
788 double VzDec = PGenVtx->position().z()+lab.Pz()/prod.Pz()*(EGenVtx->position().z()-PGenVtx->position().z());
789 double VtDec = PGenVtx->position().t()+lab.Pt()/prod.Pt()*(EGenVtx->position().t()-PGenVtx->position().t());
790 EGenVtx->set_position(HepMC::FourVector(VxDec,VyDec,VzDec,VtDec));
791 for(HepMC::GenVertex::particle_iterator dau=p->end_vertex()->particles_begin(HepMC::children); dau!=p->end_vertex()->particles_end(HepMC::children); dau++){
T getParameter(std::string const &) const
static AlgebraicMatrix initialize()
TauolappInterface(const edm::ParameterSet &)
T getUntrackedParameter(std::string const &, T const &) const
edm::ParameterSet * fPSet
std::vector< double > fScaledHadronBrRatios
bool exists(std::string const ¶meterName) const
checks if a parameter exists
HepMC::GenEvent * make_simple_tau_event(const TLorentzVector &l, int pdgid, int status)
void rmarin_(int *, int *, int *)
const HEPEUP * getHEPEUP() const
void getData(T &iHolder) const
std::vector< std::pair< int, int > > MOTHUP
void BoostProdToLabLifeTimeInDecays(HepMC::GenParticle *p, TLorentzVector &lab, TLorentzVector &prod)
std::vector< int > fLeptonModes
std::vector< FiveVector > PUP
std::vector< double > SPINUP
Abs< T >::type abs(const T &t)
void ranmar_(float *rvec, int *lenv)
How EventSelector::AcceptEvent() decides whether to accept an event for output otherwise it is excluding the probing of A single or multiple positive and the trigger will pass if any such matching triggers are PASS or EXCEPTION[A criterion thatmatches no triggers at all is detected and causes a throw.] A single negative with an expectation of appropriate bit checking in the decision and the trigger will pass if any such matching triggers are FAIL or EXCEPTION A wildcarded negative criterion that matches more than one trigger in the trigger but the state exists so we define the behavior If all triggers are the negative crieriion will lead to accepting the event(this again matches the behavior of"!*"before the partial wildcard feature was incorporated).The per-event"cost"of each negative criterion with multiple relevant triggers is about the same as!*was in the past
HepPDT::ParticleData ParticleData
void update_particles(HepMC::GenParticle *partHep, HepMC::GenEvent *theEvent, HepMC::GenParticle *p, TVector3 &boost)
void init(const edm::EventSetup &)
HepMC::GenParticle * FirstTauInChain(HepMC::GenParticle *tau)
HepMC::GenEvent * decay(HepMC::GenEvent *)
bool isLastTauInChain(const HepMC::GenParticle *tau)
std::vector< double > fScaledLeptonBrRatios
static CLHEP::HepRandomEngine * fRandomEngine
std::vector< int > fHadronModes
HepMC::GenParticle * GetMother(HepMC::GenParticle *tau)
std::pair< typename Association::data_type::first_type, double > match(Reference key, Association association, bool bestMatchByMaxValue)
Generic matching function.
edm::ESHandle< HepPDT::ParticleDataTable > fPDGTable
volatile std::atomic< bool > shutdown_flag false
double MatchedLHESpinUp(HepMC::GenParticle *tau, std::vector< HepMC::GenParticle > &p, std::vector< double > &spinup, std::vector< int > &m_idx)
void selectDecayByMDTAU()