Rivet analyses referenceMC_DALITZ_DMC analysis of Dalitz plots in three-body $D$-meson decaysExperiment: () Status: VALIDATED Authors:
Beams: * * Beam energies: ANY Run details:
Monte Carlo analysis of $D^0\to \bar{K}^0\pi^+\pi^-$, $D^0\to K^-\pi^+\pi^0$, $D^+\to K^-\pi^+\pi^+$, $D^+\to\bar{K}^0\pi^+\pi^0$, $D^+\to K^+\pi^-\pi^+$ and $D_s^+\to K^+\pi^-\pi^+$ including kinematic distributions and Dalitz plots. Source code: MC_DALITZ_D.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/DecayedParticles.hh"
4#include "Rivet/Projections/UnstableParticles.hh"
5
6namespace Rivet {
7
8
9 /// @brief Monte Carlo analysis of D-meson Dalitz decays
10 class MC_DALITZ_D : public Analysis {
11 public:
12
13 /// Constructor
14 RIVET_DEFAULT_ANALYSIS_CTOR(MC_DALITZ_D);
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==411 or
24 Cuts::abspid==421 or
25 Cuts::abspid==431);
26 declare(ufs, "UFS");
27 DecayedParticles DD(ufs);
28 DD.addStable(PID::PI0);
29 DD.addStable(PID::K0S);
30 declare(DD, "DD");
31
32 // Book histograms
33 book(_h_plus1, "h_plus1" ,200,0.,3. );
34 book(_h_minus1, "h_minus1" ,200,0.,3.2 );
35 book(_h_pipi1, "h_pipi1" ,200,0.,2. );
36 book(_h_minus2, "h_minus2" ,200,0.,3.2 );
37 book(_h_neutral2, "h_neutral2",200,0.,3.2 );
38 book(_h_pipi2, "h_pipi2" ,200,0.,2. );
39 book(_h_Kpilow3, "h_Kpilow3" ,200,0.,2. );
40 book(_h_Kpihigh3, "h_Kpihigh3" ,200,0.,3.2 );
41 book(_h_Kpiall3, "h_Kpiall3" ,200,0.,3. );
42 book(_h_pipi3, "h_pipi3" ,200,0.,2. );
43 book(_h_Kpip4, "h_Kpip4" ,200,0.,3.2 );
44 book(_h_pipi4, "h_pipi4" ,200,0.,2. );
45 book(_h_Kpi04, "h_Kpi04" ,200,0.,3.2);
46 book(_h_kppim5, "h_kppim5" ,200,0.,3. );
47 book(_h_kppip5, "h_kppip5" ,200,0.,3.1 );
48 book(_h_pippim5, "h_pippim5" ,200,0.,2. );
49 book(_h_kppim6, "h_kppim6" ,200,0.,3.5);
50 book(_h_kppip6, "h_kppip6" ,200,0.,3.5);
51 book(_h_pippim6, "h_pippim6" ,200,0.,2.5);
52 book(_h_kpkm1 ,"h_kpkm1",200,0.9,3.5);
53 book(_h_kppip7,"h_kppip7",200,0.3,3.5);
54 book(_h_kmpip1,"h_kmpip1",200,0.3,3.5);
55 book(_h_pipi5, "h_pipi5" ,200,0.,3.2 );
56 book(_h_pipi6, "h_pipi6" ,200,0.,3.2 );
57 book(_h_pipi7, "h_pipi7" ,200,0.,3.2 );
58 book(_dalitz1, "dalitz1" ,50,0.3,3.2,50,0.3,3.2);
59 book(_dalitz2, "dalitz2" ,50,0.3,3. ,50,0.3,3. );
60 book(_dalitz3, "dalitz3" ,50,0.3,2. ,50,0.07,2. );
61 book(_dalitz4, "dalitz4" ,50,0.3,3.1 ,50,0.07,2. );
62 book(_dalitz5, "dalitz5" ,50,0.,3. ,50,0.,2. );
63 book(_dalitz6, "dalitz6" ,50,0.3,3.5,50,0.07,2.5);
64 book(_dalitz7, "dalitz7" ,50,0.3,3.5,50,0.07,2.5);
65 book(_dalitz8, "dalitz8" ,50,0.,3.2,50,0.,3.2);
66 }
67
68 /// Perform the per-event analysis
69 void analyze(const Event& event) {
70 static const map<PdgId,unsigned int> & mode1 = { { 211,1}, {-211,1}, { 310,1} };
71 static const map<PdgId,unsigned int> & mode2 = { { 211,1}, {-321,1}, { 111,1} };
72 static const map<PdgId,unsigned int> & mode2CC = { {-211,1}, { 321,1}, { 111,1} };
73 static const map<PdgId,unsigned int> & mode3 = { { 211,1}, {-211,1}, { 111,1} };
74 static const map<PdgId,unsigned int> & mode4 = { { 211,2}, {-321,1} };
75 static const map<PdgId,unsigned int> & mode4CC = { {-211,2}, { 321,1} };
76 static const map<PdgId,unsigned int> & mode5 = { { 211,1}, { 111,1}, { 310,1} };
77 static const map<PdgId,unsigned int> & mode5CC = { {-211,1}, { 111,1}, { 310,1} };
78 static const map<PdgId,unsigned int> & mode6 = { { 211,1}, {-211,1}, { 321,1} };
79 static const map<PdgId,unsigned int> & mode6CC = { { 211,1}, {-211,1}, {-321,1} };
80 static const map<PdgId,unsigned int> & mode7 = { { 321,1}, {-321,1}, { 211,1} };
81 static const map<PdgId,unsigned int> & mode7CC = { { 321,1}, {-321,1}, {-211,1} };
82 DecayedParticles DD = apply<DecayedParticles>(event, "DD");
83 for (unsigned int ix=0;ix<DD.decaying().size();++ix) {
84 int sign = DD.decaying()[ix].pid()/DD.decaying()[ix].abspid();
85 if (DD.decaying()[ix].abspid()==421) {
86 if ( DD.modeMatches(ix,3,mode1 )) {
87 const Particle & K0 = DD.decayProducts()[ix].at( 310)[0];
88 const Particle & pip = DD.decayProducts()[ix].at( sign*211)[0];
89 const Particle & pim = DD.decayProducts()[ix].at(-sign*211)[0];
90 double mminus = (pim.momentum()+K0.momentum() ).mass2();
91 double mplus = (pip.momentum()+K0.momentum() ).mass2();
92 double mpipi = (pip.momentum()+pim.momentum()).mass2();
93 _h_plus1 ->fill(mplus);
94 _h_minus1 ->fill(mminus);
95 _h_pipi1 ->fill(mpipi);
96 _dalitz1 ->fill(mplus,mminus);
97 }
98 else if( ( DD.decaying()[ix].pid()>0 && DD.modeMatches(ix,3,mode2 )) ||
99 ( DD.decaying()[ix].pid()<0 && DD.modeMatches(ix,3,mode2CC))) {
100 const Particle & pi0 = DD.decayProducts()[ix].at( 111)[0];
101 const Particle & pip = DD.decayProducts()[ix].at( sign*211)[0];
102 const Particle & Km = DD.decayProducts()[ix].at(-sign*321)[0];
103 double mneut = (Km.momentum()+pip.momentum()).mass2();
104 double mminus = (Km.momentum()+pi0.momentum()).mass2();
105 double mpipi = (pip.momentum()+pi0.momentum()).mass2();
106 _h_neutral2 ->fill(mneut);
107 _h_minus2 ->fill(mminus);
108 _h_pipi2 ->fill(mpipi);
109 _dalitz2->fill(mminus,mneut);
110 }
111 else if(DD.modeMatches(ix,3,mode3)) {
112 const Particle & pi0 = DD.decayProducts()[ix].at( 111)[0];
113 const Particle & pip = DD.decayProducts()[ix].at( sign*211)[0];
114 const Particle & pim = DD.decayProducts()[ix].at(-sign*211)[0];
115 double mneut = (pim.momentum()+pip.momentum()).mass2();
116 double mminus = (pim.momentum()+pi0.momentum()).mass2();
117 double mplus = (pip.momentum()+pi0.momentum()).mass2();
118 _h_pipi5 ->fill(mplus);
119 _h_pipi6 ->fill(mminus);
120 _h_pipi7 ->fill(mneut);
121 _dalitz8 ->fill(mplus,mminus);
122 }
123 }
124 else if(DD.decaying()[ix].abspid()==411) {
125 if(DD.modeMatches(ix,3,mode4 ) || DD.modeMatches(ix,3,mode4CC)) {
126 const Particles & pip = DD.decayProducts()[ix].at( sign*211);
127 const Particle & Km = DD.decayProducts()[ix].at(-sign*321)[0];
128 double mplus = (Km.momentum() +pip[0].momentum()).mass2();
129 double mminus = (Km.momentum() +pip[1].momentum()).mass2();
130 double mpipi = (pip[0].momentum()+pip[1].momentum()).mass2();
131 if(mplus<mminus) swap(mplus,mminus);
132 _h_Kpilow3 ->fill( mminus);
133 _h_Kpihigh3->fill( mplus);
134 _h_Kpiall3 ->fill( mminus);
135 _h_Kpiall3 ->fill( mplus);
136 _h_pipi3 ->fill( mpipi);
137 _dalitz3->fill(mminus,mpipi);
138 }
139 else if(DD.modeMatches(ix,3,mode5 ) || DD.modeMatches(ix,3,mode5CC)) {
140 const Particle & pi0 = DD.decayProducts()[ix].at( 111)[0];
141 const Particle & K0 = DD.decayProducts()[ix].at( 310)[0];
142 const Particle & pip = DD.decayProducts()[ix].at( sign*211)[0];
143 double mminus = (K0.momentum()+pip.momentum()).mass2();
144 double mplus = (K0.momentum()+pi0.momentum()).mass2();
145 double mpipi = (pip.momentum()+pi0.momentum()).mass2();
146 _h_Kpip4 ->fill( mminus);
147 _h_pipi4 ->fill( mpipi );
148 _h_Kpi04 ->fill( mplus );
149 _dalitz4->fill(mplus,mpipi);
150 }
151 else if(DD.modeMatches(ix,3,mode6 ) || DD.modeMatches(ix,3,mode6CC)) {
152 const Particle & pip = DD.decayProducts()[ix].at( sign*211)[0];
153 const Particle & pim = DD.decayProducts()[ix].at(-sign*211)[0];
154 const Particle & Kp = DD.decayProducts()[ix].at( sign*321)[0];
155 double mplus = (Kp .momentum()+pip.momentum()).mass2();
156 double mminus = (Kp .momentum()+pim.momentum()).mass2();
157 double mpipi = (pip.momentum()+pim.momentum()).mass2();
158 _h_kppim5 ->fill(mminus);
159 _h_kppip5 ->fill(mplus );
160 _h_pippim5->fill(mpipi );
161 _dalitz5->fill(mminus,mpipi);
162 }
163 }
164 else if(DD.decaying()[ix].abspid()==431) {
165 if(DD.modeMatches(ix,3,mode6 ) || DD.modeMatches(ix,3,mode6CC)) {
166 const Particle & pip = DD.decayProducts()[ix].at( sign*211)[0];
167 const Particle & pim = DD.decayProducts()[ix].at(-sign*211)[0];
168 const Particle & Kp = DD.decayProducts()[ix].at( sign*321)[0];
169 double mplus = (Kp .momentum()+pip.momentum()).mass2();
170 double mminus = (Kp .momentum()+pim.momentum()).mass2();
171 double mpipi = (pip.momentum()+pim.momentum()).mass2();
172 _h_kppim6 ->fill(mminus);
173 _h_kppip6 ->fill(mplus );
174 _h_pippim6->fill(mpipi );
175 _dalitz6->fill(mminus,mpipi);
176 }
177 else if(DD.modeMatches(ix,3,mode7 ) || DD.modeMatches(ix,3,mode7CC)) {
178 const Particle & Kp = DD.decayProducts()[ix].at( sign*321)[0];
179 const Particle & Km = DD.decayProducts()[ix].at(-sign*321)[0];
180 const Particle & pip = DD.decayProducts()[ix].at( sign*211)[0];
181 double mplus = (Kp.momentum()+pip.momentum()).mass2();
182 double mminus = (Km.momentum()+pip.momentum()).mass2();
183 double mKK = (Kp.momentum()+Km .momentum()).mass2();
184 _h_kpkm1->fill(mKK);
185 _h_kppip7->fill(mplus);
186 _h_kmpip1->fill(mminus);
187 _dalitz7->fill(mKK,mminus);
188 }
189 }
190 }
191 }
192
193 /// Normalise histograms etc., after the run
194 void finalize() {
195 normalize(_h_plus1);
196 normalize(_h_minus1);
197 normalize(_h_pipi1);
198 normalize(_dalitz1);
199 normalize(_h_minus2);
200 normalize(_h_pipi2);
201 normalize(_h_neutral2);
202 normalize(_dalitz2);
203 normalize(_h_Kpilow3);
204 normalize(_h_Kpihigh3);
205 normalize(_h_Kpiall3);
206 normalize(_h_pipi3);
207 normalize(_dalitz3);
208 normalize(_h_Kpip4);
209 normalize(_h_pipi4);
210 normalize(_h_Kpi04);
211 normalize(_dalitz4);
212 normalize(_h_kppim5);
213 normalize(_h_kppip5);
214 normalize(_h_pippim5);
215 normalize(_dalitz5);
216 normalize(_h_kppim6);
217 normalize(_h_kppip6);
218 normalize(_h_pippim6);
219 normalize(_dalitz6);
220 normalize(_h_kpkm1);
221 normalize(_h_kppip7);
222 normalize(_h_kmpip1);
223 normalize(_dalitz7);
224 normalize(_h_pipi5);
225 normalize(_h_pipi6);
226 normalize(_h_pipi7);
227 normalize(_dalitz8);
228 }
229 /// @}
230
231 /// @name Histograms
232 /// @{
233 // Histograms for D^0\to \bar{K}^0\pi^+\pi^-
234 //m^2_+
235 Histo1DPtr _h_plus1;
236 //m^2_+
237 Histo1DPtr _h_minus1;
238 //m^2_{\pi\pi}
239 Histo1DPtr _h_pipi1;
240 // Dalitz plot
241 Histo2DPtr _dalitz1;
242
243 // Histograms for D^0\to K^-\pi^+\pi^0
244 // Histogram for the K^-\pi^+ mass
245 Histo1DPtr _h_minus2;
246 // Histogram for the \pi^+\pi^0 mass
247 Histo1DPtr _h_pipi2;
248 // Histogram for the K^-\pi^0 mass
249 Histo1DPtr _h_neutral2;
250 // Dalitz plot
251 Histo2DPtr _dalitz2;
252
253 // Histograms for D^+\to K^-\pi^+\pi^+
254 // Histogram for K^-\pi^+ low
255 Histo1DPtr _h_Kpilow3;
256 // Histogram for K^-\pi^+ high
257 Histo1DPtr _h_Kpihigh3;
258 // Histogram for K^-\pi^+ all
259 Histo1DPtr _h_Kpiall3;
260 // Histogram for \pi^+\pi^-
261 Histo1DPtr _h_pipi3;
262 // Dalitz plot
263 Histo2DPtr _dalitz3;
264
265 // Histograms for D^+\to\bar{K}^0\pi^+\pi^0
266 // Histogram for the \bar{K}^0\pi^+ mass
267 Histo1DPtr _h_Kpip4;
268 // Histogram for the \pi^+\pi^0 mass
269 Histo1DPtr _h_pipi4;
270 // Histogram for the \bar{K}^0\pi^0 mass
271 Histo1DPtr _h_Kpi04;
272 // Dalitz plot
273 Histo2DPtr _dalitz4;
274
275 // Histograms for D^+\to K^+\pi^-\pi^+
276 // Histogram for K^+\pi^-
277 Histo1DPtr _h_kppim5;
278 // Histogram for K^+\pi^+
279 Histo1DPtr _h_kppip5;
280 // Histogram for \pi^+\pi^-
281 Histo1DPtr _h_pippim5;
282 // Dalitz plot
283 Histo2DPtr _dalitz5;
284
285 // Histograms for D_s^+\to K^+\pi^-\pi^+
286 // Histogram for K^+\pi^-
287 Histo1DPtr _h_kppim6;
288 // Histogram for K^+\pi^+
289 Histo1DPtr _h_kppip6;
290 // Histogram for \pi^+\pi^-
291 Histo1DPtr _h_pippim6;
292 // Dalitz plot
293 Histo2DPtr _dalitz6;
294
295 // Histograms for D_s^+\to K^+K^-\pi^+
296 // Histogram for K^+K^-
297 Histo1DPtr _h_kpkm1;
298 // Histogram for K^+\pi^+
299 Histo1DPtr _h_kppip7;
300 // Histogram for K^-\pi^+
301 Histo1DPtr _h_kmpip1;
302 // Dalitz plot
303 Histo2DPtr _dalitz7;
304
305 // Histograms for D0 -> pi+pi-pi0
306 // Histogram for pi+pi0
307 Histo1DPtr _h_pipi5;
308 // Histogram for pi-pi0
309 Histo1DPtr _h_pipi6;
310 // Histogram for pi+pi-
311 Histo1DPtr _h_pipi7;
312 // Dalitz plot
313 Histo2DPtr _dalitz8;
314 /// @}
315
316 };
317
318
319 RIVET_DECLARE_PLUGIN(MC_DALITZ_D);
320
321}
|