Rivet analyses

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

namespace Rivet {


  /// @brief B -> c l nu moments
  class BELLE_2021_I1917200 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2021_I1917200);


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

    /// Book histograms and initialise projections before the run
    void init() {
      // projections
      declare(UnstableParticles(Cuts::abspid==511 || Cuts::abspid==521),"UFS");
      // histos
      for (unsigned int il=0; il<2; ++il) {
        for (unsigned int ix=0; ix<4; ++ix) {
          book(_p[il][ix], 4*il+1+ix, 1, 1);
        }
        for (unsigned int ix=4; ix<8; ++ix) {
          book(_p[il][ix], "TMP/p_"+toString(il)+"_"+toString(ix),
               refData<YODA::BinnedEstimate<string>>(4*il+1, 1, 1));
        }
      }
    }

    void findDecayProducts(const Particle& parent, Particles& em, Particles& ep,
                           Particles& nue, Particles & nueBar, bool& charm) {
      for (const Particle& p : parent.children()) {
        if(PID::isCharmHadron(p.pid())) {
          charm=true;
        }
        else if (p.pid() == PID::EMINUS || p.pid()==PID::MUON) {
          em.push_back(p);
        }
        else if (p.pid() == PID::EPLUS || p.pid()==PID::ANTIMUON) {
          ep.push_back(p);
        }
        else if (p.pid() == PID::NU_E  || p.pid()==PID::NU_MU) {
          nue.push_back(p);
        }
        else if (p.pid() == PID::NU_EBAR || p.pid()==PID::NU_MUBAR) {
          nueBar.push_back(p);
        }
        else if (PID::isBottomHadron(p.pid())) {
          findDecayProducts(p,em,ep,nue,nueBar,charm);
        }
        else if (!PID::isHadron(p.pid())) {
          findDecayProducts(p,em,ep,nue,nueBar,charm);
        }
      }
    }

    /// Perform the per-event analysis
    void analyze(const Event& event) {
      if(_edges.empty()) _edges = _p[0][0]->xEdges();
      // find and loop over Upslion(4S)
      for (const Particle& p : apply<UnstableParticles>(event, "UFS").particles()) {
        if (p.children().empty() || (p.children().size()==1 && p.children()[1].abspid()==p.abspid())) {
          continue;
        }
        // find decay products
        bool charm = false;
        Particles em,ep,nue,nueBar;
        findDecayProducts(p,em,ep,nue,nueBar,charm);
        if (!charm) continue;
        FourMomentum pl,pnu;
        unsigned int il = 0;
        if (em.size()==1 && nueBar.size()==1 && em[0].pid()+1==-nueBar[0].pid()) {
          pl  = em[0].momentum();
          pnu = nueBar[0].momentum();
          if (em[0].abspid()==13) il=1;
        }
        else if (ep.size()==1 && nue.size()==1 && nue[0].pid()==-ep[0].pid()+1) {
          pl  = ep[0].momentum();
          pnu = nue[0].momentum();
          if(ep[0].abspid()==13) il=1;
        }
        else {
          continue;
        }
        double q2 = (pl+pnu).mass2();
        vector<double> q2n(8);
        for (unsigned int ix=0; ix<8; ++ix) q2n[ix] = pow(q2,1+ix);
        if (q2<3.0) continue;
        double q2cut=3.0;
        for (unsigned int ibin=0; ibin<15; ++ibin) {
          if (q2>q2cut) {
            for(unsigned int ix=0; ix<8; ++ix) {
              _p[il][ix]->fill(_edges[ibin],q2n[ix]);
            }
          }
          else {
            break;
          }
          q2cut+=0.5;
        }
      }
    }


    /// Normalise histograms etc., after the run
    void finalize() {
      for (unsigned int il=0; il<2; ++il) {
        BinnedEstimatePtr<string> tmp[3];
        for (unsigned int ix=0;ix<3;++ix) book(tmp[ix],3*il+9+ix,1,1);
        for (unsigned int iy=0;iy<_p[il][0]->numBins();++iy) {
          const double N = _p[il][0]->bin(iy+1).effNumEntries();
          double q2  = _p[il][0]->bin(iy+1).mean(2), q4  = _p[il][1]->bin(iy+1).mean(2),
                 q6  = _p[il][2]->bin(iy+1).mean(2), q8  = _p[il][3]->bin(iy+1).mean(2),
                 q10 = _p[il][4]->bin(iy+1).mean(2), q12 = _p[il][5]->bin(iy+1).mean(2),
                 q14 = _p[il][6]->bin(iy+1).mean(2), q16 = _p[il][7]->bin(iy+1).mean(2);
          // q4 case
          if (N <= 9) continue; // q8 needs 8 events + Bessel correction
          double value = q4 - sqr(q2);
          double error = (-sqr(q4) + 4*sqr(q2)*(-sqr(q2) + 2*q4) - 4*q2*q6 + q8)/N;
          tmp[0]->bin(iy+1).set(value,sqrt(error));
          // q6 case
          if (N <= 13) continue; // q12 needs 12 events + Bessel correction
          value = q6 + q2*(2*sqr(q2) - 3*q4);
          error = (q12 - sqr(q6) - 6*q2*(q10 - 5*q4*q6) + 3*q4*(3*sqr(q4) - 2*q8) +
                   sqr(q2)*(-72*sqr(q4) + 36*sqr(q2)*(-sqr(q2) + 3*q4) - 48*q2*q6 + 21*q8))/N;
          tmp[1]->bin(iy+1).set(value,sqrt(error));
          // q8 case
          if (N <= 17) continue; // q16 needs 16 events + Bessel correction
          value = q8 + q2*(-3*q2*sqr(q2) + 6*q2*q4 - 4*q6);
          error = (q16 - 8*q6*(q10 - 2*q4*q6) - sqr(q8) - 8*q2*(q14 - 3*q4*(q10 - 4*q4*q6) - 6*q6*q8) +
                   4*sqr(q2)*(7*q12 - 28*sqr(q6) + 6*q2*(-3*q10 + 22*q4*q6) + 3*q4*(12*sqr(q4) - 11*q8) +
                                  3*sqr(q2)*(-12*sqr(q2)*(sqr(q2) - 4*q4) - 51*sqr(q4) - 28*q2*q6 + 13*q8)))/N;
          tmp[2]->bin(iy+1).set(value,sqrt(error));
        }
      }
    }

    /// @}


    /// @name Histograms
    /// @{
    BinnedProfilePtr<string> _p[2][8];
    vector<string> _edges;
    /// @}


  };


  RIVET_DECLARE_PLUGIN(BELLE_2021_I1917200);

}