1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief D_s1 decay angles
  9  class BELLE_2008_I762013 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2008_I762013);
 14
 15
 16    /// @name Analysis methods
 17    /// @{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21      // projections
 22      declare(UnstableParticles(), "UFS");
 23      // book histos
 24      book(_h_alpha,1,1,1);
 25      book(_h_beta ,2,1,1);
 26      book(_h_gamma,3,1,1);
 27    }
 28
 29    bool isK0(int id) {
 30      return id==310 || id==130 || abs(id)==311;
 31    }
 32
 33    /// Perform the per-event analysis
 34    void analyze(const Event& event) {
 35      static const int DsID = 10433;
 36      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 37      for (const Particle& p : ufs.particles(Cuts::abspid==DsID)) {
 38      	// decay angle
 39      	int sign = p.pid()/DsID;
 40      	Particle Dstar,Kaon;
 41      	if(p.children().size()!=2) continue;
 42      	if(p.children()[0].pid()==sign*413 &&
 43      	   isK0(p.children()[1].pid())) {
 44      	  Dstar = p.children()[0];
 45      	  Kaon = p.children()[1];
 46      	}
 47      	else if(p.children()[1].pid()==sign*413 &&
 48      		isK0(p.children()[0].pid())) {
 49      	  Kaon = p.children()[0];
 50      	  Dstar = p.children()[1];
 51      	}
 52      	else {
 53      	  continue;
 54      	}
 55      	// first boost to the D_s1 rest frame
 56      	LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
 57      	FourMomentum pDstar = boost1.transform(Dstar.momentum());
 58      	double cTheta = pDstar.p3().unit().dot(p.momentum().p3().unit());
 59      	_h_alpha->fill(cTheta);
 60      	if(Dstar.children().size()!=2) continue;
 61      	Particle D0,Pion;
 62      	if(Dstar.children()[0].pid()== sign*211 &&
 63      	   Dstar.children()[1].pid()== sign*421) {
 64      	  Pion = Dstar.children()[0];
 65      	  D0   = Dstar.children()[1];
 66      	}
 67      	else if(Dstar.children()[1].pid()== sign*211 &&
 68      		Dstar.children()[0].pid()== sign*421) {
 69      	  D0   = Dstar.children()[0];
 70      	  Pion = Dstar.children()[1];
 71      	}
 72      	else
 73      	  continue;
 74      	// boost to D_s frame
 75      	FourMomentum pD  = boost1.transform(D0  .momentum());
 76      	FourMomentum pK  = boost1.transform(Kaon.momentum());
 77      	FourMomentum pPi = boost1.transform(Pion.momentum());
 78      	// to D* rest frame
 79      	LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(pDstar.betaVec());
 80      	Vector3 axis = pDstar.p3().unit();
 81      	FourMomentum pp = boost2.transform(pD);
 82      	// calculate angle
 83      	double cThetap = pp.p3().unit().dot(axis);
 84      	_h_gamma->fill(cThetap);
 85	// finally beta
 86	Vector3 n1 = pD  .p3().cross(pPi.p3()).unit();
 87	Vector3 n2 = axis.cross(pK.p3()).unit();
 88	double beta = acos(n1.dot(n2));
 89	_h_beta->fill(beta);
 90      }
 91    }
 92
 93
 94    /// Normalise histograms etc., after the run
 95    void finalize() {
 96      normalize(_h_alpha);
 97      normalize(_h_beta );
 98      normalize(_h_gamma);
 99    }
100
101    /// @}
102
103
104    /// @name Histograms
105    /// @{
106    Histo1DPtr _h_alpha, _h_beta, _h_gamma;
107    /// @}
108
109
110  };
111
112
113  RIVET_DECLARE_PLUGIN(BELLE_2008_I762013);
114
115}