Added AliGenHepMCEventHeader class and put it in STEERBase library. Updated AliGenRea...
[u/mrichter/AliRoot.git] / EVGEN / AliOmegaDalitz.cxx
CommitLineData
8c90f5a2 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
31ClassImp(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//
49AliOmegaDalitz::AliOmegaDalitz():
50 AliDecayer(),
0439b2b7 51 fEPMass(0),
b8f95537 52 fMPMass(0),
53 fInit(0)
8c90f5a2 54{
55 // Constructor
56}
57
58void AliOmegaDalitz::Init()
59{
60 // Initialisation
0439b2b7 61 Int_t idecay, ibin, nbins = 1000;
8c90f5a2 62 Double_t min, max, binwidth;
0439b2b7 63 Double_t pmass, lmass, omass, emass, mmass;
8c90f5a2 64 Double_t epsilon, delta, mLL, q, kwHelp, krollWada, formFactor, weight;
65
66 // Get the particle masses
67 // electron
0439b2b7 68 emass = (TDatabasePDG::Instance()->GetParticle(kElectron))->Mass();
69 // muon
70 mmass = (TDatabasePDG::Instance()->GetParticle(kMuonPlus))->Mass();
8c90f5a2 71 // omega
72 pmass = (TDatabasePDG::Instance()->GetParticle(223)) ->Mass();
73 // pi0
74 omass = (TDatabasePDG::Instance()->GetParticle(kPi0)) ->Mass();
75
0439b2b7 76 for ( idecay = 1; idecay<=2; idecay++ )
8c90f5a2 77 {
0439b2b7 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++ )
8c90f5a2 95 {
0439b2b7 96 mLL = min + (Double_t)(ibin - 1) * binwidth + binwidth / 2.0;
97 q = (mLL / pmass) * (mLL / pmass);
98 if ( q <= 4.0 * epsilon )
99 {
8c90f5a2 100 AliFatal("Error in calculating Dalitz mass histogram binning!");
0439b2b7 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 {
8c90f5a2 107 AliFatal("Error in calculating Dalitz mass histogram binning!");
0439b2b7 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 =
8c90f5a2 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));
0439b2b7 116 weight = krollWada * formFactor;
117 if ( idecay == 1 )
118 fEPMass->AddBinContent(ibin, weight);
119 else
120 fMPMass->AddBinContent(ibin, weight);
121 }
8c90f5a2 122 }
123}
124
125
cc1f7ac3 126void AliOmegaDalitz::Decay(Int_t idlepton, TLorentzVector* pparent)
8c90f5a2 127{
128//-----------------------------------------------------------------------------
129//
130// Generate Omega Dalitz decay
131//
132//-----------------------------------------------------------------------------
133
b8f95537 134 if (!fInit) {
135 Init();
136 fInit=1;
137 }
138
8c90f5a2 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
0439b2b7 145 // lepton
146 lmass = (TDatabasePDG::Instance()->GetParticle(idlepton))->Mass();
8c90f5a2 147 // omega
d060703c 148 pmass = pparent->M();
8c90f5a2 149 // pi0
150 omass = (TDatabasePDG::Instance()->GetParticle(kPi0)) ->Mass();
151
152 // Sample the lepton pair mass from a histogram
153 for( ;; )
154 {
0439b2b7 155 if ( idlepton == kElectron )
156 lpmass = fEPMass->GetRandom();
157 else
158 lpmass = fMPMass->GetRandom();
159 if ( pmass - omass > lpmass && lpmass / 2. > lmass ) break;
8c90f5a2 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
228Int_t AliOmegaDalitz::ImportParticles(TClonesArray *particles)
229{
230 // Import TParticles in array particles
0439b2b7 231 // l+ l- pi0
8c90f5a2 232
233 TClonesArray &clonesParticles = *particles;
234
0439b2b7 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};
8c90f5a2 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
8c90f5a2 256void AliOmegaDalitz::
257Rot(Double_t pin[3], Double_t pout[3], Double_t costheta, Double_t sintheta,
271e6aae 258 Double_t cosphi, Double_t sinphi) const
8c90f5a2 259{
43ce81af 260// Perform rotation
8c90f5a2 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
05b093f3 267void 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) ||
0439b2b7 287 ((TMath::Abs(dp[1]->GetPdgCode()) != kElectron) &&
288 (TMath::Abs(dp[1]->GetPdgCode()) != kMuonPlus))) continue;
05b093f3 289 TLorentzVector omega(part->Px(), part->Py(), part->Pz(), part->Energy());
0439b2b7 290 if ( TMath::Abs(dp[1]->GetPdgCode()) == kElectron )
291 Decay(kElectron, &omega);
292 else
293 Decay(kMuonPlus, &omega);
05b093f3 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}
306AliOmegaDalitz& AliOmegaDalitz::operator=(const AliOmegaDalitz& rhs)
307{
308// Assignment operator
309 rhs.Copy(*this);
310 return *this;
311}
312
313void AliOmegaDalitz::Copy(TObject&) const
314{
315 //
316 // Copy
317 //
318 Fatal("Copy","Not implemented!\n");
319}
320
321AliOmegaDalitz::AliOmegaDalitz(const AliOmegaDalitz &dalitz)
322 : AliDecayer(),
0439b2b7 323 fEPMass(0),
b8f95537 324 fMPMass(0),
325 fInit(0)
05b093f3 326{
327 // Copy constructor
328 dalitz.Copy(*this);
329}