1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief B -> c l nu moments
  9  class BELLE_2023_I2649712 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2023_I2649712);
 14
 15
 16    /// @name Analysis methods
 17    /// @{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21      // projections
 22      declare(UnstableParticles(Cuts::abspid==511 || Cuts::abspid==521),"UFS");
 23      // histos
 24      for (unsigned int ix=0; ix<4; ++ix) {
 25        book(_p[ix  ], 1+ix, 1, 1);
 26        book(_p[ix+4], "TMP/p_"+toString(ix),refData<YODA::BinnedEstimate<string>>(1+ix, 1, 1));
 27      }
 28    }
 29
 30    void findDecayProducts(Particle parent, Particles& em, Particles& ep,
 31                           Particles& nue, Particles& nueBar, bool& charm) {
 32      for (const Particle& p : parent.children()) {
 33        if (PID::isCharmHadron(p.pid())) {
 34          charm=true;
 35        }
 36        else if (p.pid() == PID::EMINUS || p.pid()==PID::MUON) {
 37          em.push_back(p);
 38        }
 39        else if (p.pid() == PID::EPLUS || p.pid()==PID::ANTIMUON) {
 40          ep.push_back(p);
 41        }
 42        else if (p.pid() == PID::NU_E  || p.pid()==PID::NU_MU) {
 43          nue.push_back(p);
 44        }
 45        else if (p.pid() == PID::NU_EBAR || p.pid()==PID::NU_MUBAR) {
 46          nueBar.push_back(p);
 47        }
 48        else if (PID::isBottomHadron(p.pid())) {
 49          findDecayProducts(p,em,ep,nue,nueBar,charm);
 50        }
 51        else if (!PID::isHadron(p.pid())) {
 52          findDecayProducts(p,em,ep,nue,nueBar,charm);
 53        }
 54      }
 55    }
 56
 57    /// Perform the per-event analysis
 58    void analyze(const Event& event) {
 59      if(_edges.empty()) _edges=_p[0]->xEdges();
 60      // find and loop over Upslion(4S)
 61      for (const Particle& p : apply<UnstableParticles>(event, "UFS").particles()) {
 62        if (p.children().empty() || (p.children().size()==1 && p.children()[1].abspid()==p.abspid())) {
 63          continue;
 64        }
 65        // find decay products
 66        bool charm = false;
 67        Particles em,ep,nue,nueBar;
 68        findDecayProducts(p,em,ep,nue,nueBar,charm);
 69        if (!charm) continue;
 70        FourMomentum pl,pnu;
 71        if (em.size()==1 && nueBar.size()==1 && em[0].pid()+1==-nueBar[0].pid()) {
 72          pl  = em[0].mom();
 73          pnu = nueBar[0].mom();
 74        }
 75        else if(ep.size()==1 && nue.size()==1 && nue[0].pid()==-ep[0].pid()+1) {
 76          pl  = ep[0].mom();
 77          pnu = nue[0].mom();
 78       	}
 79        else {
 80          continue;
 81        }
 82        double q2 = (pl+pnu).mass2();
 83        vector<double> q2n(8);
 84        for (unsigned int ix=0;ix<8;++ix) q2n[ix] = pow(q2,1+ix);
 85        if (q2<1.5) continue;
 86        double q2cut=1.5;
 87        for (unsigned int ibin=0; ibin<15; ++ibin) {
 88          if (q2>q2cut) {
 89            for (unsigned int ix=0; ix<8; ++ix) {
 90              _p[ix]->fill(_edges[ibin],q2n[ix]);
 91            }
 92          }
 93          else break;
 94          q2cut+=0.5;
 95        }
 96      }
 97    }
 98
 99
100    /// Normalise histograms etc., after the run
101    void finalize() {
102      BinnedEstimatePtr<string> tmp[3];
103      for (unsigned int ix=0; ix<3;++ix) book(tmp[ix],5+ix,1,1);
104      for (unsigned int iy=0; iy<_p[0]->numBins();++iy) {
105        double q2  = _p[0]->bin(iy+1).mean(2)  ,  q4  = _p[1]->bin(iy+1).mean(2),
106               q6  = _p[2]->bin(iy+1).mean(2)  ,  q8  = _p[3]->bin(iy+1).mean(2),
107              q10  = _p[4]->bin(iy+1).mean(2)  ,  q12 = _p[5]->bin(iy+1).mean(2),
108          q14  = _p[6]->bin(iy+1).mean(2)  ,  q16 = _p[7]->bin(iy+1).mean(2);
109        // q4 case
110        double N = _p[0]->bin(iy+1).effNumEntries();
111        double value = q4 - sqr(q2);
112        double error = (-sqr(q4) + 4*sqr(q2)*(-sqr(q2) + 2*q4) - 4*q2*q6 + q8)/N;
113        tmp[0]->bin(iy+1).set(value,sqrt(error));
114        // q6 case
115        value = q6 + q2*(2*sqr(q2) - 3*q4);
116        error = (q12 - sqr(q6) - 6*q2*(q10 - 5*q4*q6) + 3*q4*(3*sqr(q4) - 2*q8) + 
117                 sqr(q2)*(-72*sqr(q4) + 36*sqr(q2)*(-sqr(q2) + 3*q4) - 48*q2*q6 + 21*q8))/N;
118        tmp[1]->bin(iy+1).set(value,sqrt(error));
119        // q8 case
120        value = q8 + q2*(-3*q2*sqr(q2) + 6*q2*q4 - 4*q6);
121        error = (q16 - 8*q6*(q10 - 2*q4*q6) - sqr(q8) - 8*q2*(q14 - 3*q4*(q10 - 4*q4*q6) - 6*q6*q8) + 
122                 4*sqr(q2)*(7*q12 - 28*sqr(q6) + 6*q2*(-3*q10 + 22*q4*q6) + 3*q4*(12*sqr(q4) - 11*q8) + 
123                            3*sqr(q2)*(-12*sqr(q2)*(sqr(q2) - 4*q4) - 51*sqr(q4) - 28*q2*q6 + 13*q8)))/N;
124        tmp[2]->bin(iy+1).set(value,sqrt(error));
125      }
126    }
127    /// @}
128
129
130    /// @name Histograms
131    /// @{
132    BinnedProfilePtr<string> _p[8];
133    vector<string> _edges;
134    /// @}
135
136
137  };
138
139
140  RIVET_DECLARE_PLUGIN(BELLE_2023_I2649712);
141
142}