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+ -> K+ K- pi+
10  class BELLE_2023_I2099998 : public Analysis {
11  public:
12
13    /// Constructor
14    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2023_I2099998);
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==521);
24      declare(ufs, "UFS");
25      DecayedParticles BP(ufs);
26      BP.addStable(PID::K0S);
27      BP.addStable(  411);
28      BP.addStable( -411);
29      BP.addStable(  421);
30      BP.addStable( -421);
31      BP.addStable(10441);
32      declare(BP, "BP");
33      // histograms
34      for (unsigned int ix=0; ix<3; ++ix) {
35      	book(_h[ix], 1+ix, 1, 1);
36      }
37      book(_c,"TMP/nB");
38    }
39
40
41    /// Perform the per-event analysis
42    void analyze(const Event& event) {
43      DecayedParticles BP = apply<DecayedParticles>(event, "BP");
44      // loop over particles
45      for (unsigned int ix=0; ix<BP.decaying().size(); ++ix) {
46        int sign = BP.decaying()[ix].pid()>0 ? 1 : -1;
47        _c->fill();
48        if ((sign>0 and BP.modeMatches(ix,3,mode  )) ||
49            (sign<0 and BP.modeMatches(ix,3,modeCC))) {
50          // 	  // momenta
51          FourMomentum Kp  = BP.decayProducts()[ix].at( 321*sign)[0].mom();
52          FourMomentum Km  = BP.decayProducts()[ix].at(-321*sign)[0].mom();
53          FourMomentum pip = BP.decayProducts()[ix].at( 211*sign)[0].mom();
54          // masses
55          double mKpi = (Km+pip).mass();
56          double mKK = (Kp+Km).mass();
57          _h[0]->fill( mKK);
58          _h[1]->fill(mKpi);
59          if (mKK>1.1) continue;
60          LorentzTransform boost =
61            LorentzTransform::mkFrameTransformFromBeta(BP.decaying()[ix]. mom().betaVec());
62          FourMomentum pKK = boost.transform(Kp+Km);
63          LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(pKK.betaVec());
64          double cTheta = boost2.transform(boost.transform(Kp)).p3().unit().dot(pKK.p3().unit());
65          _h[2]->fill(cTheta);
66        }
67      }
68    }
69
70
71    /// Normalise histograms etc., after the run
72    void finalize() {
73      scale(_h, 1e7/ *_c);
74      normalize(_h[2], 1.0, false);
75    }
76
77    /// @}
78
79
80    /// @name Histograms
81    /// @{
82    Histo1DPtr _h[3];
83    CounterPtr _c;
84    const map<PdgId,unsigned int> mode   = { { 321,1}, {-321,1}, { 211,1}};
85    const map<PdgId,unsigned int> modeCC = { { 321,1}, {-321,1}, {-211,1}};
86    /// @}
87
88
89  };
90
91
92  RIVET_DECLARE_PLUGIN(BELLE_2023_I2099998);
93
94}