rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

MC_DECAY_OMEGAPHIA1_3PION

MC analysis of $\omega$, $\phi$ and $a_1\to3\pi$ decays
Experiment: ()
Status: VALIDATED
Authors:
  • Peter Richardson
No references listed
Beams: * *
Beam energies: ANY
Run details:
  • Any process producing a1, omega or phi mesons

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}