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