diff --git a/[refs] b/[refs] index 01eff0daa2ef1..e2958facccd0c 100644 --- a/[refs] +++ b/[refs] @@ -1,2 +1,2 @@ --- -refs/heads/l1tmuon-upgrade-dev: 5391ca82dce0423f6ed012af12d45020d1d71263 +refs/heads/l1tmuon-upgrade-dev: 7f3a270d87d504697e3461ecb1701d7e6fdea7f2 diff --git a/trunk/GeneratorInterface/ExternalDecays/interface/TauolaInterface.h b/trunk/GeneratorInterface/ExternalDecays/interface/TauolaInterface.h index dcc97078df400..870b0d4bd730a 100644 --- a/trunk/GeneratorInterface/ExternalDecays/interface/TauolaInterface.h +++ b/trunk/GeneratorInterface/ExternalDecays/interface/TauolaInterface.h @@ -1,7 +1,14 @@ #ifndef gen_ExternalDecays_TauolaInterface_h #define gen_ExternalDecays_TauolaInterface_h +// #include "HepPDT/defs.h" +// #include "HepPDT/TableBuilder.hh" +#include "HepPDT/ParticleDataTable.hh" + #include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/Framework/interface/EventSetup.h" + namespace HepMC { @@ -19,7 +26,7 @@ namespace gen { TauolaInterface( const edm::ParameterSet& ); ~TauolaInterface(); - void init(); + void init( const edm::EventSetup& ); const std::vector& operatesOnParticles() { return fPDGs; } HepMC::GenEvent* decay( const HepMC::GenEvent* ); void statistics() ; @@ -29,6 +36,10 @@ namespace gen { // std::vector fPDGs; int fPolarization; + + edm::ESHandle fPDGTable ; + //CLHEP::HepRandomEngine* fRandomEngine ; + //CLHEP::RandFlat* fRandomGenerator; }; diff --git a/trunk/GeneratorInterface/ExternalDecays/interface/TauolaWrapper.h b/trunk/GeneratorInterface/ExternalDecays/interface/TauolaWrapper.h index 2b9ed76a4886f..af31ab8ef8fb0 100644 --- a/trunk/GeneratorInterface/ExternalDecays/interface/TauolaWrapper.h +++ b/trunk/GeneratorInterface/ExternalDecays/interface/TauolaWrapper.h @@ -26,10 +26,13 @@ //------------------------------------------------------------------------------- // -// main function of TAUOLA tau decay library +// main function(s) of TAUOLA/pretauola tau decay library extern "C" { void tauola_(int*, int*); + void tauola_srs_(int*,int*); + void taurep_(int*); + void ranmar_(float*,int*); } #define tauola tauola_ diff --git a/trunk/GeneratorInterface/ExternalDecays/src/TauolaInterface.cc b/trunk/GeneratorInterface/ExternalDecays/src/TauolaInterface.cc index 82d69697b1b55..64b6ff5607a18 100644 --- a/trunk/GeneratorInterface/ExternalDecays/src/TauolaInterface.cc +++ b/trunk/GeneratorInterface/ExternalDecays/src/TauolaInterface.cc @@ -8,6 +8,7 @@ #include "HepMC/IO_HEPEVT.h" #include "HepMC/HEPEVT_Wrapper.h" +#include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h" using namespace gen; using namespace edm; @@ -37,14 +38,20 @@ TauolaInterface::~TauolaInterface() { } -void TauolaInterface::init() +void TauolaInterface::init( const edm::EventSetup& es ) { if ( ki_taumod_.mdtau <= -1 ) // actually, need to throw exception ! return ; fPDGs.push_back( 15 ) ; - + es.getData( fPDGTable ) ; +/* + Service rng; + long seed = (long)(rng->mySeed()) ; + fRandomEngine = new HepJamesRandom(seed) ; + fRandomGenerator = new RandFlat(fRandomEngine) ; +*/ cout << "----------------------------------------------" << endl; cout << "Initializing Tauola" << endl; if ( fPolarization == 0 ) @@ -55,8 +62,11 @@ void TauolaInterface::init() { cout << "Tauola: Polarization enabled" << endl; } - int mode = -1; - tauola_( &mode, &fPolarization ); + int mode = -2; + taurep_( &mode ) ; + mode = -1; + // tauola_( &mode, &fPolarization ); + tauola_srs_( &mode, &fPolarization ); return; } @@ -82,12 +92,105 @@ HepMC::GenEvent* TauolaInterface::decay( const HepMC::GenEvent* evt ) // HepMC::HEPEVT_Wrapper::print_hepevt(); int mode = 0; - tauola_( &mode, &fPolarization ); + // tauola_( &mode, &fPolarization ); + tauola_srs_( &mode, &fPolarization ); int numPartAfterTauola = HepMC::HEPEVT_Wrapper::number_entries(); // HepMC::HEPEVT_Wrapper::print_hepevt(); HepMC::IO_HEPEVT conv; + + // before we do the conversion, we need to deal with decay vertexes + // since Tauola knows nothing about lifetimes, all decay vertexes are set to 0. + // nees to set them properly, knowing lifetime ! + // here we do it on HEPEVT record, also for consistency, although it's probably + // even easier to deal with HepMC::GenEvent record + + // find 1st "non-doc" tau + // + bool foundTau = false; + for ( int ip=1; ip<=numPartAfterTauola; ip++ ) + { + if ( std::abs( HepMC::HEPEVT_Wrapper::id( ip ) ) == 15 + && HepMC::HEPEVT_Wrapper::status( ip ) != 3 ) + { + foundTau = true; + break; + } + } + + if ( !foundTau ) + { + // no tau found + // just give up here + // + return conv.read_next_event(); + } + + std::vector PrntIndx; + + for ( int ip=numPartAfterTauola; ip>numPartBeforeTauola; ip-- ) // Fortran indexing ! + { + + // first of all, find out how many generations in decay chain + // + PrntIndx.clear(); + int Prnt = HepMC::HEPEVT_Wrapper::first_parent(ip); + ip -= (HepMC::HEPEVT_Wrapper::number_children(Prnt)-1); // such that we don't go the same part again + PrntIndx.push_back( Prnt ); + while ( abs( HepMC::HEPEVT_Wrapper::id(Prnt) ) != 15 ) // shortcut; need to loop over fPDGs... + { + int Prnt1 = HepMC::HEPEVT_Wrapper::first_parent( Prnt ); + Prnt = Prnt1; + // such that the tau always appear at the start of the list + PrntIndx.insert( PrntIndx.begin(), Prnt ); + ip -= HepMC::HEPEVT_Wrapper::number_children(Prnt); // such that we don't go the same part again + } + for ( int iprt=0; iprtparticle(HepPDT::ParticleID(abs(PartID))) ; + // + // prob = exp(-t/lifetime) ==> t = -lifetime * log(prob) + // + float prob = 0.; + int length=1; + ranmar_(&prob,&length); + double lifetime = PData->lifetime().value(); + // + // in case of Py6, this would be copied into V(5,i) + // for HEPEVT, need to check... + // + double ct = -lifetime * std::log(prob); + // + double ee = HepMC::HEPEVT_Wrapper::e( Indx ); + double px = HepMC::HEPEVT_Wrapper::px( Indx ); + double py = HepMC::HEPEVT_Wrapper::py( Indx ); + double pz = HepMC::HEPEVT_Wrapper::pz( Indx ); + double pp = std::sqrt( px*px + py*py + pz*pz ); + double mass = HepMC::HEPEVT_Wrapper::m( Indx ); + // + // this is in py6 terms: + // VDCY(J)=V(IP,J)+V(IP,5)*P(IP,J)/P(IP,5) + // + double VxDec = HepMC::HEPEVT_Wrapper::x( Indx ); + VxDec += ct * (px/mass); + double VyDec = HepMC::HEPEVT_Wrapper::y( Indx ); + VyDec += ct * (py/mass); + double VzDec = HepMC::HEPEVT_Wrapper::z( Indx ); + VzDec += ct * (pz/mass); + double VtDec = HepMC::HEPEVT_Wrapper::t( Indx ); + VtDec += ct * (ee/mass); + for ( int idau=HepMC::HEPEVT_Wrapper::first_child( Indx ); + idau<=HepMC::HEPEVT_Wrapper::last_child( Indx ); idau++ ) + { + HepMC::HEPEVT_Wrapper::set_position( idau, VxDec, VyDec, VzDec, VtDec ); + } + } + } + return conv.read_next_event(); } @@ -95,5 +198,6 @@ HepMC::GenEvent* TauolaInterface::decay( const HepMC::GenEvent* evt ) void TauolaInterface::statistics() { int mode = 1; - tauola_( &mode, &fPolarization ) ; + // tauola_( &mode, &fPolarization ) ; + tauola_srs_( &mode, &fPolarization ) ; }