1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief B -> pi,rho l+ nu
  9  class BABAR_2011_I855306 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(BABAR_2011_I855306);
 14
 15
 16    /// @name Analysis methods
 17    /// @{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21
 22      // Initialise and register projections
 23      declare(UnstableParticles(Cuts::abspid==PID::BPLUS ||
 24				Cuts::abspid==PID::B0   ), "UFS");
 25      for(unsigned int ix=0;ix<2;++ix) {
 26	book(_c[ix],"TMP/c_"+toString(ix+1));
 27	for(unsigned int iy=0;iy<3;++iy) {
 28	  book(_h[ix][iy],ix+1,1,iy+1);
 29	}
 30      }
 31    }
 32
 33    // Calculate the Q2 using mother and daughter charged lepton
 34    double q2(const Particle& B, int mesonID) {
 35      FourMomentum q = B.mom() - select(B.children(), Cuts::pid==mesonID)[0];
 36      return q*q;
 37    }
 38
 39    // Check for explicit decay into pdgids
 40    bool isSemileptonicDecay(const Particle& mother, vector<int> ids) {
 41      // Trivial check to ignore any other decays but the one in question modulo photons
 42      const Particles children = mother.children(Cuts::pid!=PID::PHOTON);
 43      if (children.size()!=ids.size()) return false;
 44      // Check for the explicit decay
 45      return all(ids, [&](int i){return count(children, hasPID(i))==1;});
 46    }
 47
 48    /// Perform the per-event analysis
 49    void analyze(const Event& event) {
 50      // Get B+ Mesons
 51      for(const Particle& p : apply<UnstableParticles>(event, "UFS").particles()) {
 52	// B0 modes
 53	if(p.abspid()==511) {
 54	  _c[0]->fill();
 55	  if (isSemileptonicDecay(p, {PID::PIMINUS, PID::POSITRON, PID::NU_E}) ||
 56	      isSemileptonicDecay(p, {PID::PIMINUS, PID::ANTIMUON, PID::NU_MU}) ) {
 57	    double qq = q2(p,PID::PIMINUS);
 58            _h[0][0]->fill(qq);
 59            _h[0][2]->fill(qq);
 60	  }
 61	  else if(isSemileptonicDecay(p, {PID::PIPLUS , PID::ELECTRON, PID::NU_EBAR}) ||
 62		  isSemileptonicDecay(p, {PID::PIPLUS , PID::MUON    , PID::NU_MUBAR})) {
 63	    double qq = q2(p,PID::PIPLUS);
 64            _h[0][0]->fill(qq);
 65            _h[0][2]->fill(qq);
 66	  }
 67	  else if (isSemileptonicDecay(p, {PID::RHOMINUS, PID::POSITRON, PID::NU_E}) ||
 68		   isSemileptonicDecay(p, {PID::RHOMINUS, PID::ANTIMUON, PID::NU_MU})) {
 69	    double qq = q2(p,PID::RHOMINUS);
 70            _h[1][0]->fill(qq);
 71            _h[1][2]->fill(qq);
 72	  }
 73	  else if( isSemileptonicDecay(p, {PID::RHOPLUS, PID::ELECTRON, PID::NU_EBAR}) ||
 74		   isSemileptonicDecay(p, {PID::RHOPLUS, PID::MUON    , PID::NU_MUBAR})) {
 75	    double qq = q2(p,PID::RHOPLUS);
 76            _h[1][0]->fill(qq);
 77            _h[1][2]->fill(qq);
 78	  }
 79	}
 80	// B+ modes
 81	else {
 82	  _c[1]->fill();
 83	  if (isSemileptonicDecay(p, {PID::PI0, PID::POSITRON, PID::NU_E}) ||
 84	      isSemileptonicDecay(p, {PID::PI0, PID::ANTIMUON, PID::NU_MU}) ||
 85	      isSemileptonicDecay(p, {PID::PI0, PID::ELECTRON, PID::NU_EBAR}) ||
 86	      isSemileptonicDecay(p, {PID::PI0, PID::MUON    , PID::NU_MUBAR})) {
 87            _h[0][1]->fill(q2(p,PID::PI0));
 88	  }
 89	  else if (isSemileptonicDecay(p, {PID::RHO0, PID::POSITRON, PID::NU_E}) ||
 90		   isSemileptonicDecay(p, {PID::RHO0, PID::ANTIMUON, PID::NU_MU}) ||
 91		   isSemileptonicDecay(p, {PID::RHO0, PID::ELECTRON, PID::NU_EBAR}) ||
 92		   isSemileptonicDecay(p, {PID::RHO0, PID::MUON    , PID::NU_MUBAR})) {
 93            _h[1][1]->fill(q2(p,PID::RHO0));
 94	  }
 95	}
 96      }
 97    }
 98
 99
100    /// Normalise histograms etc., after the run
101    void finalize() {
102      // Br in units 10^-4 and average over e/mu modes
103      for(unsigned int ix=0;ix<2;++ix) {
104	for(unsigned int iy=0;iy<3;++iy) {
105	  if(iy==1) scale(_h[ix][iy], 0.5*1e4/ *_c[1]);
106	  else      scale(_h[ix][iy], 0.5*1e4/ *_c[0]);
107	}
108      }
109    }
110
111    /// @}
112
113
114    /// @name Histograms
115    /// @{
116    Histo1DPtr _h[2][3];
117    CounterPtr _c[2];
118    /// @}
119
120
121  };
122
123
124  RIVET_DECLARE_PLUGIN(BABAR_2011_I855306);
125
126}