rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

MC_OmegaPhia1_3Pion_Decay

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