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        Scatter3D matrix = refData<YODA::Scatter3D>(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          err[iy] =val[iy]*sqrt(sqr(eff.bin(iy+1).relErrAvg()) + sqr(_h_forward[ix]->bin(iy+1).relErrW()));
126        }
127        vector<double> val2(_h_forward[ix]->numBins(),0.),err2(_h_forward[ix]->numBins(),0.);
128        for (unsigned int iy=0;iy<_h_forward[ix]->numBins();++iy) {
129          for (unsigned int iz=0;iz<_h_forward[ix]->numBins();++iz) {
130            double corr  = matrix.points()[_h_forward[ix]->numBins()*iz+iy].z()/100.;
131            double ecorr = matrix.points()[_h_forward[ix]->numBins()*iz+iy].zErrAvg()/100.;
132            val2[iy] += corr*val[iz];
133            err2[iy] += sqr(ecorr*val[iz]) + sqr(corr*err[iz]);
134          }
135          err2[iy]  = val2[iy]*sqrt(err2[iy]/sqr(val2[iy]) +  sqr(9.78/771.58));
136        }
137        for (unsigned int ibin=0; ibin<_h_forward[ix]->numBins(); ++ibin) {
138          const double dy = sqrt(err[ibin]);
139          corrected->bin(ibin+1).set(val[ibin], dy);
140        }
141      }
142    }
143
144    /// @}
145
146
147    /// @name Histograms
148    /// @{
149    Histo1DPtr _h_direct [6];
150    Histo1DPtr _h_forward[6];
151    CounterPtr _nB;
152    /// @}
153
154
155  };
156
157
158  RIVET_DECLARE_PLUGIN(BELLE_2021_I1895149);
159
160}