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