Rivet analyses

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

namespace Rivet {


  /// @brief Monte Carlo validation observables for tau polarisation
  class MC_TAUPOL : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(MC_TAUPOL);

    /// @name Analysis methods
    //@{

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

      _mode = 0; // default is ditau
      if ( getOption("MODE") == "DITAU" ) { _mode = 0; }
      if ( getOption("MODE") == "SINGLE" ) { _mode = 1; }
      _lothresh = getOption("MASSMIN", _mode? 75*GeV : 86.*GeV);
      _hithresh = getOption("MASSMAX", _mode? 85*GeV : 96*GeV);

      // Initialise and register projections
      declare(TauFinder(TauDecay::ANY), "taus");

      // Book histograms
      // tau invmass
      book(_h["m_tau"], "h_m_tau", 50,1.,3.);
      book(_h["m_tautau"], "h_m_tautau", 50,35.,420.);
      book(_h["m_B"], "h_m_B", 50,35.,420.);
      book(_h["B_pT"], "h_B_pT", 100,0.,500.);

      // longnitudal spin correlations
      book(_h["m_tau_vis"], "h_m_tau_vis",50,0.,5.);
      book(_h["m_tautau_vis"], "h_m_tautau_vis",50,0.,100.);
      book(_h["m_pipi_vis"], "h_m_pipi_vis",60,0.,100.);
      book(_h["m_pipi_vis_lm"], "h_m_pipi_vis_lm",60,0.,100.);
      book(_h["m_pipi_vis_hm"], "h_m_pipi_vis_hm",60,0.,100.);

      // transverse spin correlations
      book(_h["acoplanarity_plus"], "h_acoplanarity_plus",30,0.,6.2);
      book(_h["acop_pt40"], "h_acop_pt40",30,0.,6.2);
      book(_h["acop_pt100"], "h_acop_pt100",30,0.,6.2);
      book(_h["acop_pt200"], "h_acop_pt200",30,0.,6.2);
      book(_h["acop_pt1000"], "h_acop_pt1000",30,0.,6.2);
      book(_h["acom_pt40"], "h_acom_pt40",30,0.,6.2);
      book(_h["acom_pt100"], "h_acom_pt100",30,0.,6.2);
      book(_h["acom_pt200"], "h_acom_pt200",30,0.,6.2);
      book(_h["acom_pt1000"], "h_acom_pt1000",30,0.,6.2);
      book(_h["acoplanarity_minus"], "h_acoplanarity_minus",30,0.,6.2);

      // leptonic

      book(_h["2B_xlep"], "h_2B_xlep", 50,0.0,1.0);
      book(_h["2B_xel"], "h_2B_xel", 50,0.0,1.0);
      book(_h["2B_xmu"], "h_2B_xmu", 50,0.0,1.0);
      book(_h["xel_B"], "h_xel_B", 50,0.0,1.0);
      book(_h["xmu_B"], "h_xmu_B", 50,0.0,1.0);

      // hadronic

      book(_h["1B_xpi"], "h_1B_xpi", 30,0.0,1.0);
      book(_h["1B_xpi_lm"], "h_1B_xpi_lm", 30,0.0,1.0);
      book(_h["1B_xpi_hm"], "h_1B_xpi_hm", 30,0.0,1.0);
      book(_h["xpi_B"], "h_xpi_B", 30,0.0,1.0);
      book(_h["2B_xrho"], "h_2B_xrho", 50,0.0,1.0);

      // angles
      book(_h["t_pi"], "h_t_pi", 40, -1.0, 1.0);
      book(_h["t_rho"], "h_t_rho", 40, -1.0, 1.0);
      book(_h["pipi_theta"], "h_pipi_theta", 40, -1.0, 1.0);
      book(_ntaus, "h_n_taus", {0, 1, 2, 3, 4});

    }

    void findDecayProducts(const Particle& mother, size_t& nstable,
               Particles& ep , Particles& em , Particles& nu_e , Particles& nu_ebar,
               Particles& mup, Particles& mum, Particles& nu_mu, Particles& nu_mubar,
               Particles& pip, Particles& pim, Particles& pi0,
               Particles& Kp , Particles& Km , Particles& K0S,   Particles& K0L,
               Particles& eta, Particles& gamma) {
      for (const Particle& p : mother.children()) {
        int id = p.pid();
        if (id == PID::KPLUS) {
          Kp.push_back(p);
          ++nstable;
        }
        else if (id == PID::KMINUS ) {
          Km.push_back(p);
          ++nstable;
        }
        else if (id == PID::PIPLUS) {
          pip.push_back(p);
          ++nstable;
        }
        else if (id == PID::PIMINUS) {
          pim.push_back(p);
          ++nstable;
        }
        else if (id == PID::EPLUS) {
          ep.push_back(p);
          ++nstable;
        }
        else if (id == PID::EMINUS) {
          em.push_back(p);
          ++nstable;
        }
        else if (id == PID::NU_E) {
          nu_e.push_back(p);
          ++nstable;
        }
        else if (id == PID::NU_EBAR) {
          nu_ebar.push_back(p);
          ++nstable;
        }
        else if (id == PID::NU_MU) {
          nu_mu.push_back(p);
          ++nstable;
        }
        else if (id == PID::NU_MUBAR) {
          nu_mubar.push_back(p);
          ++nstable;
        }
        else if (id == PID::ANTIMUON) {
          mup.push_back(p);
          ++nstable;
        }
        else if (id == PID::MUON) {
          mum.push_back(p);
          ++nstable;
        }
        else if (id == PID::PI0) {
          pi0.push_back(p);
          ++nstable;
        }
        else if (id == PID::K0S) {
          K0S.push_back(p);
          ++nstable;
        }
        else if (id == PID::K0L) {
          K0L.push_back(p);
          ++nstable;
        }
        else if (id == PID::ETA) {
          eta.push_back(p);
          ++nstable;
        }
        else if (id == PID::PHOTON) {
          gamma.push_back(p);
          ++nstable;
        }
        else if ( !p.children().empty() ) {
          findDecayProducts(p, nstable, ep, em, nu_e, nu_ebar, mup, mum, nu_mu, nu_mubar,
                            pip, pim, pi0, Kp, Km, K0S, K0L, eta, gamma);
        }
        else {
          ++nstable;
        }
      }
    }


    double angleRF(const FourMomentum& p1, const FourMomentum& p2, const FourMomentum& RF) const {
       LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(RF.betaVec());
       return cos(boost(p1).angle(boost(p2)));
    }

    /// Perform the per-event analysis
    void analyze(const Event& event) {

      Particles taus, t_nu, pions;
      for (const auto& tau : apply<TauFinder>(event, "taus").particlesByPt()) {
        //remove any from hadronic decay (only want prompt taus)
        if (tau.parents()[0].isPrompt()) {
          taus += tau;
          _h["m_tau"]->fill(tau.mom().mass()/GeV);
        }
      }
      //veto ditau events without exactly 2 prompt taus
      if (taus.size() != 2 && _mode == 0) vetoEvent;
      //veto single-tau events without exactly 1 prompt tau
      if (taus.size() != 1 && _mode == 1) vetoEvent;

      FourMomentum tautaumom, taunumom, phomom, Bmom;
      if (taus.size() == 2 && _mode == 0) {
          tautaumom = taus[0].mom() + taus[1].mom();
        phomom = sum(scanAncestors(taus,22), Kin::p4, FourMomentum());
        Bmom = tautaumom + phomom;
        _h["m_tautau"]->fill(tautaumom.mass());
      }
      if (taus.size()==1 && _mode == 1) {
        t_nu = scanAncestors(taus,16);
        phomom = sum(scanAncestors(taus,22), Kin::p4, FourMomentum());
        if(t_nu.size() != 1)vetoEvent;
        if(taus[0].pid()*t_nu[0].pid() > 0) vetoEvent;
        taunumom = taus[0].mom() + t_nu[0].mom();
        Bmom = taunumom + phomom;
      }
      _ntaus->fill(int(taus.size()));
      _h["m_B"]->fill(Bmom.mass());
      _h["B_pT"]->fill(Bmom.pT());
       // Check opposite charges

      // Calculate invariant mass of tautau
      FourMomentum tt_mom_vis;
      for (const auto& t : taus) {
        FourMomentum t_mom_vis;
        for (const auto& c : t.children()) {
          if (c.isVisible()) {
            tt_mom_vis += c.mom();
            t_mom_vis += c.mom();
          }
        }
        _h["m_tau_vis"]->fill(t_mom_vis.mass());
      }

      _h["m_tautau_vis"]->fill(tt_mom_vis.mass());
      size_t OS = 2;
      if (Bmom.mass() < _lothresh)  OS = 0;
      else if (Bmom.mass() < _hithresh)  OS = 1;

      for (const Particle& tau :taus) {
        size_t nstable(0);
        Particles ep,em,nu_e,nu_ebar,mup,mum,nu_mu,nu_mubar;
        Particles pip, pim, pi0, Kp , Km, K0S, K0L, eta,gamma;
        findDecayProducts(tau, nstable,ep,em,nu_e,nu_ebar,mup,mum,nu_mu,nu_mubar,
                          pip, pim, pi0, Kp, Km, K0S, K0L,eta,gamma);
        if (tau.pid() < 0) {
          swap(pim,pip);
          swap(Kp,Km);
          swap(em,ep);
          swap(mum,mup);
          swap(nu_e ,nu_ebar );
          swap(nu_mu,nu_mubar);
        }
        // classify by number of stable decay products
        FourMomentum tmom = tau.mom();
        const LorentzTransform tboost = LorentzTransform::mkFrameTransformFromBeta(tmom.betaVec());
        if (nstable==2) {
          if (pim.size()==1) {
            const double xpi = pim[0].mom().p()/tau.mom().p();
            const FourMomentum pimom = pim[0].mom();
            const double angle = cos(tboost(pimom).angle(tmom));
            _h["xpi_B"]->fill(pim[0].mom().p()/Bmom.p());
            _h["t_pi"]->fill(angle);
            if (OS==0)      _h["1B_xpi_lm"]->fill(xpi);
            else if (OS==1) _h["1B_xpi"]->fill(xpi);
            else            _h["1B_xpi_hm"]->fill(xpi);
            pions.push_back(pim[0]);
          }
        }
        else if (nstable==3 ) {
          if (em.size()==1 && nu_ebar.size()==1) {
            const double xlep = em[0].mom().p()/tau.mom().p();
            _h["2B_xlep"]->fill(xlep);
            _h["2B_xel"]->fill(xlep);
            _h["xel_B"]->fill(em[0].mom().p()/Bmom.p());
          }
          else if (mum.size()==1 && nu_mubar.size()==1) {
            const double xlep = mum[0].mom().p()/tau.mom().p();
            _h["2B_xlep"]->fill(xlep);
            _h["2B_xmu"]->fill(xlep);
            _h["xmu_B"]->fill(mum[0].mom().p()/Bmom.p());
          }
          else if (pim.size()==1 && pi0.size()==1) {
            const FourMomentum ptot = pim[0].mom()+pi0[0].mom();
            const double xrho = ptot.p()/tau.mom().p();
            _h["2B_xrho"]->fill(xrho);
            const double angle = cos(tboost(ptot).angle(tmom));
            _h["t_rho"]->fill(angle);
          }
        } //end nstable = 3

      } //end tau loop

      // TAU LONGNITUDAL SPIN EFFECTS PART
      //fill tautau->pipi inv mass
      if (pions.size() == 2){
        FourMomentum pi1mom = pions[0].mom();
        FourMomentum pi2mom = pions[1].mom();
        FourMomentum pipi = pions[0].mom() + pions[1].mom();
        if (OS==0)      _h["m_pipi_vis_lm"]->fill(pipi.mass()/GeV);
        else if (OS==1) _h["m_pipi_vis"]->fill(pipi.mass()/GeV);
        else            _h["m_pipi_vis_hm"]->fill(pipi.mass()/GeV);
        _h["pipi_theta"]->fill(angleRF(pi1mom,pi2mom,Bmom));
      }

        // TAU TRANSVERSE SPIN EFFECTS PART
      if (pions.size() == 2 && _mode != 1 && OS == 1){

        FourMomentum tautau = taus[0].mom() + taus[1].mom();

        // calculate angles between (tau, tau_pi) and (tau2, tau2_pi) planes
        const LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(tautau.betaVec());
          taus[0] = Particle(PID::ANY, boost.transform(taus[0].mom()));
          pions[0] = Particle(PID::ANY, boost.transform(pions[0].mom()));
        taus[1] = Particle(PID::ANY, boost.transform(taus[1].mom()));
          pions[1] = Particle(PID::ANY, boost.transform(pions[1].mom()));

        // calculate the angle
        Vector3 t0xp0 = taus[0].p3().cross(pions[0].p3()).unit();
        Vector3 t1xp1 = taus[1].p3().cross(pions[1].p3()).unit();
        double phi_star_sign = acos(t0xp0.dot(t1xp1));

        // extend acoplanarity angle to range (0, 2pi)
        Vector3 p0xp1 = pions[0].p3().cross(pions[1].p3()).unit();
        double sign = p0xp1.dot(taus[0].p3().unit());
        if (sign > 0)  phi_star_sign = 2* M_PI - phi_star_sign;


        // alpha angles in the lab frame
        Vector3  tlv3_Ez(0.0, 0.0, 1.0 );
        Vector3 cross1_pi = tlv3_Ez.cross(taus[0].p3()).unit();
        Vector3 cross2_pi = pions[0].p3().cross(taus[0].p3()).unit();
        const double alpha_pi = acos(abs(cross1_pi.dot(cross2_pi)));

        if (alpha_pi > M_PI/4.) {
          _h["acoplanarity_plus"]->fill(phi_star_sign);
          if (Bmom.pT() <  40*GeV)       _h["acop_pt40"]->fill(phi_star_sign);
          else if (Bmom.pT() < 100*GeV)  _h["acop_pt100"]->fill(phi_star_sign);
          else if (Bmom.pT() < 200*GeV)  _h["acop_pt200"]->fill(phi_star_sign);
          else                           _h["acop_pt1000"]->fill(phi_star_sign);
        }

        if (alpha_pi < M_PI/4.) {
          _h["acoplanarity_minus"]->fill(phi_star_sign);
          if (Bmom.pT() <  40*GeV)       _h["acom_pt40"]->fill(phi_star_sign);
          else if (Bmom.pT() < 100*GeV)  _h["acom_pt100"]->fill(phi_star_sign);
          else if (Bmom.pT() < 200*GeV)  _h["acom_pt200"]->fill(phi_star_sign);
          else                           _h["acom_pt1000"]->fill(phi_star_sign);
        }

      } // end of tau spin corr

    } //end void function


    /// Normalise histograms etc., after the run
    void finalize() {
      normalize(_h);
      normalize(_ntaus);
    }

    //@}

    Particles scanAncestors(const Particles& taus, int id){
      Particles rtn;
      for (const Particle& tau : taus) {
        Particle thistau = tau;
        while (true) {
          // check the tau has parents first else end loop
          if (thistau.parents().empty())  break;
          //loop over taus siblings
          for (const Particle& sibling : thistau.parents()[0].children()) {
            // find ID
            if (sibling.abspid() == id) {
              // add particles to rtn vector
              bool duplicate = true;
              for (const auto& p : rtn) {
                if (sibling.isSame(p)) duplicate = false;
              }
              if (duplicate) rtn += sibling;
            }
          }
          // want to loop over taus parent's siblings next
          thistau = thistau.parents()[0];
          // if taus parent is not a tau end the loop
          if (thistau.abspid() != PID::TAU)  break;
        }
      }
      return rtn;
    }

    /// @name Histograms
    //@{

    // histograms for leptonic decay
    map<string,Histo1DPtr> _h;

    BinnedHistoPtr<int> _ntaus;

    size_t _mode;

    double _lothresh, _hithresh;

    //@}

  };


  RIVET_DECLARE_PLUGIN(MC_TAUPOL);
}