Rivet analyses

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/TauFinder.hh"
#include "Rivet/Tools/ParticleUtils.hh"

namespace Rivet {


  class PDG_TAUS : public Analysis {
  public:

    /// Constructor
    PDG_TAUS()
        : Analysis("PDG_TAUS") { }


    /// @name Analysis methods
    /// @{

    /// Book histograms and initialise projections before the run
    void init() {

      TauFinder tauleptonic(TauDecay::LEPTONIC); // open cuts, leptonic decays
      declare(tauleptonic, "TauLeptonic");

      TauFinder tauhadronic(TauDecay::HADRONIC); // open cuts, hadronic decays
      declare(tauhadronic, "TauHadronic");

      populateDecayMap();

      book(_h_ratio_mu, 1, 1, 1);
      book(_h_ratio_el, 1, 1, 2);
      book(_h_1prong_pinu, 2, 1, 1);
      book(_h_1prong_Kpnu, 2, 1, 2);
      book(_h_1prong_pipinu, 2, 1, 3);
      book(_h_1prong_Kppinu, 2, 1, 4);
      book(_h_1prong_pipipinu, 2, 1, 5);
      book(_h_1prong_Knpinu, 2, 1, 6);
      book(_h_3prong_pipipinu, 2, 2, 1);
      book(_h_5prong, 2, 3, 1);

      book(_weights_had, "TMP/weights_had");
      book(_weights_mu, "TMP/weights_mu");
      book(_weights_el, "TMP/weights_el");
    }


    /// Perform the per-event analysis
    void analyze(const Event& e) {
      const TauFinder& taulep = apply<TauFinder>(e, "TauLeptonic");
      const TauFinder& tauhad = apply<TauFinder>(e, "TauHadronic");

      // Hadronic tau decays --- prong decays
      for (const Particle& tau : tauhad.taus()) {
        _weights_had->fill();
        int prongs = countProngs(tau); // number of charged particles among decay products
        // Only do 1 prong decays here
        if (prongs == 1) {
          ////// Exclusive decay modes "1-prong"
          if (analyzeDecay(tau, decay_pids["pinu"], true)) _h_1prong_pinu->fill(1);
          if (analyzeDecay(tau, decay_pids["Kpnu"], true)) _h_1prong_Kpnu->fill(1);
          if (analyzeDecay(tau, decay_pids["pipinu"], true)) _h_1prong_pipinu->fill(1);
          if (analyzeDecay(tau, decay_pids["Kppinu"], true)) _h_1prong_Kppinu->fill(1);
          if (analyzeDecay(tau, decay_pids["pipipinu"], true)) _h_1prong_pipipinu->fill(1);
          // Kshort, Klong --- (twice) filling the K0 labelled PDG histo
          if (analyzeDecay(tau, decay_pids["KSpinu"], true)) _h_1prong_Knpinu->fill(1);
          if (analyzeDecay(tau, decay_pids["KLpinu"], true)) _h_1prong_Knpinu->fill(1);
        }
        else if (prongs == 3) {
          if (analyzeDecay(tau, decay_pids["3pipipinu"], true)) _h_3prong_pipipinu->fill(1);
        }
        else if (prongs == 5 && !any(tau.stableDescendants(), HasAbsPID(310)))
          _h_5prong->fill(1);
      }

      // Leptonic tau decays --- look for radiative and non-radiative 1 prong decays
      for (const Particle& tau : taulep.taus()) {
        int prongs = countProngs(tau); // number of charged particles among decay products
        // Only do 1 prong decays here
        if (prongs == 1) {
          analyzeRadiativeDecay(tau, decay_pids["muids"], _weights_mu, true, _h_ratio_mu);
          analyzeRadiativeDecay(tau, decay_pids["elids"], _weights_el, true, _h_ratio_el);
        }
      }
    }


    /// Normalise histograms etc., after the run
    void finalize() {
      scale(_h_ratio_mu, 1. / *_weights_mu);
      scale(_h_ratio_el, 1. / *_weights_el);

      const YODA::Counter norm = *_weights_had + *_weights_mu + *_weights_el;
      scale(_h_1prong_pinu, 1. / norm);
      scale(_h_1prong_Kpnu, 1. / norm);
      scale(_h_1prong_pipinu, 1. / norm);
      scale(_h_1prong_Kppinu, 1. / norm);
      scale(_h_1prong_pipipinu, 1. / norm);
      scale(_h_1prong_Knpinu, 1. / norm);
      scale(_h_3prong_pipipinu, 1. / norm);
      scale(_h_5prong, 1. / norm);
    }


    // Count charged decay products
    int countProngs(Particle mother) {
      int n_prongs = 0;
      for (Particle p : mother.stableDescendants())
        if (p.charge3() != 0) ++n_prongs;
      return n_prongs;
    }


    // Set up a lookup table for decays
    void populateDecayMap() {
      decay_pids["muids"] = {{13, 14, 16}};
      decay_pids["elids"] = {{11, 12, 16}};
      decay_pids["pinu"] = {{211, 16}};
      decay_pids["Kpnu"] = {{321, 16}};
      decay_pids["pipinu"] = {{111, 211, 16}};
      decay_pids["Kppinu"] = {{111, 321, 16}};
      decay_pids["pipipinu"] = {{111, 111, 211, 16}};
      decay_pids["KSpinu"] = {{211, 310, 16}};
      decay_pids["KLpinu"] = {{211, 130, 16}};
      decay_pids["3pipipinu"] = {{211, 211, 211, 16}};
    }


    bool analyzeDecay(Particle mother, const vector<int>& ids, bool absolute) {
      const Particles parts = {mother};
      return cascadeContains(parts, ids, absolute, true);
    }


    // Look for radiative (and non-radiative) tau decays to fill a ratio histo
    void analyzeRadiativeDecay(Particle mother,
                               vector<int> ids,
                               CounterPtr& w_incl,
                               bool absolute,
                               Histo1DPtr h_ratio) {
      // w_incl   ... reference to a global weight counter for all leptonic tau decays
      // h_ratio  ... pointer to ratio histo

      // There is no point in looking for decays with less particles than to be analysed
      const Particles& descendants = mother.stableDescendants();
      if (descendants.size() >= ids.size()) {
        const Particles parts = {mother};
        bool decayfound = cascadeContains(parts, ids, absolute, true);
        // Do not increment counters if the specified decay products were not found
        if (decayfound) {
          w_incl->fill(); // the (global) weight counter for leptonic decays
          bool radiative = any(descendants, HasPID(PID::PHOTON));

          // Only fill the histo if there is a radiative decay
          if (radiative) {
            // Iterate over decay products to find photon with 5 MeV energy
            for (const Particle& son : mother.stableDescendants()) {
              if (son.pid() == PID::PHOTON) {
                // Require photons to have at least 5 MeV energy in the rest frame of the tau
                // boosted taus
                if (!mother.momentum().betaVec().isZero()) {
                  LorentzTransform cms_boost = LorentzTransform::mkFrameTransformFromBeta(
                      mother.momentum().betaVec());
                  if (cms_boost.transform(son.momentum())[0] / MeV > 5.) {
                    h_ratio->fill(1);
                    break;
                  }
                }
                // not boosted taus
                else {
                  if (son.momentum()[0] / MeV > 5.) {
                    h_ratio->fill(1);
                    break;
                  }
                }
              }
            } // end loop over decay products
          } // end of radiative
        } // end of decayfound
      } // end of first if
    }


  private:

    /// @name Histograms
    /// @{
    Histo1DPtr _h_ratio_mu, _h_ratio_el;
    Histo1DPtr _h_1prong_pinu, _h_1prong_Kpnu, _h_1prong_Kppinu, _h_1prong_pipinu, _h_1prong_pipipinu,
        _h_1prong_Knpinu;
    Histo1DPtr _h_3prong_pipipinu;
    Histo1DPtr _h_5prong;
    /// @}

    CounterPtr _weights_had, _weights_mu, _weights_el;
    map<string, vector<int>> decay_pids;
  };


  RIVET_DECLARE_PLUGIN(PDG_TAUS);

}