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