Adding histos for the LEGO train
[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         fInit(0)
54 {
55     // Constructor
56 }
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   emass = (TDatabasePDG::Instance()->GetParticle(kElectron))->Mass();
69   // muon
70   mmass = (TDatabasePDG::Instance()->GetParticle(kMuonPlus))->Mass();
71   // omega
72   pmass = (TDatabasePDG::Instance()->GetParticle(223))      ->Mass();
73   // pi0
74   omass = (TDatabasePDG::Instance()->GetParticle(kPi0))     ->Mass();
75
76   for ( idecay = 1; idecay<=2; idecay++ )
77   {
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++ )
95     {
96       mLL = min + (Double_t)(ibin - 1) * binwidth + binwidth / 2.0;
97       q    = (mLL / pmass) * (mLL / pmass);
98       if ( q <= 4.0 * epsilon )
99       {
100         AliFatal("Error in calculating Dalitz mass histogram binning!");
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       {
107         AliFatal("Error in calculating Dalitz mass histogram binning!");
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 = 
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));
116       weight = krollWada * formFactor;   
117       if ( idecay == 1 ) 
118         fEPMass->AddBinContent(ibin, weight);
119       else
120         fMPMass->AddBinContent(ibin, weight);
121     }
122   }
123 }
124
125
126 void AliOmegaDalitz::Decay(Int_t idlepton, TLorentzVector* pparent)
127 {
128 //-----------------------------------------------------------------------------
129 //
130 //  Generate Omega Dalitz decay
131 //
132 //-----------------------------------------------------------------------------
133
134     if (!fInit) {
135         Init();
136         fInit=1;
137     }
138     
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
145   // lepton
146   lmass = (TDatabasePDG::Instance()->GetParticle(idlepton))->Mass();
147   // omega
148   pmass = pparent->M();
149   // pi0
150   omass = (TDatabasePDG::Instance()->GetParticle(kPi0))     ->Mass();
151
152   // Sample the lepton pair mass from a histogram
153   for( ;; ) 
154   {
155     if ( idlepton == kElectron )
156       lpmass = fEPMass->GetRandom();
157     else
158       lpmass = fMPMass->GetRandom();
159     if ( pmass - omass > lpmass && lpmass / 2. > lmass ) break;
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
228 Int_t AliOmegaDalitz::ImportParticles(TClonesArray *particles)
229 {
230     // Import TParticles in array particles
231     // l+ l- pi0
232     
233     TClonesArray &clonesParticles = *particles;
234
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};
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
256 void AliOmegaDalitz::
257 Rot(Double_t pin[3], Double_t pout[3], Double_t costheta, Double_t sintheta,
258     Double_t cosphi, Double_t sinphi) const
259 {
260 // Perform rotation
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
267 void 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) ||
287         ((TMath::Abs(dp[1]->GetPdgCode()) != kElectron) &&
288          (TMath::Abs(dp[1]->GetPdgCode()) != kMuonPlus))) continue;
289     TLorentzVector omega(part->Px(), part->Py(), part->Pz(), part->Energy());
290     if ( TMath::Abs(dp[1]->GetPdgCode()) == kElectron ) 
291       Decay(kElectron, &omega);
292     else
293       Decay(kMuonPlus, &omega);
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 }
306 AliOmegaDalitz& AliOmegaDalitz::operator=(const  AliOmegaDalitz& rhs)
307 {
308 // Assignment operator
309     rhs.Copy(*this);
310     return *this;
311 }
312
313 void AliOmegaDalitz::Copy(TObject&) const
314 {
315     //
316     // Copy 
317     //
318     Fatal("Copy","Not implemented!\n");
319 }
320
321 AliOmegaDalitz::AliOmegaDalitz(const AliOmegaDalitz &dalitz)
322   : AliDecayer(),
323     fEPMass(0),
324     fMPMass(0),
325     fInit(0)
326 {
327   // Copy constructor
328   dalitz.Copy(*this);
329 }