Update timestamps for new AMANDA simulation (17/02/2015)
[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 costheta, sintheta, cosphi, sinphi, phi;
142   
143   // Get the particle masses
144   // lepton
145   lmass = (TDatabasePDG::Instance()->GetParticle(idlepton))->Mass();
146   // omega
147   pmass = pparent->M();
148   // pi0
149   omass = (TDatabasePDG::Instance()->GetParticle(kPi0))     ->Mass();
150
151   // Sample the lepton pair mass from a histogram
152   for( ;; ) 
153   {
154     if ( idlepton == kElectron )
155       lpmass = fEPMass->GetRandom();
156     else
157       lpmass = fMPMass->GetRandom();
158     if ( pmass - omass > lpmass && lpmass / 2. > lmass ) break;
159   }
160   
161   // lepton pair kinematics in virtual photon rest frame
162   e1 = lpmass / 2.;
163   p1 = TMath::Sqrt((e1 + lmass) * (e1 - lmass));
164   // betaSquare = 1.0 - 4.0 * (lmass * lmass) / (lpmass * lpmass);
165   // lambda      = betaSquare / (2.0 - betaSquare);
166   costheta = (2.0 * gRandom->Rndm()) - 1.;
167   sintheta = TMath::Sqrt((1. + costheta) * (1. - costheta));
168   phi      = 2.0 * TMath::ACos(-1.) * gRandom->Rndm();
169   sinphi   = TMath::Sin(phi);
170   cosphi   = TMath::Cos(phi);
171   // momentum vectors of leptons in virtual photon rest frame
172   Double_t pProd1[3] = {p1 * sintheta * cosphi, 
173                        p1 * sintheta * sinphi, 
174                        p1 * costheta};
175   Double_t pProd2[3] = {-1.0 * p1 * sintheta * cosphi, 
176                        -1.0 * p1 * sintheta * sinphi, 
177                        -1.0 * p1 * costheta};
178   
179   // pizero kinematics in omega rest frame
180   e3       = (pmass * pmass + omass * omass - lpmass * lpmass)/(2. * pmass);
181   p3       = TMath::Sqrt((e3 + omass)  * (e3 - omass));
182   costheta = (2.0 * gRandom->Rndm()) - 1.;
183   sintheta = TMath::Sqrt((1. + costheta) * (1. - costheta));
184   phi      = 2.0 * TMath::ACos(-1.) * gRandom->Rndm();
185   sinphi   = TMath::Sin(phi);
186   cosphi   = TMath::Cos(phi);       
187   // pizero 4-vector in omega rest frame
188   fProducts[2].SetPx(p3 * sintheta * cosphi);
189   fProducts[2].SetPy(p3 * sintheta * sinphi);
190   fProducts[2].SetPz(p3 * costheta);
191   fProducts[2].SetE(e3);
192
193   // lepton 4-vectors in properly rotated virtual photon rest frame
194   Double_t pRot1[3] = {0.};
195   Rot(pProd1, pRot1, costheta, -sintheta, -cosphi, -sinphi);
196   Double_t pRot2[3] = {0.};
197   Rot(pProd2, pRot2, costheta, -sintheta, -cosphi, -sinphi); 
198   fProducts[0].SetPx(pRot1[0]);
199   fProducts[0].SetPy(pRot1[1]);
200   fProducts[0].SetPz(pRot1[2]);
201   fProducts[0].SetE(e1);
202   fProducts[1].SetPx(pRot2[0]);
203   fProducts[1].SetPy(pRot2[1]);
204   fProducts[1].SetPz(pRot2[2]);
205   fProducts[1].SetE(e1);
206
207   // boost the dilepton into the omega's rest frame 
208   Double_t eLPparent = TMath::Sqrt(p3 * p3 + lpmass * lpmass);
209   TVector3 boostPair( -1.0 * fProducts[2].Px() / eLPparent, 
210                       -1.0 * fProducts[2].Py() / eLPparent,
211                       -1.0 * fProducts[2].Pz() / eLPparent);
212   fProducts[0].Boost(boostPair);
213   fProducts[1].Boost(boostPair);
214
215   // boost all decay products into the lab frame 
216   TVector3 boostLab( pparent->Px() / pparent->E(), 
217                      pparent->Py() / pparent->E(),
218                      pparent->Pz() / pparent->E());
219   fProducts[0].Boost(boostLab);
220   fProducts[1].Boost(boostLab);
221   fProducts[2].Boost(boostLab);
222   
223   return;
224
225 }
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     Int_t pdg   [3] = {kElectron, -kElectron, kPi0};
235     if ( fProducts[1].M() > 0.001 ) 
236     {  
237       pdg[0] =  kMuonPlus;
238       pdg[1] = -kMuonPlus;
239     }
240     Int_t parent[3] = {0, 0, -1};
241     Int_t d1    [3] = {-1, -1, 1};
242     Int_t d2    [3] = {-1, -1, 2};
243     for (Int_t i = 2; i > -1; i--) {
244         Double_t px = fProducts[i].Px();
245         Double_t py = fProducts[i].Py();
246         Double_t pz = fProducts[i].Pz();
247         Double_t e  = fProducts[i].E();
248         //
249         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     return (3);
252 }
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   pout[0] = pin[0]*costheta*cosphi-pin[1]*sinphi+pin[2]*sintheta*cosphi;
261   pout[1] = pin[0]*costheta*sinphi+pin[1]*cosphi+pin[2]*sintheta*sinphi;
262   pout[2] = -1.0  * pin[0] * sintheta + pin[2] * costheta;
263   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   printf("-->Correcting Dalitz \n");
272   Int_t nt = array->GetEntriesFast();
273   TParticle* dp[3];
274   for (Int_t i = 0; i < nt; i++) {
275     TParticle* part = (TParticle*) (array->At(i));
276     if (part->GetPdgCode() != 223)                     continue;
277     
278     Int_t fd = part->GetFirstDaughter() - 1;
279     Int_t ld = part->GetLastDaughter()  - 1;
280
281     if (fd < 0)                                        continue;
282     if ((ld - fd) != 2)                                continue;
283
284     for (Int_t j = 0; j < 3; j++) dp[j] = (TParticle*) (array->At(fd+j));
285     if ((dp[0]->GetPdgCode() != kPi0) ||
286         ((TMath::Abs(dp[1]->GetPdgCode()) != kElectron) &&
287          (TMath::Abs(dp[1]->GetPdgCode()) != kMuonPlus))) continue;
288     TLorentzVector omega(part->Px(), part->Py(), part->Pz(), part->Energy());
289     if ( TMath::Abs(dp[1]->GetPdgCode()) == kElectron ) 
290       Decay(kElectron, &omega);
291     else
292       Decay(kMuonPlus, &omega);
293     for (Int_t j = 0; j < 3; j++) dp[j]->SetMomentum(fProducts[2-j]);
294     printf("original %13.3f %13.3f %13.3f %13.3f \n", 
295            part->Px(), part->Py(), part->Pz(), part->Energy());
296     printf("new       %13.3f %13.3f %13.3f %13.3f \n", 
297            dp[0]->Px()+dp[1]->Px()+dp[2]->Px(), 
298            dp[0]->Py()+dp[1]->Py()+dp[2]->Py(), 
299            dp[0]->Pz()+dp[1]->Pz()+dp[2]->Pz(), 
300            dp[0]->Energy()+dp[1]->Energy()+dp[2]->Energy());
301
302
303   }
304 }
305 AliOmegaDalitz& AliOmegaDalitz::operator=(const  AliOmegaDalitz& rhs)
306 {
307 // Assignment operator
308     rhs.Copy(*this);
309     return *this;
310 }
311
312 void AliOmegaDalitz::Copy(TObject&) const
313 {
314     //
315     // Copy 
316     //
317     Fatal("Copy","Not implemented!\n");
318 }
319
320 AliOmegaDalitz::AliOmegaDalitz(const AliOmegaDalitz &dalitz)
321   : AliDecayer(),
322     fEPMass(0),
323     fMPMass(0),
324     fInit(0)
325 {
326   // Copy constructor
327   dalitz.Copy(*this);
328 }