Coverity fix.
[u/mrichter/AliRoot.git] / 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 }