Pythia8Interface/src/ResonanceDecayFilterHook.cc

0001 #include "Pythia8/Pythia.h"
0002 #include "GeneratorInterface/Pythia8Interface/interface/ResonanceDecayFilterHook.h"
0003
0004 using namespace Pythia8;
0005
0006 //--------------------------------------------------------------------------
0007 bool ResonanceDecayFilterHook::initAfterBeams() {
0008   filter_ = settingsPtr->flag("ResonanceDecayFilter:filter");
0009   exclusive_ = settingsPtr->flag("ResonanceDecayFilter:exclusive");
0010   eMuAsEquivalent_ = settingsPtr->flag("ResonanceDecayFilter:eMuAsEquivalent");
0011   eMuTauAsEquivalent_ = settingsPtr->flag("ResonanceDecayFilter:eMuTauAsEquivalent");
0012   allNuAsEquivalent_ = settingsPtr->flag("ResonanceDecayFilter:allNuAsEquivalent");
0013   udscAsEquivalent_ = settingsPtr->flag("ResonanceDecayFilter:udscAsEquivalent");
0014   udscbAsEquivalent_ = settingsPtr->flag("ResonanceDecayFilter:udscbAsEquivalent");
0015   wzAsEquivalent_ = settingsPtr->flag("ResonanceDecayFilter:wzAsEquivalent");
0016   auto mothers = settingsPtr->mvec("ResonanceDecayFilter:mothers");
0017   mothers_.clear();
0018   mothers_.insert(mothers.begin(), mothers.end());
0019   daughters_ = settingsPtr->mvec("ResonanceDecayFilter:daughters");
0020
0021   requestedDaughters_.clear();
0022
0023   for (int id : daughters_) {
0024     int did = std::abs(id);
0025     if (did == 13 && (eMuAsEquivalent_ || eMuTauAsEquivalent_)) {
0026       did = 11;
0027     }
0028     if (did == 15 && eMuTauAsEquivalent_) {
0029       did = 11;
0030     }
0031     if ((did == 14 || did == 16) && allNuAsEquivalent_) {
0032       did = 12;
0033     }
0034     if ((did == 2 || did == 3 || did == 4) && udscAsEquivalent_) {
0035       did = 1;
0036     }
0037     if ((did == 2 || did == 3 || did == 4 || did == 5) && udscbAsEquivalent_) {
0038       did = 1;
0039     }
0040     if ((did == 23 || did == 24) && wzAsEquivalent_) {
0041       did = 23;
0042     }
0043
0044     ++requestedDaughters_[std::abs(did)];
0045   }
0046
0047   return true;
0048 }
0049
0050 //--------------------------------------------------------------------------
0051 bool ResonanceDecayFilterHook::checkVetoResonanceDecays(const Event &process) {
0052   if (!filter_)
0053     return false;
0054
0055   observedDaughters_.clear();
0056
0057   //count decay products
0058   for (int i = 0; i < process.size(); ++i) {
0059     const Particle &p = process[i];
0060
0061     int did = std::abs(p.id());
0062
0063     if (did == 13 && (eMuAsEquivalent_ || eMuTauAsEquivalent_)) {
0064       did = 11;
0065     }
0066     if (did == 15 && eMuTauAsEquivalent_) {
0067       did = 11;
0068     }
0069     if ((did == 14 || did == 16) && allNuAsEquivalent_) {
0070       did = 12;
0071     }
0072     if ((did == 2 || did == 3 || did == 4) && udscAsEquivalent_) {
0073       did = 1;
0074     }
0075     if ((did == 2 || did == 3 || did == 4 || did == 5) && udscbAsEquivalent_) {
0076       did = 1;
0077     }
0078     if ((did == 23 || did == 24) && wzAsEquivalent_) {
0079       did = 23;
0080     }
0081
0082     int mid = p.mother1() > 0 ? std::abs(process[p.mother1()].id()) : 0;
0083
0084     //if no list of mothers is provided, then all particles
0085     //in hard process and resonance decays are counted together
0086     if (mothers_.empty() || mothers_.count(mid) || mothers_.count(-mid))
0087       ++observedDaughters_[did];
0088   }
0089
0090   //check if criteria is satisfied
0091   //inclusive mode: at least as many decay products as requested
0092   //exclusive mode: exactly as many decay products as requested
0093   //(but additional particle types not appearing in the list of requested daughter id's are ignored)
0094   for (const auto &reqpair : requestedDaughters_) {
0095     int reqid = reqpair.first;
0096     int reqcount = reqpair.second;
0097
0098     int obscount = 0;
0099     for (const auto &obspair : observedDaughters_) {
0100       int obsid = obspair.first;
0101
0102       if (obsid == reqid) {
0103         obscount = obspair.second;
0104         break;
0105       }
0106     }
0107
0108     //inclusive criteria not satisfied, veto event
0109     if (obscount < reqcount)
0110       return true;
0111
0112     //exclusive criteria not satisfied, veto event
0113     if (exclusive_ && obscount > reqcount)
0114       return true;
0115   }
0116
0117   //all criteria satisfied, don't veto
0118   return false;
0119 }