Rivet analyses referenceMC_DECAY_OMEGAPHIA1_3PIONMC analysis of $\omega$, $\phi$ and $a_1\to3\pi$ decaysExperiment: () Status: VALIDATED Authors:
Beams: * * Beam energies: ANY Run details:
Analysis of the mass distributions and Dalitz plots in $a_1\to3\pi$ and $\omega,\phi\to\pi^+\pi^-\pi^0$ decays. Based on old Herwig++ internal analysis. Source code: MC_DECAY_OMEGAPHIA1_3PION.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/UnstableParticles.hh"
4
5namespace Rivet {
6
7
8 class MC_DECAY_OMEGAPHIA1_3PION : public Analysis {
9 public:
10
11 /// Constructor
12 RIVET_DEFAULT_ANALYSIS_CTOR(MC_DECAY_OMEGAPHIA1_3PION);
13
14
15 /// @name Analysis methods
16 /// @{
17
18 /// Book histograms and initialise projections before the run
19 void init() {
20
21 // Initialise and register projections
22 declare(UnstableParticles(), "UFS");
23 // book histograms a_1
24 // Histograms for a_10 -> pi0pi0pi0
25 book(_hist0, "hist0",200,0.2,1.5);
26 // // dalitz plot
27 book(_dalitz0, "dalitz0",50,0.2,1.5,50,0.2,1.5);
28 // Histograms for a_1+ -> pi0pi0pi+
29 // Mass of the pi0pi0 pair
30 book(_hist1A, "hist1A",200,0.2,1.5);
31 // Mass of the pi0pi+ pair
32 book(_hist1B, "hist1B",200,0.2,1.5);
33 // dalitz plot
34 book(_dalitz1, "dalitz1",50,0.2,1.5,50,0.2,1.5);
35 // Histograms for a_10 -> pi+pi-pi0
36 // Mass of the pi+pi- pair
37 book(_hist2A, "hist2A",200,0.2,1.5);
38 // Mass of the pi+pi0 pair
39 book(_hist2B, "hist2B",200,0.2,1.5);
40 // Mass of the pi-pi0 pair
41 book(_hist2C, "hist2C",200,0.2,1.5);
42 // dalitz plot
43 book(_dalitz2, "dalitz2",50,0.2,1.5,50,0.2,1.5);
44 // Histograms for a_1+ -> pi+pi+pi-
45 // Mass of the pi+pi+ pair
46 book(_hist3A, "hist3A",200,0.2,1.5);
47 // Mass of the pi+pi- pair
48 book(_hist3B, "hist3B",200,0.2,1.5);
49 // dalitz plot
50 book(_dalitz3, "dalitz3",50,0.2,1.5,50,0.2,1.5);
51
52 // Book histograms omega/phi
53 for(unsigned int ix=0;ix<2;++ix) {
54 double mmax = ix==0 ? 0.8 : 1.0;
55 std::ostringstream title1; title1 << "xhist_" << ix+1;
56 _h_xhist .push_back(Histo1DPtr());
57 book(_h_xhist.back(), title1.str(),200,-300.,300. );
58 std::ostringstream title2; title2 << "yhist_" << ix+1;
59 _h_yhist .push_back(Histo1DPtr());
60 book(_h_yhist.back(), title2.str(),200,0. ,400. );
61 std::ostringstream title3; title3 << "mplus_" << ix+1;
62 _h_mplus .push_back(Histo1DPtr());
63 book(_h_mplus.back(), title3.str(),200,200.,mmax*1000.);
64 std::ostringstream title4; title4 << "mminus_" << ix+1;
65 _h_mminus .push_back(Histo1DPtr());
66 book(_h_mminus.back(), title4.str(),200,200.,mmax*1000.);
67 std::ostringstream title5; title5 << "m0_" << ix+1;
68 _h_m0 .push_back(Histo1DPtr());
69 book(_h_m0.back(), title5.str(),200,200.,mmax*1000.);
70 std::ostringstream title6; title6 << "dalitz_" << ix+1;
71 _h_dalitz.push_back(Histo2DPtr());
72 book(_h_dalitz.back(), title6.str(),50,0.2,mmax,50,0.2,mmax);
73 }
74 }
75
76
77 void findDecayProducts(const Particle & mother, unsigned int & nstable,
78 Particles & pip , Particles & pim , Particles & pi0) {
79 for(const Particle & p : mother.children()) {
80 int id = p.pid();
81 if (id == PID::PIPLUS) {
82 pip.push_back(p);
83 ++nstable;
84 }
85 else if (id == PID::PIMINUS) {
86 pim.push_back(p);
87 ++nstable;
88 }
89 else if (id == PID::PI0) {
90 pi0.push_back(p);
91 ++nstable;
92 }
93 else if ( !p.children().empty() ) {
94 findDecayProducts(p, nstable, pip, pim, pi0);
95 }
96 else
97 ++nstable;
98 }
99 }
100
101 /// Perform the per-event analysis
102 void analyze(const Event& event) {
103 for(const Particle& meson :
104 apply<UnstableParticles>(event, "UFS").particles(Cuts::pid==PID::PHI || Cuts::pid==PID::OMEGA ||
105 Cuts::abspid==20213 || Cuts::pid==20113 )) {
106 unsigned int nstable(0);
107 Particles pip, pim, pi0;
108 findDecayProducts(meson, nstable, pip, pim, pi0);
109 if(nstable !=3) continue;
110 if(meson.pid()<0) {
111 swap(pim,pip);
112 }
113 if(meson.pid()== PID::PHI || meson.pid()==PID::OMEGA) {
114 if(pip.size()!=1 || pim.size()!=1 || pi0.size()!=1) continue;
115 unsigned int iloc = meson.pid() == PID::OMEGA ? 0 : 1;
116 LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(meson.momentum().betaVec());
117 FourMomentum pp = boost.transform(pip[0].momentum());
118 FourMomentum pm = boost.transform(pim[0].momentum());
119 FourMomentum p0 = boost.transform(pi0[0].momentum());
120 double mp = (pp+p0).mass(), mm = (pm+pp).mass();
121 _h_mplus [iloc]->fill(mp/MeV);
122 _h_mminus[iloc]->fill((pm+p0).mass()/MeV);
123 _h_m0 [iloc]->fill(mm/MeV);
124 double x = pp.t()-pm.t();
125 double y = p0.t()-p0.mass();
126 _h_xhist[iloc]->fill(x/MeV);
127 _h_yhist[iloc]->fill(y/MeV);
128 _h_dalitz[iloc]->fill(mp,mm);
129 }
130 else {
131 // a_1+ -> pi+pi+pi-
132 if(pip.size()==2&&pim.size()==1) {
133 _hist3A->fill((pip[0].momentum()+pip[1].momentum()).mass());
134 _hist3B->fill((pip[0].momentum()+pim[0].momentum()).mass());
135 _hist3B->fill((pip[1].momentum()+pim[0].momentum()).mass());
136 _dalitz3->fill((pip[0].momentum()+pim[0].momentum()).mass(),(pip[1].momentum()+pim[0].momentum()).mass());
137 _dalitz3->fill((pip[1].momentum()+pim[0].momentum()).mass(),(pip[0].momentum()+pim[0].momentum()).mass());
138 }
139 // a_1+ -> pi0pi0pi+
140 else if(pip.size()==1&&pi0.size()==2) {
141 _hist1A->fill((pi0[0].momentum()+pi0[1].momentum()).mass());
142 _hist1B->fill((pip[0].momentum()+pi0[0].momentum()).mass());
143 _hist1B->fill((pip[0].momentum()+pi0[1].momentum()).mass());
144 _dalitz1->fill((pip[0].momentum()+pi0[0].momentum()).mass(),(pip[0].momentum()+pi0[1].momentum()).mass());
145 _dalitz1->fill((pip[0].momentum()+pi0[1].momentum()).mass(),(pip[0].momentum()+pi0[0].momentum()).mass());
146 }
147 // a_10 -> pi0pi0pi0
148 else if(pi0.size()==3) {
149 _hist0->fill((pi0[0].momentum()+pi0[1].momentum()).mass());
150 _hist0->fill((pi0[0].momentum()+pi0[2].momentum()).mass());
151 _hist0->fill((pi0[1].momentum()+pi0[2].momentum()).mass());
152 _dalitz0->fill((pi0[0].momentum()+pi0[1].momentum()).mass(),(pi0[0].momentum()+pi0[2].momentum()).mass());
153 _dalitz0->fill((pi0[0].momentum()+pi0[1].momentum()).mass(),(pi0[1].momentum()+pi0[2].momentum()).mass());
154 _dalitz0->fill((pi0[0].momentum()+pi0[2].momentum()).mass(),(pi0[1].momentum()+pi0[2].momentum()).mass());
155 _dalitz0->fill((pi0[0].momentum()+pi0[2].momentum()).mass(),(pi0[0].momentum()+pi0[1].momentum()).mass());
156 _dalitz0->fill((pi0[1].momentum()+pi0[2].momentum()).mass(),(pi0[0].momentum()+pi0[1].momentum()).mass());
157 _dalitz0->fill((pi0[1].momentum()+pi0[2].momentum()).mass(),(pi0[0].momentum()+pi0[2].momentum()).mass());
158 }
159 // a_10 -> pi+pi-pi0
160 else if(pi0.size()==1&&pip.size()==1&&pim.size()==1) {
161 _hist2A->fill((pim[0].momentum()+pip[0].momentum()).mass());
162 _hist2B->fill((pip[0].momentum()+pi0[0].momentum()).mass());
163 _hist2C->fill((pim[0].momentum()+pi0[0].momentum()).mass());
164 _dalitz2->fill((pim[0].momentum()+pi0[0].momentum()).mass(),(pip[0].momentum()+pi0[0].momentum()).mass());
165 }
166 }
167 }
168 }
169
170
171 /// Normalise histograms etc., after the run
172 void finalize() {
173 // a_1
174 normalize(_hist0);
175 normalize(_dalitz0);
176 normalize(_hist1A);
177 normalize(_hist1B);
178 normalize(_dalitz1);
179 normalize(_hist2A);
180 normalize(_hist2B);
181 normalize(_hist2C);
182 normalize(_dalitz2);
183 normalize(_hist3A);
184 normalize(_hist3B);
185 normalize(_dalitz3);
186 // omega/phi
187 for(unsigned int ix=0;ix<2;++ix) {
188 normalize(_h_xhist [ix]);
189 normalize(_h_yhist [ix]);
190 normalize(_h_mplus [ix]);
191 normalize(_h_mminus[ix]);
192 normalize(_h_m0 [ix]);
193 normalize(_h_dalitz[ix]);
194 }
195 }
196
197 /// @}
198
199 /// @name Histograms a_1
200 /// @{
201 // Histograms for a_10 -> pi0pi0pi0
202 Histo1DPtr _hist0;
203 // dalitz plot
204 Histo2DPtr _dalitz0;
205 // Histograms for a_1+ -> pi0pi0pi+
206 // Mass of the pi0pi0 pair
207 Histo1DPtr _hist1A;
208 // Mass of the pi0pi+ pair
209 Histo1DPtr _hist1B;
210 // dalitz plot
211 Histo2DPtr _dalitz1;
212 // Histograms for a_10 -> pi+pi-pi0
213 // Mass of the pi+pi- pair
214 Histo1DPtr _hist2A;
215 // Mass of the pi+pi0 pair
216 Histo1DPtr _hist2B;
217 // Mass of the pi-pi0 pair
218 Histo1DPtr _hist2C;
219 // dalitz plot
220 Histo2DPtr _dalitz2;
221 // Histograms for a_1+ -> pi+pi+pi-
222 // Mass of the pi+pi+ pair
223 Histo1DPtr _hist3A;
224 // Mass of the pi+pi- pair
225 Histo1DPtr _hist3B;
226 // dalitz plot
227 Histo2DPtr _dalitz3;
228 /// @}
229
230 /// @name Histograms omega/phi
231 /// @{
232 // Histogram for the x-values
233 vector<Histo1DPtr> _h_xhist;
234 // Histogram for the y-values
235 vector<Histo1DPtr> _h_yhist;
236 // The mass of the \rho^+
237 vector<Histo1DPtr> _h_mplus;
238 // The mass of the \rho^-
239 vector<Histo1DPtr> _h_mminus;
240 // The mass of the \rho^0
241 vector<Histo1DPtr> _h_m0;
242 // Dalitz plot
243 vector<Histo2DPtr> _h_dalitz;
244 /// @}
245 };
246
247
248 RIVET_DECLARE_PLUGIN(MC_DECAY_OMEGAPHIA1_3PION);
249
250}
|