]>
Commit | Line | Data |
---|---|---|
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 | ||
30 | ClassImp(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 | 48 | AliDecayerExodus::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 | ||
66 | void 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 | ||
307 | Double_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 | 336 | Double_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 | 350 | void 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 | ||
1013 | void 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 | ||
1024 | Int_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 | ||
1169 | void 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 | ||
1354 | AliDecayerExodus& AliDecayerExodus::operator=(const AliDecayerExodus& rhs) | |
1355 | { | |
1356 | // Assignment operator | |
1357 | rhs.Copy(*this); | |
1358 | return *this; | |
1359 | } | |
1360 | ||
1361 | void AliDecayerExodus::Copy(TObject&) const | |
1362 | { | |
1363 | // | |
1364 | // Copy | |
1365 | // | |
1366 | Fatal("Copy","Not implemented!\n"); | |
1367 | } | |
1368 | ||
1369 | ||
1370 | AliDecayerExodus::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 |