Pythia8Interface/src/BiasedTauDecayer.cc

0001 #include "Pythia8/Pythia.h"
0002 #include "GeneratorInterface/Pythia8Interface/interface/BiasedTauDecayer.h"
0003 using namespace Pythia8;
0004
0005 //==========================================================================
0006
0007 // Specialized decayer for resonance decays to taus to allowing biasing to
0008 // leptonic decays
0009 //
0010
0011 BiasedTauDecayer::BiasedTauDecayer(Info* infoPtr, Settings* settingsPtr) {
0012   decayer = TauDecays();
0013   decayer.initInfoPtr(*infoPtr);
0014   decayer.init();
0015   filter_ = settingsPtr->flag("BiasedTauDecayer:filter");
0016   eDecays_ = settingsPtr->flag("BiasedTauDecayer:eDecays");
0017   muDecays_ = settingsPtr->flag("BiasedTauDecayer:muDecays");
0018 }
0019
0020 bool BiasedTauDecayer::decay(
0021     std::vector<int>& idProd, std::vector<double>& mProd, std::vector<Vec4>& pProd, int iDec, const Event& event) {
0022   if (!filter_)
0023     return false;
0024   if (idProd[0] != 15 && idProd[0] != -15)
0025     return false;
0026   int iStart = event[iDec].iTopCopyId();
0027   int iMom = event[iStart].mother1();
0028   int idMom = event[iMom].idAbs();
0029   if (idMom != 23 && idMom != 24 && idMom != 25)
0030     return false;
0031   int iDau1 = event[iMom].daughter1();
0032   int iDau2 = event[iMom].daughter2();
0033   int iBot1 = event[iDau1].iBotCopyId();
0034   int iBot2 = event[iDau2].iBotCopyId();
0035   int iDecSis = (iDec == iBot1) ? iBot2 : iBot1;
0036   // Check if sister has been decayed
0037   // Since taus decays are correlated, use one decay, store the other
0038   bool notDecayed = event[iDecSis].status() > 0 ? true : false;
0039   if (notDecayed) {
0040     // bias for leptonic decays
0041     bool hasLepton = (eDecays_ || muDecays_) ? false : true;
0042     Event decay;
0043     int i1 = -1;
0044     int i2 = -1;
0045     while (!hasLepton) {
0046       decay = event;
0047       decayer.decay(iDec, decay);
0048       // check for lepton in first decay
0049       i1 = decay[iDec].daughter1();
0050       i2 = decay[iDec].daughter2();
0051       for (int i = i1; i < i2 + 1; ++i) {
0052         if (decay[i].isLepton() && decay[i].isCharged()) {
0053           if ((eDecays_ && abs(decay[i].id()) == 11) || (muDecays_ && abs(decay[i].id()) == 13)) {
0054             hasLepton = true;
0055             break;
0056           }
0057         }
0058       }
0059       if (hasLepton)
0060         break;
0061       // check for lepton in second decay
0062       i1 = decay[iDecSis].daughter1();
0063       i2 = decay[iDecSis].daughter2();
0064       for (int i = i1; i < i2 + 1; ++i) {
0065         if (decay[i].isLepton() && decay[i].isCharged()) {
0066           if ((eDecays_ && abs(decay[i].id()) == 11) || (muDecays_ && abs(decay[i].id()) == 13)) {
0067             hasLepton = true;
0068             break;
0069           }
0070         }
0071       }
0072     }
0073     // Return decay products
0074     i1 = decay[iDec].daughter1();
0075     i2 = decay[iDec].daughter2();
0076     for (int i = i1; i < i2 + 1; ++i) {
0077       idProd.push_back(decay[i].id());
0078       mProd.push_back(decay[i].m());
0079       pProd.push_back(decay[i].p());
0080     }
0081     // Store correlated decay products
0082     i1 = decay[iDecSis].daughter1();
0083     i2 = decay[iDecSis].daughter2();
0084     idProdSave.clear();
0085     mProdSave.clear();
0086     pProdSave.clear();
0087     for (int i = i1; i < i2 + 1; ++i) {
0088       idProdSave.push_back(decay[i].id());
0089       mProdSave.push_back(decay[i].m());
0090       pProdSave.push_back(decay[i].p());
0091     }
0092   } else {
0093     // Return stored decay products
0094     for (size_t i = 0; i < idProdSave.size(); ++i) {
0095       idProd.push_back(idProdSave[i]);
0096       mProd.push_back(mProdSave[i]);
0097       pProd.push_back(pProdSave[i]);
0098     }
0099   }
0100
0101   return true;
0102 }