Rivet analyses referenceBESII_2004_I652399Mass and angular distributions in $J/\psi\to\omega \pi^+\pi^-$Experiment: BESII (BEPC) Inspire ID: 652399 Status: VALIDATED NOHEPDATA Authors:
Beam energies: (1.6, 1.6) GeV Run details:
Mass and angular distributions in $J/\psi\to\omega \pi^+\pi^-$. The data were read from the figures in the paper and are not corrected for acceptance. Source code: BESII_2004_I652399.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/Beam.hh"
4#include "Rivet/Projections/UnstableParticles.hh"
5#include "Rivet/Projections/DecayedParticles.hh"
6
7namespace Rivet {
8
9
10 /// @brief J/psi -> omega pi+ pi-
11 class BESII_2004_I652399 : public Analysis {
12 public:
13
14 /// Constructor
15 RIVET_DEFAULT_ANALYSIS_CTOR(BESII_2004_I652399);
16
17
18 /// @name Analysis methods
19 /// @{
20
21 /// Book histograms and initialise projections before the run
22 void init() {
23 // Initialise and register projections
24 UnstableParticles ufs = UnstableParticles(Cuts::abspid==443);
25 declare(ufs, "UFS");
26 DecayedParticles PSI(ufs);
27 PSI.addStable(223);
28 declare(PSI, "PSI");
29 declare(Beam(), "Beams");
30 // histogram
31 for (unsigned int ix=0; ix<2; ++ix) {
32 book(_h_mass[ix],1,1,1+ix);
33 for (unsigned int iy=0; iy<4; ++iy) {
34 book(_h_angle[iy][ix],2,1+iy,1+ix);
35 }
36 }
37 }
38
39 void findChildren(const Particle& p, Particles& pim, Particles& pip,
40 Particles& pi0, unsigned int &ncount) {
41 for (const Particle& child : p.children()) {
42 if (child.pid()==PID::PIPLUS) {
43 pip.push_back(child);
44 ncount+=1;
45 }
46 else if (child.pid()==PID::PIMINUS) {
47 pim.push_back(child);
48 ncount+=1;
49 }
50 else if (child.pid()==PID::PI0) {
51 pi0.push_back(child);
52 ncount+=1;
53 }
54 else if (child.children().empty()) {
55 ncount+=1;
56 }
57 else {
58 findChildren(child,pim,pip,pi0,ncount);
59 }
60 }
61 }
62
63 /// Perform the per-event analysis
64 void analyze(const Event& event) {
65 // get the axis, direction of incoming electron
66 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
67 Vector3 axis;
68 if (beams.first.pid()>0) {
69 axis = beams.first .mom().p3().unit();
70 }
71 else {
72 axis = beams.second.mom().p3().unit();
73 }
74 // J/psi decay products
75 DecayedParticles PSI = apply<DecayedParticles>(event, "PSI");
76 if ( PSI.decaying().size()!=1) vetoEvent;
77 if (!PSI.modeMatches(0,3,mode)) vetoEvent;
78 const Particle& Pip = PSI.decayProducts()[0].at( 211)[0];
79 const Particle& Pim = PSI.decayProducts()[0].at(-211)[0];
80 const Particle& omega = PSI.decayProducts()[0].at( 223)[0];
81 // mass histograms
82 FourMomentum pPiPi = Pip.mom()+Pim .mom();
83 const double mPiPi = pPiPi.mass();
84 _h_mass[0]->fill(mPiPi/GeV);
85 _h_mass[1]->fill((Pip.mom()+omega.mom()).mass()/GeV);
86 _h_mass[1]->fill((Pim.mom()+omega.mom()).mass()/GeV);
87 unsigned int iplot=0;
88 if (mPiPi>.3 && mPiPi<.5) iplot=0;
89 else if (mPiPi>.5 && mPiPi<.7) iplot=1;
90 else return;
91 if (abs(Pip.p3().unit().dot(axis))>0.84) return;
92 if (abs(Pim.p3().unit().dot(axis))>0.84) return;
93 LorentzTransform boost;
94 if (PSI.decaying()[0].mom().p3().mod() > 1*MeV) {
95 boost = LorentzTransform::mkFrameTransformFromBeta(PSI.decaying()[0].mom().betaVec());
96 }
97 pPiPi = boost.transform(pPiPi);
98 // omega decay
99 unsigned int ncount=0;
100 Particles pip, pim, pi0;
101 findChildren(omega,pim,pip,pi0,ncount);
102 if ( ncount!=3 || !(pim.size()==1 && pip.size()==1 && pi0.size()==1)) return;
103 if (abs(pip[0].p3().unit().dot(axis))>0.84) return;
104 if (abs(pim[0].p3().unit().dot(axis))>0.84) return;
105 FourMomentum pOmega = boost.transform(omega.mom());
106 Vector3 e1Z = pOmega.p3().unit();
107 Vector3 e1Y = e1Z.cross(axis).unit();
108 Vector3 e1X = e1Y.cross(e1Z).unit();
109 const double cOmega = e1Z.dot(axis);
110 _h_angle[1][iplot]->fill(cOmega);
111 // pion angles
112 LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(pPiPi.betaVec());
113 FourMomentum pPi = boost2.transform(boost.transform(Pip.mom()));
114 const double cPi = pPi.p3().unit().dot(e1Z);
115 _h_angle[2][iplot]->fill(cPi);
116 Vector3 trans1 = pPi.p3().unit()-cPi*e1Z;
117 // omega angles
118 LorentzTransform boost3 = LorentzTransform::mkFrameTransformFromBeta(pOmega.betaVec());
119 FourMomentum ppip = boost3.transform(boost.transform(pip[0].mom()));
120 FourMomentum ppim = boost3.transform(boost.transform(pim[0].mom()));
121 Vector3 nW = ppip.p3().cross(ppim.p3()).unit();
122 const double bW = nW.dot(e1X);
123 _h_angle[3][iplot]->fill(bW);
124 Vector3 trans2 = ppip.p3()-ppip.p3().dot(e1Z)*e1Z;
125 const double chi = abs(atan2(trans1.cross(trans2).dot(e1Z),trans1.dot(trans2)));
126 _h_angle[0][iplot]->fill(chi/M_PI*180.);
127 }
128
129
130 /// Normalise histograms etc., after the run
131 void finalize() {
132 normalize(_h_mass, 1.0, false);
133 for (unsigned int iy=0; iy<4; ++iy) {
134 normalize(_h_angle[iy], 1.0, false);
135 }
136 }
137
138 /// @}
139
140
141 /// @name Histograms
142 /// @{
143 Histo1DPtr _h_mass[2],_h_angle[4][2];
144 const map<PdgId,unsigned int> mode = { { 223,1}, { 211,1},{-211,1}};
145 /// @}
146
147
148 };
149
150
151 RIVET_DECLARE_PLUGIN(BESII_2004_I652399);
152
153}
|