1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief tau -> K+ KS nu_tau
  9  class BABAR_2018_I1679886 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(BABAR_2018_I1679886);
 14
 15
 16    /// @name Analysis methods
 17    /// @{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21
 22      declare(UnstableParticles(), "UFS");
 23      book(_h_KK, 1, 1, 1);
 24      book(_v_KK, 1, 1, 2);
 25      book(_c,"TMP/c");
 26    }
 27
 28    void findDecayProducts(const Particle & mother, unsigned int & nstable,
 29                           unsigned int & nK0, unsigned int & nKp,
 30			   unsigned int & nKm, FourMomentum & ptot) {
 31      for(const Particle & p : mother.children()) {
 32        int id = p.pid();
 33        if ( id == PID::KPLUS ) {
 34	  ++nKp;
 35          ++nstable;
 36	  ptot += p.momentum();
 37	}
 38        else if (id == PID::KMINUS ) {
 39	  ++nKm;
 40          ++nstable;
 41	  ptot += p.momentum();
 42	}
 43        else if (id == PID::K0S) {
 44          ++nK0;
 45          ++nstable;
 46	  ptot += p.momentum();
 47        }
 48        else if (id == PID::PI0 || id == PID::PIPLUS || id == PID::PIMINUS) {
 49          ++nstable;
 50        }
 51        else if ( !p.children().empty() ) {
 52          findDecayProducts(p, nstable, nK0, nKp, nKm, ptot);
 53        }
 54        else
 55          ++nstable;
 56      }
 57    }
 58
 59    /// Perform the per-event analysis
 60    void analyze(const Event& event) {
 61
 62      // Loop over taus
 63      for(const Particle& tau : apply<UnstableParticles>(event, "UFS").particles(Cuts::abspid==PID::TAU)) {
 64        unsigned int nstable(0),nK0(0),nKp(0),nKm(0);
 65      	FourMomentum p_tot(0,0,0,0);
 66        findDecayProducts(tau, nstable, nK0, nKp, nKm, p_tot);
 67        if (tau.pid() < 0) {
 68      	  swap(nKp,nKm);
 69      	}
 70       	if(nstable!=3) continue;
 71      	if(nKm==1 && nK0==1 ) {
 72          double q = p_tot.mass();
 73          _h_KK->fill(q);
 74          double mtau = tau.mass();
 75          double Cq   = q*sqr(sqr(mtau)-sqr(q))*(sqr(mtau)+2.*sqr(q));
 76          double fact = pow(mtau,8)/Cq;
 77          _v_KK->fill(q,fact);
 78          _c->fill();
 79        }
 80      }
 81
 82
 83    }
 84
 85
 86    /// Normalise histograms etc., after the run
 87    void finalize() {
 88      normalize(_h_KK);
 89      // values from PDG 2023 for constants, see eqn 1 of paper
 90      double brKK = 0.5*1.486e-3; // PDG is K0 divide by 2 to get KS0
 91      double brE  = 0.1782;
 92      double Vud  = .97367;
 93      scale(_v_KK, 1000.*brKK/brE/sqr(Vud)/12./M_PI/ *_c);
 94    }
 95
 96    /// @}
 97
 98
 99    /// @name Histograms
100    /// @{
101    Histo1DPtr _h_KK,_v_KK;
102    CounterPtr _c;
103    /// @}
104
105
106  };
107
108
109  RIVET_DECLARE_PLUGIN(BABAR_2018_I1679886);
110
111
112}