1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4#include "Rivet/Projections/DecayedParticles.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief B -> K* l+ l-
 10  class BELLE_2019_I1728173 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2019_I1728173);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22      // Initialise and register projections
 23      UnstableParticles ufs = UnstableParticles(Cuts::abspid==511 or
 24						Cuts::abspid==521);
 25      declare(ufs, "UFS");
 26      DecayedParticles BB(ufs);
 27      BB.addStable(   443);
 28      BB.addStable(100443);
 29      BB.addStable( 313);
 30      BB.addStable( 323);
 31      BB.addStable(-313);
 32      BB.addStable(-323);
 33      declare(BB, "BB");
 34      // book histograms
 35      for(unsigned int ix=0;ix<4;++ix)
 36      	for(unsigned int iy=0;iy<3;++iy) {
 37      	  book(_h_br[ix][iy],"TMP/h_br_"+toString(ix+1)+"_"+toString(iy+1),refData(1,1+iy,1));
 38      	}
 39      for(unsigned int ix=0;ix<2;++ix)
 40      	book(_c[ix],"TMP/nB_"+toString(ix+1));
 41    }
 42
 43
 44    /// Perform the per-event analysis
 45    void analyze(const Event& event) {
 46      static const map<PdgId,unsigned int> & mode1   = { { 323,1},{ 13,1}, {-13,1}};
 47      static const map<PdgId,unsigned int> & mode1CC = { {-323,1},{ 13,1}, {-13,1}};
 48      static const map<PdgId,unsigned int> & mode2   = { { 313,1},{ 13,1}, {-13,1}};
 49      static const map<PdgId,unsigned int> & mode2CC = { {-313,1},{ 13,1}, {-13,1}};
 50      static const map<PdgId,unsigned int> & mode3   = { { 323,1},{ 11,1}, {-11,1}};
 51      static const map<PdgId,unsigned int> & mode3CC = { {-323,1},{ 11,1}, {-11,1}};
 52      static const map<PdgId,unsigned int> & mode4   = { { 313,1},{ 11,1}, {-11,1}};
 53      static const map<PdgId,unsigned int> & mode4CC = { {-313,1},{ 11,1}, {-11,1}};
 54      DecayedParticles BB = apply<DecayedParticles>(event, "BB");
 55      // loop over particles
 56      for(unsigned int ix=0;ix<BB.decaying().size();++ix) {
 57      	if(BB.decaying()[ix].abspid()==521) _c[0]->fill();
 58      	else                                _c[1]->fill();
 59      	int imode=0;
 60      	if ((BB.decaying()[ix].pid()>0 && BB.modeMatches(ix,3,mode1)) ||
 61	    (BB.decaying()[ix].pid()<0 && BB.modeMatches(ix,3,mode1CC)))       imode=0;
 62	else if ((BB.decaying()[ix].pid()>0 && BB.modeMatches(ix,3,mode2)) ||
 63		 (BB.decaying()[ix].pid()<0 && BB.modeMatches(ix,3,mode2CC)))  imode=1;
 64      	else if ((BB.decaying()[ix].pid()>0 && BB.modeMatches(ix,3,mode3)) ||
 65		 (BB.decaying()[ix].pid()<0 && BB.modeMatches(ix,3,mode3CC)))  imode=2;
 66      	else if ((BB.decaying()[ix].pid()>0 && BB.modeMatches(ix,3,mode4)) ||
 67		 (BB.decaying()[ix].pid()<0 && BB.modeMatches(ix,3,mode4CC)))  imode=3;
 68      	else continue;
 69      	int il = imode<2 ? 13 : 11;
 70      	const Particle & lp = BB.decayProducts()[ix].at(-il)[0];
 71      	const Particle & lm = BB.decayProducts()[ix].at( il)[0];
 72      	double qq = (lp.momentum()+lm.momentum()).mass2();
 73      	for(unsigned int iy=0;iy<3;++iy) {
 74      	  _h_br[imode][iy]->fill(qq);
 75      	}
 76      }
 77    }
 78
 79
 80    /// Normalise histograms etc., after the run
 81    void finalize() {
 82      // ratio of lifetimes
 83      double rLife = 1.078;
 84      // normalize BR plots
 85      for(unsigned int ix=0;ix<4;++ix) {
 86      	for(unsigned int iy=0;iy<3;++iy) {
 87      	  if(ix%2==0) {
 88      	     scale(_h_br[ix][iy],1e7/ *_c[0]);
 89      	  }
 90      	  else {
 91      	     scale(_h_br[ix][iy],1e7*rLife/ *_c[1]);
 92      	  }
 93      	}
 94      }
 95      // RK*
 96      for (unsigned int ix=0;ix<3;++ix) {
 97      	Estimate1DPtr RK;
 98      	book(RK,1,1+ix,3);
 99      	divide(_h_br[0][ix],_h_br[2][ix],RK);
100      	book(RK,1,1+ix,2);
101      	divide(_h_br[1][ix],_h_br[3][ix],RK);
102      	book(RK,1,1+ix,1);
103        size_t iloc=1;
104      	for (size_t ibin=0; ibin < _h_br[1][ix]->numBins(); ++ibin) {
105          while(_h_br[0][ix]->bin(iloc).isMasked()) ++iloc;
106      	  double num     = _h_br[0][ix]->bin(iloc).sumW()   +_h_br[1][ix]->bin(iloc).sumW();
107      	  double numErr2 = sqr(_h_br[0][ix]->bin(iloc).errW())+sqr(_h_br[1][ix]->bin(iloc).errW());
108      	  double den     = _h_br[2][ix]->bin(iloc).sumW()   +_h_br[3][ix]->bin(iloc).sumW();
109      	  double denErr2 = sqr(_h_br[2][ix]->bin(iloc).errW())+sqr(_h_br[3][ix]->bin(iloc).errW());
110      	  double val(0.),err(0.);
111      	  if(num>0. && den>0.) {
112      	    val = num/den;
113      	    err = val*(numErr2/sqr(num)+denErr2/sqr(den));
114      	  }
115      	  RK->bin(iloc).set(val, err);
116          ++iloc;
117      	}
118      }
119    }
120
121    /// @}
122
123
124    /// @name Histograms
125    /// @{
126    CounterPtr _c[2];
127    Histo1DPtr _h_br[4][3];
128    /// @}
129
130  };
131
132
133  RIVET_DECLARE_PLUGIN(BELLE_2019_I1728173);
134
135}