TauolappInterface.cc

Go to the documentation of this file.

/* this is the code for the new Tauola++ */

#include <iostream>

#include "GeneratorInterface/TauolaInterface/interface/TauolappInterface.h"

#include "Tauola/Tauola.h"
#include "Tauola/TauolaHepMCEvent.h"
#include "Tauola/Log.h"
#include "Tauola/TauolaHepMCParticle.h"
#include "Tauola/TauolaParticle.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "FWCore/Utilities/interface/RandomNumberGenerator.h"
#include "FWCore/Utilities/interface/Exception.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "CLHEP/Random/RandomEngine.h"

#include "HepMC/GenEvent.h"
#include "HepMC/IO_HEPEVT.h"
#include "HepMC/HEPEVT_Wrapper.h"

#include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"

// LHE Run
#include "SimDataFormats/GeneratorProducts/interface/LHERunInfoProduct.h"
#include "GeneratorInterface/LHEInterface/interface/LHERunInfo.h"

// LHE Event
#include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h"
#include "GeneratorInterface/LHEInterface/interface/LHEEvent.h"

using namespace gen;
using namespace edm;
using namespace std;

CLHEP::HepRandomEngine*  TauolappInterface::fRandomEngine = nullptr;

extern "C" {

  void gen::ranmar_( float *rvec, int *lenv )
  {
    for(int i = 0; i < *lenv; i++)
      *rvec++ = TauolappInterface::flat();
    return;
  }
  
  void gen::rmarin_( int*, int*, int* )
  {
    return;
  }
  
}

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

TauolappInterface::~TauolappInterface(){if ( fPSet != nullptr ) delete fPSet;}

void TauolappInterface::init( const edm::EventSetup& es ){
  if ( fIsInitialized ) return; // do init only once
  if ( fPSet == nullptr ) 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;
}

double TauolappInterface::flat()
{
   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::GenEvent* TauolappInterface::decay( HepMC::GenEvent* evt ){
   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!=nullptr){
     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.empty() ){
       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::statistics()
{
   return;
}

void TauolappInterface::decodeMDTAU( int mdtau )
{

   // 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 TauolappInterface::selectDecayByMDTAU()
{

   
   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::selectLeptonic()
{
   
   float prob = flat();
   
   if ( prob > 0. && prob <= fScaledLeptonBrRatios[0] ) 
   {
      return 1;
   }
   else if ( prob > fScaledLeptonBrRatios[1] && prob <=1. )
   {
      return 2;
   }
      
   return 0;
}

int TauolappInterface::selectHadronic()
{

   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;

}

HepMC::GenEvent* TauolappInterface::make_simple_tau_event(const TLorentzVector &l,int pdgid,int status){
  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;
}

void TauolappInterface::update_particles(HepMC::GenParticle* partHep,HepMC::GenEvent* theEvent,HepMC::GenParticle* p,TVector3 &boost){
  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);
      }
    }
  }
}

bool TauolappInterface::isLastTauInChain(const HepMC::GenParticle* tau){
  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;
}

double TauolappInterface::MatchedLHESpinUp(HepMC::GenParticle* tau, std::vector<HepMC::GenParticle> &p, std::vector<double> &spinup,
                       std::vector<int> &m_idx){
  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;
}

HepMC::GenParticle* TauolappInterface::FirstTauInChain(HepMC::GenParticle* tau){
  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;
}

HepMC::GenParticle* TauolappInterface::GetMother(HepMC::GenParticle* tau){
  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::BoostProdToLabLifeTimeInDecays(HepMC::GenParticle* p,TLorentzVector &lab, TLorentzVector &prod){
  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 
    } 
  } 
}