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 omega -> pi+pi-pi0
10  class BESIII_2018_I1703033 : public Analysis {
11  public:
12
13    /// Constructor
14    RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2018_I1703033);
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::pid== 223);
24      declare(ufs, "UFS");
25      DecayedParticles OMEGA(ufs);
26      OMEGA.addStable(PID::PI0);
27      declare(OMEGA,"OMEGA");
28      // histograms
29      book(_h_z  ,1,1,1);
30      book(_h_phi,1,1,2);
31      book(_h_cos,1,1,3);
32      book(_h_ppi,1,1,4);
33      book(_dalitz,"dalitz",50,-1,1,50,-1.1,0.9);
34    }
35
36    /// Perform the per-event analysis
37    void analyze(const Event& event) {
38      static const map<PdgId,unsigned int> & mode = { { 211,1}, {-211,1}, { 111,1}};
39      static const std::complex<double> ii(0.,1.);
40      DecayedParticles OMEGA = apply<DecayedParticles>(event, "OMEGA");
41      // loop over particles
42      for(unsigned int ix=0;ix<OMEGA.decaying().size();++ix) {
43	// pi+ pi- pi0
44	if (!OMEGA.modeMatches(ix,3,mode)) continue;
45	// apply omega mass cut
46	if(abs(OMEGA.decaying()[ix].mass()-0.78265)>0.04) continue;
47	// decay products
48	const Particle & pip = OMEGA.decayProducts()[ix].at( 211)[0];
49	const Particle & pim = OMEGA.decayProducts()[ix].at(-211)[0];
50	const Particle & pi0 = OMEGA.decayProducts()[ix].at( 111)[0];
51	LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(OMEGA.decaying()[ix].momentum().betaVec());
52	FourMomentum pp = boost.transform(pip.momentum());
53	FourMomentum pm = boost.transform(pim.momentum());
54	FourMomentum p0 = boost.transform(pi0.momentum());
55	// variables from eqn 1 in paper (use version from Eqn 2.1 of 1010.3946)
56	double Qc = OMEGA.decaying()[ix].mass()-pi0.mass()-pim.mass()-pip.mass();
57	double x = sqrt(3.)*(pm.E()-pp.E())/Qc;
58	double y = 3./Qc*(p0.E()-pi0.mass())-1.;
59	_dalitz->fill(x,y);
60	// plot arg and phase of variable
61	std::complex<double> zz = x+ii*y;
62	_h_z->fill(norm(zz));
63	double phi = arg(zz);
64	if(phi<0.) phi+=2.*M_PI;
65	_h_phi->fill(phi);
66	double ctheta = -pp.p3().unit().dot(p0.p3().unit());
67	_h_cos->fill(ctheta);
68	_h_ppi->fill(p0.p3().mod());
69      }
70    }
71
72
73    /// Normalise histograms etc., after the run
74    void finalize() {
75      normalize(_h_z  );
76      normalize(_h_phi);
77      normalize(_h_cos);
78      normalize(_h_ppi);
79      normalize(_dalitz);
80    }
81
82    /// @}
83
84
85    /// @name Histograms
86    /// @{
87    Histo1DPtr _h_z,_h_phi,_h_cos,_h_ppi;
88    Histo2DPtr _dalitz;
89    /// @}
90
91
92  };
93
94
95  RIVET_DECLARE_PLUGIN(BESIII_2018_I1703033);
96
97}