small update
[u/mrichter/AliRoot.git] / EVGEN / AliDecayerExodus.cxx
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
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
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;
95     mass_min = 0.;
96     mass_max = 0.;
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); 
296      weight_jpsi    = (Float_t)Lorentz(mass_bin,vmass_jpsi,vwidth_jpsi);
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
310 // of rho, omega and phi
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;
318  
319   if(mass>pimass){
320   corr = vwidth*(vmass/mass)*exp(1.5*log((mass*mass/4.0-pimass*pimass)
321          /(vmass*vmass/4.0-pimass*pimass)));
322   }
323
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
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
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{
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;
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