Line data Source code
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 : #include <TF1.h>
29 :
30 :
31 6 : ClassImp(AliDecayerExodus)
32 :
33 : //---------------------------------------------------------------------------------------------------
34 : //
35 : // Generate electron-pair mass distributions for Dalitz decays according
36 : // to the Kroll-Wada parametrization: N. Kroll, W. Wada: Phys. Rev 98(1955)1355
37 : // and generate electron-pair mass distributions for resonances according
38 : // to the Gounaris-Sakurai parametrization: G.J. Gounaris, J.J. Sakurai: Phys.Rev.Lett. 21(1968)244
39 : //
40 : // For the electromagnetic form factor the parameterization from
41 : // Lepton-G is used: L.G. Landsberg et al.: Phys. Rep. 128(1985)301
42 : //
43 : // Ralf Averbeck (R.Averbeck@gsi.de)
44 : // Irem Erdemir (irem.erdemir@cern.ch)
45 : //
46 : //---------------------------------------------------------------------------------------------------
47 :
48 :
49 55 : AliDecayerExodus::AliDecayerExodus():
50 1 : AliDecayer(),
51 1 : fEPMassPion(0),
52 1 : fEPMassEta(0),
53 1 : fEPMassEtaPrime(0),
54 1 : fEPMassRho(0),
55 1 : fEPMassOmega(0),
56 1 : fEPMassOmegaDalitz(0),
57 1 : fEPMassPhi(0),
58 1 : fEPMassPhiDalitz(0),
59 1 : fEPMassJPsi(0),
60 3 : fPol(new TF1("dsigdcostheta","1.+[0]*x*x",-1.,1.)), /* Polarization Function for resonances */
61 1 : fInit(0)
62 :
63 5 : {
64 : // Constructor
65 2 : }
66 :
67 :
68 : void AliDecayerExodus::Init()
69 : {
70 :
71 : // Initialisation
72 : //
73 : Int_t ibin, nbins;
74 : Double_t min, maxpion, maxeta, maxomega, maxetaprime, maxphi, binwidth_pion, binwidth_eta, binwidth_omega, binwidth_etaprime, binwidth_phi;
75 : Double_t pionmass, etamass, omegamass, etaprimemass, phimass, emass, proton_mass, omasspion, omasseta, omassgamma;
76 : Double_t epsilon_pion, epsilon_eta, epsilon_omega, epsilon_etaprime, epsilon_phi;
77 : Double_t delta_pion, delta_eta, delta_omega, delta_etaprime, delta_phi;
78 : Double_t mLL_pion, mLL_eta, mLL_omega, mLL_etaprime, mLL_phi;
79 : Double_t q_pion, q_eta, q_omega, q_etaprime, q_phi;
80 : Double_t kwHelp_pion, kwHelp_eta, kwHelp_omega, kwHelp_etaprime, kwHelp_phi;
81 : Double_t krollWada_pion, krollWada_eta, krollWada_omega, krollWada_etaprime, krollWada_phi;
82 : Double_t formFactor_pion, formFactor_eta, formFactor_omega, formFactor_etaprime, formFactor_phi;
83 : Double_t weight_pion, weight_eta, weight_omega_dalitz, weight_etaprime, weight_phi_dalitz;
84 :
85 : Float_t binwidth;
86 : Float_t mass_bin, mass_min, mass_max;
87 : Double_t vmass_rho, vmass_omega, vmass_phi, vmass_jpsi, vwidth_rho, vwidth_omega, vwidth_phi, vwidth_jpsi;
88 : Double_t weight_rho, weight_omega, weight_phi, weight_jpsi;
89 :
90 : //================================================================================//
91 : // Create electron pair mass histograms from dalitz decays //
92 : //================================================================================//
93 :
94 : // Get the particle masses
95 : // parent
96 : nbins = 1000;
97 : mass_min = 0.;
98 : mass_max = 0.;
99 0 : pionmass = (TDatabasePDG::Instance()->GetParticle(111))->Mass();
100 0 : etamass = (TDatabasePDG::Instance()->GetParticle(221))->Mass();
101 0 : omegamass = (TDatabasePDG::Instance()->GetParticle(223))->Mass();
102 0 : etaprimemass = (TDatabasePDG::Instance()->GetParticle(331))->Mass();
103 0 : phimass = (TDatabasePDG::Instance()->GetParticle(333))->Mass();
104 : // child - electron
105 0 : emass = (TDatabasePDG::Instance()->GetParticle(11))->Mass();
106 : // child - other : third childs from Dalitz decays
107 : omasspion = pionmass;
108 : omasseta = etamass;
109 : omassgamma = 0.;
110 :
111 0 : min = 2.0 * emass;
112 : maxpion = pionmass - omassgamma;
113 : maxeta = etamass - omassgamma;
114 0 : maxomega = omegamass - pionmass;
115 : maxetaprime = etaprimemass - omassgamma;
116 0 : maxphi = phimass - omasseta;
117 :
118 0 : binwidth_pion = (maxpion - min) / (Double_t)nbins;
119 0 : binwidth_eta = (maxeta - min) / (Double_t)nbins;
120 0 : binwidth_omega = (maxomega - min) / (Double_t)nbins;
121 0 : binwidth_etaprime = (maxetaprime - min) / (Double_t)nbins;
122 0 : binwidth_phi = (maxphi - min) / (Double_t)nbins;
123 :
124 :
125 0 : epsilon_pion = (emass / pionmass) * (emass / pionmass);
126 0 : epsilon_eta = (emass / etamass) * (emass / etamass);
127 0 : epsilon_omega = (emass / omegamass) * (emass / omegamass);
128 0 : epsilon_etaprime = (emass / etaprimemass) * (emass / etaprimemass);
129 0 : epsilon_phi = (emass / phimass) * (emass / phimass);
130 :
131 :
132 0 : delta_pion = (omassgamma / pionmass) * (omassgamma / pionmass);
133 0 : delta_eta = (omassgamma / etamass) * (omassgamma / etamass);
134 0 : delta_omega = (omasspion / omegamass) * (omasspion / omegamass);
135 0 : delta_etaprime = (omassgamma / etaprimemass) * (omassgamma / etaprimemass);
136 0 : delta_phi = (omasseta / phimass) * (omasseta / phimass);
137 :
138 : // create pair mass histograms for Dalitz decays of pi0, eta, omega, eta' and phi
139 0 : if (!fEPMassPion) {delete fEPMassPion; fEPMassPion = new TH1F("fEPMassPion", "Dalitz electron pair from pion", nbins, min, maxpion); }
140 0 : if (!fEPMassEta) {delete fEPMassEta; fEPMassEta = new TH1F("fEPMassEta", "Dalitz electron pair from eta", nbins, min, maxeta);}
141 0 : if (!fEPMassOmegaDalitz) {delete fEPMassOmegaDalitz; fEPMassOmegaDalitz = new TH1F("fEPMassOmegaDalitz", "Dalitz electron pair from omega ", nbins, min, maxomega);}
142 0 : if (!fEPMassEtaPrime) {delete fEPMassEtaPrime; fEPMassEtaPrime = new TH1F("fEPMassEtaPrime", "Dalitz electron pair from etaprime", nbins, min, maxetaprime);}
143 0 : 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 0 : for (ibin = 1; ibin <= nbins; ibin++ )
149 : {
150 0 : mLL_pion = min + (Double_t)(ibin - 1) * binwidth_pion + binwidth_pion / 2.0;
151 0 : mLL_eta = min + (Double_t)(ibin - 1) * binwidth_eta + binwidth_eta / 2.0;
152 0 : mLL_omega = min + (Double_t)(ibin - 1) * binwidth_omega + binwidth_omega / 2.0;
153 0 : mLL_etaprime = min + (Double_t)(ibin - 1) * binwidth_etaprime + binwidth_etaprime / 2.0;
154 0 : mLL_phi = min + (Double_t)(ibin - 1) * binwidth_phi + binwidth_phi / 2.0;
155 :
156 :
157 0 : q_pion = (mLL_pion / pionmass) * (mLL_pion / pionmass);
158 0 : q_eta = (mLL_eta / etamass) * (mLL_eta / etamass);
159 0 : q_omega = (mLL_omega / omegamass)*(mLL_omega / omegamass);
160 0 : q_etaprime = (mLL_etaprime / etaprimemass) * (mLL_etaprime / etaprimemass);
161 0 : q_phi = (mLL_phi / phimass) * (mLL_phi / phimass);
162 :
163 0 : 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 0 : AliFatal("Error in calculating Dalitz mass histogram binning!");
166 0 : }
167 :
168 :
169 0 : kwHelp_pion = (1.0 + q_pion / (1.0 - delta_pion)) * (1.0 + q_pion / (1.0 - delta_pion))
170 0 : - 4.0 * q_pion / ((1.0 - delta_pion) * (1.0 - delta_pion));
171 :
172 0 : kwHelp_eta = (1.0 + q_eta / (1.0 - delta_eta)) * (1.0 + q_eta / (1.0 - delta_eta))
173 0 : - 4.0 * q_eta / ((1.0 - delta_eta) * (1.0 - delta_eta));
174 :
175 0 : kwHelp_omega = (1.0 + q_omega / (1.0 - delta_omega)) * (1.0 + q_omega / (1.0 - delta_omega))
176 0 : - 4.0 * q_omega / ((1.0 - delta_omega) * (1.0 - delta_omega));
177 :
178 0 : kwHelp_etaprime = (1.0 + q_etaprime / (1.0 - delta_etaprime)) * (1.0 + q_etaprime / (1.0 - delta_etaprime))
179 0 : - 4.0 * q_etaprime / ((1.0 - delta_etaprime) * (1.0 - delta_etaprime));
180 :
181 0 : kwHelp_phi = (1.0 + q_phi / (1.0 - delta_phi)) * (1.0 + q_phi / (1.0 - delta_phi))
182 0 : - 4.0 * q_phi / ((1.0 - delta_phi) * (1.0 - delta_phi));
183 :
184 :
185 :
186 :
187 0 : if ( kwHelp_pion <= 0.0 || kwHelp_eta <= 0.0 || kwHelp_omega <= 0.0 || kwHelp_etaprime <= 0.0 || kwHelp_phi <= 0.0 )
188 : {
189 0 : AliFatal("Error in calculating Dalitz mass histogram binning!");
190 :
191 0 : }
192 :
193 :
194 : // Invariant mass distributions of electron pairs from Dalitz decays
195 : // using Kroll-Wada function
196 :
197 0 : krollWada_pion = (2.0 / mLL_pion) * TMath::Exp(1.5 * TMath::Log(kwHelp_pion))
198 0 : * TMath::Sqrt(1.0 - 4.0 * epsilon_pion / q_pion)
199 0 : * (1.0 + 2.0 * epsilon_pion / q_pion);
200 :
201 :
202 0 : krollWada_eta = (2.0 / mLL_eta) * TMath::Exp(1.5 * TMath::Log(kwHelp_eta))
203 0 : * TMath::Sqrt(1.0 - 4.0 * epsilon_eta / q_eta)
204 0 : * (1.0 + 2.0 * epsilon_eta / q_eta);
205 :
206 :
207 0 : krollWada_omega = (2.0 / mLL_omega) * TMath::Exp(1.5 * TMath::Log(kwHelp_omega))
208 0 : * TMath::Sqrt(1.0 - 4.0 * epsilon_omega / q_omega)
209 0 : * (1.0 + 2.0 * epsilon_omega / q_omega);
210 :
211 :
212 0 : krollWada_etaprime = (2.0 / mLL_etaprime) * TMath::Exp(1.5 * TMath::Log(kwHelp_etaprime))
213 0 : * TMath::Sqrt(1.0 - 4.0 * epsilon_etaprime / q_etaprime)
214 0 : * (1.0 + 2.0 * epsilon_etaprime / q_etaprime);
215 :
216 0 : krollWada_phi = (2.0 / mLL_phi) * TMath::Exp(1.5 * TMath::Log(kwHelp_phi))
217 0 : * TMath::Sqrt(1.0 - 4.0 * epsilon_phi / q_phi)
218 0 : * (1.0 + 2.0 * epsilon_phi / q_phi);
219 :
220 :
221 :
222 : // Form factors from Lepton-G
223 0 : formFactor_pion = 1.0/(1.0-5.5*mLL_pion*mLL_pion);
224 0 : formFactor_eta = 1.0/(1.0-1.9*mLL_eta*mLL_eta);
225 0 : formFactor_omega = (TMath::Power(TMath::Power(0.6519,2),2))
226 0 : /(TMath::Power(TMath::Power(0.6519,2)-TMath::Power(mLL_omega, 2), 2)
227 0 : + TMath::Power(0.04198, 2)*TMath::Power(0.6519, 2));
228 0 : formFactor_etaprime = (TMath::Power(TMath::Power(0.764,2),2))
229 0 : /(TMath::Power(TMath::Power(0.764,2)-TMath::Power(mLL_etaprime, 2), 2)
230 0 : + TMath::Power(0.1020, 2)*TMath::Power(0.764, 2));
231 : formFactor_phi = 1.0;
232 :
233 :
234 :
235 :
236 0 : weight_pion = krollWada_pion * formFactor_pion * formFactor_pion;
237 0 : weight_eta = krollWada_eta * formFactor_eta * formFactor_eta;
238 0 : weight_omega_dalitz = krollWada_omega * formFactor_omega;
239 0 : 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 0 : fEPMassPion ->AddBinContent(ibin, weight_pion);
245 0 : fEPMassEta ->AddBinContent(ibin, weight_eta);
246 0 : fEPMassOmegaDalitz->AddBinContent(ibin, weight_omega_dalitz);
247 0 : fEPMassEtaPrime ->AddBinContent(ibin, weight_etaprime);
248 0 : 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 0 : vmass_rho = (TDatabasePDG::Instance()->GetParticle(113))->Mass();
263 0 : vmass_omega = (TDatabasePDG::Instance()->GetParticle(223))->Mass();
264 0 : vmass_phi = (TDatabasePDG::Instance()->GetParticle(333))->Mass();
265 0 : vmass_jpsi = (TDatabasePDG::Instance()->GetParticle(443))->Mass();
266 : // Get the particle widths
267 : // parent
268 0 : vwidth_rho = (TDatabasePDG::Instance()->GetParticle(113))->Width();
269 0 : vwidth_omega = (TDatabasePDG::Instance()->GetParticle(223))->Width();
270 0 : vwidth_phi = (TDatabasePDG::Instance()->GetParticle(333))->Width();
271 0 : vwidth_jpsi = (TDatabasePDG::Instance()->GetParticle(443))->Width();
272 :
273 :
274 0 : if ( mass_min == 0. && mass_max == 0. )
275 : {
276 : mass_min = 2.*pimass;
277 : mass_max = 5.;
278 0 : }
279 :
280 0 : binwidth = (mass_max-mass_min)/(Double_t)nbins;
281 :
282 : // create pair mass histograms for resonances of rho, omega, phi and jpsi
283 0 : if (!fEPMassRho) {delete fEPMassRho; fEPMassRho = new TH1F("fEPMassRho","mass rho",nbins,mass_min,mass_max);}
284 0 : if (!fEPMassOmega) {delete fEPMassOmega; fEPMassOmega = new TH1F("fEPMassOmega","mass omega",nbins,mass_min,mass_max);}
285 0 : if (!fEPMassPhi) {delete fEPMassPhi; fEPMassPhi = new TH1F("fEPMassPhi","mass phi",nbins,mass_min,mass_max);}
286 0 : if (!fEPMassJPsi) {delete fEPMassJPsi; fEPMassJPsi = new TH1F("fEPMassJPsi","mass jpsi",nbins,mass_min,mass_max);}
287 :
288 0 : for (ibin=1; ibin<=nbins; ibin++ )
289 : {
290 0 : mass_bin = mass_min+(Double_t)(ibin-1)*binwidth+binwidth/2.0;
291 :
292 0 : weight_rho = (Float_t)GounarisSakurai(mass_bin,vmass_rho,vwidth_rho,emass);
293 0 : weight_omega = (Float_t)GounarisSakurai(mass_bin,vmass_omega,vwidth_omega,emass);
294 0 : weight_phi = (Float_t)GounarisSakurai(mass_bin,vmass_phi,vwidth_phi,emass);
295 0 : weight_jpsi = (Float_t)Lorentz(mass_bin,vmass_jpsi,vwidth_jpsi);
296 :
297 : // Fill histograms of electron pair masses from resonance decays
298 0 : fEPMassRho ->AddBinContent(ibin,weight_rho);
299 0 : fEPMassOmega->AddBinContent(ibin,weight_omega);
300 0 : fEPMassPhi ->AddBinContent(ibin,weight_phi);
301 0 : fEPMassJPsi ->AddBinContent(ibin,weight_jpsi);
302 : }
303 :
304 0 : }
305 :
306 : Double_t AliDecayerExodus::GounarisSakurai(Float_t mass, Double_t vmass, Double_t vwidth, Double_t emass)
307 : {
308 : // Invariant mass distributions of electron pairs from resonance decays
309 : // of rho, omega and phi
310 : // using Gounaris-Sakurai function
311 :
312 : Double_t corr = 0.;
313 : Double_t epsilon = 0.;
314 : Double_t weight = 0.;
315 :
316 : Double_t pimass = 0.13956995;
317 :
318 0 : if(mass>pimass){
319 0 : corr = vwidth*(vmass/mass)*exp(1.5*log((mass*mass/4.0-pimass*pimass)
320 0 : /(vmass*vmass/4.0-pimass*pimass)));
321 0 : }
322 :
323 0 : epsilon = (emass/mass)*(emass/mass);
324 :
325 0 : if ( 1.0-4.0*epsilon>=0.0 )
326 : {
327 0 : weight = sqrt(1.0-4.0*epsilon)*(1.0+2.0*epsilon)/
328 0 : ((vmass*vmass-mass*mass)*(vmass*vmass-mass*mass)+
329 0 : (vmass*corr)*(vmass*corr));
330 0 : }
331 0 : return weight;
332 : }
333 :
334 :
335 : Double_t AliDecayerExodus::Lorentz(Float_t mass, Double_t vmass, Double_t vwidth)
336 : {
337 : // Invariant mass distributions of electron pairs from resonance decay
338 : // of jpsi (and it can also be used for other particles except rho, omega and phi)
339 : // using Lorentz function
340 :
341 : Double_t weight;
342 :
343 0 : weight = (vwidth*vwidth/4.0)/(vwidth*vwidth/4.0+(vmass-mass)*(vmass-mass));
344 :
345 0 : return weight;
346 :
347 : }
348 :
349 : void AliDecayerExodus::Decay(Int_t idpart, TLorentzVector* pparent)
350 : {
351 :
352 0 : if (!fInit) {
353 0 : Init();
354 0 : fInit=1;
355 0 : }
356 :
357 : //local variables for dalitz/2-body decay:
358 : Double_t pmass, epmass, realp_mass, e1, p1, e3, p3;
359 : Double_t wp_res, mp_res, md_res, epmass_res, Ed_res, pd_res;
360 : Double_t PolPar;
361 0 : TLorentzVector fProducts_res[2], fProducts_dalitz[3];
362 : Int_t idRho=113;
363 : Int_t idOmega=223;
364 : Int_t idPhi=333;
365 : Int_t idJPsi=443;
366 : Int_t idPi0=111;
367 : Int_t idEta=221;
368 : Int_t idEtaPrime=331;
369 :
370 : // Get the particle masses of daughters
371 : Double_t emass, proton_mass, omass_pion, omass_eta, omass_gamma;
372 0 : emass = (TDatabasePDG::Instance()->GetParticle(11)) ->Mass();
373 0 : proton_mass = (TDatabasePDG::Instance()->GetParticle(2212)) ->Mass();
374 0 : omass_pion = (TDatabasePDG::Instance()->GetParticle(111))->Mass();
375 0 : omass_eta = (TDatabasePDG::Instance()->GetParticle(221))->Mass();
376 0 : omass_gamma = (TDatabasePDG::Instance()->GetParticle(22)) ->Mass();
377 :
378 : //flat angular distributions
379 : Double_t costheta, sintheta, cosphi, sinphi, phi;
380 : Double_t beta_square, lambda;
381 0 : costheta = (2.0 * gRandom->Rndm()) - 1.;
382 0 : sintheta = TMath::Sqrt((1. + costheta) * (1. - costheta));
383 0 : phi = 2.0 * TMath::ACos(-1.) * gRandom->Rndm();
384 0 : sinphi = TMath::Sin(phi);
385 0 : cosphi = TMath::Cos(phi);
386 :
387 :
388 : //-----------------------------------------------------------------------------//
389 : // Generate Dalitz decays: Pi0/Eta/Omega/EtaPrime/Phi //
390 : //-----------------------------------------------------------------------------//
391 :
392 0 : if(idpart==idPi0||idpart==idEta||idpart==idOmega||idpart==idEtaPrime||idpart==idPhi){
393 :
394 : //get the parent mass
395 0 : pmass = pparent->M();
396 :
397 : // Sample the electron pair mass from a histogram
398 0 : for(;;){
399 0 : if(idpart==idPi0){
400 0 : epmass = fEPMassPion->GetRandom();
401 : realp_mass=omass_gamma;
402 0 : }else if(idpart==idEta){
403 0 : epmass = fEPMassEta->GetRandom();
404 : realp_mass=omass_gamma;
405 0 : }else if(idpart==idOmega){
406 0 : epmass = fEPMassOmegaDalitz->GetRandom();
407 : realp_mass=omass_pion;
408 0 : }else if(idpart==idEtaPrime){
409 0 : epmass = fEPMassEtaPrime->GetRandom();
410 : realp_mass=omass_gamma;
411 0 : }else if(idpart==idPhi){
412 0 : epmass = fEPMassPhiDalitz->GetRandom();
413 : realp_mass=omass_eta;
414 0 : }else{ printf(" Exodus ERROR: Dalitz mass parametrization not found \n");
415 0 : return;
416 : }
417 0 : if(pmass-realp_mass>epmass && epmass/2.>emass) break;
418 : }
419 :
420 : // electron pair kinematics in virtual photon rest frame
421 : e1 = epmass / 2.;
422 0 : p1 = TMath::Sqrt((e1 + emass) * (e1 - emass));
423 :
424 :
425 : //Polarization parameters (lambda) for Dalitz:
426 0 : if ( realp_mass<0.01 ){
427 0 : beta_square = 1.0 - 4.0*(emass*emass)/(epmass*epmass);
428 0 : lambda = beta_square/(2.0-beta_square);
429 0 : do{
430 0 : costheta = (2.0*gRandom->Rndm())-1.;
431 0 : }
432 0 : while ( (1.0+lambda*costheta*costheta)<(2.0*gRandom->Rndm()) );
433 0 : sintheta = TMath::Sqrt((1. + costheta) * (1. - costheta));
434 0 : phi = 2.0 * TMath::ACos(-1.) * gRandom->Rndm();
435 0 : sinphi = TMath::Sin(phi);
436 0 : cosphi = TMath::Cos(phi);
437 0 : }
438 :
439 : // momentum vectors of electrons in virtual photon rest frame
440 0 : Double_t pProd1[3] = {p1 * sintheta * cosphi,
441 0 : p1 * sintheta * sinphi,
442 0 : p1 * costheta};
443 0 : Double_t pProd2[3] = {-1.0 * p1 * sintheta * cosphi,
444 0 : -1.0 * p1 * sintheta * sinphi,
445 0 : -1.0 * p1 * costheta};
446 0 : fProducts_dalitz[0].SetPx(pProd1[0]);
447 0 : fProducts_dalitz[0].SetPy(pProd1[1]);
448 0 : fProducts_dalitz[0].SetPz(pProd1[2]);
449 0 : fProducts_dalitz[0].SetE(e1);
450 0 : fProducts_dalitz[1].SetPx(pProd2[0]);
451 0 : fProducts_dalitz[1].SetPy(pProd2[1]);
452 0 : fProducts_dalitz[1].SetPz(pProd2[2]);
453 0 : fProducts_dalitz[1].SetE(e1);
454 :
455 : // third child kinematics in parent meson rest frame
456 0 : e3 = (pmass*pmass + realp_mass*realp_mass - epmass*epmass)/(2. * pmass);
457 0 : p3 = TMath::Sqrt((e3+realp_mass) * (e3-realp_mass));
458 :
459 : // third child 4-vector in parent meson rest frame
460 0 : costheta = (2.0 * gRandom->Rndm()) - 1.;
461 0 : sintheta = TMath::Sqrt((1. + costheta) * (1. - costheta));
462 0 : phi = 2.0 * TMath::ACos(-1.) * gRandom->Rndm();
463 0 : sinphi = TMath::Sin(phi);
464 0 : cosphi = TMath::Cos(phi);
465 0 : fProducts_dalitz[2].SetPx(p3 * sintheta * cosphi);
466 0 : fProducts_dalitz[2].SetPy(p3 * sintheta * sinphi);
467 0 : fProducts_dalitz[2].SetPz(p3 * costheta);
468 0 : fProducts_dalitz[2].SetE(e3);
469 :
470 : // boost the dielectron into the parent meson's rest frame
471 0 : Double_t eLPparent = TMath::Sqrt(p3*p3 + epmass*epmass);
472 0 : TVector3 boostPair( -1.0 * fProducts_dalitz[2].Px() / eLPparent,
473 0 : -1.0 * fProducts_dalitz[2].Py() / eLPparent,
474 0 : -1.0 * fProducts_dalitz[2].Pz() / eLPparent);
475 0 : fProducts_dalitz[0].Boost(boostPair);
476 0 : fProducts_dalitz[1].Boost(boostPair);
477 :
478 : // boost all decay products into the lab frame
479 0 : TVector3 boostLab(pparent->Px() / pparent->E(),
480 0 : pparent->Py() / pparent->E(),
481 0 : pparent->Pz() / pparent->E());
482 0 : fProducts_dalitz[0].Boost(boostLab);
483 0 : fProducts_dalitz[1].Boost(boostLab);
484 0 : fProducts_dalitz[2].Boost(boostLab);
485 :
486 0 : if(idpart==idPi0) {
487 0 : fProducts_pion[0]=fProducts_dalitz[0];
488 0 : fProducts_pion[1]=fProducts_dalitz[1];
489 0 : fProducts_pion[2]=fProducts_dalitz[2];
490 0 : }else if(idpart==idEta){
491 0 : fProducts_eta[0]=fProducts_dalitz[0];
492 0 : fProducts_eta[1]=fProducts_dalitz[1];
493 0 : fProducts_eta[2]=fProducts_dalitz[2];
494 0 : }else if(idpart==idOmega){
495 0 : fProducts_omega_dalitz[0]=fProducts_dalitz[0];
496 0 : fProducts_omega_dalitz[1]=fProducts_dalitz[1];
497 0 : fProducts_omega_dalitz[2]=fProducts_dalitz[2];
498 0 : }else if(idpart==idEtaPrime){
499 0 : fProducts_etaprime[0]=fProducts_dalitz[0];
500 0 : fProducts_etaprime[1]=fProducts_dalitz[1];
501 0 : fProducts_etaprime[2]=fProducts_dalitz[2];
502 0 : }else if(idpart==idPhi){
503 0 : fProducts_phi_dalitz[0]=fProducts_dalitz[0];
504 0 : fProducts_phi_dalitz[1]=fProducts_dalitz[1];
505 0 : fProducts_phi_dalitz[2]=fProducts_dalitz[2];
506 : }
507 :
508 0 : }
509 :
510 :
511 : //-----------------------------------------------------------------------------//
512 : // Generate 2-body resonance decays: Rho/Omega/Phi/JPsi //
513 : //-----------------------------------------------------------------------------//
514 :
515 0 : if(idpart==idRho||idpart==idOmega||idpart==idPhi||idpart==idJPsi){
516 :
517 :
518 : //get the parent mass
519 0 : mp_res = pparent->M();
520 :
521 : //check daughters mass
522 : md_res=emass;
523 0 : if ( mp_res < 2.*md_res ){
524 0 : printf("res into ee Decay kinematically impossible! \n");
525 0 : return;
526 : }
527 :
528 : // Sample the electron pair mass from a histogram and set Polarization
529 : for( ;; ) {
530 0 : if(idpart==idRho){
531 0 : epmass_res = fEPMassRho->GetRandom();
532 : PolPar=0.;
533 0 : }else if(idpart==idOmega){
534 0 : epmass_res = fEPMassOmega->GetRandom();
535 : PolPar=0.;
536 0 : }else if(idpart==idPhi){
537 0 : epmass_res = fEPMassPhi->GetRandom();
538 : PolPar=0.;
539 0 : }else if(idpart==idPhi){
540 0 : epmass_res = fEPMassPhi->GetRandom();
541 : PolPar=0.;
542 0 : }else if(idpart==idJPsi){
543 0 : epmass_res = fEPMassJPsi->GetRandom();
544 : PolPar=0.;
545 0 : }else{ printf(" Exodus ERROR: Resonance mass G-S parametrization not found \n");
546 0 : return;
547 : }
548 0 : if ( mp_res < 2.*epmass_res ) break;
549 : }
550 :
551 : // electron pair kinematics in virtual photon rest frame
552 0 : Ed_res = epmass_res/2.;
553 0 : pd_res = TMath::Sqrt((Ed_res+md_res)*(Ed_res-md_res));
554 :
555 : // momentum vectors of electrons in virtual photon rest frame
556 0 : fPol->SetParameter(0,PolPar);
557 0 : costheta = fPol->GetRandom();
558 0 : sintheta = TMath::Sqrt((1. + costheta)*(1. - costheta));
559 0 : fProducts_res[0].SetPx(pd_res * sintheta * cosphi);
560 0 : fProducts_res[0].SetPy(pd_res * sintheta * sinphi);
561 0 : fProducts_res[0].SetPz(pd_res * costheta);
562 0 : fProducts_res[0].SetE(Ed_res);
563 0 : fProducts_res[1].SetPx(-1.0 * pd_res * sintheta * cosphi);
564 0 : fProducts_res[1].SetPy(-1.0 * pd_res * sintheta * sinphi);
565 0 : fProducts_res[1].SetPz(-1.0 * pd_res * costheta);
566 0 : fProducts_res[1].SetE(Ed_res);
567 :
568 : // Beam parameters in LAB frame
569 0 : TLorentzVector pProj, pTarg;
570 : Double_t BeamE=3500.;
571 0 : pProj.SetPxPyPzE(0.,0.,-1.*BeamE,TMath::Sqrt(BeamE*BeamE+proton_mass*proton_mass)); // Beam 1
572 0 : pTarg.SetPxPyPzE(0.,0.,BeamE,TMath::Sqrt(BeamE*BeamE+proton_mass*proton_mass)); // Beam 2
573 :
574 : //re-build parent with G-S mass
575 0 : TLorentzVector pparent_corr;
576 0 : pparent_corr.SetPx(pparent->Px());
577 0 : pparent_corr.SetPy(pparent->Py());
578 0 : pparent_corr.SetPz(pparent->Pz());
579 0 : pparent_corr.SetE(sqrt(pow(pparent->P(),2)+pow(epmass_res,2)));
580 :
581 : //Boost Beam from CM to Resonance rest frame
582 0 : TVector3 betaResCM;
583 0 : betaResCM = (-1./pparent_corr.E()*pparent_corr.Vect());
584 0 : pProj.Boost(betaResCM);
585 0 : pTarg.Boost(betaResCM);
586 :
587 : //Define Zaxis in C-S frame and rotate legs to it
588 0 : TVector3 zaxisCS;
589 0 : zaxisCS=(((pProj.Vect()).Unit())-((pTarg.Vect()).Unit())).Unit();
590 0 : fProducts_res[0].RotateUz(zaxisCS);
591 0 : fProducts_res[1].RotateUz(zaxisCS);
592 :
593 : // boost decay products into the lab frame
594 0 : TVector3 boostLab_res_corr(pparent_corr.Px() / pparent_corr.E(),
595 0 : pparent_corr.Py() / pparent_corr.E(),
596 0 : pparent_corr.Pz() / pparent_corr.E());
597 0 : fProducts_res[0].Boost(boostLab_res_corr);
598 0 : fProducts_res[1].Boost(boostLab_res_corr);
599 :
600 0 : if(idpart==idRho) {
601 0 : fProducts_rho[0]=fProducts_res[0];
602 0 : fProducts_rho[1]=fProducts_res[1];
603 0 : }else if(idpart==idOmega){
604 0 : fProducts_omega[0]=fProducts_res[0];
605 0 : fProducts_omega[1]=fProducts_res[1];
606 0 : }else if(idpart==idPhi){
607 0 : fProducts_phi[0]=fProducts_res[0];
608 0 : fProducts_phi[1]=fProducts_res[1];
609 0 : }else if(idpart==idJPsi){
610 0 : fProducts_jpsi[0]=fProducts_res[0];
611 0 : fProducts_jpsi[1]=fProducts_res[1];
612 : }
613 :
614 0 : }
615 :
616 0 : return;
617 0 : }
618 :
619 : void AliDecayerExodus::Rot(Double_t pin[3], Double_t pout[3], Double_t costheta, Double_t sintheta,
620 : Double_t cosphi, Double_t sinphi) const
621 : {
622 : // Perform rotation
623 0 : pout[0] = pin[0]*costheta*cosphi-pin[1]*sinphi+pin[2]*sintheta*cosphi;
624 0 : pout[1] = pin[0]*costheta*sinphi+pin[1]*cosphi+pin[2]*sintheta*sinphi;
625 0 : pout[2] = -1.0 * pin[0] * sintheta + pin[2] * costheta;
626 0 : return;
627 : }
628 :
629 :
630 : Int_t AliDecayerExodus::ImportParticles(TClonesArray *particles)
631 : {
632 : //
633 : // Import particles for Dalitz and resonance decays
634 : //
635 :
636 : TClonesArray &clonesParticles = *particles;
637 :
638 : Int_t i, k;
639 : Double_t px, py, pz, e;
640 :
641 0 : Int_t pdgD [3][3] = { {kElectron, -kElectron, 22}, // pizero, eta, etaprime
642 : {kElectron, -kElectron, 111}, // omega dalitz
643 : {kElectron, -kElectron, 221} }; // phi dalitz
644 :
645 0 : Int_t pdgR [2] = {kElectron, -kElectron}; // rho, omega, phi, jpsi
646 :
647 :
648 :
649 0 : Int_t parentD[3] = { 0, 0, -1};
650 0 : Int_t dauD1 [3] = {-1, -1, 1};
651 0 : Int_t dauD2 [3] = {-1, -1, 2};
652 :
653 0 : Int_t parentR[2] = { 0, 0};
654 0 : Int_t dauR1 [2] = { -1, -1};
655 0 : Int_t dauR2 [2] = { -1, -1};
656 :
657 0 : for (Int_t j = 0; j < 9; j++){
658 :
659 : // pizero
660 0 : if(j==0){
661 0 : for (i = 2; i > -1; i--) {
662 0 : px = fProducts_pion[i].Px();
663 0 : py = fProducts_pion[i].Py();
664 0 : pz = fProducts_pion[i].Pz();
665 0 : e = fProducts_pion[i].E();
666 0 : new(clonesParticles[2 - i]) TParticle(pdgD[0][i], 1, parentD[i], -1, dauD1[i], dauD2[i], px, py, pz, e, 0., 0., 0., 0.);
667 : }
668 0 : return (3);
669 : }
670 :
671 : // rho
672 0 : else if(j==1){
673 0 : for (k = 1; k > -1; k--) {
674 0 : px = fProducts_rho[k].Px();
675 0 : py = fProducts_rho[k].Py();
676 0 : pz = fProducts_rho[k].Pz();
677 0 : e = fProducts_rho[k].E();
678 0 : new(clonesParticles[1 - k]) TParticle(pdgR[k], 1, parentR[k], -1, dauR1[k], dauR2[k], px, py, pz, e, 0., 0., 0., 0.);
679 : }
680 0 : return (2);
681 : }
682 :
683 : // eta
684 0 : else if(j==2){
685 0 : for (i = 2; i > -1; i--) {
686 0 : px = fProducts_eta[i].Px();
687 0 : py = fProducts_eta[i].Py();
688 0 : pz = fProducts_eta[i].Pz();
689 0 : e = fProducts_eta[i].E();
690 0 : new(clonesParticles[2 - i]) TParticle(pdgD[0][i], 1, parentD[i], -1, dauD1[i], dauD2[i], px, py, pz, e, 0., 0., 0., 0.);
691 : }
692 0 : return (3);
693 : }
694 :
695 : // omega dalitz
696 0 : else if(j==3){
697 0 : for (i = 2; i > -1; i--) {
698 0 : px = fProducts_omega_dalitz[i].Px();
699 0 : py = fProducts_omega_dalitz[i].Py();
700 0 : pz = fProducts_omega_dalitz[i].Pz();
701 0 : e = fProducts_omega_dalitz[i].E();
702 0 : new(clonesParticles[2 - i]) TParticle(pdgD[1][i], 1, parentD[i], -1, dauD1[i], dauD2[i], px, py, pz, e, 0., 0., 0., 0.);
703 : }
704 0 : return (3);
705 : }
706 :
707 : // omega direct
708 0 : else if(j==4){
709 0 : for (k = 1; k > -1; k--) {
710 0 : px = fProducts_rho[k].Px();
711 0 : py = fProducts_rho[k].Py();
712 0 : pz = fProducts_rho[k].Pz();
713 0 : e = fProducts_rho[k].E();
714 0 : new(clonesParticles[1 - k]) TParticle(pdgR[k], 1, parentR[k], -1, dauR1[k], dauR2[k], px, py, pz, e, 0., 0., 0., 0.);
715 : }
716 0 : return (2);
717 : }
718 :
719 : // etaprime
720 0 : else if(j==5){
721 0 : for (i = 2; i > -1; i--) {
722 0 : px = fProducts_etaprime[i].Px();
723 0 : py = fProducts_etaprime[i].Py();
724 0 : pz = fProducts_etaprime[i].Pz();
725 0 : e = fProducts_etaprime[i].E();
726 0 : new(clonesParticles[2 - i]) TParticle(pdgD[0][i], 1, parentD[i], -1, dauD1[i], dauD2[i], px, py, pz, e, 0., 0., 0., 0.);
727 : }
728 0 : return (3);
729 : }
730 :
731 : // phi dalitz
732 0 : else if(j==6){
733 0 : for (i = 2; i > -1; i--) {
734 0 : px = fProducts_phi_dalitz[i].Px();
735 0 : py = fProducts_phi_dalitz[i].Py();
736 0 : pz = fProducts_phi_dalitz[i].Pz();
737 0 : e = fProducts_phi_dalitz[i].E();
738 0 : new(clonesParticles[2 - i]) TParticle(pdgD[2][i], 1, parentD[i], -1, dauD1[i], dauD2[i], px, py, pz, e, 0., 0., 0., 0.);
739 : }
740 0 : return (3);
741 : }
742 :
743 :
744 : // phi direct
745 0 : else if(j==7){
746 0 : for (k = 1; k > -1; k--) {
747 0 : px = fProducts_phi[k].Px();
748 0 : py = fProducts_phi[k].Py();
749 0 : pz = fProducts_phi[k].Pz();
750 0 : e = fProducts_phi[k].E();
751 0 : new(clonesParticles[1 - k]) TParticle(pdgR[k], 1, parentR[k], -1, dauR1[k], dauR2[k], px, py, pz, e, 0., 0., 0., 0.);
752 : }
753 0 : return (2);
754 : }
755 :
756 : // jpsi direct
757 0 : else if(j==8){
758 0 : for (k = 1; k > -1; k--) {
759 0 : px = fProducts_jpsi[k].Px();
760 0 : py = fProducts_jpsi[k].Py();
761 0 : pz = fProducts_jpsi[k].Pz();
762 0 : e = fProducts_jpsi[k].E();
763 0 : new(clonesParticles[1 - k]) TParticle(pdgR[k], 1, parentR[k], -1, dauR1[k], dauR2[k], px, py, pz, e, 0., 0., 0., 0.);
764 : }
765 0 : return (2);
766 : }
767 :
768 : }
769 :
770 0 : return particles->GetEntries();
771 :
772 0 : }
773 :
774 :
775 : void AliDecayerExodus::Decay(TClonesArray *array)
776 : {
777 : // Replace all Dalitz(pi0,eta,omega,eta',phi) and resonance(rho,omega,phi,jpsi) decays with the correct matrix element decays
778 : // for di-electron cocktail calculations
779 :
780 :
781 0 : Int_t nt = array->GetEntriesFast();
782 0 : TParticle* dp3[3];
783 0 : TParticle* dp2[2];
784 : Int_t fd3, ld3, fd2, ld2, fd, ld;
785 : Int_t j, k;
786 :
787 0 : for (Int_t i = 0; i < nt; i++) {
788 0 : TParticle* part = (TParticle*) (array->At(i));
789 0 : if (part->GetPdgCode() != 111 || part->GetPdgCode() != 221 || part->GetPdgCode() != 223 || part->GetPdgCode() != 331 || part->GetPdgCode() != 333 || part->GetPdgCode() != 443 ) continue;
790 :
791 : //
792 : // Pizero Dalitz
793 : //
794 0 : if(part->GetPdgCode() == 111){
795 :
796 0 : fd3 = part->GetFirstDaughter() - 1;
797 0 : ld3 = part->GetLastDaughter() - 1;
798 :
799 0 : if (fd3 < 0) continue;
800 0 : if ((ld3 - fd3) != 2) continue;
801 :
802 0 : for (j = 0; j < 3; j++) dp3[j] = (TParticle*) (array->At(fd3+j));
803 :
804 0 : if((dp3[0]->GetPdgCode() != 22) && (TMath::Abs(dp3[1]->GetPdgCode()) != 11)) continue;
805 :
806 0 : TLorentzVector Pizero(part->Px(), part->Py(), part->Pz(), part->Energy());
807 0 : Decay(111, &Pizero);
808 0 : for (j = 0; j < 3; j++) dp3[j]->SetMomentum(fProducts_pion[2-j]);
809 0 : }
810 :
811 :
812 : //
813 : // Eta Dalitz
814 : //
815 :
816 0 : if(part->GetPdgCode() == 221){
817 :
818 0 : fd3 = part->GetFirstDaughter() - 1;
819 0 : ld3 = part->GetLastDaughter() - 1;
820 :
821 0 : if (fd3 < 0) continue;
822 0 : if ((ld3 - fd3) != 2) continue;
823 :
824 0 : for (j = 0; j < 3; j++) dp3[j] = (TParticle*) (array->At(fd3+j));
825 :
826 0 : if((dp3[0]->GetPdgCode() != 22) && ((TMath::Abs(dp3[1]->GetPdgCode()) != 11))) continue;
827 :
828 0 : TLorentzVector Eta(part->Px(), part->Py(), part->Pz(), part->Energy());
829 0 : Decay(221, &Eta);
830 0 : for (j = 0; j < 3; j++) dp3[j]->SetMomentum(fProducts_eta[2-j]);
831 0 : }
832 :
833 : //
834 : // Rho
835 : //
836 :
837 0 : if(part->GetPdgCode() == 113){
838 :
839 0 : fd2 = part->GetFirstDaughter() - 1;
840 0 : ld2 = part->GetLastDaughter() - 1;
841 :
842 0 : if (fd2 < 0) continue;
843 0 : if ((ld2 - fd2) != 1) continue;
844 :
845 0 : for (k = 0; k < 2; k++) dp2[k] = (TParticle*) (array->At(fd2+k));
846 :
847 0 : if((dp2[0]->GetPdgCode() != 11) && ((TMath::Abs(dp2[1]->GetPdgCode()) != 11))) continue;
848 :
849 0 : TLorentzVector Rho(part->Px(), part->Py(), part->Pz(), part->Energy());
850 0 : Decay(113, &Rho);
851 0 : for (k = 0; k < 2; k++) dp2[k]->SetMomentum(fProducts_rho[1-k]);
852 0 : }
853 :
854 : //
855 : // Omega dalitz and direct
856 : //
857 :
858 0 : if(part->GetPdgCode() == 223){
859 :
860 0 : fd = part->GetFirstDaughter() - 1;
861 0 : ld = part->GetLastDaughter() - 1;
862 :
863 0 : if (fd < 0) continue;
864 :
865 0 : if ((ld - fd) == 2){
866 :
867 0 : for (j = 0; j < 3; j++) dp3[j] = (TParticle*) (array->At(fd+j));
868 0 : if( dp3[0]->GetPdgCode() != 111 && (TMath::Abs(dp3[1]->GetPdgCode()) != 11)) continue;
869 :
870 0 : TLorentzVector Omegadalitz(part->Px(), part->Py(), part->Pz(), part->Energy());
871 0 : Decay(223, &Omegadalitz);
872 0 : for (j = 0; j < 3; j++) dp3[j]->SetMomentum(fProducts_omega_dalitz[2-j]);
873 0 : }
874 :
875 0 : else if ((ld - fd) == 1) {
876 :
877 0 : for (k = 0; k < 2; k++) dp2[k] = (TParticle*) (array->At(fd+k));
878 0 : if( dp2[0]->GetPdgCode() != 11 && (TMath::Abs(dp2[1]->GetPdgCode()) != 11)) continue;
879 :
880 0 : TLorentzVector Omega(part->Px(), part->Py(), part->Pz(), part->Energy());
881 0 : Decay(223, &Omega);
882 0 : for (k = 0; k < 2; k++) dp2[k]->SetMomentum(fProducts_omega[1-k]);
883 0 : }
884 : }
885 :
886 : //
887 : // Etaprime dalitz
888 : //
889 :
890 0 : if(part->GetPdgCode() == 331){
891 :
892 0 : fd3 = part->GetFirstDaughter() - 1;
893 0 : ld3 = part->GetLastDaughter() - 1;
894 :
895 0 : if (fd3 < 0) continue;
896 0 : if ((ld3 - fd3) != 2) continue;
897 :
898 0 : for (j = 0; j < 3; j++) dp3[j] = (TParticle*) (array->At(fd3+j));
899 :
900 0 : if((dp3[0]->GetPdgCode() != 22) && ((TMath::Abs(dp3[1]->GetPdgCode()) != 11))) continue;
901 :
902 0 : TLorentzVector Etaprime(part->Px(), part->Py(), part->Pz(), part->Energy());
903 0 : Decay(331, &Etaprime);
904 0 : for (j = 0; j < 3; j++) dp3[j]->SetMomentum(fProducts_etaprime[2-j]);
905 0 : }
906 :
907 : //
908 : // Phi dalitz and direct
909 : //
910 0 : if(part->GetPdgCode() == 333){
911 :
912 0 : fd = part->GetFirstDaughter() - 1;
913 0 : ld = part->GetLastDaughter() - 1;
914 :
915 0 : if (fd < 0) continue;
916 0 : if ((ld - fd) == 2){
917 0 : for (j = 0; j < 3; j++) dp3[j] = (TParticle*) (array->At(fd+j));
918 0 : if( dp3[0]->GetPdgCode() != 221 && (TMath::Abs(dp3[1]->GetPdgCode()) != 11)) continue;
919 :
920 0 : TLorentzVector Phidalitz(part->Px(), part->Py(), part->Pz(), part->Energy());
921 0 : Decay(333, &Phidalitz);
922 0 : for (j = 0; j < 3; j++) dp3[j]->SetMomentum(fProducts_phi_dalitz[2-j]);
923 0 : }
924 :
925 0 : else if ((ld - fd) == 1) {
926 0 : for (k = 0; k < 2; k++) dp2[k] = (TParticle*) (array->At(fd+k));
927 0 : if( dp2[0]->GetPdgCode() != 11 && (TMath::Abs(dp2[1]->GetPdgCode()) != 11)) continue;
928 :
929 0 : TLorentzVector Phi(part->Px(), part->Py(), part->Pz(), part->Energy());
930 0 : Decay(333, &Phi);
931 0 : for (k = 0; k < 2; k++) dp2[k]->SetMomentum(fProducts_phi[1-k]);
932 0 : }
933 : }
934 :
935 : //
936 : // JPsi
937 : //
938 :
939 0 : if(part->GetPdgCode() == 443){
940 :
941 0 : fd2 = part->GetFirstDaughter() - 1;
942 0 : ld2 = part->GetLastDaughter() - 1;
943 :
944 0 : if (fd2 < 0) continue;
945 0 : if ((ld2 - fd2) != 1) continue;
946 :
947 0 : for (k = 0; k < 2; k++) dp2[k] = (TParticle*) (array->At(fd2+k));
948 :
949 0 : if((dp2[0]->GetPdgCode() != 11) && ((TMath::Abs(dp2[1]->GetPdgCode()) != 11))) continue;
950 :
951 0 : TLorentzVector JPsi(part->Px(), part->Py(), part->Pz(), part->Energy());
952 0 : Decay(443, &JPsi);
953 0 : for (k = 0; k < 2; k++) dp2[k]->SetMomentum(fProducts_jpsi[1-k]);
954 0 : }
955 :
956 0 : }
957 0 : }
958 :
959 :
960 : AliDecayerExodus& AliDecayerExodus::operator=(const AliDecayerExodus& rhs)
961 : {
962 : // Assignment operator
963 0 : rhs.Copy(*this);
964 0 : return *this;
965 : }
966 :
967 : void AliDecayerExodus::Copy(TObject&) const
968 : {
969 : //
970 : // Copy
971 : //
972 0 : Fatal("Copy","Not implemented!\n");
973 0 : }
974 :
975 :
976 0 : AliDecayerExodus::AliDecayerExodus(const AliDecayerExodus &decayer)
977 0 : : AliDecayer(),
978 0 : fEPMassPion(0),
979 0 : fEPMassEta(0),
980 0 : fEPMassEtaPrime(0),
981 0 : fEPMassRho(0),
982 0 : fEPMassOmega(0),
983 0 : fEPMassOmegaDalitz(0),
984 0 : fEPMassPhi(0),
985 0 : fEPMassPhiDalitz(0),
986 0 : fEPMassJPsi(0),
987 0 : fInit(0)
988 0 : {
989 : // Copy Constructor
990 0 : decayer.Copy(*this);
991 0 : }
992 :
993 :
994 :
|
