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+ -> eta' rho+
10  class CLEOII_1998_I466173 : public Analysis {
11  public:
12
13    /// Constructor
14    RIVET_DEFAULT_ANALYSIS_CTOR(CLEOII_1998_I466173);
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==431);
24      declare(ufs, "UFS");
25      DecayedParticles DS(ufs);
26      DS.addStable(PID::RHOPLUS);
27      DS.addStable(PID::ETAPRIME);
28      declare(DS,"DS");
29      // histos
30      book(_h, 1, 1, 1);
31    }
32
33
34    /// Perform the per-event analysis
35    void analyze(const Event& event) {
36      DecayedParticles DS = apply<DecayedParticles>(event, "DS");
37      // loop over particles
38      for (unsigned int ix=0; ix < DS.decaying().size(); ++ix) {
39        int sign = 1;
40        if (DS.decaying()[ix].pid()>0 && DS.modeMatches(ix,2,mode))         sign = 1;
41        else if  (DS.decaying()[ix].pid()<0 && DS.modeMatches(ix,2,modeCC)) sign =-1;
42        else continue;
43        const Particle& rhoP = DS.decayProducts()[ix].at( sign*213)[0];
44        // rho decay products
45        Particle piP;
46        if (rhoP.children().size()!=2) vetoEvent;
47        if (rhoP.children()[0].pid()==211 &&
48            rhoP.children()[1].pid()==111) {
49          piP = rhoP.children()[0];
50        }
51        else if (rhoP.children()[1].pid()==211 &&
52                 rhoP.children()[0].pid()==111) {
53          piP = rhoP.children()[1];
54        }
55        else  vetoEvent;
56        // compute the helicity angle
57        LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(DS.decaying()[ix].mom().betaVec());
58        FourMomentum prho=boost1.transform(rhoP.mom());
59        LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(prho.betaVec());
60        Vector3 axis = boost2.transform(boost1.transform(piP.mom())).p3().unit();
61        _h->fill(axis.dot(prho.p3().unit()));
62      }
63    }
64
65
66    /// Normalise histograms etc., after the run
67    void finalize() {
68      normalize(_h, 1.0, false);
69    }
70
71    /// @}
72
73
74    /// @name Histograms
75    /// @{
76    Histo1DPtr _h;
77    const map<PdgId,unsigned int> mode   = { { 213,1}, { 331,1}};
78    const map<PdgId,unsigned int> modeCC = { {-213,1}, { 331,1}};
79    /// @}
80
81
82  };
83
84
85  RIVET_DECLARE_PLUGIN(CLEOII_1998_I466173);
86
87}