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 D_s+ -> K+K- e+ nu_e
 10  class BABAR_2008_I790461 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BABAR_2008_I790461);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22
 23      // Initialise and register projections
 24      UnstableParticles ufs = UnstableParticles(Cuts::pid==431);
 25      declare(ufs, "UFS");
 26      DecayedParticles DS(ufs);
 27      DS.addStable(PID::PI0);
 28      DS.addStable(PID::K0S);
 29      DS.addStable(PID::ETA);
 30      DS.addStable(PID::ETAPRIME);
 31      declare(DS, "DS");
 32      
 33      // Book histograms
 34      book(_h[0],1,1,1);
 35      for(unsigned int ix=0;ix<4;++ix)
 36	book(_h[ix+1],2,1,1+ix);
 37    }
 38
 39
 40    /// Perform the per-event analysis
 41    void analyze(const Event& event) {
 42      static const map<PdgId,unsigned int> & mode = { { 321,1}, {-321,1}, {-11,1}, { 12,1}};
 43      DecayedParticles DS = apply<DecayedParticles>(event, "DS");
 44      // loop over particles
 45      for(unsigned int ix=0;ix<DS.decaying().size();++ix) {
 46	if(!DS.modeMatches(ix,4,mode)) continue;
 47	const Particle & Kp = DS.decayProducts()[ix].at( 321)[0];
 48	const Particle & Km = DS.decayProducts()[ix].at(-321)[0];
 49	const Particle & ep = DS.decayProducts()[ix].at( -11)[0];
 50	const Particle & nue= DS.decayProducts()[ix].at(  12)[0];
 51	FourMomentum pPhi = Kp.momentum()+Km.momentum(); 
 52	_h[0]->fill(pPhi.mass());
 53	FourMomentum qq = DS.decaying()[ix].momentum()-pPhi;
 54	_h[1]->fill(qq.mass2());
 55	// boost momenta to DS rest frame
 56	LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(DS.decaying()[ix].momentum().betaVec());
 57	FourMomentum pPHI = boost.transform(pPhi);
 58	Matrix3 ptoz(-pPHI.p3().unit(), Vector3(0,0,1));
 59	boost.preMult(ptoz);
 60	// the momenta in frane to W along z
 61	FourMomentum pD  = boost.transform(DS.decaying()[ix].momentum());
 62	FourMomentum pKp = boost.transform(Kp .momentum());
 63	FourMomentum pKm = boost.transform(Km .momentum());
 64	FourMomentum pe  = boost.transform(ep .momentum());
 65	FourMomentum pnu = boost.transform(nue.momentum());
 66	pPhi = pKp+pKm;
 67	qq = pD-pPhi;
 68	LorentzTransform boostK = LorentzTransform::mkFrameTransformFromBeta(pPhi.betaVec());
 69	Vector3 axisK = boostK.transform(pKp).p3().unit();
 70	_h[3]->fill(axisK.dot(pPhi.p3().unit()));
 71	LorentzTransform boostW = LorentzTransform::mkFrameTransformFromBeta(  qq.betaVec());
 72	Vector3 axisE = boostW.transform(pe).p3().unit();
 73	_h[2]->fill(axisE.dot(qq.p3().unit()));
 74	axisK.setZ(0.);
 75	axisE.setZ(0.);
 76	double chi = atan2(axisE.cross(axisK).dot(qq.p3().unit()), axisE.dot(axisK));
 77	_h[4]->fill(chi);
 78      }
 79    }
 80
 81
 82    /// Normalise histograms etc., after the run
 83    void finalize() {
 84      for(unsigned int ix=0;ix<5;++ix)
 85	normalize(_h[ix]);
 86    }
 87
 88    /// @}
 89
 90
 91    /// @name Histograms
 92    /// @{
 93    Histo1DPtr _h[5];
 94    /// @}
 95
 96
 97  };
 98
 99
100  RIVET_DECLARE_PLUGIN(BABAR_2008_I790461);
101
102}