]> git.uio.no Git - u/mrichter/AliRoot.git/blame - EVGEN/AliOmegaDalitz.cxx
Adding more histograms
[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(),
51 fLPMass(0)
52{
53 // Constructor
54}
55
56void 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;
8c90f5a2 105 fLPMass->AddBinContent(ibin, weight);
106 }
107}
108
109
110void 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
d060703c 127 pmass = pparent->M();
8c90f5a2 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
204Int_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
228void AliOmegaDalitz::
229Rot(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
05b093f3 238void 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}
273AliOmegaDalitz& AliOmegaDalitz::operator=(const AliOmegaDalitz& rhs)
274{
275// Assignment operator
276 rhs.Copy(*this);
277 return *this;
278}
279
280void AliOmegaDalitz::Copy(TObject&) const
281{
282 //
283 // Copy
284 //
285 Fatal("Copy","Not implemented!\n");
286}
287
288AliOmegaDalitz::AliOmegaDalitz(const AliOmegaDalitz &dalitz)
289 : AliDecayer(),
290 fLPMass(0)
291{
292 // Copy constructor
293 dalitz.Copy(*this);
294}