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