]> git.uio.no Git - u/mrichter/AliRoot.git/blame - EVGEN/AliGenMUONLMR.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Adding track references also in AliZDCv5
[u/mrichter/AliRoot.git] / EVGEN / AliGenMUONLMR.cxx
CommitLineData
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
12ClassImp(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 47void 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
310AliGenMUONLMR::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
335AliGenMUONLMR& 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
366AliGenMUONLMR::~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
391void 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 404Double_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 424Double_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
433void 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 582void 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
625void 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 667void 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
802Double_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 847Double_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