CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/GeneratorInterface/ExternalDecays/src/TauolaInterface.cc

Go to the documentation of this file.
00001 /* for old tauola27 */
00002 #include <iostream>
00003 
00004 #include "GeneratorInterface/Pythia6Interface/interface/Pythia6Service.h"
00005 
00006 #include "GeneratorInterface/ExternalDecays/interface/TauolaInterface.h"
00007 #include "GeneratorInterface/ExternalDecays/interface/TauolaWrapper.h"
00008 
00009 #include "FWCore/ServiceRegistry/interface/Service.h"
00010 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
00011 
00012 #include "HepMC/GenEvent.h"
00013 #include "HepMC/IO_HEPEVT.h"
00014 #include "HepMC/HEPEVT_Wrapper.h"
00015 
00016 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
00017 
00018 #include "GeneratorInterface/ExternalDecays/interface/DecayRandomEngine.h"
00019 
00020 using namespace gen;
00021 using namespace edm;
00022 using namespace std;
00023 
00024 extern "C" {
00025 
00026   void ranmar_( float *rvec, int *lenv )
00027   {
00028 
00029       for(int i = 0; i < *lenv; i++)
00030          *rvec++ = decayRandomEngine->flat();
00031 
00032       return;
00033 
00034   }
00035   
00036   void rmarin_( int*, int*, int* )
00037   {
00038 
00039      return;
00040 
00041   }
00042 
00043 }
00044 
00045 //
00046 //   General Note: While there're no explicit calls or otherwise "links" to Pythia6 anywhere,
00047 //   we're using Pythia6Service here because we run pretauola rather than "core" tauola;
00048 //   pretauola is an extension on top of tauola, which is tied to Pythia6 via several routines;
00049 //   most heavily use one is PYR - we can't avoid it (other Pythia6-tied routines we avoid)
00050 //
00051 
00052 TauolaInterface::TauolaInterface( const ParameterSet& pset )
00053    : fIsInitialized(false)
00054 {
00055    fPy6Service = new Pythia6Service;
00056 
00057    fPolarization = pset.getParameter<bool>("UseTauolaPolarization") ? 1 : 0 ;
00058 
00059    // set Tauola defaults
00060    //
00061    ki_taumod_.pjak1 = -1;
00062    ki_taumod_.pjak2 = -1;
00063    ki_taumod_.mdtau = -1;
00064    
00065    // read tau decay mode switches
00066    //
00067    ParameterSet cards = pset.getParameter< ParameterSet >("InputCards");
00068    ki_taumod_.pjak1 = cards.getParameter< int >( "pjak1" );
00069    ki_taumod_.pjak2 = cards.getParameter< int >( "pjak2" );
00070    ki_taumod_.mdtau = cards.getParameter< int >( "mdtau" );
00071    
00072 } 
00073 
00074 TauolaInterface::~TauolaInterface()
00075 {
00076    delete fPy6Service;
00077 }
00078 
00079 void TauolaInterface::init( const edm::EventSetup& es )
00080 {
00081    
00082    if ( fIsInitialized ) return; // do init only once
00083    
00084    if ( ki_taumod_.mdtau <= -1 ) // actually, need to throw exception !
00085       return ;
00086    
00087    fPDGs.push_back( 15 ) ;
00088    es.getData( fPDGTable ) ;
00089 
00090         cout << "----------------------------------------------" << endl;
00091         cout << "Initializing Tauola" << endl;
00092         if ( fPolarization == 0 )
00093         {
00094            cout << "Tauola: Polarization disabled" << endl;
00095         } 
00096         else if ( fPolarization == 1 )
00097         {
00098            cout << "Tauola: Polarization enabled" << endl;
00099         }
00100 
00101 // FIXME !!!
00102 // This is a temporary hack - we're re-using master generator's seed to init RANMAR
00103 // FIXME !!!
00104 //   This is now off because ranmar has been overriden (see code above) to use pyr_(...)
00105 //   - this way we're using guaranteed initialized rndm generator... BUT !!! in the long
00106 //   run we may want a separate random stream for tauola...
00107 
00108 //   Service<RandomNumberGenerator> rng;
00109 //   int seed = rng->mySeed() ;
00110 //   int ntot=0, ntot2=0;
00111 //   rmarin_( &seed, &ntot, &ntot2 );
00112 
00113    int mode = -2;
00114    taurep_( &mode ) ;
00115    mode = -1;
00116    // tauola_( &mode, &fPolarization );
00117    // tauola_srs_( &mode, &fPolarization );
00118    //
00119    // We're using the call(...) method here because it'll make sure that Py6 
00120    // is initialized, and that's done only once, and will grab exatly that instance
00121    //
00122    fPy6Service->call( tauola_srs_, &mode, &fPolarization ); 
00123    
00124    fIsInitialized = true;
00125    
00126    return;
00127 }
00128 
00129 HepMC::GenEvent* TauolaInterface::decay( HepMC::GenEvent* evt )
00130 {
00131    
00132    // event record convertor
00133    //
00134    HepMC::IO_HEPEVT conv;
00135    
00136    if ( !fIsInitialized ) return conv.read_next_event();
00137    
00138    // We are using random numbers, we are fetched through Pythia6Service
00139    // (through ranmar_ below) -> so grab the instance during decay()
00140 
00141    Pythia6Service::InstanceWrapper pythia6InstanceGuard( fPy6Service );
00142 
00143    // fill up HEPEVT common block
00144    //
00145    // IDEALLY, this should be the way to go
00146    // BUT !!! this utility fills it up in the "reshuffled" order,
00147    // and later on Tauola chocks on it 
00148    //
00149    // Needs to be sorted out, eith in HepMC, or in Tauola, or both !!!
00150    // 
00151    // At present, this thing blindly relies on the assumption that
00152    // HEPEVT is always there - which wont be the case with Py8 or Hwg++
00153    //
00154    //HepMC::IO_HEPEVT conv;
00155    //conv.write_event( evt ) ;
00156    
00157    int numPartBeforeTauola = HepMC::HEPEVT_Wrapper::number_entries();
00158    // HepMC::HEPEVT_Wrapper::print_hepevt();
00159    
00160    int mode = 0;
00161    // tauola_( &mode, &fPolarization );
00162    fPy6Service->call( tauola_srs_, &mode, &fPolarization );
00163    
00164    int numPartAfterTauola = HepMC::HEPEVT_Wrapper::number_entries();
00165    // HepMC::HEPEVT_Wrapper::print_hepevt();
00166    
00167    // before we do the conversion, we need to deal with decay vertexes
00168    // since Tauola knows nothing about lifetimes, all decay vertexes are set to 0. 
00169    // nees to set them properly, knowing lifetime !
00170    // here we do it on HEPEVT record, also for consistency, although it's probably
00171    // even easier to deal with HepMC::GenEvent record  
00172    
00173    // find 1st "non-doc" tau
00174    //
00175    bool foundTau = false;
00176    for ( int ip=1; ip<=numPartAfterTauola; ip++ )
00177    {
00178       if ( std::abs( HepMC::HEPEVT_Wrapper::id( ip ) ) == 15
00179            && HepMC::HEPEVT_Wrapper::status( ip ) != 3 )
00180       {
00181          foundTau = true;
00182          break;
00183       }
00184    }
00185    
00186    if ( !foundTau )
00187    {
00188       // no tau found
00189       // just give up here
00190       //
00191       return conv.read_next_event();
00192    }
00193    
00194    std::vector<int> PrntIndx;
00195    
00196    for ( int ip=numPartAfterTauola; ip>numPartBeforeTauola; ip-- ) // Fortran indexing !
00197    {
00198       
00199       // first of all, find out how many generations in decay chain
00200       //
00201       PrntIndx.clear();
00202       int Prnt = HepMC::HEPEVT_Wrapper::first_parent(ip);
00203       ip -= (HepMC::HEPEVT_Wrapper::number_children(Prnt)-1); // such that we don't go the same part again
00204       PrntIndx.push_back( Prnt );
00205       while ( abs( HepMC::HEPEVT_Wrapper::id(Prnt) ) != 15 ) // shortcut; need to loop over fPDGs...
00206       {
00207          int Prnt1 = HepMC::HEPEVT_Wrapper::first_parent( Prnt );
00208          Prnt = Prnt1;
00209          // such that the tau always appear at the start of the list
00210          PrntIndx.insert( PrntIndx.begin(), Prnt );
00211          ip -= HepMC::HEPEVT_Wrapper::number_children(Prnt); // such that we don't go the same part again
00212       }
00213       for ( size_t iprt=0; iprt<PrntIndx.size(); iprt++ )
00214       {  
00215           int Indx = PrntIndx[iprt];
00216           int PartID = HepMC::HEPEVT_Wrapper::id( Indx );
00217           const HepPDT::ParticleData* 
00218              PData = fPDGTable->particle(HepPDT::ParticleID(abs(PartID))) ;
00219          //
00220          // prob = exp(-t/lifetime) ==> t = -lifetime * log(prob)
00221          //
00222          float prob = 0.;
00223          int length=1;
00224          ranmar_(&prob,&length);
00225          double lifetime = PData->lifetime().value();
00226          //
00227          // in case of Py6, this would be copied into V(5,i)
00228          // for HEPEVT, need to check...
00229          //
00230          double ct = -lifetime * std::log(prob);
00231          //
00232          double ee = HepMC::HEPEVT_Wrapper::e( Indx );
00233          double px = HepMC::HEPEVT_Wrapper::px( Indx );
00234          double py = HepMC::HEPEVT_Wrapper::py( Indx );
00235          double pz = HepMC::HEPEVT_Wrapper::pz( Indx );
00236          // double pp = std::sqrt( px*px + py*py + pz*pz );
00237          double mass = HepMC::HEPEVT_Wrapper::m( Indx );
00238          //
00239          // this is in py6 terms:
00240          //  VDCY(J)=V(IP,J)+V(IP,5)*P(IP,J)/P(IP,5)
00241          //
00242          double VxDec = HepMC::HEPEVT_Wrapper::x( Indx );
00243          VxDec += ct * (px/mass);
00244          double VyDec = HepMC::HEPEVT_Wrapper::y( Indx );
00245          VyDec += ct * (py/mass);
00246          double VzDec = HepMC::HEPEVT_Wrapper::z( Indx );
00247          VzDec += ct * (pz/mass);
00248          double VtDec = HepMC::HEPEVT_Wrapper::t( Indx );
00249          VtDec += ct * (ee/mass);
00250          for ( int idau=HepMC::HEPEVT_Wrapper::first_child( Indx );
00251                    idau<=HepMC::HEPEVT_Wrapper::last_child( Indx ); idau++ )
00252          {
00253             HepMC::HEPEVT_Wrapper::set_position( idau, VxDec, VyDec, VzDec, VtDec );
00254          }
00255       }
00256    }
00257    
00258    return conv.read_next_event();
00259       
00260 }
00261 
00262 void TauolaInterface::statistics()
00263 {
00264    int mode = 1;
00265    // tauola_( &mode, &fPolarization ) ;
00266    // tauola_srs_( &mode, &fPolarization ) ;
00267    fPy6Service->call( tauola_srs_, &mode, &fPolarization );
00268    return;
00269 }
00270 
00271 /* */
00272 
00273 /* this is the code for the new Tauola++ 
00274 
00275 #include <iostream>
00276 
00277 #include "GeneratorInterface/ExternalDecays/interface/TauolaInterface.h"
00278 
00279 #include "Tauola.h"
00280 #include "TauolaHepMCEvent.h"
00281 #include "Log.h"
00282 
00283 #include "FWCore/ServiceRegistry/interface/Service.h"
00284 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
00285 #include "FWCore/Utilities/interface/Exception.h"
00286 
00287 #include "CLHEP/Random/RandomEngine.h"
00288 
00289 #include "HepMC/GenEvent.h"
00290 #include "HepMC/IO_HEPEVT.h"
00291 #include "HepMC/HEPEVT_Wrapper.h"
00292 
00293 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
00294 
00295 // #include "GeneratorInterface/ExternalDecays/interface/DecayRandomEngine.h"
00296 
00297 
00298 extern "C" {
00299 
00300   void gen::ranmar_( float *rvec, int *lenv )
00301   {
00302       TauolaInterface* instance = TauolaInterface::getInstance();
00303       for(int i = 0; i < *lenv; i++)
00304          // *rvec++ = decayRandomEngine->flat();
00305          *rvec++ = instance->flat();
00306       return;
00307   }
00308   
00309   void gen::rmarin_( int*, int*, int* )
00310   {
00311      return;
00312   }
00313 
00314 }
00315 
00316 using namespace gen;
00317 using namespace edm;
00318 using namespace std;
00319 
00320 TauolaInterface* TauolaInterface::fInstance = 0;
00321 
00322 
00323 TauolaInterface::TauolaInterface()
00324    : fPolarization(false), fPSet(0), fIsInitialized(false)
00325 {
00326    
00327    Service<RandomNumberGenerator> rng;
00328    if(!rng.isAvailable()) {
00329     throw cms::Exception("Configuration")
00330        << "The RandomNumberProducer module requires the RandomNumberGeneratorService\n"
00331           "which appears to be absent.  Please add that service to your configuration\n"
00332           "or remove the modules that require it." << std::endl;
00333    }
00334    
00335    fRandomEngine = &rng->getEngine();
00336 
00337 }
00338 
00339 
00340 //TauolaInterface::TauolaInterface( const ParameterSet& pset )
00341 //   : fIsInitialized(false)
00342 //{
00343 //
00344 //   Tauola::setDecayingParticle(15);
00345 //   // --> ??? Tauola::setRadiation(false);
00346 //
00347 //   // polarization switch 
00348 //   //
00349 //   // fPolarization = pset.getParameter<bool>("UseTauolaPolarization") ? 1 : 0 ;
00350 //   fPolarization = pset.getParameter<bool>("UseTauolaPolarization");
00351 //   
00352 //   // read tau decay mode switches
00353 //   //
00354 //   ParameterSet cards = pset.getParameter< ParameterSet >("InputCards");
00355 //   Tauola::setSameParticleDecayMode( cards.getParameter< int >( "pjak1" ) ) ;
00356 //   Tauola::setOppositeParticleDecayMode( cards.getParameter< int >( "pjak2" ) ) ;
00357 //
00358 //   Tauola::setTauLifetime(0.0);
00359 //   Tauola::spin_correlation.setAll(fPolarization);
00360 //
00361 //   // some more options, copied over from an example 
00362 //   // - maybe will use later...
00363 //   //
00364 //   //Tauola::setEtaK0sPi(0,0,0); // switches to decay eta K0_S and pi0 1/0 on/off. 
00365 //   //
00366 //
00367 //} 
00368 
00369 
00370 TauolaInterface* TauolaInterface::getInstance()
00371 {
00372 
00373    if ( fInstance == 0 ) fInstance = new TauolaInterface() ;
00374    return fInstance;
00375 
00376 }
00377 
00378 
00379 TauolaInterface::~TauolaInterface()
00380 {
00381 
00382    if ( fPSet != 0 ) delete fPSet;
00383    if ( fInstance == this ) fInstance = 0;
00384 
00385 }
00386 
00387 void TauolaInterface::setPSet( const ParameterSet& pset )
00388 {
00389 
00390    if ( fPSet != 0 ) 
00391    {
00392       throw cms::Exception("TauolaInterfaceError")
00393          << "Attempt to override Tauola an existing ParameterSet\n"
00394          << std::endl;   
00395    }
00396    
00397    fPSet = new ParameterSet(pset);
00398    
00399    return;
00400 
00401 }
00402 
00403 void TauolaInterface::init( const edm::EventSetup& es )
00404 {
00405    
00406    if ( fIsInitialized ) return; // do init only once
00407    
00408    if ( fPSet == 0 ) 
00409    {
00410 
00411       throw cms::Exception("TauolaInterfaceError")
00412          << "Attempt to initialize Tauola with an empty ParameterSet\n"
00413          << std::endl;   
00414    }
00415    
00416    fIsInitialized = true;
00417       
00418    es.getData( fPDGTable ) ;
00419 
00420    Tauola::setDecayingParticle(15);
00421    // --> ??? Tauola::setRadiation(false);
00422 
00423    // polarization switch 
00424    //
00425    // fPolarization = fPSet->getParameter<bool>("UseTauolaPolarization") ? 1 : 0 ;
00426    fPolarization = fPSet->getParameter<bool>("UseTauolaPolarization");
00427    
00428    // read tau decay mode switches
00429    //
00430    ParameterSet cards = fPSet->getParameter< ParameterSet >("InputCards");
00431    Tauola::setSameParticleDecayMode( cards.getParameter< int >( "pjak1" ) ) ;
00432    Tauola::setOppositeParticleDecayMode( cards.getParameter< int >( "pjak2" ) ) ;
00433 
00434    Tauola::setTauLifetime(0.0);
00435    Tauola::spin_correlation.setAll(fPolarization);
00436 
00437    // some more options, copied over from an example 
00438    // - maybe will use later...
00439    //
00440    //Tauola::setEtaK0sPi(0,0,0); // switches to decay eta K0_S and pi0 1/0 on/off. 
00441    //
00442 
00443 //
00444 //   const HepPDT::ParticleData* 
00445 //         PData = fPDGTable->particle(HepPDT::ParticleID( abs(Tauola::getDecayingParticle()) )) ;
00446 //   double lifetime = PData->lifetime().value();
00447 //   Tauola::setTauLifetime( lifetime );
00448 
00449    fPDGs.push_back( Tauola::getDecayingParticle() );
00450          
00451    Tauola::initialise();
00452    Log::LogWarning(false);
00453    
00454    return;
00455 }
00456 
00457 float TauolaInterface::flat()
00458 {
00459 
00460    if ( !fPSet )
00461    {
00462       // throw
00463       throw cms::Exception("TauolaInterfaceError")
00464          << "Attempt to run random number generator of un-initialized Tauola\n"
00465          << std::endl;   
00466    }
00467    
00468    if ( !fIsInitialized ) 
00469    {
00470       // throw
00471       throw cms::Exception("TauolaInterfaceError")
00472          << "Attempt to run random number generator of un-initialized Tauola\n"
00473          << std::endl;   
00474    }
00475    
00476    return fRandomEngine->flat();
00477 
00478 }
00479 
00480 HepMC::GenEvent* TauolaInterface::decay( HepMC::GenEvent* evt )
00481 {
00482       
00483    if ( !fIsInitialized ) return evt;
00484    
00485    int NPartBefore = evt->particles_size();
00486    int NVtxBefore  = evt->vertices_size();
00487    
00488    // what do we do if Hep::GenEvent size is larger than 10K ???
00489    // Tauola (& Photos, BTW) can only handle up to 10K via HEPEVT,
00490    // and in case of CMS, it's only up to 4K !!!
00491    //
00492    // if ( NPartBefore > 10000 ) return evt;
00493    //
00494    
00495     //construct tmp TAUOLA event
00496     //
00497     TauolaHepMCEvent * t_event = new TauolaHepMCEvent(evt);
00498    
00499     // another option: if one lets Pythia or another master gen to decay taus, 
00500     //                 we have to undecay them first
00501     // t_event->undecayTaus();
00502     
00503     // run Tauola on the tmp event - HepMC::GenEvernt will be MODIFIED !!!
00504     //
00505     t_event->decayTaus();
00506     
00507     // delet tmp Tauola event
00508     //
00509     delete t_event; 
00510     
00511     // fix barcodes of the newly added particles
00512     //
00513     for(HepMC::GenEvent::particle_const_iterator it = evt->particles_begin(); 
00514         it != evt->particles_end(); ++it ) 
00515     {
00516        //HepMC::GenParticle* GenPrt = (*it);
00517        if ( (*it)->barcode() > 10000 )
00518        {
00519           int NewBarcode = ((*it)->barcode()-10000) + NPartBefore;
00520           (*it)->suggest_barcode( NewBarcode ); 
00521        }
00522     }
00523     
00524 
00525     // do we also need to apply the lifetime and vtx position shift ??? 
00526     // (see TauolaInterface, for example)
00527     //
00528 
00529     for ( int iv=NVtxBefore+1; iv<=evt->vertices_size(); iv++ )
00530     {
00531        HepMC::GenVertex* GenVtx = evt->barcode_to_vertex(-iv);
00532        HepMC::GenParticle* GenPart = *(GenVtx->particles_in_const_begin());
00533        HepMC::GenVertex* ProdVtx = GenPart->production_vertex();
00534        HepMC::FourVector PMom = GenPart->momentum();
00535        double mass = GenPart->generated_mass();
00536        const HepPDT::ParticleData* 
00537              PData = fPDGTable->particle(HepPDT::ParticleID(abs(GenPart->pdg_id()))) ;
00538        double lifetime = PData->lifetime().value();
00539        float prob = 0.;
00540        int length=1;
00541        ranmar_(&prob,&length);
00542        double ct = -lifetime * std::log(prob);
00543        double VxDec = GenVtx->position().x();
00544        VxDec += ct * (PMom.px()/mass);
00545        VxDec += ProdVtx->position().x();
00546        double VyDec = GenVtx->position().y();
00547        VyDec += ct * (PMom.py()/mass);
00548        VyDec += ProdVtx->position().y();
00549        double VzDec = GenVtx->position().z();
00550        VzDec += ct * (PMom.pz()/mass);
00551        VzDec += ProdVtx->position().z();
00552        double VtDec = GenVtx->position().t();
00553        VtDec += ct * (PMom.e()/mass);
00554        VtDec += ProdVtx->position().t();
00555        GenVtx->set_position( HepMC::FourVector(VxDec,VyDec,VzDec,VtDec) );       
00556     }
00557     
00558     return evt;
00559       
00560 }
00561 
00562 void TauolaInterface::statistics()
00563 {
00564    return;
00565 }
00566 
00567 */