]>
Commit | Line | Data |
---|---|---|
368fc61e | 1 | //#include <TClonesArray.h> |
2 | ||
3 | #include <TDatabasePDG.h> | |
4 | #include <TFile.h> | |
5 | #include "AliConst.h" | |
6 | #include "AliGenMUONLMR.h" | |
7 | #include "AliMC.h" | |
8 | #include "AliRun.h" | |
c279681a | 9 | #include "AliLog.h" |
368fc61e | 10 | #include "AliGenEventHeader.h" |
11 | ||
12 | ClassImp(AliGenMUONLMR) | |
13 | ||
0d19baad | 14 | AliGenMUONLMR::AliGenMUONLMR () : AliGenMC(), |
15 | fNMuMin(2), | |
16 | fCMSEnergy(kNCMSEnergies), | |
17 | fGenSingleProc(-1), | |
18 | fYCM(0), | |
19 | fCosTheta(0x0), | |
20 | fRhoLineShape(0x0), | |
21 | fHMultMu(0x0), | |
22 | fHNProc(0x0) { | |
23 | // | |
24 | // default constructor | |
25 | // | |
26 | for (int i=0; i<fgkNpart; i++) { | |
27 | fPDG[i] = 0; | |
28 | fScaleMult[i] = 1.; | |
29 | fPt[i] = NULL; | |
30 | fY[i] = NULL; | |
31 | fMult[i] = NULL; | |
0d19baad | 32 | fParticle[i] = NULL; |
33 | } | |
14467d5d | 34 | for (int i=0; i<2; i++) { |
35 | fMu[i] = NULL; | |
36 | fDecay[i] = NULL; | |
37 | } | |
38 | ||
39 | for (int i=0; i<3; i++) { | |
40 | fDalitz[i] = NULL; | |
41 | } | |
42 | ||
0d19baad | 43 | } |
35e2c739 | 44 | |
c279681a | 45 | //----------------------------------------------------------- |
368fc61e | 46 | |
c279681a | 47 | void AliGenMUONLMR::SetCMSEnergy(CMSEnergies energy){ |
0d19baad | 48 | fCMSEnergy = energy; |
49 | // initialize pt and y distributions according to a fit to | |
50 | // Pythia simulation at sqrt(s) = 7 TeV | |
51 | for (Int_t ipart=0; ipart < fgkNpart; ipart++) fScaleMult[ipart] = 1; | |
52 | fScaleMult[kPionLMR] = 0; // set pion multiplicity to zero | |
53 | fScaleMult[kKaonLMR] = 0; // set kaon multiplicity to zero | |
54 | const char* fdname[2] = {"fDecPion","fDecKaon"}; | |
55 | Double_t ctau[2] = {7.8045, 3.712}; | |
56 | Int_t pdg[7] = {211, 321, 221, 113, 223, 333, 331}; | |
57 | const char* fptname[7] = {"fPtPion","fPtKaon","fPtEta","fPtRho","fPtOmega","fPtPhi","fPtEtaPrime"}; | |
58 | const char* fyname[7] = {"fYPion","fYKaon","fYEta","fYRho","fYOmega","fYPhi","fYEtaPrime"}; | |
59 | const char* fnname[7] = {"fMultPion","fMultKaon","fMultEta","fMultRho","fMultOmega","fMultPhi","fMultEtaPrime"}; | |
60 | Double_t ptparam[7][9]; | |
61 | Double_t yparam[7][9]; | |
62 | Double_t nparam[7][9]; | |
1f18dff1 | 63 | |
64 | // parameters for 8 TeV generation | |
65 | if (fCMSEnergy==kCMS8000GeV) { | |
66 | AliInfo ("Using pp parameterization at 8 TeV\n"); | |
67 | ||
68 | // Parameters of transverse momentum spectra | |
69 | Double_t ptparam8000[7][9] = {{1,0.427,2.52,0,0,0,0,0,0}, // pions from 7 TeV section of code | |
70 | {1,0.58,2.57,0,0,0,0,0,0}, // kaons from 7 TeV section of code | |
71 | {1,0.657,2.685,0,0,0,0,0,0}, // eta from PYTHIA6.4 ATLAS-CSC at 8 TeV | |
72 | {1,1.44,3.16,0,0,0,0,0,0}, // rho+omega from 7 TeV section of code | |
73 | {1,1.44,3.16,0,0,0,0,0,0}, // rho+omega from 7 TeV section of code | |
74 | {1,1.16,2.74,0,0,0,0,0,0}, // phi from 7 TeV section of code | |
75 | {1,0.755,2.578,0,0,0,0,0,0}}; // etaPrime from PYTHIA 6.4 ATLAS-CSC at 8 TeV | |
76 | ||
77 | // Parameters of rapidity spectra | |
78 | Double_t yparam8000[7][9] = {{1,0.8251,3.657,0,0,0,0,0,0}, // pions from 7 TeV section of code | |
79 | {1,1.83,2.698,0,0,0,0,0,0}, // kaons from 7 TeV section of code | |
80 | {1,0.0509,3.96,0,0,0,0,0,0}, // eta from PYTHIA6.4 ATLAS-CSC at 8 TeV | |
81 | {1,0.0489,3.961,0,0,0,0,0,0}, // rho from PYTHIA6.4 ATLAS-CSC at 8 TeV | |
82 | {1,0.0650,3.966,0,0,0,0,0,0}, // omega from PYTHIA 6.4 ATLAS-CSC at 8 TeV | |
83 | {1,1.279,2.745,0,0,0,0,0,0}, // phi from from PYTHIA6.4 ATLAS-CSC at 8 TeV | |
84 | {1,0.1627,3.883,0,0,0,0,0,0}}; // eta prime from PYTHIA6.4 ATLAS-CSC at 8 TeV | |
85 | ||
86 | // Parameters of multiplicity spectra | |
87 | Double_t nparam8000[7][9] = {{353.582, 6.76263, 1.66979, 998.445, 9.73281, 12.6704, 175.187, 29.08, 40.2531}, //pions from 7 TeV section of code | |
88 | {1.e4, 0.2841, 0,0,0,0,0,0,0}, // kaons from 7 TeV section of code | |
89 | {2.279e4, 0.2622, 0,0,0,0,0,0,0}, // eta from PYTHIA6.4 ATLAS-CSC at 8 TeV | |
90 | {1.564e4, 0.1713, 0,0,0,0,0,0,0}, // rho from PYTHIA6.4 ATLAS-CSC at 8 TeV | |
91 | {1.662e4, 0.183, 0,0,0,0,0,0,0}, // omega from PYTHIA6.4 ATLAS-CSC at 8 TeV | |
92 | {6.723e4, 1.121, 0,0,0,0,0,0,0}, // phi from PYTHIA6.4 ATLAS-CSC at 8 TeV | |
93 | {5.005e4, 0.6971, 0,0,0,0,0,0,0}}; // eta prime from PYTHIA6.4 ATLAS-CSC at 8 TeV | |
94 | ||
95 | for (Int_t i=0; i<fgkNpart; i++) { | |
96 | for (Int_t j=0; j<9; j++) { | |
97 | ptparam[i][j] = ptparam8000[i][j]; | |
98 | yparam[i][j] = yparam8000[i][j]; | |
99 | nparam[i][j] = nparam8000[i][j]; | |
100 | } | |
101 | } | |
102 | } | |
c279681a | 103 | |
0d19baad | 104 | // parameters for 7 TeV generation |
1f18dff1 | 105 | else if (fCMSEnergy==kCMS7000GeV) { |
0d19baad | 106 | AliInfo ("Using pp parameterization at 7 TeV\n"); |
107 | Double_t ptparam7000[7][9] = {{1,0.427,2.52,0,0,0,0,0,0}, // pions from Pythia | |
108 | {1,0.58,2.57,0,0,0,0,0,0}, // kaons from Pythia | |
109 | {1,0.641,2.62,0,0,0,0,0,0}, // eta from Pythia | |
110 | {1,1.44,3.16,0,0,0,0,0,0}, // rho+omega from ALICE muon | |
111 | {1,1.44,3.16,0,0,0,0,0,0}, // rho+omega from ALICE muon | |
112 | {1,1.16,2.74,0,0,0,0,0,0}, // phi from ALICE muon | |
113 | {1,0.72,2.5,0,0,0,0,0,0}}; // etaPrime from Pythia | |
c279681a | 114 | |
0d19baad | 115 | Double_t yparam7000[7][9] = {{1,0.8251,3.657,0,0,0,0,0,0}, // pions from pythia |
116 | {1,1.83,2.698,0,0,0,0,0,0}, // kaons from pythia | |
117 | {1,1.169,3.282,0,0,0,0,0,0}, // eta from pythia | |
118 | {1,1.234,3.264,0,0,0,0,0,0}, // rho from pythia | |
119 | {1,1.311,3.223,0,0,0,0,0,0}, // omega from pythia | |
120 | {1,2.388,2.129,0,0,0,0,0,0}, // phi from pythia | |
121 | {1,1.13,3.3,0,0,0,0,0,0}}; // eta prime from pythia | |
c279681a | 122 | |
0d19baad | 123 | // multiplicity parameters from pythia |
124 | Double_t nparam7000[7][9] = {{353.582, 6.76263, 1.66979, 998.445, 9.73281, 12.6704, 175.187, 29.08, 40.2531}, | |
125 | {1.e4, 0.2841, 0,0,0,0,0,0,0}, | |
126 | {1.e4, 0.2647, 0,0,0,0,0,0,0}, | |
127 | {7055, 0.1786, 0,0,0,0,0,0,0}, | |
128 | {7500, 0.1896, 0,0,0,0,0,0,0}, | |
129 | {5.e4, 1.167, 0,0,0,0,0,0,0}, | |
130 | {2.9e4, 0.714, 0,0,0,0,0,0,0}}; | |
c279681a | 131 | |
0d19baad | 132 | for (Int_t i=0; i<fgkNpart; i++) { |
133 | for (Int_t j=0; j<9; j++) { | |
134 | ptparam[i][j] = ptparam7000[i][j]; | |
135 | yparam[i][j] = yparam7000[i][j]; | |
136 | nparam[i][j] = nparam7000[i][j]; | |
137 | } | |
138 | } | |
139 | } | |
1f18dff1 | 140 | else if (fCMSEnergy==kCMS5020GeVpPb || fCMSEnergy==kCMS5020GeVPbp) { |
0d19baad | 141 | AliInfo ("Using pPb parameterization at 5.02 TeV\n"); |
142 | Double_t ptparam5020[7][9] = {{1,0.427,2.52,0,0,0,0,0,0}, // pions from Pythia at 7 TeV | |
143 | {1,0.58,2.57,0,0,0,0,0,0}, // kaons from Pythia at 7 TeV | |
144 | {1,0.665,2.796,0,0,0,0,0,0}, // eta from Pythia at 5.02 TeV | |
145 | {1,1.66,3.12,0,0,0,0,0,0}, // rho+omega from ALICE muon | |
146 | {1,1.66,3.12,0,0,0,0,0,0}, // rho+omega from ALICE muon | |
147 | {1,2.03,3.13,0,0,0,0,0,0}, // phi from ALICE muon | |
148 | {1,0.767,2.713,0,0,0,0,0,0}}; // etaPrime from Pythia at 5.02 TeV | |
c279681a | 149 | |
0d19baad | 150 | Double_t yparam5020[7][9] = {{1,0.8251,3.657,0,0,0,0,0,0}, // pions from pythia at 7 TeV |
151 | {1,1.83,2.698,0,0,0,0,0,0}, // kaons from pythia at 7 TeV | |
152 | {1,1.169,3.282,0,0,0,0,0,0}, // eta from pythia at 7 TeV | |
153 | {1,1.234,3.264,0,0,0,0,0,0}, // rho from pythia at 7 TeV | |
154 | {1,1.311,3.223,0,0,0,0,0,0}, // omega from pythia at 7 TeV | |
155 | {1,2.388,2.129,0,0,0,0,0,0}, // phi from pythia at 7 TeV | |
156 | {1,1.13,3.3,0,0,0,0,0,0}}; // eta prime from pythia at 7 TeV | |
c279681a | 157 | |
0d19baad | 158 | // multiplicity parameters from pythia at 7 TeV |
159 | Double_t nparam5020[7][9] = {{353.582, 6.76263, 1.66979, 998.445, 9.73281, 12.6704, 175.187, 29.08, 40.2531}, | |
160 | {1.e4, 0.2841, 0,0,0,0,0,0,0}, | |
161 | {1.e4, 0.2647, 0,0,0,0,0,0,0}, | |
162 | {7055, 0.1786, 0,0,0,0,0,0,0}, | |
163 | {7500, 0.1896, 0,0,0,0,0,0,0}, | |
164 | {5.e4, 1.167, 0,0,0,0,0,0,0}, | |
165 | {2.9e4, 0.714, 0,0,0,0,0,0,0}}; | |
c279681a | 166 | |
c279681a | 167 | for (Int_t i=0; i<fgkNpart; i++) { |
0d19baad | 168 | for (Int_t j=0; j<9; j++) { |
169 | ptparam[i][j] = ptparam5020[i][j]; | |
170 | yparam[i][j] = yparam5020[i][j]; | |
171 | nparam[i][j] = nparam5020[i][j]; | |
172 | } | |
173 | } | |
174 | if (fCMSEnergy==kCMS5020GeVpPb) fYCM = -0.4654; | |
175 | else fYCM = 0.4654; | |
176 | } | |
177 | else if (fCMSEnergy==kCMS2760GeV){ | |
178 | // parameters for 2.76 generation | |
179 | // pt params has been determined as <pt>ALICE_2.76 = <pt>ALICE_7 * <pt>PYTHIA_2.76 / <pt>PYTHIA_7 | |
180 | AliInfo ("Using pp parameterization at 2.76 TeV\n"); | |
181 | Double_t yparam2760[7][9] = {{1,0.8251,3.657,0,0,0,0,0,0}, // pions from pythia | |
182 | {1,1.83,2.698,0,0,0,0,0,0}, // kaons from pythia | |
183 | {1,0.011,3.474,0,0,0,0,0,0}, // eta from pythia | |
184 | {1,-0.01,3.409,0,0,0,0,0,0}, // rho from pythia | |
185 | {1,-0.037,3.294,0,0,0,0,0,0}, // omega from pythia | |
186 | {1,-0.016,2.717,0,0,0,0,0,0}, // phi from pythia | |
187 | {1,-0.010,3.312,0,0,0,0,0,0}}; // eta prime from pythia | |
c279681a | 188 | |
0d19baad | 189 | Double_t ptparam2760[7][9] = {{1,0.1665,8.878,0,0,0,0,0,0}, // pions from Pythia |
190 | {1,0.1657,8.591,0,0,0,0,0,0}, // kaons from Pythia | |
191 | {1,0.641,2.62,0,0,0,0,0,0}, // eta from ALICE 7 TeV | |
192 | {1,1.3551,3.16,0,0,0,0,0,0}, // rho with <pt> scaled | |
193 | {1,1.3551,3.16,0,0,0,0,0,0}, // omega with <pt> scaled | |
194 | {1,1.0811,2.74,0,0,0,0,0,0}, // phi with <pt> scaled | |
195 | {1,0.72,2.5,0,0,0,0,0,0}}; // etaPrime from ALICE 7 TeV | |
196 | ||
197 | Double_t nparam2760[7][9] = {{9752,-2.693,3.023,9.5e9,-84.68,16.75,-14.06,635.3,-423.2}, // pions | |
198 | {1.e5, 1.538, 0,0,0,0,0,0,0}, // kaons | |
199 | {1.e4, 0.351, 0,0,0,0,0,0,0}, // eta | |
200 | {1.e4, 0.2471, 0,0,0,0,0,0,0}, // rho | |
201 | {1.e4, 0.2583, 0,0,0,0,0,0,0}, // omega | |
202 | {1.e5, 1.393, 0,0,0,0,0,0,0}, // phi | |
203 | {1.e4, 0.9005, 0,0,0,0,0,0,0}}; // etaPrime | |
204 | ||
205 | for (Int_t i=0; i<fgkNpart; i++) { | |
206 | for (Int_t j=0; j<9; j++) { | |
207 | ptparam[i][j] = ptparam2760[i][j]; | |
208 | yparam[i][j] = yparam2760[i][j]; | |
209 | nparam[i][j] = nparam2760[i][j]; | |
210 | } | |
368fc61e | 211 | } |
0d19baad | 212 | } |
213 | else AliFatal("Energy not correctly defined"); | |
c279681a | 214 | |
0d19baad | 215 | for (Int_t i=0; i<fgkNpart; i++) { |
216 | fPDG[i] = pdg[i]; | |
217 | if (i!=0) { | |
218 | fMult[i] = new TF1(fnname[i],"[0]*exp(-[1]*x)",0,30); | |
219 | fMult[i]->SetParameters(nparam[i][0],nparam[i][1]); | |
220 | } | |
221 | else { | |
222 | fMult[i] = new TF1(fnname[i],"gaus(0)+gaus(3)+gaus(6)",0,150); | |
223 | for (Int_t j=0; j<9; j++) fMult[i]->SetParameter(j,nparam[i][j]); | |
224 | } | |
225 | ||
226 | fPt[i] = new TF1(fptname[i],AliGenMUONLMR::PtDistr,0,20,3); | |
227 | fPt[i]->SetParameters(ptparam[i][0], ptparam[i][1], ptparam[i][2]); | |
228 | fY[i] = new TF1(fyname[i],AliGenMUONLMR::YDistr,-10,10,4); | |
229 | fY[i]->SetParameters(yparam[i][0], yparam[i][1], yparam[i][2],fYCM); | |
230 | } | |
c279681a | 231 | |
0d19baad | 232 | for(Int_t i = 0; i<2; i++){ |
233 | fDecay[i] = new TF1(fdname[i],"exp(-x/[0])",0,150); | |
234 | fDecay[i]->SetParameter(0,ctau[i]); | |
235 | } | |
236 | ||
237 | for (Int_t ipart = 0; ipart < fgkNpart; ipart++) { | |
238 | fParticle[ipart] = new TParticle(); | |
239 | fParticle[ipart]->SetPdgCode(fPDG[ipart]); | |
240 | } | |
c279681a | 241 | |
0d19baad | 242 | TDatabasePDG *pdgdb = TDatabasePDG::Instance(); |
243 | Double_t mumass = pdgdb->GetParticle(13)->Mass(); | |
244 | fMu[0] = new TParticle(); | |
245 | fMu[0]->SetPdgCode(-13); | |
246 | fMu[0]->SetCalcMass(mumass); | |
247 | fMu[1] = new TParticle(); | |
248 | fMu[1]->SetPdgCode(13); | |
249 | fMu[1]->SetCalcMass(mumass); | |
c279681a | 250 | |
0d19baad | 251 | // function for polarized theta distributions |
252 | fCosTheta = new TF1 ("fCosTheta","1+[0]*x*x",-1,1); | |
253 | fCosTheta->SetParameter(0,1); | |
c279681a | 254 | |
0d19baad | 255 | // Dalitz decays |
256 | Int_t nbins = 1000; | |
257 | Double_t xmin = 0, xmax = 2; | |
258 | fDalitz[0] = new TH1F("hDalitzEta","",nbins,xmin,xmax); | |
259 | fDalitz[1] = new TH1F("hDalitzOmega","",nbins,xmin,xmax); | |
260 | fDalitz[2] = new TH1F("hDalitzEtaPrime","",nbins,xmin,xmax); | |
c279681a | 261 | |
0d19baad | 262 | Double_t meta = pdgdb->GetParticle("eta")->Mass(); |
263 | Double_t momega = pdgdb->GetParticle("omega")->Mass(); | |
264 | Double_t metaPrime = pdgdb->GetParticle("eta'")->Mass(); | |
265 | Double_t mpi0 = pdgdb->GetParticle("pi0")->Mass(); | |
266 | Double_t md3 = 0, mres = 0; | |
c279681a | 267 | |
0d19baad | 268 | for (Int_t index = 0; index < 3; index++) { |
269 | if (index == 0) { | |
270 | mres = meta; | |
271 | md3 = 0; | |
272 | } | |
273 | else if (index == 1) { | |
274 | mres = momega; | |
275 | md3 = mpi0; | |
276 | } | |
277 | else if (index == 2) { | |
278 | mres = metaPrime; | |
279 | md3 = 0; | |
280 | } | |
281 | Double_t delta = md3 * md3 / (mres * mres); | |
282 | Double_t epsilon = mumass * mumass / (mres * mres); | |
283 | nbins = fDalitz[index]->GetNbinsX(); | |
284 | xmin = fDalitz[index]->GetXaxis()->GetXmin(); | |
285 | Double_t deltax = fDalitz[index]->GetBinWidth(1); | |
286 | Double_t xd = xmin - deltax/2.; | |
287 | for (Int_t ibin = 0; ibin< nbins; ibin++) { | |
288 | Double_t dalval = 0; | |
289 | xd += deltax; | |
290 | if (xd > 4. *epsilon) { | |
291 | Double_t bracket = TMath::Power(1. + xd/(1. - delta),2) | |
292 | - 4. * xd / ((1. - delta) * (1. - delta)); | |
293 | if (bracket > 0) { | |
294 | dalval = TMath::Power(bracket,1.5) /xd * | |
295 | TMath::Sqrt(1 - 4 * epsilon / xd) * (1 + 2 * epsilon / xd) * | |
296 | FormFactor(xd * mres * mres, index); | |
297 | fDalitz[index]->Fill(xd,dalval); | |
c279681a | 298 | } |
0d19baad | 299 | } |
300 | } | |
301 | } | |
c279681a | 302 | |
0d19baad | 303 | fRhoLineShape = new TF1("fRhoLineShape",RhoLineShapeNew,0,2,2); |
304 | fHMultMu = new TH1D("fHMultMu","Muon multiplicity",20,-0.5,19.5); | |
305 | fHNProc = new TH1D("fHNProc","Number of gen. evts. per process in 4 pi",9,-0.5,8.5); | |
c279681a | 306 | |
368fc61e | 307 | } |
35e2c739 | 308 | //----------------------------------------------------------- |
309 | ||
310 | AliGenMUONLMR::AliGenMUONLMR (AliGenMUONLMR &gen) : AliGenMC(), | |
311 | fNMuMin(gen.fNMuMin), | |
0d19baad | 312 | fCMSEnergy(gen.fCMSEnergy), |
313 | fGenSingleProc(gen.fGenSingleProc), | |
314 | fYCM(gen.fYCM), | |
35e2c739 | 315 | fCosTheta(gen.fCosTheta), |
316 | fRhoLineShape(gen.fRhoLineShape), | |
317 | fHMultMu(gen.fHMultMu), | |
318 | fHNProc(gen.fHNProc) { | |
319 | for (Int_t i=0; i < fgkNpart; i++) { | |
320 | fPDG[i] = gen.fPDG[i]; | |
321 | fScaleMult[i] = gen.fScaleMult[i]; | |
322 | fPt[i] = (TF1*) gen.fPt[i]->Clone(); | |
323 | fY[i] = (TF1*) gen.fY[i]->Clone(); | |
324 | fMult[i] = (TF1*) gen.fMult[i]->Clone(); | |
325 | fParticle[i] = (TParticle*) gen.fParticle[i]->Clone(); | |
326 | } | |
327 | ||
328 | for(Int_t i = 0; i<2; i++) fDecay[i] = (TF1*) gen.fDecay[i]->Clone(); | |
329 | for(Int_t i = 0; i<3; i++) fDalitz[i] = (TH1F*) gen.fDalitz[i]->Clone(); | |
330 | for(Int_t i = 0; i<2; i++) fMu[i] = (TParticle*) gen.fMu[i]->Clone(); | |
331 | } | |
332 | ||
333 | //----------------------------------------------------------- | |
334 | ||
335 | AliGenMUONLMR& AliGenMUONLMR::operator=(const AliGenMUONLMR &gen) { | |
0d19baad | 336 | // Assignment operator |
a59583d3 | 337 | if (this!=&gen) { |
338 | fNMuMin = gen.fNMuMin; | |
0d19baad | 339 | fCMSEnergy = gen.fCMSEnergy; |
a59583d3 | 340 | fGenSingleProc = gen.fGenSingleProc; |
0d19baad | 341 | fYCM = gen.fYCM; |
a59583d3 | 342 | fCosTheta = (TF1*) gen.fCosTheta->Clone(); |
343 | fRhoLineShape = (TF1*) gen.fRhoLineShape->Clone(); | |
344 | fHMultMu = (TH1D*) gen.fHMultMu->Clone(); | |
345 | fHNProc = (TH1D*) gen.fHNProc->Clone(); | |
346 | ||
347 | for (Int_t i=0; i < fgkNpart; i++) { | |
348 | fPDG[i] = gen.fPDG[i]; | |
349 | fScaleMult[i] = gen.fScaleMult[i]; | |
350 | fPt[i] = (TF1*) gen.fPt[i]->Clone(); | |
351 | fY[i] = (TF1*) gen.fY[i]->Clone(); | |
352 | fMult[i] = (TF1*) gen.fMult[i]->Clone(); | |
353 | fParticle[i] = (TParticle*) gen.fParticle[i]->Clone(); | |
354 | } | |
355 | ||
356 | for(Int_t i = 0; i<2; i++) fDecay[i] = (TF1*) gen.fDecay[i]->Clone(); | |
357 | for(Int_t i = 0; i<3; i++) fDalitz[i] = (TH1F*) gen.fDalitz[i]->Clone(); | |
358 | for(Int_t i = 0; i<2; i++) fMu[i] = (TParticle*) gen.fMu[i]->Clone(); | |
35e2c739 | 359 | } |
35e2c739 | 360 | return *this; |
361 | } | |
362 | ||
363 | ||
368fc61e | 364 | //----------------------------------------------------------- |
365 | ||
366 | AliGenMUONLMR::~AliGenMUONLMR() | |
367 | { | |
368 | // Default destructor | |
369 | for (Int_t i=0; i<7; i++) { | |
370 | delete fPt[i]; | |
371 | delete fY[i]; | |
372 | delete fMult[i]; | |
373 | delete fParticle[i]; | |
374 | } | |
375 | ||
376 | for (Int_t i=0; i<2; i++) { | |
377 | delete fDecay[i]; | |
378 | delete fMu[i]; | |
379 | } | |
380 | ||
381 | for (Int_t i=0; i<3; i++) delete fDalitz[i]; | |
382 | ||
383 | delete fCosTheta; fCosTheta = 0; | |
384 | delete fRhoLineShape; fRhoLineShape = 0; | |
385 | delete fHMultMu; fHMultMu = 0; | |
386 | delete fHNProc; fHNProc = 0; | |
387 | } | |
388 | ||
389 | //----------------------------------------------------------- | |
390 | ||
391 | void AliGenMUONLMR::FinishRun(){ | |
392 | // save some histograms to an output file | |
393 | Int_t nbins = fHNProc->GetNbinsX(); | |
c279681a | 394 | for (Int_t ibin=1; ibin <= nbins; ibin++) AliInfo (Form("ibin = %d nEvProc = %g", |
0d19baad | 395 | ibin,fHNProc->GetBinContent(ibin))); |
396 | TFile *fout = new TFile("AliGenMUONLMR_histos.root","recreate"); | |
368fc61e | 397 | fHMultMu->Write(); |
398 | fHNProc->Write(); | |
399 | fout->Close(); | |
400 | } | |
401 | ||
402 | //----------------------------------------------------------- | |
403 | ||
35e2c739 | 404 | Double_t AliGenMUONLMR::YDistr(Double_t *px, Double_t *par){ |
368fc61e | 405 | // function for rapidity distribution: plateau at par[0] + |
406 | // gaussian tails centered at par[1] and with par[2]=sigma | |
0d19baad | 407 | Double_t ylab = px[0]; |
408 | Double_t y0 = par[3]; // center of mass rapidity | |
368fc61e | 409 | Double_t func = 0; |
0d19baad | 410 | if (ylab<y0+par[1] && ylab>y0-par[1]) func = par[0]; |
411 | else if (ylab>y0+par[1]) { | |
412 | Double_t z = (ylab-(par[1]+y0) )/(par[2]); | |
413 | func = par[0] * TMath::Exp(-0.5 * z * z); | |
414 | } | |
415 | else { | |
416 | Double_t z = (ylab-(-par[1]+y0) )/(par[2]); | |
368fc61e | 417 | func = par[0] * TMath::Exp(-0.5 * z * z); |
418 | } | |
419 | return func; | |
420 | } | |
421 | ||
422 | //----------------------------------------------------------- | |
423 | ||
35e2c739 | 424 | Double_t AliGenMUONLMR::PtDistr(Double_t *px, Double_t *par){ |
368fc61e | 425 | // pt distribution: power law |
426 | Double_t x = px[0]; | |
427 | Double_t func = par[0] * x / TMath::Power((1+(x/par[1])*(x/par[1])),par[2]); | |
428 | return func; | |
429 | } | |
430 | ||
431 | //----------------------------------------------------------- | |
432 | ||
433 | void AliGenMUONLMR::Generate() { | |
434 | // | |
435 | // generate the low mass resonances and their decays according to | |
436 | // the multiplicity parameterized by pythia and BR from PDG | |
437 | // rapidity distributions parametrized from pythia | |
438 | // pt distributions from data (or pythia for etaprime) | |
439 | // | |
440 | Double_t pxPushed[100], pyPushed[100], pzPushed[100], ePushed[100]; | |
441 | Int_t nmuons = -1, npartPushed = 0, pdgPushed[100]; | |
442 | Double_t polar[3]= {0,0,0}; // Polarisation of the parent particle (for GEANT tracking) | |
443 | Double_t origin0[3]; // Origin of the generated parent particle (for GEANT tracking) | |
444 | // Calculating vertex position per event | |
445 | for (Int_t j=0;j<3;j++) origin0[j]=fOrigin[j]; | |
446 | if(fVertexSmear==kPerEvent) { | |
447 | Vertex(); | |
448 | for (Int_t j=0;j<3;j++) origin0[j]=fVertex[j]; | |
449 | } | |
450 | ||
368fc61e | 451 | TParticle *mother; |
452 | TDatabasePDG* pdg = TDatabasePDG::Instance(); | |
453 | ||
454 | Double_t pt, y, phi, mass, px, py, pz, ene, mt; | |
455 | ||
456 | const Int_t nproc = 9; | |
457 | Int_t idRes[nproc] = {kEtaLMR, kEtaLMR, kRhoLMR, kOmegaLMR, kOmegaLMR, kPhiLMR, kEtaPrimeLMR, kPionLMR, kKaonLMR}; | |
35e2c739 | 458 | Double_t BR[nproc] = {5.8e-6, 3.1e-4, 4.55e-5, 7.28e-5, 1.3e-4, 2.86e-4, 1.04e-4, 1, 0.6344}; |
368fc61e | 459 | // Double_t BR[nproc] = {1, 1, 1, 1, 1, 1, 1, 1, 1}; |
460 | Int_t idDec[nproc] = {0, 1, 0, 0, 1, 0, 1, 2, 2}; // 0:2body, 1:Dalitz, 2:pi/K | |
461 | Int_t mult[nproc] = {0,0,0,0,0,0,0,0,0}; | |
462 | ||
463 | while (nmuons < fNMuMin) { | |
464 | ||
465 | nmuons = 0; | |
466 | npartPushed = 0; | |
467 | for (Int_t iproc=0; iproc<nproc; iproc++) { | |
468 | if (fGenSingleProc == -1) { | |
469 | mult[iproc] = Int_t(fMult[idRes[iproc]]->GetRandom()*fScaleMult[idRes[iproc]]); | |
470 | } | |
471 | else { | |
472 | if (iproc==fGenSingleProc) { | |
473 | mult[iproc] = 1; | |
474 | BR[iproc] = 1; | |
475 | } | |
476 | else { | |
477 | mult[iproc] = 0; | |
478 | BR[iproc] = 0; | |
479 | } | |
480 | } | |
481 | } | |
482 | ||
483 | if (fGenSingleProc == -1) { | |
484 | mult[1] = mult[0]; | |
485 | mult[4] = mult[3]; | |
486 | } | |
487 | ||
488 | for (Int_t iproc = 0; iproc < nproc; iproc++) { | |
0d19baad | 489 | // printf ("Multiplicity for process %d is %d\n",iproc,mult[iproc]); |
368fc61e | 490 | for (Int_t imult=0; imult<mult[iproc]; imult++) { |
491 | if (gRandom->Rndm() < BR[iproc]) { | |
492 | fHNProc->Fill(iproc); | |
493 | Int_t ipart = idRes[iproc]; | |
494 | pt = fPt[ipart]->GetRandom(); | |
495 | y = fY[ipart]->GetRandom(); | |
496 | phi = gRandom->Rndm() * 2 * TMath::Pi(); | |
497 | mass = pdg->GetParticle(fPDG[ipart])->Mass(); | |
498 | px = pt * TMath::Cos(phi); | |
499 | py = pt * TMath::Sin(phi); | |
500 | mt = TMath::Sqrt(pt * pt + mass * mass); | |
501 | pz = mt * TMath::SinH(y); | |
502 | ene = mt * TMath::CosH(y); | |
503 | ||
504 | mother = fParticle[ipart]; | |
505 | mother->SetMomentum(px,py,pz,ene); | |
506 | mother->SetCalcMass(mass); | |
507 | if (!KinematicSelection(mother,0)) continue; | |
508 | ||
509 | Bool_t hasDecayed = kTRUE; | |
510 | if (idDec[iproc] == 0) Decay2Body(mother); | |
511 | else if (idDec[iproc] == 1) DalitzDecay(mother); | |
512 | else DecayPiK(mother,hasDecayed); | |
513 | if (!hasDecayed) continue; | |
514 | Bool_t isMu0Acc = KinematicSelection(fMu[0],1); | |
515 | Bool_t isMu1Acc = KinematicSelection(fMu[1],1); | |
516 | Bool_t isMuFromPiKAcc = kTRUE; | |
517 | ||
518 | if (idDec[iproc] == 2) isMuFromPiKAcc = (mother->GetPdgCode()>0) ? isMu0Acc : isMu1Acc; | |
519 | // mother | |
520 | if ((idDec[iproc] < 2 && (isMu0Acc || isMu1Acc)) || | |
521 | (idDec[iproc] == 2 && isMuFromPiKAcc)) { | |
522 | pdgPushed[npartPushed] = mother->GetPdgCode(); | |
523 | pxPushed[npartPushed] = mother->Px(); | |
524 | pyPushed[npartPushed] = mother->Py(); | |
525 | pzPushed[npartPushed] = mother->Pz(); | |
526 | ePushed[npartPushed] = mother->Energy(); | |
527 | npartPushed++; | |
528 | if (isMu0Acc && (idDec[iproc] < 2 || mother->GetPdgCode() > 0)) { | |
529 | pdgPushed[npartPushed] = fMu[0]->GetPdgCode(); | |
530 | pxPushed[npartPushed] = fMu[0]->Px(); | |
531 | pyPushed[npartPushed] = fMu[0]->Py(); | |
532 | pzPushed[npartPushed] = fMu[0]->Pz(); | |
533 | ePushed[npartPushed] = fMu[0]->Energy(); | |
534 | npartPushed++; | |
535 | nmuons++; | |
536 | } | |
537 | ||
538 | if (isMu1Acc && (idDec[iproc] < 2 || mother->GetPdgCode() < 0)) { | |
539 | pdgPushed[npartPushed] = fMu[1]->GetPdgCode(); | |
540 | pxPushed[npartPushed] = fMu[1]->Px(); | |
541 | pyPushed[npartPushed] = fMu[1]->Py(); | |
542 | pzPushed[npartPushed] = fMu[1]->Pz(); | |
543 | ePushed[npartPushed] = fMu[1]->Energy(); | |
544 | npartPushed++; | |
545 | nmuons++; | |
546 | } | |
547 | } | |
548 | } // end if BR | |
549 | } // end loop on multiplicity | |
550 | } // end loop on process | |
551 | fHMultMu->Fill(nmuons); | |
552 | } // keep on generating until at least a muon is created in the event | |
553 | ||
554 | Int_t ntmother = 0, ntchild =0; | |
555 | for (Int_t ipart = 0; ipart < npartPushed; ipart++) { | |
556 | if (TMath::Abs(pdgPushed[ipart]) != 13) { // particle is not a muon, hence it's a mother | |
557 | PushTrack(0,-1,pdgPushed[ipart], | |
558 | pxPushed[ipart],pyPushed[ipart],pzPushed[ipart],ePushed[ipart], | |
35e2c739 | 559 | origin0[0],origin0[1],origin0[2],0., |
368fc61e | 560 | polar[0],polar[1],polar[2], |
561 | kPPrimary,ntmother,1,11); | |
562 | KeepTrack(ntmother); | |
563 | } | |
564 | else { | |
565 | PushTrack(1,ntmother,pdgPushed[ipart], | |
566 | pxPushed[ipart],pyPushed[ipart],pzPushed[ipart],ePushed[ipart], | |
35e2c739 | 567 | origin0[0],origin0[1],origin0[2],0., |
368fc61e | 568 | polar[0],polar[1],polar[2], |
569 | kPDecay,ntchild,1,1); | |
570 | KeepTrack(ntchild); | |
571 | } | |
572 | } | |
573 | SetHighWaterMark(ntchild); | |
574 | AliGenEventHeader* header = new AliGenEventHeader("LMR"); | |
575 | header->SetPrimaryVertex(fVertex); | |
576 | header->SetNProduced(fNprimaries); | |
577 | AddHeader(header); | |
578 | } | |
579 | ||
580 | //------------------------------------------------------------------ | |
581 | ||
35e2c739 | 582 | void AliGenMUONLMR::Decay2Body(TParticle *mother){ |
368fc61e | 583 | // performs decay in two muons of the low mass resonances |
584 | Double_t md1 = fMu[0]->GetMass(); | |
585 | Int_t pdg = mother->GetPdgCode(); | |
586 | Double_t mres =0; | |
587 | // if mother is a rho, extract the mass from its line shape | |
588 | // otherwise consider the resonance mass | |
589 | if (pdg == 113) mres = fRhoLineShape->GetRandom(); | |
590 | else mres = mother->GetCalcMass(); | |
591 | // while (mres < md1 + md2) mres = fDsigmaDm[res]->GetRandom(); | |
592 | // energies and momenta in rest frame | |
593 | Double_t e1 = mres / 2.; | |
594 | Double_t p1 = TMath::Sqrt((e1 + md1)*(e1 - md1)); | |
595 | // orientation in decaying particle rest frame | |
596 | Double_t costheta = gRandom->Rndm() * 2 - 1; | |
597 | Double_t sintheta = TMath::Sqrt((1. + costheta)*(1. - costheta)); | |
598 | Double_t phi = 2. * TMath::Pi() * gRandom->Rndm(); | |
599 | Double_t px1 = p1 * sintheta * TMath::Cos(phi); | |
600 | Double_t py1 = p1 * sintheta * TMath::Sin(phi); | |
601 | Double_t pz1 = p1 * costheta; | |
602 | ||
603 | // boost muons into lab frame | |
604 | ||
605 | TLorentzVector vmother, v1, v2; | |
606 | // TLorentzVector boosted1, boosted2; | |
607 | vmother.SetPxPyPzE(mother->Px(),mother->Py(),mother->Pz(),mother->Energy()); | |
608 | v1.SetPxPyPzE(px1,py1,pz1,e1); | |
609 | v2.SetPxPyPzE(-px1,-py1,-pz1,e1); | |
610 | ||
611 | TVector3 betaParent = (1./vmother.E())*vmother.Vect(); // beta = p/E | |
612 | v1.Boost(betaParent); | |
613 | v2.Boost(betaParent); | |
614 | ||
0d19baad | 615 | // TLorentzVector vtot = v1 + v2; |
616 | // printf ("mother: %g %g %g %g\n",vmother.Px(), vmother.Py(), vmother.Pz(), vmother.E()); | |
617 | // printf ("vtot : %g %g %g %g\n",vtot.Px(), vtot.Py(), vtot.Pz(), vtot.E()); | |
368fc61e | 618 | |
619 | fMu[0]->SetMomentum(v1.Px(),v1.Py(),v1.Pz(),v1.E()); | |
620 | fMu[1]->SetMomentum(v2.Px(),v2.Py(),v2.Pz(),v2.E()); | |
621 | } | |
622 | ||
623 | //------------------------------------------------------------------ | |
624 | ||
625 | void AliGenMUONLMR::DecayPiK(TParticle *mother, Bool_t &hasDecayed){ | |
626 | // performs decays of pions and kaons | |
627 | Double_t md1 = fMu[0]->GetMass(); | |
628 | // extract the mass from the resonance's line shape | |
629 | Double_t mres = mother->GetMass(); | |
630 | // choose the pi/k sign, assuming 50% probabilities for both signs | |
631 | Int_t sign = (gRandom->Rndm() > 0.5) ? 1 : -1; | |
632 | mother->SetPdgCode(sign * TMath::Abs(mother->GetPdgCode())); | |
633 | ||
634 | // energies and momenta in rest frame | |
635 | Double_t e1 = (mres*mres + md1*md1)/(2*mres); | |
636 | Double_t p1 = TMath::Sqrt((e1 + md1)*(e1 - md1)); | |
637 | // orientation in decaying particle rest frame | |
638 | Double_t costheta = gRandom->Rndm() * 2 - 1; | |
639 | Double_t sintheta = TMath::Sqrt((1. + costheta)*(1. - costheta)); | |
640 | Double_t phi = 2. * TMath::Pi() * gRandom->Rndm(); | |
641 | Double_t px1 = p1 * sintheta * TMath::Cos(phi); | |
642 | Double_t py1 = p1 * sintheta * TMath::Sin(phi); | |
643 | Double_t pz1 = p1 * costheta; | |
644 | ||
645 | // boost muons into lab frame | |
646 | TLorentzVector vmother, v1; | |
647 | vmother.SetPxPyPzE(mother->Px(),mother->Py(),mother->Pz(),mother->Energy()); | |
648 | v1.SetPxPyPzE(px1,py1,pz1,e1); | |
649 | ||
650 | TVector3 betaParent = (1./vmother.E())*vmother.Vect(); // beta = p/E | |
651 | v1.Boost(betaParent); | |
652 | if (mother->GetPdgCode()>0) fMu[0]->SetMomentum(v1.Px(),v1.Py(),v1.Pz(),v1.E()); | |
653 | else fMu[1]->SetMomentum(v1.Px(),v1.Py(),v1.Pz(),v1.E()); | |
654 | ||
a59583d3 | 655 | Int_t idmother = 0; |
35e2c739 | 656 | if (TMath::Abs(mother->GetPdgCode())== 211) idmother = 0; |
657 | if (TMath::Abs(mother->GetPdgCode())== 321) idmother = 1; | |
368fc61e | 658 | Double_t gammaRes = mother->Energy()/mres; |
659 | Double_t zResCM = fDecay[idmother]->GetRandom(); | |
660 | Double_t zResLab = gammaRes*zResCM; | |
661 | if(zResLab > 0.938) hasDecayed = 0; // 0.938: distance from IP to absorber + lambda_i | |
662 | else hasDecayed = 1; | |
368fc61e | 663 | } |
664 | ||
665 | //------------------------------------------------------------------- | |
666 | ||
35e2c739 | 667 | void AliGenMUONLMR::DalitzDecay(TParticle *mother){ |
368fc61e | 668 | // |
669 | // perform dalitz decays of eta, omega and etaprime | |
670 | // | |
671 | //in the rest frame of the virtual photon: | |
672 | Double_t mres = mother->GetCalcMass(); | |
673 | Double_t mumass = fMu[0]->GetMass(); | |
674 | Double_t md3 = 0; // unless differently specified, third particle is a photon | |
675 | if (mother->GetPdgCode() == 223) md3 = 0.134977; // if mother is an omega, third particle is a pi0 | |
a59583d3 | 676 | Int_t index = 0; |
35e2c739 | 677 | if (mother->GetPdgCode() == 221) index = 0; // eta |
678 | else if (mother->GetPdgCode() == 223) index = 1; // omega | |
679 | else if (mother->GetPdgCode() == 331) index = 2; // etaPrime | |
368fc61e | 680 | Int_t flag = 0; |
681 | Double_t xd=0, mvirt2=0; | |
682 | Double_t countIt = 0; | |
683 | while (flag==0) { | |
684 | xd = fDalitz[index]->GetRandom(); | |
685 | mvirt2 = xd * mres * mres; // mass of virtual photon | |
686 | // check kinematics | |
687 | if (mres - md3 > TMath::Sqrt(mvirt2) && TMath::Sqrt(mvirt2)/2. > mumass) flag=1; | |
688 | if (++countIt>1E11) { | |
689 | mvirt2 = mres * mres * 0.998; | |
690 | break; | |
691 | } | |
692 | } | |
693 | ||
694 | // | |
695 | // Generate muons in virtual photon rest frame. | |
696 | // z axis is the virt. photon direction (before boost) | |
697 | // | |
698 | ||
699 | Double_t e1 = TMath::Sqrt(mvirt2)/2.; // energy of mu1 in the virtual photon frame | |
700 | Double_t psquare = (e1 + mumass)*(e1 - mumass); | |
701 | if (psquare<0) { | |
c279681a | 702 | AliError(Form("sqrt of psquare = %f put to 0\n",psquare)); |
368fc61e | 703 | psquare = 0; |
704 | } | |
705 | Double_t p1 = TMath::Sqrt(psquare); | |
706 | //theta angle between the pos. muon and the virtual photon | |
707 | Double_t costheta = fCosTheta->GetRandom(); | |
708 | if (costheta>1) costheta = 1; | |
709 | if (costheta<-1) costheta = -1; | |
710 | Double_t sintheta = TMath::Sqrt((1. + costheta)*(1. - costheta)); | |
711 | Double_t phi = 2 * TMath::Pi() * gRandom->Rndm(); | |
712 | Double_t sinphi = TMath::Sin(phi); | |
713 | Double_t cosphi = TMath::Cos(phi); | |
714 | ||
715 | // fill 4-vectors of leptons in the virtual photon frame | |
716 | ||
717 | Double_t px1 = p1*sintheta*cosphi; | |
718 | Double_t py1 = p1*sintheta*sinphi; | |
719 | Double_t pz1 = p1*costheta; | |
720 | Double_t px2 = -p1*sintheta*cosphi; | |
721 | Double_t py2 = -p1*sintheta*sinphi; | |
722 | Double_t pz2 = -p1*costheta; | |
723 | Double_t e2 = e1; | |
724 | ||
725 | fMu[0]->SetMomentum(px1,py1,pz1,e1); | |
726 | fMu[1]->SetMomentum(px2,py2,pz2,e2); | |
727 | ||
728 | // calculate components of non-dilepton in CMS of parent resonance | |
729 | ||
730 | Double_t e3 = (mres * mres + md3 * md3 - mvirt2) / (2.*mres); | |
731 | Double_t psquare3 = (e3 + md3)*(e3 - md3); | |
732 | if (psquare3<0) { | |
c279681a | 733 | AliError(Form("Sqrt of psquare3 = %f put to 0\n",psquare3)); |
368fc61e | 734 | psquare3 = 0; |
735 | } | |
736 | Double_t p3 = TMath::Sqrt(psquare3); | |
737 | Double_t costheta2 = 2.* gRandom->Rndm() - 1.; // angle between virtual photon and resonance | |
738 | if (costheta2>1) costheta2 = 1; | |
739 | if (costheta2<-1) costheta2 = -1; | |
740 | Double_t sintheta2 = TMath::Sqrt((1. + costheta2)*(1. - costheta2)); | |
741 | Double_t phi2 = 2 * TMath::Pi() * gRandom->Rndm(); | |
742 | Double_t sinphi2 = TMath::Sin(phi2); | |
743 | Double_t cosphi2 = TMath::Cos(phi2); | |
744 | Double_t px3 = p3*sintheta2*cosphi2; | |
745 | Double_t py3 = p3*sintheta2*sinphi2; | |
746 | Double_t pz3 = p3*costheta2; | |
747 | TLorentzVector v3(px3,py3,pz3,e3); | |
748 | ||
749 | sintheta2 = -sintheta2; | |
750 | cosphi2 = -cosphi2; | |
751 | sinphi2 = -sinphi2; | |
752 | ||
753 | Double_t px1new = px1*costheta2*cosphi2 - py1*sinphi2 + pz1*sintheta2*cosphi2; | |
754 | Double_t py1new = px1*costheta2*sinphi2 + py1*cosphi2 + pz1*sintheta2*sinphi2; | |
755 | Double_t pz1new =-px1*sintheta2 + pz1*costheta2; | |
756 | Double_t px2new = px2*costheta2*cosphi2 - py2*sinphi2 + pz2*sintheta2*cosphi2; | |
757 | Double_t py2new = px2*costheta2*sinphi2 + py2*cosphi2 + pz2*sintheta2*sinphi2; | |
758 | Double_t pz2new =-px2*sintheta2 + pz2*costheta2; | |
759 | ||
760 | fMu[0]->SetMomentum(px1new,py1new,pz1new,e1); | |
761 | fMu[1]->SetMomentum(px2new,py2new,pz2new,e2); | |
762 | ||
763 | Double_t evirt = mres - e3; | |
764 | Double_t pxvirt = -px3; | |
765 | Double_t pyvirt = -py3; | |
766 | Double_t pzvirt = -pz3; | |
767 | TLorentzVector vvirt(pxvirt,pyvirt,pzvirt,evirt); | |
768 | TVector3 betaVirt = (1./evirt) * vvirt.Vect(); // virtual photon beta in res frame | |
769 | ||
770 | TLorentzVector v1(px1,py1,pz1,e1); | |
771 | TLorentzVector v2(px2,py2,pz2,e2); | |
772 | ||
773 | // boost the muons in the frame where the resonance is at rest | |
774 | ||
775 | v1.Boost(betaVirt); | |
776 | v2.Boost(betaVirt); | |
777 | ||
778 | // boost muons and third particle in lab frame | |
779 | ||
780 | TLorentzVector vmother(mother->Px(), mother->Py(), mother->Pz(), mother->Energy()); | |
781 | TVector3 resBetaLab = (1./vmother.E())*vmother.Vect(); // eta beta in lab frame | |
782 | v1.Boost(resBetaLab); | |
783 | v2.Boost(resBetaLab); | |
784 | v3.Boost(resBetaLab); | |
785 | vvirt.Boost(resBetaLab); | |
786 | ||
787 | fMu[0]->SetMomentum(v1.Px(),v1.Py(),v1.Pz(),v1.E()); | |
788 | fMu[1]->SetMomentum(v2.Px(),v2.Py(),v2.Pz(),v2.E()); | |
0d19baad | 789 | // part3->SetMomentum(v3.Px(),v3.Py(),v3.Pz(),v3.E()); |
368fc61e | 790 | |
0d19baad | 791 | // TLorentzVector vtot = v1 + v2 + v3; |
792 | // TLorentzVector vdimu = v1 + v2; | |
793 | // printf ("mother: %g %g %g %g\n",vmother.Px(), vmother.Py(), vmother.Pz(), vmother.E()); | |
794 | // printf ("vtot : %g %g %g %g\n",vtot.Px(), vtot.Py(), vtot.Pz(), vtot.E()); | |
795 | // printf ("vvirt : %g %g %g %g\n",vvirt.Px(), vvirt.Py(), vvirt.Pz(), vvirt.E()); | |
796 | // printf ("vdimu : %g %g %g %g\n",vdimu.Px(), vdimu.Py(), vdimu.Pz(), vdimu.E()); | |
368fc61e | 797 | |
798 | } | |
799 | ||
800 | //------------------------------------------------------------------ | |
801 | ||
802 | Double_t AliGenMUONLMR::FormFactor(Double_t q2, Int_t decay){ | |
803 | // Calculates the form factor for Dalitz decays A->B+l+l | |
804 | // Returns: |F(q^2)|^2 | |
805 | // | |
806 | // References: L.G. Landsberg, Physics Reports 128 No.6 (1985) 301-376. | |
807 | ||
808 | Double_t ff2, mass2; | |
809 | Double_t n2, n4, m2; | |
810 | // Lepton-G | |
811 | ||
812 | Double_t lambda2inv = 0; | |
813 | switch (decay) { | |
814 | case 0: // eta -> mu mu gamma | |
0d19baad | 815 | // eta -> l+ l- gamma: pole approximation |
368fc61e | 816 | lambda2inv = 1.95; |
817 | mass2 = fParticle[kEtaLMR]->GetMass() * fParticle[kEtaLMR]->GetMass(); | |
818 | if (q2 < mass2) ff2 = TMath::Power(1./(1.-lambda2inv*q2),2); | |
819 | else ff2 = 0; | |
820 | break; | |
821 | case 1: // omega -> mu mu pi0 | |
822 | // omega -> l+ l- pi0: pole approximation | |
823 | mass2 = fParticle[kOmegaLMR]->GetMass() * fParticle[kOmegaLMR]->GetMass(); | |
824 | lambda2inv = 2.26; | |
825 | if (q2 < mass2) ff2 = TMath::Power(1./(1.-lambda2inv*q2),2); | |
826 | else ff2 = 0; | |
827 | break; | |
828 | case 2: // etaPrime -> mu mu gamma | |
829 | mass2 = fParticle[kEtaPrimeLMR]->GetMass() * fParticle[kEtaPrimeLMR]->GetMass(); | |
830 | // eta' -> l+ l- gamma: Breit-Wigner fitted to data | |
831 | n2 = 0.764 * 0.764; | |
832 | n4 = n2 * n2; | |
833 | m2 = 0.1020 * 0.1020; | |
834 | if (q2 < mass2) ff2 = n4 / (TMath::Power(n2-q2,2) + m2 * n2); | |
835 | else ff2 = 0; | |
836 | break; | |
837 | default: | |
c279681a | 838 | AliError ("FormFactor: Decay not found"); |
368fc61e | 839 | return 0; |
840 | break; | |
841 | } | |
842 | return ff2; | |
843 | } | |
844 | ||
845 | //____________________________________________________________ | |
846 | ||
35e2c739 | 847 | Double_t AliGenMUONLMR::RhoLineShapeNew(Double_t *x, Double_t */*para*/){ |
368fc61e | 848 | //new parameterization implemented by Hiroyuki Sako (GSI) |
849 | Double_t mass = *x; | |
850 | double r, GammaTot; | |
851 | Double_t mRho = TDatabasePDG::Instance()->GetParticle("rho0")->Mass(); | |
852 | Double_t mPi = TDatabasePDG::Instance()->GetParticle("pi0")->Mass(); | |
853 | Double_t mMu = TDatabasePDG::Instance()->GetParticle("mu-")->Mass(); | |
854 | Double_t Gamma0 = TDatabasePDG::Instance()->GetParticle("rho0")->Width(); | |
855 | ||
856 | const double Norm = 0.0744416*1.01; | |
35e2c739 | 857 | |
368fc61e | 858 | // 0.0744416 at m = 0.72297 |
859 | // is the max number with Norm=1 (for rho) | |
860 | ||
861 | double mThreshold = 2.*mPi; | |
862 | ||
863 | const double T = 0.170; // Assumption of pi+ temperature [GeV/c^2] | |
864 | //const double T = 0.11; // Taken from fit to pi+ temperature [GeV/c^2] | |
865 | // with Reference: LEBC-EHS collab., Z. Phys. C 50 (1991) 405 | |
866 | ||
867 | if (mass < mThreshold) { | |
868 | r = 0.; | |
869 | return r; | |
870 | } | |
871 | ||
872 | double k = sqrt(0.25*mass*mass-(mThreshold/2)*(mThreshold/2)); | |
873 | double k0 = sqrt(0.25*mRho*mRho-(mThreshold/2)*(mThreshold/2)); | |
874 | ||
875 | GammaTot = (k/k0)*(k/k0)*(k/k0)*(mRho/mass)*(mRho/mass)*Gamma0; | |
876 | ||
877 | double FormFactor2 = 1/((mass*mass-mRho*mRho)*(mass*mass-mRho*mRho)+ | |
878 | mass*mass*GammaTot*GammaTot); | |
879 | ||
880 | r = pow(mass,1.5)*pow((1-mThreshold*mThreshold/(mass*mass)),1.5)* | |
881 | ((mass*mass+2*mMu*mMu)/(mass*mass))*(pow((mass*mass-4*mMu*mMu),0.5)/mass)*FormFactor2 | |
882 | *exp(-mass/T)/Norm; | |
883 | ||
884 | return r; | |
885 | } | |
c279681a | 886 |