Class for Omega Dalitz Decays (R. Averbeck, AM)
[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         fLPMass(0)
52 {
53     // Constructor
54 }
55
56 void 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;   
105     printf(" %5d %13.3f \n", ibin, weight);
106     
107     fLPMass->AddBinContent(ibin, weight);
108   }
109 }
110
111
112 void AliOmegaDalitz::Decay(Int_t /*idpart*/, TLorentzVector* pparent)
113 {
114 //-----------------------------------------------------------------------------
115 //
116 //  Generate Omega Dalitz decay
117 //
118 //-----------------------------------------------------------------------------
119
120   Double_t pmass, lmass, omass, lpmass;
121   Double_t e1, p1, e3, p3;
122   Double_t betaSquare, lambda;
123   Double_t costheta, sintheta, cosphi, sinphi, phi;
124   
125   // Get the particle masses
126   // electron
127   lmass = (TDatabasePDG::Instance()->GetParticle(kElectron))->Mass();
128   // omega
129   pmass = (TDatabasePDG::Instance()->GetParticle(223))      ->Mass();
130   // pi0
131   omass = (TDatabasePDG::Instance()->GetParticle(kPi0))     ->Mass();
132
133   // Sample the lepton pair mass from a histogram
134   for( ;; ) 
135   {
136       lpmass = fLPMass->GetRandom();
137       if ( pmass - omass > lpmass && lpmass / 2. > lmass ) break;
138   }
139   
140   // lepton pair kinematics in virtual photon rest frame
141   e1 = lpmass / 2.;
142   p1 = TMath::Sqrt((e1 + lmass) * (e1 - lmass));
143   betaSquare = 1.0 - 4.0 * (lmass * lmass) / (lpmass * lpmass);
144   lambda      = betaSquare / (2.0 - betaSquare);
145   costheta = (2.0 * gRandom->Rndm()) - 1.;
146   sintheta = TMath::Sqrt((1. + costheta) * (1. - costheta));
147   phi      = 2.0 * TMath::ACos(-1.) * gRandom->Rndm();
148   sinphi   = TMath::Sin(phi);
149   cosphi   = TMath::Cos(phi);
150   // momentum vectors of leptons in virtual photon rest frame
151   Double_t pProd1[3] = {p1 * sintheta * cosphi, 
152                        p1 * sintheta * sinphi, 
153                        p1 * costheta};
154   Double_t pProd2[3] = {-1.0 * p1 * sintheta * cosphi, 
155                        -1.0 * p1 * sintheta * sinphi, 
156                        -1.0 * p1 * costheta};
157   
158   // pizero kinematics in omega rest frame
159   e3       = (pmass * pmass + omass * omass - lpmass * lpmass)/(2. * pmass);
160   p3       = TMath::Sqrt((e3 + omass)  * (e3 - omass));
161   costheta = (2.0 * gRandom->Rndm()) - 1.;
162   sintheta = TMath::Sqrt((1. + costheta) * (1. - costheta));
163   phi      = 2.0 * TMath::ACos(-1.) * gRandom->Rndm();
164   sinphi   = TMath::Sin(phi);
165   cosphi   = TMath::Cos(phi);       
166   // pizero 4-vector in omega rest frame
167   fProducts[2].SetPx(p3 * sintheta * cosphi);
168   fProducts[2].SetPy(p3 * sintheta * sinphi);
169   fProducts[2].SetPz(p3 * costheta);
170   fProducts[2].SetE(e3);
171
172   // lepton 4-vectors in properly rotated virtual photon rest frame
173   Double_t pRot1[3] = {0.};
174   Rot(pProd1, pRot1, costheta, -sintheta, -cosphi, -sinphi);
175   Double_t pRot2[3] = {0.};
176   Rot(pProd2, pRot2, costheta, -sintheta, -cosphi, -sinphi); 
177   fProducts[0].SetPx(pRot1[0]);
178   fProducts[0].SetPy(pRot1[1]);
179   fProducts[0].SetPz(pRot1[2]);
180   fProducts[0].SetE(e1);
181   fProducts[1].SetPx(pRot2[0]);
182   fProducts[1].SetPy(pRot2[1]);
183   fProducts[1].SetPz(pRot2[2]);
184   fProducts[1].SetE(e1);
185
186   // boost the dilepton into the omega's rest frame 
187   Double_t eLPparent = TMath::Sqrt(p3 * p3 + lpmass * lpmass);
188   TVector3 boostPair( -1.0 * fProducts[2].Px() / eLPparent, 
189                       -1.0 * fProducts[2].Py() / eLPparent,
190                       -1.0 * fProducts[2].Pz() / eLPparent);
191   fProducts[0].Boost(boostPair);
192   fProducts[1].Boost(boostPair);
193
194   // boost all decay products into the lab frame 
195   TVector3 boostLab( pparent->Px() / pparent->E(), 
196                      pparent->Py() / pparent->E(),
197                      pparent->Pz() / pparent->E());
198   fProducts[0].Boost(boostLab);
199   fProducts[1].Boost(boostLab);
200   fProducts[2].Boost(boostLab);
201   
202   return;
203
204 }
205
206 Int_t AliOmegaDalitz::ImportParticles(TClonesArray *particles)
207 {
208     // Import TParticles in array particles
209     // e+ e- pi0
210     
211     TClonesArray &clonesParticles = *particles;
212
213     Int_t pdg   [3] = {kPi0, kElectron, -kElectron};
214     Int_t parent[3] = {-1, 0, 0};
215     Int_t d1    [3] = {1, -1, -1};
216     Int_t d2    [3] = {2, -1, -1};
217     for (Int_t i = 2; i > -1; i--) {
218         Double_t px = fProducts[i].Px();
219         Double_t py = fProducts[i].Py();
220         Double_t pz = fProducts[i].Pz();
221         Double_t e  = fProducts[i].E();
222         //
223         new(clonesParticles[2 - i]) TParticle(pdg[i], 1, parent[i], -1, d1[i], d2[i], px, py, pz, e, 0., 0., 0., 0.);
224     }
225     return (3);
226 }
227
228
229
230 void AliOmegaDalitz::
231 Rot(Double_t pin[3], Double_t pout[3], Double_t costheta, Double_t sintheta,
232     Double_t cosphi, Double_t sinphi)
233 {
234   pout[0] = pin[0]*costheta*cosphi-pin[1]*sinphi+pin[2]*sintheta*cosphi;
235   pout[1] = pin[0]*costheta*sinphi+pin[1]*cosphi+pin[2]*sintheta*sinphi;
236   pout[2] = -1.0  * pin[0] * sintheta + pin[2] * costheta;
237   return;
238 }
239