#include <TauolappInterface.h>

Inheritance diagram for gen::TauolappInterface:

Public Member Functions
HepMC::GenEvent *	decay (HepMC::GenEvent *)

void	disablePolarization ()

void	enablePolarization ()

void	init (const edm::EventSetup &)

const std::vector< int > &	operatesOnParticles ()

virtual void	SetLHE (lhef::LHEEvent *l)

void	setRandomEngine (CLHEP::HepRandomEngine *v)

void	statistics ()

	TauolappInterface (const edm::ParameterSet &)

	~TauolappInterface ()

Public Member Functions inherited from gen::TauolaInterfaceBase
virtual void	SetDecayRandomEngine (CLHEP::HepRandomEngine *decayRandomEngine)

	TauolaInterfaceBase ()

	TauolaInterfaceBase (const edm::ParameterSet &)

virtual	~TauolaInterfaceBase ()

Static Public Member Functions
static double	flat ()

Private Member Functions
void	BoostProdToLabLifeTimeInDecays (HepMC::GenParticle *p, TLorentzVector &lab, TLorentzVector &prod)

void	decodeMDTAU (int)

HepMC::GenParticle *	FirstTauInChain (HepMC::GenParticle *tau)

HepMC::GenParticle *	GetMother (HepMC::GenParticle *tau)

bool	isLastTauInChain (const HepMC::GenParticle *tau)

HepMC::GenEvent *	make_simple_tau_event (const TLorentzVector &l, int pdgid, int status)

double	MatchedLHESpinUp (HepMC::GenParticle *tau, std::vector< HepMC::GenParticle > &p, std::vector< double > &spinup, std::vector< int > &m_idx)

void	selectDecayByMDTAU ()

int	selectHadronic ()

int	selectLeptonic ()

void	update_particles (HepMC::GenParticle partHep, HepMC::GenEvent theEvent, HepMC::GenParticle *p, TVector3 &boost)

Private Attributes
double	dmMatch

bool	dolhe

bool	dolheBosonCorr

std::vector< int >	fHadronModes

bool	fIsInitialized

std::vector< int >	fLeptonModes

int	fMDTAU

std::vector< int >	fPDGs

edm::ESHandle < HepPDT::ParticleDataTable >	fPDGTable

bool	fPolarization

edm::ParameterSet *	fPSet

std::vector< double >	fScaledHadronBrRatios

std::vector< double >	fScaledLeptonBrRatios

bool	fSelectDecayByEvent

lhef::LHEEvent *	lhe

double	lifetime

int	ntries

Static Private Attributes
static CLHEP::HepRandomEngine *	fRandomEngine = 0

Additional Inherited Members
Protected Attributes inherited from gen::TauolaInterfaceBase
std::vector< int >	fPDGs

Detailed Description

Definition at line 30 of file TauolappInterface.h.

Constructor & Destructor Documentation

TauolappInterface::TauolappInterface ( const edm::ParameterSet & pset )

Definition at line 55 of file TauolappInterface.cc.

References edm::hlt::Exception, and fPSet.

                                                                 :
   fPolarization(false),
   fPSet(0),
   fIsInitialized(false),
   fMDTAU(-1),
   fSelectDecayByEvent(false),
   lhe(NULL),
   dmMatch(0.5),
   dolhe(false),
   dolheBosonCorr(false),
   ntries(10)
 {
   if ( fPSet != 0 ) throw cms::Exception("TauolappInterfaceError") << "Attempt to override Tauola an existing ParameterSet\n" << std::endl;
   fPSet = new ParameterSet(pset);
 }

TauolappInterface::~TauolappInterface ( )

Definition at line 71 of file TauolappInterface.cc.

References fPSet.

71 {if ( fPSet != 0 ) delete fPSet;}

gen::TauolappInterface::fPSet

edm::ParameterSet * fPSet

Definition: TauolappInterface.h:67

Member Function Documentation

void TauolappInterface::BoostProdToLabLifeTimeInDecays	(	HepMC::GenParticle *	p,
		TLorentzVector &	lab,
		TLorentzVector &	prod
	)

private

Definition at line 789 of file TauolappInterface.cc.

Referenced by decay().

                                                                                                                    {
   if(p->end_vertex() && p->production_vertex()){ 
     HepMC::GenVertex* PGenVtx=p->production_vertex();
     HepMC::GenVertex* EGenVtx=p->end_vertex(); 
     double VxDec = PGenVtx->position().x()+lab.Px()/prod.Px()*(EGenVtx->position().x()-PGenVtx->position().x()); 
     double VyDec = PGenVtx->position().y()+lab.Py()/prod.Py()*(EGenVtx->position().y()-PGenVtx->position().y()); 
     double VzDec = PGenVtx->position().z()+lab.Pz()/prod.Pz()*(EGenVtx->position().z()-PGenVtx->position().z()); 
     double VtDec = PGenVtx->position().t()+lab.Pt()/prod.Pt()*(EGenVtx->position().t()-PGenVtx->position().t()); 
     EGenVtx->set_position(HepMC::FourVector(VxDec,VyDec,VzDec,VtDec)); 
     for(HepMC::GenVertex::particle_iterator dau=p->end_vertex()->particles_begin(HepMC::children); dau!=p->end_vertex()->particles_end(HepMC::children); dau++){ 
       BoostProdToLabLifeTimeInDecays((*dau),lab,prod); //recursively modify everything in the decay chain 
     } 
   } 
 }  

HepMC::GenEvent * TauolappInterface::decay ( HepMC::GenEvent * evt )

virtual

Reimplemented from gen::TauolaInterfaceBase.

Definition at line 201 of file TauolappInterface.cc.

References funct::abs(), BoostProdToLabLifeTimeInDecays(), dolhe, dolheBosonCorr, fIsInitialized, flat(), fSelectDecayByEvent, GenParticle::GenParticle, lhef::LHEEvent::getHEPEUP(), GetMother(), i, lhef::HEPEUP::IDUP, isLastTauInChain(), lhef::HEPEUP::ISTUP, getDQMSummary::iter, j, gen::k, lhe, visualization-live-secondInstance_cfg::m, make_simple_tau_event(), match(), MatchedLHESpinUp(), lhef::HEPEUP::MOTHUP, ntries, NULL, gen::p, p1, parents, lhef::HEPEUP::PUP, selectDecayByMDTAU(), lhef::HEPEUP::SPINUP, ntuplemaker::status, and update_particles().

                                                            {
    if ( !fIsInitialized ) return evt;
    Tauolapp::Tauola::setRandomGenerator(gen::TauolappInterface::flat);   // rest tauola++ random number incase other modules use tauola++
    int NPartBefore = evt->particles_size();
    int NVtxBefore  = evt->vertices_size();
    
    // what do we do if Hep::GenEvent size is larger than 10K ???
    // Tauola (& Photos, BTW) can only handle up to 10K via HEPEVT,
    // and in case of CMS, it's only up to 4K !!!
    // override decay mode if needs be
 
    if(fSelectDecayByEvent){
      selectDecayByMDTAU();
    }
    if(dolhe && lhe!=NULL){
      std::vector<HepMC::GenParticle> particles;
      std::vector<int> m_idx; 
      std::vector<double> spinup=lhe->getHEPEUP()->SPINUP;
      std::vector<int> pdg =lhe->getHEPEUP()->IDUP;
      for(unsigned int i=0;i<spinup.size();i++){
        particles.push_back(HepMC::GenParticle(HepMC::FourVector(lhe->getHEPEUP()->PUP.at(i)[0],lhe->getHEPEUP()->PUP.at(i)[1],lhe->getHEPEUP()->PUP.at(i)[2],lhe->getHEPEUP()->PUP.at(i)[3]),lhe->getHEPEUP()->IDUP.at(i)));
        int status = lhe->getHEPEUP()->ISTUP.at(i);
      particles.at(particles.size()-1).set_generated_mass(lhe->getHEPEUP()->PUP.at(i)[4]);
      particles.at(particles.size()-1).set_status(status > 0 ? (status == 2 ? 3 : status) : 3);
      m_idx.push_back(lhe->getHEPEUP()->MOTHUP.at(i).first-1);// correct for fortran index offset
      }
      // match to taus in hepmc and identify mother of taus  
      bool hastaus(false);
      std::vector<HepMC::GenParticle*> match;
      for(HepMC::GenEvent::particle_const_iterator iter = evt->particles_begin(); iter != evt->particles_end(); iter++) {
        if(abs((*iter)->pdg_id())==15){
      hastaus=true;
      int mother_pid(0);
      // check imediate parent to avoid parent tau ie tau->taugamma
      for(HepMC::GenVertex::particle_iterator mother=(*iter)->production_vertex()->particles_begin(HepMC::parents); mother!=(*iter)->production_vertex()->particles_end(HepMC::parents); mother++) {
        mother_pid = (*mother)->pdg_id();
        if(mother_pid!=(*iter)->pdg_id()){
          // match against lhe record
          if(abs(mother_pid) == 24 || // W
         abs(mother_pid) == 37 || // H+/-
         abs(mother_pid) == 23 || // Z
         abs(mother_pid) == 22 || // gamma
         abs(mother_pid) == 25 || // H0 SM
         abs(mother_pid) == 35 || // H0
         abs(mother_pid) == 36    // A0
         ){
            bool isfound=false;
            for(unsigned int k=0;k<match.size();k++){
          if((*mother)==match.at(k))isfound=true;
            }
            if(!isfound) match.push_back(*mother);
          }
        }
      }
        }
      }
      if(hastaus){
        // if is single gauge boson decay and match helicities
        if(match.size()==1 && dolheBosonCorr){
      for(int i=0; i<ntries;i++){
        // re-decay taus then check if helicities match 
        auto * t_event = new Tauolapp::TauolaHepMCEvent(evt);
        t_event->undecayTaus();
        t_event->decayTaus();
        bool ismatch=true;
        for(unsigned int j=0;j<spinup.size();j++){
          if(abs(pdg.at(j))==15){
            double diffhelminus=(-1.0*(double)Tauolapp::Tauola::getHelMinus()-spinup.at(j)); // -1.0 to correct for tauola feature 
            double diffhelplus=((double)Tauolapp::Tauola::getHelPlus()-spinup.at(j));
            if(pdg.at(j)==15  && diffhelminus>0.5) ismatch=false;
            if(pdg.at(j)==-15 && diffhelplus>0.5)  ismatch=false;
          }
        }
        delete t_event;
        if(ismatch) break;
      }
        }
        else{
      // If the event does not contain a single gauge boson the code will be run with 
      // remove all tau decays
      auto * t_event = new Tauolapp::TauolaHepMCEvent(evt);
      t_event->undecayTaus();
      delete t_event;
      // decay all taus manually based on the helicity
      for(HepMC::GenEvent::particle_const_iterator iter = evt->particles_begin(); iter != evt->particles_end(); iter++) {
        if(abs((*iter)->pdg_id())==15 && isLastTauInChain(*iter)){
          TLorentzVector ltau((*iter)->momentum().px(),(*iter)->momentum().py(),(*iter)->momentum().pz(),(*iter)->momentum().e());
          HepMC::GenParticle *m=GetMother(*iter);
          TLorentzVector mother(m->momentum().px(),m->momentum().py(),m->momentum().pz(),m->momentum().e());
          TVector3 boost=-1.0*mother.BoostVector(); // boost into mother's CM frame
          TLorentzVector ltau_lab=ltau;
          ltau.Boost(boost);
          mother.Boost(boost);
          HepMC::GenEvent *tauevt=make_simple_tau_event(ltau,(*iter)->pdg_id(),(*iter)->status());
          HepMC::GenParticle *p=(*(tauevt->particles_begin()));
          Tauolapp::TauolaParticle *tp = new Tauolapp::TauolaHepMCParticle(p);
          double helicity=MatchedLHESpinUp(*iter,particles,spinup,m_idx); // get helicity from lhe
          if((*iter)->pdg_id()==15) helicity*=-1.0;
          tp->undecay();
          // use |S_{tau}|=0.999999 to avoid issues with numerical roundoff 
          Tauolapp::Tauola::decayOne(tp,true,0,0,((double)helicity)*0.999999);
          boost*=-1.0; // boost back to lab frame
          mother.Boost(boost);
          update_particles((*iter),evt,p,boost);
          //correct tau liftetime for boost (change rest frame from mothers to taus)
          BoostProdToLabLifeTimeInDecays((*iter),ltau_lab,ltau);
          delete tauevt;
        }
      }
        }
      }
    }
    else{
      //construct tmp TAUOLA event
      auto * t_event = new Tauolapp::TauolaHepMCEvent(evt);
      //t_event->undecayTaus();
      t_event->decayTaus();
      delete t_event; 
    }
 
    for ( int iv=NVtxBefore+1; iv<=evt->vertices_size(); iv++ ){
      HepMC::GenVertex* GenVtx = evt->barcode_to_vertex(-iv);
      //
      // now find decay products with funky barcode, weed out and replace with clones of sensible barcode
      // we can NOT change the barcode while iterating, because iterators do depend on the barcoding
      // thus we have to take a 2-step procedure
      //
      std::vector<int> BCodes;
      BCodes.clear();
      for (HepMC::GenVertex::particle_iterator pitr= GenVtx->particles_begin(HepMC::children);
       pitr != GenVtx->particles_end(HepMC::children); ++pitr){
        if ( (*pitr)->barcode() > 10000 ){
      BCodes.push_back( (*pitr)->barcode() );
        }
      }
      if ( BCodes.size() > 0 ){
        for ( size_t ibc=0; ibc<BCodes.size(); ibc++ ){
      HepMC::GenParticle* p1 = evt->barcode_to_particle( BCodes[ibc] );
      int nbc = p1->barcode() - 10000 + NPartBefore;
      p1->suggest_barcode( nbc );
        }
      }             
    }
 
    for (HepMC::GenEvent::particle_const_iterator p= evt->particles_begin(); p != evt->particles_end(); ++p){
      if((*p)->end_vertex() && (*p)->status() == 1)(*p)->set_status(2);
      if((*p)->end_vertex() && (*p)->end_vertex()->particles_out_size()==0) edm::LogWarning("TauolappInterface::decay error: empty end vertex!");
    }
 
    return evt;
 }

void TauolappInterface::decodeMDTAU ( int mdtau )

private

Definition at line 358 of file TauolappInterface.cc.

References fHadronModes, fLeptonModes, fScaledHadronBrRatios, fScaledLeptonBrRatios, fSelectDecayByEvent, gamprt, i, mathSSE::return(), and taubra_.

Referenced by init().

 {
 
    // Note-1:
    // I have to hack the common block directly because set<...>DecayMode(...)
    // only changes it in the Tauola++ instance but does NOT passes it over
    // to the Fortran core - this it does only one, via initialize() stuff...
    //
    // So I'll do both ways of settings, just for consistency...
    // but I probably need to communicate it to the Tauola(++) team...
    //
    
    // Note-2: 
    // originally, the 1xx settings are meant for tau's from hard event,
    // and the 2xx settings are for any tau in the event record;
    //
    // later one, we'll have to take this into account...
    // but first I'll have to sort out what happens in the 1xx case
    // to tau's coming outside of hard event (if any in the record)   
    //
    
    if ( mdtau == 101 || mdtau == 201 )
    {
       // override with electron mode for both tau's
       //
       Tauolapp::jaki_.jak1 = 1;
       Tauolapp::jaki_.jak2 = 1;
       Tauolapp::Tauola::setSameParticleDecayMode( 1 ) ;
       Tauolapp::Tauola::setOppositeParticleDecayMode( 1 ) ;
       return;
    }
    
    if ( mdtau == 102 || mdtau == 202 )
    {
       // override with muon mode for both tau's
       //
       Tauolapp::jaki_.jak1 = 2;
       Tauolapp::jaki_.jak2 = 2;
       Tauolapp::Tauola::setSameParticleDecayMode( 2 ) ;
       Tauolapp::Tauola::setOppositeParticleDecayMode( 2 ) ;
       return;
    }
 
    if ( mdtau == 111 || mdtau == 211 )
    {
       // override with electron mode for 1st tau 
       // and any mode for 2nd tau
       //
       Tauolapp::jaki_.jak1 = 1;
       Tauolapp::jaki_.jak2 = 0;
       Tauolapp::Tauola::setSameParticleDecayMode( 1 ) ;
       Tauolapp::Tauola::setOppositeParticleDecayMode( 0 ) ;
       return;
    }
 
    if ( mdtau == 112 || mdtau == 212 )
    {
       // override with muon mode for the 1st tau 
       // and any mode for the 2nd tau
       //
       Tauolapp::jaki_.jak1 = 2;
       Tauolapp::jaki_.jak2 = 0;
       Tauolapp::Tauola::setSameParticleDecayMode( 2 ) ;
       Tauolapp::Tauola::setOppositeParticleDecayMode( 0 ) ;
       return;
    }
    
    if ( mdtau == 121 || mdtau == 221 )
    {
       // override with any mode for the 1st tau 
       // and electron mode for the 2nd tau
       //
       Tauolapp::jaki_.jak1 = 0;
       Tauolapp::jaki_.jak2 = 1;
       Tauolapp::Tauola::setSameParticleDecayMode( 0 ) ;
       Tauolapp::Tauola::setOppositeParticleDecayMode( 1 ) ;
       return;
    }
    
    if ( mdtau == 122 || mdtau == 222 )
    {
       // override with any mode for the 1st tau 
       // and muon mode for the 2nd tau
       //
       Tauolapp::jaki_.jak1 = 0;
       Tauolapp::jaki_.jak2 = 2;
       Tauolapp::Tauola::setSameParticleDecayMode( 0 ) ;
       Tauolapp::Tauola::setOppositeParticleDecayMode( 2 ) ;
       return;
    }
 
    if ( mdtau == 140 || mdtau == 240 )
    {
       // override with pi+/- nutau mode for both tau's 
       //
       Tauolapp::jaki_.jak1 = 3;
       Tauolapp::jaki_.jak2 = 3;
       Tauolapp::Tauola::setSameParticleDecayMode( 3 ) ;
       Tauolapp::Tauola::setOppositeParticleDecayMode( 3 ) ;
       return;
    }
 
    if ( mdtau == 141 || mdtau == 241 )
    {
       // override with pi+/- nutau mode for the 1st tau 
       // and any mode for the 2nd tau
       //
       Tauolapp::jaki_.jak1 = 3;
       Tauolapp::jaki_.jak2 = 0;
       Tauolapp::Tauola::setSameParticleDecayMode( 3 ) ;
       Tauolapp::Tauola::setOppositeParticleDecayMode( 0 ) ;
       return;
    }
 
    if ( mdtau == 142 || mdtau == 242 )
    {
       // override with any mode for the 1st tau 
       // and pi+/- nutau mode for 2nd tau
       //
       Tauolapp::jaki_.jak1 = 0;
       Tauolapp::jaki_.jak2 = 3;
       Tauolapp::Tauola::setSameParticleDecayMode( 0 ) ;
       Tauolapp::Tauola::setOppositeParticleDecayMode( 3 ) ;
       return;
    }
    
    // OK, we come here for semi-inclusive modes
    //
    
    // First of all, leptons and hadron modes sums
    //
    // re-scale branching ratios, just in case...
    //
    double sumBra = 0;
    
    // the number of decay modes is hardcoded at 22 because that's what it is right now in Tauola
    // in the future, perhaps an asscess method would be useful - communicate to Tauola team...
    //
    
    for ( int i=0; i<22; i++ )
    {
       sumBra += Tauolapp::taubra_.gamprt[i];
    }
    if ( sumBra == 0. ) return ; // perhaps need to throw ?
    for ( int i=0; i<22; i++ )
    {
       double newBra = Tauolapp::taubra_.gamprt[i] / sumBra;
       Tauolapp::Tauola::setTauBr( i+1, newBra ); 
    }
    sumBra = 1.0;
    
    double sumLeptonBra = Tauolapp::taubra_.gamprt[0] + Tauolapp::taubra_.gamprt[1];
    double sumHadronBra = sumBra - sumLeptonBra;
    
    for ( int i=0; i<2; i++ )
    {
       fLeptonModes.push_back( i+1 );
       fScaledLeptonBrRatios.push_back( (Tauolapp::taubra_.gamprt[i]/sumLeptonBra) );  
    }
    for ( int i=2; i<22; i++ )
    {
       fHadronModes.push_back( i+1 );
       fScaledHadronBrRatios.push_back( (Tauolapp::taubra_.gamprt[i]/sumHadronBra) ); 
    }
 
    fSelectDecayByEvent = true;
    return;
       
 }

void gen::TauolappInterface::disablePolarization ( )

inlinevirtual

Reimplemented from gen::TauolaInterfaceBase.

Definition at line 38 of file TauolappInterface.h.

References fPolarization.

38 { fPolarization = false; return; }

gen::TauolappInterface::fPolarization

bool fPolarization

Definition: TauolappInterface.h:65

void gen::TauolappInterface::enablePolarization ( )

inlinevirtual

Reimplemented from gen::TauolaInterfaceBase.

Definition at line 37 of file TauolappInterface.h.

References fPolarization.

37 { fPolarization = true; return; }

gen::TauolappInterface::fPolarization

bool fPolarization

Definition: TauolappInterface.h:65

HepMC::GenParticle * TauolappInterface::FirstTauInChain ( HepMC::GenParticle * tau )

private

Definition at line 768 of file TauolappInterface.cc.

References parents, and metsig::tau.

Referenced by MatchedLHESpinUp().

                                                                          {
   if ( tau->production_vertex() ) {
     HepMC::GenVertex::particle_iterator mother;
     for (mother = tau->production_vertex()->particles_begin(HepMC::parents); mother!= tau->production_vertex()->particles_end(HepMC::parents); mother++ ) {
       if((*mother)->pdg_id()==tau->pdg_id()) return FirstTauInChain(*mother); // recursive call to get mother with different pdgid
     }
   }
   return tau;
 }

double TauolappInterface::flat ( void )

static

Definition at line 190 of file TauolappInterface.cc.

References edm::hlt::Exception, and fRandomEngine.

Referenced by decay(), init(), gen::ranmar_(), and selectLeptonic().

 {
    if ( !fRandomEngine ) {
      throw cms::Exception("LogicError")
        << "TauolaInterface::flat: Attempt to generate random number when engine pointer is null\n"
        << "This might mean that the code was modified to generate a random number outside the\n"
        << "event and beginLuminosityBlock methods, which is not allowed.\n";
    }
    return fRandomEngine->flat();
 }

HepMC::GenParticle * TauolappInterface::GetMother ( HepMC::GenParticle * tau )

private

Definition at line 778 of file TauolappInterface.cc.

References parents, and metsig::tau.

Referenced by decay(), and MatchedLHESpinUp().

                                                                    {
   if ( tau->production_vertex() ) {
     HepMC::GenVertex::particle_iterator mother;
     for (mother = tau->production_vertex()->particles_begin(HepMC::parents); mother!= tau->production_vertex()->particles_end(HepMC::parents); mother++ ) {
       if((*mother)->pdg_id() == tau->pdg_id()) return GetMother(*mother); // recursive call to get mother with different pdgid 
       return (*mother);
     }
   }
   return tau;
 }

void TauolappInterface::init ( const edm::EventSetup & es )

virtual

Reimplemented from gen::TauolaInterfaceBase.

Definition at line 73 of file TauolappInterface.cc.

References funct::abs(), looper::cfg, gather_cfg::cout, decodeMDTAU(), dmMatch, dolhe, dolheBosonCorr, edm::hlt::Exception, edm::ParameterSet::exists(), fIsInitialized, flat(), fMDTAU, fPDGs, fPDGTable, fPolarization, fPSet, edm::EventSetup::getData(), edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), i, initialize(), lifetime, ntries, RecoTau_DiTaus_pt_20-420_cfg::ParticleID, and AlCaHLTBitMon_QueryRunRegistry::string.

                                                      {
   if ( fIsInitialized ) return; // do init only once
   if ( fPSet == 0 ) throw cms::Exception("TauolappInterfaceError") << "Attempt to initialize Tauola with an empty ParameterSet\n" << std::endl;   
   
   fIsInitialized = true;
       
    es.getData( fPDGTable ) ;
 
    Tauolapp::Tauola::setDecayingParticle(15);
 
    // LHE Information
    dmMatch=fPSet->getUntrackedParameter<double>("dmMatch",0.5);
    dolhe=fPSet->getUntrackedParameter<bool>("dolhe",false);
    dolheBosonCorr=fPSet->getUntrackedParameter<bool>("dolheBosonCorr",true);
    ntries=fPSet->getUntrackedParameter<int>("ntries",10);
 
    // polarization switch 
    // fPolarization = fPSet->getParameter<bool>("UseTauolaPolarization") ? 1 : 0 ;
    fPolarization = fPSet->getParameter<bool>("UseTauolaPolarization");
    
    // read tau decay mode switches
    //
    ParameterSet cards = fPSet->getParameter< ParameterSet >("InputCards");
    
    fMDTAU = cards.getParameter< int >( "mdtau" );
 
    if ( fMDTAU == 0 || fMDTAU == 1 )
    {
       Tauolapp::Tauola::setSameParticleDecayMode( cards.getParameter< int >( "pjak1" ) ) ;
       Tauolapp::Tauola::setOppositeParticleDecayMode( cards.getParameter< int >( "pjak2" ) ) ;
    }
 
    Tauolapp::Tauola::spin_correlation.setAll(fPolarization);
 
    // some more options, copied over from an example 
    // Default values
    //Tauola::setEtaK0sPi(0,0,0); // switches to decay eta K0_S and pi0 1/0 on/off. 
  
    const HepPDT::ParticleData* PData = fPDGTable->particle(HepPDT::ParticleID( abs(Tauolapp::Tauola::getDecayingParticle()) )) ;
    double lifetime = PData->lifetime().value();
    Tauolapp::Tauola::setTauLifetime( lifetime );
 
    fPDGs.push_back( Tauolapp::Tauola::getDecayingParticle() );
 
    Tauolapp::Tauola::setRandomGenerator(gen::TauolappInterface::flat);         
 
    if (fPSet->exists("parameterSets")){
      std::vector<std::string> par = fPSet->getParameter< std::vector<std::string> >("parameterSets");
      for (unsigned int ip=0; ip<par.size(); ++ip ){
        std::string curSet = par[ip];
        if(curSet=="setNewCurrents") Tauolapp::Tauola::setNewCurrents(fPSet->getParameter<int>(curSet));
      }
    }
 
    Tauolapp::Tauola::initialize();
 
    Tauolapp::Tauola::spin_correlation.setAll(fPolarization);// Tauola switches this on during Tauola::initialise(); so we add this here to keep it on/off
 
    if (fPSet->exists("parameterSets")){
      std::vector<std::string> par = fPSet->getParameter< std::vector<std::string> >("parameterSets");
      for (unsigned int ip=0; ip<par.size(); ++ip ){
        std::string curSet = par[ip];
        if(curSet=="spinCorrelationSetAll") Tauolapp::Tauola::spin_correlation.setAll(fPSet->getParameter<bool>(curSet));
        if(curSet=="spinCorrelationGAMMA") Tauolapp::Tauola::spin_correlation.GAMMA=fPSet->getParameter<bool>(curSet);
        if(curSet=="spinCorrelationZ0") Tauolapp::Tauola::spin_correlation.Z0=fPSet->getParameter<bool>(curSet);
        if(curSet=="spinCorrelationHIGGS") Tauolapp::Tauola::spin_correlation.HIGGS=fPSet->getParameter<bool>(curSet);
        if(curSet=="spinCorrelationHIGGSH") Tauolapp::Tauola::spin_correlation.HIGGS_H=fPSet->getParameter<bool>(curSet);
        if(curSet=="spinCorrelationHIGGSA") Tauolapp::Tauola::spin_correlation.HIGGS_A=fPSet->getParameter<bool>(curSet);
        if(curSet=="spinCorrelationHIGGSPLUS") Tauolapp::Tauola::spin_correlation.HIGGS_PLUS=fPSet->getParameter<bool>(curSet);
        if(curSet=="spinCorrelationHIGGSMINUS") Tauolapp::Tauola::spin_correlation.HIGGS_MINUS=fPSet->getParameter<bool>(curSet);
        if(curSet=="spinCorrelationWPLUS") Tauolapp::Tauola::spin_correlation.W_PLUS=fPSet->getParameter<bool>(curSet);
        if(curSet=="spinCorrelationWMINUS") Tauolapp::Tauola::spin_correlation.W_MINUS=fPSet->getParameter<bool>(curSet);
        
        if(curSet=="setHiggsScalarPseudoscalarPDG") Tauolapp::Tauola::setHiggsScalarPseudoscalarPDG(fPSet->getParameter<int>(curSet));
        if(curSet=="setHiggsScalarPseudoscalarMixingAngle") Tauolapp::Tauola::setHiggsScalarPseudoscalarMixingAngle(fPSet->getParameter<double>(curSet));
        
        if(curSet=="setRadiation") Tauolapp::Tauola::setRadiation(fPSet->getParameter<bool>(curSet));
        if(curSet=="setRadiationCutOff") Tauolapp::Tauola::setRadiationCutOff(fPSet->getParameter<double>(curSet));
        
        if(curSet=="setEtaK0sPi"){
      std::vector<int> vpar = fPSet->getParameter<std::vector<int> >(curSet);
      if(vpar.size()==3) Tauolapp::Tauola::setEtaK0sPi(vpar[0],vpar[1],vpar[2]);
      else {std::cout << "WARNING invalid size for setEtaK0sPi: " << vpar.size() << " Require 3 elements " << std::endl;}
        }
        
        if(curSet=="setTaukle"){
      std::vector<double> vpar = fPSet->getParameter<std::vector<double> >(curSet);
      if(vpar.size()==4) Tauolapp::Tauola::setTaukle(vpar[0], vpar[1], vpar[2], vpar[3]);
      else {std::cout << "WARNING invalid size for setTaukle: " << vpar.size() << " Require 4 elements " << std::endl;}
        }
        
        if(curSet=="setTauBr"){
      edm::ParameterSet cfg = fPSet->getParameter<edm::ParameterSet>(curSet);
      std::vector<int> vJAK = cfg.getParameter<std::vector<int> >("JAK");
      std::vector<double> vBR = cfg.getParameter<std::vector<double> >("BR");
      if(vJAK.size() == vBR.size()){
        for(unsigned int i=0;i<vJAK.size();i++) Tauolapp::Tauola::setTauBr(vJAK[i],vBR[i]);
      }
      else {std::cout << "WARNING invalid size for setTauBr - JAK: " << vJAK.size() << " BR: " << vBR.size() << std::endl;}
        }
      }
    }
 
    // override decay modes if needs be
    //
    // we have to do it AFTER init because otherwises branching ratios are NOT filled in
    //
    if ( fMDTAU != 0 && fMDTAU != 1 )
    {
       decodeMDTAU( fMDTAU );
    }
 
    Tauolapp::Log::LogWarning(false);
    
    return;
 }

bool TauolappInterface::isLastTauInChain ( const HepMC::GenParticle * tau )

private

Definition at line 740 of file TauolappInterface.cc.

Referenced by decay().

                                                                    {
   if ( tau->end_vertex() ) {
     HepMC::GenVertex::particle_iterator dau;
     for (dau = tau->end_vertex()->particles_begin(HepMC::children); dau!= tau->end_vertex()->particles_end(HepMC::children); dau++ ) {
       int dau_pid = (*dau)->pdg_id();
       if(dau_pid == tau->pdg_id()) return false;
     }
   }
   return true;
 }

HepMC::GenEvent * TauolappInterface::make_simple_tau_event	(	const TLorentzVector &	l,
		int	pdgid,
		int	status
	)

private

Definition at line 704 of file TauolappInterface.cc.

References event(), and GenParticle::GenParticle.

Referenced by decay().

                                                                                                    {
   HepMC::GenEvent *event = new HepMC::GenEvent();
   // make tau's four vector
   HepMC::FourVector momentum_tau1(l.Px(),l.Py(),l.Pz(),l.E());
   // make particles
   HepMC::GenParticle *tau1 = new HepMC::GenParticle(momentum_tau1,pdgid,status);
   // make the vertex
   HepMC::GenVertex *vertex = new HepMC::GenVertex();
   vertex->add_particle_out(tau1);
   event->add_vertex(vertex);
   return event;
 }

double TauolappInterface::MatchedLHESpinUp	(	HepMC::GenParticle *	tau,
		std::vector< HepMC::GenParticle > &	p,
		std::vector< double > &	spinup,
		std::vector< int > &	m_idx
	)

private

Definition at line 751 of file TauolappInterface.cc.

References dmMatch, FirstTauInChain(), GenParticle::GenParticle, GetMother(), i, visualization-live-secondInstance_cfg::m, and edmStreamStallGrapher::t.

Referenced by decay().

                                              {
   HepMC::GenParticle* Tau=FirstTauInChain(tau);
   HepMC::GenParticle* mother=GetMother(Tau);
   TLorentzVector t(tau->momentum().px(),tau->momentum().py(),tau->momentum().pz(),tau->momentum().e());
   TLorentzVector m(mother->momentum().px(),mother->momentum().py(),mother->momentum().pz(),mother->momentum().e());
   for(unsigned int i=0;i<p.size();i++){
     if(tau->pdg_id()==p.at(i).pdg_id()){
       if(mother->pdg_id()==p.at(m_idx.at(i)).pdg_id()){
     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());
     if(fabs(m.M()-pm.M())<dmMatch)return spinup.at(i);
       }
     }
   }
   return 0;
 }

const std::vector<int>& gen::TauolappInterface::operatesOnParticles ( )

inlinevirtual

Reimplemented from gen::TauolaInterfaceBase.

Definition at line 40 of file TauolappInterface.h.

References fPDGs.

40 { return fPDGs; }

gen::TauolappInterface::fPDGs

std::vector< int > fPDGs

Definition: TauolappInterface.h:64

void TauolappInterface::selectDecayByMDTAU ( )

private

Definition at line 528 of file TauolappInterface.cc.

References fMDTAU, jak2, alignBH_cfg::mode, mathSSE::return(), selectHadronic(), and selectLeptonic().

Referenced by decay().

 {
 
    
    if ( fMDTAU == 100 || fMDTAU == 200 )
    {
       int mode = selectLeptonic();
       Tauolapp::jaki_.jak1 = mode;
       Tauolapp::Tauola::setSameParticleDecayMode( mode );
       mode = selectLeptonic();
       Tauolapp::jaki_.jak2 = mode;
       Tauolapp::Tauola::setOppositeParticleDecayMode( mode );
       return ;
    }
    
    int modeL = selectLeptonic();
    int modeH = selectHadronic();
    
    if ( fMDTAU == 110 || fMDTAU == 210 )
    {
       Tauolapp::jaki_.jak1 = modeL;
       Tauolapp::jaki_.jak2 = 0;
       Tauolapp::Tauola::setSameParticleDecayMode( modeL );
       Tauolapp::Tauola::setOppositeParticleDecayMode( 0 );
       return ;
    }
    
    if ( fMDTAU == 120 || fMDTAU == 22 )
    {
       Tauolapp::jaki_.jak1 = 0;
       Tauolapp::jaki_.jak2 = modeL;
       Tauolapp::Tauola::setSameParticleDecayMode( 0 );
       Tauolapp::Tauola::setOppositeParticleDecayMode( modeL );
       return;      
    }
    
    if ( fMDTAU == 114 || fMDTAU == 214 )
    {
       Tauolapp::jaki_.jak1 = modeL;
       Tauolapp::jaki_.jak2 = modeH;
       Tauolapp::Tauola::setSameParticleDecayMode( modeL );
       Tauolapp::Tauola::setOppositeParticleDecayMode( modeH );
       return;      
    }
 
    if ( fMDTAU == 124 || fMDTAU == 224 )
    {
       Tauolapp::jaki_.jak1 = modeH;
       Tauolapp::jaki_.jak2 = modeL;
       Tauolapp::Tauola::setSameParticleDecayMode( modeH );
       Tauolapp::Tauola::setOppositeParticleDecayMode( modeL );
       return;      
    }
 
    if ( fMDTAU == 115 || fMDTAU == 215 )
    {
       Tauolapp::jaki_.jak1 = 1;
       Tauolapp::jaki_.jak2 = modeH;
       Tauolapp::Tauola::setSameParticleDecayMode( 1 );
       Tauolapp::Tauola::setOppositeParticleDecayMode( modeH );
       return;      
    }
 
    if ( fMDTAU == 125 || fMDTAU == 225 )
    {
       Tauolapp::jaki_.jak1 = modeH;
       Tauolapp::jaki_.jak2 = 1;
       Tauolapp::Tauola::setSameParticleDecayMode( modeH );
       Tauolapp::Tauola::setOppositeParticleDecayMode( 1 );
       return;      
    }
 
    if ( fMDTAU == 116 || fMDTAU == 216 )
    {
       Tauolapp::jaki_.jak1 = 2;
       Tauolapp::jaki_.jak2 = modeH;
       Tauolapp::Tauola::setSameParticleDecayMode( 2 );
       Tauolapp::Tauola::setOppositeParticleDecayMode( modeH );
       return;      
    }
 
    if ( fMDTAU == 126 || fMDTAU == 226 )
    {
       Tauolapp::jaki_.jak1 = modeH;
       Tauolapp::jaki_.jak2 = 2;
       Tauolapp::Tauola::setSameParticleDecayMode( modeH );
       Tauolapp::Tauola::setOppositeParticleDecayMode( 2 );
       return;      
    }
 
    if ( fMDTAU == 130 || fMDTAU == 230 )
    {
       Tauolapp::jaki_.jak1 = modeH;
       Tauolapp::jaki_.jak2 = selectHadronic();
       Tauolapp::Tauola::setSameParticleDecayMode( modeH );
       Tauolapp::Tauola::setOppositeParticleDecayMode( Tauolapp::jaki_.jak2 );
       return;      
    }
 
    if ( fMDTAU == 131 || fMDTAU == 231 )
    {
       Tauolapp::jaki_.jak1 = modeH;
       Tauolapp::jaki_.jak2 = 0;
       Tauolapp::Tauola::setSameParticleDecayMode( modeH );
       Tauolapp::Tauola::setOppositeParticleDecayMode( 0 );
       return;      
    }
 
    if ( fMDTAU == 132 || fMDTAU == 232 )
    {
       Tauolapp::jaki_.jak1 = 0;
       Tauolapp::jaki_.jak2 = modeH;
       Tauolapp::Tauola::setSameParticleDecayMode( 0 );
       Tauolapp::Tauola::setOppositeParticleDecayMode( modeH );
       return;      
    }
    
    // unlikely that we get here on unknown mdtau 
    // - there's a protection earlier
    // but if we do, just set defaults
    // probably need to spit a warning...
    //
    Tauolapp::Tauola::setSameParticleDecayMode( 0 );
    Tauolapp::Tauola::setOppositeParticleDecayMode( 0 );
       
    return;
    
 
 }

int TauolappInterface::selectHadronic ( )

private

Definition at line 675 of file TauolappInterface.cc.

References fHadronModes, fScaledHadronBrRatios, i, and gen::ranmar_().

Referenced by selectDecayByMDTAU().

 {
 
    float prob = 0.;
    int len = 1;
    ranmar_(&prob,&len);
    
    double sumBra = fScaledHadronBrRatios[0];
    if ( prob > 0. && prob <= sumBra ) 
    {
       return fHadronModes[0];
    }
    else
    {
       int NN = fScaledHadronBrRatios.size();
       for ( int i=1; i<NN; i++ )
       {
          if ( prob > sumBra && prob <= (sumBra+fScaledHadronBrRatios[i]) ) 
      {
         return fHadronModes[i];
      }
      sumBra += fScaledHadronBrRatios[i];
       }
    }
    
    return 0;
 
 }

int TauolappInterface::selectLeptonic ( )

private

Definition at line 658 of file TauolappInterface.cc.

References flat(), and fScaledLeptonBrRatios.

Referenced by selectDecayByMDTAU().

 {
    
    float prob = flat();
    
    if ( prob > 0. && prob <= fScaledLeptonBrRatios[0] ) 
    {
       return 1;
    }
    else if ( prob > fScaledLeptonBrRatios[1] && prob <=1. )
    {
       return 2;
    }
       
    return 0;
 }

virtual void gen::TauolappInterface::SetLHE ( lhef::LHEEvent * l )

inlinevirtual

Reimplemented from gen::TauolaInterfaceBase.

Definition at line 43 of file TauolappInterface.h.

References prof2calltree::l, and lhe.

43 {lhe=l;}

gen::TauolappInterface::lhe

lhef::LHEEvent * lhe

Definition: TauolappInterface.h:76

prof2calltree.l

tuple l

Definition: prof2calltree.py:116

void gen::TauolappInterface::setRandomEngine ( CLHEP::HepRandomEngine * v )

inlinevirtual

Implements gen::TauolaInterfaceBase.

Definition at line 44 of file TauolappInterface.h.

References fRandomEngine, and gen::v.

44 { fRandomEngine = v; }

gen::v

double v[5][pyjets_maxn]

Definition: Cascade2Hadronizer.cc:80

gen::TauolappInterface::fRandomEngine

static CLHEP::HepRandomEngine * fRandomEngine

Definition: TauolappInterface.h:63

void TauolappInterface::statistics ( )

virtual

Reimplemented from gen::TauolaInterfaceBase.

Definition at line 353 of file TauolappInterface.cc.

 {
    return;
 }

void TauolappInterface::update_particles	(	HepMC::GenParticle *	partHep,
		HepMC::GenEvent *	theEvent,
		HepMC::GenParticle *	p,
		TVector3 &	boost
	)

private

Definition at line 717 of file TauolappInterface.cc.

References ztail::d, GenParticle::GenParticle, and prof2calltree::l.

Referenced by decay().

                                                                                                                              {
   partHep->set_status(p->status());
   if(p->end_vertex()){
     if(!partHep->end_vertex()){
       HepMC::GenVertex* vtx = new HepMC::GenVertex(p->end_vertex()->position());
       theEvent->add_vertex(vtx);
       vtx->add_particle_in(partHep);
     }
     if(p->end_vertex()->particles_out_size()!=0){
       for(HepMC::GenVertex::particles_out_const_iterator d=p->end_vertex()->particles_out_const_begin(); d!=p->end_vertex()->particles_out_const_end();d++){
     // Create daughter and add to event
     TLorentzVector l((*d)->momentum().px(),(*d)->momentum().py(),(*d)->momentum().pz(),(*d)->momentum().e());
     l.Boost(boost);
     HepMC::FourVector momentum(l.Px(),l.Py(),l.Pz(),l.E());
     HepMC::GenParticle *daughter = new HepMC::GenParticle(momentum,(*d)->pdg_id(),(*d)->status());
     daughter->suggest_barcode(theEvent->particles_size()+1);
     partHep->end_vertex()->add_particle_out(daughter);
     if((*d)->end_vertex()) update_particles(daughter,theEvent,(*d),boost);
       }
     }
   }
 }

Member Data Documentation

double gen::TauolappInterface::dmMatch

private

Definition at line 78 of file TauolappInterface.h.

Referenced by init(), and MatchedLHESpinUp().

bool gen::TauolappInterface::dolhe

private

Definition at line 79 of file TauolappInterface.h.

Referenced by decay(), and init().

bool gen::TauolappInterface::dolheBosonCorr

private

Definition at line 80 of file TauolappInterface.h.

Referenced by decay(), and init().

std::vector<int> gen::TauolappInterface::fHadronModes

private

Definition at line 73 of file TauolappInterface.h.

Referenced by decodeMDTAU(), and selectHadronic().

bool gen::TauolappInterface::fIsInitialized

private

Definition at line 68 of file TauolappInterface.h.

Referenced by decay(), and init().

std::vector<int> gen::TauolappInterface::fLeptonModes

private

Definition at line 72 of file TauolappInterface.h.

Referenced by decodeMDTAU().

int gen::TauolappInterface::fMDTAU

private

Definition at line 70 of file TauolappInterface.h.

Referenced by init(), and selectDecayByMDTAU().

std::vector<int> gen::TauolappInterface::fPDGs

private

Definition at line 64 of file TauolappInterface.h.

Referenced by init(), and operatesOnParticles().

edm::ESHandle<HepPDT::ParticleDataTable> gen::TauolappInterface::fPDGTable

private

Definition at line 66 of file TauolappInterface.h.

Referenced by init().

bool gen::TauolappInterface::fPolarization

private

Definition at line 65 of file TauolappInterface.h.

Referenced by disablePolarization(), enablePolarization(), and init().

edm::ParameterSet* gen::TauolappInterface::fPSet

private

Definition at line 67 of file TauolappInterface.h.

Referenced by init(), TauolappInterface(), and ~TauolappInterface().

CLHEP::HepRandomEngine * TauolappInterface::fRandomEngine = 0

staticprivate

Definition at line 63 of file TauolappInterface.h.

Referenced by flat(), and setRandomEngine().

std::vector<double> gen::TauolappInterface::fScaledHadronBrRatios

private

Definition at line 75 of file TauolappInterface.h.

Referenced by decodeMDTAU(), and selectHadronic().

std::vector<double> gen::TauolappInterface::fScaledLeptonBrRatios

private

Definition at line 74 of file TauolappInterface.h.

Referenced by decodeMDTAU(), and selectLeptonic().

bool gen::TauolappInterface::fSelectDecayByEvent

private

Definition at line 71 of file TauolappInterface.h.

Referenced by decay(), and decodeMDTAU().

lhef::LHEEvent* gen::TauolappInterface::lhe

private

Definition at line 76 of file TauolappInterface.h.

Referenced by decay(), and SetLHE().

double gen::TauolappInterface::lifetime

private

Definition at line 82 of file TauolappInterface.h.

Referenced by init().

int gen::TauolappInterface::ntries

private

Definition at line 81 of file TauolappInterface.h.

Referenced by decay(), and init().

Public Member Functions

Static Public Member Functions

Private Member Functions

Private Attributes

Static Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation