1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief B -> X_u l nu
  9  class BELLE_2021_I1895149 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2021_I1895149);
 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(),"UFS");
 23      // histograms
 24      for(unsigned int ix=0;ix<6;++ix) {
 25	book(_h_direct[ix],1+ix,1,1);
 26	book(_h_forward[ix],"TMP/h_"+toString(ix+1),refData(7+ix,1,1));
 27      }
 28      book(_nB,"/TMP/nB");
 29    }
 30
 31    void findDecayProducts(Particle parent, Particles & em, Particles & ep,
 32			   Particles & nue, Particles & nueBar, bool & charm) {
 33      for(const Particle & p : parent.children()) {
 34	if(PID::isCharmHadron(p.pid())) {
 35	  charm=true;
 36	}
 37	else if(p.pid() == PID::EMINUS) {
 38	  em.push_back(p);
 39	}
 40	else if(p.pid() == PID::EPLUS) {
 41	  ep.push_back(p);
 42	}
 43	else if(p.pid() == PID::NU_E || p.pid()==PID::NU_MU) {
 44	  nue.push_back(p);
 45	}
 46	else if(p.pid() == PID::NU_EBAR || p.pid()==PID::NU_MUBAR) {
 47	  nueBar.push_back(p);
 48	}
 49	else if(PID::isBottomHadron(p.pid())) {
 50	  findDecayProducts(p,em,ep,nue,nueBar,charm);
 51	}
 52	else if(!PID::isHadron(p.pid())) {
 53	  findDecayProducts(p,em,ep,nue,nueBar,charm);
 54	}
 55      }
 56    }
 57
 58    /// Perform the per-event analysis
 59    void analyze(const Event& event) {
 60      // find and loop over Upslion(4S)
 61      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 62      for (const Particle& p : ufs.particles(Cuts::pid==300553)) {
 63      	for(const Particle & p2 : p.children()) {
 64      	  if(p2.abspid()!=511 && p2.abspid()!=521) continue;
 65	  _nB->fill();
 66	  bool charm = false;
 67	  Particles em,ep,nue,nueBar;
 68	  findDecayProducts(p2,em,ep,nue,nueBar,charm);
 69	  if(charm) continue;
 70	  FourMomentum pl,pnu;
 71	  if(em.size()==1 && nueBar.size()==1) {
 72	    pl  = em[0].momentum();
 73	    pnu = nueBar[0].momentum();
 74	  }
 75	  else if(ep.size()==1 && nue.size()==1) {
 76	    pl  = ep[0].momentum();
 77	    pnu = nue[0].momentum();
 78	  }
 79	  else
 80	    continue;
 81      	  LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(p2.momentum().betaVec());
 82	  pl  = boost.transform(pl );
 83	  pnu = boost.transform(pnu);
 84	  FourMomentum pB = boost.transform(p2.momentum());
 85	  FourMomentum q = pl+pnu;
 86	  FourMomentum pX = pB-q;
 87	  double p3 = pX.p();
 88	  _h_forward[0]->fill(pl.E());
 89	  _h_forward[1]->fill(q.mass2());
 90	  _h_forward[2]->fill(pX.mass());
 91	  _h_forward[3]->fill(pX.mass2());
 92	  _h_forward[4]->fill(pX.E()-p3);
 93	  _h_forward[5]->fill(pX.E()+p3);
 94	  if(pl.E()>1) {
 95	    _h_direct[0]->fill(pl.E());
 96	    _h_direct[1]->fill(q.mass2());
 97	    _h_direct[2]->fill(pX.mass());
 98	    _h_direct[3]->fill(pX.mass2());
 99	    _h_direct[4]->fill(pX.E()-p3);
100	    _h_direct[5]->fill(pX.E()+p3);
101	  }
102	}
103      }
104    }
105
106
107    /// Normalise histograms etc., after the run
108    void finalize() {
109      for(unsigned int ix=0;ix<6;++ix) {
110        // unfolded dist, scale by 1/2 /no of B's (2 as using e and mu modes)
111        scale(_h_direct[ix], 0.5/ *_nB);
112        // forward folding scale to BELLE no of B's
113        scale(_h_forward[ix], 2.*771.58e6/ *_nB);
114        // get the efficiency product and divide by it
115        unsigned int iloc = ix<2 ? 3+ix : (ix<4 ? ix-1 : ix+1);
116        Estimate1D eff    = refData<YODA::Estimate1D>(iloc+24,1,1);
117        Estimate2D matrix = refData<YODA::Estimate2D>(  19+ix,1,1);
118        // scatter for the result
119        Estimate1DPtr corrected;
120        book(corrected,ix+7,1,1);
121        vector<double> val(_h_forward[ix]->numBins(),0.),err(_h_forward[ix]->numBins(),0.);
122        // first divide by eff
123        for (unsigned int iy=0;iy<_h_forward[ix]->numBins();++iy) {
124          val[iy] = _h_forward[ix]->bin(iy+1).sumW()/eff.bin(iy+1).val();
125          double relE = eff.bin(iy+1).totalErrAvg()/eff.bin(iy+1).val();
126          err[iy] =val[iy]*sqrt(sqr(relE) + sqr(_h_forward[ix]->bin(iy+1).relErrW()));
127        }
128        vector<double> val2(_h_forward[ix]->numBins(),0.),err2(_h_forward[ix]->numBins(),0.);
129        for (unsigned int iy=0;iy<_h_forward[ix]->numBins();++iy) {
130          for (unsigned int iz=0;iz<_h_forward[ix]->numBins();++iz) {
131            double corr  = matrix.bin((_h_forward[ix]->numBins()+2)*(iz+1)+iy+1).val()/100.;
132            double ecorr = matrix.bin((_h_forward[ix]->numBins()+2)*(iz+1)+iy+1).totalErrAvg()/100.;
133            val2[iy] += corr*val[iz];
134            err2[iy] += sqr(ecorr*val[iz]) + sqr(corr*err[iz]);
135          }
136          err2[iy]  = val2[iy]*sqrt(err2[iy]/sqr(val2[iy]) +  sqr(9.78/771.58));
137        }
138        for (unsigned int ibin=0; ibin<_h_forward[ix]->numBins(); ++ibin) {
139          const double dy = sqrt(err[ibin]);
140          corrected->bin(ibin+1).set(val[ibin], dy);
141        }
142      }
143    }
144
145    /// @}
146
147
148    /// @name Histograms
149    /// @{
150    Histo1DPtr _h_direct [6];
151    Histo1DPtr _h_forward[6];
152    CounterPtr _nB;
153    /// @}
154
155
156  };
157
158
159  RIVET_DECLARE_PLUGIN(BELLE_2021_I1895149);
160
161}