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+ -> pi+ pi0 eta'
10  class BESIII_2022_I2030993 : public Analysis {
11  public:
12
13    /// Constructor
14    RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2022_I2030993);
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::PI0);
27      DS.addStable(PID::K0S);
28      DS.addStable(PID::ETAPRIME);
29      declare(DS,"DS");
30      // histograms
31      book(_h_etapip,1,1,1);
32      book(_h_etapi0,1,1,2);
33      book(_h_pipi , 1,1,3);
34      book(_dalitz, "dalitz",50,0.,1.1,50,1,3.5);
35    }
36
37    /// Perform the per-event analysis
38    void analyze(const Event& event) {
39      static const map<PdgId,unsigned int> & mode   = { { 211,1}, { 111,1}, { 331,1}};
40      static const map<PdgId,unsigned int> & modeCC = { {-211,1}, { 111,1}, { 331,1}};
41      DecayedParticles DS = apply<DecayedParticles>(event, "DS");
42      // loop over particles
43      for(unsigned int ix=0;ix<DS.decaying().size();++ix) {
44	int sign = 1;
45	if (DS.decaying()[ix].pid()>0 && DS.modeMatches(ix,3,mode)) {
46	  sign=1;
47	}
48	else if  (DS.decaying()[ix].pid()<0 && DS.modeMatches(ix,3,modeCC)) {
49	  sign=-1;
50	}
51	else
52	  continue;
53	const Particle  & pip = DS.decayProducts()[ix].at( sign*211)[0];
54	const Particle  & pi0 = DS.decayProducts()[ix].at(      111)[0];
55	const Particle  & eta = DS.decayProducts()[ix].at(      331)[0];
56	double mplus = (eta.momentum()+pip.momentum()).mass2();
57	double mzero = (eta.momentum()+pi0.momentum()).mass2();
58	double mpipi = (pip.momentum()+pi0.momentum()).mass2();
59	_dalitz  ->fill(mpipi,mplus);
60	_h_pipi  ->fill(sqrt(mpipi));
61	_h_etapip->fill(sqrt(mplus));
62	_h_etapi0->fill(sqrt(mzero));
63      }
64    }
65
66
67    /// Normalise histograms etc., after the run
68    void finalize() {
69      normalize(_h_pipi  ,1.,false);
70      normalize(_h_etapip,1.,false);
71      normalize(_h_etapi0,1.,false);
72      normalize(_dalitz  ,1.,false);
73    }
74
75    /// @}
76
77
78    /// @name Histograms
79    /// @{
80    Histo1DPtr _h_pipi, _h_etapip,_h_etapi0;
81    Histo2DPtr _dalitz;
82    /// @}
83
84
85  };
86
87
88  RIVET_DECLARE_PLUGIN(BESIII_2022_I2030993);
89
90}