Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-12-13 23:50:10 

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   Tauolapp::Log::IgnoreRedirection(true);
0210 
0211   return;
0212 }
0213 
0214 double TauolappInterface::flat() {
0215   if (!fRandomEngine) {
0216     throw cms::Exception("LogicError")
0217         << "TauolaInterface::flat: Attempt to generate random number when engine pointer is null\n"
0218         << "This might mean that the code was modified to generate a random number outside the\n"
0219         << "event and beginLuminosityBlock methods, which is not allowed.\n";
0220   }
0221   return fRandomEngine->flat();
0222 }
0223 
0224 HepMC::GenEvent* TauolappInterface::decay(HepMC::GenEvent* evt) {
0225   if (!fIsInitialized)
0226     return evt;
0227   Tauolapp::Tauola::setRandomGenerator(
0228       gen::TauolappInterface::flat);  // rest tauola++ random number incase other modules use tauola++
0229   int NPartBefore = evt->particles_size();
0230   int NVtxBefore = evt->vertices_size();
0231 
0232   // what do we do if Hep::GenEvent size is larger than 10K ???
0233   // Tauola (& Photos, BTW) can only handle up to 10K via HEPEVT,
0234   // and in case of CMS, it's only up to 4K !!!
0235   // override decay mode if needs be
0236 
0237   if (fSelectDecayByEvent) {
0238     selectDecayByMDTAU();
0239   }
0240   if (dolhe && lhe != nullptr) {
0241     std::vector<HepMC::GenParticle> particles;
0242     std::vector<int> m_idx;
0243     std::vector<double> spinup = lhe->getHEPEUP()->SPINUP;
0244     std::vector<int> pdg = lhe->getHEPEUP()->IDUP;
0245     for (unsigned int i = 0; i < spinup.size(); i++) {
0246       particles.push_back(HepMC::GenParticle(HepMC::FourVector(lhe->getHEPEUP()->PUP.at(i)[0],
0247                                                                lhe->getHEPEUP()->PUP.at(i)[1],
0248                                                                lhe->getHEPEUP()->PUP.at(i)[2],
0249                                                                lhe->getHEPEUP()->PUP.at(i)[3]),
0250                                              lhe->getHEPEUP()->IDUP.at(i)));
0251       int status = lhe->getHEPEUP()->ISTUP.at(i);
0252       particles.at(particles.size() - 1).set_generated_mass(lhe->getHEPEUP()->PUP.at(i)[4]);
0253       particles.at(particles.size() - 1).set_status(status > 0 ? (status == 2 ? 3 : status) : 3);
0254       m_idx.push_back(lhe->getHEPEUP()->MOTHUP.at(i).first - 1);  // correct for fortran index offset
0255     }
0256     // match to taus in hepmc and identify mother of taus
0257     bool hastaus(false);
0258     std::vector<HepMC::GenParticle*> match;
0259     for (HepMC::GenEvent::particle_const_iterator iter = evt->particles_begin(); iter != evt->particles_end(); iter++) {
0260       if (abs((*iter)->pdg_id()) == 15) {
0261         hastaus = true;
0262         int mother_pid(0);
0263         // check imediate parent to avoid parent tau ie tau->taugamma
0264         for (HepMC::GenVertex::particle_iterator mother = (*iter)->production_vertex()->particles_begin(HepMC::parents);
0265              mother != (*iter)->production_vertex()->particles_end(HepMC::parents);
0266              mother++) {
0267           mother_pid = (*mother)->pdg_id();
0268           if (mother_pid != (*iter)->pdg_id()) {
0269             // match against lhe record
0270             if (abs(mother_pid) == 24 ||  // W
0271                 abs(mother_pid) == 37 ||  // H+/-
0272                 abs(mother_pid) == 23 ||  // Z
0273                 abs(mother_pid) == 22 ||  // gamma
0274                 abs(mother_pid) == 25 ||  // H0 SM
0275                 abs(mother_pid) == 35 ||  // H0
0276                 abs(mother_pid) == 36     // A0
0277             ) {
0278               bool isfound = false;
0279               for (unsigned int k = 0; k < match.size(); k++) {
0280                 if ((*mother) == match.at(k))
0281                   isfound = true;
0282               }
0283               if (!isfound)
0284                 match.push_back(*mother);
0285             }
0286           }
0287         }
0288       }
0289     }
0290     if (hastaus) {
0291       // if is single gauge boson decay and match helicities
0292       if (match.size() == 1 && dolheBosonCorr) {
0293         for (int i = 0; i < ntries; i++) {
0294           // re-decay taus then check if helicities match
0295           auto* t_event = new Tauolapp::TauolaHepMCEvent(evt);
0296           t_event->undecayTaus();
0297           t_event->decayTaus();
0298           bool ismatch = true;
0299           for (unsigned int j = 0; j < spinup.size(); j++) {
0300             if (abs(pdg.at(j)) == 15) {
0301               double diffhelminus = (-1.0 * (double)Tauolapp::Tauola::getHelMinus() -
0302                                      spinup.at(j));  // -1.0 to correct for tauola feature
0303               double diffhelplus = ((double)Tauolapp::Tauola::getHelPlus() - spinup.at(j));
0304               if (pdg.at(j) == 15 && diffhelminus > 0.5)
0305                 ismatch = false;
0306               if (pdg.at(j) == -15 && diffhelplus > 0.5)
0307                 ismatch = false;
0308             }
0309           }
0310           delete t_event;
0311           if (ismatch)
0312             break;
0313         }
0314       } else {
0315         // If the event does not contain a single gauge boson the code will be run with
0316         // remove all tau decays
0317         auto* t_event = new Tauolapp::TauolaHepMCEvent(evt);
0318         t_event->undecayTaus();
0319         delete t_event;
0320         // decay all taus manually based on the helicity
0321         for (HepMC::GenEvent::particle_const_iterator iter = evt->particles_begin(); iter != evt->particles_end();
0322              iter++) {
0323           if (abs((*iter)->pdg_id()) == 15 && isLastTauInChain(*iter)) {
0324             TLorentzVector ltau(
0325                 (*iter)->momentum().px(), (*iter)->momentum().py(), (*iter)->momentum().pz(), (*iter)->momentum().e());
0326             HepMC::GenParticle* m = GetMother(*iter);
0327             TLorentzVector mother(m->momentum().px(), m->momentum().py(), m->momentum().pz(), m->momentum().e());
0328             TVector3 boost = -1.0 * mother.BoostVector();  // boost into mother's CM frame
0329             TLorentzVector ltau_lab = ltau;
0330             ltau.Boost(boost);
0331             mother.Boost(boost);
0332             HepMC::GenEvent* tauevt = make_simple_tau_event(ltau, (*iter)->pdg_id(), (*iter)->status());
0333             HepMC::GenParticle* p = (*(tauevt->particles_begin()));
0334             Tauolapp::TauolaParticle* tp = new Tauolapp::TauolaHepMCParticle(p);
0335             double helicity = MatchedLHESpinUp(*iter, particles, spinup, m_idx);  // get helicity from lhe
0336             if ((*iter)->pdg_id() == 15)
0337               helicity *= -1.0;
0338             tp->undecay();
0339             // use |S_{tau}|=0.999999 to avoid issues with numerical roundoff
0340             Tauolapp::Tauola::decayOne(tp, true, 0, 0, ((double)helicity) * 0.999999);
0341             boost *= -1.0;  // boost back to lab frame
0342             mother.Boost(boost);
0343             update_particles((*iter), evt, p, boost);
0344             //correct tau liftetime for boost (change rest frame from mothers to taus)
0345             BoostProdToLabLifeTimeInDecays((*iter), ltau_lab, ltau);
0346             delete tauevt;
0347           }
0348         }
0349       }
0350     }
0351   } else {
0352     //construct tmp TAUOLA event
0353     auto* t_event = new Tauolapp::TauolaHepMCEvent(evt);
0354     //t_event->undecayTaus();
0355     t_event->decayTaus();
0356     delete t_event;
0357   }
0358 
0359   for (int iv = NVtxBefore + 1; iv <= evt->vertices_size(); iv++) {
0360     HepMC::GenVertex* GenVtx = evt->barcode_to_vertex(-iv);
0361     //
0362     // now find decay products with funky barcode, weed out and replace with clones of sensible barcode
0363     // we can NOT change the barcode while iterating, because iterators do depend on the barcoding
0364     // thus we have to take a 2-step procedure
0365     //
0366     std::vector<int> BCodes;
0367     BCodes.clear();
0368     for (HepMC::GenVertex::particle_iterator pitr = GenVtx->particles_begin(HepMC::children);
0369          pitr != GenVtx->particles_end(HepMC::children);
0370          ++pitr) {
0371       if ((*pitr)->barcode() > 10000) {
0372         BCodes.push_back((*pitr)->barcode());
0373       }
0374     }
0375     if (!BCodes.empty()) {
0376       for (size_t ibc = 0; ibc < BCodes.size(); ibc++) {
0377         HepMC::GenParticle* p1 = evt->barcode_to_particle(BCodes[ibc]);
0378         int nbc = p1->barcode() - 10000 + NPartBefore;
0379         p1->suggest_barcode(nbc);
0380       }
0381     }
0382   }
0383 
0384   for (HepMC::GenEvent::particle_const_iterator p = evt->particles_begin(); p != evt->particles_end(); ++p) {
0385     if ((*p)->end_vertex() && (*p)->status() == 1)
0386       (*p)->set_status(2);
0387     if ((*p)->end_vertex() && (*p)->end_vertex()->particles_out_size() == 0)
0388       edm::LogWarning("TauolappInterface::decay error: empty end vertex!");
0389   }
0390 
0391   return evt;
0392 }
0393 
0394 void TauolappInterface::statistics() { return; }
0395 
0396 void TauolappInterface::decodeMDTAU(int mdtau) {
0397   // Note-1:
0398   // I have to hack the common block directly because set<...>DecayMode(...)
0399   // only changes it in the Tauola++ instance but does NOT passes it over
0400   // to the Fortran core - this it does only one, via initialize() stuff...
0401   //
0402   // So I'll do both ways of settings, just for consistency...
0403   // but I probably need to communicate it to the Tauola(++) team...
0404   //
0405 
0406   // Note-2:
0407   // originally, the 1xx settings are meant for tau's from hard event,
0408   // and the 2xx settings are for any tau in the event record;
0409   //
0410   // later one, we'll have to take this into account...
0411   // but first I'll have to sort out what happens in the 1xx case
0412   // to tau's coming outside of hard event (if any in the record)
0413   //
0414 
0415   if (mdtau == 101 || mdtau == 201) {
0416     // override with electron mode for both tau's
0417     //
0418     Tauolapp::jaki_.jak1 = 1;
0419     Tauolapp::jaki_.jak2 = 1;
0420     Tauolapp::Tauola::setSameParticleDecayMode(1);
0421     Tauolapp::Tauola::setOppositeParticleDecayMode(1);
0422     return;
0423   }
0424 
0425   if (mdtau == 102 || mdtau == 202) {
0426     // override with muon mode for both tau's
0427     //
0428     Tauolapp::jaki_.jak1 = 2;
0429     Tauolapp::jaki_.jak2 = 2;
0430     Tauolapp::Tauola::setSameParticleDecayMode(2);
0431     Tauolapp::Tauola::setOppositeParticleDecayMode(2);
0432     return;
0433   }
0434 
0435   if (mdtau == 111 || mdtau == 211) {
0436     // override with electron mode for 1st tau
0437     // and any mode for 2nd tau
0438     //
0439     Tauolapp::jaki_.jak1 = 1;
0440     Tauolapp::jaki_.jak2 = 0;
0441     Tauolapp::Tauola::setSameParticleDecayMode(1);
0442     Tauolapp::Tauola::setOppositeParticleDecayMode(0);
0443     return;
0444   }
0445 
0446   if (mdtau == 112 || mdtau == 212) {
0447     // override with muon mode for the 1st tau
0448     // and any mode for the 2nd tau
0449     //
0450     Tauolapp::jaki_.jak1 = 2;
0451     Tauolapp::jaki_.jak2 = 0;
0452     Tauolapp::Tauola::setSameParticleDecayMode(2);
0453     Tauolapp::Tauola::setOppositeParticleDecayMode(0);
0454     return;
0455   }
0456 
0457   if (mdtau == 121 || mdtau == 221) {
0458     // override with any mode for the 1st tau
0459     // and electron mode for the 2nd tau
0460     //
0461     Tauolapp::jaki_.jak1 = 0;
0462     Tauolapp::jaki_.jak2 = 1;
0463     Tauolapp::Tauola::setSameParticleDecayMode(0);
0464     Tauolapp::Tauola::setOppositeParticleDecayMode(1);
0465     return;
0466   }
0467 
0468   if (mdtau == 122 || mdtau == 222) {
0469     // override with any mode for the 1st tau
0470     // and muon mode for the 2nd tau
0471     //
0472     Tauolapp::jaki_.jak1 = 0;
0473     Tauolapp::jaki_.jak2 = 2;
0474     Tauolapp::Tauola::setSameParticleDecayMode(0);
0475     Tauolapp::Tauola::setOppositeParticleDecayMode(2);
0476     return;
0477   }
0478 
0479   if (mdtau == 140 || mdtau == 240) {
0480     // override with pi+/- nutau mode for both tau's
0481     //
0482     Tauolapp::jaki_.jak1 = 3;
0483     Tauolapp::jaki_.jak2 = 3;
0484     Tauolapp::Tauola::setSameParticleDecayMode(3);
0485     Tauolapp::Tauola::setOppositeParticleDecayMode(3);
0486     return;
0487   }
0488 
0489   if (mdtau == 141 || mdtau == 241) {
0490     // override with pi+/- nutau mode for the 1st tau
0491     // and any mode for the 2nd tau
0492     //
0493     Tauolapp::jaki_.jak1 = 3;
0494     Tauolapp::jaki_.jak2 = 0;
0495     Tauolapp::Tauola::setSameParticleDecayMode(3);
0496     Tauolapp::Tauola::setOppositeParticleDecayMode(0);
0497     return;
0498   }
0499 
0500   if (mdtau == 142 || mdtau == 242) {
0501     // override with any mode for the 1st tau
0502     // and pi+/- nutau mode for 2nd tau
0503     //
0504     Tauolapp::jaki_.jak1 = 0;
0505     Tauolapp::jaki_.jak2 = 3;
0506     Tauolapp::Tauola::setSameParticleDecayMode(0);
0507     Tauolapp::Tauola::setOppositeParticleDecayMode(3);
0508     return;
0509   }
0510 
0511   // OK, we come here for semi-inclusive modes
0512   //
0513 
0514   // First of all, leptons and hadron modes sums
0515   //
0516   // re-scale branching ratios, just in case...
0517   //
0518   double sumBra = 0;
0519 
0520   // the number of decay modes is hardcoded at 22 because that's what it is right now in Tauola
0521   // in the future, perhaps an asscess method would be useful - communicate to Tauola team...
0522   //
0523 
0524   for (int i = 0; i < 22; i++) {
0525     sumBra += Tauolapp::taubra_.gamprt[i];
0526   }
0527   if (sumBra == 0.)
0528     return;  // perhaps need to throw ?
0529   for (int i = 0; i < 22; i++) {
0530     double newBra = Tauolapp::taubra_.gamprt[i] / sumBra;
0531     Tauolapp::Tauola::setTauBr(i + 1, newBra);
0532   }
0533   sumBra = 1.0;
0534 
0535   double sumLeptonBra = Tauolapp::taubra_.gamprt[0] + Tauolapp::taubra_.gamprt[1];
0536   double sumHadronBra = sumBra - sumLeptonBra;
0537 
0538   for (int i = 0; i < 2; i++) {
0539     fLeptonModes.push_back(i + 1);
0540     fScaledLeptonBrRatios.push_back((Tauolapp::taubra_.gamprt[i] / sumLeptonBra));
0541   }
0542   for (int i = 2; i < 22; i++) {
0543     fHadronModes.push_back(i + 1);
0544     fScaledHadronBrRatios.push_back((Tauolapp::taubra_.gamprt[i] / sumHadronBra));
0545   }
0546 
0547   fSelectDecayByEvent = true;
0548   return;
0549 }
0550 
0551 void TauolappInterface::selectDecayByMDTAU() {
0552   if (fMDTAU == 100 || fMDTAU == 200) {
0553     int mode = selectLeptonic();
0554     Tauolapp::jaki_.jak1 = mode;
0555     Tauolapp::Tauola::setSameParticleDecayMode(mode);
0556     mode = selectLeptonic();
0557     Tauolapp::jaki_.jak2 = mode;
0558     Tauolapp::Tauola::setOppositeParticleDecayMode(mode);
0559     return;
0560   }
0561 
0562   int modeL = selectLeptonic();
0563   int modeH = selectHadronic();
0564 
0565   if (fMDTAU == 110 || fMDTAU == 210) {
0566     Tauolapp::jaki_.jak1 = modeL;
0567     Tauolapp::jaki_.jak2 = 0;
0568     Tauolapp::Tauola::setSameParticleDecayMode(modeL);
0569     Tauolapp::Tauola::setOppositeParticleDecayMode(0);
0570     return;
0571   }
0572 
0573   if (fMDTAU == 120 || fMDTAU == 22) {
0574     Tauolapp::jaki_.jak1 = 0;
0575     Tauolapp::jaki_.jak2 = modeL;
0576     Tauolapp::Tauola::setSameParticleDecayMode(0);
0577     Tauolapp::Tauola::setOppositeParticleDecayMode(modeL);
0578     return;
0579   }
0580 
0581   if (fMDTAU == 114 || fMDTAU == 214) {
0582     Tauolapp::jaki_.jak1 = modeL;
0583     Tauolapp::jaki_.jak2 = modeH;
0584     Tauolapp::Tauola::setSameParticleDecayMode(modeL);
0585     Tauolapp::Tauola::setOppositeParticleDecayMode(modeH);
0586     return;
0587   }
0588 
0589   if (fMDTAU == 124 || fMDTAU == 224) {
0590     Tauolapp::jaki_.jak1 = modeH;
0591     Tauolapp::jaki_.jak2 = modeL;
0592     Tauolapp::Tauola::setSameParticleDecayMode(modeH);
0593     Tauolapp::Tauola::setOppositeParticleDecayMode(modeL);
0594     return;
0595   }
0596 
0597   if (fMDTAU == 115 || fMDTAU == 215) {
0598     Tauolapp::jaki_.jak1 = 1;
0599     Tauolapp::jaki_.jak2 = modeH;
0600     Tauolapp::Tauola::setSameParticleDecayMode(1);
0601     Tauolapp::Tauola::setOppositeParticleDecayMode(modeH);
0602     return;
0603   }
0604 
0605   if (fMDTAU == 125 || fMDTAU == 225) {
0606     Tauolapp::jaki_.jak1 = modeH;
0607     Tauolapp::jaki_.jak2 = 1;
0608     Tauolapp::Tauola::setSameParticleDecayMode(modeH);
0609     Tauolapp::Tauola::setOppositeParticleDecayMode(1);
0610     return;
0611   }
0612 
0613   if (fMDTAU == 116 || fMDTAU == 216) {
0614     Tauolapp::jaki_.jak1 = 2;
0615     Tauolapp::jaki_.jak2 = modeH;
0616     Tauolapp::Tauola::setSameParticleDecayMode(2);
0617     Tauolapp::Tauola::setOppositeParticleDecayMode(modeH);
0618     return;
0619   }
0620 
0621   if (fMDTAU == 126 || fMDTAU == 226) {
0622     Tauolapp::jaki_.jak1 = modeH;
0623     Tauolapp::jaki_.jak2 = 2;
0624     Tauolapp::Tauola::setSameParticleDecayMode(modeH);
0625     Tauolapp::Tauola::setOppositeParticleDecayMode(2);
0626     return;
0627   }
0628 
0629   if (fMDTAU == 130 || fMDTAU == 230) {
0630     Tauolapp::jaki_.jak1 = modeH;
0631     Tauolapp::jaki_.jak2 = selectHadronic();
0632     Tauolapp::Tauola::setSameParticleDecayMode(modeH);
0633     Tauolapp::Tauola::setOppositeParticleDecayMode(Tauolapp::jaki_.jak2);
0634     return;
0635   }
0636 
0637   if (fMDTAU == 131 || fMDTAU == 231) {
0638     Tauolapp::jaki_.jak1 = modeH;
0639     Tauolapp::jaki_.jak2 = 0;
0640     Tauolapp::Tauola::setSameParticleDecayMode(modeH);
0641     Tauolapp::Tauola::setOppositeParticleDecayMode(0);
0642     return;
0643   }
0644 
0645   if (fMDTAU == 132 || fMDTAU == 232) {
0646     Tauolapp::jaki_.jak1 = 0;
0647     Tauolapp::jaki_.jak2 = modeH;
0648     Tauolapp::Tauola::setSameParticleDecayMode(0);
0649     Tauolapp::Tauola::setOppositeParticleDecayMode(modeH);
0650     return;
0651   }
0652 
0653   // unlikely that we get here on unknown mdtau
0654   // - there's a protection earlier
0655   // but if we do, just set defaults
0656   // probably need to spit a warning...
0657   //
0658   Tauolapp::Tauola::setSameParticleDecayMode(0);
0659   Tauolapp::Tauola::setOppositeParticleDecayMode(0);
0660 
0661   return;
0662 }
0663 
0664 int TauolappInterface::selectLeptonic() {
0665   float prob = flat();
0666 
0667   if (prob > 0. && prob <= fScaledLeptonBrRatios[0]) {
0668     return 1;
0669   } else if (prob > fScaledLeptonBrRatios[1] && prob <= 1.) {
0670     return 2;
0671   }
0672 
0673   return 0;
0674 }
0675 
0676 int TauolappInterface::selectHadronic() {
0677   float prob = 0.;
0678   int len = 1;
0679   ranmar_(&prob, &len);
0680 
0681   double sumBra = fScaledHadronBrRatios[0];
0682   if (prob > 0. && prob <= sumBra) {
0683     return fHadronModes[0];
0684   } else {
0685     int NN = fScaledHadronBrRatios.size();
0686     for (int i = 1; i < NN; i++) {
0687       if (prob > sumBra && prob <= (sumBra + fScaledHadronBrRatios[i])) {
0688         return fHadronModes[i];
0689       }
0690       sumBra += fScaledHadronBrRatios[i];
0691     }
0692   }
0693 
0694   return 0;
0695 }
0696 
0697 HepMC::GenEvent* TauolappInterface::make_simple_tau_event(const TLorentzVector& l, int pdgid, int status) {
0698   HepMC::GenEvent* event = new HepMC::GenEvent();
0699   // make tau's four vector
0700   HepMC::FourVector momentum_tau1(l.Px(), l.Py(), l.Pz(), l.E());
0701   // make particles
0702   HepMC::GenParticle* tau1 = new HepMC::GenParticle(momentum_tau1, pdgid, status);
0703   // make the vertex
0704   HepMC::GenVertex* vertex = new HepMC::GenVertex();
0705   vertex->add_particle_out(tau1);
0706   event->add_vertex(vertex);
0707   return event;
0708 }
0709 
0710 void TauolappInterface::update_particles(HepMC::GenParticle* partHep,
0711                                          HepMC::GenEvent* theEvent,
0712                                          HepMC::GenParticle* p,
0713                                          TVector3& boost) {
0714   partHep->set_status(p->status());
0715   if (p->end_vertex()) {
0716     if (!partHep->end_vertex()) {
0717       HepMC::GenVertex* vtx = new HepMC::GenVertex(p->end_vertex()->position());
0718       theEvent->add_vertex(vtx);
0719       vtx->add_particle_in(partHep);
0720     }
0721     if (p->end_vertex()->particles_out_size() != 0) {
0722       for (HepMC::GenVertex::particles_out_const_iterator d = p->end_vertex()->particles_out_const_begin();
0723            d != p->end_vertex()->particles_out_const_end();
0724            d++) {
0725         // Create daughter and add to event
0726         TLorentzVector l((*d)->momentum().px(), (*d)->momentum().py(), (*d)->momentum().pz(), (*d)->momentum().e());
0727         l.Boost(boost);
0728         HepMC::FourVector momentum(l.Px(), l.Py(), l.Pz(), l.E());
0729         HepMC::GenParticle* daughter = new HepMC::GenParticle(momentum, (*d)->pdg_id(), (*d)->status());
0730         daughter->suggest_barcode(theEvent->particles_size() + 1);
0731         partHep->end_vertex()->add_particle_out(daughter);
0732         if ((*d)->end_vertex())
0733           update_particles(daughter, theEvent, (*d), boost);
0734       }
0735     }
0736   }
0737 }
0738 
0739 bool TauolappInterface::isLastTauInChain(const HepMC::GenParticle* tau) {
0740   if (tau->end_vertex()) {
0741     HepMC::GenVertex::particle_iterator dau;
0742     for (dau = tau->end_vertex()->particles_begin(HepMC::children);
0743          dau != tau->end_vertex()->particles_end(HepMC::children);
0744          dau++) {
0745       int dau_pid = (*dau)->pdg_id();
0746       if (dau_pid == tau->pdg_id())
0747         return false;
0748     }
0749   }
0750   return true;
0751 }
0752 
0753 double TauolappInterface::MatchedLHESpinUp(HepMC::GenParticle* tau,
0754                                            std::vector<HepMC::GenParticle>& p,
0755                                            std::vector<double>& spinup,
0756                                            std::vector<int>& m_idx) {
0757   HepMC::GenParticle* Tau = FirstTauInChain(tau);
0758   HepMC::GenParticle* mother = GetMother(Tau);
0759   TLorentzVector t(tau->momentum().px(), tau->momentum().py(), tau->momentum().pz(), tau->momentum().e());
0760   TLorentzVector m(mother->momentum().px(), mother->momentum().py(), mother->momentum().pz(), mother->momentum().e());
0761   for (unsigned int i = 0; i < p.size(); i++) {
0762     if (tau->pdg_id() == p.at(i).pdg_id()) {
0763       if (mother->pdg_id() == p.at(m_idx.at(i)).pdg_id()) {
0764         TLorentzVector pm(p.at(m_idx.at(i)).momentum().px(),
0765                           p.at(m_idx.at(i)).momentum().py(),
0766                           p.at(m_idx.at(i)).momentum().pz(),
0767                           p.at(m_idx.at(i)).momentum().e());
0768         if (fabs(m.M() - pm.M()) < dmMatch)
0769           return spinup.at(i);
0770       }
0771     }
0772   }
0773   return 0;
0774 }
0775 
0776 HepMC::GenParticle* TauolappInterface::FirstTauInChain(HepMC::GenParticle* tau) {
0777   if (tau->production_vertex()) {
0778     HepMC::GenVertex::particle_iterator mother;
0779     for (mother = tau->production_vertex()->particles_begin(HepMC::parents);
0780          mother != tau->production_vertex()->particles_end(HepMC::parents);
0781          mother++) {
0782       if ((*mother)->pdg_id() == tau->pdg_id())
0783         return FirstTauInChain(*mother);  // recursive call to get mother with different pdgid
0784     }
0785   }
0786   return tau;
0787 }
0788 
0789 HepMC::GenParticle* TauolappInterface::GetMother(HepMC::GenParticle* tau) {
0790   if (tau->production_vertex()) {
0791     HepMC::GenVertex::particle_iterator mother;
0792     for (mother = tau->production_vertex()->particles_begin(HepMC::parents);
0793          mother != tau->production_vertex()->particles_end(HepMC::parents);
0794          mother++) {
0795       if ((*mother)->pdg_id() == tau->pdg_id())
0796         return GetMother(*mother);  // recursive call to get mother with different pdgid
0797       return (*mother);
0798     }
0799   }
0800   return tau;
0801 }
0802 
0803 void TauolappInterface::BoostProdToLabLifeTimeInDecays(HepMC::GenParticle* p,
0804                                                        TLorentzVector& lab,
0805                                                        TLorentzVector& prod) {
0806   if (p->end_vertex() && p->production_vertex()) {
0807     HepMC::GenVertex* PGenVtx = p->production_vertex();
0808     HepMC::GenVertex* EGenVtx = p->end_vertex();
0809     double VxDec = PGenVtx->position().x() + lab.Px() / prod.Px() * (EGenVtx->position().x() - PGenVtx->position().x());
0810     double VyDec = PGenVtx->position().y() + lab.Py() / prod.Py() * (EGenVtx->position().y() - PGenVtx->position().y());
0811     double VzDec = PGenVtx->position().z() + lab.Pz() / prod.Pz() * (EGenVtx->position().z() - PGenVtx->position().z());
0812     double VtDec = PGenVtx->position().t() + lab.Pt() / prod.Pt() * (EGenVtx->position().t() - PGenVtx->position().t());
0813     EGenVtx->set_position(HepMC::FourVector(VxDec, VyDec, VzDec, VtDec));
0814     for (HepMC::GenVertex::particle_iterator dau = p->end_vertex()->particles_begin(HepMC::children);
0815          dau != p->end_vertex()->particles_end(HepMC::children);
0816          dau++) {
0817       BoostProdToLabLifeTimeInDecays((*dau), lab, prod);  //recursively modify everything in the decay chain
0818     }
0819   }
0820 }
[ 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!