Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-10-19 10:28:30 

0001 /* this is the code for the new Tauola++ */
0002 
0003 #include <iostream>
0004 
0005 #include "GeneratorInterface/TauolaInterface/interface/TauolappInterface.h"
0006 
0007 #include "Tauola/Tauola.h"
0008 #include "Tauola/TauolaHepMCEvent.h"
0009 #include "Tauola/Log.h"
0010 #include "Tauola/TauolaHepMCParticle.h"
0011 #include "Tauola/TauolaParticle.h"
0012 #include "FWCore/ServiceRegistry/interface/Service.h"
0013 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0014 #include "FWCore/Utilities/interface/Exception.h"
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016 
0017 #include "CLHEP/Random/RandomEngine.h"
0018 
0019 #include "HepMC/GenEvent.h"
0020 #include "HepMC/IO_HEPEVT.h"
0021 #include "HepMC/HEPEVT_Wrapper.h"
0022 
0023 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
0024 
0025 // LHE Run
0026 #include "SimDataFormats/GeneratorProducts/interface/LHERunInfoProduct.h"
0027 #include "GeneratorInterface/LHEInterface/interface/LHERunInfo.h"
0028 
0029 // LHE Event
0030 #include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h"
0031 #include "GeneratorInterface/LHEInterface/interface/LHEEvent.h"
0032 
0033 using namespace gen;
0034 using namespace edm;
0035 using namespace std;
0036 
0037 CLHEP::HepRandomEngine* TauolappInterface::fRandomEngine = nullptr;
0038 
0039 extern "C" {
0040 
0041 void gen::ranmar_(float* rvec, int* lenv) {
0042   for (int i = 0; i < *lenv; i++)
0043     *rvec++ = TauolappInterface::flat();
0044   return;
0045 }
0046 
0047 void gen::rmarin_(int*, int*, int*) { return; }
0048 }
0049 
0050 TauolappInterface::TauolappInterface(const edm::ParameterSet& pset, edm::ConsumesCollector iCollector)
0051     : fPolarization(false),
0052       fPDGTableToken(iCollector.esConsumes<edm::Transition::BeginLuminosityBlock>()),
0053       fPSet(nullptr),
0054       fIsInitialized(false),
0055       fMDTAU(-1),
0056       fSelectDecayByEvent(false),
0057       lhe(nullptr),
0058       dmMatch(0.5),
0059       dolhe(false),
0060       dolheBosonCorr(false),
0061       ntries(10) {
0062   fPSet = new ParameterSet(pset);
0063 }
0064 
0065 TauolappInterface::~TauolappInterface() {
0066   if (fPSet != nullptr)
0067     delete fPSet;
0068 }
0069 
0070 void TauolappInterface::init(const edm::EventSetup& es) {
0071   if (fIsInitialized)
0072     return;  // do init only once
0073   if (fPSet == nullptr)
0074     throw cms::Exception("TauolappInterfaceError") << "Attempt to initialize Tauola with an empty ParameterSet\n"
0075                                                    << std::endl;
0076 
0077   fIsInitialized = true;
0078 
0079   fPDGTable = es.getHandle(fPDGTableToken);
0080 
0081   Tauolapp::Tauola::setDecayingParticle(15);
0082 
0083   // LHE Information
0084   dmMatch = fPSet->getUntrackedParameter<double>("dmMatch", 0.5);
0085   dolhe = fPSet->getUntrackedParameter<bool>("dolhe", false);
0086   dolheBosonCorr = fPSet->getUntrackedParameter<bool>("dolheBosonCorr", true);
0087   ntries = fPSet->getUntrackedParameter<int>("ntries", 10);
0088 
0089   // polarization switch
0090   // fPolarization = fPSet->getParameter<bool>("UseTauolaPolarization") ? 1 : 0 ;
0091   fPolarization = fPSet->getParameter<bool>("UseTauolaPolarization");
0092 
0093   // read tau decay mode switches
0094   //
0095   ParameterSet cards = fPSet->getParameter<ParameterSet>("InputCards");
0096 
0097   fMDTAU = cards.getParameter<int>("mdtau");
0098 
0099   if (fMDTAU == 0 || fMDTAU == 1) {
0100     Tauolapp::Tauola::setSameParticleDecayMode(cards.getParameter<int>("pjak1"));
0101     Tauolapp::Tauola::setOppositeParticleDecayMode(cards.getParameter<int>("pjak2"));
0102   }
0103 
0104   Tauolapp::Tauola::spin_correlation.setAll(fPolarization);
0105 
0106   // some more options, copied over from an example
0107   // Default values
0108   //Tauola::setEtaK0sPi(0,0,0); // switches to decay eta K0_S and pi0 1/0 on/off.
0109 
0110   const HepPDT::ParticleData* PData =
0111       fPDGTable->particle(HepPDT::ParticleID(abs(Tauolapp::Tauola::getDecayingParticle())));
0112   double lifetime = PData->lifetime().value();
0113   Tauolapp::Tauola::setTauLifetime(lifetime);
0114 
0115   fPDGs.push_back(Tauolapp::Tauola::getDecayingParticle());
0116 
0117   Tauolapp::Tauola::setRandomGenerator(gen::TauolappInterface::flat);
0118 
0119   if (fPSet->exists("parameterSets")) {
0120     std::vector<std::string> par = fPSet->getParameter<std::vector<std::string> >("parameterSets");
0121     for (unsigned int ip = 0; ip < par.size(); ++ip) {
0122       std::string curSet = par[ip];
0123       if (curSet == "setNewCurrents")
0124         Tauolapp::Tauola::setNewCurrents(fPSet->getParameter<int>(curSet));
0125     }
0126   }
0127 
0128   Tauolapp::Tauola::initialize();
0129 
0130   Tauolapp::Tauola::spin_correlation.setAll(
0131       fPolarization);  // Tauola switches this on during Tauola::initialise(); so we add this here to keep it on/off
0132 
0133   if (fPSet->exists("parameterSets")) {
0134     std::vector<std::string> par = fPSet->getParameter<std::vector<std::string> >("parameterSets");
0135     for (unsigned int ip = 0; ip < par.size(); ++ip) {
0136       std::string curSet = par[ip];
0137       if (curSet == "spinCorrelationSetAll")
0138         Tauolapp::Tauola::spin_correlation.setAll(fPSet->getParameter<bool>(curSet));
0139       if (curSet == "spinCorrelationGAMMA")
0140         Tauolapp::Tauola::spin_correlation.GAMMA = fPSet->getParameter<bool>(curSet);
0141       if (curSet == "spinCorrelationZ0")
0142         Tauolapp::Tauola::spin_correlation.Z0 = fPSet->getParameter<bool>(curSet);
0143       if (curSet == "spinCorrelationHIGGS")
0144         Tauolapp::Tauola::spin_correlation.HIGGS = fPSet->getParameter<bool>(curSet);
0145       if (curSet == "spinCorrelationHIGGSH")
0146         Tauolapp::Tauola::spin_correlation.HIGGS_H = fPSet->getParameter<bool>(curSet);
0147       if (curSet == "spinCorrelationHIGGSA")
0148         Tauolapp::Tauola::spin_correlation.HIGGS_A = fPSet->getParameter<bool>(curSet);
0149       if (curSet == "spinCorrelationHIGGSPLUS")
0150         Tauolapp::Tauola::spin_correlation.HIGGS_PLUS = fPSet->getParameter<bool>(curSet);
0151       if (curSet == "spinCorrelationHIGGSMINUS")
0152         Tauolapp::Tauola::spin_correlation.HIGGS_MINUS = fPSet->getParameter<bool>(curSet);
0153       if (curSet == "spinCorrelationWPLUS")
0154         Tauolapp::Tauola::spin_correlation.W_PLUS = fPSet->getParameter<bool>(curSet);
0155       if (curSet == "spinCorrelationWMINUS")
0156         Tauolapp::Tauola::spin_correlation.W_MINUS = fPSet->getParameter<bool>(curSet);
0157 
0158       if (curSet == "setHiggsScalarPseudoscalarPDG")
0159         Tauolapp::Tauola::setHiggsScalarPseudoscalarPDG(fPSet->getParameter<int>(curSet));
0160       if (curSet == "setHiggsScalarPseudoscalarMixingAngle")
0161         Tauolapp::Tauola::setHiggsScalarPseudoscalarMixingAngle(fPSet->getParameter<double>(curSet));
0162 
0163       if (curSet == "setRadiation")
0164         Tauolapp::Tauola::setRadiation(fPSet->getParameter<bool>(curSet));
0165       if (curSet == "setRadiationCutOff")
0166         Tauolapp::Tauola::setRadiationCutOff(fPSet->getParameter<double>(curSet));
0167 
0168       if (curSet == "setEtaK0sPi") {
0169         std::vector<int> vpar = fPSet->getParameter<std::vector<int> >(curSet);
0170         if (vpar.size() == 3)
0171           Tauolapp::Tauola::setEtaK0sPi(vpar[0], vpar[1], vpar[2]);
0172         else {
0173           std::cout << "WARNING invalid size for setEtaK0sPi: " << vpar.size() << " Require 3 elements " << std::endl;
0174         }
0175       }
0176 
0177       if (curSet == "setTaukle") {
0178         std::vector<double> vpar = fPSet->getParameter<std::vector<double> >(curSet);
0179         if (vpar.size() == 4)
0180           Tauolapp::Tauola::setTaukle(vpar[0], vpar[1], vpar[2], vpar[3]);
0181         else {
0182           std::cout << "WARNING invalid size for setTaukle: " << vpar.size() << " Require 4 elements " << std::endl;
0183         }
0184       }
0185 
0186       if (curSet == "setTauBr") {
0187         edm::ParameterSet cfg = fPSet->getParameter<edm::ParameterSet>(curSet);
0188         std::vector<int> vJAK = cfg.getParameter<std::vector<int> >("JAK");
0189         std::vector<double> vBR = cfg.getParameter<std::vector<double> >("BR");
0190         if (vJAK.size() == vBR.size()) {
0191           for (unsigned int i = 0; i < vJAK.size(); i++)
0192             Tauolapp::Tauola::setTauBr(vJAK[i], vBR[i]);
0193         } else {
0194           std::cout << "WARNING invalid size for setTauBr - JAK: " << vJAK.size() << " BR: " << vBR.size() << std::endl;
0195         }
0196       }
0197     }
0198   }
0199 
0200   // override decay modes if needs be
0201   //
0202   // we have to do it AFTER init because otherwises branching ratios are NOT filled in
0203   //
0204   if (fMDTAU != 0 && fMDTAU != 1) {
0205     decodeMDTAU(fMDTAU);
0206   }
0207 
0208   Tauolapp::Log::LogWarning(false);
0209 
0210   return;
0211 }
0212 
0213 double TauolappInterface::flat() {
0214   if (!fRandomEngine) {
0215     throw cms::Exception("LogicError")
0216         << "TauolaInterface::flat: Attempt to generate random number when engine pointer is null\n"
0217         << "This might mean that the code was modified to generate a random number outside the\n"
0218         << "event and beginLuminosityBlock methods, which is not allowed.\n";
0219   }
0220   return fRandomEngine->flat();
0221 }
0222 
0223 HepMC::GenEvent* TauolappInterface::decay(HepMC::GenEvent* evt) {
0224   if (!fIsInitialized)
0225     return evt;
0226   Tauolapp::Tauola::setRandomGenerator(
0227       gen::TauolappInterface::flat);  // rest tauola++ random number incase other modules use tauola++
0228   int NPartBefore = evt->particles_size();
0229   int NVtxBefore = evt->vertices_size();
0230 
0231   // what do we do if Hep::GenEvent size is larger than 10K ???
0232   // Tauola (& Photos, BTW) can only handle up to 10K via HEPEVT,
0233   // and in case of CMS, it's only up to 4K !!!
0234   // override decay mode if needs be
0235 
0236   if (fSelectDecayByEvent) {
0237     selectDecayByMDTAU();
0238   }
0239   if (dolhe && lhe != nullptr) {
0240     std::vector<HepMC::GenParticle> particles;
0241     std::vector<int> m_idx;
0242     std::vector<double> spinup = lhe->getHEPEUP()->SPINUP;
0243     std::vector<int> pdg = lhe->getHEPEUP()->IDUP;
0244     for (unsigned int i = 0; i < spinup.size(); i++) {
0245       particles.push_back(HepMC::GenParticle(HepMC::FourVector(lhe->getHEPEUP()->PUP.at(i)[0],
0246                                                                lhe->getHEPEUP()->PUP.at(i)[1],
0247                                                                lhe->getHEPEUP()->PUP.at(i)[2],
0248                                                                lhe->getHEPEUP()->PUP.at(i)[3]),
0249                                              lhe->getHEPEUP()->IDUP.at(i)));
0250       int status = lhe->getHEPEUP()->ISTUP.at(i);
0251       particles.at(particles.size() - 1).set_generated_mass(lhe->getHEPEUP()->PUP.at(i)[4]);
0252       particles.at(particles.size() - 1).set_status(status > 0 ? (status == 2 ? 3 : status) : 3);
0253       m_idx.push_back(lhe->getHEPEUP()->MOTHUP.at(i).first - 1);  // correct for fortran index offset
0254     }
0255     // match to taus in hepmc and identify mother of taus
0256     bool hastaus(false);
0257     std::vector<HepMC::GenParticle*> match;
0258     for (HepMC::GenEvent::particle_const_iterator iter = evt->particles_begin(); iter != evt->particles_end(); iter++) {
0259       if (abs((*iter)->pdg_id()) == 15) {
0260         hastaus = true;
0261         int mother_pid(0);
0262         // check imediate parent to avoid parent tau ie tau->taugamma
0263         for (HepMC::GenVertex::particle_iterator mother = (*iter)->production_vertex()->particles_begin(HepMC::parents);
0264              mother != (*iter)->production_vertex()->particles_end(HepMC::parents);
0265              mother++) {
0266           mother_pid = (*mother)->pdg_id();
0267           if (mother_pid != (*iter)->pdg_id()) {
0268             // match against lhe record
0269             if (abs(mother_pid) == 24 ||  // W
0270                 abs(mother_pid) == 37 ||  // H+/-
0271                 abs(mother_pid) == 23 ||  // Z
0272                 abs(mother_pid) == 22 ||  // gamma
0273                 abs(mother_pid) == 25 ||  // H0 SM
0274                 abs(mother_pid) == 35 ||  // H0
0275                 abs(mother_pid) == 36     // A0
0276             ) {
0277               bool isfound = false;
0278               for (unsigned int k = 0; k < match.size(); k++) {
0279                 if ((*mother) == match.at(k))
0280                   isfound = true;
0281               }
0282               if (!isfound)
0283                 match.push_back(*mother);
0284             }
0285           }
0286         }
0287       }
0288     }
0289     if (hastaus) {
0290       // if is single gauge boson decay and match helicities
0291       if (match.size() == 1 && dolheBosonCorr) {
0292         for (int i = 0; i < ntries; i++) {
0293           // re-decay taus then check if helicities match
0294           auto* t_event = new Tauolapp::TauolaHepMCEvent(evt);
0295           t_event->undecayTaus();
0296           t_event->decayTaus();
0297           bool ismatch = true;
0298           for (unsigned int j = 0; j < spinup.size(); j++) {
0299             if (abs(pdg.at(j)) == 15) {
0300               double diffhelminus = (-1.0 * (double)Tauolapp::Tauola::getHelMinus() -
0301                                      spinup.at(j));  // -1.0 to correct for tauola feature
0302               double diffhelplus = ((double)Tauolapp::Tauola::getHelPlus() - spinup.at(j));
0303               if (pdg.at(j) == 15 && diffhelminus > 0.5)
0304                 ismatch = false;
0305               if (pdg.at(j) == -15 && diffhelplus > 0.5)
0306                 ismatch = false;
0307             }
0308           }
0309           delete t_event;
0310           if (ismatch)
0311             break;
0312         }
0313       } else {
0314         // If the event does not contain a single gauge boson the code will be run with
0315         // remove all tau decays
0316         auto* t_event = new Tauolapp::TauolaHepMCEvent(evt);
0317         t_event->undecayTaus();
0318         delete t_event;
0319         // decay all taus manually based on the helicity
0320         for (HepMC::GenEvent::particle_const_iterator iter = evt->particles_begin(); iter != evt->particles_end();
0321              iter++) {
0322           if (abs((*iter)->pdg_id()) == 15 && isLastTauInChain(*iter)) {
0323             TLorentzVector ltau(
0324                 (*iter)->momentum().px(), (*iter)->momentum().py(), (*iter)->momentum().pz(), (*iter)->momentum().e());
0325             HepMC::GenParticle* m = GetMother(*iter);
0326             TLorentzVector mother(m->momentum().px(), m->momentum().py(), m->momentum().pz(), m->momentum().e());
0327             TVector3 boost = -1.0 * mother.BoostVector();  // boost into mother's CM frame
0328             TLorentzVector ltau_lab = ltau;
0329             ltau.Boost(boost);
0330             mother.Boost(boost);
0331             HepMC::GenEvent* tauevt = make_simple_tau_event(ltau, (*iter)->pdg_id(), (*iter)->status());
0332             HepMC::GenParticle* p = (*(tauevt->particles_begin()));
0333             Tauolapp::TauolaParticle* tp = new Tauolapp::TauolaHepMCParticle(p);
0334             double helicity = MatchedLHESpinUp(*iter, particles, spinup, m_idx);  // get helicity from lhe
0335             if ((*iter)->pdg_id() == 15)
0336               helicity *= -1.0;
0337             tp->undecay();
0338             // use |S_{tau}|=0.999999 to avoid issues with numerical roundoff
0339             Tauolapp::Tauola::decayOne(tp, true, 0, 0, ((double)helicity) * 0.999999);
0340             boost *= -1.0;  // boost back to lab frame
0341             mother.Boost(boost);
0342             update_particles((*iter), evt, p, boost);
0343             //correct tau liftetime for boost (change rest frame from mothers to taus)
0344             BoostProdToLabLifeTimeInDecays((*iter), ltau_lab, ltau);
0345             delete tauevt;
0346           }
0347         }
0348       }
0349     }
0350   } else {
0351     //construct tmp TAUOLA event
0352     auto* t_event = new Tauolapp::TauolaHepMCEvent(evt);
0353     //t_event->undecayTaus();
0354     t_event->decayTaus();
0355     delete t_event;
0356   }
0357 
0358   for (int iv = NVtxBefore + 1; iv <= evt->vertices_size(); iv++) {
0359     HepMC::GenVertex* GenVtx = evt->barcode_to_vertex(-iv);
0360     //
0361     // now find decay products with funky barcode, weed out and replace with clones of sensible barcode
0362     // we can NOT change the barcode while iterating, because iterators do depend on the barcoding
0363     // thus we have to take a 2-step procedure
0364     //
0365     std::vector<int> BCodes;
0366     BCodes.clear();
0367     for (HepMC::GenVertex::particle_iterator pitr = GenVtx->particles_begin(HepMC::children);
0368          pitr != GenVtx->particles_end(HepMC::children);
0369          ++pitr) {
0370       if ((*pitr)->barcode() > 10000) {
0371         BCodes.push_back((*pitr)->barcode());
0372       }
0373     }
0374     if (!BCodes.empty()) {
0375       for (size_t ibc = 0; ibc < BCodes.size(); ibc++) {
0376         HepMC::GenParticle* p1 = evt->barcode_to_particle(BCodes[ibc]);
0377         int nbc = p1->barcode() - 10000 + NPartBefore;
0378         p1->suggest_barcode(nbc);
0379       }
0380     }
0381   }
0382 
0383   for (HepMC::GenEvent::particle_const_iterator p = evt->particles_begin(); p != evt->particles_end(); ++p) {
0384     if ((*p)->end_vertex() && (*p)->status() == 1)
0385       (*p)->set_status(2);
0386     if ((*p)->end_vertex() && (*p)->end_vertex()->particles_out_size() == 0)
0387       edm::LogWarning("TauolappInterface::decay error: empty end vertex!");
0388   }
0389 
0390   return evt;
0391 }
0392 
0393 void TauolappInterface::statistics() { return; }
0394 
0395 void TauolappInterface::decodeMDTAU(int mdtau) {
0396   // Note-1:
0397   // I have to hack the common block directly because set<...>DecayMode(...)
0398   // only changes it in the Tauola++ instance but does NOT passes it over
0399   // to the Fortran core - this it does only one, via initialize() stuff...
0400   //
0401   // So I'll do both ways of settings, just for consistency...
0402   // but I probably need to communicate it to the Tauola(++) team...
0403   //
0404 
0405   // Note-2:
0406   // originally, the 1xx settings are meant for tau's from hard event,
0407   // and the 2xx settings are for any tau in the event record;
0408   //
0409   // later one, we'll have to take this into account...
0410   // but first I'll have to sort out what happens in the 1xx case
0411   // to tau's coming outside of hard event (if any in the record)
0412   //
0413 
0414   if (mdtau == 101 || mdtau == 201) {
0415     // override with electron mode for both tau's
0416     //
0417     Tauolapp::jaki_.jak1 = 1;
0418     Tauolapp::jaki_.jak2 = 1;
0419     Tauolapp::Tauola::setSameParticleDecayMode(1);
0420     Tauolapp::Tauola::setOppositeParticleDecayMode(1);
0421     return;
0422   }
0423 
0424   if (mdtau == 102 || mdtau == 202) {
0425     // override with muon mode for both tau's
0426     //
0427     Tauolapp::jaki_.jak1 = 2;
0428     Tauolapp::jaki_.jak2 = 2;
0429     Tauolapp::Tauola::setSameParticleDecayMode(2);
0430     Tauolapp::Tauola::setOppositeParticleDecayMode(2);
0431     return;
0432   }
0433 
0434   if (mdtau == 111 || mdtau == 211) {
0435     // override with electron mode for 1st tau
0436     // and any mode for 2nd tau
0437     //
0438     Tauolapp::jaki_.jak1 = 1;
0439     Tauolapp::jaki_.jak2 = 0;
0440     Tauolapp::Tauola::setSameParticleDecayMode(1);
0441     Tauolapp::Tauola::setOppositeParticleDecayMode(0);
0442     return;
0443   }
0444 
0445   if (mdtau == 112 || mdtau == 212) {
0446     // override with muon mode for the 1st tau
0447     // and any mode for the 2nd tau
0448     //
0449     Tauolapp::jaki_.jak1 = 2;
0450     Tauolapp::jaki_.jak2 = 0;
0451     Tauolapp::Tauola::setSameParticleDecayMode(2);
0452     Tauolapp::Tauola::setOppositeParticleDecayMode(0);
0453     return;
0454   }
0455 
0456   if (mdtau == 121 || mdtau == 221) {
0457     // override with any mode for the 1st tau
0458     // and electron mode for the 2nd tau
0459     //
0460     Tauolapp::jaki_.jak1 = 0;
0461     Tauolapp::jaki_.jak2 = 1;
0462     Tauolapp::Tauola::setSameParticleDecayMode(0);
0463     Tauolapp::Tauola::setOppositeParticleDecayMode(1);
0464     return;
0465   }
0466 
0467   if (mdtau == 122 || mdtau == 222) {
0468     // override with any mode for the 1st tau
0469     // and muon mode for the 2nd tau
0470     //
0471     Tauolapp::jaki_.jak1 = 0;
0472     Tauolapp::jaki_.jak2 = 2;
0473     Tauolapp::Tauola::setSameParticleDecayMode(0);
0474     Tauolapp::Tauola::setOppositeParticleDecayMode(2);
0475     return;
0476   }
0477 
0478   if (mdtau == 140 || mdtau == 240) {
0479     // override with pi+/- nutau mode for both tau's
0480     //
0481     Tauolapp::jaki_.jak1 = 3;
0482     Tauolapp::jaki_.jak2 = 3;
0483     Tauolapp::Tauola::setSameParticleDecayMode(3);
0484     Tauolapp::Tauola::setOppositeParticleDecayMode(3);
0485     return;
0486   }
0487 
0488   if (mdtau == 141 || mdtau == 241) {
0489     // override with pi+/- nutau mode for the 1st tau
0490     // and any mode for the 2nd tau
0491     //
0492     Tauolapp::jaki_.jak1 = 3;
0493     Tauolapp::jaki_.jak2 = 0;
0494     Tauolapp::Tauola::setSameParticleDecayMode(3);
0495     Tauolapp::Tauola::setOppositeParticleDecayMode(0);
0496     return;
0497   }
0498 
0499   if (mdtau == 142 || mdtau == 242) {
0500     // override with any mode for the 1st tau
0501     // and pi+/- nutau mode for 2nd tau
0502     //
0503     Tauolapp::jaki_.jak1 = 0;
0504     Tauolapp::jaki_.jak2 = 3;
0505     Tauolapp::Tauola::setSameParticleDecayMode(0);
0506     Tauolapp::Tauola::setOppositeParticleDecayMode(3);
0507     return;
0508   }
0509 
0510   // OK, we come here for semi-inclusive modes
0511   //
0512 
0513   // First of all, leptons and hadron modes sums
0514   //
0515   // re-scale branching ratios, just in case...
0516   //
0517   double sumBra = 0;
0518 
0519   // the number of decay modes is hardcoded at 22 because that's what it is right now in Tauola
0520   // in the future, perhaps an asscess method would be useful - communicate to Tauola team...
0521   //
0522 
0523   for (int i = 0; i < 22; i++) {
0524     sumBra += Tauolapp::taubra_.gamprt[i];
0525   }
0526   if (sumBra == 0.)
0527     return;  // perhaps need to throw ?
0528   for (int i = 0; i < 22; i++) {
0529     double newBra = Tauolapp::taubra_.gamprt[i] / sumBra;
0530     Tauolapp::Tauola::setTauBr(i + 1, newBra);
0531   }
0532   sumBra = 1.0;
0533 
0534   double sumLeptonBra = Tauolapp::taubra_.gamprt[0] + Tauolapp::taubra_.gamprt[1];
0535   double sumHadronBra = sumBra - sumLeptonBra;
0536 
0537   for (int i = 0; i < 2; i++) {
0538     fLeptonModes.push_back(i + 1);
0539     fScaledLeptonBrRatios.push_back((Tauolapp::taubra_.gamprt[i] / sumLeptonBra));
0540   }
0541   for (int i = 2; i < 22; i++) {
0542     fHadronModes.push_back(i + 1);
0543     fScaledHadronBrRatios.push_back((Tauolapp::taubra_.gamprt[i] / sumHadronBra));
0544   }
0545 
0546   fSelectDecayByEvent = true;
0547   return;
0548 }
0549 
0550 void TauolappInterface::selectDecayByMDTAU() {
0551   if (fMDTAU == 100 || fMDTAU == 200) {
0552     int mode = selectLeptonic();
0553     Tauolapp::jaki_.jak1 = mode;
0554     Tauolapp::Tauola::setSameParticleDecayMode(mode);
0555     mode = selectLeptonic();
0556     Tauolapp::jaki_.jak2 = mode;
0557     Tauolapp::Tauola::setOppositeParticleDecayMode(mode);
0558     return;
0559   }
0560 
0561   int modeL = selectLeptonic();
0562   int modeH = selectHadronic();
0563 
0564   if (fMDTAU == 110 || fMDTAU == 210) {
0565     Tauolapp::jaki_.jak1 = modeL;
0566     Tauolapp::jaki_.jak2 = 0;
0567     Tauolapp::Tauola::setSameParticleDecayMode(modeL);
0568     Tauolapp::Tauola::setOppositeParticleDecayMode(0);
0569     return;
0570   }
0571 
0572   if (fMDTAU == 120 || fMDTAU == 22) {
0573     Tauolapp::jaki_.jak1 = 0;
0574     Tauolapp::jaki_.jak2 = modeL;
0575     Tauolapp::Tauola::setSameParticleDecayMode(0);
0576     Tauolapp::Tauola::setOppositeParticleDecayMode(modeL);
0577     return;
0578   }
0579 
0580   if (fMDTAU == 114 || fMDTAU == 214) {
0581     Tauolapp::jaki_.jak1 = modeL;
0582     Tauolapp::jaki_.jak2 = modeH;
0583     Tauolapp::Tauola::setSameParticleDecayMode(modeL);
0584     Tauolapp::Tauola::setOppositeParticleDecayMode(modeH);
0585     return;
0586   }
0587 
0588   if (fMDTAU == 124 || fMDTAU == 224) {
0589     Tauolapp::jaki_.jak1 = modeH;
0590     Tauolapp::jaki_.jak2 = modeL;
0591     Tauolapp::Tauola::setSameParticleDecayMode(modeH);
0592     Tauolapp::Tauola::setOppositeParticleDecayMode(modeL);
0593     return;
0594   }
0595 
0596   if (fMDTAU == 115 || fMDTAU == 215) {
0597     Tauolapp::jaki_.jak1 = 1;
0598     Tauolapp::jaki_.jak2 = modeH;
0599     Tauolapp::Tauola::setSameParticleDecayMode(1);
0600     Tauolapp::Tauola::setOppositeParticleDecayMode(modeH);
0601     return;
0602   }
0603 
0604   if (fMDTAU == 125 || fMDTAU == 225) {
0605     Tauolapp::jaki_.jak1 = modeH;
0606     Tauolapp::jaki_.jak2 = 1;
0607     Tauolapp::Tauola::setSameParticleDecayMode(modeH);
0608     Tauolapp::Tauola::setOppositeParticleDecayMode(1);
0609     return;
0610   }
0611 
0612   if (fMDTAU == 116 || fMDTAU == 216) {
0613     Tauolapp::jaki_.jak1 = 2;
0614     Tauolapp::jaki_.jak2 = modeH;
0615     Tauolapp::Tauola::setSameParticleDecayMode(2);
0616     Tauolapp::Tauola::setOppositeParticleDecayMode(modeH);
0617     return;
0618   }
0619 
0620   if (fMDTAU == 126 || fMDTAU == 226) {
0621     Tauolapp::jaki_.jak1 = modeH;
0622     Tauolapp::jaki_.jak2 = 2;
0623     Tauolapp::Tauola::setSameParticleDecayMode(modeH);
0624     Tauolapp::Tauola::setOppositeParticleDecayMode(2);
0625     return;
0626   }
0627 
0628   if (fMDTAU == 130 || fMDTAU == 230) {
0629     Tauolapp::jaki_.jak1 = modeH;
0630     Tauolapp::jaki_.jak2 = selectHadronic();
0631     Tauolapp::Tauola::setSameParticleDecayMode(modeH);
0632     Tauolapp::Tauola::setOppositeParticleDecayMode(Tauolapp::jaki_.jak2);
0633     return;
0634   }
0635 
0636   if (fMDTAU == 131 || fMDTAU == 231) {
0637     Tauolapp::jaki_.jak1 = modeH;
0638     Tauolapp::jaki_.jak2 = 0;
0639     Tauolapp::Tauola::setSameParticleDecayMode(modeH);
0640     Tauolapp::Tauola::setOppositeParticleDecayMode(0);
0641     return;
0642   }
0643 
0644   if (fMDTAU == 132 || fMDTAU == 232) {
0645     Tauolapp::jaki_.jak1 = 0;
0646     Tauolapp::jaki_.jak2 = modeH;
0647     Tauolapp::Tauola::setSameParticleDecayMode(0);
0648     Tauolapp::Tauola::setOppositeParticleDecayMode(modeH);
0649     return;
0650   }
0651 
0652   // unlikely that we get here on unknown mdtau
0653   // - there's a protection earlier
0654   // but if we do, just set defaults
0655   // probably need to spit a warning...
0656   //
0657   Tauolapp::Tauola::setSameParticleDecayMode(0);
0658   Tauolapp::Tauola::setOppositeParticleDecayMode(0);
0659 
0660   return;
0661 }
0662 
0663 int TauolappInterface::selectLeptonic() {
0664   float prob = flat();
0665 
0666   if (prob > 0. && prob <= fScaledLeptonBrRatios[0]) {
0667     return 1;
0668   } else if (prob > fScaledLeptonBrRatios[1] && prob <= 1.) {
0669     return 2;
0670   }
0671 
0672   return 0;
0673 }
0674 
0675 int TauolappInterface::selectHadronic() {
0676   float prob = 0.;
0677   int len = 1;
0678   ranmar_(&prob, &len);
0679 
0680   double sumBra = fScaledHadronBrRatios[0];
0681   if (prob > 0. && prob <= sumBra) {
0682     return fHadronModes[0];
0683   } else {
0684     int NN = fScaledHadronBrRatios.size();
0685     for (int i = 1; i < NN; i++) {
0686       if (prob > sumBra && prob <= (sumBra + fScaledHadronBrRatios[i])) {
0687         return fHadronModes[i];
0688       }
0689       sumBra += fScaledHadronBrRatios[i];
0690     }
0691   }
0692 
0693   return 0;
0694 }
0695 
0696 HepMC::GenEvent* TauolappInterface::make_simple_tau_event(const TLorentzVector& l, int pdgid, int status) {
0697   HepMC::GenEvent* event = new HepMC::GenEvent();
0698   // make tau's four vector
0699   HepMC::FourVector momentum_tau1(l.Px(), l.Py(), l.Pz(), l.E());
0700   // make particles
0701   HepMC::GenParticle* tau1 = new HepMC::GenParticle(momentum_tau1, pdgid, status);
0702   // make the vertex
0703   HepMC::GenVertex* vertex = new HepMC::GenVertex();
0704   vertex->add_particle_out(tau1);
0705   event->add_vertex(vertex);
0706   return event;
0707 }
0708 
0709 void TauolappInterface::update_particles(HepMC::GenParticle* partHep,
0710                                          HepMC::GenEvent* theEvent,
0711                                          HepMC::GenParticle* p,
0712                                          TVector3& boost) {
0713   partHep->set_status(p->status());
0714   if (p->end_vertex()) {
0715     if (!partHep->end_vertex()) {
0716       HepMC::GenVertex* vtx = new HepMC::GenVertex(p->end_vertex()->position());
0717       theEvent->add_vertex(vtx);
0718       vtx->add_particle_in(partHep);
0719     }
0720     if (p->end_vertex()->particles_out_size() != 0) {
0721       for (HepMC::GenVertex::particles_out_const_iterator d = p->end_vertex()->particles_out_const_begin();
0722            d != p->end_vertex()->particles_out_const_end();
0723            d++) {
0724         // Create daughter and add to event
0725         TLorentzVector l((*d)->momentum().px(), (*d)->momentum().py(), (*d)->momentum().pz(), (*d)->momentum().e());
0726         l.Boost(boost);
0727         HepMC::FourVector momentum(l.Px(), l.Py(), l.Pz(), l.E());
0728         HepMC::GenParticle* daughter = new HepMC::GenParticle(momentum, (*d)->pdg_id(), (*d)->status());
0729         daughter->suggest_barcode(theEvent->particles_size() + 1);
0730         partHep->end_vertex()->add_particle_out(daughter);
0731         if ((*d)->end_vertex())
0732           update_particles(daughter, theEvent, (*d), boost);
0733       }
0734     }
0735   }
0736 }
0737 
0738 bool TauolappInterface::isLastTauInChain(const HepMC::GenParticle* tau) {
0739   if (tau->end_vertex()) {
0740     HepMC::GenVertex::particle_iterator dau;
0741     for (dau = tau->end_vertex()->particles_begin(HepMC::children);
0742          dau != tau->end_vertex()->particles_end(HepMC::children);
0743          dau++) {
0744       int dau_pid = (*dau)->pdg_id();
0745       if (dau_pid == tau->pdg_id())
0746         return false;
0747     }
0748   }
0749   return true;
0750 }
0751 
0752 double TauolappInterface::MatchedLHESpinUp(HepMC::GenParticle* tau,
0753                                            std::vector<HepMC::GenParticle>& p,
0754                                            std::vector<double>& spinup,
0755                                            std::vector<int>& m_idx) {
0756   HepMC::GenParticle* Tau = FirstTauInChain(tau);
0757   HepMC::GenParticle* mother = GetMother(Tau);
0758   TLorentzVector t(tau->momentum().px(), tau->momentum().py(), tau->momentum().pz(), tau->momentum().e());
0759   TLorentzVector m(mother->momentum().px(), mother->momentum().py(), mother->momentum().pz(), mother->momentum().e());
0760   for (unsigned int i = 0; i < p.size(); i++) {
0761     if (tau->pdg_id() == p.at(i).pdg_id()) {
0762       if (mother->pdg_id() == p.at(m_idx.at(i)).pdg_id()) {
0763         TLorentzVector pm(p.at(m_idx.at(i)).momentum().px(),
0764                           p.at(m_idx.at(i)).momentum().py(),
0765                           p.at(m_idx.at(i)).momentum().pz(),
0766                           p.at(m_idx.at(i)).momentum().e());
0767         if (fabs(m.M() - pm.M()) < dmMatch)
0768           return spinup.at(i);
0769       }
0770     }
0771   }
0772   return 0;
0773 }
0774 
0775 HepMC::GenParticle* TauolappInterface::FirstTauInChain(HepMC::GenParticle* tau) {
0776   if (tau->production_vertex()) {
0777     HepMC::GenVertex::particle_iterator mother;
0778     for (mother = tau->production_vertex()->particles_begin(HepMC::parents);
0779          mother != tau->production_vertex()->particles_end(HepMC::parents);
0780          mother++) {
0781       if ((*mother)->pdg_id() == tau->pdg_id())
0782         return FirstTauInChain(*mother);  // recursive call to get mother with different pdgid
0783     }
0784   }
0785   return tau;
0786 }
0787 
0788 HepMC::GenParticle* TauolappInterface::GetMother(HepMC::GenParticle* tau) {
0789   if (tau->production_vertex()) {
0790     HepMC::GenVertex::particle_iterator mother;
0791     for (mother = tau->production_vertex()->particles_begin(HepMC::parents);
0792          mother != tau->production_vertex()->particles_end(HepMC::parents);
0793          mother++) {
0794       if ((*mother)->pdg_id() == tau->pdg_id())
0795         return GetMother(*mother);  // recursive call to get mother with different pdgid
0796       return (*mother);
0797     }
0798   }
0799   return tau;
0800 }
0801 
0802 void TauolappInterface::BoostProdToLabLifeTimeInDecays(HepMC::GenParticle* p,
0803                                                        TLorentzVector& lab,
0804                                                        TLorentzVector& prod) {
0805   if (p->end_vertex() && p->production_vertex()) {
0806     HepMC::GenVertex* PGenVtx = p->production_vertex();
0807     HepMC::GenVertex* EGenVtx = p->end_vertex();
0808     double VxDec = PGenVtx->position().x() + lab.Px() / prod.Px() * (EGenVtx->position().x() - PGenVtx->position().x());
0809     double VyDec = PGenVtx->position().y() + lab.Py() / prod.Py() * (EGenVtx->position().y() - PGenVtx->position().y());
0810     double VzDec = PGenVtx->position().z() + lab.Pz() / prod.Pz() * (EGenVtx->position().z() - PGenVtx->position().z());
0811     double VtDec = PGenVtx->position().t() + lab.Pt() / prod.Pt() * (EGenVtx->position().t() - PGenVtx->position().t());
0812     EGenVtx->set_position(HepMC::FourVector(VxDec, VyDec, VzDec, VtDec));
0813     for (HepMC::GenVertex::particle_iterator dau = p->end_vertex()->particles_begin(HepMC::children);
0814          dau != p->end_vertex()->particles_end(HepMC::children);
0815          dau++) {
0816       BoostProdToLabLifeTimeInDecays((*dau), lab, prod);  //recursively modify everything in the decay chain
0817     }
0818   }
0819 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!