rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

MC_TAU_Decay

Analysis of Kinematic distributions for $\tau$ lepton decays
Experiment: ()
Status: VALIDATED
Authors:
  • Peter Richardson
No references listed
Beams: * *
Beam energies: ANY
Run details:
  • Any type of process producing tau leptons

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_TAU_Decay.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_TAU_Decay : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(MC_TAU_Decay);
 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_TAU_Decay);
663
664}