1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Tools/ParticleUtils.hh"
  4#include "Rivet/Projections/TauFinder.hh"
  5
  6namespace Rivet {
  7
  8
  9  class PDG_TAUS : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    PDG_TAUS()
 14      : Analysis("PDG_TAUS")
 15    {   }
 16
 17
 18    /// @name Analysis methods
 19    /// @{
 20
 21    /// Book histograms and initialise projections before the run
 22    void init() {
 23
 24      TauFinder tauleptonic(TauDecay::LEPTONIC); // open cuts, leptonic decays
 25      declare(tauleptonic, "TauLeptonic");
 26
 27      TauFinder tauhadronic(TauDecay::HADRONIC); // open cuts, hadronic decays
 28      declare(tauhadronic, "TauHadronic");
 29
 30      populateDecayMap();
 31
 32      book(_h_ratio_mu        ,1, 1, 1);
 33      book(_h_ratio_el        ,1, 1, 2);
 34      book(_h_1prong_pinu     ,2, 1, 1);
 35      book(_h_1prong_Kpnu     ,2, 1, 2);
 36      book(_h_1prong_pipinu   ,2, 1, 3);
 37      book(_h_1prong_Kppinu   ,2, 1, 4);
 38      book(_h_1prong_pipipinu ,2, 1, 5);
 39      book(_h_1prong_Knpinu   ,2, 1, 6);
 40      book(_h_3prong_pipipinu ,2, 2, 1);
 41      book(_h_5prong          ,2, 3, 1);
 42
 43      book(_weights_had, "TMP/weights_had");
 44      book(_weights_mu, "TMP/weights_mu");
 45      book(_weights_el, "TMP/weights_el");
 46    }
 47
 48
 49    /// Perform the per-event analysis
 50    void analyze(const Event& e) {
 51      const TauFinder& taulep = apply<TauFinder>(e, "TauLeptonic");
 52      const TauFinder& tauhad = apply<TauFinder>(e, "TauHadronic");
 53
 54      // Hadronic tau decays --- prong decays
 55      for(const Particle& tau : tauhad.taus()) {
 56        _weights_had->fill();
 57        int prongs = countProngs(tau); // number of charged particles among decay products
 58        // Only do 1 prong decays here
 59        if (prongs == 1) {
 60          ////// Exclusive decay modes "1-prong"
 61          if (analyzeDecay(tau,   decay_pids["pinu"], true))     _h_1prong_pinu->fill(1);
 62          if (analyzeDecay(tau,   decay_pids["Kpnu"], true))     _h_1prong_Kpnu->fill(1);
 63          if (analyzeDecay(tau, decay_pids["pipinu"], true))     _h_1prong_pipinu->fill(1);
 64          if (analyzeDecay(tau, decay_pids["Kppinu"]  , true))   _h_1prong_Kppinu->fill(1);
 65          if (analyzeDecay(tau, decay_pids["pipipinu"], true))   _h_1prong_pipipinu->fill(1);
 66          // Kshort, Klong --- (twice) filling the K0 labelled PDG histo
 67          if (analyzeDecay(tau, decay_pids["KSpinu"]  , true))   _h_1prong_Knpinu->fill(1);
 68          if (analyzeDecay(tau, decay_pids["KLpinu"]  , true))   _h_1prong_Knpinu->fill(1);
 69        }
 70        else if (prongs == 3) {
 71          if (analyzeDecay(tau, decay_pids["3pipipinu"], true))  _h_3prong_pipipinu->fill(1);
 72        }
 73        else if (prongs == 5 && !any(tau.stableDescendants(), HasAbsPID(310))) _h_5prong->fill(1);
 74      }
 75
 76      // Leptonic tau decays --- look for radiative and non-radiative 1 prong decays
 77      for(const Particle& tau : taulep.taus()) {
 78        int prongs = countProngs(tau); // number of charged particles among decay products
 79        // Only do 1 prong decays here
 80        if (prongs == 1) {
 81          analyzeRadiativeDecay(tau, decay_pids["muids"], _weights_mu, true, _h_ratio_mu);
 82          analyzeRadiativeDecay(tau, decay_pids["elids"], _weights_el, true, _h_ratio_el);
 83        }
 84      }
 85    }
 86
 87
 88    /// Normalise histograms etc., after the run
 89    void finalize() {
 90      scale(_h_ratio_mu, 1. / *_weights_mu);
 91      scale(_h_ratio_el, 1. / *_weights_el);
 92
 93      const YODA::Counter norm = *_weights_had + *_weights_mu + *_weights_el;
 94      scale(_h_1prong_pinu,     1./norm);
 95      scale(_h_1prong_Kpnu,     1./norm);
 96      scale(_h_1prong_pipinu,   1./norm);
 97      scale(_h_1prong_Kppinu,   1./norm);
 98      scale(_h_1prong_pipipinu, 1./norm);
 99      scale(_h_1prong_Knpinu,   1./norm);
100      scale(_h_3prong_pipipinu, 1./norm);
101      scale(_h_5prong,          1./norm);
102    }
103
104
105    // Count charged decay products
106    int countProngs(Particle mother) {
107      int n_prongs = 0;
108      for(Particle p : mother.stableDescendants())
109        if (p.charge3()!=0) ++n_prongs;
110      return n_prongs;
111    }
112
113
114    // Set up a lookup table for decays
115    void populateDecayMap() {
116      decay_pids["muids"]     = {{ 13,14,16 }};
117      decay_pids["elids"]     = {{ 11,12,16 }};
118      decay_pids["pinu"]      = {{ 211,16 }};
119      decay_pids["Kpnu"]      = {{ 321,16 }};
120      decay_pids["pipinu"]    = {{ 111,211,16 }};
121      decay_pids["Kppinu"]    = {{ 111,321,16 }};
122      decay_pids["pipipinu"]  = {{ 111,111,211,16 }};
123      decay_pids["KSpinu"]    = {{ 211,310,16 }};
124      decay_pids["KLpinu"]    = {{ 211,130,16 }};
125      decay_pids["3pipipinu"] = {{ 211,211,211,16 }};
126    }
127
128
129    bool analyzeDecay(Particle mother, const vector<int>& ids, bool absolute) {
130      const Particles parts = { mother };
131      return cascadeContains(parts, ids, absolute, true);
132    }
133
134
135    // Look for radiative (and non-radiative) tau decays to fill a ratio histo
136    void analyzeRadiativeDecay(Particle mother, vector<int> ids, CounterPtr &w_incl, bool absolute, Histo1DPtr h_ratio) {
137      // w_incl   ... reference to a global weight counter for all leptonic tau decays
138      // h_ratio  ... pointer to ratio histo
139
140      // There is no point in looking for decays with less particles than to be analysed
141      const Particles& descendants = mother.stableDescendants();
142      if (descendants.size() >= ids.size()) {
143        const Particles parts = { mother };
144        bool decayfound = cascadeContains(parts, ids, absolute, true);
145        // Do not increment counters if the specified decay products were not found
146        if (decayfound) {
147          w_incl->fill(); // the (global) weight counter for leptonic decays
148          bool radiative = any(descendants, HasPID(PID::PHOTON));
149
150          // Only fill the histo if there is a radiative decay
151          if (radiative) {
152            // Iterate over decay products to find photon with 5 MeV energy
153            for (const Particle& son : mother.stableDescendants()) {
154              if (son.pid() == PID::PHOTON) {
155                // Require photons to have at least 5 MeV energy in the rest frame of the tau
156                // boosted taus
157                if (!mother.momentum().betaVec().isZero()) {
158                  LorentzTransform cms_boost = LorentzTransform::mkFrameTransformFromBeta(mother.momentum().betaVec());
159                  if (cms_boost.transform(son.momentum())[0]/MeV > 5.) {
160                    h_ratio->fill(1);
161                    break;
162                  }
163                }
164                // not boosted taus
165                else {
166                  if (son.momentum()[0]/MeV > 5.) {
167                    h_ratio->fill(1);
168                    break;
169                  }
170                }
171              }
172            } // end loop over decay products
173          } // end of radiative
174        } // end of decayfound
175      } // end of first if
176    }
177
178
179  private:
180
181    /// @name Histograms
182    /// @{
183    Histo1DPtr _h_ratio_mu, _h_ratio_el;
184    Histo1DPtr _h_1prong_pinu, _h_1prong_Kpnu, _h_1prong_Kppinu, _h_1prong_pipinu, _h_1prong_pipipinu, _h_1prong_Knpinu;
185    Histo1DPtr _h_3prong_pipipinu;
186    Histo1DPtr _h_5prong;
187    /// @}
188
189    CounterPtr _weights_had, _weights_mu, _weights_el;
190    map<string, vector<int> > decay_pids;
191
192  };
193
194
195  RIVET_DECLARE_PLUGIN(PDG_TAUS);
196
197}