Rivet analyses referenceMC_DECAY_TAUAnalysis of Kinematic distributions for $\tau$ lepton decaysExperiment: () Status: VALIDATED Authors:
Beams: * * Beam energies: ANY Run details:
Simple analysis of kinematic distributions in tau lepton decays. This includes the mass distribution of the hadronic decay products of the $\tau$ in the 2, 3, 4 and 5 hadron decays. The two hadron modes included are $\tau^-\to\nu_\tau\{\pi^-\pi^0,K^-\pi^0,K^0\pi^-,K^-\eta,K^-K^0\}$. The three hadron modes included are $\tau^-\to\nu_\tau \pi^+\pi^-\pi^-$, $\tau^-\to\nu_\tau \pi^0\pi^0\pi^-$, $\tau^-\to\nu_\tau K^-K^+\pi^-$, $\tau^-\to\nu_\tau K^0\bar{K}^0\pi^-$, $\tau^-\to\nu_\tau K^-K^0\pi^0$, $\tau^-\to\nu_\tau \pi^0-\pi^0K^-$,$\tau^-\to\nu_\tau K^-\pi^-\pi^+$, $\tau^-\to\nu_\tau \pi^-K^0\pi^0$, $\tau^-\to\nu_\tau \pi^-\pi^0\eta$, $\tau^-\to\nu_\tau \pi^-\pi^0\gamma$. The mass distributions in the four and five pion decays are included. The leptonic modes are also included. Charge conjugate modes are combined. This is based on a number of old Herwig internal analyses. Source code: MC_DECAY_TAU.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalState.hh"
4#include "Rivet/Projections/FastJets.hh"
5#include "Rivet/Projections/UnstableParticles.hh"
6
7namespace Rivet {
8
9
10 /// @brief Tau-lepton decay observables
11 class MC_DECAY_TAU : public Analysis {
12 public:
13
14 /// Constructor
15 RIVET_DEFAULT_ANALYSIS_CTOR(MC_DECAY_TAU);
16
17
18 /// @name Analysis methods
19 /// @{
20
21 /// Book histograms and initialise projections before the run
22 void init() {
23
24 // Initialise and register projections
25 declare(UnstableParticles(),"UFS");
26
27 // Book histograms
28
29 // Leptonic
30 book(_h_2B_m2enu, "h_2B_m2enu", 200,0.,3.15);
31 book(_h_2B_menu, "h_2B_menu" , 200,0.,1.8 );
32
33 // 1 hadron
34 book(_h_1B_xpi, "h_1B_xpi", 25,0.0,1.0);
35
36 // 2 hadrons
37 book(_h_2B_m2pipi, "h_2B_m2pipi", 200,0.,3.15);
38 book(_h_2B_mpipi, "h_2B_mpipi" , 200,0.,1.8 );
39 book(_h_2B_m2munu, "h_2B_m2munu", 200,0.,3.15);
40 book(_h_2B_mmunu, "h_2B_mmunu" , 200,0.,1.8 );
41 book(_h_2B_m2KpiA, "h_2B_m2KpiA", 200,0.,3.15);
42 book(_h_2B_mKpiA, "h_2B_mKpiA" , 200,0.,1.8 );
43 book(_h_2B_m2KpiB, "h_2B_m2KpiB", 200,0.,3.15);
44 book(_h_2B_mKpiB, "h_2B_mKpiB" , 200,0.,1.8 );
45 book(_h_2B_m2Keta, "h_2B_m2Keta", 200,0.,3.15);
46 book(_h_2B_mKeta, "h_2B_mKeta" , 200,0.,1.8 );
47 book(_h_2B_m2KK, "h_2B_m2KK" , 200,0.,3.15);
48 book(_h_2B_mKK, "h_2B_mKK" , 200,0.,1.8 );
49
50 // 3 hadrons
51 Histo1DPtr dummy;
52 for (size_t ix = 0; ix < 4; ++ix) {
53 if (ix < 3) {
54 book(dummy, strcat("h_3B_pippimpim_", ix+1), 200,0.,1.8);
55 _h_3B_pippimpim.push_back(dummy);
56 book(dummy, strcat("h_3B_pi0pi0pim_", ix+1), 200,0.,1.8);
57 _h_3B_pi0pi0pim.push_back(dummy);
58 book(dummy, strcat("h_3B_pi0pi0km_", ix+1), 200,0.,1.8);
59 _h_3B_pi0pi0km.push_back(dummy);
60 book(dummy, strcat("h_3B_kspimks_", ix+1), 200,0.,1.8);
61 _h_3B_kspimks.push_back(dummy);
62 book(dummy, strcat("h_3B_klpimkl_", ix+1), 200,0.,1.8);
63 _h_3B_klpimkl.push_back(dummy);
64 }
65 book(dummy, strcat("h_3B_kmpimkp_", ix+1), 200,0.,1.8);
66 _h_3B_kmpimkp.push_back(dummy);
67 book(dummy, strcat("h_3B_kmpi0k0_", ix+1), 200,0.,1.8);
68 _h_3B_kmpi0k0.push_back(dummy);
69 book(dummy, strcat("h_3B_kmpimpip_", ix+1), 200,0.,1.8);
70 _h_3B_kmpimpip.push_back(dummy);
71 book(dummy, strcat("h_3B_pimk0pi0_", ix+1), 200,0.,1.8);
72 _h_3B_pimk0pi0.push_back(dummy);
73 book(dummy, strcat("h_3B_pimpi0eta_", ix+1), 200,0.,1.8);
74 _h_3B_pimpi0eta.push_back(dummy);
75 book(dummy, strcat("h_3B_pimpi0gamma_", ix+1), 200,0.,1.8);
76 _h_3B_pimpi0gamma.push_back(dummy);
77 book(dummy, strcat("h_3B_kspimkl_", ix+1), 200,0.,1.8);
78 _h_3B_kspimkl.push_back(dummy);
79 }
80 // 4 pion decays
81 for (size_t ix=0;ix<5;++ix) {
82 book(dummy, strcat("h_4B_pipi_", ix+1), 200,0.,1.8);
83 _h_4B_pipi.push_back(dummy);
84 book(dummy, strcat("h_4B_pipipi_", ix+1), 200,0.,1.8);
85 _h_4B_pipipi.push_back(dummy);
86 }
87 book(dummy, "h_4B_pipi_6", 200,0.,1.8);
88 _h_4B_pipi.push_back(dummy);
89 for (size_t ix=0;ix<2;++ix) {
90 book(dummy, strcat("h_4B_pipipipi_", ix+1), 200,0.,1.8);
91 _h_4B_pipipipi.push_back(dummy);
92 }
93 // 5 pion decays
94 // 2 pi0 2pi- pi+
95 book(_h_5B_q1, "h_5B_q1",200,0.,1.8);
96 for (size_t ix=0;ix<5;++ix) {
97 book(dummy, strcat("h_5B_pipi1_", ix+1), 200,0.,1.8);
98 _h_5B_pipi1.push_back(dummy);
99 }
100 for (size_t ix=0;ix<5;++ix) {
101 book(dummy, strcat("h_5B_pipipi1_", ix+1), 200,0.,1.8);
102 _h_5B_pipipi1.push_back(dummy);
103 }
104 for (size_t ix=0;ix<3;++ix) {
105 book(dummy, strcat("h_5B_pipipipi1_", ix+1), 200,0.,1.8);
106 _h_5B_pipipipi1.push_back(dummy);
107 }
108 // 4 pi0 pi-
109 book(_h_5B_q2, "h_5B_q2",200,0.,1.8);
110 for (size_t ix=0;ix<2;++ix) {
111 book(dummy, strcat("h_5B_pipi2_", ix+1), 200,0.,1.8);
112 _h_5B_pipi2.push_back(dummy);
113 }
114 for (size_t ix=0;ix<2;++ix) {
115 book(dummy, strcat("h_5B_pipipi2_", ix+1), 200,0.,1.8);
116 _h_5B_pipipi2.push_back(dummy);
117 }
118 for (size_t ix=0;ix<2;++ix) {
119 book(dummy, strcat("h_5B_pipipipi2_", ix+1), 200,0.,1.8);
120 _h_5B_pipipipi2.push_back(dummy);
121 }
122 // 3 pi- 2 pi+
123 book(_h_5B_q3, "h_5B_q3",200,0.,1.8);
124 for (size_t ix=0;ix<3;++ix) {
125 book(dummy, strcat("h_5B_pipi3_", ix+1), 200,0.,1.8);
126 _h_5B_pipi3.push_back(dummy);
127 }
128 for (size_t ix=0;ix<3;++ix) {
129 book(dummy, strcat("h_5B_pipipi3_", ix+1), 200,0.,1.8);
130 _h_5B_pipipi3.push_back(dummy);
131 }
132 for (size_t ix=0;ix<2;++ix) {
133 book(dummy, strcat("h_5B_pipipipi3_", ix+1), 200,0.,1.8);
134 _h_5B_pipipipi3.push_back(dummy);
135 }
136 }
137
138
139 void findDecayProducts(const Particle & mother, size_t & nstable,
140 Particles & ep , Particles & em , Particles & nu_e , Particles & nu_ebar,
141 Particles & mup , Particles & mum , Particles & nu_mu, Particles & nu_mubar,
142 Particles & pip , Particles & pim , Particles & pi0 ,
143 Particles & Kp , Particles & Km , Particles & K0S , Particles & K0L,
144 Particles & eta , Particles & gamma) {
145 for (const Particle & p : mother.children()) {
146 int id = p.pid();
147 if ( id == PID::KPLUS ) {
148 Kp.push_back(p);
149 ++nstable;
150 }
151 else if (id == PID::KMINUS ) {
152 Km.push_back(p);
153 ++nstable;
154 }
155 else if (id == PID::PIPLUS) {
156 pip.push_back(p);
157 ++nstable;
158 }
159 else if (id == PID::PIMINUS) {
160 pim.push_back(p);
161 ++nstable;
162 }
163 else if (id == PID::EPLUS) {
164 ep.push_back(p);
165 ++nstable;
166 }
167 else if (id == PID::EMINUS) {
168 em.push_back(p);
169 ++nstable;
170 }
171 else if (id == PID::NU_E) {
172 nu_e.push_back(p);
173 ++nstable;
174 }
175 else if (id == PID::NU_EBAR) {
176 nu_ebar.push_back(p);
177 ++nstable;
178 }
179 else if (id == PID::NU_MU) {
180 nu_mu.push_back(p);
181 ++nstable;
182 }
183 else if (id == PID::NU_MUBAR) {
184 nu_mubar.push_back(p);
185 ++nstable;
186 }
187 else if (id == PID::ANTIMUON) {
188 mup.push_back(p);
189 ++nstable;
190 }
191 else if (id == PID::MUON) {
192 mum.push_back(p);
193 ++nstable;
194 }
195 else if (id == PID::PI0) {
196 pi0.push_back(p);
197 ++nstable;
198 }
199 else if (id == PID::K0S) {
200 K0S.push_back(p);
201 ++nstable;
202 }
203 else if (id == PID::K0L) {
204 K0L.push_back(p);
205 ++nstable;
206 }
207 else if (id == PID::ETA) {
208 eta.push_back(p);
209 ++nstable;
210 }
211 else if (id == PID::PHOTON) {
212 gamma.push_back(p);
213 ++nstable;
214 }
215 else if ( !p.children().empty() ) {
216 findDecayProducts(p, nstable,ep,em,nu_e,nu_ebar,mup,mum,nu_mu,nu_mubar,
217 pip, pim, pi0,Kp , Km, K0S, K0L,eta,gamma);
218 }
219 else
220 ++nstable;
221 }
222 }
223
224
225 /// Perform the per-event analysis
226 void analyze(const Event& event) {
227 for (const Particle& tau : apply<UnstableParticles>(event, "UFS").particles(Cuts::abspid==PID::TAU)) {
228 size_t nstable(0);
229 Particles ep,em,nu_e,nu_ebar,mup,mum,nu_mu,nu_mubar;
230 Particles pip, pim, pi0, Kp , Km, K0S, K0L, eta,gamma;
231 findDecayProducts(tau, nstable,ep,em,nu_e,nu_ebar,mup,mum,nu_mu,nu_mubar,
232 pip, pim, pi0,Kp , Km, K0S, K0L,eta,gamma);
233 if(tau.pid()<0) {
234 swap(pim,pip);
235 swap(Kp,Km);
236 swap(em,ep);
237 swap(mum,mup);
238 swap(nu_e ,nu_ebar );
239 swap(nu_mu,nu_mubar);
240 }
241 // cerr << "testing before loop " << nstable << " "
242 // << pip.size() << " " << pim.size() << " " << pi0.size() << " "
243 // << Kp.size() << " " << Km.size() << " " << K0S.size() << " " << K0L.size() << "\n";
244 // 2 hadrons
245 if (nstable==2) {
246 if (pim.size()==1) {
247 double xpi = pim[0].mom().p()/tau.mom().p();
248 _h_1B_xpi->fill(xpi);
249 }
250 }
251 else if (nstable==3 ) {
252 if (em.size()==1 && nu_ebar.size()==1) {
253 FourMomentum ptot = em[0].mom()+nu_ebar[0].mom();
254 double mass2 = ptot.mass2();
255 _h_2B_m2enu->fill(mass2/GeV2);
256 _h_2B_menu ->fill(sqrt(mass2)/GeV);
257 }
258 else if (mum.size()==1 && nu_mubar.size()==1) {
259 FourMomentum ptot = mum[0].mom()+nu_mubar[0].mom();
260 double mass2 = ptot.mass2();
261 _h_2B_m2munu->fill(mass2/GeV2);
262 _h_2B_mmunu ->fill(sqrt(mass2)/GeV);
263 }
264 else if (pim.size()==1 && pi0.size()==1) {
265 FourMomentum ptot = pim[0].mom()+pi0[0].mom();
266 double mass2 = ptot.mass2();
267 _h_2B_m2pipi->fill(mass2/GeV2);
268 _h_2B_mpipi ->fill(sqrt(mass2)/GeV);
269 }
270 else if (Km.size()==1 && pi0.size()==1) {
271 FourMomentum ptot = Km[0].mom()+pi0[0].mom();
272 double mass2 = ptot.mass2();
273 _h_2B_m2KpiA->fill(mass2/GeV2);
274 _h_2B_mKpiA ->fill(sqrt(mass2)/GeV);
275 }
276 else if (K0S.size()==1 && pim.size()==1) {
277 FourMomentum ptot = K0S[0].mom()+pim[0].mom();
278 double mass2 = ptot.mass2();
279 _h_2B_m2KpiB->fill(mass2/GeV2);
280 _h_2B_mKpiB ->fill(sqrt(mass2)/GeV);
281 }
282 else if (K0L.size()==1 && pim.size()==1) {
283 FourMomentum ptot = K0L[0].mom()+pim[0].mom();
284 double mass2 = ptot.mass2();
285 _h_2B_m2KpiB->fill(mass2/GeV2);
286 _h_2B_mKpiB ->fill(sqrt(mass2)/GeV);
287 }
288 else if (K0S.size()==1 && Km.size()==1) {
289 FourMomentum ptot = K0S[0].mom()+Km[0].mom();
290 double mass2 = ptot.mass2();
291 _h_2B_m2KK->fill(mass2/GeV2);
292 _h_2B_mKK ->fill(sqrt(mass2)/GeV);
293 }
294 else if (K0L.size()==1 && Km.size()==1) {
295 FourMomentum ptot = K0L[0].mom()+Km[0].mom();
296 double mass2 = ptot.mass2();
297 _h_2B_m2KK->fill(mass2/GeV2);
298 _h_2B_mKK ->fill(sqrt(mass2)/GeV);
299 }
300 else if (eta.size()==1 && Km.size()==1) {
301 FourMomentum ptot = eta[0].mom()+Km[0].mom();
302 double mass2 = ptot.mass2();
303 _h_2B_m2Keta->fill(mass2/GeV2);
304 _h_2B_mKeta ->fill(sqrt(mass2)/GeV);
305 }
306 }
307 else if (nstable==4) {
308 if (pim.size()==2 && pip.size()==1) {
309 _h_3B_pippimpim[0]->fill((pim[0].mom()+pim[1].mom()+pip[0].mom()).mass()/GeV);
310 _h_3B_pippimpim[1]->fill((pim[0].mom()+pim[1].mom()).mass()/GeV);
311 _h_3B_pippimpim[2]->fill((pim[0].mom()+pip[0].mom()).mass()/GeV);
312 _h_3B_pippimpim[2]->fill((pim[1].mom()+pip[0].mom()).mass()/GeV);
313 }
314 else if (pim.size()==1 && pi0.size()==2) {
315 _h_3B_pi0pi0pim[0]->fill((pi0[0].mom()+pi0[1].mom()+pim[0].mom()).mass()/GeV);
316 _h_3B_pi0pi0pim[1]->fill((pi0[0].mom()+pi0[1].mom()).mass()/GeV);
317 _h_3B_pi0pi0pim[2]->fill((pi0[0].mom()+pim[0].mom()).mass()/GeV);
318 _h_3B_pi0pi0pim[2]->fill((pi0[1].mom()+pim[0].mom()).mass()/GeV);
319 }
320 else if (Km.size()==1 && Kp.size()==1 && pim.size()==1) {
321 _h_3B_kmpimkp[0]->fill((Km[0].mom()+pim[0].mom()+Kp[0].mom()).mass()/GeV);
322 _h_3B_kmpimkp[1]->fill((Km[0].mom()+pim[0].mom()).mass()/GeV);
323 _h_3B_kmpimkp[2]->fill((Km[0].mom()+ Kp[0].mom()).mass()/GeV);
324 _h_3B_kmpimkp[3]->fill((Kp[0].mom()+pim[0].mom()).mass()/GeV);
325 }
326 else if ((K0S.size()==1||K0L.size()==1) && Km.size()==1 && pi0.size()==1) {
327 FourMomentum pk = K0L.size()==1 ? K0L[0].mom() : K0S[0].mom();
328 _h_3B_kmpi0k0[0]->fill((Km[0].mom()+pi0[0].mom()+pk).mass()/GeV);
329 _h_3B_kmpi0k0[1]->fill((Km[0].mom()+pi0[0].mom()).mass()/GeV);
330 _h_3B_kmpi0k0[2]->fill((Km[0].mom()+pk ).mass()/GeV);
331 _h_3B_kmpi0k0[3]->fill((pk+pi0[0].mom()).mass()/GeV);
332 }
333 else if (pi0.size()==2 && Km.size()==1) {
334 _h_3B_pi0pi0km[0]->fill((pi0[0].mom()+pi0[1].mom()+Km[0].mom()).mass()/GeV);
335 _h_3B_pi0pi0km[1]->fill((pi0[0].mom()+pi0[1].mom()).mass()/GeV);
336 _h_3B_pi0pi0km[2]->fill((pi0[0].mom()+Km[0].mom() ).mass()/GeV);
337 _h_3B_pi0pi0km[2]->fill((pi0[1].mom()+Km[0].mom() ).mass()/GeV);
338 }
339 else if (Km.size()==1 && pim.size()==1 && pip.size()==1) {
340 _h_3B_kmpimpip[0]->fill((pip[0].mom()+pim[0].mom()+Km[0].mom()).mass()/GeV);
341 _h_3B_kmpimpip[1]->fill((Km[0].mom()+pim[0].mom()).mass()/GeV);
342 _h_3B_kmpimpip[2]->fill((Km[0].mom()+pip[0].mom() ).mass()/GeV);
343 _h_3B_kmpimpip[3]->fill((pip[0].mom()+pim[0].mom() ).mass()/GeV);
344 }
345 else if (pim.size()==1 && (K0S.size()==1||K0L.size()==1) && pi0.size()==1) {
346 FourMomentum pk = K0L.size()==1 ? K0L[0].mom() : K0S[0].mom();
347 _h_3B_pimk0pi0[0]->fill((pim[0].mom()+pi0[0].mom()+pk).mass()/GeV);
348 _h_3B_pimk0pi0[1]->fill((pim[0].mom()+pk).mass()/GeV);
349 _h_3B_pimk0pi0[2]->fill((pim[0].mom()+pi0[0].mom() ).mass()/GeV);
350 _h_3B_pimk0pi0[3]->fill((pk+pi0[0].mom()).mass()/GeV);
351 }
352 else if (pim.size()==1 && pi0.size()==1 && eta.size()==1) {
353 _h_3B_pimpi0eta[0]->fill((pim[0].mom()+pi0[0].mom()+eta[0].mom()).mass()/GeV);
354 _h_3B_pimpi0eta[1]->fill((pim[0].mom()+pi0[0].mom()).mass()/GeV);
355 _h_3B_pimpi0eta[2]->fill((pim[0].mom()+eta[0].mom()).mass()/GeV);
356 _h_3B_pimpi0eta[3]->fill((pi0[0].mom()+eta[0].mom()).mass()/GeV);
357 }
358 else if (pim.size()==1 && pi0.size()==1 && gamma.size()==1) {
359 _h_3B_pimpi0gamma[0]->fill((pim[0].mom()+pi0[0].mom()+gamma[0].mom()).mass()/GeV);
360 _h_3B_pimpi0gamma[1]->fill((pim[0].mom()+pi0[0].mom()).mass()/GeV);
361 _h_3B_pimpi0gamma[2]->fill((pim[0].mom()+gamma[0].mom()).mass()/GeV);
362 _h_3B_pimpi0gamma[3]->fill((pi0[0].mom()+gamma[0].mom()).mass()/GeV);
363 }
364 else if (K0S.size()==2 && pim.size()==1) {
365 _h_3B_kspimks[0]->fill((pim[0].mom()+K0S[0].mom()+K0S[1].mom()).mass()/GeV);
366 _h_3B_kspimks[1]->fill((pim[0].mom()+K0S[0].mom()).mass()/GeV);
367 _h_3B_kspimks[1]->fill((pim[0].mom()+K0S[1].mom()).mass()/GeV);
368 _h_3B_kspimks[2]->fill((K0S [0].mom()+K0S[1].mom()).mass()/GeV);
369 }
370 else if (K0L.size()==2 && pim.size()==1) {
371 _h_3B_klpimkl[0]->fill((pim[0].mom()+K0L[0].mom()+K0L[1].mom()).mass()/GeV);
372 _h_3B_klpimkl[1]->fill((pim[0].mom()+K0L[0].mom()).mass()/GeV);
373 _h_3B_klpimkl[1]->fill((pim[0].mom()+K0L[1].mom()).mass()/GeV);
374 _h_3B_klpimkl[2]->fill((K0L [0].mom()+K0L[1].mom()).mass()/GeV);
375 }
376 else if (K0S.size()==1 && K0L.size()==1 && pim.size()==1) {
377 _h_3B_kspimkl[0]->fill((pim[0].mom()+K0S[0].mom()+K0L[0].mom()).mass()/GeV);
378 _h_3B_kspimkl[1]->fill((pim[0].mom()+K0S[0].mom()).mass()/GeV);
379 _h_3B_kspimkl[2]->fill((K0S[0].mom() +K0L[0].mom()).mass()/GeV);
380 _h_3B_kspimkl[3]->fill((pim[0].mom()+K0L[0].mom()).mass()/GeV);
381 }
382 }
383 else if (nstable==5) {
384 if (pi0.size()==3 && pim.size()==1) {
385 _h_4B_pipi[0] ->fill( (pi0[0].mom()+pim[0].mom()).mass()/GeV);
386 _h_4B_pipi[0] ->fill( (pi0[1].mom()+pim[0].mom()).mass()/GeV);
387 _h_4B_pipi[0] ->fill( (pi0[2].mom()+pim[0].mom()).mass()/GeV);
388 _h_4B_pipi[1] ->fill( (pi0[0].mom()+pi0[1].mom()).mass()/GeV);
389 _h_4B_pipi[1] ->fill( (pi0[0].mom()+pi0[2].mom()).mass()/GeV);
390 _h_4B_pipi[1] ->fill( (pi0[1].mom()+pi0[2].mom()).mass()/GeV);
391 _h_4B_pipipi[0] ->fill( (pi0[0].mom()+pi0[1].mom()+pi0[2].mom()).mass()/GeV);
392 _h_4B_pipipi[1] ->fill( (pi0[0].mom()+pi0[1].mom()+pim[0].mom()).mass()/GeV);
393 _h_4B_pipipi[1] ->fill( (pi0[0].mom()+pi0[2].mom()+pim[0].mom()).mass()/GeV);
394 _h_4B_pipipi[1] ->fill( (pi0[1].mom()+pi0[2].mom()+pim[0].mom()).mass()/GeV);
395 _h_4B_pipipipi[0] ->fill( (pi0[0].mom()+pi0[1].mom()+pi0[2].mom()+pim[0].mom()).mass()/GeV);
396 }
397 else if (pi0.size()==1 && pip.size()==1 && pim.size()==2) {
398 _h_4B_pipi[2] ->fill((pi0[0].mom()+pip[0].mom()).mass()/GeV);
399 _h_4B_pipi[3] ->fill((pi0[0].mom()+pim[0].mom()).mass()/GeV);
400 _h_4B_pipi[3] ->fill((pi0[0].mom()+pim[1].mom()).mass()/GeV);
401 _h_4B_pipi[4] ->fill((pip[0].mom()+pim[0].mom()).mass()/GeV);
402 _h_4B_pipi[4] ->fill((pip[0].mom()+pim[1].mom()).mass()/GeV);
403 _h_4B_pipi[5] ->fill((pim[0].mom()+pim[1].mom()).mass()/GeV);
404 _h_4B_pipipi[2] ->fill( (pi0[0].mom()+pip[0].mom()+pim[0].mom()).mass()/GeV);
405 _h_4B_pipipi[2] ->fill( (pi0[0].mom()+pip[0].mom()+pim[1].mom()).mass()/GeV);
406 _h_4B_pipipi[3] ->fill( (pip[0].mom()+pim[0].mom()+pim[1].mom()).mass()/GeV);
407 _h_4B_pipipi[4] ->fill( (pi0[0].mom()+pim[0].mom()+pim[1].mom()).mass()/GeV);
408 _h_4B_pipipipi[1] ->fill( (pi0[0].mom()+pip[0].mom()+pim[0].mom()+pim[1].mom()).mass()/GeV);
409 }
410 }
411 else if (nstable==6) {
412 // 2 pi0 2pi- pi+
413 if (pi0.size()==2 && pim.size()==2 && pip.size()==1) {
414 FourMomentum ptotal = pim[0].mom()+pim[1].mom() + pip[0].mom()+pi0[0].mom()+pi0[1].mom();
415 _h_5B_pipi1[0]->fill((pim[0].mom()+pim[1].mom()).mass()/GeV);
416 _h_5B_pipi1[1]->fill((pim[0].mom()+pip[0].mom()).mass()/GeV);
417 _h_5B_pipi1[1]->fill((pim[1].mom()+pip[0].mom()).mass()/GeV);
418 _h_5B_pipi1[2]->fill((pim[0].mom()+pi0[0].mom()).mass()/GeV);
419 _h_5B_pipi1[2]->fill((pim[0].mom()+pi0[1].mom()).mass()/GeV);
420 _h_5B_pipi1[2]->fill((pim[1].mom()+pi0[0].mom()).mass()/GeV);
421 _h_5B_pipi1[2]->fill((pim[1].mom()+pi0[1].mom()).mass()/GeV);
422 _h_5B_pipi1[3]->fill((pip[0].mom()+pi0[0].mom()).mass()/GeV);
423 _h_5B_pipi1[3]->fill((pip[0].mom()+pi0[1].mom()).mass()/GeV);
424 _h_5B_pipi1[4]->fill((pi0[0].mom()+pi0[1].mom()).mass()/GeV);
425 _h_5B_pipipi1[0]->fill((pim[0].mom()+pim[1].mom()-ptotal).mass()/GeV);
426 _h_5B_pipipi1[1]->fill((pim[0].mom()+pip[0].mom()-ptotal).mass()/GeV);
427 _h_5B_pipipi1[1]->fill((pim[1].mom()+pip[0].mom()-ptotal).mass()/GeV);
428 _h_5B_pipipi1[2]->fill((pim[0].mom()+pi0[0].mom()-ptotal).mass()/GeV);
429 _h_5B_pipipi1[2]->fill((pim[0].mom()+pi0[1].mom()-ptotal).mass()/GeV);
430 _h_5B_pipipi1[2]->fill((pim[1].mom()+pi0[0].mom()-ptotal).mass()/GeV);
431 _h_5B_pipipi1[2]->fill((pim[1].mom()+pi0[1].mom()-ptotal).mass()/GeV);
432 _h_5B_pipipi1[3]->fill((pip[0].mom()+pi0[0].mom()-ptotal).mass()/GeV);
433 _h_5B_pipipi1[3]->fill((pip[0].mom()+pi0[1].mom()-ptotal).mass()/GeV);
434 _h_5B_pipipi1[4]->fill((pi0[0].mom()+pi0[1].mom()-ptotal).mass()/GeV);
435 _h_5B_pipipipi1[0]->fill((ptotal-pim[0].mom()).mass()/GeV);
436 _h_5B_pipipipi1[0]->fill((ptotal-pim[1].mom()).mass()/GeV);
437 _h_5B_pipipipi1[1]->fill((ptotal-pip[0].mom()).mass()/GeV);
438 _h_5B_pipipipi1[2]->fill((ptotal-pi0[0].mom()).mass()/GeV);
439 _h_5B_pipipipi1[2]->fill((ptotal-pi0[1].mom()).mass()/GeV);
440 _h_5B_q1->fill(ptotal.mass()/GeV);
441 }
442 // 4 pi0 pi-
443 else if (pi0.size()==4 && pim.size()==1) {
444 FourMomentum ptotal = pi0[0].mom()+pi0[1].mom()+pi0[2].mom() + pi0[3].mom()+pim[0].mom();
445 _h_5B_pipi2[0]->fill((pim[0].mom()+pi0[0].mom()).mass()/GeV);
446 _h_5B_pipi2[0]->fill((pim[0].mom()+pi0[1].mom()).mass()/GeV);
447 _h_5B_pipi2[0]->fill((pim[0].mom()+pi0[2].mom()).mass()/GeV);
448 _h_5B_pipi2[0]->fill((pim[0].mom()+pi0[3].mom()).mass()/GeV);
449 _h_5B_pipi2[1]->fill((pi0[0].mom()+pi0[1].mom()).mass()/GeV);
450 _h_5B_pipi2[1]->fill((pi0[0].mom()+pi0[2].mom()).mass()/GeV);
451 _h_5B_pipi2[1]->fill((pi0[0].mom()+pi0[3].mom()).mass()/GeV);
452 _h_5B_pipi2[1]->fill((pi0[1].mom()+pi0[2].mom()).mass()/GeV);
453 _h_5B_pipi2[1]->fill((pi0[1].mom()+pi0[3].mom()).mass()/GeV);
454 _h_5B_pipi2[1]->fill((pi0[2].mom()+pi0[3].mom()).mass()/GeV);
455 _h_5B_pipipi2[0]->fill((pim[0].mom()+pi0[0].mom()-ptotal).mass()/GeV);
456 _h_5B_pipipi2[0]->fill((pim[0].mom()+pi0[1].mom()-ptotal).mass()/GeV);
457 _h_5B_pipipi2[0]->fill((pim[0].mom()+pi0[2].mom()-ptotal).mass()/GeV);
458 _h_5B_pipipi2[0]->fill((pim[0].mom()+pi0[3].mom()-ptotal).mass()/GeV);
459 _h_5B_pipipi2[1]->fill((pi0[0].mom()+pi0[1].mom()-ptotal).mass()/GeV);
460 _h_5B_pipipi2[1]->fill((pi0[0].mom()+pi0[2].mom()-ptotal).mass()/GeV);
461 _h_5B_pipipi2[1]->fill((pi0[0].mom()+pi0[3].mom()-ptotal).mass()/GeV);
462 _h_5B_pipipi2[1]->fill((pi0[1].mom()+pi0[2].mom()-ptotal).mass()/GeV);
463 _h_5B_pipipi2[1]->fill((pi0[1].mom()+pi0[3].mom()-ptotal).mass()/GeV);
464 _h_5B_pipipi2[1]->fill((pi0[2].mom()+pi0[3].mom()-ptotal).mass()/GeV);
465 _h_5B_pipipipi2[0]->fill((ptotal-pim[0].mom()).mass()/GeV);
466 _h_5B_pipipipi2[1]->fill((ptotal-pi0[0].mom()).mass()/GeV);
467 _h_5B_pipipipi2[1]->fill((ptotal-pi0[1].mom()).mass()/GeV);
468 _h_5B_pipipipi2[1]->fill((ptotal-pi0[2].mom()).mass()/GeV);
469 _h_5B_pipipipi2[1]->fill((ptotal-pi0[3].mom()).mass()/GeV);
470 _h_5B_q2->fill(ptotal.mass()/GeV);
471 }
472 // 3 pi- 2pi+
473 else if (pim.size()==3 && pip.size()==2) {
474 FourMomentum ptotal = pim[0].mom()+pim[1].mom() + pim[2].mom()+pip[0].mom()+pip[1].mom();
475 _h_5B_pipi3[0]->fill((pip[0].mom()+pip[1].mom()).mass()/GeV);
476 _h_5B_pipi3[1]->fill((pim[0].mom()+pip[0].mom()).mass()/GeV);
477 _h_5B_pipi3[1]->fill((pim[0].mom()+pip[1].mom()).mass()/GeV);
478 _h_5B_pipi3[1]->fill((pim[1].mom()+pip[0].mom()).mass()/GeV);
479 _h_5B_pipi3[1]->fill((pim[1].mom()+pip[1].mom()).mass()/GeV);
480 _h_5B_pipi3[1]->fill((pim[2].mom()+pip[0].mom()).mass()/GeV);
481 _h_5B_pipi3[1]->fill((pim[2].mom()+pip[1].mom()).mass()/GeV);
482 _h_5B_pipi3[2]->fill((pim[0].mom()+pim[1].mom()).mass()/GeV);
483 _h_5B_pipi3[2]->fill((pim[0].mom()+pim[2].mom()).mass()/GeV);
484 _h_5B_pipi3[2]->fill((pim[1].mom()+pim[2].mom()).mass()/GeV);
485 _h_5B_pipipi3[0]->fill((pip[0].mom()+pip[1].mom()-ptotal).mass()/GeV);
486 _h_5B_pipipi3[1]->fill((pim[0].mom()+pip[0].mom()-ptotal).mass()/GeV);
487 _h_5B_pipipi3[1]->fill((pim[0].mom()+pip[1].mom()-ptotal).mass()/GeV);
488 _h_5B_pipipi3[1]->fill((pim[1].mom()+pip[0].mom()-ptotal).mass()/GeV);
489 _h_5B_pipipi3[1]->fill((pim[1].mom()+pip[1].mom()-ptotal).mass()/GeV);
490 _h_5B_pipipi3[1]->fill((pim[2].mom()+pip[0].mom()-ptotal).mass()/GeV);
491 _h_5B_pipipi3[1]->fill((pim[2].mom()+pip[1].mom()-ptotal).mass()/GeV);
492 _h_5B_pipipi3[2]->fill((pim[0].mom()+pim[1].mom()-ptotal).mass()/GeV);
493 _h_5B_pipipi3[2]->fill((pim[0].mom()+pim[2].mom()-ptotal).mass()/GeV);
494 _h_5B_pipipi3[2]->fill((pim[1].mom()+pim[2].mom()-ptotal).mass()/GeV);
495 _h_5B_pipipipi3[0]->fill((ptotal-pim[0].mom()).mass()/GeV);
496 _h_5B_pipipipi3[0]->fill((ptotal-pim[1].mom()).mass()/GeV);
497 _h_5B_pipipipi3[0]->fill((ptotal-pim[2].mom()).mass()/GeV);
498 _h_5B_pipipipi3[1]->fill((ptotal-pip[0].mom()).mass()/GeV);
499 _h_5B_pipipipi3[1]->fill((ptotal-pip[1].mom()).mass()/GeV);
500 _h_5B_q3->fill(ptotal.mass()/GeV);
501 }
502 }
503 }
504 }
505
506
507 /// Normalise histograms etc., after the run
508 void finalize() {
509
510 // Leptonic
511 normalize(_h_2B_m2enu);
512 normalize(_h_2B_menu );
513
514 // 1 hadron
515 normalize(_h_1B_xpi);
516
517 // 2 hadrons
518 normalize(_h_2B_m2pipi);
519 normalize(_h_2B_mpipi );
520 normalize(_h_2B_m2munu);
521 normalize(_h_2B_mmunu );
522 normalize(_h_2B_m2KpiA);
523 normalize(_h_2B_mKpiA );
524 normalize(_h_2B_m2KpiB);
525 normalize(_h_2B_mKpiB );
526 normalize(_h_2B_m2Keta);
527 normalize(_h_2B_mKeta );
528 normalize(_h_2B_m2KK );
529 normalize(_h_2B_mKK );
530
531 // 3 hadrons
532 for (size_t ix = 0; ix < 4; ++ix) {
533 if (ix < 3) {
534 normalize(_h_3B_pippimpim [ix]);
535 normalize(_h_3B_pi0pi0pim [ix]);
536 normalize(_h_3B_pi0pi0km [ix]);
537 normalize(_h_3B_kspimks [ix]);
538 normalize(_h_3B_klpimkl [ix]);
539 }
540 normalize(_h_3B_kmpimkp [ix]);
541 normalize(_h_3B_kmpi0k0 [ix]);
542 normalize(_h_3B_kmpimpip [ix]);
543 normalize(_h_3B_pimk0pi0 [ix]);
544 normalize(_h_3B_pimpi0eta [ix]);
545 normalize(_h_3B_pimpi0gamma[ix]);
546 normalize(_h_3B_kspimkl [ix]);
547 }
548
549 // 4 pion decays
550 for (size_t ix=0;ix<5;++ix) {
551 normalize(_h_4B_pipi [ix]);
552 normalize(_h_4B_pipipi[ix]);
553 }
554 normalize(_h_4B_pipi[5]);
555 for (size_t ix=0;ix<2;++ix) {
556 normalize(_h_4B_pipipipi[ix]);
557 }
558
559 // 5 pions
560 normalize(_h_5B_q1);
561 for (size_t ix=0;ix<5;++ix) {
562 normalize(_h_5B_pipi1);
563 normalize(_h_5B_pipipi1);
564 }
565 for (size_t ix=0;ix<3;++ix) {
566 normalize(_h_5B_pipipipi1);
567 }
568
569 // 4 pi0 pi-
570 normalize(_h_5B_q2);
571 for (size_t ix=0;ix<2;++ix) {
572 normalize(_h_5B_pipi2);
573 normalize(_h_5B_pipipi2);
574 normalize(_h_5B_pipipipi2);
575 }
576
577 // 3 pi- 2 pi+
578 normalize(_h_5B_q3);
579 for (size_t ix=0;ix<3;++ix) {
580 normalize(_h_5B_pipi3);
581 normalize(_h_5B_pipipi3);
582 }
583 for (size_t ix=0;ix<2;++ix) {
584 normalize(_h_5B_pipipipi3);
585 }
586 }
587
588 /// @}
589
590
591 /// @name Histograms
592 /// @{
593
594 /// Histograms for leptonic decay
595 Histo1DPtr _h_2B_m2enu,_h_2B_menu;
596 Histo1DPtr _h_2B_m2munu,_h_2B_mmunu;
597
598 /// Histograms for 1 hadron decay
599 Histo1DPtr _h_1B_xpi;
600
601 /// Histograms for 2 hadron decay
602 Histo1DPtr _h_2B_m2pipi,_h_2B_mpipi;
603 Histo1DPtr _h_2B_m2KpiA,_h_2B_m2KpiB,_h_2B_mKpiA,_h_2B_mKpiB;
604 Histo1DPtr _h_2B_m2Keta,_h_2B_mKeta;
605 Histo1DPtr _h_2B_m2KK,_h_2B_mKK;
606
607 // Histograms for 3 hadron decay
608 /// Histograms for tau^- -> nu_tau pi^+pi^-pi^-
609 vector<Histo1DPtr> _h_3B_pippimpim;
610 /// Histograms for tau^- -> nu_tau pi^0pi^0pi^-
611 vector<Histo1DPtr> _h_3B_pi0pi0pim;
612 /// Histograms for tau^- -> nu_tau K^-K^+pi^-
613 vector<Histo1DPtr> _h_3B_kmpimkp;
614 /// Histograms for tau^- -> nu_tau K^-K^0pi^0
615 vector<Histo1DPtr> _h_3B_kmpi0k0;
616 /// Histograms for tau^- -> nu_tau pi^0pi^0K^-
617 vector<Histo1DPtr> _h_3B_pi0pi0km;
618 /// Histograms for tau^- -> nu_tau K^-pi^-pi^+
619 vector<Histo1DPtr> _h_3B_kmpimpip;
620 /// Histograms for tau^- -> nu_tau pi^-K^0pi^0
621 vector<Histo1DPtr> _h_3B_pimk0pi0;
622 /// Histograms for tau^- -> nu_tau pi^-pi^0eta
623 vector<Histo1DPtr> _h_3B_pimpi0eta;
624 /// Histograms for tau^- -> nu_tau pi^-pi^0gamma
625 vector<Histo1DPtr> _h_3B_pimpi0gamma;
626 /// Histograms for tau^- -> nu_tau K^0_SK^0_Spi^-
627 vector<Histo1DPtr> _h_3B_kspimks;
628 /// Histograms for tau^- -> nu_tau K^0_LK^0_Lpi^-
629 vector<Histo1DPtr> _h_3B_klpimkl;
630 /// Histograms for tau^- -> nu_tau K^0_SK^0_Lpi^-
631 vector<Histo1DPtr> _h_3B_kspimkl;
632
633 // Histograms for 4 pion decay
634 /// Histograms for the pipi mass distributions
635 vector<Histo1DPtr> _h_4B_pipi;
636 /// Histograms for the pipipi mass distributions
637 vector<Histo1DPtr> _h_4B_pipipi;
638 /// Histograms for the pipipipi mass distributions
639 vector<Histo1DPtr> _h_4B_pipipipi;
640
641 // Histograms for 5 pion decay
642 /// 2 pi0 2 pi- pi+
643 Histo1DPtr _h_5B_q1;
644 vector<Histo1DPtr> _h_5B_pipi1;
645 vector<Histo1DPtr> _h_5B_pipipi1;
646 vector<Histo1DPtr> _h_5B_pipipipi1;
647 /// 4 pi0 pi-
648 Histo1DPtr _h_5B_q2;
649 vector<Histo1DPtr> _h_5B_pipi2;
650 vector<Histo1DPtr> _h_5B_pipipi2;
651 vector<Histo1DPtr> _h_5B_pipipipi2;
652 /// 3 pi- 2 pi+
653 Histo1DPtr _h_5B_q3;
654 vector<Histo1DPtr> _h_5B_pipi3;
655 vector<Histo1DPtr> _h_5B_pipipi3;
656 vector<Histo1DPtr> _h_5B_pipipipi3;
657 /// @}
658
659 };
660
661
662 RIVET_DECLARE_PLUGIN(MC_DECAY_TAU);
663
664}
|