1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief B -> Xs l+l-
  9  class BABAR_2013_I1272843 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(BABAR_2013_I1272843);
 14
 15
 16    /// @name Analysis methods
 17    /// @{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21      // Initialise and register projections
 22      declare(UnstableParticles(Cuts::abspid==521 || Cuts::abspid==511), "UFS");
 23      // Book histograms
 24      for(unsigned int ix=0;ix<2;++ix) book(_p_A[ix],1,1+ix,4);
 25      for(unsigned int iy=0;iy<3;++iy) {
 26	book(_h_mX[iy],2,1,1+iy);
 27	for(unsigned int ix=0;ix<2;++ix)
 28	  book(_h_q2[ix][iy],1,1+ix,1+iy);
 29      }
 30      book(_nBottom, "TMP/BottomCounter");
 31    }
 32
 33    void findDecayProducts(bool & charm, const Particle& mother,
 34                           unsigned int& nK0, unsigned int& nKp,
 35			   unsigned int& nKm, Particles & lp, Particles & lm) {
 36      for (const Particle & p : mother.children()) {
 37        int id = p.pid();
 38	if(PID::isCharmHadron(p.pid())) charm = true;
 39	else if ( id == PID::POSITRON || id == PID::ANTIMUON) lp.push_back(p);
 40	else if ( id == PID::ELECTRON || id == PID::MUON    ) lm.push_back(p);
 41        else if ( id == PID::KPLUS )    ++nKp;
 42        else if (id == PID::KMINUS )    ++nKm;
 43        else if (id == PID::K0S)        ++nK0;
 44        else if (id == PID::PI0 || id == PID::PIPLUS || id == PID::PIMINUS) {
 45	  continue;
 46        }
 47        else if ( !p.children().empty() ) {
 48          findDecayProducts(charm, p, nK0, nKp, nKm,lp,lm);
 49        }
 50      }
 51    }
 52
 53    /// Perform the per-event analysis
 54    void analyze(const Event& event) {
 55      // Loop over bottoms
 56      for (const Particle& bottom : apply<UnstableParticles>(event, "UFS").particles()) {
 57       	// remove mixing entries etc
 58	if(bottom.children()[0].abspid()==bottom.abspid()) continue;
 59	_nBottom->fill();
 60	bool charm = false;
 61	Particles lp,lm;
 62	unsigned int nK0(0),nKp(0),nKm(0);
 63	findDecayProducts(charm,bottom, nK0, nKp, nKm,lp,lm);
 64	if(charm) continue;
 65        unsigned int nk = nKp-nKm+nK0;
 66	if( nk % 2 == 0) continue;
 67	if (lp.size()!=1 || lm.size()!=1 || lp[0].pid()!=-lm[0].pid()) continue;
 68	FourMomentum pl = (lp[0].momentum()+lm[0].momentum());
 69	FourMomentum pX=bottom.momentum()-pl;
 70	double q2 = pl.mass2();
 71	double mX = pX.mass();
 72	if(lm[0].pid()==11) {
 73	  _h_q2[0][0]->fill(q2);
 74	  _h_q2[1][0]->fill(q2);
 75	  _h_mX[0]->fill(mX);
 76	}
 77	else {
 78	  _h_q2[0][1]->fill(q2);
 79	  _h_q2[1][1]->fill(q2);
 80	  _h_mX[1]->fill(mX);
 81	}
 82	_h_q2[0][2]->fill(q2);
 83	_h_q2[1][2]->fill(q2);
 84	_h_mX[2]->fill(mX);
 85	// A_CP doesn't include charmonium regions
 86	if(q2<6.8 || (q2>10.1 && q2<12.9) || q2>14.2) {
 87	  double wgt = bottom.pid()>0 ? -1 : 1;
 88	  _p_A[0]->fill(q2,wgt);
 89	  _p_A[1]->fill(q2,wgt);
 90	}
 91      }
 92    }
 93
 94
 95    /// Normalise histograms etc., after the run
 96    void finalize() {
 97      for(unsigned int iy=0;iy<3;++iy) {
 98	scale(_h_mX[iy],1e6/ *_nBottom);
 99	for(unsigned int ix=0;ix<2;++ix)
100	  scale(_h_q2[ix][iy],1e6/ *_nBottom);
101      }
102    }
103    /// @}
104
105
106    /// @name Histograms
107    /// @{
108    Histo1DPtr _h_q2[2][3],_h_mX[3];
109    Profile1DPtr _p_A[2];
110    CounterPtr _nBottom;
111    /// @}
112
113
114  };
115
116
117  RIVET_DECLARE_PLUGIN(BABAR_2013_I1272843);
118
119}