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