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 -> phi phi K
 10  class BABAR_2011_I901433 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BABAR_2011_I901433);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22      UnstableParticles ufs = UnstableParticles(Cuts::abspid==511 or
 23						Cuts::abspid==521);
 24      declare(ufs, "UFS");
 25      DecayedParticles BB(ufs);
 26      BB.addStable(PID::PHI);
 27      BB.addStable(PID::K0S);
 28      declare(BB, "BB");
 29      // histograms
 30      for(unsigned int ix=0;ix<6;++ix)
 31	book(_h[ix],1,1,1+ix);
 32    }
 33
 34    /// Perform the per-event analysis
 35    void analyze(const Event& event) {
 36      DecayedParticles BB = apply<DecayedParticles>(event, "BB");
 37      static const map<PdgId,unsigned int> & mode1   = { { 333,2},{ 321,1}};
 38      static const map<PdgId,unsigned int> & mode1CC = { { 333,2},{-321,1}};
 39      static const map<PdgId,unsigned int> & mode2   = { { 333,2},{ 310,1}};
 40      for(unsigned int ix=0;ix<BB.decaying().size();++ix) {
 41	if(BB.modeMatches(ix,3,mode1) || BB.modeMatches(ix,3,mode1CC) ||
 42	   BB.modeMatches(ix,3,mode2)) {
 43	  // phi mesons
 44	  const Particles & phi = BB.decayProducts()[ix].at(333);
 45	  // bost to B rest frane
 46	  LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(BB.decaying()[ix].momentum().betaVec());
 47	  FourMomentum pPhiPhi= boost1.transform(phi[0].momentum()+phi[1].momentum());
 48	  double mPhiPhi = pPhiPhi.mass();
 49	  int iloc=-1;
 50	  if(mPhiPhi>2.94 && mPhiPhi<3.02)
 51	    iloc=0;
 52	  else if(mPhiPhi<2.85)
 53	    iloc=3;
 54	  else continue;
 55	  // cos theta phi phi
 56	  Vector3 axis1 = pPhiPhi.p3().unit();
 57	  LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(pPhiPhi.betaVec());
 58	  Vector3 axis2 = boost2.transform(boost1.transform(phi[0].momentum())).p3().unit();
 59	  _h[iloc+2]->fill(abs(axis1.dot(axis2)));
 60	  // now for the phi decays
 61	  Vector3 Trans[2];
 62	  bool foundPhi=true;
 63	  for(unsigned int ix=0;ix<2;++ix) {
 64	    if(phi[ix].children().size()!=2|| phi[ix].children()[0].pid()!=-phi[ix].children()[1].pid() ||
 65	       phi[ix].children()[0].abspid()!=321) {
 66	      foundPhi = false;
 67	      break;
 68	    }
 69	    Particle Km = phi[ix].children()[0];
 70	    Particle Kp = phi[ix].children()[1];
 71	    if(Kp.pid()<0) swap(Km,Kp);
 72	    FourMomentum pKp  = boost2.transform(boost1.transform(Kp.momentum()));
 73	    FourMomentum pPhi = boost2.transform(boost1.transform(phi[ix].momentum()));
 74	    LorentzTransform boost3 = LorentzTransform::mkFrameTransformFromBeta(pPhi.betaVec());
 75	    pKp = boost3.transform(pKp);
 76	    double cK = axis2.dot(pKp.p3().unit());
 77	    _h[iloc+1]->fill(cK);
 78	    Trans[ix] = pKp.p3() - cK*pKp.p3().mod()*axis2;
 79	  }
 80	  if(!foundPhi) continue;
 81	  double chi = atan2(Trans[0].cross(Trans[1]).dot(axis2),Trans[0].dot(Trans[1]));
 82	  _h[iloc]->fill(abs(chi));
 83	}
 84      }
 85    }
 86
 87
 88    /// Normalise histograms etc., after the run
 89    void finalize() {
 90      for(unsigned int ix=0;ix<6;++ix)
 91	normalize(_h[ix],1.,false);
 92    }
 93
 94    /// @}
 95
 96
 97    /// @name Histograms
 98    /// @{
 99    Histo1DPtr _h[6];
100    /// @}
101
102
103  };
104
105
106  RIVET_DECLARE_PLUGIN(BABAR_2011_I901433);
107
108}