1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 #include "AliDecayerExodus.h"
18 #include <Riostream.h>
23 #include <TParticle.h>
24 #include <TDatabasePDG.h>
26 #include <TLorentzVector.h>
27 #include <TClonesArray.h>
30 ClassImp(AliDecayerExodus)
32 //---------------------------------------------------------------------------------------------------
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
39 // For the electromagnetic form factor the parameterization from
40 // Lepton-G is used: L.G. Landsberg et al.: Phys. Rep. 128(1985)301
42 // Ralf Averbeck (R.Averbeck@gsi.de)
43 // Irem Erdemir (irem.erdemir@cern.ch)
45 //---------------------------------------------------------------------------------------------------
48 AliDecayerExodus::AliDecayerExodus():
55 fEPMassOmegaDalitz(0),
66 void AliDecayerExodus::Init()
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;
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;
88 //================================================================================//
89 // Create electron pair mass histograms from dalitz decays //
90 //================================================================================//
92 // Get the particle masses
96 pionmass = (TDatabasePDG::Instance()->GetParticle(111))->Mass();
97 etamass = (TDatabasePDG::Instance()->GetParticle(221))->Mass();
98 omegamass = (TDatabasePDG::Instance()->GetParticle(223))->Mass();
99 etaprimemass = (TDatabasePDG::Instance()->GetParticle(331))->Mass();
100 phimass = (TDatabasePDG::Instance()->GetParticle(333))->Mass();
102 emass = (TDatabasePDG::Instance()->GetParticle(11))->Mass();
103 // child - other : third childs from Dalitz decays
104 omasspion = pionmass;
109 maxpion = pionmass - omassgamma;
110 maxeta = etamass - omassgamma;
111 maxomega = omegamass - pionmass;
112 maxetaprime = etaprimemass - omassgamma;
113 maxphi = phimass - omasseta;
115 binwidth_pion = (maxpion - min) / (Double_t)nbins;
116 binwidth_eta = (maxeta - min) / (Double_t)nbins;
117 binwidth_omega = (maxomega - min) / (Double_t)nbins;
118 binwidth_etaprime = (maxetaprime - min) / (Double_t)nbins;
119 binwidth_phi = (maxphi - min) / (Double_t)nbins;
122 epsilon_pion = (emass / pionmass) * (emass / pionmass);
123 epsilon_eta = (emass / etamass) * (emass / etamass);
124 epsilon_omega = (emass / omegamass) * (emass / omegamass);
125 epsilon_etaprime = (emass / etaprimemass) * (emass / etaprimemass);
126 epsilon_phi = (emass / phimass) * (emass / phimass);
129 delta_pion = (omassgamma / pionmass) * (omassgamma / pionmass);
130 delta_eta = (omassgamma / etamass) * (omassgamma / etamass);
131 delta_omega = (omasspion / omegamass) * (omasspion / omegamass);
132 delta_etaprime = (omassgamma / etaprimemass) * (omassgamma / etaprimemass);
133 delta_phi = (omasseta / phimass) * (omasseta / phimass);
137 // create pair mass histograms for Dalitz decays of pi0, eta, omega, eta' and phi
138 if (!fEPMassPion) {delete fEPMassPion; fEPMassPion = new TH1F("fEPMassPion", "Dalitz electron pair from pion", nbins, min, maxpion); }
139 if (!fEPMassEta) {delete fEPMassEta; fEPMassEta = new TH1F("fEPMassEta", "Dalitz electron pair from eta", nbins, min, maxeta);}
140 if (!fEPMassOmegaDalitz) {delete fEPMassOmegaDalitz; fEPMassOmegaDalitz = new TH1F("fEPMassOmegaDalitz", "Dalitz electron pair from omega ", nbins, min, maxomega);}
141 if (!fEPMassEtaPrime) {delete fEPMassEtaPrime; fEPMassEtaPrime = new TH1F("fEPMassEtaPrime", "Dalitz electron pair from etaprime", nbins, min, maxetaprime);}
142 if (!fEPMassPhiDalitz) {delete fEPMassPhiDalitz; fEPMassPhiDalitz = new TH1F("fEPMassPhiDalitz", "Dalitz electron pair from phi ", nbins, min, maxphi);}
145 mLL_pion = mLL_eta = mLL_omega = mLL_etaprime = mLL_phi = 0.;
147 for (ibin = 1; ibin <= nbins; ibin++ )
149 mLL_pion = min + (Double_t)(ibin - 1) * binwidth_pion + binwidth_pion / 2.0;
150 mLL_eta = min + (Double_t)(ibin - 1) * binwidth_eta + binwidth_eta / 2.0;
151 mLL_omega = min + (Double_t)(ibin - 1) * binwidth_omega + binwidth_omega / 2.0;
152 mLL_etaprime = min + (Double_t)(ibin - 1) * binwidth_etaprime + binwidth_etaprime / 2.0;
153 mLL_phi = min + (Double_t)(ibin - 1) * binwidth_phi + binwidth_phi / 2.0;
156 q_pion = (mLL_pion / pionmass) * (mLL_pion / pionmass);
157 q_eta = (mLL_eta / etamass) * (mLL_eta / etamass);
158 q_omega = (mLL_omega / omegamass)*(mLL_omega / omegamass);
159 q_etaprime = (mLL_etaprime / etaprimemass) * (mLL_etaprime / etaprimemass);
160 q_phi = (mLL_phi / phimass) * (mLL_phi / phimass);
162 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 AliFatal("Error in calculating Dalitz mass histogram binning!");
168 kwHelp_pion = (1.0 + q_pion / (1.0 - delta_pion)) * (1.0 + q_pion / (1.0 - delta_pion))
169 - 4.0 * q_pion / ((1.0 - delta_pion) * (1.0 - delta_pion));
171 kwHelp_eta = (1.0 + q_eta / (1.0 - delta_eta)) * (1.0 + q_eta / (1.0 - delta_eta))
172 - 4.0 * q_eta / ((1.0 - delta_eta) * (1.0 - delta_eta));
174 kwHelp_omega = (1.0 + q_omega / (1.0 - delta_omega)) * (1.0 + q_omega / (1.0 - delta_omega))
175 - 4.0 * q_omega / ((1.0 - delta_omega) * (1.0 - delta_omega));
177 kwHelp_etaprime = (1.0 + q_etaprime / (1.0 - delta_etaprime)) * (1.0 + q_etaprime / (1.0 - delta_etaprime))
178 - 4.0 * q_etaprime / ((1.0 - delta_etaprime) * (1.0 - delta_etaprime));
180 kwHelp_phi = (1.0 + q_phi / (1.0 - delta_phi)) * (1.0 + q_phi / (1.0 - delta_phi))
181 - 4.0 * q_phi / ((1.0 - delta_phi) * (1.0 - delta_phi));
186 if ( kwHelp_pion <= 0.0 || kwHelp_eta <= 0.0 || kwHelp_omega <= 0.0 || kwHelp_etaprime <= 0.0 || kwHelp_phi <= 0.0 )
188 AliFatal("Error in calculating Dalitz mass histogram binning!");
193 // Invariant mass distributions of electron pairs from Dalitz decays
194 // using Kroll-Wada function
196 krollWada_pion = (2.0 / mLL_pion) * TMath::Exp(1.5 * TMath::Log(kwHelp_pion))
197 * TMath::Sqrt(1.0 - 4.0 * epsilon_pion / q_pion)
198 * (1.0 + 2.0 * epsilon_pion / q_pion);
201 krollWada_eta = (2.0 / mLL_eta) * TMath::Exp(1.5 * TMath::Log(kwHelp_eta))
202 * TMath::Sqrt(1.0 - 4.0 * epsilon_eta / q_eta)
203 * (1.0 + 2.0 * epsilon_eta / q_eta);
206 krollWada_omega = (2.0 / mLL_omega) * TMath::Exp(1.5 * TMath::Log(kwHelp_omega))
207 * TMath::Sqrt(1.0 - 4.0 * epsilon_omega / q_omega)
208 * (1.0 + 2.0 * epsilon_omega / q_omega);
211 krollWada_etaprime = (2.0 / mLL_etaprime) * TMath::Exp(1.5 * TMath::Log(kwHelp_etaprime))
212 * TMath::Sqrt(1.0 - 4.0 * epsilon_etaprime / q_etaprime)
213 * (1.0 + 2.0 * epsilon_etaprime / q_etaprime);
215 krollWada_phi = (2.0 / mLL_phi) * TMath::Exp(1.5 * TMath::Log(kwHelp_phi))
216 * TMath::Sqrt(1.0 - 4.0 * epsilon_phi / q_phi)
217 * (1.0 + 2.0 * epsilon_phi / q_phi);
221 // Form factors from Lepton-G
222 formFactor_pion = 1.0/(1.0-5.5*mLL_pion*mLL_pion);
223 formFactor_eta = 1.0/(1.0-1.9*mLL_eta*mLL_eta);
224 formFactor_omega = (TMath::Power(TMath::Power(0.6519,2),2))
225 /(TMath::Power(TMath::Power(0.6519,2)-TMath::Power(mLL_omega, 2), 2)
226 + TMath::Power(0.04198, 2)*TMath::Power(0.6519, 2));
227 formFactor_etaprime = (TMath::Power(TMath::Power(0.764,2),2))
228 /(TMath::Power(TMath::Power(0.764,2)-TMath::Power(mLL_etaprime, 2), 2)
229 + TMath::Power(0.1020, 2)*TMath::Power(0.764, 2));
230 formFactor_phi = 1.0;
235 weight_pion = krollWada_pion * formFactor_pion * formFactor_pion;
236 weight_eta = krollWada_eta * formFactor_eta * formFactor_eta;
237 weight_omega_dalitz = krollWada_omega * formFactor_omega;
238 weight_etaprime = krollWada_etaprime * formFactor_etaprime;
239 weight_phi_dalitz = krollWada_phi * formFactor_phi * formFactor_phi;
242 // Fill histograms of electron pair masses from dalitz decays
243 fEPMassPion ->AddBinContent(ibin, weight_pion);
244 fEPMassEta ->AddBinContent(ibin, weight_eta);
245 fEPMassOmegaDalitz->AddBinContent(ibin, weight_omega_dalitz);
246 fEPMassEtaPrime ->AddBinContent(ibin, weight_etaprime);
247 fEPMassPhiDalitz ->AddBinContent(ibin, weight_phi_dalitz);
253 //===================================================================================//
254 // Create electron pair mass histograms from resonance decays //
255 //===================================================================================//
257 Double_t pimass = 0.13956995;
259 // Get the particle masses
261 vmass_rho = (TDatabasePDG::Instance()->GetParticle(113))->Mass();
262 vmass_omega = (TDatabasePDG::Instance()->GetParticle(223))->Mass();
263 vmass_phi = (TDatabasePDG::Instance()->GetParticle(333))->Mass();
264 vmass_jpsi = (TDatabasePDG::Instance()->GetParticle(443))->Mass();
265 // Get the particle widths
267 vwidth_rho = (TDatabasePDG::Instance()->GetParticle(113))->Width();
268 vwidth_omega = (TDatabasePDG::Instance()->GetParticle(223))->Width();
269 vwidth_phi = (TDatabasePDG::Instance()->GetParticle(333))->Width();
270 vwidth_jpsi = (TDatabasePDG::Instance()->GetParticle(443))->Width();
273 if ( mass_min == 0. && mass_max == 0. )
275 mass_min = 2.*pimass;
279 binwidth = (mass_max-mass_min)/(Double_t)nbins;
281 // create pair mass histograms for resonances of rho, omega, phi and jpsi
282 if (!fEPMassRho) {delete fEPMassRho; fEPMassRho = new TH1F("fEPMassRho","mass rho",nbins,mass_min,mass_max);}
283 if (!fEPMassOmega) {delete fEPMassOmega; fEPMassOmega = new TH1F("fEPMassOmega","mass omega",nbins,mass_min,mass_max);}
284 if (!fEPMassPhi) {delete fEPMassPhi; fEPMassPhi = new TH1F("fEPMassPhi","mass phi",nbins,mass_min,mass_max);}
285 if (!fEPMassJPsi) {delete fEPMassJPsi; fEPMassJPsi = new TH1F("fEPMassJPsi","mass jpsi",nbins,mass_min,mass_max);}
288 for (ibin=1; ibin<=nbins; ibin++ )
290 mass_bin = mass_min+(Double_t)(ibin-1)*binwidth+binwidth/2.0;
292 weight_rho = (Float_t)GounarisSakurai(mass_bin,vmass_rho,vwidth_rho,emass);
293 weight_omega = (Float_t)GounarisSakurai(mass_bin,vmass_omega,vwidth_omega,emass);
294 weight_phi = (Float_t)GounarisSakurai(mass_bin,vmass_phi,vwidth_phi,emass);
295 weight_jpsi = (Float_t)Lorentz(mass_bin,vmass_jpsi,vwidth_jpsi);
297 // Fill histograms of electron pair masses from resonance decays
298 fEPMassRho ->AddBinContent(ibin,weight_rho);
299 fEPMassOmega->AddBinContent(ibin,weight_omega);
300 fEPMassPhi ->AddBinContent(ibin,weight_phi);
301 fEPMassJPsi ->AddBinContent(ibin,weight_jpsi);
306 Double_t AliDecayerExodus::GounarisSakurai(Float_t mass, Double_t vmass, Double_t vwidth, Double_t emass)
308 // Invariant mass distributions of electron pairs from resonance decays
309 // of rho, omega and phi
310 // using Gounaris-Sakurai function
313 Double_t epsilon = 0.;
314 Double_t weight = 0.;
316 Double_t pimass = 0.13956995;
318 corr = vwidth*(vmass/mass)*exp(1.5*log((mass*mass/4.0-pimass*pimass)
319 /(vmass*vmass/4.0-pimass*pimass)));
321 epsilon = (emass/mass)*(emass/mass);
323 if ( 1.0-4.0*epsilon>=0.0 )
325 weight = sqrt(1.0-4.0*epsilon)*(1.0+2.0*epsilon)/
326 ((vmass*vmass-mass*mass)*(vmass*vmass-mass*mass)+
327 (vmass*corr)*(vmass*corr));
333 Double_t AliDecayerExodus::Lorentz(Float_t mass, Double_t vmass, Double_t vwidth)
335 // Invariant mass distributions of electron pairs from resonance decay
336 // of jpsi (and it can also be used for other particles except rho, omega and phi)
337 // using Lorentz function
341 weight = (vwidth*vwidth/4.0)/(vwidth*vwidth/4.0+(vmass-mass)*(vmass-mass));
347 void AliDecayerExodus::Decay(Int_t idpart, TLorentzVector* pparent)
356 Double_t pmass_pion, pmass_eta, pmass_omega_dalitz, pmass_etaprime, pmass_phi_dalitz;
357 Double_t emass, omass_pion, omass_eta, omass_gamma, epmass_pion, epmass_eta, epmass_omega_dalitz, epmass_etaprime, epmass_phi_dalitz;
358 Double_t e1_pion, e1_eta, e1_omega, e1_etaprime, e1_phi;
359 Double_t p1_pion, p1_eta, p1_omega, p1_etaprime, p1_phi;
360 Double_t e3_gamma_pion, e3_gamma_eta, e3_pion, e3_gamma_etaprime, e3_eta;
361 Double_t p3_gamma_pion, p3_gamma_eta, p3_pion, p3_gamma_etaprime, p3_eta;
363 Double_t wp_rho, wp_omega, wp_phi, wp_jpsi, epmass_rho, epmass_omega, epmass_phi, epmass_jpsi;
364 Double_t mp_rho, mp_omega, mp_phi, mp_jpsi, md_rho, md_omega, md_phi, md_jpsi;
365 Double_t Ed_rho, Ed_omega, Ed_phi, Ed_jpsi, pd_rho, pd_omega, pd_phi, pd_jpsi;
368 md_rho = md_omega = md_phi = md_jpsi = 0.;
371 Double_t costheta, sintheta, cosphi, sinphi, phi;
373 // Get the particle masses of daughters
374 emass = (TDatabasePDG::Instance()->GetParticle(11)) ->Mass();
375 omass_pion = (TDatabasePDG::Instance()->GetParticle(111))->Mass();
376 omass_eta = (TDatabasePDG::Instance()->GetParticle(221))->Mass();
377 omass_gamma = (TDatabasePDG::Instance()->GetParticle(22)) ->Mass();
379 // Get the particle widths of mothers for resonances
380 wp_rho = (TDatabasePDG::Instance()->GetParticle(113))->Width();
381 wp_omega = (TDatabasePDG::Instance()->GetParticle(223))->Width();
382 wp_phi = (TDatabasePDG::Instance()->GetParticle(333))->Width();
383 wp_jpsi = (TDatabasePDG::Instance()->GetParticle(443))->Width();
385 costheta = (2.0 * gRandom->Rndm()) - 1.;
386 sintheta = TMath::Sqrt((1. + costheta) * (1. - costheta));
387 phi = 2.0 * TMath::ACos(-1.) * gRandom->Rndm();
388 sinphi = TMath::Sin(phi);
389 cosphi = TMath::Cos(phi);
392 //-----------------------------------------------------------------------------//
393 // Generate Pizero Dalitz decay //
394 //-----------------------------------------------------------------------------//
397 pmass_pion = pparent->M();
400 // Sample the electron pair mass from a histogram
401 epmass_pion = fEPMassPion->GetRandom();
402 if (pmass_pion-omass_gamma>epmass_pion && epmass_pion/2.>emass) break;
405 // electron pair kinematics in virtual photon rest frame
406 e1_pion = epmass_pion / 2.;
407 p1_pion = TMath::Sqrt((e1_pion + emass) * (e1_pion - emass));
409 // momentum vectors of electrons in virtual photon rest frame
410 Double_t pProd1_pion[3] = {p1_pion * sintheta * cosphi,
411 p1_pion * sintheta * sinphi,
414 Double_t pProd2_pion[3] = {-1.0 * p1_pion * sintheta * cosphi,
415 -1.0 * p1_pion * sintheta * sinphi,
416 -1.0 * p1_pion * costheta};
419 // third child kinematics in parent meson rest frame
420 e3_gamma_pion = (pmass_pion * pmass_pion - epmass_pion * epmass_pion)/(2. * pmass_pion);
421 p3_gamma_pion = TMath::Sqrt((e3_gamma_pion * e3_gamma_pion));
424 // third child 4-vector in parent meson rest frame
425 fProducts_pion[2].SetPx(p3_gamma_pion * sintheta * cosphi);
426 fProducts_pion[2].SetPy(p3_gamma_pion * sintheta * sinphi);
427 fProducts_pion[2].SetPz(p3_gamma_pion * costheta);
428 fProducts_pion[2].SetE(e3_gamma_pion);
431 // electron 4-vectors in properly rotated virtual photon rest frame
432 Double_t pRot1_pion[3] = {0.};
433 Rot(pProd1_pion, pRot1_pion, costheta, -sintheta, -cosphi, -sinphi);
434 Double_t pRot2_pion[3] = {0.};
435 Rot(pProd2_pion, pRot2_pion, costheta, -sintheta, -cosphi, -sinphi);
436 fProducts_pion[0].SetPx(pRot1_pion[0]);
437 fProducts_pion[0].SetPy(pRot1_pion[1]);
438 fProducts_pion[0].SetPz(pRot1_pion[2]);
439 fProducts_pion[0].SetE(e1_pion);
440 fProducts_pion[1].SetPx(pRot2_pion[0]);
441 fProducts_pion[1].SetPy(pRot2_pion[1]);
442 fProducts_pion[1].SetPz(pRot2_pion[2]);
443 fProducts_pion[1].SetE(e1_pion);
445 // boost the dielectron into the parent meson's rest frame
446 Double_t eLPparent_pion = TMath::Sqrt(p3_gamma_pion * p3_gamma_pion + epmass_pion * epmass_pion);
447 TVector3 boostPair_pion( -1.0 * fProducts_pion[2].Px() / eLPparent_pion,
448 -1.0 * fProducts_pion[2].Py() / eLPparent_pion,
449 -1.0 * fProducts_pion[2].Pz() / eLPparent_pion);
450 fProducts_pion[0].Boost(boostPair_pion);
451 fProducts_pion[1].Boost(boostPair_pion);
453 // boost all decay products into the lab frame
454 TVector3 boostLab_pion(pparent->Px() / pparent->E(),
455 pparent->Py() / pparent->E(),
456 pparent->Pz() / pparent->E());
458 fProducts_pion[0].Boost(boostLab_pion);
459 fProducts_pion[1].Boost(boostLab_pion);
460 fProducts_pion[2].Boost(boostLab_pion);
465 //-----------------------------------------------------------------------------//
466 // Generate Rho resonance decay //
467 //-----------------------------------------------------------------------------//
469 else if(idpart==113){
470 // calculate rho mass
472 mp_rho = pparent->M();
475 Double_t x_rho=pparent->Px(); Double_t y_rho=pparent->Py(); Double_t z_rho=pparent->Pz();
476 Double_t t_rho=pparent->E();
477 Double_t p_rho=x_rho*x_rho+y_rho*y_rho+z_rho*z_rho;
478 Double_t Q2_rho=abs((t_rho*t_rho)-(p_rho*p_rho));
479 mp_rho = sqrt(Q2_rho);
482 if ( mp_rho < 2.*md_rho )
484 printf("Rho into ee Decay kinematically impossible! \n");
489 // Sample the electron pair mass from a histogram
490 epmass_rho = fEPMassRho->GetRandom();
491 if ( mp_rho < 2.*epmass_rho ) break;
494 // electron pair kinematics in virtual photon rest frame
495 Ed_rho = epmass_rho/2.;
496 pd_rho = TMath::Sqrt((Ed_rho+md_rho)*(Ed_rho-md_rho));
498 // momentum vectors of electrons in virtual photon rest frame
499 Double_t pProd1_rho[3] = {pd_rho * sintheta * cosphi,
500 pd_rho * sintheta * sinphi,
503 Double_t pProd2_rho[3] = {-1.0 * pd_rho * sintheta * cosphi,
504 -1.0 * pd_rho * sintheta * sinphi,
505 -1.0 * pd_rho * costheta};
508 // electron 4 vectors in properly rotated virtual photon rest frame
509 Double_t pRot1_rho[3] = {0.};
510 Rot(pProd1_rho, pRot1_rho, costheta, -sintheta, -cosphi, -sinphi);
511 Double_t pRot2_rho[3] = {0.};
512 Rot(pProd2_rho, pRot2_rho, costheta, -sintheta, -cosphi, -sinphi);
513 fProducts_rho[0].SetPx(pRot1_rho[0]);
514 fProducts_rho[0].SetPy(pRot1_rho[1]);
515 fProducts_rho[0].SetPz(pRot1_rho[2]);
516 fProducts_rho[0].SetE(Ed_rho);
517 fProducts_rho[1].SetPx(pRot2_rho[0]);
518 fProducts_rho[1].SetPy(pRot2_rho[1]);
519 fProducts_rho[1].SetPz(pRot2_rho[2]);
520 fProducts_rho[1].SetE(Ed_rho);
523 // boost decay products into the lab frame
524 TVector3 boostLab_rho(pparent->Px() / pparent->E(),
525 pparent->Py() / pparent->E(),
526 pparent->Pz() / pparent->E());
528 fProducts_rho[0].Boost(boostLab_rho);
529 fProducts_rho[1].Boost(boostLab_rho);
533 //-----------------------------------------------------------------------------//
534 // Generate Eta Dalitz decay //
535 //-----------------------------------------------------------------------------//
537 else if(idpart==221){
538 pmass_eta = pparent->M();
541 // Sample the electron pair mass from a histogram
542 epmass_eta = fEPMassEta->GetRandom();
543 if(pmass_eta-omass_gamma>epmass_eta && epmass_eta/2.>emass)
547 // electron pair kinematics in virtual photon rest frame
548 e1_eta = epmass_eta / 2.;
549 p1_eta = TMath::Sqrt((e1_eta + emass) * (e1_eta - emass));
551 // momentum vectors of electrons in virtual photon rest frame
552 Double_t pProd1_eta[3] = {p1_eta * sintheta * cosphi,
553 p1_eta * sintheta * sinphi,
555 Double_t pProd2_eta[3] = {-1.0 * p1_eta * sintheta * cosphi,
556 -1.0 * p1_eta * sintheta * sinphi,
557 -1.0 * p1_eta * costheta};
559 // third child kinematics in parent meson rest frame
560 e3_gamma_eta = (pmass_eta * pmass_eta - epmass_eta * epmass_eta)/(2. * pmass_eta);
561 p3_gamma_eta = TMath::Sqrt((e3_gamma_eta * e3_gamma_eta));
564 // third child 4-vector in parent meson rest frame
565 fProducts_eta[2].SetPx(p3_gamma_eta * sintheta * cosphi);
566 fProducts_eta[2].SetPy(p3_gamma_eta * sintheta * sinphi);
567 fProducts_eta[2].SetPz(p3_gamma_eta * costheta);
568 fProducts_eta[2].SetE(e3_gamma_eta);
570 // electron 4-vectors in properly rotated virtual photon rest frame
571 Double_t pRot1_eta[3] = {0.};
572 Rot(pProd1_eta, pRot1_eta, costheta, -sintheta, -cosphi, -sinphi);
573 Double_t pRot2_eta[3] = {0.};
574 Rot(pProd2_eta, pRot2_eta, costheta, -sintheta, -cosphi, -sinphi);
575 fProducts_eta[0].SetPx(pRot1_eta[0]);
576 fProducts_eta[0].SetPy(pRot1_eta[1]);
577 fProducts_eta[0].SetPz(pRot1_eta[2]);
578 fProducts_eta[0].SetE(e1_eta);
579 fProducts_eta[1].SetPx(pRot2_eta[0]);
580 fProducts_eta[1].SetPy(pRot2_eta[1]);
581 fProducts_eta[1].SetPz(pRot2_eta[2]);
582 fProducts_eta[1].SetE(e1_eta);
584 // boost the dielectron into the parent meson's rest frame
585 Double_t eLPparent_eta = TMath::Sqrt(p3_gamma_eta * p3_gamma_eta + epmass_eta * epmass_eta);
586 TVector3 boostPair_eta( -1.0 * fProducts_eta[2].Px() / eLPparent_eta,
587 -1.0 * fProducts_eta[2].Py() / eLPparent_eta,
588 -1.0 * fProducts_eta[2].Pz() / eLPparent_eta);
589 fProducts_eta[0].Boost(boostPair_eta);
590 fProducts_eta[1].Boost(boostPair_eta);
592 // boost all decay products into the lab frame
593 TVector3 boostLab_eta(pparent->Px() / pparent->E(),
594 pparent->Py() / pparent->E(),
595 pparent->Pz() / pparent->E());
597 fProducts_eta[0].Boost(boostLab_eta);
598 fProducts_eta[1].Boost(boostLab_eta);
599 fProducts_eta[2].Boost(boostLab_eta);
604 //-----------------------------------------------------------------------------//
605 // Generate Omega Dalitz decay //
606 //-----------------------------------------------------------------------------//
608 else if(idpart==223){
609 pmass_omega_dalitz = pparent->M();
611 // Sample the electron pair mass from a histogram
612 epmass_omega_dalitz = fEPMassOmegaDalitz->GetRandom();
613 if(pmass_omega_dalitz-omass_pion>epmass_omega_dalitz && epmass_omega_dalitz/2.>emass)
616 // electron pair kinematics in virtual photon rest frame
617 e1_omega = epmass_omega_dalitz / 2.;
618 p1_omega = TMath::Sqrt((e1_omega + emass) * (e1_omega - emass));
620 // momentum vectors of electrons in virtual photon rest frame
621 Double_t pProd1_omega_dalitz[3] = {p1_omega * sintheta * cosphi,
622 p1_omega * sintheta * sinphi,
623 p1_omega * costheta};
624 Double_t pProd2_omega_dalitz[3] = {-1.0 * p1_omega * sintheta * cosphi,
625 -1.0 * p1_omega * sintheta * sinphi,
626 -1.0 * p1_omega * costheta};
628 // third child kinematics in parent meson rest frame
629 e3_pion = (pmass_omega_dalitz * pmass_omega_dalitz + omass_pion * omass_pion - epmass_omega_dalitz * epmass_omega_dalitz)/(2. * pmass_omega_dalitz);
630 p3_pion = TMath::Sqrt((e3_pion + omass_pion) * (e3_pion - omass_pion));
632 // third child 4-vector in parent meson rest frame
633 fProducts_omega_dalitz[2].SetPx(p3_pion * sintheta * cosphi);
634 fProducts_omega_dalitz[2].SetPy(p3_pion * sintheta * sinphi);
635 fProducts_omega_dalitz[2].SetPz(p3_pion * costheta);
636 fProducts_omega_dalitz[2].SetE(e3_pion);
638 // lepton 4-vectors in properly rotated virtual photon rest frame
639 Double_t pRot1_omega_dalitz[3] = {0.};
640 Rot(pProd1_omega_dalitz, pRot1_omega_dalitz, costheta, -sintheta, -cosphi, -sinphi);
641 Double_t pRot2_omega_dalitz[3] = {0.};
642 Rot(pProd2_omega_dalitz, pRot2_omega_dalitz, costheta, -sintheta, -cosphi, -sinphi);
643 fProducts_omega_dalitz[0].SetPx(pRot1_omega_dalitz[0]);
644 fProducts_omega_dalitz[0].SetPy(pRot1_omega_dalitz[1]);
645 fProducts_omega_dalitz[0].SetPz(pRot1_omega_dalitz[2]);
646 fProducts_omega_dalitz[0].SetE(e1_omega);
647 fProducts_omega_dalitz[1].SetPx(pRot2_omega_dalitz[0]);
648 fProducts_omega_dalitz[1].SetPy(pRot2_omega_dalitz[1]);
649 fProducts_omega_dalitz[1].SetPz(pRot2_omega_dalitz[2]);
650 fProducts_omega_dalitz[1].SetE(e1_omega);
652 // boost the dielectron into the parent meson's rest frame
653 Double_t eLPparent_omega = TMath::Sqrt(p3_pion * p3_pion + epmass_omega_dalitz * epmass_omega_dalitz);
654 TVector3 boostPair_omega( -1.0 * fProducts_omega_dalitz[2].Px() / eLPparent_omega,
655 -1.0 * fProducts_omega_dalitz[2].Py() / eLPparent_omega,
656 -1.0 * fProducts_omega_dalitz[2].Pz() / eLPparent_omega);
657 fProducts_omega_dalitz[0].Boost(boostPair_omega);
658 fProducts_omega_dalitz[1].Boost(boostPair_omega);
660 // boost all decay products into the lab frame
661 TVector3 boostLab_omega_dalitz(pparent->Px() / pparent->E(),
662 pparent->Py() / pparent->E(),
663 pparent->Pz() / pparent->E());
665 fProducts_omega_dalitz[0].Boost(boostLab_omega_dalitz);
666 fProducts_omega_dalitz[1].Boost(boostLab_omega_dalitz);
667 fProducts_omega_dalitz[2].Boost(boostLab_omega_dalitz);
670 //-----------------------------------------------------------------------------//
671 // Generate Omega resonance decay //
672 //-----------------------------------------------------------------------------//
675 // calculate omega mass
676 mp_omega = pparent->M();
679 Double_t x_omega=pparent->Px(); Double_t y_omega=pparent->Py(); Double_t z_omega=pparent->Pz();
680 Double_t t_omega=pparent->E();
681 Double_t p_omega=x_omega*x_omega+y_omega*y_omega+z_omega*z_omega;
682 Double_t Q2_omega= abs((t_omega*t_omega)-(p_omega*p_omega));
683 mp_omega = sqrt(Q2_omega);
687 if ( mp_omega< 2.*md_omega )
689 printf("Omega into ee Decay kinematically impossible! \n");
694 // Sample the electron pair mass from a histogram
695 epmass_omega = fEPMassOmega->GetRandom();
696 if( mp_omega < 2.*epmass_omega ) break;
699 // electron pair kinematics in virtual photon rest frame
700 Ed_omega = epmass_omega/2.;
701 pd_omega = TMath::Sqrt((Ed_omega+md_omega)*(Ed_omega-md_omega));
703 // momentum vectors of electrons in virtual photon rest frame
704 Double_t pProd1_omega[3] = {pd_omega * sintheta * cosphi,
705 pd_omega * sintheta * sinphi,
706 pd_omega * costheta};
708 Double_t pProd2_omega[3] = {-1.0 * pd_omega * sintheta * cosphi,
709 -1.0 * pd_omega * sintheta * sinphi,
710 -1.0 * pd_omega * costheta};
713 // lepton 4 vectors in properly rotated virtual photon rest frame
714 Double_t pRot1_omega[3] = {0.};
715 Rot(pProd1_omega, pRot1_omega, costheta, -sintheta, -cosphi, -sinphi);
716 Double_t pRot2_omega[3] = {0.};
717 Rot(pProd2_omega, pRot2_omega, costheta, -sintheta, -cosphi, -sinphi);
718 fProducts_omega[0].SetPx(pRot1_omega[0]);
719 fProducts_omega[0].SetPy(pRot1_omega[1]);
720 fProducts_omega[0].SetPz(pRot1_omega[2]);
721 fProducts_omega[0].SetE(Ed_omega);
722 fProducts_omega[1].SetPx(pRot2_omega[0]);
723 fProducts_omega[1].SetPy(pRot2_omega[1]);
724 fProducts_omega[1].SetPz(pRot2_omega[2]);
725 fProducts_omega[1].SetE(Ed_omega);
727 // boost decay products into the lab frame
728 TVector3 boostLab_omega(pparent->Px() / pparent->E(),
729 pparent->Py() / pparent->E(),
730 pparent->Pz() / pparent->E());
732 fProducts_omega[0].Boost(boostLab_omega);
733 fProducts_omega[1].Boost(boostLab_omega);
737 //-----------------------------------------------------------------------------//
738 // Generate Etaprime Dalitz decay //
739 //-----------------------------------------------------------------------------//
741 else if(idpart==331){
742 pmass_etaprime = pparent->M();
744 // Sample the electron pair mass from a histogram
745 epmass_etaprime = fEPMassEtaPrime->GetRandom();
746 if(pmass_etaprime-omass_gamma>epmass_etaprime && epmass_etaprime/2.>emass)
749 // electron pair kinematics in virtual photon rest frame
750 e1_etaprime = epmass_etaprime / 2.;
751 p1_etaprime = TMath::Sqrt((e1_etaprime + emass) * (e1_etaprime - emass));
753 // momentum vectors of electrons in virtual photon rest frame
754 Double_t pProd1_etaprime[3] = {p1_etaprime * sintheta * cosphi,
755 p1_etaprime * sintheta * sinphi,
756 p1_etaprime * costheta};
757 Double_t pProd2_etaprime[3] = {-1.0 * p1_etaprime * sintheta * cosphi,
758 -1.0 * p1_etaprime * sintheta * sinphi,
759 -1.0 * p1_etaprime * costheta};
761 // third child kinematics in parent meson rest frame
762 e3_gamma_etaprime = (pmass_etaprime * pmass_etaprime + omass_gamma * omass_gamma - epmass_etaprime * epmass_etaprime)/(2. * pmass_etaprime);
763 p3_gamma_etaprime = TMath::Sqrt((e3_gamma_etaprime + omass_gamma) * (e3_gamma_etaprime - omass_gamma));
765 // third child 4-vector in parent meson rest frame
766 fProducts_etaprime[2].SetPx(p3_gamma_etaprime * sintheta * cosphi);
767 fProducts_etaprime[2].SetPy(p3_gamma_etaprime * sintheta * sinphi);
768 fProducts_etaprime[2].SetPz(p3_gamma_etaprime * costheta);
769 fProducts_etaprime[2].SetE(e3_gamma_etaprime);
771 // electron 4-vectors in properly rotated virtual photon rest frame
772 Double_t pRot1_etaprime[3] = {0.};
773 Rot(pProd1_etaprime, pRot1_etaprime, costheta, -sintheta, -cosphi, -sinphi);
774 Double_t pRot2_etaprime[3] = {0.};
775 Rot(pProd2_etaprime, pRot2_etaprime, costheta, -sintheta, -cosphi, -sinphi);
776 fProducts_etaprime[0].SetPx(pRot1_etaprime[0]);
777 fProducts_etaprime[0].SetPy(pRot1_etaprime[1]);
778 fProducts_etaprime[0].SetPz(pRot1_etaprime[2]);
779 fProducts_etaprime[0].SetE(e1_etaprime);
780 fProducts_etaprime[1].SetPx(pRot2_etaprime[0]);
781 fProducts_etaprime[1].SetPy(pRot2_etaprime[1]);
782 fProducts_etaprime[1].SetPz(pRot2_etaprime[2]);
783 fProducts_etaprime[1].SetE(e1_etaprime);
785 // boost the dielectron into the parent meson's rest frame
786 Double_t eLPparent_etaprime = TMath::Sqrt(p3_gamma_etaprime * p3_gamma_etaprime + epmass_etaprime * epmass_etaprime);
787 TVector3 boostPair_etaprime( -1.0 * fProducts_etaprime[2].Px() / eLPparent_etaprime,
788 -1.0 * fProducts_etaprime[2].Py() / eLPparent_etaprime,
789 -1.0 * fProducts_etaprime[2].Pz() / eLPparent_etaprime);
790 fProducts_etaprime[0].Boost(boostPair_etaprime);
791 fProducts_etaprime[1].Boost(boostPair_etaprime);
793 // boost all decay products into the lab frame
794 TVector3 boostLab_etaprime(pparent->Px() / pparent->E(),
795 pparent->Py() / pparent->E(),
796 pparent->Pz() / pparent->E());
798 fProducts_etaprime[0].Boost(boostLab_etaprime);
799 fProducts_etaprime[1].Boost(boostLab_etaprime);
800 fProducts_etaprime[2].Boost(boostLab_etaprime);
804 //-----------------------------------------------------------------------------//
805 // Generate Phi Dalitz decay //
806 //-----------------------------------------------------------------------------//
808 else if(idpart==333){
809 pmass_phi_dalitz = pparent->M();
811 // Sample the electron pair mass from a histogram
812 epmass_phi_dalitz = fEPMassPhiDalitz->GetRandom();
813 if(pmass_phi_dalitz-omass_eta>epmass_phi_dalitz && epmass_phi_dalitz/2.>emass)
816 // electron pair kinematics in virtual photon rest frame
817 e1_phi = epmass_phi_dalitz / 2.;
818 p1_phi = TMath::Sqrt((e1_phi + emass) * (e1_phi - emass));
820 // momentum vectors of electrons in virtual photon rest frame
821 Double_t pProd1_phi_dalitz[3] = {p1_phi * sintheta * cosphi,
822 p1_phi * sintheta * sinphi,
824 Double_t pProd2_phi_dalitz[3] = {-1.0 * p1_phi * sintheta * cosphi,
825 -1.0 * p1_phi * sintheta * sinphi,
826 -1.0 * p1_phi * costheta};
828 // third child kinematics in parent meson rest frame
829 e3_eta = (pmass_phi_dalitz * pmass_phi_dalitz + omass_eta * omass_eta - epmass_phi_dalitz * epmass_phi_dalitz)/(2. * pmass_phi_dalitz);
830 p3_eta = TMath::Sqrt((e3_eta + omass_eta) * (e3_eta - omass_eta));
832 // third child 4-vector in parent meson rest frame
833 fProducts_phi_dalitz[2].SetPx(p3_eta * sintheta * cosphi);
834 fProducts_phi_dalitz[2].SetPy(p3_eta * sintheta * sinphi);
835 fProducts_phi_dalitz[2].SetPz(p3_eta * costheta);
836 fProducts_phi_dalitz[2].SetE(e3_eta);
838 // electron 4-vectors in properly rotated virtual photon rest frame
839 Double_t pRot1_phi_dalitz[3] = {0.};
840 Rot(pProd1_phi_dalitz, pRot1_phi_dalitz, costheta, -sintheta, -cosphi, -sinphi);
841 Double_t pRot2_phi_dalitz[3] = {0.};
842 Rot(pProd2_phi_dalitz, pRot2_phi_dalitz, costheta, -sintheta, -cosphi, -sinphi);
843 fProducts_phi_dalitz[0].SetPx(pRot1_phi_dalitz[0]);
844 fProducts_phi_dalitz[0].SetPy(pRot1_phi_dalitz[1]);
845 fProducts_phi_dalitz[0].SetPz(pRot1_phi_dalitz[2]);
846 fProducts_phi_dalitz[0].SetE(e1_phi);
847 fProducts_phi_dalitz[1].SetPx(pRot2_phi_dalitz[0]);
848 fProducts_phi_dalitz[1].SetPy(pRot2_phi_dalitz[1]);
849 fProducts_phi_dalitz[1].SetPz(pRot2_phi_dalitz[2]);
850 fProducts_phi_dalitz[1].SetE(e1_phi);
852 // boost the dielectron into the parent meson's rest frame
853 Double_t eLPparent_phi = TMath::Sqrt(p3_eta * p3_eta + epmass_phi_dalitz * epmass_phi_dalitz);
854 TVector3 boostPair_phi( -1.0 * fProducts_phi_dalitz[2].Px() / eLPparent_phi,
855 -1.0 * fProducts_phi_dalitz[2].Py() / eLPparent_phi,
856 -1.0 * fProducts_phi_dalitz[2].Pz() / eLPparent_phi);
857 fProducts_phi_dalitz[0].Boost(boostPair_phi);
858 fProducts_phi_dalitz[1].Boost(boostPair_phi);
860 // boost all decay products into the lab frame
861 TVector3 boostLab_phi_dalitz(pparent->Px() / pparent->E(),
862 pparent->Py() / pparent->E(),
863 pparent->Pz() / pparent->E());
865 fProducts_phi_dalitz[0].Boost(boostLab_phi_dalitz);
866 fProducts_phi_dalitz[1].Boost(boostLab_phi_dalitz);
867 fProducts_phi_dalitz[2].Boost(boostLab_phi_dalitz);
870 //-----------------------------------------------------------------------------//
871 // Generate Phi resonance decay //
872 //-----------------------------------------------------------------------------//
875 // calculate phi mass
876 mp_phi = pparent->M();
879 Double_t x_phi=pparent->Px(); Double_t y_phi=pparent->Py(); Double_t z_phi=pparent->Pz();
880 Double_t t_phi=pparent->E();
881 Double_t p_phi=x_phi*x_phi+y_phi*y_phi+z_phi*z_phi;
882 Double_t Q2_phi= abs((t_phi*t_phi)-(p_phi*p_phi));
883 mp_phi = sqrt(Q2_phi);
886 if ( mp_phi< 2.*md_phi )
888 printf("Phi into ee Decay kinematically impossible! \n");
893 // Sample the electron pair mass from a histogram
894 epmass_phi = fEPMassPhi->GetRandom();
895 if(mp_phi < 2.*epmass_phi) break;
898 // electron pair kinematics in virtual photon rest frame
899 Ed_phi = epmass_phi/2.;
900 pd_phi = TMath::Sqrt((Ed_phi+md_phi)*(Ed_phi-md_phi));
902 // momentum vectors of electrons in virtual photon rest frame
903 Double_t pProd1_phi[3] = {pd_phi * sintheta * cosphi,
904 pd_phi * sintheta * sinphi,
906 Double_t pProd2_phi[3] = {-1.0 * pd_phi * sintheta * cosphi,
907 -1.0 * pd_phi * sintheta * sinphi,
908 -1.0 * pd_phi * costheta};
910 // electron 4 vectors in properly rotated virtual photon rest frame
911 Double_t pRot1_phi[3] = {0.};
912 Rot(pProd1_phi, pRot1_phi, costheta, -sintheta, -cosphi, -sinphi);
913 Double_t pRot2_phi[3] = {0.};
914 Rot(pProd2_phi, pRot2_phi, costheta, -sintheta, -cosphi, -sinphi);
915 fProducts_phi[0].SetPx(pRot1_phi[0]);
916 fProducts_phi[0].SetPy(pRot1_phi[1]);
917 fProducts_phi[0].SetPz(pRot1_phi[2]);
918 fProducts_phi[0].SetE(Ed_phi);
919 fProducts_phi[1].SetPx(pRot2_phi[0]);
920 fProducts_phi[1].SetPy(pRot2_phi[1]);
921 fProducts_phi[1].SetPz(pRot2_phi[2]);
922 fProducts_phi[1].SetE(Ed_phi);
924 // boost decay products into the lab frame
925 TVector3 boostLab_phi(pparent->Px() / pparent->E(),
926 pparent->Py() / pparent->E(),
927 pparent->Pz() / pparent->E());
929 fProducts_phi[0].Boost(boostLab_phi);
930 fProducts_phi[1].Boost(boostLab_phi);
934 //-----------------------------------------------------------------------------//
935 // Generate Jpsi resonance decay //
936 //-----------------------------------------------------------------------------//
938 else if(idpart==443){
939 // calculate jpsi mass
941 mp_jpsi = pparent->M();
944 /*Double_t x_jpsi=pparent->Px();
945 Double_t y_jpsi=pparent->Py();
946 Double_t z_jpsi=pparent->Pz();
947 Double_t t_jpsi=pparent->E();
948 Double_t p_jpsi=x_jpsi*x_jpsi+y_jpsi*y_jpsi+z_jpsi*z_jpsi;
949 Double_t Q2_jpsi= abs((t_jpsi*t_jpsi)-(p_jpsi*p_jpsi));
950 mp_jpsi = sqrt(Q2_jpsi);*/
957 if ( mp_jpsi < 2.*md_jpsi )
959 printf("JPsi into ee Decay kinematically impossible! \n");
964 // Sample the electron pair mass from a histogram
965 epmass_jpsi = fEPMassJPsi->GetRandom();
966 if ( mp_jpsi < 2.*epmass_jpsi ) break;
968 // electron pair kinematics in virtual photon rest frame
969 Ed_jpsi = epmass_jpsi/2.;
970 pd_jpsi = TMath::Sqrt((Ed_jpsi+md_jpsi)*(Ed_jpsi-md_jpsi));
972 // momentum vectors of electrons in virtual photon rest frame
973 Double_t pProd1_jpsi[3] = {pd_jpsi * sintheta * cosphi,
974 pd_jpsi * sintheta * sinphi,
977 Double_t pProd2_jpsi[3] = {-1.0 * pd_jpsi * sintheta * cosphi,
978 -1.0 * pd_jpsi * sintheta * sinphi,
979 -1.0 * pd_jpsi * costheta};
982 // electron 4 vectors in properly rotated virtual photon rest frame
983 Double_t pRot1_jpsi[3] = {0.};
984 Rot(pProd1_jpsi, pRot1_jpsi, costheta, -sintheta, -cosphi, -sinphi);
985 Double_t pRot2_jpsi[3] = {0.};
986 Rot(pProd2_jpsi, pRot2_jpsi, costheta, -sintheta, -cosphi, -sinphi);
987 fProducts_jpsi[0].SetPx(pRot1_jpsi[0]);
988 fProducts_jpsi[0].SetPy(pRot1_jpsi[1]);
989 fProducts_jpsi[0].SetPz(pRot1_jpsi[2]);
990 fProducts_jpsi[0].SetE(Ed_jpsi);
991 fProducts_jpsi[1].SetPx(pRot2_jpsi[0]);
992 fProducts_jpsi[1].SetPy(pRot2_jpsi[1]);
993 fProducts_jpsi[1].SetPz(pRot2_jpsi[2]);
994 fProducts_jpsi[1].SetE(Ed_jpsi);
997 // boost decay products into the lab frame
998 TVector3 boostLab_jpsi(pparent->Px() / pparent->E(),
999 pparent->Py() / pparent->E(),
1000 pparent->Pz() / pparent->E());
1002 fProducts_jpsi[0].Boost(boostLab_jpsi);
1003 fProducts_jpsi[1].Boost(boostLab_jpsi);
1010 void AliDecayerExodus::Rot(Double_t pin[3], Double_t pout[3], Double_t costheta, Double_t sintheta,
1011 Double_t cosphi, Double_t sinphi) const
1014 pout[0] = pin[0]*costheta*cosphi-pin[1]*sinphi+pin[2]*sintheta*cosphi;
1015 pout[1] = pin[0]*costheta*sinphi+pin[1]*cosphi+pin[2]*sintheta*sinphi;
1016 pout[2] = -1.0 * pin[0] * sintheta + pin[2] * costheta;
1021 Int_t AliDecayerExodus::ImportParticles(TClonesArray *particles)
1024 // Import particles for Dalitz and resonance decays
1027 TClonesArray &clonesParticles = *particles;
1030 Double_t px, py, pz, e;
1032 Int_t pdgD [3][3] = { {kElectron, -kElectron, 22}, // pizero, eta, etaprime
1033 {kElectron, -kElectron, 111}, // omega dalitz
1034 {kElectron, -kElectron, 221} }; // phi dalitz
1036 Int_t pdgR [2] = {kElectron, -kElectron}; // rho, omega, phi, jpsi
1040 Int_t parentD[3] = { 0, 0, -1};
1041 Int_t dauD1 [3] = {-1, -1, 1};
1042 Int_t dauD2 [3] = {-1, -1, 2};
1044 Int_t parentR[2] = { 0, 0};
1045 Int_t dauR1 [2] = { -1, -1};
1046 Int_t dauR2 [2] = { -1, -1};
1048 for (Int_t j = 0; j < 9; j++){
1052 for (i = 2; i > -1; i--) {
1053 px = fProducts_pion[i].Px();
1054 py = fProducts_pion[i].Py();
1055 pz = fProducts_pion[i].Pz();
1056 e = fProducts_pion[i].E();
1057 new(clonesParticles[2 - i]) TParticle(pdgD[0][i], 1, parentD[i], -1, dauD1[i], dauD2[i], px, py, pz, e, 0., 0., 0., 0.);
1064 for (k = 1; k > -1; k--) {
1065 px = fProducts_rho[k].Px();
1066 py = fProducts_rho[k].Py();
1067 pz = fProducts_rho[k].Pz();
1068 e = fProducts_rho[k].E();
1069 new(clonesParticles[1 - k]) TParticle(pdgR[k], 1, parentR[k], -1, dauR1[k], dauR2[k], px, py, pz, e, 0., 0., 0., 0.);
1076 for (i = 2; i > -1; i--) {
1077 px = fProducts_eta[i].Px();
1078 py = fProducts_eta[i].Py();
1079 pz = fProducts_eta[i].Pz();
1080 e = fProducts_eta[i].E();
1081 new(clonesParticles[2 - i]) TParticle(pdgD[0][i], 1, parentD[i], -1, dauD1[i], dauD2[i], px, py, pz, e, 0., 0., 0., 0.);
1088 for (i = 2; i > -1; i--) {
1089 px = fProducts_omega_dalitz[i].Px();
1090 py = fProducts_omega_dalitz[i].Py();
1091 pz = fProducts_omega_dalitz[i].Pz();
1092 e = fProducts_omega_dalitz[i].E();
1093 new(clonesParticles[2 - i]) TParticle(pdgD[1][i], 1, parentD[i], -1, dauD1[i], dauD2[i], px, py, pz, e, 0., 0., 0., 0.);
1100 for (k = 1; k > -1; k--) {
1101 px = fProducts_rho[k].Px();
1102 py = fProducts_rho[k].Py();
1103 pz = fProducts_rho[k].Pz();
1104 e = fProducts_rho[k].E();
1105 new(clonesParticles[1 - k]) TParticle(pdgR[k], 1, parentR[k], -1, dauR1[k], dauR2[k], px, py, pz, e, 0., 0., 0., 0.);
1112 for (i = 2; i > -1; i--) {
1113 px = fProducts_etaprime[i].Px();
1114 py = fProducts_etaprime[i].Py();
1115 pz = fProducts_etaprime[i].Pz();
1116 e = fProducts_etaprime[i].E();
1117 new(clonesParticles[2 - i]) TParticle(pdgD[0][i], 1, parentD[i], -1, dauD1[i], dauD2[i], px, py, pz, e, 0., 0., 0., 0.);
1124 for (i = 2; i > -1; i--) {
1125 px = fProducts_phi_dalitz[i].Px();
1126 py = fProducts_phi_dalitz[i].Py();
1127 pz = fProducts_phi_dalitz[i].Pz();
1128 e = fProducts_phi_dalitz[i].E();
1129 new(clonesParticles[2 - i]) TParticle(pdgD[2][i], 1, parentD[i], -1, dauD1[i], dauD2[i], px, py, pz, e, 0., 0., 0., 0.);
1137 for (k = 1; k > -1; k--) {
1138 px = fProducts_phi[k].Px();
1139 py = fProducts_phi[k].Py();
1140 pz = fProducts_phi[k].Pz();
1141 e = fProducts_phi[k].E();
1142 new(clonesParticles[1 - k]) TParticle(pdgR[k], 1, parentR[k], -1, dauR1[k], dauR2[k], px, py, pz, e, 0., 0., 0., 0.);
1149 for (k = 1; k > -1; k--) {
1150 px = fProducts_jpsi[k].Px();
1151 py = fProducts_jpsi[k].Py();
1152 pz = fProducts_jpsi[k].Pz();
1153 e = fProducts_jpsi[k].E();
1154 new(clonesParticles[1 - k]) TParticle(pdgR[k], 1, parentR[k], -1, dauR1[k], dauR2[k], px, py, pz, e, 0., 0., 0., 0.);
1161 return particles->GetEntries();
1166 void AliDecayerExodus::Decay(TClonesArray *array)
1168 // Replace all Dalitz(pi0,eta,omega,eta',phi) and resonance(rho,omega,phi,jpsi) decays with the correct matrix element decays
1169 // for di-electron cocktail calculations
1172 Int_t nt = array->GetEntriesFast();
1175 Int_t fd3, ld3, fd2, ld2, fd, ld;
1178 for (Int_t i = 0; i < nt; i++) {
1179 TParticle* part = (TParticle*) (array->At(i));
1180 if (part->GetPdgCode() != 111 || part->GetPdgCode() != 221 || part->GetPdgCode() != 223 || part->GetPdgCode() != 331 || part->GetPdgCode() != 333 || part->GetPdgCode() != 443 ) continue;
1185 if(part->GetPdgCode() == 111){
1187 fd3 = part->GetFirstDaughter() - 1;
1188 ld3 = part->GetLastDaughter() - 1;
1190 if (fd3 < 0) continue;
1191 if ((ld3 - fd3) != 2) continue;
1193 for (j = 0; j < 3; j++) dp3[j] = (TParticle*) (array->At(fd3+j));
1195 if((dp3[0]->GetPdgCode() != 22) && (TMath::Abs(dp3[1]->GetPdgCode()) != 11)) continue;
1197 TLorentzVector Pizero(part->Px(), part->Py(), part->Pz(), part->Energy());
1198 Decay(111, &Pizero);
1199 for (j = 0; j < 3; j++) dp3[j]->SetMomentum(fProducts_pion[2-j]);
1207 if(part->GetPdgCode() == 221){
1209 fd3 = part->GetFirstDaughter() - 1;
1210 ld3 = part->GetLastDaughter() - 1;
1212 if (fd3 < 0) continue;
1213 if ((ld3 - fd3) != 2) continue;
1215 for (j = 0; j < 3; j++) dp3[j] = (TParticle*) (array->At(fd3+j));
1217 if((dp3[0]->GetPdgCode() != 22) && ((TMath::Abs(dp3[1]->GetPdgCode()) != 11))) continue;
1219 TLorentzVector Eta(part->Px(), part->Py(), part->Pz(), part->Energy());
1221 for (j = 0; j < 3; j++) dp3[j]->SetMomentum(fProducts_eta[2-j]);
1228 if(part->GetPdgCode() == 113){
1230 fd2 = part->GetFirstDaughter() - 1;
1231 ld2 = part->GetLastDaughter() - 1;
1233 if (fd2 < 0) continue;
1234 if ((ld2 - fd2) != 1) continue;
1236 for (k = 0; k < 2; k++) dp2[k] = (TParticle*) (array->At(fd2+k));
1238 if((dp2[0]->GetPdgCode() != 11) && ((TMath::Abs(dp2[1]->GetPdgCode()) != 11))) continue;
1240 TLorentzVector Rho(part->Px(), part->Py(), part->Pz(), part->Energy());
1242 for (k = 0; k < 2; k++) dp2[k]->SetMomentum(fProducts_rho[1-k]);
1246 // Omega dalitz and direct
1249 if(part->GetPdgCode() == 223){
1251 fd = part->GetFirstDaughter() - 1;
1252 ld = part->GetLastDaughter() - 1;
1254 if (fd < 0) continue;
1256 if ((ld - fd) == 2){
1258 for (j = 0; j < 3; j++) dp3[j] = (TParticle*) (array->At(fd+j));
1259 if( dp3[0]->GetPdgCode() != 111 && (TMath::Abs(dp3[1]->GetPdgCode()) != 11)) continue;
1261 TLorentzVector Omegadalitz(part->Px(), part->Py(), part->Pz(), part->Energy());
1262 Decay(223, &Omegadalitz);
1263 for (j = 0; j < 3; j++) dp3[j]->SetMomentum(fProducts_omega_dalitz[2-j]);
1266 else if ((ld - fd) == 1) {
1268 for (k = 0; k < 2; k++) dp2[k] = (TParticle*) (array->At(fd+k));
1269 if( dp2[0]->GetPdgCode() != 11 && (TMath::Abs(dp2[1]->GetPdgCode()) != 11)) continue;
1271 TLorentzVector Omega(part->Px(), part->Py(), part->Pz(), part->Energy());
1273 for (k = 0; k < 2; k++) dp2[k]->SetMomentum(fProducts_omega[1-k]);
1281 if(part->GetPdgCode() == 331){
1283 fd3 = part->GetFirstDaughter() - 1;
1284 ld3 = part->GetLastDaughter() - 1;
1286 if (fd3 < 0) continue;
1287 if ((ld3 - fd3) != 2) continue;
1289 for (j = 0; j < 3; j++) dp3[j] = (TParticle*) (array->At(fd3+j));
1291 if((dp3[0]->GetPdgCode() != 22) && ((TMath::Abs(dp3[1]->GetPdgCode()) != 11))) continue;
1293 TLorentzVector Etaprime(part->Px(), part->Py(), part->Pz(), part->Energy());
1294 Decay(331, &Etaprime);
1295 for (j = 0; j < 3; j++) dp3[j]->SetMomentum(fProducts_etaprime[2-j]);
1299 // Phi dalitz and direct
1301 if(part->GetPdgCode() == 333){
1303 fd = part->GetFirstDaughter() - 1;
1304 ld = part->GetLastDaughter() - 1;
1306 if (fd < 0) continue;
1307 if ((ld - fd) == 2){
1308 for (j = 0; j < 3; j++) dp3[j] = (TParticle*) (array->At(fd+j));
1309 if( dp3[0]->GetPdgCode() != 221 && (TMath::Abs(dp3[1]->GetPdgCode()) != 11)) continue;
1311 TLorentzVector Phidalitz(part->Px(), part->Py(), part->Pz(), part->Energy());
1312 Decay(333, &Phidalitz);
1313 for (j = 0; j < 3; j++) dp3[j]->SetMomentum(fProducts_phi_dalitz[2-j]);
1316 else if ((ld - fd) == 1) {
1317 for (k = 0; k < 2; k++) dp2[k] = (TParticle*) (array->At(fd+k));
1318 if( dp2[0]->GetPdgCode() != 11 && (TMath::Abs(dp2[1]->GetPdgCode()) != 11)) continue;
1320 TLorentzVector Phi(part->Px(), part->Py(), part->Pz(), part->Energy());
1322 for (k = 0; k < 2; k++) dp2[k]->SetMomentum(fProducts_phi[1-k]);
1330 if(part->GetPdgCode() == 443){
1332 fd2 = part->GetFirstDaughter() - 1;
1333 ld2 = part->GetLastDaughter() - 1;
1335 if (fd2 < 0) continue;
1336 if ((ld2 - fd2) != 1) continue;
1338 for (k = 0; k < 2; k++) dp2[k] = (TParticle*) (array->At(fd2+k));
1340 if((dp2[0]->GetPdgCode() != 11) && ((TMath::Abs(dp2[1]->GetPdgCode()) != 11))) continue;
1342 TLorentzVector JPsi(part->Px(), part->Py(), part->Pz(), part->Energy());
1344 for (k = 0; k < 2; k++) dp2[k]->SetMomentum(fProducts_jpsi[1-k]);
1351 AliDecayerExodus& AliDecayerExodus::operator=(const AliDecayerExodus& rhs)
1353 // Assignment operator
1358 void AliDecayerExodus::Copy(TObject&) const
1363 Fatal("Copy","Not implemented!\n");
1367 AliDecayerExodus::AliDecayerExodus(const AliDecayerExodus &decayer)
1374 fEPMassOmegaDalitz(0),
1376 fEPMassPhiDalitz(0),
1381 decayer.Copy(*this);