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 Lambda_c+ -> Lambda Kbar0 K+
 10  class BABAR_2006_I722213 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BABAR_2006_I722213);
 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==4122);
 24      declare(ufs, "UFS");
 25      DecayedParticles LAMBDAC(ufs);
 26      LAMBDAC.addStable(PID::PI0);
 27      LAMBDAC.addStable(PID::K0S);
 28      LAMBDAC.addStable(PID::ETA);
 29      LAMBDAC.addStable(PID::LAMBDA);
 30      LAMBDAC.addStable(-PID::LAMBDA);
 31      declare(LAMBDAC, "LAMBDAC");
 32      // histograms
 33      for(unsigned int ix=0;ix<3;++ix)
 34	book(_h_mass[ix],1,1,1+ix);
 35      for(unsigned int ix=0;ix<2;++ix)
 36	book(_h_ctheta[ix],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   = { { PID::LAMBDA,1}, { 321,1}, { 310,1}};
 43      static const map<PdgId,unsigned int> & modeCC = { {-PID::LAMBDA,1}, {-321,1}, { 310,1}};
 44      DecayedParticles LAMBDAC = apply<DecayedParticles>(event, "LAMBDAC");
 45      // loop over particles
 46      for(unsigned int ix=0;ix<LAMBDAC.decaying().size();++ix) {
 47	int sign = 1;
 48	if (LAMBDAC.decaying()[ix].pid()>0 && LAMBDAC.modeMatches(ix,3,mode)) {
 49	  sign=1;
 50	}
 51	else if  (LAMBDAC.decaying()[ix].pid()<0 && LAMBDAC.modeMatches(ix,3,modeCC)) {
 52	  sign=-1;
 53	}
 54	else
 55	  continue;
 56	const Particle & lam = LAMBDAC.decayProducts()[ix].at( sign*PID::LAMBDA)[0];
 57	const Particle & K0  = LAMBDAC.decayProducts()[ix].at( 310)[0];
 58	const Particle & Kp  = LAMBDAC.decayProducts()[ix].at( sign*321)[0];
 59	double mLamK = (lam.momentum()+K0.momentum()).mass();
 60	_h_mass[0]->fill(mLamK);
 61	_h_mass[1]->fill((Kp .momentum()+K0.momentum()).mass()-Kp.mass()-K0.mass());
 62	_h_mass[2]->fill((lam.momentum()+Kp.momentum()).mass()-Kp.mass()-lam.mass());
 63	// take Xi(1690) to be any resonance in mass region
 64	if(mLamK<1.6725 || mLamK>1.6975) continue;
 65	if(LAMBDAC.decaying()[ix].children().size()!=2) continue;
 66	if(LAMBDAC.decaying()[ix].children()[0].abspid()!=321 &&
 67	   LAMBDAC.decaying()[ix].children()[1].abspid()!=321) continue;
 68	LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(LAMBDAC.decaying()[ix].momentum().betaVec());
 69	FourMomentum pbaryon1 = boost1.transform(LAMBDAC.decaying()[ix].momentum());
 70	FourMomentum pbaryon2 = boost1.transform(lam.momentum());
 71	LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(pbaryon1.betaVec());
 72	Vector3 axis = pbaryon1.p3().unit();
 73	FourMomentum pp = boost2.transform(pbaryon2);
 74	// calculate angle
 75	double cTheta = pp.p3().unit().dot(axis);
 76	_h_ctheta[0]->fill(cTheta);
 77	_h_ctheta[1]->fill(cTheta);
 78      }
 79    }
 80
 81
 82    /// Normalise histograms etc., after the run
 83    void finalize() {
 84      for(unsigned int ix=0;ix<3;++ix)
 85	normalize(_h_mass  [ix]);
 86      for(unsigned int ix=0;ix<2;++ix)
 87	normalize(_h_ctheta[ix]);
 88    }
 89
 90    /// @}
 91
 92
 93    /// @name Histograms
 94    /// @{
 95    Histo1DPtr _h_mass  [3];
 96    Histo1DPtr _h_ctheta[2];
 97    /// @}
 98
 99
100  };
101
102
103  RIVET_DECLARE_PLUGIN(BABAR_2006_I722213);
104
105}