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 pi pi
 10  class BABAR_2008_I782405 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BABAR_2008_I782405);
 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==521);
 24      declare(ufs, "UFS");
 25      DecayedParticles BP(ufs);
 26      BP.addStable(PID::K0S);
 27      declare(BP, "BP");
 28      // histos
 29      for(unsigned int iy=0;iy<3;++iy) {
 30	book(_h_mass[iy],1,1,1+iy);
 31	for(unsigned int ix=0;ix<2;++ix)
 32	  book(_h_mass2[ix][iy],2,1+iy,1+ix);
 33      }
 34    }
 35
 36
 37    /// Perform the per-event analysis
 38    void analyze(const Event& event) {
 39      static const map<PdgId,unsigned int> & mode   = { { 321,1}, { 211,1}, {-211,1}};
 40      static const map<PdgId,unsigned int> & modeCC = { {-321,1}, { 211,1}, {-211,1}};
 41      DecayedParticles BP = apply<DecayedParticles>(event, "BP");
 42      // loop over particles
 43      for(unsigned int ix=0;ix<BP.decaying().size();++ix) {
 44      	int sign = BP.decaying()[ix].pid()>0 ? 1 : -1;
 45       	if((sign>0 and BP.modeMatches(ix,3,mode  )) ||
 46	   (sign<0 and BP.modeMatches(ix,3,modeCC))) {
 47	  LorentzTransform boost =
 48	    LorentzTransform::mkFrameTransformFromBeta(BP.decaying()[ix]. momentum().betaVec());
 49	  // momenta
 50	  FourMomentum pip  = boost.transform(BP.decayProducts()[ix].at( 211*sign)[0].momentum());
 51	  FourMomentum pim  = boost.transform(BP.decayProducts()[ix].at(-211*sign)[0].momentum());
 52	  FourMomentum Kp   = boost.transform(BP.decayProducts()[ix].at( 321*sign)[0].momentum());
 53	  // mass distributions
 54	  double mKpi = (Kp+pim).mass();
 55	  double mpipi = (pip+pim).mass();
 56	  // vetos
 57	  if(mKpi>1.756 && mKpi<1.931) continue;
 58	  if(mpipi>1.660 && mpipi<1.800) continue;
 59	  if(mpipi>3.019 && mpipi<3.179) continue;
 60	  if(mpipi>3.627 && mpipi<3.747) continue;
 61	  // mass distributions
 62	  if(mpipi>2.) _h_mass[0]->fill(mKpi);
 63	  if(mKpi>2.) {
 64	    _h_mass[1]->fill(mpipi);
 65	    _h_mass[2]->fill(mpipi);
 66	  }
 67	  FourMomentum ppipi = pim+pip;
 68	  LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(ppipi.betaVec());
 69	  double cTheta = boost2.transform(pim).p3().unit().dot(Kp.p3().unit());
 70	  _h_mass2[(sign+1)/2][0]->fill(mpipi);
 71	  if(cTheta>0)
 72	    _h_mass2[(sign+1)/2][1]->fill(mpipi);
 73	  else
 74	    _h_mass2[(sign+1)/2][2]->fill(mpipi);
 75	}
 76      }
 77    }
 78
 79
 80    /// Normalise histograms etc., after the run
 81    void finalize() {
 82      for(unsigned int iy=0;iy<3;++iy) {
 83	normalize(_h_mass[iy],1.,false);
 84	for(unsigned int ix=0;ix<2;++ix)
 85	  normalize(_h_mass2[ix][iy],1.,false);
 86      }
 87    }
 88
 89    /// @}
 90
 91
 92    /// @name Histograms
 93    /// @{
 94    Histo1DPtr _h_mass[3],_h_mass2[2][3];
 95    /// @}
 96
 97
 98  };
 99
100
101  RIVET_DECLARE_PLUGIN(BABAR_2008_I782405);
102
103}