Line data Source code
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 6 : 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 7 : AliOmegaDalitz::AliOmegaDalitz():
50 1 : AliDecayer(),
51 1 : fEPMass(0),
52 1 : fMPMass(0),
53 1 : fInit(0)
54 5 : {
55 : // Constructor
56 2 : }
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 0 : emass = (TDatabasePDG::Instance()->GetParticle(kElectron))->Mass();
69 : // muon
70 0 : mmass = (TDatabasePDG::Instance()->GetParticle(kMuonPlus))->Mass();
71 : // omega
72 0 : pmass = (TDatabasePDG::Instance()->GetParticle(223)) ->Mass();
73 : // pi0
74 0 : omass = (TDatabasePDG::Instance()->GetParticle(kPi0)) ->Mass();
75 :
76 0 : for ( idecay = 1; idecay<=2; idecay++ )
77 : {
78 0 : if ( idecay == 1 )
79 0 : lmass = emass;
80 : else
81 : lmass = mmass;
82 :
83 0 : min = 2.0 * lmass;
84 0 : max = pmass - omass;
85 0 : binwidth = (max - min) / (Double_t)nbins;
86 0 : if ( idecay == 1 )
87 0 : fEPMass = new TH1F("fEPMass", "Dalitz electron pair", nbins, min, max);
88 : else
89 0 : fMPMass = new TH1F("fMPMass", "Dalitz muon pair", nbins, min, max);
90 :
91 0 : epsilon = (lmass / pmass) * (lmass / pmass);
92 0 : delta = (omass / pmass) * (omass / pmass);
93 :
94 0 : for ( ibin = 1; ibin <= nbins; ibin++ )
95 : {
96 0 : mLL = min + (Double_t)(ibin - 1) * binwidth + binwidth / 2.0;
97 0 : q = (mLL / pmass) * (mLL / pmass);
98 0 : if ( q <= 4.0 * epsilon )
99 : {
100 0 : AliFatal("Error in calculating Dalitz mass histogram binning!");
101 0 : }
102 :
103 0 : kwHelp = (1.0 + q / (1.0 - delta)) * (1.0 + q / (1.0 - delta))
104 0 : - 4.0 * q / ((1.0 - delta) * (1.0 - delta));
105 0 : if ( kwHelp <= 0.0 )
106 : {
107 0 : AliFatal("Error in calculating Dalitz mass histogram binning!");
108 0 : }
109 0 : krollWada = (2.0 / mLL) * TMath::Exp(1.5 * TMath::Log(kwHelp))
110 0 : * TMath::Sqrt(1.0 - 4.0 * epsilon / q)
111 0 : * (1.0 + 2.0 * epsilon / q);
112 : formFactor =
113 0 : (TMath::Power(TMath::Power(0.6519,2),2))
114 0 : / (TMath::Power(TMath::Power(0.6519,2)-TMath::Power(mLL, 2), 2)
115 0 : + TMath::Power(0.04198, 2)*TMath::Power(0.6519, 2));
116 0 : weight = krollWada * formFactor;
117 0 : if ( idecay == 1 )
118 0 : fEPMass->AddBinContent(ibin, weight);
119 : else
120 0 : fMPMass->AddBinContent(ibin, weight);
121 : }
122 : }
123 0 : }
124 :
125 :
126 : void AliOmegaDalitz::Decay(Int_t idlepton, TLorentzVector* pparent)
127 : {
128 : //-----------------------------------------------------------------------------
129 : //
130 : // Generate Omega Dalitz decay
131 : //
132 : //-----------------------------------------------------------------------------
133 :
134 0 : if (!fInit) {
135 0 : Init();
136 0 : fInit=1;
137 0 : }
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 0 : lmass = (TDatabasePDG::Instance()->GetParticle(idlepton))->Mass();
146 : // omega
147 0 : pmass = pparent->M();
148 : // pi0
149 0 : omass = (TDatabasePDG::Instance()->GetParticle(kPi0)) ->Mass();
150 :
151 : // Sample the lepton pair mass from a histogram
152 0 : for( ;; )
153 : {
154 0 : if ( idlepton == kElectron )
155 0 : lpmass = fEPMass->GetRandom();
156 : else
157 0 : lpmass = fMPMass->GetRandom();
158 0 : if ( pmass - omass > lpmass && lpmass / 2. > lmass ) break;
159 : }
160 :
161 : // lepton pair kinematics in virtual photon rest frame
162 : e1 = lpmass / 2.;
163 0 : p1 = TMath::Sqrt((e1 + lmass) * (e1 - lmass));
164 : // betaSquare = 1.0 - 4.0 * (lmass * lmass) / (lpmass * lpmass);
165 : // lambda = betaSquare / (2.0 - betaSquare);
166 0 : costheta = (2.0 * gRandom->Rndm()) - 1.;
167 0 : sintheta = TMath::Sqrt((1. + costheta) * (1. - costheta));
168 0 : phi = 2.0 * TMath::ACos(-1.) * gRandom->Rndm();
169 0 : sinphi = TMath::Sin(phi);
170 0 : cosphi = TMath::Cos(phi);
171 : // momentum vectors of leptons in virtual photon rest frame
172 0 : Double_t pProd1[3] = {p1 * sintheta * cosphi,
173 0 : p1 * sintheta * sinphi,
174 0 : p1 * costheta};
175 0 : Double_t pProd2[3] = {-1.0 * p1 * sintheta * cosphi,
176 0 : -1.0 * p1 * sintheta * sinphi,
177 0 : -1.0 * p1 * costheta};
178 :
179 : // pizero kinematics in omega rest frame
180 0 : e3 = (pmass * pmass + omass * omass - lpmass * lpmass)/(2. * pmass);
181 0 : p3 = TMath::Sqrt((e3 + omass) * (e3 - omass));
182 0 : costheta = (2.0 * gRandom->Rndm()) - 1.;
183 0 : sintheta = TMath::Sqrt((1. + costheta) * (1. - costheta));
184 0 : phi = 2.0 * TMath::ACos(-1.) * gRandom->Rndm();
185 0 : sinphi = TMath::Sin(phi);
186 0 : cosphi = TMath::Cos(phi);
187 : // pizero 4-vector in omega rest frame
188 0 : fProducts[2].SetPx(p3 * sintheta * cosphi);
189 0 : fProducts[2].SetPy(p3 * sintheta * sinphi);
190 0 : fProducts[2].SetPz(p3 * costheta);
191 0 : fProducts[2].SetE(e3);
192 :
193 : // lepton 4-vectors in properly rotated virtual photon rest frame
194 0 : Double_t pRot1[3] = {0.};
195 0 : Rot(pProd1, pRot1, costheta, -sintheta, -cosphi, -sinphi);
196 0 : Double_t pRot2[3] = {0.};
197 0 : Rot(pProd2, pRot2, costheta, -sintheta, -cosphi, -sinphi);
198 0 : fProducts[0].SetPx(pRot1[0]);
199 0 : fProducts[0].SetPy(pRot1[1]);
200 0 : fProducts[0].SetPz(pRot1[2]);
201 0 : fProducts[0].SetE(e1);
202 0 : fProducts[1].SetPx(pRot2[0]);
203 0 : fProducts[1].SetPy(pRot2[1]);
204 0 : fProducts[1].SetPz(pRot2[2]);
205 0 : fProducts[1].SetE(e1);
206 :
207 : // boost the dilepton into the omega's rest frame
208 0 : Double_t eLPparent = TMath::Sqrt(p3 * p3 + lpmass * lpmass);
209 0 : TVector3 boostPair( -1.0 * fProducts[2].Px() / eLPparent,
210 0 : -1.0 * fProducts[2].Py() / eLPparent,
211 0 : -1.0 * fProducts[2].Pz() / eLPparent);
212 0 : fProducts[0].Boost(boostPair);
213 0 : fProducts[1].Boost(boostPair);
214 :
215 : // boost all decay products into the lab frame
216 0 : TVector3 boostLab( pparent->Px() / pparent->E(),
217 0 : pparent->Py() / pparent->E(),
218 0 : pparent->Pz() / pparent->E());
219 0 : fProducts[0].Boost(boostLab);
220 0 : fProducts[1].Boost(boostLab);
221 0 : fProducts[2].Boost(boostLab);
222 :
223 : return;
224 :
225 0 : }
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 0 : Int_t pdg [3] = {kElectron, -kElectron, kPi0};
235 0 : if ( fProducts[1].M() > 0.001 )
236 : {
237 0 : pdg[0] = kMuonPlus;
238 0 : pdg[1] = -kMuonPlus;
239 0 : }
240 0 : Int_t parent[3] = {0, 0, -1};
241 0 : Int_t d1 [3] = {-1, -1, 1};
242 0 : Int_t d2 [3] = {-1, -1, 2};
243 0 : for (Int_t i = 2; i > -1; i--) {
244 0 : Double_t px = fProducts[i].Px();
245 0 : Double_t py = fProducts[i].Py();
246 0 : Double_t pz = fProducts[i].Pz();
247 0 : Double_t e = fProducts[i].E();
248 : //
249 0 : 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 0 : return (3);
252 0 : }
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 0 : pout[0] = pin[0]*costheta*cosphi-pin[1]*sinphi+pin[2]*sintheta*cosphi;
261 0 : pout[1] = pin[0]*costheta*sinphi+pin[1]*cosphi+pin[2]*sintheta*sinphi;
262 0 : pout[2] = -1.0 * pin[0] * sintheta + pin[2] * costheta;
263 0 : 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 0 : printf("-->Correcting Dalitz \n");
272 0 : Int_t nt = array->GetEntriesFast();
273 0 : TParticle* dp[3];
274 0 : for (Int_t i = 0; i < nt; i++) {
275 0 : TParticle* part = (TParticle*) (array->At(i));
276 0 : if (part->GetPdgCode() != 223) continue;
277 :
278 0 : Int_t fd = part->GetFirstDaughter() - 1;
279 0 : Int_t ld = part->GetLastDaughter() - 1;
280 :
281 0 : if (fd < 0) continue;
282 0 : if ((ld - fd) != 2) continue;
283 :
284 0 : for (Int_t j = 0; j < 3; j++) dp[j] = (TParticle*) (array->At(fd+j));
285 0 : if ((dp[0]->GetPdgCode() != kPi0) ||
286 0 : ((TMath::Abs(dp[1]->GetPdgCode()) != kElectron) &&
287 0 : (TMath::Abs(dp[1]->GetPdgCode()) != kMuonPlus))) continue;
288 0 : TLorentzVector omega(part->Px(), part->Py(), part->Pz(), part->Energy());
289 0 : if ( TMath::Abs(dp[1]->GetPdgCode()) == kElectron )
290 0 : Decay(kElectron, &omega);
291 : else
292 0 : Decay(kMuonPlus, &omega);
293 0 : for (Int_t j = 0; j < 3; j++) dp[j]->SetMomentum(fProducts[2-j]);
294 0 : printf("original %13.3f %13.3f %13.3f %13.3f \n",
295 0 : part->Px(), part->Py(), part->Pz(), part->Energy());
296 0 : printf("new %13.3f %13.3f %13.3f %13.3f \n",
297 0 : dp[0]->Px()+dp[1]->Px()+dp[2]->Px(),
298 0 : dp[0]->Py()+dp[1]->Py()+dp[2]->Py(),
299 0 : dp[0]->Pz()+dp[1]->Pz()+dp[2]->Pz(),
300 0 : dp[0]->Energy()+dp[1]->Energy()+dp[2]->Energy());
301 :
302 :
303 0 : }
304 0 : }
305 : AliOmegaDalitz& AliOmegaDalitz::operator=(const AliOmegaDalitz& rhs)
306 : {
307 : // Assignment operator
308 0 : rhs.Copy(*this);
309 0 : return *this;
310 : }
311 :
312 : void AliOmegaDalitz::Copy(TObject&) const
313 : {
314 : //
315 : // Copy
316 : //
317 0 : Fatal("Copy","Not implemented!\n");
318 0 : }
319 :
320 0 : AliOmegaDalitz::AliOmegaDalitz(const AliOmegaDalitz &dalitz)
321 0 : : AliDecayer(),
322 0 : fEPMass(0),
323 0 : fMPMass(0),
324 0 : fInit(0)
325 0 : {
326 : // Copy constructor
327 0 : dalitz.Copy(*this);
328 0 : }
|
