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