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