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