Rivet analyses referenceMC_DECAY_MESON_MESON_LEPTONSMC analyis of $V\to P \ell^+\ell^-$ and $P\to V \ell^+\ell^-$ decaysExperiment: () Status: VALIDATED Authors:
Beams: * * Beam energies: ANY Run details:
A Monte Carlo analysis for the decay of vector mesons to a pseudoscalar meson, or a pseudoscalar meson to a vector meson, and an $e^+e^-$ or $\mu^+\mu^-$ pair. All such decays are automatically identified. Based on old Herwig++ internal analysis. Source code: MC_DECAY_MESON_MESON_LEPTONS.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/UnstableParticles.hh"
4#include "Rivet/Tools/ParticleIdUtils.hh"
5
6namespace Rivet {
7
8
9 /// @brief MC decay M -> M l+ l-
10 class MC_DECAY_MESON_MESON_LEPTONS : public Analysis {
11 public:
12
13 /// Constructor
14 RIVET_DEFAULT_ANALYSIS_CTOR(MC_DECAY_MESON_MESON_LEPTONS);
15
16
17 /// @name Analysis methods
18 /// @{
19
20 /// Book histograms and initialise projections before the run
21 void init() {
22
23 // Initialise and register projections
24 declare(UnstableParticles(),"UFS");
25
26 // book the Histos
27 // pi0 dalitz
28 bookHistos(111,22,11,0.140);
29 // eta dalitz
30 bookHistos(221,22,11,0.55);
31 bookHistos(221,22,13,0.55);
32 // eta' dalitz
33 bookHistos(331,22,11,0.96);
34 bookHistos(331,22,13,0.96);
35 // omega -> pi0
36 bookHistos(223,111,11,0.8);
37 bookHistos(223,111,13,0.8);
38 // phi -> pi0
39 bookHistos(333,111,11,1.1);
40 bookHistos(333,111,13,1.1);
41 // phi -> eta
42 bookHistos(333,221,11,1.1);
43 bookHistos(333,221,13,1.1);
44 // J/psi dalitz
45 bookHistos(443,22,11,3.1);
46 bookHistos(443,22,13,3.1);
47 // B -> s gamma
48 bookHistos(511,313,11,5.3);
49 bookHistos(511,313,13,5.3);
50 }
51
52 void bookHistos(int id1, int id2, int il, double dM) {
53 if(abs(id2)%10==3 || id2==22) {
54 _incoming_P.push_back(id1);
55 _outgoingV.push_back(id2);
56 _outgoingf_P.push_back(il);
57 std::ostringstream title;
58 title << "h2_" << abs(id1);
59 if(id1>0) title << "p";
60 else title << "m";
61 title << "_" << abs(id2);
62 if(id2>0) title << "p";
63 else title << "m";
64 title << "_" << il << "_";
65 _mff_P .push_back(Histo1DPtr());
66 book(_mff_P.back(), title.str()+"mff" , 100, 0., dM);
67 _mVf .push_back(Histo1DPtr());
68 book(_mVf.back(), title.str()+"mVf" , 100, 0., dM);
69 _mVfbar.push_back(Histo1DPtr());
70 book(_mVfbar.back(), title.str()+"mVfbar", 100, 0., dM);
71 }
72 else {
73 _incomingV.push_back(id1);
74 _outgoingP.push_back(id2);
75 _outgoingf_V.push_back(il);
76 std::ostringstream title;
77 title << "h_" << abs(id1);
78 if(id1>0) title << "p";
79 else title << "m";
80 title << "_" << abs(id2);
81 if(id2>0) title << "p";
82 else title << "m";
83 title << "_" << il << "_";
84 _mff_V .push_back(Histo1DPtr());
85 book(_mff_V.back(), title.str()+"mff" , 100, 0., dM);
86 _mPf .push_back(Histo1DPtr());
87 book(_mPf.back(), title.str()+"mPf" , 100, 0., dM);
88 _mPfbar.push_back(Histo1DPtr());
89 book(_mPfbar.back(), title.str()+"mPfbar", 100, 0., dM);
90 }
91 }
92
93 void findDecayProducts(const Particle & mother,
94 unsigned int & nstable,
95 Particles& lp, Particles& lm,
96 Particles& scalar,
97 Particles& vector) {
98 for(const Particle & p : mother.children()) {
99 int id = p.pid();
100 if ( id == PID::EMINUS || id == PID::MUON ) {
101 lm.push_back(p);
102 ++nstable;
103 }
104 else if (id == PID::EPLUS || id == PID::ANTIMUON) {
105 lp.push_back(p);
106 ++nstable;
107 }
108 else if (abs(id)%10==1 && PID::isMeson(id)) {
109 scalar.push_back(p);
110 ++nstable;
111 }
112 else if ((abs(id)%10==3 && PID::isMeson(id)) ||
113 id==PID::PHOTON ) {
114 vector.push_back(p);
115 ++nstable;
116 }
117 else if ( !p.children().empty() ) {
118 findDecayProducts(p,nstable,lp,lm,scalar,vector);
119 }
120 else
121 ++nstable;
122 }
123 }
124
125 /// Perform the per-event analysis
126 void analyze(const Event& event) {
127 // loop over unstable particles
128 for(const Particle& iMeson : apply<UnstableParticles>(event, "UFS").particles()) {
129 // only consider scalar/vector mesons
130 long pid = iMeson.pid();
131 if(!PID::isMeson(pid)) continue;
132 if(abs(pid)%10!=3 and abs(pid)%10!=1 ) continue;
133 Particles lp,lm,scalar,vector;
134 unsigned int nstable(0);
135 findDecayProducts(iMeson,nstable,lp,lm,scalar,vector);
136 if(nstable!=3 || lp.size()!=1 || lm.size()!=1 || lp[0].pid()!=-lm[0].pid()) continue;
137 if(scalar.size()==1) {
138 // check if we already have this decay
139 unsigned int ix=0; bool found(false);
140 while(!found&&ix<_incomingV.size()) {
141 if(_incomingV[ix]==pid && _outgoingP[ix]==scalar[0].pid() &&
142 _outgoingf_V[ix]==lm[0].pid()) {
143 found=true;
144 }
145 else {
146 ++ix;
147 }
148 }
149 // create a new graph if needed
150 if(!found) {
151 MSG_WARNING("MC_DECAY_MESON_MESON_LEPTONS S" << iMeson.pid() << " " << scalar[0].pid() << " "
152 << iMeson.mass() << "\n");
153 continue;
154 }
155 // add the results to the histogram
156 _mff_V [ix]->fill((lm [0].momentum()+lp[0].momentum()).mass());
157 _mPf [ix]->fill((scalar[0].momentum()+lm[0].momentum()).mass());
158 _mPfbar[ix]->fill((scalar[0].momentum()+lp[0].momentum()).mass());
159 }
160 else if(vector.size()==1) {
161 // check if we already have this decay
162 unsigned int ix=0; bool found(false);
163 while(!found&&ix<_incoming_P.size()) {
164 if(_incoming_P[ix]==pid && _outgoingV[ix]==vector[0].pid() &&
165 _outgoingf_P[ix]==lm[0].pid()) {
166 found=true;
167 }
168 else {
169 ++ix;
170 }
171 }
172 // create a new graph if needed
173 if(!found) {
174 MSG_WARNING("MC_DECAY_MESON_MESON_LEPTONS V" << iMeson.pid() << " " << vector[0].pid() << " "
175 << iMeson.mass() << "\n");
176 continue;
177 }
178 // add the results to the histogram
179 _mff_P [ix]->fill((lm [0].momentum()+lp[0].momentum()).mass());
180 _mVf [ix]->fill((vector[0].momentum()+lm[0].momentum()).mass());
181 _mVfbar[ix]->fill((vector[0].momentum()+lp[0].momentum()).mass());
182 }
183 }
184 }
185
186 /// Normalise histograms etc., after the run
187 void finalize() {
188
189 // normalize to unity V->P
190 for(unsigned int ix=0;ix<_mff_V.size();++ix) {
191 normalize(_mff_V);
192 normalize(_mPf);
193 normalize(_mPfbar);
194 }
195 // normalize to unity P->V
196 for(unsigned int ix=0;ix<_mff_P.size();++ix) {
197 normalize(_mff_P);
198 normalize(_mVf);
199 normalize(_mVfbar);
200 }
201
202 }
203
204 /// @}
205
206
207
208 /// @name Histograms for V -> P
209 /// @{
210 /**
211 * PDG codes of the incoming particles
212 */
213 vector<long> _incomingV;
214
215 /**
216 * PDG codes of the outgoing pseudoscalar mesons
217 */
218 vector<long> _outgoingP;
219
220 /**
221 * PDG codes of the outgoing fermion
222 */
223 vector<long> _outgoingf_V;
224
225 /**
226 * Histograms for the mass of the fermion-antifermion pair
227 */
228 vector<Histo1DPtr> _mff_V;
229
230 /**
231 * Histograms for the masses of the pseudoscalar and the fermion
232 */
233 vector<Histo1DPtr> _mPf;
234
235 /**
236 * Histograms for the masses of the pseudoscalar and the antifermion
237 */
238 vector<Histo1DPtr> _mPfbar;
239 /// @}
240
241 /// @name Histograms P->V
242 /// @{
243 /**
244 * PDG codes of the incoming_P particles
245 */
246 vector<long> _incoming_P;
247
248 /**
249 * PDG codes of the outgoing vector mesons
250 */
251 vector<long> _outgoingV;
252
253 /**
254 * PDG codes of the outgoing fermion
255 */
256 vector<long> _outgoingf_P;
257
258 /**
259 * Histograms for the mass of the fermion-antifermion pair
260 */
261 vector<Histo1DPtr> _mff_P;
262
263 /**
264 * Histograms for the masses of the vector and the fermion
265 */
266 vector<Histo1DPtr> _mVf;
267
268 /**
269 * Histograms for the masses of the vector and the antifermion
270 */
271 vector<Histo1DPtr> _mVfbar;
272 /// @}
273
274
275 };
276
277
278 RIVET_DECLARE_PLUGIN(MC_DECAY_MESON_MESON_LEPTONS);
279
280}
|