1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/UnstableParticles.hh"
 4
 5namespace Rivet {
 6
 7
 8  /// @brief B -> s gamma
 9  class BELLE_2022_I2167323 : public Analysis {
10  public:
11
12    /// Constructor
13    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2022_I2167323);
14
15
16    /// @name Analysis methods
17    /// @{
18
19    /// Book histograms and initialise projections before the run
20    void init() {
21      // Initialise and register projections
22      declare(UnstableParticles(Cuts::abspid==521 || Cuts::abspid==511), "UFS");
23      // Book histograms
24      book(_h_spectrum, 1, 1, 1);
25      book(_nBottom, "TMP/BottomCounter");
26    }
27
28    void findDecayProducts(const Particle& mother,
29                           unsigned int& nK0, unsigned int& nKp,
30			   unsigned int& nKm) {
31      for (const Particle & p : mother.children()) {
32        int id = p.pid();
33        if ( id == PID::KPLUS )      ++nKp;
34        else if (id == PID::KMINUS ) ++nKm;
35        else if (id == PID::K0S)     ++nK0;
36        else if (id == PID::PI0 || id == PID::PIPLUS || id == PID::PIMINUS) {
37	  continue;
38        }
39        else if ( !p.children().empty() ) {
40          findDecayProducts(p, nK0, nKp, nKm);
41        }
42      }
43    }
44
45    /// Perform the per-event analysis
46    void analyze(const Event& event) {
47      // Loop over bottoms
48      for (const Particle& bottom : apply<UnstableParticles>(event, "UFS").particles()) {
49	// remove mixing entries etc
50	if(bottom.children()[0].abspid()==bottom.abspid()) continue;
51        _nBottom->fill();
52	FourMomentum pgamma(0.,0.,0.,0.);
53	unsigned int ngamma = 0;
54        for (const Particle & child : bottom.children()) {
55	  if (child.pid() == PID::PHOTON) {
56            ngamma += 1;
57            pgamma += child.momentum();
58          }
59	}
60	if (ngamma != 1) continue;
61        unsigned int nK0(0),nKp(0),nKm(0);
62        FourMomentum p_tot(0,0,0,0);
63        findDecayProducts(bottom, nK0, nKp, nKm);
64        unsigned int nk = nKp-nKm+nK0;
65        if (nk % 2 == 1) {
66	  const LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(bottom.momentum().betaVec());
67	  double eGamma = boost.transform(pgamma).E();
68	  _h_spectrum->fill(eGamma);
69	}
70      }
71    }
72
73
74    /// Normalise histograms etc., after the run
75    void finalize() {
76      scale(_h_spectrum, 1e4/_nBottom->sumW());
77    }
78
79    /// @}
80
81
82    /// @name Histograms
83    /// @{
84    Histo1DPtr _h_spectrum;
85    CounterPtr _nBottom;
86    /// @}
87
88
89  };
90
91
92  RIVET_DECLARE_PLUGIN(BELLE_2022_I2167323);
93
94}