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