EVGEN/AliOmegaDalitz.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 /* $Id$ */
  18
  19
  20 #include "AliOmegaDalitz.h"
  21 #include <TMath.h>
  22 #include <AliLog.h>
  23 #include <TH1.h>
  24 #include <TPDGCode.h>
  25 #include <TRandom.h>
  26 #include <TParticle.h>
  27 #include <TDatabasePDG.h>
  28 #include <TLorentzVector.h>
  29 #include <TClonesArray.h>
  30
  31 ClassImp(AliOmegaDalitz)
  32
  33
  34 //-----------------------------------------------------------------------------
  35 //
  36 // Generate lepton-pair mass distributions for Dalitz decays according
  37 // to the Kroll-Wada parametrization: N. Kroll, W. Wada: Phys. Rev 98(1955)1355
  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 //-----------------------------------------------------------------------------
  43 //
  44 // Adaption for AliRoot
  45 //
  46 // R. Averbeck
  47 // A. Morsch
  48 //
  49 AliOmegaDalitz::AliOmegaDalitz():
  50         AliDecayer(),
  51         fEPMass(0),
  52         fMPMass(0),
  53         fInit(0)
  54 {
  55     // Constructor
  56 }
  57
  58 void AliOmegaDalitz::Init()
  59 {
  60   // Initialisation
  61   Int_t    idecay, ibin, nbins = 1000;
  62   Double_t min, max, binwidth;
  63   Double_t pmass, lmass, omass, emass, mmass;
  64   Double_t epsilon, delta, mLL, q, kwHelp, krollWada, formFactor, weight;
  65
  66   // Get the particle masses
  67   // electron
  68   emass = (TDatabasePDG::Instance()->GetParticle(kElectron))->Mass();
  69   // muon
  70   mmass = (TDatabasePDG::Instance()->GetParticle(kMuonPlus))->Mass();
  71   // omega
  72   pmass = (TDatabasePDG::Instance()->GetParticle(223))      ->Mass();
  73   // pi0
  74   omass = (TDatabasePDG::Instance()->GetParticle(kPi0))     ->Mass();
  75
  76   for ( idecay = 1; idecay<=2; idecay++ )
  77   {
  78     if ( idecay == 1 )
  79       lmass = emass;
  80     else
  81       lmass = mmass;
  82
  83     min = 2.0 * lmass;
  84     max = pmass - omass;
  85     binwidth = (max - min) / (Double_t)nbins;
  86     if ( idecay == 1 )
  87       fEPMass = new TH1F("fEPMass", "Dalitz electron pair", nbins, min, max);
  88     else
  89       fMPMass = new TH1F("fMPMass", "Dalitz muon pair", nbins, min, max);
  90
  91     epsilon = (lmass / pmass) * (lmass / pmass);
  92     delta   = (omass / pmass) * (omass / pmass);
  93
  94     for ( ibin = 1; ibin <= nbins; ibin++ )
  95     {
  96       mLL = min + (Double_t)(ibin - 1) * binwidth + binwidth / 2.0;
  97       q    = (mLL / pmass) * (mLL / pmass);
  98       if ( q <= 4.0 * epsilon )
  99       {
 100         AliFatal("Error in calculating Dalitz mass histogram binning!");
 101       }
 102
 103       kwHelp = (1.0 + q /  (1.0 - delta)) * (1.0 + q / (1.0 - delta))
 104               - 4.0 * q / ((1.0 - delta) * (1.0 - delta));
 105       if ( kwHelp <= 0.0 )
 106       {
 107         AliFatal("Error in calculating Dalitz mass histogram binning!");
 108       }
 109       krollWada = (2.0 / mLL) * TMath::Exp(1.5 * TMath::Log(kwHelp))
 110                               * TMath::Sqrt(1.0 - 4.0 * epsilon / q)
 111                               * (1.0 + 2.0 * epsilon / q);
 112       formFactor =
 113         (TMath::Power(TMath::Power(0.6519,2),2))
 114         / (TMath::Power(TMath::Power(0.6519,2)-TMath::Power(mLL, 2), 2)
 115            + TMath::Power(0.04198, 2)*TMath::Power(0.6519, 2));
 116       weight = krollWada * formFactor;
 117       if ( idecay == 1 )
 118         fEPMass->AddBinContent(ibin, weight);
 119       else
 120         fMPMass->AddBinContent(ibin, weight);
 121     }
 122   }
 123 }
 124
 125
 126 void AliOmegaDalitz::Decay(Int_t idlepton, TLorentzVector* pparent)
 127 {
 128 //-----------------------------------------------------------------------------
 129 //
 130 //  Generate Omega Dalitz decay
 131 //
 132 //-----------------------------------------------------------------------------
 133
 134     if (!fInit) {
 135         Init();
 136         fInit=1;
 137     }
 138
 139   Double_t pmass, lmass, omass, lpmass;
 140   Double_t e1, p1, e3, p3;
 141   Double_t betaSquare, lambda;
 142   Double_t costheta, sintheta, cosphi, sinphi, phi;
 143
 144   // Get the particle masses
 145   // lepton
 146   lmass = (TDatabasePDG::Instance()->GetParticle(idlepton))->Mass();
 147   // omega
 148   pmass = pparent->M();
 149   // pi0
 150   omass = (TDatabasePDG::Instance()->GetParticle(kPi0))     ->Mass();
 151
 152   // Sample the lepton pair mass from a histogram
 153   for( ;; )
 154   {
 155     if ( idlepton == kElectron )
 156       lpmass = fEPMass->GetRandom();
 157     else
 158       lpmass = fMPMass->GetRandom();
 159     if ( pmass - omass > lpmass && lpmass / 2. > lmass ) break;
 160   }
 161
 162   // lepton pair kinematics in virtual photon rest frame
 163   e1 = lpmass / 2.;
 164   p1 = TMath::Sqrt((e1 + lmass) * (e1 - lmass));
 165   betaSquare = 1.0 - 4.0 * (lmass * lmass) / (lpmass * lpmass);
 166   lambda      = betaSquare / (2.0 - betaSquare);
 167   costheta = (2.0 * gRandom->Rndm()) - 1.;
 168   sintheta = TMath::Sqrt((1. + costheta) * (1. - costheta));
 169   phi      = 2.0 * TMath::ACos(-1.) * gRandom->Rndm();
 170   sinphi   = TMath::Sin(phi);
 171   cosphi   = TMath::Cos(phi);
 172   // momentum vectors of leptons in virtual photon rest frame
 173   Double_t pProd1[3] = {p1 * sintheta * cosphi,
 174                        p1 * sintheta * sinphi,
 175                        p1 * costheta};
 176   Double_t pProd2[3] = {-1.0 * p1 * sintheta * cosphi,
 177                        -1.0 * p1 * sintheta * sinphi,
 178                        -1.0 * p1 * costheta};
 179
 180   // pizero kinematics in omega rest frame
 181   e3       = (pmass * pmass + omass * omass - lpmass * lpmass)/(2. * pmass);
 182   p3       = TMath::Sqrt((e3 + omass)  * (e3 - omass));
 183   costheta = (2.0 * gRandom->Rndm()) - 1.;
 184   sintheta = TMath::Sqrt((1. + costheta) * (1. - costheta));
 185   phi      = 2.0 * TMath::ACos(-1.) * gRandom->Rndm();
 186   sinphi   = TMath::Sin(phi);
 187   cosphi   = TMath::Cos(phi);
 188   // pizero 4-vector in omega rest frame
 189   fProducts[2].SetPx(p3 * sintheta * cosphi);
 190   fProducts[2].SetPy(p3 * sintheta * sinphi);
 191   fProducts[2].SetPz(p3 * costheta);
 192   fProducts[2].SetE(e3);
 193
 194   // lepton 4-vectors in properly rotated virtual photon rest frame
 195   Double_t pRot1[3] = {0.};
 196   Rot(pProd1, pRot1, costheta, -sintheta, -cosphi, -sinphi);
 197   Double_t pRot2[3] = {0.};
 198   Rot(pProd2, pRot2, costheta, -sintheta, -cosphi, -sinphi);
 199   fProducts[0].SetPx(pRot1[0]);
 200   fProducts[0].SetPy(pRot1[1]);
 201   fProducts[0].SetPz(pRot1[2]);
 202   fProducts[0].SetE(e1);
 203   fProducts[1].SetPx(pRot2[0]);
 204   fProducts[1].SetPy(pRot2[1]);
 205   fProducts[1].SetPz(pRot2[2]);
 206   fProducts[1].SetE(e1);
 207
 208   // boost the dilepton into the omega's rest frame
 209   Double_t eLPparent = TMath::Sqrt(p3 * p3 + lpmass * lpmass);
 210   TVector3 boostPair( -1.0 * fProducts[2].Px() / eLPparent,
 211                       -1.0 * fProducts[2].Py() / eLPparent,
 212                       -1.0 * fProducts[2].Pz() / eLPparent);
 213   fProducts[0].Boost(boostPair);
 214   fProducts[1].Boost(boostPair);
 215
 216   // boost all decay products into the lab frame
 217   TVector3 boostLab( pparent->Px() / pparent->E(),
 218                      pparent->Py() / pparent->E(),
 219                      pparent->Pz() / pparent->E());
 220   fProducts[0].Boost(boostLab);
 221   fProducts[1].Boost(boostLab);
 222   fProducts[2].Boost(boostLab);
 223
 224   return;
 225
 226 }
 227
 228 Int_t AliOmegaDalitz::ImportParticles(TClonesArray *particles)
 229 {
 230     // Import TParticles in array particles
 231     // l+ l- pi0
 232
 233     TClonesArray &clonesParticles = *particles;
 234
 235     Int_t pdg   [3] = {kElectron, -kElectron, kPi0};
 236     if ( fProducts[1].M() > 0.001 )
 237     {
 238       pdg[0] =  kMuonPlus;
 239       pdg[1] = -kMuonPlus;
 240     }
 241     Int_t parent[3] = {0, 0, -1};
 242     Int_t d1    [3] = {-1, -1, 1};
 243     Int_t d2    [3] = {-1, -1, 2};
 244     for (Int_t i = 2; i > -1; i--) {
 245         Double_t px = fProducts[i].Px();
 246         Double_t py = fProducts[i].Py();
 247         Double_t pz = fProducts[i].Pz();
 248         Double_t e  = fProducts[i].E();
 249         //
 250         new(clonesParticles[2 - i]) TParticle(pdg[i], 1, parent[i], -1, d1[i], d2[i], px, py, pz, e, 0., 0., 0., 0.);
 251     }
 252     return (3);
 253 }
 254
 255
 256 void AliOmegaDalitz::
 257 Rot(Double_t pin[3], Double_t pout[3], Double_t costheta, Double_t sintheta,
 258     Double_t cosphi, Double_t sinphi) const
 259 {
 260 // Perform rotation
 261   pout[0] = pin[0]*costheta*cosphi-pin[1]*sinphi+pin[2]*sintheta*cosphi;
 262   pout[1] = pin[0]*costheta*sinphi+pin[1]*cosphi+pin[2]*sintheta*sinphi;
 263   pout[2] = -1.0  * pin[0] * sintheta + pin[2] * costheta;
 264   return;
 265 }
 266
 267 void AliOmegaDalitz::Decay(TClonesArray* array)
 268 {
 269   //
 270   // Replace all omega dalitz decays with the correct matrix element decays
 271   //
 272   printf("-->Correcting Dalitz \n");
 273   Int_t nt = array->GetEntriesFast();
 274   TParticle* dp[3];
 275   for (Int_t i = 0; i < nt; i++) {
 276     TParticle* part = (TParticle*) (array->At(i));
 277     if (part->GetPdgCode() != 223)                     continue;
 278
 279     Int_t fd = part->GetFirstDaughter() - 1;
 280     Int_t ld = part->GetLastDaughter()  - 1;
 281
 282     if (fd < 0)                                        continue;
 283     if ((ld - fd) != 2)                                continue;
 284
 285     for (Int_t j = 0; j < 3; j++) dp[j] = (TParticle*) (array->At(fd+j));
 286     if ((dp[0]->GetPdgCode() != kPi0) ||
 287         ((TMath::Abs(dp[1]->GetPdgCode()) != kElectron) &&
 288          (TMath::Abs(dp[1]->GetPdgCode()) != kMuonPlus))) continue;
 289     TLorentzVector omega(part->Px(), part->Py(), part->Pz(), part->Energy());
 290     if ( TMath::Abs(dp[1]->GetPdgCode()) == kElectron )
 291       Decay(kElectron, &omega);
 292     else
 293       Decay(kMuonPlus, &omega);
 294     for (Int_t j = 0; j < 3; j++) dp[j]->SetMomentum(fProducts[2-j]);
 295     printf("original %13.3f %13.3f %13.3f %13.3f \n",
 296            part->Px(), part->Py(), part->Pz(), part->Energy());
 297     printf("new       %13.3f %13.3f %13.3f %13.3f \n",
 298            dp[0]->Px()+dp[1]->Px()+dp[2]->Px(),
 299            dp[0]->Py()+dp[1]->Py()+dp[2]->Py(),
 300            dp[0]->Pz()+dp[1]->Pz()+dp[2]->Pz(),
 301            dp[0]->Energy()+dp[1]->Energy()+dp[2]->Energy());
 302
 303
 304   }
 305 }
 306 AliOmegaDalitz& AliOmegaDalitz::operator=(const  AliOmegaDalitz& rhs)
 307 {
 308 // Assignment operator
 309     rhs.Copy(*this);
 310     return *this;
 311 }
 312
 313 void AliOmegaDalitz::Copy(TObject&) const
 314 {
 315     //
 316     // Copy
 317     //
 318     Fatal("Copy","Not implemented!\n");
 319 }
 320
 321 AliOmegaDalitz::AliOmegaDalitz(const AliOmegaDalitz &dalitz)
 322   : AliDecayer(),
 323     fEPMass(0),
 324     fMPMass(0),
 325     fInit(0)
 326 {
 327   // Copy constructor
 328   dalitz.Copy(*this);
 329 }