gen::TauolappInterface Class Reference

#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 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.

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

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, 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, 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.

{ fPolarization = false; return; }

void gen::TauolappInterface::enablePolarization ( )

inlinevirtual

Reimplemented from gen::TauolaInterfaceBase.

Definition at line 37 of file TauolappInterface.h.

References fPolarization.

{ fPolarization = true; return; }

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 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, 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 tree::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.

{ return fPDGs; }

void TauolappInterface::selectDecayByMDTAU ( )

private

Definition at line 528 of file TauolappInterface.cc.

References fMDTAU, 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.

{lhe=l;}

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

inlinevirtual

Implements gen::TauolaInterfaceBase.

Definition at line 44 of file TauolappInterface.h.

References fRandomEngine, and gen::v.

{ fRandomEngine = v; }

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().