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+ eta gamma
10  class BELLE_2018_I1663447 : public Analysis {
11  public:
12
13    /// Constructor
14    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2018_I1663447);
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::ETA);
27      declare(BP, "BP");
28      for(unsigned int ix=0;ix<2;++ix)
29	for(unsigned int iy=0;iy<2;++iy)
30	  book(_h[ix][iy],1+ix,1,1+iy);
31    }
32
33
34    /// Perform the per-event analysis
35    void analyze(const Event& event) {
36      static const map<PdgId,unsigned int> & mode   = { { 221,1},{ 321,1}, {22,1}};
37      static const map<PdgId,unsigned int> & modeCC = { { 221,1},{-321,1}, {22,1}};
38      DecayedParticles BP = apply<DecayedParticles>(event, "BP");
39      // loop over particles
40      for(unsigned int ix=0;ix<BP.decaying().size();++ix) {
41      	int sign = 1;
42      	if (BP.decaying()[ix].pid()>0 && BP.modeMatches(ix,3,mode)) {
43      	  sign=1;
44      	}
45      	else if  (BP.decaying()[ix].pid()<0 && BP.modeMatches(ix,3,modeCC)) {
46      	  sign=-1;
47      	}
48      	else
49      	  continue;
50	const Particle & Kp    = BP.decayProducts()[ix].at( sign*321)[0];
51	const Particle & eta   = BP.decayProducts()[ix].at(      221)[0];
52	const Particle & gamma = BP.decayProducts()[ix].at(       22)[0];
53	FourMomentum pKeta = Kp.momentum()+eta.momentum();
54	double mass = pKeta.mass();
55	LorentzTransform boostB = LorentzTransform::mkFrameTransformFromBeta(BP.decaying()[ix].momentum().betaVec());
56	pKeta = boostB.transform(pKeta);
57	LorentzTransform boostKeta = LorentzTransform::mkFrameTransformFromBeta(pKeta.betaVec());
58	FourMomentum pK     = boostKeta.transform(boostB.transform(Kp   .momentum()));
59	FourMomentum pGamma = boostB.transform(gamma.momentum());
60	double cTheta = pK.p3().unit().dot(pGamma.p3().unit());
61	for(unsigned int iy=0;iy<2;++iy) {
62	  _h[0][iy]->fill(cTheta);
63	  _h[1][iy]->fill(mass );
64	}
65      }
66    }
67
68
69    /// Normalise histograms etc., after the run
70    void finalize() {
71      for(unsigned int ix=0;ix<2;++ix)
72	for(unsigned int iy=0;iy<2;++iy)
73	  normalize(_h[ix][iy],1.,false);
74    }
75
76    /// @}
77
78
79    /// @name Histograms
80    /// @{
81    Histo1DPtr _h[2][2];
82    /// @}
83
84
85  };
86
87
88  RIVET_DECLARE_PLUGIN(BELLE_2018_I1663447);
89
90}