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 costheta, sintheta, cosphi, sinphi, phi;
 142
 143   // Get the particle masses
 144   // lepton
 145   lmass = (TDatabasePDG::Instance()->GetParticle(idlepton))->Mass();
 146   // omega
 147   pmass = pparent->M();
 148   // pi0
 149   omass = (TDatabasePDG::Instance()->GetParticle(kPi0))     ->Mass();
 150
 151   // Sample the lepton pair mass from a histogram
 152   for( ;; )
 153   {
 154     if ( idlepton == kElectron )
 155       lpmass = fEPMass->GetRandom();
 156     else
 157       lpmass = fMPMass->GetRandom();
 158     if ( pmass - omass > lpmass && lpmass / 2. > lmass ) break;
 159   }
 160
 161   // lepton pair kinematics in virtual photon rest frame
 162   e1 = lpmass / 2.;
 163   p1 = TMath::Sqrt((e1 + lmass) * (e1 - lmass));
 164   // betaSquare = 1.0 - 4.0 * (lmass * lmass) / (lpmass * lpmass);
 165   // lambda      = betaSquare / (2.0 - betaSquare);
 166   costheta = (2.0 * gRandom->Rndm()) - 1.;
 167   sintheta = TMath::Sqrt((1. + costheta) * (1. - costheta));
 168   phi      = 2.0 * TMath::ACos(-1.) * gRandom->Rndm();
 169   sinphi   = TMath::Sin(phi);
 170   cosphi   = TMath::Cos(phi);
 171   // momentum vectors of leptons in virtual photon rest frame
 172   Double_t pProd1[3] = {p1 * sintheta * cosphi,
 173                        p1 * sintheta * sinphi,
 174                        p1 * costheta};
 175   Double_t pProd2[3] = {-1.0 * p1 * sintheta * cosphi,
 176                        -1.0 * p1 * sintheta * sinphi,
 177                        -1.0 * p1 * costheta};
 178
 179   // pizero kinematics in omega rest frame
 180   e3       = (pmass * pmass + omass * omass - lpmass * lpmass)/(2. * pmass);
 181   p3       = TMath::Sqrt((e3 + omass)  * (e3 - omass));
 182   costheta = (2.0 * gRandom->Rndm()) - 1.;
 183   sintheta = TMath::Sqrt((1. + costheta) * (1. - costheta));
 184   phi      = 2.0 * TMath::ACos(-1.) * gRandom->Rndm();
 185   sinphi   = TMath::Sin(phi);
 186   cosphi   = TMath::Cos(phi);
 187   // pizero 4-vector in omega rest frame
 188   fProducts[2].SetPx(p3 * sintheta * cosphi);
 189   fProducts[2].SetPy(p3 * sintheta * sinphi);
 190   fProducts[2].SetPz(p3 * costheta);
 191   fProducts[2].SetE(e3);
 192
 193   // lepton 4-vectors in properly rotated virtual photon rest frame
 194   Double_t pRot1[3] = {0.};
 195   Rot(pProd1, pRot1, costheta, -sintheta, -cosphi, -sinphi);
 196   Double_t pRot2[3] = {0.};
 197   Rot(pProd2, pRot2, costheta, -sintheta, -cosphi, -sinphi);
 198   fProducts[0].SetPx(pRot1[0]);
 199   fProducts[0].SetPy(pRot1[1]);
 200   fProducts[0].SetPz(pRot1[2]);
 201   fProducts[0].SetE(e1);
 202   fProducts[1].SetPx(pRot2[0]);
 203   fProducts[1].SetPy(pRot2[1]);
 204   fProducts[1].SetPz(pRot2[2]);
 205   fProducts[1].SetE(e1);
 206
 207   // boost the dilepton into the omega's rest frame
 208   Double_t eLPparent = TMath::Sqrt(p3 * p3 + lpmass * lpmass);
 209   TVector3 boostPair( -1.0 * fProducts[2].Px() / eLPparent,
 210                       -1.0 * fProducts[2].Py() / eLPparent,
 211                       -1.0 * fProducts[2].Pz() / eLPparent);
 212   fProducts[0].Boost(boostPair);
 213   fProducts[1].Boost(boostPair);
 214
 215   // boost all decay products into the lab frame
 216   TVector3 boostLab( pparent->Px() / pparent->E(),
 217                      pparent->Py() / pparent->E(),
 218                      pparent->Pz() / pparent->E());
 219   fProducts[0].Boost(boostLab);
 220   fProducts[1].Boost(boostLab);
 221   fProducts[2].Boost(boostLab);
 222
 223   return;
 224
 225 }
 226
 227 Int_t AliOmegaDalitz::ImportParticles(TClonesArray *particles)
 228 {
 229     // Import TParticles in array particles
 230     // l+ l- pi0
 231
 232     TClonesArray &clonesParticles = *particles;
 233
 234     Int_t pdg   [3] = {kElectron, -kElectron, kPi0};
 235     if ( fProducts[1].M() > 0.001 )
 236     {
 237       pdg[0] =  kMuonPlus;
 238       pdg[1] = -kMuonPlus;
 239     }
 240     Int_t parent[3] = {0, 0, -1};
 241     Int_t d1    [3] = {-1, -1, 1};
 242     Int_t d2    [3] = {-1, -1, 2};
 243     for (Int_t i = 2; i > -1; i--) {
 244         Double_t px = fProducts[i].Px();
 245         Double_t py = fProducts[i].Py();
 246         Double_t pz = fProducts[i].Pz();
 247         Double_t e  = fProducts[i].E();
 248         //
 249         new(clonesParticles[2 - i]) TParticle(pdg[i], 1, parent[i], -1, d1[i], d2[i], px, py, pz, e, 0., 0., 0., 0.);
 250     }
 251     return (3);
 252 }
 253
 254
 255 void AliOmegaDalitz::
 256 Rot(Double_t pin[3], Double_t pout[3], Double_t costheta, Double_t sintheta,
 257     Double_t cosphi, Double_t sinphi) const
 258 {
 259 // Perform rotation
 260   pout[0] = pin[0]*costheta*cosphi-pin[1]*sinphi+pin[2]*sintheta*cosphi;
 261   pout[1] = pin[0]*costheta*sinphi+pin[1]*cosphi+pin[2]*sintheta*sinphi;
 262   pout[2] = -1.0  * pin[0] * sintheta + pin[2] * costheta;
 263   return;
 264 }
 265
 266 void AliOmegaDalitz::Decay(TClonesArray* array)
 267 {
 268   //
 269   // Replace all omega dalitz decays with the correct matrix element decays
 270   //
 271   printf("-->Correcting Dalitz \n");
 272   Int_t nt = array->GetEntriesFast();
 273   TParticle* dp[3];
 274   for (Int_t i = 0; i < nt; i++) {
 275     TParticle* part = (TParticle*) (array->At(i));
 276     if (part->GetPdgCode() != 223)                     continue;
 277
 278     Int_t fd = part->GetFirstDaughter() - 1;
 279     Int_t ld = part->GetLastDaughter()  - 1;
 280
 281     if (fd < 0)                                        continue;
 282     if ((ld - fd) != 2)                                continue;
 283
 284     for (Int_t j = 0; j < 3; j++) dp[j] = (TParticle*) (array->At(fd+j));
 285     if ((dp[0]->GetPdgCode() != kPi0) ||
 286         ((TMath::Abs(dp[1]->GetPdgCode()) != kElectron) &&
 287          (TMath::Abs(dp[1]->GetPdgCode()) != kMuonPlus))) continue;
 288     TLorentzVector omega(part->Px(), part->Py(), part->Pz(), part->Energy());
 289     if ( TMath::Abs(dp[1]->GetPdgCode()) == kElectron )
 290       Decay(kElectron, &omega);
 291     else
 292       Decay(kMuonPlus, &omega);
 293     for (Int_t j = 0; j < 3; j++) dp[j]->SetMomentum(fProducts[2-j]);
 294     printf("original %13.3f %13.3f %13.3f %13.3f \n",
 295            part->Px(), part->Py(), part->Pz(), part->Energy());
 296     printf("new       %13.3f %13.3f %13.3f %13.3f \n",
 297            dp[0]->Px()+dp[1]->Px()+dp[2]->Px(),
 298            dp[0]->Py()+dp[1]->Py()+dp[2]->Py(),
 299            dp[0]->Pz()+dp[1]->Pz()+dp[2]->Pz(),
 300            dp[0]->Energy()+dp[1]->Energy()+dp[2]->Energy());
 301
 302
 303   }
 304 }
 305 AliOmegaDalitz& AliOmegaDalitz::operator=(const  AliOmegaDalitz& rhs)
 306 {
 307 // Assignment operator
 308     rhs.Copy(*this);
 309     return *this;
 310 }
 311
 312 void AliOmegaDalitz::Copy(TObject&) const
 313 {
 314     //
 315     // Copy
 316     //
 317     Fatal("Copy","Not implemented!\n");
 318 }
 319
 320 AliOmegaDalitz::AliOmegaDalitz(const AliOmegaDalitz &dalitz)
 321   : AliDecayer(),
 322     fEPMass(0),
 323     fMPMass(0),
 324     fInit(0)
 325 {
 326   // Copy constructor
 327   dalitz.Copy(*this);
 328 }