]> git.uio.no Git - u/mrichter/AliRoot.git/blame - EVGEN/AliDecayerExodus.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Moving required CMake version from 2.8.4 to 2.8.8
[u/mrichter/AliRoot.git] / EVGEN / AliDecayerExodus.cxx
CommitLineData
b5c4afc6 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16
17#include "AliDecayerExodus.h"
18#include <Riostream.h>
19#include <TMath.h>
20#include <AliLog.h>
21#include <TH1.h>
22#include <TRandom.h>
23#include <TParticle.h>
24#include <TDatabasePDG.h>
25#include <TPDGCode.h>
26#include <TLorentzVector.h>
27#include <TClonesArray.h>
28
29
30ClassImp(AliDecayerExodus)
31
76e6e4c5 32//---------------------------------------------------------------------------------------------------
33//
34// Generate electron-pair mass distributions for Dalitz decays according
35// to the Kroll-Wada parametrization: N. Kroll, W. Wada: Phys. Rev 98(1955)1355
36// and generate electron-pair mass distributions for resonances according
37// to the Gounaris-Sakurai parametrization: G.J. Gounaris, J.J. Sakurai: Phys.Rev.Lett. 21(1968)244
38//
39// For the electromagnetic form factor the parameterization from
40// Lepton-G is used: L.G. Landsberg et al.: Phys. Rep. 128(1985)301
41//
42// Ralf Averbeck (R.Averbeck@gsi.de)
43// Irem Erdemir (irem.erdemir@cern.ch)
44//
45//---------------------------------------------------------------------------------------------------
46
47
b5c4afc6 48AliDecayerExodus::AliDecayerExodus():
49 AliDecayer(),
50 fEPMassPion(0),
51 fEPMassEta(0),
52 fEPMassEtaPrime(0),
53 fEPMassRho(0),
54 fEPMassOmega(0),
55 fEPMassOmegaDalitz(0),
56 fEPMassPhi(0),
57 fEPMassPhiDalitz(0),
58 fEPMassJPsi(0),
59 fInit(0)
60
61{
62// Constructor
63}
64
65
66void AliDecayerExodus::Init()
67{
68
69// Initialisation
70//
71 Int_t ibin, nbins;
72 Double_t min, maxpion, maxeta, maxomega, maxetaprime, maxphi, binwidth_pion, binwidth_eta, binwidth_omega, binwidth_etaprime, binwidth_phi;
73 Double_t pionmass, etamass, omegamass, etaprimemass, phimass, emass, omasspion, omasseta, omassgamma;
74 Double_t epsilon_pion, epsilon_eta, epsilon_omega, epsilon_etaprime, epsilon_phi;
75 Double_t delta_pion, delta_eta, delta_omega, delta_etaprime, delta_phi;
76 Double_t mLL_pion, mLL_eta, mLL_omega, mLL_etaprime, mLL_phi;
77 Double_t q_pion, q_eta, q_omega, q_etaprime, q_phi;
78 Double_t kwHelp_pion, kwHelp_eta, kwHelp_omega, kwHelp_etaprime, kwHelp_phi;
79 Double_t krollWada_pion, krollWada_eta, krollWada_omega, krollWada_etaprime, krollWada_phi;
80 Double_t formFactor_pion, formFactor_eta, formFactor_omega, formFactor_etaprime, formFactor_phi;
81 Double_t weight_pion, weight_eta, weight_omega_dalitz, weight_etaprime, weight_phi_dalitz;
82
83 Float_t binwidth;
84 Float_t mass_bin, mass_min, mass_max;
85 Double_t vmass_rho, vmass_omega, vmass_phi, vmass_jpsi, vwidth_rho, vwidth_omega, vwidth_phi, vwidth_jpsi;
86 Double_t weight_rho, weight_omega, weight_phi, weight_jpsi;
87
88//================================================================================//
89// Create electron pair mass histograms from dalitz decays //
90//================================================================================//
91
92 // Get the particle masses
93 // parent
94 nbins = 1000;
df2fcbcb 95 mass_min = 0.;
96 mass_max = 0.;
b5c4afc6 97 pionmass = (TDatabasePDG::Instance()->GetParticle(111))->Mass();
98 etamass = (TDatabasePDG::Instance()->GetParticle(221))->Mass();
99 omegamass = (TDatabasePDG::Instance()->GetParticle(223))->Mass();
100 etaprimemass = (TDatabasePDG::Instance()->GetParticle(331))->Mass();
101 phimass = (TDatabasePDG::Instance()->GetParticle(333))->Mass();
102 // child - electron
103 emass = (TDatabasePDG::Instance()->GetParticle(11))->Mass();
104 // child - other : third childs from Dalitz decays
105 omasspion = pionmass;
106 omasseta = etamass;
107 omassgamma = 0.;
108
109 min = 2.0 * emass;
110 maxpion = pionmass - omassgamma;
111 maxeta = etamass - omassgamma;
112 maxomega = omegamass - pionmass;
113 maxetaprime = etaprimemass - omassgamma;
114 maxphi = phimass - omasseta;
115
116 binwidth_pion = (maxpion - min) / (Double_t)nbins;
117 binwidth_eta = (maxeta - min) / (Double_t)nbins;
118 binwidth_omega = (maxomega - min) / (Double_t)nbins;
119 binwidth_etaprime = (maxetaprime - min) / (Double_t)nbins;
120 binwidth_phi = (maxphi - min) / (Double_t)nbins;
121
122
123 epsilon_pion = (emass / pionmass) * (emass / pionmass);
124 epsilon_eta = (emass / etamass) * (emass / etamass);
125 epsilon_omega = (emass / omegamass) * (emass / omegamass);
126 epsilon_etaprime = (emass / etaprimemass) * (emass / etaprimemass);
127 epsilon_phi = (emass / phimass) * (emass / phimass);
128
129
130 delta_pion = (omassgamma / pionmass) * (omassgamma / pionmass);
131 delta_eta = (omassgamma / etamass) * (omassgamma / etamass);
132 delta_omega = (omasspion / omegamass) * (omasspion / omegamass);
133 delta_etaprime = (omassgamma / etaprimemass) * (omassgamma / etaprimemass);
134 delta_phi = (omasseta / phimass) * (omasseta / phimass);
135
136
137
138 // create pair mass histograms for Dalitz decays of pi0, eta, omega, eta' and phi
139 if (!fEPMassPion) {delete fEPMassPion; fEPMassPion = new TH1F("fEPMassPion", "Dalitz electron pair from pion", nbins, min, maxpion); }
140 if (!fEPMassEta) {delete fEPMassEta; fEPMassEta = new TH1F("fEPMassEta", "Dalitz electron pair from eta", nbins, min, maxeta);}
141 if (!fEPMassOmegaDalitz) {delete fEPMassOmegaDalitz; fEPMassOmegaDalitz = new TH1F("fEPMassOmegaDalitz", "Dalitz electron pair from omega ", nbins, min, maxomega);}
142 if (!fEPMassEtaPrime) {delete fEPMassEtaPrime; fEPMassEtaPrime = new TH1F("fEPMassEtaPrime", "Dalitz electron pair from etaprime", nbins, min, maxetaprime);}
143 if (!fEPMassPhiDalitz) {delete fEPMassPhiDalitz; fEPMassPhiDalitz = new TH1F("fEPMassPhiDalitz", "Dalitz electron pair from phi ", nbins, min, maxphi);}
144
145
146 mLL_pion = mLL_eta = mLL_omega = mLL_etaprime = mLL_phi = 0.;
147
148 for (ibin = 1; ibin <= nbins; ibin++ )
149 {
150 mLL_pion = min + (Double_t)(ibin - 1) * binwidth_pion + binwidth_pion / 2.0;
151 mLL_eta = min + (Double_t)(ibin - 1) * binwidth_eta + binwidth_eta / 2.0;
152 mLL_omega = min + (Double_t)(ibin - 1) * binwidth_omega + binwidth_omega / 2.0;
153 mLL_etaprime = min + (Double_t)(ibin - 1) * binwidth_etaprime + binwidth_etaprime / 2.0;
154 mLL_phi = min + (Double_t)(ibin - 1) * binwidth_phi + binwidth_phi / 2.0;
155
156
157 q_pion = (mLL_pion / pionmass) * (mLL_pion / pionmass);
158 q_eta = (mLL_eta / etamass) * (mLL_eta / etamass);
159 q_omega = (mLL_omega / omegamass)*(mLL_omega / omegamass);
160 q_etaprime = (mLL_etaprime / etaprimemass) * (mLL_etaprime / etaprimemass);
161 q_phi = (mLL_phi / phimass) * (mLL_phi / phimass);
162
163 if ( q_pion <= 4.0 * epsilon_pion || q_eta <= 4.0 * epsilon_eta || q_omega <= 4.0 * epsilon_omega || q_etaprime <= 4.0 * epsilon_etaprime || q_phi <= 4.0 * epsilon_phi )
164 {
165 AliFatal("Error in calculating Dalitz mass histogram binning!");
166 }
167
168
169 kwHelp_pion = (1.0 + q_pion / (1.0 - delta_pion)) * (1.0 + q_pion / (1.0 - delta_pion))
170 - 4.0 * q_pion / ((1.0 - delta_pion) * (1.0 - delta_pion));
171
172 kwHelp_eta = (1.0 + q_eta / (1.0 - delta_eta)) * (1.0 + q_eta / (1.0 - delta_eta))
173 - 4.0 * q_eta / ((1.0 - delta_eta) * (1.0 - delta_eta));
174
175 kwHelp_omega = (1.0 + q_omega / (1.0 - delta_omega)) * (1.0 + q_omega / (1.0 - delta_omega))
176 - 4.0 * q_omega / ((1.0 - delta_omega) * (1.0 - delta_omega));
177
178 kwHelp_etaprime = (1.0 + q_etaprime / (1.0 - delta_etaprime)) * (1.0 + q_etaprime / (1.0 - delta_etaprime))
179 - 4.0 * q_etaprime / ((1.0 - delta_etaprime) * (1.0 - delta_etaprime));
180
181 kwHelp_phi = (1.0 + q_phi / (1.0 - delta_phi)) * (1.0 + q_phi / (1.0 - delta_phi))
182 - 4.0 * q_phi / ((1.0 - delta_phi) * (1.0 - delta_phi));
183
184
185
186
187 if ( kwHelp_pion <= 0.0 || kwHelp_eta <= 0.0 || kwHelp_omega <= 0.0 || kwHelp_etaprime <= 0.0 || kwHelp_phi <= 0.0 )
188 {
189 AliFatal("Error in calculating Dalitz mass histogram binning!");
190
191 }
192
193
194 // Invariant mass distributions of electron pairs from Dalitz decays
195 // using Kroll-Wada function
196
197 krollWada_pion = (2.0 / mLL_pion) * TMath::Exp(1.5 * TMath::Log(kwHelp_pion))
198 * TMath::Sqrt(1.0 - 4.0 * epsilon_pion / q_pion)
199 * (1.0 + 2.0 * epsilon_pion / q_pion);
200
201
202 krollWada_eta = (2.0 / mLL_eta) * TMath::Exp(1.5 * TMath::Log(kwHelp_eta))
203 * TMath::Sqrt(1.0 - 4.0 * epsilon_eta / q_eta)
204 * (1.0 + 2.0 * epsilon_eta / q_eta);
205
206
207 krollWada_omega = (2.0 / mLL_omega) * TMath::Exp(1.5 * TMath::Log(kwHelp_omega))
208 * TMath::Sqrt(1.0 - 4.0 * epsilon_omega / q_omega)
209 * (1.0 + 2.0 * epsilon_omega / q_omega);
210
211
212 krollWada_etaprime = (2.0 / mLL_etaprime) * TMath::Exp(1.5 * TMath::Log(kwHelp_etaprime))
213 * TMath::Sqrt(1.0 - 4.0 * epsilon_etaprime / q_etaprime)
214 * (1.0 + 2.0 * epsilon_etaprime / q_etaprime);
215
216 krollWada_phi = (2.0 / mLL_phi) * TMath::Exp(1.5 * TMath::Log(kwHelp_phi))
217 * TMath::Sqrt(1.0 - 4.0 * epsilon_phi / q_phi)
218 * (1.0 + 2.0 * epsilon_phi / q_phi);
219
220
221
222 // Form factors from Lepton-G
223 formFactor_pion = 1.0/(1.0-5.5*mLL_pion*mLL_pion);
224 formFactor_eta = 1.0/(1.0-1.9*mLL_eta*mLL_eta);
225 formFactor_omega = (TMath::Power(TMath::Power(0.6519,2),2))
226 /(TMath::Power(TMath::Power(0.6519,2)-TMath::Power(mLL_omega, 2), 2)
227 + TMath::Power(0.04198, 2)*TMath::Power(0.6519, 2));
228 formFactor_etaprime = (TMath::Power(TMath::Power(0.764,2),2))
229 /(TMath::Power(TMath::Power(0.764,2)-TMath::Power(mLL_etaprime, 2), 2)
230 + TMath::Power(0.1020, 2)*TMath::Power(0.764, 2));
231 formFactor_phi = 1.0;
232
233
234
235
236 weight_pion = krollWada_pion * formFactor_pion * formFactor_pion;
237 weight_eta = krollWada_eta * formFactor_eta * formFactor_eta;
238 weight_omega_dalitz = krollWada_omega * formFactor_omega;
239 weight_etaprime = krollWada_etaprime * formFactor_etaprime;
240 weight_phi_dalitz = krollWada_phi * formFactor_phi * formFactor_phi;
241
242
243 // Fill histograms of electron pair masses from dalitz decays
244 fEPMassPion ->AddBinContent(ibin, weight_pion);
245 fEPMassEta ->AddBinContent(ibin, weight_eta);
246 fEPMassOmegaDalitz->AddBinContent(ibin, weight_omega_dalitz);
247 fEPMassEtaPrime ->AddBinContent(ibin, weight_etaprime);
248 fEPMassPhiDalitz ->AddBinContent(ibin, weight_phi_dalitz);
249 }
250
251
252
253
254//===================================================================================//
255// Create electron pair mass histograms from resonance decays //
256//===================================================================================//
257
258 Double_t pimass = 0.13956995;
259
260 // Get the particle masses
261 // parent
262 vmass_rho = (TDatabasePDG::Instance()->GetParticle(113))->Mass();
263 vmass_omega = (TDatabasePDG::Instance()->GetParticle(223))->Mass();
264 vmass_phi = (TDatabasePDG::Instance()->GetParticle(333))->Mass();
265 vmass_jpsi = (TDatabasePDG::Instance()->GetParticle(443))->Mass();
266 // Get the particle widths
267 // parent
268 vwidth_rho = (TDatabasePDG::Instance()->GetParticle(113))->Width();
269 vwidth_omega = (TDatabasePDG::Instance()->GetParticle(223))->Width();
270 vwidth_phi = (TDatabasePDG::Instance()->GetParticle(333))->Width();
271 vwidth_jpsi = (TDatabasePDG::Instance()->GetParticle(443))->Width();
272
273
274 if ( mass_min == 0. && mass_max == 0. )
275 {
276 mass_min = 2.*pimass;
277 mass_max = 5.;
278 }
279
280 binwidth = (mass_max-mass_min)/(Double_t)nbins;
281
282 // create pair mass histograms for resonances of rho, omega, phi and jpsi
283 if (!fEPMassRho) {delete fEPMassRho; fEPMassRho = new TH1F("fEPMassRho","mass rho",nbins,mass_min,mass_max);}
284 if (!fEPMassOmega) {delete fEPMassOmega; fEPMassOmega = new TH1F("fEPMassOmega","mass omega",nbins,mass_min,mass_max);}
285 if (!fEPMassPhi) {delete fEPMassPhi; fEPMassPhi = new TH1F("fEPMassPhi","mass phi",nbins,mass_min,mass_max);}
286 if (!fEPMassJPsi) {delete fEPMassJPsi; fEPMassJPsi = new TH1F("fEPMassJPsi","mass jpsi",nbins,mass_min,mass_max);}
287
288
289 for (ibin=1; ibin<=nbins; ibin++ )
290 {
291 mass_bin = mass_min+(Double_t)(ibin-1)*binwidth+binwidth/2.0;
292
293 weight_rho = (Float_t)GounarisSakurai(mass_bin,vmass_rho,vwidth_rho,emass);
294 weight_omega = (Float_t)GounarisSakurai(mass_bin,vmass_omega,vwidth_omega,emass);
295 weight_phi = (Float_t)GounarisSakurai(mass_bin,vmass_phi,vwidth_phi,emass);
76e6e4c5 296 weight_jpsi = (Float_t)Lorentz(mass_bin,vmass_jpsi,vwidth_jpsi);
b5c4afc6 297
298 // Fill histograms of electron pair masses from resonance decays
299 fEPMassRho ->AddBinContent(ibin,weight_rho);
300 fEPMassOmega->AddBinContent(ibin,weight_omega);
301 fEPMassPhi ->AddBinContent(ibin,weight_phi);
302 fEPMassJPsi ->AddBinContent(ibin,weight_jpsi);
303 }
304
305}
306
307Double_t AliDecayerExodus::GounarisSakurai(Float_t mass, Double_t vmass, Double_t vwidth, Double_t emass)
308{
309// Invariant mass distributions of electron pairs from resonance decays
76e6e4c5 310// of rho, omega and phi
b5c4afc6 311// using Gounaris-Sakurai function
312
313 Double_t corr = 0.;
314 Double_t epsilon = 0.;
315 Double_t weight = 0.;
316
317 Double_t pimass = 0.13956995;
76e6e4c5 318
5abb8148 319 if(mass>pimass){
76e6e4c5 320 corr = vwidth*(vmass/mass)*exp(1.5*log((mass*mass/4.0-pimass*pimass)
321 /(vmass*vmass/4.0-pimass*pimass)));
5abb8148 322 }
b5c4afc6 323
b5c4afc6 324 epsilon = (emass/mass)*(emass/mass);
325
326 if ( 1.0-4.0*epsilon>=0.0 )
327 {
328 weight = sqrt(1.0-4.0*epsilon)*(1.0+2.0*epsilon)/
329 ((vmass*vmass-mass*mass)*(vmass*vmass-mass*mass)+
330 (vmass*corr)*(vmass*corr));
331 }
332 return weight;
333}
334
335
76e6e4c5 336Double_t AliDecayerExodus::Lorentz(Float_t mass, Double_t vmass, Double_t vwidth)
337{
338// Invariant mass distributions of electron pairs from resonance decay
339// of jpsi (and it can also be used for other particles except rho, omega and phi)
340// using Lorentz function
341
342 Double_t weight;
343
344 weight = (vwidth*vwidth/4.0)/(vwidth*vwidth/4.0+(vmass-mass)*(vmass-mass));
345
346 return weight;
347
348}
349
b5c4afc6 350void AliDecayerExodus::Decay(Int_t idpart, TLorentzVector* pparent)
351{
352
353 if (!fInit) {
354 Init();
355 fInit=1;
356 }
357
358
359 Double_t pmass_pion, pmass_eta, pmass_omega_dalitz, pmass_etaprime, pmass_phi_dalitz;
360 Double_t emass, omass_pion, omass_eta, omass_gamma, epmass_pion, epmass_eta, epmass_omega_dalitz, epmass_etaprime, epmass_phi_dalitz;
361 Double_t e1_pion, e1_eta, e1_omega, e1_etaprime, e1_phi;
362 Double_t p1_pion, p1_eta, p1_omega, p1_etaprime, p1_phi;
363 Double_t e3_gamma_pion, e3_gamma_eta, e3_pion, e3_gamma_etaprime, e3_eta;
364 Double_t p3_gamma_pion, p3_gamma_eta, p3_pion, p3_gamma_etaprime, p3_eta;
365
366 Double_t wp_rho, wp_omega, wp_phi, wp_jpsi, epmass_rho, epmass_omega, epmass_phi, epmass_jpsi;
367 Double_t mp_rho, mp_omega, mp_phi, mp_jpsi, md_rho, md_omega, md_phi, md_jpsi;
368 Double_t Ed_rho, Ed_omega, Ed_phi, Ed_jpsi, pd_rho, pd_omega, pd_phi, pd_jpsi;
369
370
371 md_rho = md_omega = md_phi = md_jpsi = 0.;
372
373
374 Double_t costheta, sintheta, cosphi, sinphi, phi;
375
376 // Get the particle masses of daughters
377 emass = (TDatabasePDG::Instance()->GetParticle(11)) ->Mass();
378 omass_pion = (TDatabasePDG::Instance()->GetParticle(111))->Mass();
379 omass_eta = (TDatabasePDG::Instance()->GetParticle(221))->Mass();
380 omass_gamma = (TDatabasePDG::Instance()->GetParticle(22)) ->Mass();
381
382 // Get the particle widths of mothers for resonances
383 wp_rho = (TDatabasePDG::Instance()->GetParticle(113))->Width();
384 wp_omega = (TDatabasePDG::Instance()->GetParticle(223))->Width();
385 wp_phi = (TDatabasePDG::Instance()->GetParticle(333))->Width();
386 wp_jpsi = (TDatabasePDG::Instance()->GetParticle(443))->Width();
387
388 costheta = (2.0 * gRandom->Rndm()) - 1.;
389 sintheta = TMath::Sqrt((1. + costheta) * (1. - costheta));
390 phi = 2.0 * TMath::ACos(-1.) * gRandom->Rndm();
391 sinphi = TMath::Sin(phi);
392 cosphi = TMath::Cos(phi);
393
394
395//-----------------------------------------------------------------------------//
396// Generate Pizero Dalitz decay //
397//-----------------------------------------------------------------------------//
398
399 if(idpart==111){
400 pmass_pion = pparent->M();
401
402 for(;;){
403 // Sample the electron pair mass from a histogram
404 epmass_pion = fEPMassPion->GetRandom();
405 if (pmass_pion-omass_gamma>epmass_pion && epmass_pion/2.>emass) break;
406 }
407
408 // electron pair kinematics in virtual photon rest frame
409 e1_pion = epmass_pion / 2.;
410 p1_pion = TMath::Sqrt((e1_pion + emass) * (e1_pion - emass));
411
412 // momentum vectors of electrons in virtual photon rest frame
413 Double_t pProd1_pion[3] = {p1_pion * sintheta * cosphi,
414 p1_pion * sintheta * sinphi,
415 p1_pion * costheta};
416
417 Double_t pProd2_pion[3] = {-1.0 * p1_pion * sintheta * cosphi,
418 -1.0 * p1_pion * sintheta * sinphi,
419 -1.0 * p1_pion * costheta};
420
421
422 // third child kinematics in parent meson rest frame
423 e3_gamma_pion = (pmass_pion * pmass_pion - epmass_pion * epmass_pion)/(2. * pmass_pion);
424 p3_gamma_pion = TMath::Sqrt((e3_gamma_pion * e3_gamma_pion));
425
426
427 // third child 4-vector in parent meson rest frame
428 fProducts_pion[2].SetPx(p3_gamma_pion * sintheta * cosphi);
429 fProducts_pion[2].SetPy(p3_gamma_pion * sintheta * sinphi);
430 fProducts_pion[2].SetPz(p3_gamma_pion * costheta);
431 fProducts_pion[2].SetE(e3_gamma_pion);
432
433
434 // electron 4-vectors in properly rotated virtual photon rest frame
435 Double_t pRot1_pion[3] = {0.};
436 Rot(pProd1_pion, pRot1_pion, costheta, -sintheta, -cosphi, -sinphi);
437 Double_t pRot2_pion[3] = {0.};
438 Rot(pProd2_pion, pRot2_pion, costheta, -sintheta, -cosphi, -sinphi);
439 fProducts_pion[0].SetPx(pRot1_pion[0]);
440 fProducts_pion[0].SetPy(pRot1_pion[1]);
441 fProducts_pion[0].SetPz(pRot1_pion[2]);
442 fProducts_pion[0].SetE(e1_pion);
443 fProducts_pion[1].SetPx(pRot2_pion[0]);
444 fProducts_pion[1].SetPy(pRot2_pion[1]);
445 fProducts_pion[1].SetPz(pRot2_pion[2]);
446 fProducts_pion[1].SetE(e1_pion);
447
448 // boost the dielectron into the parent meson's rest frame
449 Double_t eLPparent_pion = TMath::Sqrt(p3_gamma_pion * p3_gamma_pion + epmass_pion * epmass_pion);
450 TVector3 boostPair_pion( -1.0 * fProducts_pion[2].Px() / eLPparent_pion,
451 -1.0 * fProducts_pion[2].Py() / eLPparent_pion,
452 -1.0 * fProducts_pion[2].Pz() / eLPparent_pion);
453 fProducts_pion[0].Boost(boostPair_pion);
454 fProducts_pion[1].Boost(boostPair_pion);
455
456 // boost all decay products into the lab frame
457 TVector3 boostLab_pion(pparent->Px() / pparent->E(),
458 pparent->Py() / pparent->E(),
459 pparent->Pz() / pparent->E());
460
461 fProducts_pion[0].Boost(boostLab_pion);
462 fProducts_pion[1].Boost(boostLab_pion);
463 fProducts_pion[2].Boost(boostLab_pion);
464
465 }
466
467
468//-----------------------------------------------------------------------------//
469// Generate Rho resonance decay //
470//-----------------------------------------------------------------------------//
471
472 else if(idpart==113){
473 // calculate rho mass
474 if(wp_rho!=0.0){
475 mp_rho = pparent->M();
476 }
477 else{
478 Double_t x_rho=pparent->Px(); Double_t y_rho=pparent->Py(); Double_t z_rho=pparent->Pz();
479 Double_t t_rho=pparent->E();
480 Double_t p_rho=x_rho*x_rho+y_rho*y_rho+z_rho*z_rho;
481 Double_t Q2_rho=abs((t_rho*t_rho)-(p_rho*p_rho));
482 mp_rho = sqrt(Q2_rho);
483 }
484 // daughter
485 if ( mp_rho < 2.*md_rho )
486 {
487 printf("Rho into ee Decay kinematically impossible! \n");
488 return;
489 }
490
491 for( ;; ) {
492 // Sample the electron pair mass from a histogram
493 epmass_rho = fEPMassRho->GetRandom();
494 if ( mp_rho < 2.*epmass_rho ) break;
495 }
496
497 // electron pair kinematics in virtual photon rest frame
498 Ed_rho = epmass_rho/2.;
499 pd_rho = TMath::Sqrt((Ed_rho+md_rho)*(Ed_rho-md_rho));
500
501 // momentum vectors of electrons in virtual photon rest frame
502 Double_t pProd1_rho[3] = {pd_rho * sintheta * cosphi,
503 pd_rho * sintheta * sinphi,
504 pd_rho * costheta};
505
506 Double_t pProd2_rho[3] = {-1.0 * pd_rho * sintheta * cosphi,
507 -1.0 * pd_rho * sintheta * sinphi,
508 -1.0 * pd_rho * costheta};
509
510
511 // electron 4 vectors in properly rotated virtual photon rest frame
512 Double_t pRot1_rho[3] = {0.};
513 Rot(pProd1_rho, pRot1_rho, costheta, -sintheta, -cosphi, -sinphi);
514 Double_t pRot2_rho[3] = {0.};
515 Rot(pProd2_rho, pRot2_rho, costheta, -sintheta, -cosphi, -sinphi);
516 fProducts_rho[0].SetPx(pRot1_rho[0]);
517 fProducts_rho[0].SetPy(pRot1_rho[1]);
518 fProducts_rho[0].SetPz(pRot1_rho[2]);
519 fProducts_rho[0].SetE(Ed_rho);
520 fProducts_rho[1].SetPx(pRot2_rho[0]);
521 fProducts_rho[1].SetPy(pRot2_rho[1]);
522 fProducts_rho[1].SetPz(pRot2_rho[2]);
523 fProducts_rho[1].SetE(Ed_rho);
524
525
526 // boost decay products into the lab frame
527 TVector3 boostLab_rho(pparent->Px() / pparent->E(),
528 pparent->Py() / pparent->E(),
529 pparent->Pz() / pparent->E());
530
531 fProducts_rho[0].Boost(boostLab_rho);
532 fProducts_rho[1].Boost(boostLab_rho);
533 }
534
535
536//-----------------------------------------------------------------------------//
537// Generate Eta Dalitz decay //
538//-----------------------------------------------------------------------------//
539
540 else if(idpart==221){
541 pmass_eta = pparent->M();
542
543 for(;;){
544 // Sample the electron pair mass from a histogram
545 epmass_eta = fEPMassEta->GetRandom();
546 if(pmass_eta-omass_gamma>epmass_eta && epmass_eta/2.>emass)
547 break;
548 }
549
550 // electron pair kinematics in virtual photon rest frame
551 e1_eta = epmass_eta / 2.;
552 p1_eta = TMath::Sqrt((e1_eta + emass) * (e1_eta - emass));
553
554 // momentum vectors of electrons in virtual photon rest frame
555 Double_t pProd1_eta[3] = {p1_eta * sintheta * cosphi,
556 p1_eta * sintheta * sinphi,
557 p1_eta * costheta};
558 Double_t pProd2_eta[3] = {-1.0 * p1_eta * sintheta * cosphi,
559 -1.0 * p1_eta * sintheta * sinphi,
560 -1.0 * p1_eta * costheta};
561
562 // third child kinematics in parent meson rest frame
563 e3_gamma_eta = (pmass_eta * pmass_eta - epmass_eta * epmass_eta)/(2. * pmass_eta);
564 p3_gamma_eta = TMath::Sqrt((e3_gamma_eta * e3_gamma_eta));
565
566
567 // third child 4-vector in parent meson rest frame
568 fProducts_eta[2].SetPx(p3_gamma_eta * sintheta * cosphi);
569 fProducts_eta[2].SetPy(p3_gamma_eta * sintheta * sinphi);
570 fProducts_eta[2].SetPz(p3_gamma_eta * costheta);
571 fProducts_eta[2].SetE(e3_gamma_eta);
572
573 // electron 4-vectors in properly rotated virtual photon rest frame
574 Double_t pRot1_eta[3] = {0.};
575 Rot(pProd1_eta, pRot1_eta, costheta, -sintheta, -cosphi, -sinphi);
576 Double_t pRot2_eta[3] = {0.};
577 Rot(pProd2_eta, pRot2_eta, costheta, -sintheta, -cosphi, -sinphi);
578 fProducts_eta[0].SetPx(pRot1_eta[0]);
579 fProducts_eta[0].SetPy(pRot1_eta[1]);
580 fProducts_eta[0].SetPz(pRot1_eta[2]);
581 fProducts_eta[0].SetE(e1_eta);
582 fProducts_eta[1].SetPx(pRot2_eta[0]);
583 fProducts_eta[1].SetPy(pRot2_eta[1]);
584 fProducts_eta[1].SetPz(pRot2_eta[2]);
585 fProducts_eta[1].SetE(e1_eta);
586
587 // boost the dielectron into the parent meson's rest frame
588 Double_t eLPparent_eta = TMath::Sqrt(p3_gamma_eta * p3_gamma_eta + epmass_eta * epmass_eta);
589 TVector3 boostPair_eta( -1.0 * fProducts_eta[2].Px() / eLPparent_eta,
590 -1.0 * fProducts_eta[2].Py() / eLPparent_eta,
591 -1.0 * fProducts_eta[2].Pz() / eLPparent_eta);
592 fProducts_eta[0].Boost(boostPair_eta);
593 fProducts_eta[1].Boost(boostPair_eta);
594
595 // boost all decay products into the lab frame
596 TVector3 boostLab_eta(pparent->Px() / pparent->E(),
597 pparent->Py() / pparent->E(),
598 pparent->Pz() / pparent->E());
599
600 fProducts_eta[0].Boost(boostLab_eta);
601 fProducts_eta[1].Boost(boostLab_eta);
602 fProducts_eta[2].Boost(boostLab_eta);
603
604 }
605
606
607//-----------------------------------------------------------------------------//
608// Generate Omega Dalitz decay //
609//-----------------------------------------------------------------------------//
610
611 else if(idpart==223){
612 pmass_omega_dalitz = pparent->M();
613 for(;;){
614 // Sample the electron pair mass from a histogram
615 epmass_omega_dalitz = fEPMassOmegaDalitz->GetRandom();
616 if(pmass_omega_dalitz-omass_pion>epmass_omega_dalitz && epmass_omega_dalitz/2.>emass)
617 break;}
618
619 // electron pair kinematics in virtual photon rest frame
620 e1_omega = epmass_omega_dalitz / 2.;
621 p1_omega = TMath::Sqrt((e1_omega + emass) * (e1_omega - emass));
622
623 // momentum vectors of electrons in virtual photon rest frame
624 Double_t pProd1_omega_dalitz[3] = {p1_omega * sintheta * cosphi,
625 p1_omega * sintheta * sinphi,
626 p1_omega * costheta};
627 Double_t pProd2_omega_dalitz[3] = {-1.0 * p1_omega * sintheta * cosphi,
628 -1.0 * p1_omega * sintheta * sinphi,
629 -1.0 * p1_omega * costheta};
630
631 // third child kinematics in parent meson rest frame
632 e3_pion = (pmass_omega_dalitz * pmass_omega_dalitz + omass_pion * omass_pion - epmass_omega_dalitz * epmass_omega_dalitz)/(2. * pmass_omega_dalitz);
633 p3_pion = TMath::Sqrt((e3_pion + omass_pion) * (e3_pion - omass_pion));
634
635 // third child 4-vector in parent meson rest frame
636 fProducts_omega_dalitz[2].SetPx(p3_pion * sintheta * cosphi);
637 fProducts_omega_dalitz[2].SetPy(p3_pion * sintheta * sinphi);
638 fProducts_omega_dalitz[2].SetPz(p3_pion * costheta);
639 fProducts_omega_dalitz[2].SetE(e3_pion);
640
641 // lepton 4-vectors in properly rotated virtual photon rest frame
642 Double_t pRot1_omega_dalitz[3] = {0.};
643 Rot(pProd1_omega_dalitz, pRot1_omega_dalitz, costheta, -sintheta, -cosphi, -sinphi);
644 Double_t pRot2_omega_dalitz[3] = {0.};
645 Rot(pProd2_omega_dalitz, pRot2_omega_dalitz, costheta, -sintheta, -cosphi, -sinphi);
646 fProducts_omega_dalitz[0].SetPx(pRot1_omega_dalitz[0]);
647 fProducts_omega_dalitz[0].SetPy(pRot1_omega_dalitz[1]);
648 fProducts_omega_dalitz[0].SetPz(pRot1_omega_dalitz[2]);
649 fProducts_omega_dalitz[0].SetE(e1_omega);
650 fProducts_omega_dalitz[1].SetPx(pRot2_omega_dalitz[0]);
651 fProducts_omega_dalitz[1].SetPy(pRot2_omega_dalitz[1]);
652 fProducts_omega_dalitz[1].SetPz(pRot2_omega_dalitz[2]);
653 fProducts_omega_dalitz[1].SetE(e1_omega);
654
655 // boost the dielectron into the parent meson's rest frame
656 Double_t eLPparent_omega = TMath::Sqrt(p3_pion * p3_pion + epmass_omega_dalitz * epmass_omega_dalitz);
657 TVector3 boostPair_omega( -1.0 * fProducts_omega_dalitz[2].Px() / eLPparent_omega,
658 -1.0 * fProducts_omega_dalitz[2].Py() / eLPparent_omega,
659 -1.0 * fProducts_omega_dalitz[2].Pz() / eLPparent_omega);
660 fProducts_omega_dalitz[0].Boost(boostPair_omega);
661 fProducts_omega_dalitz[1].Boost(boostPair_omega);
662
663 // boost all decay products into the lab frame
664 TVector3 boostLab_omega_dalitz(pparent->Px() / pparent->E(),
665 pparent->Py() / pparent->E(),
666 pparent->Pz() / pparent->E());
667
668 fProducts_omega_dalitz[0].Boost(boostLab_omega_dalitz);
669 fProducts_omega_dalitz[1].Boost(boostLab_omega_dalitz);
670 fProducts_omega_dalitz[2].Boost(boostLab_omega_dalitz);
671
672
673//-----------------------------------------------------------------------------//
674// Generate Omega resonance decay //
675//-----------------------------------------------------------------------------//
676
677 if(wp_omega!=0.0){
678 // calculate omega mass
679 mp_omega = pparent->M();
680 }
681 else{
682 Double_t x_omega=pparent->Px(); Double_t y_omega=pparent->Py(); Double_t z_omega=pparent->Pz();
683 Double_t t_omega=pparent->E();
684 Double_t p_omega=x_omega*x_omega+y_omega*y_omega+z_omega*z_omega;
685 Double_t Q2_omega= abs((t_omega*t_omega)-(p_omega*p_omega));
686 mp_omega = sqrt(Q2_omega);
687 }
688
689 // daughter
690 if ( mp_omega< 2.*md_omega )
691 {
692 printf("Omega into ee Decay kinematically impossible! \n");
693 return;
694 }
695
696 for( ;; ) {
697 // Sample the electron pair mass from a histogram
698 epmass_omega = fEPMassOmega->GetRandom();
699 if( mp_omega < 2.*epmass_omega ) break;
700 }
701
702 // electron pair kinematics in virtual photon rest frame
703 Ed_omega = epmass_omega/2.;
704 pd_omega = TMath::Sqrt((Ed_omega+md_omega)*(Ed_omega-md_omega));
705
706 // momentum vectors of electrons in virtual photon rest frame
707 Double_t pProd1_omega[3] = {pd_omega * sintheta * cosphi,
708 pd_omega * sintheta * sinphi,
709 pd_omega * costheta};
710
711 Double_t pProd2_omega[3] = {-1.0 * pd_omega * sintheta * cosphi,
712 -1.0 * pd_omega * sintheta * sinphi,
713 -1.0 * pd_omega * costheta};
714
715
716 // lepton 4 vectors in properly rotated virtual photon rest frame
717 Double_t pRot1_omega[3] = {0.};
718 Rot(pProd1_omega, pRot1_omega, costheta, -sintheta, -cosphi, -sinphi);
719 Double_t pRot2_omega[3] = {0.};
720 Rot(pProd2_omega, pRot2_omega, costheta, -sintheta, -cosphi, -sinphi);
721 fProducts_omega[0].SetPx(pRot1_omega[0]);
722 fProducts_omega[0].SetPy(pRot1_omega[1]);
723 fProducts_omega[0].SetPz(pRot1_omega[2]);
724 fProducts_omega[0].SetE(Ed_omega);
725 fProducts_omega[1].SetPx(pRot2_omega[0]);
726 fProducts_omega[1].SetPy(pRot2_omega[1]);
727 fProducts_omega[1].SetPz(pRot2_omega[2]);
728 fProducts_omega[1].SetE(Ed_omega);
729
730 // boost decay products into the lab frame
731 TVector3 boostLab_omega(pparent->Px() / pparent->E(),
732 pparent->Py() / pparent->E(),
733 pparent->Pz() / pparent->E());
734
735 fProducts_omega[0].Boost(boostLab_omega);
736 fProducts_omega[1].Boost(boostLab_omega);
737
738 }
739
740//-----------------------------------------------------------------------------//
741// Generate Etaprime Dalitz decay //
742//-----------------------------------------------------------------------------//
743
744 else if(idpart==331){
745 pmass_etaprime = pparent->M();
746 for(;;){
747 // Sample the electron pair mass from a histogram
748 epmass_etaprime = fEPMassEtaPrime->GetRandom();
749 if(pmass_etaprime-omass_gamma>epmass_etaprime && epmass_etaprime/2.>emass)
750 break;}
751
752 // electron pair kinematics in virtual photon rest frame
753 e1_etaprime = epmass_etaprime / 2.;
754 p1_etaprime = TMath::Sqrt((e1_etaprime + emass) * (e1_etaprime - emass));
755
756 // momentum vectors of electrons in virtual photon rest frame
757 Double_t pProd1_etaprime[3] = {p1_etaprime * sintheta * cosphi,
758 p1_etaprime * sintheta * sinphi,
759 p1_etaprime * costheta};
760 Double_t pProd2_etaprime[3] = {-1.0 * p1_etaprime * sintheta * cosphi,
761 -1.0 * p1_etaprime * sintheta * sinphi,
762 -1.0 * p1_etaprime * costheta};
763
764 // third child kinematics in parent meson rest frame
765 e3_gamma_etaprime = (pmass_etaprime * pmass_etaprime + omass_gamma * omass_gamma - epmass_etaprime * epmass_etaprime)/(2. * pmass_etaprime);
766 p3_gamma_etaprime = TMath::Sqrt((e3_gamma_etaprime + omass_gamma) * (e3_gamma_etaprime - omass_gamma));
767
768 // third child 4-vector in parent meson rest frame
769 fProducts_etaprime[2].SetPx(p3_gamma_etaprime * sintheta * cosphi);
770 fProducts_etaprime[2].SetPy(p3_gamma_etaprime * sintheta * sinphi);
771 fProducts_etaprime[2].SetPz(p3_gamma_etaprime * costheta);
772 fProducts_etaprime[2].SetE(e3_gamma_etaprime);
773
774 // electron 4-vectors in properly rotated virtual photon rest frame
775 Double_t pRot1_etaprime[3] = {0.};
776 Rot(pProd1_etaprime, pRot1_etaprime, costheta, -sintheta, -cosphi, -sinphi);
777 Double_t pRot2_etaprime[3] = {0.};
778 Rot(pProd2_etaprime, pRot2_etaprime, costheta, -sintheta, -cosphi, -sinphi);
779 fProducts_etaprime[0].SetPx(pRot1_etaprime[0]);
780 fProducts_etaprime[0].SetPy(pRot1_etaprime[1]);
781 fProducts_etaprime[0].SetPz(pRot1_etaprime[2]);
782 fProducts_etaprime[0].SetE(e1_etaprime);
783 fProducts_etaprime[1].SetPx(pRot2_etaprime[0]);
784 fProducts_etaprime[1].SetPy(pRot2_etaprime[1]);
785 fProducts_etaprime[1].SetPz(pRot2_etaprime[2]);
786 fProducts_etaprime[1].SetE(e1_etaprime);
787
788 // boost the dielectron into the parent meson's rest frame
789 Double_t eLPparent_etaprime = TMath::Sqrt(p3_gamma_etaprime * p3_gamma_etaprime + epmass_etaprime * epmass_etaprime);
790 TVector3 boostPair_etaprime( -1.0 * fProducts_etaprime[2].Px() / eLPparent_etaprime,
791 -1.0 * fProducts_etaprime[2].Py() / eLPparent_etaprime,
792 -1.0 * fProducts_etaprime[2].Pz() / eLPparent_etaprime);
793 fProducts_etaprime[0].Boost(boostPair_etaprime);
794 fProducts_etaprime[1].Boost(boostPair_etaprime);
795
796 // boost all decay products into the lab frame
797 TVector3 boostLab_etaprime(pparent->Px() / pparent->E(),
798 pparent->Py() / pparent->E(),
799 pparent->Pz() / pparent->E());
800
801 fProducts_etaprime[0].Boost(boostLab_etaprime);
802 fProducts_etaprime[1].Boost(boostLab_etaprime);
803 fProducts_etaprime[2].Boost(boostLab_etaprime);
804
805 }
806
807//-----------------------------------------------------------------------------//
808// Generate Phi Dalitz decay //
809//-----------------------------------------------------------------------------//
810
811 else if(idpart==333){
812 pmass_phi_dalitz = pparent->M();
813 for(;;){
814 // Sample the electron pair mass from a histogram
815 epmass_phi_dalitz = fEPMassPhiDalitz->GetRandom();
816 if(pmass_phi_dalitz-omass_eta>epmass_phi_dalitz && epmass_phi_dalitz/2.>emass)
817 break;}
818
819 // electron pair kinematics in virtual photon rest frame
820 e1_phi = epmass_phi_dalitz / 2.;
821 p1_phi = TMath::Sqrt((e1_phi + emass) * (e1_phi - emass));
822
823 // momentum vectors of electrons in virtual photon rest frame
824 Double_t pProd1_phi_dalitz[3] = {p1_phi * sintheta * cosphi,
825 p1_phi * sintheta * sinphi,
826 p1_phi * costheta};
827 Double_t pProd2_phi_dalitz[3] = {-1.0 * p1_phi * sintheta * cosphi,
828 -1.0 * p1_phi * sintheta * sinphi,
829 -1.0 * p1_phi * costheta};
830
831 // third child kinematics in parent meson rest frame
832 e3_eta = (pmass_phi_dalitz * pmass_phi_dalitz + omass_eta * omass_eta - epmass_phi_dalitz * epmass_phi_dalitz)/(2. * pmass_phi_dalitz);
833 p3_eta = TMath::Sqrt((e3_eta + omass_eta) * (e3_eta - omass_eta));
834
835 // third child 4-vector in parent meson rest frame
836 fProducts_phi_dalitz[2].SetPx(p3_eta * sintheta * cosphi);
837 fProducts_phi_dalitz[2].SetPy(p3_eta * sintheta * sinphi);
838 fProducts_phi_dalitz[2].SetPz(p3_eta * costheta);
839 fProducts_phi_dalitz[2].SetE(e3_eta);
840
841 // electron 4-vectors in properly rotated virtual photon rest frame
842 Double_t pRot1_phi_dalitz[3] = {0.};
843 Rot(pProd1_phi_dalitz, pRot1_phi_dalitz, costheta, -sintheta, -cosphi, -sinphi);
844 Double_t pRot2_phi_dalitz[3] = {0.};
845 Rot(pProd2_phi_dalitz, pRot2_phi_dalitz, costheta, -sintheta, -cosphi, -sinphi);
846 fProducts_phi_dalitz[0].SetPx(pRot1_phi_dalitz[0]);
847 fProducts_phi_dalitz[0].SetPy(pRot1_phi_dalitz[1]);
848 fProducts_phi_dalitz[0].SetPz(pRot1_phi_dalitz[2]);
849 fProducts_phi_dalitz[0].SetE(e1_phi);
850 fProducts_phi_dalitz[1].SetPx(pRot2_phi_dalitz[0]);
851 fProducts_phi_dalitz[1].SetPy(pRot2_phi_dalitz[1]);
852 fProducts_phi_dalitz[1].SetPz(pRot2_phi_dalitz[2]);
853 fProducts_phi_dalitz[1].SetE(e1_phi);
854
855 // boost the dielectron into the parent meson's rest frame
856 Double_t eLPparent_phi = TMath::Sqrt(p3_eta * p3_eta + epmass_phi_dalitz * epmass_phi_dalitz);
857 TVector3 boostPair_phi( -1.0 * fProducts_phi_dalitz[2].Px() / eLPparent_phi,
858 -1.0 * fProducts_phi_dalitz[2].Py() / eLPparent_phi,
859 -1.0 * fProducts_phi_dalitz[2].Pz() / eLPparent_phi);
860 fProducts_phi_dalitz[0].Boost(boostPair_phi);
861 fProducts_phi_dalitz[1].Boost(boostPair_phi);
862
863 // boost all decay products into the lab frame
864 TVector3 boostLab_phi_dalitz(pparent->Px() / pparent->E(),
865 pparent->Py() / pparent->E(),
866 pparent->Pz() / pparent->E());
867
868 fProducts_phi_dalitz[0].Boost(boostLab_phi_dalitz);
869 fProducts_phi_dalitz[1].Boost(boostLab_phi_dalitz);
870 fProducts_phi_dalitz[2].Boost(boostLab_phi_dalitz);
871
872
873//-----------------------------------------------------------------------------//
874// Generate Phi resonance decay //
875//-----------------------------------------------------------------------------//
876
877 if(wp_phi!=0.0){
878 // calculate phi mass
879 mp_phi = pparent->M();
880 }
881 else{
882 Double_t x_phi=pparent->Px(); Double_t y_phi=pparent->Py(); Double_t z_phi=pparent->Pz();
883 Double_t t_phi=pparent->E();
884 Double_t p_phi=x_phi*x_phi+y_phi*y_phi+z_phi*z_phi;
885 Double_t Q2_phi= abs((t_phi*t_phi)-(p_phi*p_phi));
886 mp_phi = sqrt(Q2_phi);
887 }
888
889 if ( mp_phi< 2.*md_phi )
890 {
891 printf("Phi into ee Decay kinematically impossible! \n");
892 return;
893 }
894
895 for( ;; ) {
896 // Sample the electron pair mass from a histogram
897 epmass_phi = fEPMassPhi->GetRandom();
898 if(mp_phi < 2.*epmass_phi) break;
899 }
900
901 // electron pair kinematics in virtual photon rest frame
902 Ed_phi = epmass_phi/2.;
903 pd_phi = TMath::Sqrt((Ed_phi+md_phi)*(Ed_phi-md_phi));
904
905 // momentum vectors of electrons in virtual photon rest frame
906 Double_t pProd1_phi[3] = {pd_phi * sintheta * cosphi,
907 pd_phi * sintheta * sinphi,
908 pd_phi * costheta};
909 Double_t pProd2_phi[3] = {-1.0 * pd_phi * sintheta * cosphi,
910 -1.0 * pd_phi * sintheta * sinphi,
911 -1.0 * pd_phi * costheta};
912
913 // electron 4 vectors in properly rotated virtual photon rest frame
914 Double_t pRot1_phi[3] = {0.};
915 Rot(pProd1_phi, pRot1_phi, costheta, -sintheta, -cosphi, -sinphi);
916 Double_t pRot2_phi[3] = {0.};
917 Rot(pProd2_phi, pRot2_phi, costheta, -sintheta, -cosphi, -sinphi);
918 fProducts_phi[0].SetPx(pRot1_phi[0]);
919 fProducts_phi[0].SetPy(pRot1_phi[1]);
920 fProducts_phi[0].SetPz(pRot1_phi[2]);
921 fProducts_phi[0].SetE(Ed_phi);
922 fProducts_phi[1].SetPx(pRot2_phi[0]);
923 fProducts_phi[1].SetPy(pRot2_phi[1]);
924 fProducts_phi[1].SetPz(pRot2_phi[2]);
925 fProducts_phi[1].SetE(Ed_phi);
926
927 // boost decay products into the lab frame
928 TVector3 boostLab_phi(pparent->Px() / pparent->E(),
929 pparent->Py() / pparent->E(),
930 pparent->Pz() / pparent->E());
931
932 fProducts_phi[0].Boost(boostLab_phi);
933 fProducts_phi[1].Boost(boostLab_phi);
934
935 }
936
937//-----------------------------------------------------------------------------//
938// Generate Jpsi resonance decay //
939//-----------------------------------------------------------------------------//
940
941 else if(idpart==443){
942 // calculate jpsi mass
943 if(wp_jpsi!=0.0){
944 mp_jpsi = pparent->M();
945 }
946 else{
76e6e4c5 947 /*Double_t x_jpsi=pparent->Px();
948 Double_t y_jpsi=pparent->Py();
949 Double_t z_jpsi=pparent->Pz();
950 Double_t t_jpsi=pparent->E();
951 Double_t p_jpsi=x_jpsi*x_jpsi+y_jpsi*y_jpsi+z_jpsi*z_jpsi;
952 Double_t Q2_jpsi= abs((t_jpsi*t_jpsi)-(p_jpsi*p_jpsi));
953 mp_jpsi = sqrt(Q2_jpsi);*/
954
955 mp_jpsi = 3.096;
b5c4afc6 956
957 }
958
959 // daughter
960 if ( mp_jpsi < 2.*md_jpsi )
961 {
962 printf("JPsi into ee Decay kinematically impossible! \n");
963 return;
964 }
965
966 for( ;; ) {
967 // Sample the electron pair mass from a histogram
968 epmass_jpsi = fEPMassJPsi->GetRandom();
969 if ( mp_jpsi < 2.*epmass_jpsi ) break;
970 }
971 // electron pair kinematics in virtual photon rest frame
972 Ed_jpsi = epmass_jpsi/2.;
973 pd_jpsi = TMath::Sqrt((Ed_jpsi+md_jpsi)*(Ed_jpsi-md_jpsi));
974
975 // momentum vectors of electrons in virtual photon rest frame
976 Double_t pProd1_jpsi[3] = {pd_jpsi * sintheta * cosphi,
977 pd_jpsi * sintheta * sinphi,
978 pd_jpsi * costheta};
979
980 Double_t pProd2_jpsi[3] = {-1.0 * pd_jpsi * sintheta * cosphi,
981 -1.0 * pd_jpsi * sintheta * sinphi,
982 -1.0 * pd_jpsi * costheta};
983
984
985 // electron 4 vectors in properly rotated virtual photon rest frame
986 Double_t pRot1_jpsi[3] = {0.};
987 Rot(pProd1_jpsi, pRot1_jpsi, costheta, -sintheta, -cosphi, -sinphi);
988 Double_t pRot2_jpsi[3] = {0.};
989 Rot(pProd2_jpsi, pRot2_jpsi, costheta, -sintheta, -cosphi, -sinphi);
990 fProducts_jpsi[0].SetPx(pRot1_jpsi[0]);
991 fProducts_jpsi[0].SetPy(pRot1_jpsi[1]);
992 fProducts_jpsi[0].SetPz(pRot1_jpsi[2]);
993 fProducts_jpsi[0].SetE(Ed_jpsi);
994 fProducts_jpsi[1].SetPx(pRot2_jpsi[0]);
995 fProducts_jpsi[1].SetPy(pRot2_jpsi[1]);
996 fProducts_jpsi[1].SetPz(pRot2_jpsi[2]);
997 fProducts_jpsi[1].SetE(Ed_jpsi);
998
999
1000 // boost decay products into the lab frame
1001 TVector3 boostLab_jpsi(pparent->Px() / pparent->E(),
1002 pparent->Py() / pparent->E(),
1003 pparent->Pz() / pparent->E());
1004
1005 fProducts_jpsi[0].Boost(boostLab_jpsi);
1006 fProducts_jpsi[1].Boost(boostLab_jpsi);
1007
1008 }
1009
1010 return;
1011}
1012
1013void AliDecayerExodus::Rot(Double_t pin[3], Double_t pout[3], Double_t costheta, Double_t sintheta,
1014 Double_t cosphi, Double_t sinphi) const
1015{
1016// Perform rotation
1017 pout[0] = pin[0]*costheta*cosphi-pin[1]*sinphi+pin[2]*sintheta*cosphi;
1018 pout[1] = pin[0]*costheta*sinphi+pin[1]*cosphi+pin[2]*sintheta*sinphi;
1019 pout[2] = -1.0 * pin[0] * sintheta + pin[2] * costheta;
1020 return;
1021}
1022
1023
1024Int_t AliDecayerExodus::ImportParticles(TClonesArray *particles)
1025{
1026//
1027// Import particles for Dalitz and resonance decays
1028//
1029
1030 TClonesArray &clonesParticles = *particles;
1031
1032 Int_t i, k;
1033 Double_t px, py, pz, e;
1034
1035 Int_t pdgD [3][3] = { {kElectron, -kElectron, 22}, // pizero, eta, etaprime
1036 {kElectron, -kElectron, 111}, // omega dalitz
1037 {kElectron, -kElectron, 221} }; // phi dalitz
1038
1039 Int_t pdgR [2] = {kElectron, -kElectron}; // rho, omega, phi, jpsi
1040
1041
1042
1043 Int_t parentD[3] = { 0, 0, -1};
1044 Int_t dauD1 [3] = {-1, -1, 1};
1045 Int_t dauD2 [3] = {-1, -1, 2};
1046
1047 Int_t parentR[2] = { 0, 0};
1048 Int_t dauR1 [2] = { -1, -1};
1049 Int_t dauR2 [2] = { -1, -1};
1050
1051 for (Int_t j = 0; j < 9; j++){
1052
1053 // pizero
1054 if(j==0){
1055 for (i = 2; i > -1; i--) {
1056 px = fProducts_pion[i].Px();
1057 py = fProducts_pion[i].Py();
1058 pz = fProducts_pion[i].Pz();
1059 e = fProducts_pion[i].E();
1060 new(clonesParticles[2 - i]) TParticle(pdgD[0][i], 1, parentD[i], -1, dauD1[i], dauD2[i], px, py, pz, e, 0., 0., 0., 0.);
1061 }
1062 return (3);
1063 }
1064
1065 // rho
1066 else if(j==1){
1067 for (k = 1; k > -1; k--) {
1068 px = fProducts_rho[k].Px();
1069 py = fProducts_rho[k].Py();
1070 pz = fProducts_rho[k].Pz();
1071 e = fProducts_rho[k].E();
1072 new(clonesParticles[1 - k]) TParticle(pdgR[k], 1, parentR[k], -1, dauR1[k], dauR2[k], px, py, pz, e, 0., 0., 0., 0.);
1073 }
1074 return (2);
1075 }
1076
1077 // eta
1078 else if(j==2){
1079 for (i = 2; i > -1; i--) {
1080 px = fProducts_eta[i].Px();
1081 py = fProducts_eta[i].Py();
1082 pz = fProducts_eta[i].Pz();
1083 e = fProducts_eta[i].E();
1084 new(clonesParticles[2 - i]) TParticle(pdgD[0][i], 1, parentD[i], -1, dauD1[i], dauD2[i], px, py, pz, e, 0., 0., 0., 0.);
1085 }
1086 return (3);
1087 }
1088
1089 // omega dalitz
1090 else if(j==3){
1091 for (i = 2; i > -1; i--) {
1092 px = fProducts_omega_dalitz[i].Px();
1093 py = fProducts_omega_dalitz[i].Py();
1094 pz = fProducts_omega_dalitz[i].Pz();
1095 e = fProducts_omega_dalitz[i].E();
1096 new(clonesParticles[2 - i]) TParticle(pdgD[1][i], 1, parentD[i], -1, dauD1[i], dauD2[i], px, py, pz, e, 0., 0., 0., 0.);
1097 }
1098 return (3);
1099 }
1100
1101 // omega direct
1102 else if(j==4){
1103 for (k = 1; k > -1; k--) {
1104 px = fProducts_rho[k].Px();
1105 py = fProducts_rho[k].Py();
1106 pz = fProducts_rho[k].Pz();
1107 e = fProducts_rho[k].E();
1108 new(clonesParticles[1 - k]) TParticle(pdgR[k], 1, parentR[k], -1, dauR1[k], dauR2[k], px, py, pz, e, 0., 0., 0., 0.);
1109 }
1110 return (2);
1111 }
1112
1113 // etaprime
1114 else if(j==5){
1115 for (i = 2; i > -1; i--) {
1116 px = fProducts_etaprime[i].Px();
1117 py = fProducts_etaprime[i].Py();
1118 pz = fProducts_etaprime[i].Pz();
1119 e = fProducts_etaprime[i].E();
1120 new(clonesParticles[2 - i]) TParticle(pdgD[0][i], 1, parentD[i], -1, dauD1[i], dauD2[i], px, py, pz, e, 0., 0., 0., 0.);
1121 }
1122 return (3);
1123 }
1124
1125 // phi dalitz
1126 else if(j==6){
1127 for (i = 2; i > -1; i--) {
1128 px = fProducts_phi_dalitz[i].Px();
1129 py = fProducts_phi_dalitz[i].Py();
1130 pz = fProducts_phi_dalitz[i].Pz();
1131 e = fProducts_phi_dalitz[i].E();
1132 new(clonesParticles[2 - i]) TParticle(pdgD[2][i], 1, parentD[i], -1, dauD1[i], dauD2[i], px, py, pz, e, 0., 0., 0., 0.);
1133 }
1134 return (3);
1135 }
1136
1137
1138 // phi direct
1139 else if(j==7){
1140 for (k = 1; k > -1; k--) {
1141 px = fProducts_phi[k].Px();
1142 py = fProducts_phi[k].Py();
1143 pz = fProducts_phi[k].Pz();
1144 e = fProducts_phi[k].E();
1145 new(clonesParticles[1 - k]) TParticle(pdgR[k], 1, parentR[k], -1, dauR1[k], dauR2[k], px, py, pz, e, 0., 0., 0., 0.);
1146 }
1147 return (2);
1148 }
1149
1150 // jpsi direct
1151 else if(j==8){
1152 for (k = 1; k > -1; k--) {
1153 px = fProducts_jpsi[k].Px();
1154 py = fProducts_jpsi[k].Py();
1155 pz = fProducts_jpsi[k].Pz();
1156 e = fProducts_jpsi[k].E();
1157 new(clonesParticles[1 - k]) TParticle(pdgR[k], 1, parentR[k], -1, dauR1[k], dauR2[k], px, py, pz, e, 0., 0., 0., 0.);
1158 }
1159 return (2);
1160 }
1161
1162 }
1163
1164 return particles->GetEntries();
1165
1166}
1167
1168
1169void AliDecayerExodus::Decay(TClonesArray *array)
1170{
1171 // Replace all Dalitz(pi0,eta,omega,eta',phi) and resonance(rho,omega,phi,jpsi) decays with the correct matrix element decays
1172 // for di-electron cocktail calculations
1173
1174
1175 Int_t nt = array->GetEntriesFast();
1176 TParticle* dp3[3];
1177 TParticle* dp2[2];
1178 Int_t fd3, ld3, fd2, ld2, fd, ld;
1179 Int_t j, k;
1180
1181 for (Int_t i = 0; i < nt; i++) {
1182 TParticle* part = (TParticle*) (array->At(i));
1183 if (part->GetPdgCode() != 111 || part->GetPdgCode() != 221 || part->GetPdgCode() != 223 || part->GetPdgCode() != 331 || part->GetPdgCode() != 333 || part->GetPdgCode() != 443 ) continue;
1184
1185 //
1186 // Pizero Dalitz
1187 //
1188 if(part->GetPdgCode() == 111){
1189
1190 fd3 = part->GetFirstDaughter() - 1;
1191 ld3 = part->GetLastDaughter() - 1;
1192
1193 if (fd3 < 0) continue;
1194 if ((ld3 - fd3) != 2) continue;
1195
1196 for (j = 0; j < 3; j++) dp3[j] = (TParticle*) (array->At(fd3+j));
1197
1198 if((dp3[0]->GetPdgCode() != 22) && (TMath::Abs(dp3[1]->GetPdgCode()) != 11)) continue;
1199
1200 TLorentzVector Pizero(part->Px(), part->Py(), part->Pz(), part->Energy());
1201 Decay(111, &Pizero);
1202 for (j = 0; j < 3; j++) dp3[j]->SetMomentum(fProducts_pion[2-j]);
1203 }
1204
1205
1206 //
1207 // Eta Dalitz
1208 //
1209
1210 if(part->GetPdgCode() == 221){
1211
1212 fd3 = part->GetFirstDaughter() - 1;
1213 ld3 = part->GetLastDaughter() - 1;
1214
1215 if (fd3 < 0) continue;
1216 if ((ld3 - fd3) != 2) continue;
1217
1218 for (j = 0; j < 3; j++) dp3[j] = (TParticle*) (array->At(fd3+j));
1219
1220 if((dp3[0]->GetPdgCode() != 22) && ((TMath::Abs(dp3[1]->GetPdgCode()) != 11))) continue;
1221
1222 TLorentzVector Eta(part->Px(), part->Py(), part->Pz(), part->Energy());
1223 Decay(221, &Eta);
1224 for (j = 0; j < 3; j++) dp3[j]->SetMomentum(fProducts_eta[2-j]);
1225 }
1226
1227 //
1228 // Rho
1229 //
1230
1231 if(part->GetPdgCode() == 113){
1232
1233 fd2 = part->GetFirstDaughter() - 1;
1234 ld2 = part->GetLastDaughter() - 1;
1235
1236 if (fd2 < 0) continue;
1237 if ((ld2 - fd2) != 1) continue;
1238
1239 for (k = 0; k < 2; k++) dp2[k] = (TParticle*) (array->At(fd2+k));
1240
1241 if((dp2[0]->GetPdgCode() != 11) && ((TMath::Abs(dp2[1]->GetPdgCode()) != 11))) continue;
1242
1243 TLorentzVector Rho(part->Px(), part->Py(), part->Pz(), part->Energy());
1244 Decay(113, &Rho);
1245 for (k = 0; k < 2; k++) dp2[k]->SetMomentum(fProducts_rho[1-k]);
1246 }
1247
1248 //
1249 // Omega dalitz and direct
1250 //
1251
1252 if(part->GetPdgCode() == 223){
1253
1254 fd = part->GetFirstDaughter() - 1;
1255 ld = part->GetLastDaughter() - 1;
1256
1257 if (fd < 0) continue;
1258
1259 if ((ld - fd) == 2){
1260
1261 for (j = 0; j < 3; j++) dp3[j] = (TParticle*) (array->At(fd+j));
1262 if( dp3[0]->GetPdgCode() != 111 && (TMath::Abs(dp3[1]->GetPdgCode()) != 11)) continue;
1263
1264 TLorentzVector Omegadalitz(part->Px(), part->Py(), part->Pz(), part->Energy());
1265 Decay(223, &Omegadalitz);
1266 for (j = 0; j < 3; j++) dp3[j]->SetMomentum(fProducts_omega_dalitz[2-j]);
1267 }
1268
1269 else if ((ld - fd) == 1) {
1270
1271 for (k = 0; k < 2; k++) dp2[k] = (TParticle*) (array->At(fd+k));
1272 if( dp2[0]->GetPdgCode() != 11 && (TMath::Abs(dp2[1]->GetPdgCode()) != 11)) continue;
1273
1274 TLorentzVector Omega(part->Px(), part->Py(), part->Pz(), part->Energy());
1275 Decay(223, &Omega);
1276 for (k = 0; k < 2; k++) dp2[k]->SetMomentum(fProducts_omega[1-k]);
1277 }
1278 }
1279
1280 //
1281 // Etaprime dalitz
1282 //
1283
1284 if(part->GetPdgCode() == 331){
1285
1286 fd3 = part->GetFirstDaughter() - 1;
1287 ld3 = part->GetLastDaughter() - 1;
1288
1289 if (fd3 < 0) continue;
1290 if ((ld3 - fd3) != 2) continue;
1291
1292 for (j = 0; j < 3; j++) dp3[j] = (TParticle*) (array->At(fd3+j));
1293
1294 if((dp3[0]->GetPdgCode() != 22) && ((TMath::Abs(dp3[1]->GetPdgCode()) != 11))) continue;
1295
1296 TLorentzVector Etaprime(part->Px(), part->Py(), part->Pz(), part->Energy());
1297 Decay(331, &Etaprime);
1298 for (j = 0; j < 3; j++) dp3[j]->SetMomentum(fProducts_etaprime[2-j]);
1299 }
1300
1301 //
1302 // Phi dalitz and direct
1303 //
1304 if(part->GetPdgCode() == 333){
1305
1306 fd = part->GetFirstDaughter() - 1;
1307 ld = part->GetLastDaughter() - 1;
1308
1309 if (fd < 0) continue;
1310 if ((ld - fd) == 2){
1311 for (j = 0; j < 3; j++) dp3[j] = (TParticle*) (array->At(fd+j));
1312 if( dp3[0]->GetPdgCode() != 221 && (TMath::Abs(dp3[1]->GetPdgCode()) != 11)) continue;
1313
1314 TLorentzVector Phidalitz(part->Px(), part->Py(), part->Pz(), part->Energy());
1315 Decay(333, &Phidalitz);
1316 for (j = 0; j < 3; j++) dp3[j]->SetMomentum(fProducts_phi_dalitz[2-j]);
1317 }
1318
1319 else if ((ld - fd) == 1) {
1320 for (k = 0; k < 2; k++) dp2[k] = (TParticle*) (array->At(fd+k));
1321 if( dp2[0]->GetPdgCode() != 11 && (TMath::Abs(dp2[1]->GetPdgCode()) != 11)) continue;
1322
1323 TLorentzVector Phi(part->Px(), part->Py(), part->Pz(), part->Energy());
1324 Decay(333, &Phi);
1325 for (k = 0; k < 2; k++) dp2[k]->SetMomentum(fProducts_phi[1-k]);
1326 }
1327 }
1328
1329 //
1330 // JPsi
1331 //
1332
1333 if(part->GetPdgCode() == 443){
1334
1335 fd2 = part->GetFirstDaughter() - 1;
1336 ld2 = part->GetLastDaughter() - 1;
1337
1338 if (fd2 < 0) continue;
1339 if ((ld2 - fd2) != 1) continue;
1340
1341 for (k = 0; k < 2; k++) dp2[k] = (TParticle*) (array->At(fd2+k));
1342
1343 if((dp2[0]->GetPdgCode() != 11) && ((TMath::Abs(dp2[1]->GetPdgCode()) != 11))) continue;
1344
1345 TLorentzVector JPsi(part->Px(), part->Py(), part->Pz(), part->Energy());
1346 Decay(443, &JPsi);
1347 for (k = 0; k < 2; k++) dp2[k]->SetMomentum(fProducts_jpsi[1-k]);
1348 }
1349
1350 }
1351}
1352
1353
1354AliDecayerExodus& AliDecayerExodus::operator=(const AliDecayerExodus& rhs)
1355{
1356 // Assignment operator
1357 rhs.Copy(*this);
1358 return *this;
1359}
1360
1361void AliDecayerExodus::Copy(TObject&) const
1362{
1363 //
1364 // Copy
1365 //
1366 Fatal("Copy","Not implemented!\n");
1367}
1368
1369
1370AliDecayerExodus::AliDecayerExodus(const AliDecayerExodus &decayer)
1371 : AliDecayer(),
1372 fEPMassPion(0),
1373 fEPMassEta(0),
1374 fEPMassEtaPrime(0),
1375 fEPMassRho(0),
1376 fEPMassOmega(0),
1377 fEPMassOmegaDalitz(0),
1378 fEPMassPhi(0),
1379 fEPMassPhiDalitz(0),
1380 fEPMassJPsi(0),
1381 fInit(0)
1382{
1383 // Copy Constructor
1384 decayer.Copy(*this);
1385}
1386
1387
1388