1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/TauFinder.hh"
  4#include "Rivet/Projections/IdentifiedFinalState.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief MC validation analysis for tau decays involving photons
 10  class MC_TAUS_PHOTONS : public Analysis {
 11    public:
 12
 13      RIVET_DEFAULT_ANALYSIS_CTOR(MC_TAUS_PHOTONS);
 14
 15      /// Book projections and histograms
 16      void init() {
 17        TauFinder taus(TauDecay::ANY, Cuts::pT > 500*MeV);
 18        declare(taus, "Taus");
 19
 20        IdentifiedFinalState photons(Cuts::pT > 500*MeV, PID::PHOTON);
 21        declare(photons, "Photons");
 22
 23        book(_nPhotonsEl, "nPhotonsEl", 20, 0, 20);
 24        book(_nPhotonsMu, "nPhotonsMu", 20, 0, 20);
 25        book(_nPhotonsHad, "nPhotonsHad", 20, 0, 20);
 26
 27        book(_tauMassEl, "tauMassEl", 50, 0, 5);
 28        book(_tauMassMu, "tauMassMu", 50, 0, 5);
 29        book(_tauMassHad, "tauMassHad", 50, 0, 5);
 30
 31        book(_pFracPhotonsEl, "pFracPhotonsEl", 50, 0, 2);
 32        book(_pFracPhotonsMu, "pFracPhotonsMu", 50, 0, 2);
 33        book(_pFracPhotonsHad, "pFracPhotonsHad", 50, 0, 2);
 34
 35        book(_logPFracPhotonsEl, "logPFracPhotonsEl", 50, -5, 0);
 36        book(_logPFracPhotonsMu, "logPFracPhotonsMu", 50, -5, 0);
 37        book(_logPFracPhotonsHad, "logPFracPhotonsHad", 50, -5, 0);
 38
 39        book(_logPFracNotPhotonsEl, "logPFracNotPhotonsEl", 50, -5, 0);
 40        book(_logPFracNotPhotonsMu, "logPFracNotPhotonsMu", 50, -5, 0);
 41        book(_logPFracNotPhotonsHad, "logPFracNotPhotonsHad", 50, -5, 0);
 42
 43        book(_restFramePhotonsEnergyEl, "RestFramePhotonsEnergyEl", 50, 0, 3);
 44        book(_restFramePhotonsEnergyMu, "RestFramePhotonsEnergyMu", 50, 0, 3);
 45        book(_restFramePhotonsEnergyHad, "RestFramePhotonsEnergyHad", 50, 0, 3);
 46
 47        return;
 48      }
 49
 50
 51      /// Per-event analysis
 52      void analyze(const Event& event) {
 53        const Particles taus = apply<TauFinder>(event, "Taus").particlesByPt();
 54        const Particles photons = apply<IdentifiedFinalState>(event, "Photons").particlesByPt();
 55
 56        for (const auto& tau : taus) {
 57          Particles descendants = tau.stableDescendants();
 58          bool hasHad = false, hasEl = false, hasMu = false;
 59
 60          for (const auto& descendant : descendants) {
 61            if (isHadron(descendant)) {
 62              hasHad = true;
 63            } else if (abs(descendant.pid()) == PID::ELECTRON) {
 64              hasEl = true;
 65            } else if (abs(descendant.pid()) == PID::MUON) {
 66              hasMu = true;
 67            }
 68
 69            continue;
 70          }
 71
 72          // decaymode: hadronic = 0, electron = 1, muon = 2
 73          int decaymode = -1;
 74          if (hasHad) decaymode = 0;
 75          else if (hasEl) decaymode = 1;
 76          else if (hasMu) decaymode = 2;
 77          assert(decaymode >= 0);
 78
 79          FourMomentum taumom;
 80          FourMomentum photonmom;
 81          int nphotons = 0;
 82          for (const auto& descendant : descendants) {
 83            if (descendant.pid() == PID::PHOTON) {
 84              photonmom += descendant.mom();
 85              nphotons++;
 86            }
 87
 88            taumom += descendant.mom();
 89
 90            continue;
 91          }
 92
 93          LorentzTransform boost =
 94            LorentzTransform::mkFrameTransformFromBeta(taumom.betaVec());
 95
 96          float taumass = taumom.mass();
 97
 98          float photonfrac = photonmom.p() / tau.p();
 99          if (photonfrac >= 1.0)
100            photonfrac = 1.0 - 1e-6;
101
102          if (decaymode == 0) {
103            _nPhotonsHad->fill(nphotons);
104            _tauMassHad->fill(taumass / GeV);
105
106            _pFracPhotonsHad->fill(photonfrac);
107            _logPFracPhotonsHad->fill(log(photonfrac));
108            _logPFracNotPhotonsHad->fill(log(1.0 - photonfrac));
109            _restFramePhotonsEnergyHad->fill(boost(photonmom).E() / GeV);
110
111          } else if (decaymode == 1) {
112            _nPhotonsEl->fill(nphotons);
113            _tauMassEl->fill(taumass / GeV);
114
115            _pFracPhotonsEl->fill(photonfrac);
116            _logPFracPhotonsEl->fill(log(photonfrac));
117            _logPFracNotPhotonsEl->fill(log(1.0 - photonfrac));
118            _restFramePhotonsEnergyEl->fill(boost(photonmom).E() / GeV);
119
120          } else if (decaymode == 2) {
121            _nPhotonsMu->fill(nphotons);
122            _tauMassMu->fill(taumass / GeV);
123
124            _pFracPhotonsMu->fill(photonfrac);
125            _logPFracPhotonsMu->fill(log(photonfrac));
126            _logPFracNotPhotonsMu->fill(log(1.0 - photonfrac));
127            _restFramePhotonsEnergyMu->fill(boost(photonmom).E() / GeV);
128
129          }
130
131          continue;
132        }
133
134        return;
135      }
136
137
138      /// Normalisations etc.
139      void finalize() {
140        _nPhotonsEl->normalize();
141        _nPhotonsMu->normalize();
142        _nPhotonsHad->normalize();
143
144        _tauMassEl->normalize();
145        _tauMassMu->normalize();
146        _tauMassHad->normalize();
147
148        _pFracPhotonsEl->normalize();
149        _pFracPhotonsMu->normalize();
150        _pFracPhotonsHad->normalize();
151
152        _logPFracPhotonsEl->normalize();
153        _logPFracPhotonsMu->normalize();
154        _logPFracPhotonsHad->normalize();
155
156        _logPFracNotPhotonsEl->normalize();
157        _logPFracNotPhotonsMu->normalize();
158        _logPFracNotPhotonsHad->normalize();
159
160        _restFramePhotonsEnergyEl->normalize();
161        _restFramePhotonsEnergyMu->normalize();
162        _restFramePhotonsEnergyHad->normalize();
163
164        return;
165      }
166
167    private:
168
169      Histo1DPtr
170          _nPhotonsEl
171        , _nPhotonsMu
172        , _nPhotonsHad
173        , _tauMassEl
174        , _tauMassMu
175        , _tauMassHad
176        , _pFracPhotonsEl
177        , _pFracPhotonsMu
178        , _pFracPhotonsHad
179        , _logPFracPhotonsEl
180        , _logPFracPhotonsMu
181        , _logPFracPhotonsHad
182        , _logPFracNotPhotonsEl
183        , _logPFracNotPhotonsMu
184        , _logPFracNotPhotonsHad
185        , _restFramePhotonsEnergyEl
186        , _restFramePhotonsEnergyMu
187        , _restFramePhotonsEnergyHad
188        ;
189
190  };
191
192
193  RIVET_DECLARE_PLUGIN(MC_TAUS_PHOTONS);
194
195}