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