d046a3ce70dde3fcbc9c06e6e3684407feeb3425
[u/mrichter/AliRoot.git] / EVGEN / AliGenMUONLMR.cxx
1 //#include <TClonesArray.h> 
2
3 #include <TDatabasePDG.h>
4 #include <TFile.h>
5 #include "AliConst.h"
6 #include "AliLog.h"
7 #include "AliGenMUONLMR.h" 
8 #include "AliMC.h" 
9 #include "AliRun.h" 
10 #include "AliGenEventHeader.h"
11
12 ClassImp(AliGenMUONLMR)
13
14
15 AliGenMUONLMR::AliGenMUONLMR () : AliGenMC(), fNMuMin(2), fGenSingleProc(-1), fCosTheta(0), fRhoLineShape(0), fHMultMu(0), fHNProc(0) { 
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}; 
26   char* fptname[7] = {"fPtPion","fPtKaon","fPtEta","fPtRho","fPtOmega","fPtPhi","fPtEtaPrime"};
27   char* fyname[7] = {"fYPion","fYKaon","fYEta","fYRho","fYOmega","fYPhi","fYEtaPrime"}; 
28   char* fnname[7] = {"fMultPion","fMultKaon","fMultEta","fMultRho","fMultOmega","fMultPhi","fMultEtaPrime"};
29   char* fdname[2] = {"fDecPion","fDecKaon"};
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);
125     Int_t nbins0 = fDalitz[index]->GetNbinsX();
126     Double_t xmin0 = fDalitz[index]->GetXaxis()->GetXmin(); 
127     Double_t deltax =  fDalitz[index]->GetBinWidth(1);
128     Double_t xd = xmin0 - deltax/2.; 
129     for (Int_t ibin = 0; ibin< nbins0; ibin++) { 
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
152 AliGenMUONLMR::~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
177 void 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
190 Double_t AliGenMUONLMR::YDistr(const Double_t *px, const Double_t *par){ 
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
205 Double_t AliGenMUONLMR::PtDistr(const Double_t *px, const Double_t *par){
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
214 void 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)
225   // Calculating vertex position per event
226   for (Int_t j=0;j<3;j++) origin0[j]=fOrigin[j];
227   if(fVertexSmear==kPerEvent) {
228     Vertex();
229     for (Int_t j=0;j<3;j++) origin0[j]=fVertex[j];
230   }
231   
232   printf ("In Generate()\n");
233   TParticle *mother; 
234   TDatabasePDG* pdg = TDatabasePDG::Instance();
235
236   Double_t pt, y, phi, mass, px, py, pz, ene, mt; 
237
238   const Int_t nproc = 9; 
239   Int_t idRes[nproc] = {kEtaLMR, kEtaLMR, kRhoLMR, kOmegaLMR, kOmegaLMR, kPhiLMR, kEtaPrimeLMR, kPionLMR, kKaonLMR}; 
240   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};
241   //  Double_t BR[nproc] = {1, 1, 1, 1, 1, 1, 1, 1, 1};
242   Int_t idDec[nproc] = {0, 1, 0, 0, 1, 0, 1, 2, 2}; // 0:2body, 1:Dalitz, 2:pi/K 
243   Int_t mult[nproc] = {0,0,0,0,0,0,0,0,0}; 
244
245   while (nmuons < fNMuMin) { 
246
247     nmuons = 0; 
248     npartPushed = 0; 
249     for (Int_t iproc=0; iproc<nproc; iproc++) { 
250       if (fGenSingleProc == -1) { 
251         mult[iproc] = Int_t(fMult[idRes[iproc]]->GetRandom()*fScaleMult[idRes[iproc]]); 
252       }
253       else { 
254         if (iproc==fGenSingleProc) { 
255           mult[iproc] = 1; 
256           BR[iproc] = 1;
257         } 
258         else { 
259           mult[iproc] = 0; 
260           BR[iproc] = 0;
261         }
262       }
263     }
264     
265     if (fGenSingleProc == -1) { 
266       mult[1] = mult[0]; 
267       mult[4] = mult[3]; 
268     }
269     
270     for (Int_t iproc = 0; iproc < nproc; iproc++) { 
271 //       printf ("Multiplicity for process %d is %d\n",iproc,mult[iproc]); 
272       for (Int_t imult=0; imult<mult[iproc]; imult++) { 
273         if (gRandom->Rndm() < BR[iproc]) { 
274           fHNProc->Fill(iproc); 
275           Int_t ipart = idRes[iproc]; 
276           pt  = fPt[ipart]->GetRandom(); 
277           y   = fY[ipart]->GetRandom(); 
278           phi = gRandom->Rndm() * 2 * TMath::Pi(); 
279           mass = pdg->GetParticle(fPDG[ipart])->Mass(); 
280           px  = pt * TMath::Cos(phi); 
281           py  = pt * TMath::Sin(phi); 
282           mt  = TMath::Sqrt(pt * pt + mass * mass);
283           pz  = mt * TMath::SinH(y); 
284           ene = mt * TMath::CosH(y); 
285         
286           mother = fParticle[ipart]; 
287           mother->SetMomentum(px,py,pz,ene); 
288           mother->SetCalcMass(mass);
289           if (!KinematicSelection(mother,0)) continue; 
290
291           Bool_t hasDecayed = kTRUE;
292           if (idDec[iproc] == 0) Decay2Body(mother);
293           else if (idDec[iproc] == 1) DalitzDecay(mother); 
294           else DecayPiK(mother,hasDecayed); 
295           if (!hasDecayed) continue; 
296           Bool_t isMu0Acc = KinematicSelection(fMu[0],1); 
297           Bool_t isMu1Acc = KinematicSelection(fMu[1],1); 
298           Bool_t isMuFromPiKAcc = kTRUE; 
299
300           if (idDec[iproc] == 2) isMuFromPiKAcc = (mother->GetPdgCode()>0) ? isMu0Acc : isMu1Acc;
301           // mother 
302           if ((idDec[iproc]  < 2 && (isMu0Acc || isMu1Acc)) || 
303               (idDec[iproc] == 2 && isMuFromPiKAcc)) { 
304             pdgPushed[npartPushed] = mother->GetPdgCode(); 
305             pxPushed[npartPushed] = mother->Px(); 
306             pyPushed[npartPushed] = mother->Py(); 
307             pzPushed[npartPushed] = mother->Pz();
308             ePushed[npartPushed] = mother->Energy(); 
309             npartPushed++; 
310             if (isMu0Acc && (idDec[iproc] < 2 || mother->GetPdgCode() > 0)) { 
311               pdgPushed[npartPushed] = fMu[0]->GetPdgCode(); 
312               pxPushed[npartPushed] = fMu[0]->Px(); 
313               pyPushed[npartPushed] = fMu[0]->Py(); 
314               pzPushed[npartPushed] = fMu[0]->Pz();
315               ePushed[npartPushed] = fMu[0]->Energy(); 
316               npartPushed++; 
317               nmuons++; 
318             }
319             
320             if (isMu1Acc && (idDec[iproc] < 2 || mother->GetPdgCode() < 0)) { 
321               pdgPushed[npartPushed] = fMu[1]->GetPdgCode(); 
322               pxPushed[npartPushed] = fMu[1]->Px(); 
323               pyPushed[npartPushed] = fMu[1]->Py(); 
324               pzPushed[npartPushed] = fMu[1]->Pz();
325               ePushed[npartPushed] = fMu[1]->Energy(); 
326               npartPushed++; 
327               nmuons++; 
328             }
329           } 
330         } // end if BR
331       } // end loop on multiplicity 
332     }  // end loop on process 
333     fHMultMu->Fill(nmuons); 
334   } // keep on generating until at least a muon is created in the event
335   
336   Int_t ntmother = 0, ntchild =0; 
337   for (Int_t ipart = 0; ipart < npartPushed; ipart++) { 
338     if (TMath::Abs(pdgPushed[ipart]) != 13) { // particle is not a muon, hence it's a mother
339       PushTrack(0,-1,pdgPushed[ipart],
340                 pxPushed[ipart],pyPushed[ipart],pzPushed[ipart],ePushed[ipart],
341                 origin0[0],origin0[1],origin0[2],0.,
342                 polar[0],polar[1],polar[2],
343                 kPPrimary,ntmother,1,11);
344       KeepTrack(ntmother); 
345     }
346     else { 
347       PushTrack(1,ntmother,pdgPushed[ipart],
348                 pxPushed[ipart],pyPushed[ipart],pzPushed[ipart],ePushed[ipart],
349                 origin0[0],origin0[1],origin0[2],0.,
350                 polar[0],polar[1],polar[2],
351                 kPDecay,ntchild,1,1);
352       KeepTrack(ntchild); 
353     }
354   }
355   SetHighWaterMark(ntchild); 
356   AliGenEventHeader* header = new AliGenEventHeader("LMR");
357   header->SetPrimaryVertex(fVertex);
358   header->SetNProduced(fNprimaries);
359   AddHeader(header); 
360 }
361
362 //------------------------------------------------------------------
363
364 void AliGenMUONLMR::Decay2Body(const TParticle *mother){ 
365   // performs decay in two muons of the low mass resonances
366   Double_t md1 = fMu[0]->GetMass(); 
367   Int_t pdg = mother->GetPdgCode(); 
368   Double_t mres =0; 
369   // if mother is a rho, extract the mass from its line shape
370   // otherwise consider the resonance mass 
371   if (pdg == 113) mres = fRhoLineShape->GetRandom(); 
372   else mres = mother->GetCalcMass(); 
373   //  while (mres < md1 + md2) mres =  fDsigmaDm[res]->GetRandom();
374   // energies and momenta in rest frame 
375   Double_t e1 = mres / 2.;
376   Double_t p1 = TMath::Sqrt((e1 + md1)*(e1 - md1)); 
377   // orientation in decaying particle rest frame
378   Double_t costheta = gRandom->Rndm() * 2 - 1;
379   Double_t sintheta = TMath::Sqrt((1. + costheta)*(1. - costheta));
380   Double_t phi      = 2. * TMath::Pi() * gRandom->Rndm(); 
381   Double_t px1      = p1 * sintheta * TMath::Cos(phi); 
382   Double_t py1      = p1 * sintheta * TMath::Sin(phi); 
383   Double_t pz1      = p1 * costheta; 
384
385   // boost muons into lab frame 
386
387   TLorentzVector vmother, v1, v2;
388   //  TLorentzVector boosted1, boosted2;   
389   vmother.SetPxPyPzE(mother->Px(),mother->Py(),mother->Pz(),mother->Energy());
390   v1.SetPxPyPzE(px1,py1,pz1,e1); 
391   v2.SetPxPyPzE(-px1,-py1,-pz1,e1); 
392
393   TVector3 betaParent = (1./vmother.E())*vmother.Vect(); // beta = p/E
394   v1.Boost(betaParent);
395   v2.Boost(betaParent);
396
397 //   TLorentzVector vtot = v1 + v2; 
398 //   printf ("mother: %g   %g    %g     %g\n",vmother.Px(), vmother.Py(), vmother.Pz(), vmother.E());
399 //   printf ("vtot  : %g   %g    %g     %g\n",vtot.Px(), vtot.Py(), vtot.Pz(), vtot.E());
400
401   fMu[0]->SetMomentum(v1.Px(),v1.Py(),v1.Pz(),v1.E());
402   fMu[1]->SetMomentum(v2.Px(),v2.Py(),v2.Pz(),v2.E());
403
404
405 //------------------------------------------------------------------
406
407 void AliGenMUONLMR::DecayPiK(TParticle *mother, Bool_t &hasDecayed){ 
408   // performs decays of pions and kaons
409   Double_t md1 = fMu[0]->GetMass(); 
410   // extract the mass from the resonance's line shape
411   Double_t mres = mother->GetMass(); 
412   // choose the pi/k sign, assuming 50% probabilities for both signs
413   Int_t sign = (gRandom->Rndm() > 0.5) ? 1 : -1;
414   mother->SetPdgCode(sign * TMath::Abs(mother->GetPdgCode())); 
415
416   // energies and momenta in rest frame 
417   Double_t e1 = (mres*mres + md1*md1)/(2*mres);
418   Double_t p1 = TMath::Sqrt((e1 + md1)*(e1 - md1)); 
419   // orientation in decaying particle rest frame
420   Double_t costheta = gRandom->Rndm() * 2 - 1;
421   Double_t sintheta = TMath::Sqrt((1. + costheta)*(1. - costheta));
422   Double_t phi      = 2. * TMath::Pi() * gRandom->Rndm(); 
423   Double_t px1      = p1 * sintheta * TMath::Cos(phi); 
424   Double_t py1      = p1 * sintheta * TMath::Sin(phi); 
425   Double_t pz1      = p1 * costheta;  
426
427   // boost muons into lab frame 
428   TLorentzVector vmother, v1;
429   vmother.SetPxPyPzE(mother->Px(),mother->Py(),mother->Pz(),mother->Energy());
430   v1.SetPxPyPzE(px1,py1,pz1,e1); 
431
432   TVector3 betaParent = (1./vmother.E())*vmother.Vect(); // beta = p/E
433   v1.Boost(betaParent);  
434   if (mother->GetPdgCode()>0) fMu[0]->SetMomentum(v1.Px(),v1.Py(),v1.Pz(),v1.E());
435   else fMu[1]->SetMomentum(v1.Px(),v1.Py(),v1.Pz(),v1.E());
436
437   Int_t idmother = 0; 
438   if (TMath::Abs(mother->GetPdgCode())== 211) {
439       idmother = 0; 
440   } else if (TMath::Abs(mother->GetPdgCode())== 321) {
441       idmother = 1; 
442   } else {
443       AliWarning("Mother particle is not a pion or kaon \n");
444   }
445   
446   Double_t gammaRes = mother->Energy()/mres;
447   Double_t zResCM = fDecay[idmother]->GetRandom();
448   Double_t zResLab = gammaRes*zResCM;  
449   if(zResLab > 0.938) hasDecayed = 0; // 0.938: distance from IP to absorber + lambda_i
450   else hasDecayed = 1;
451
452
453
454 //-------------------------------------------------------------------
455
456 void AliGenMUONLMR::DalitzDecay(const TParticle *mother){
457   //
458   // perform dalitz decays of eta, omega and etaprime 
459   //
460   //in the rest frame of the virtual photon:
461   Double_t mres = mother->GetCalcMass(); 
462   Double_t mumass  = fMu[0]->GetMass(); 
463   Double_t md3  = 0;  // unless differently specified, third particle is a photon 
464   if (mother->GetPdgCode() == 223) md3 = 0.134977; // if mother is an omega, third particle is a pi0
465
466   Int_t index = 0; 
467   if (TMath::Abs(mother->GetPdgCode())== 221) {
468     // eta
469       index = 0; 
470   } else if (TMath::Abs(mother->GetPdgCode())== 223) {
471     // omega
472       index = 1; 
473   } else if (mother->GetPdgCode() == 331) {
474     // eta'
475       index = 2; 
476   } else {
477       AliWarning("Mother particle is not a eta, eta' or omega \n");
478   }
479
480   Int_t flag = 0; 
481   Double_t xd=0, mvirt2=0; 
482   Double_t countIt = 0;
483   while (flag==0) {  
484     xd       = fDalitz[index]->GetRandom(); 
485     mvirt2   = xd * mres * mres;   // mass of virtual photon 
486     // check kinematics 
487     if (mres - md3 > TMath::Sqrt(mvirt2) && TMath::Sqrt(mvirt2)/2. > mumass) flag=1;
488     if (++countIt>1E11) {
489       mvirt2 =  mres * mres * 0.998; 
490       break;
491     }
492   }  
493
494   //
495   //        Generate muons in virtual photon rest frame. 
496   //        z axis is the virt. photon direction (before boost)  
497   //
498
499   Double_t e1 = TMath::Sqrt(mvirt2)/2.; // energy of mu1 in the virtual photon frame
500   Double_t psquare = (e1 + mumass)*(e1 - mumass); 
501   if (psquare<0) {
502     printf("Error in AliGenMUONLMR::DalitzDecay: sqrt of psquare = %f put to 0\n",psquare); 
503     psquare = 0;
504   }
505   Double_t p1 = TMath::Sqrt(psquare);
506   //theta angle between the pos. muon and the virtual photon 
507   Double_t costheta = fCosTheta->GetRandom();
508   if (costheta>1)  costheta = 1; 
509   if (costheta<-1) costheta = -1; 
510   Double_t sintheta = TMath::Sqrt((1. + costheta)*(1. - costheta));
511   Double_t phi      = 2 * TMath::Pi() * gRandom->Rndm();
512   Double_t sinphi   = TMath::Sin(phi);
513   Double_t cosphi   = TMath::Cos(phi);
514
515   // fill 4-vectors of leptons in the virtual photon frame
516
517   Double_t px1 = p1*sintheta*cosphi; 
518   Double_t py1 = p1*sintheta*sinphi; 
519   Double_t pz1 = p1*costheta; 
520   Double_t px2 = -p1*sintheta*cosphi; 
521   Double_t py2 = -p1*sintheta*sinphi; 
522   Double_t pz2 = -p1*costheta; 
523   Double_t e2  = e1; 
524
525   fMu[0]->SetMomentum(px1,py1,pz1,e1); 
526   fMu[1]->SetMomentum(px2,py2,pz2,e2); 
527
528   // calculate components of non-dilepton in CMS of parent resonance 
529
530   Double_t e3 = (mres * mres + md3 * md3 - mvirt2) / (2.*mres);
531   Double_t psquare3 = (e3 + md3)*(e3 - md3); 
532   if (psquare3<0) {
533     printf("Error in AliGenMUONLMR::DalitzDecay: sqrt of psquare3 = %f put to 0\n",psquare3); 
534     psquare3 = 0;
535   }
536   Double_t p3 = TMath::Sqrt(psquare3);
537   Double_t costheta2 = 2.* gRandom->Rndm() - 1.;   // angle between virtual photon and resonance
538   if (costheta2>1)  costheta2 = 1; 
539   if (costheta2<-1) costheta2 = -1; 
540   Double_t sintheta2 = TMath::Sqrt((1. + costheta2)*(1. - costheta2));
541   Double_t phi2      = 2 * TMath::Pi() * gRandom->Rndm();
542   Double_t sinphi2   = TMath::Sin(phi2);
543   Double_t cosphi2   = TMath::Cos(phi2);
544   Double_t px3 = p3*sintheta2*cosphi2; 
545   Double_t py3 = p3*sintheta2*sinphi2; 
546   Double_t pz3 = p3*costheta2; 
547   TLorentzVector v3(px3,py3,pz3,e3); 
548
549   sintheta2 = -sintheta2;
550   cosphi2   = -cosphi2;
551   sinphi2   = -sinphi2;
552
553   Double_t px1new = px1*costheta2*cosphi2 - py1*sinphi2 + pz1*sintheta2*cosphi2; 
554   Double_t py1new = px1*costheta2*sinphi2 + py1*cosphi2 + pz1*sintheta2*sinphi2; 
555   Double_t pz1new =-px1*sintheta2                       + pz1*costheta2; 
556   Double_t px2new = px2*costheta2*cosphi2 - py2*sinphi2 + pz2*sintheta2*cosphi2; 
557   Double_t py2new = px2*costheta2*sinphi2 + py2*cosphi2 + pz2*sintheta2*sinphi2; 
558   Double_t pz2new =-px2*sintheta2                       + pz2*costheta2; 
559
560   fMu[0]->SetMomentum(px1new,py1new,pz1new,e1); 
561   fMu[1]->SetMomentum(px2new,py2new,pz2new,e2); 
562
563   Double_t evirt = mres - e3; 
564   Double_t pxvirt = -px3;
565   Double_t pyvirt = -py3;
566   Double_t pzvirt = -pz3;
567   TLorentzVector vvirt(pxvirt,pyvirt,pzvirt,evirt); 
568   TVector3 betaVirt = (1./evirt) * vvirt.Vect(); // virtual photon beta in res frame
569
570   TLorentzVector v1(px1,py1,pz1,e1); 
571   TLorentzVector v2(px2,py2,pz2,e2);
572
573   // boost the muons in the frame where the resonance is at rest 
574
575   v1.Boost(betaVirt); 
576   v2.Boost(betaVirt); 
577
578   // boost muons and third particle in lab frame
579
580   TLorentzVector vmother(mother->Px(), mother->Py(), mother->Pz(), mother->Energy());  
581   TVector3 resBetaLab = (1./vmother.E())*vmother.Vect(); // eta beta in lab frame
582   v1.Boost(resBetaLab); 
583   v2.Boost(resBetaLab); 
584   v3.Boost(resBetaLab); 
585   vvirt.Boost(resBetaLab); 
586
587   fMu[0]->SetMomentum(v1.Px(),v1.Py(),v1.Pz(),v1.E());
588   fMu[1]->SetMomentum(v2.Px(),v2.Py(),v2.Pz(),v2.E());
589 //   part3->SetMomentum(v3.Px(),v3.Py(),v3.Pz(),v3.E());
590
591 //   TLorentzVector vtot = v1 + v2 + v3; 
592 //   TLorentzVector vdimu = v1 + v2; 
593 //   printf ("mother: %g   %g    %g     %g\n",vmother.Px(), vmother.Py(), vmother.Pz(), vmother.E());
594 //   printf ("vtot  : %g   %g    %g     %g\n",vtot.Px(), vtot.Py(), vtot.Pz(), vtot.E());
595 //   printf ("vvirt : %g   %g    %g     %g\n",vvirt.Px(), vvirt.Py(), vvirt.Pz(), vvirt.E());
596 //   printf ("vdimu : %g   %g    %g     %g\n",vdimu.Px(), vdimu.Py(), vdimu.Pz(), vdimu.E());
597
598 }
599
600 //------------------------------------------------------------------
601
602 Double_t AliGenMUONLMR::FormFactor(Double_t q2, Int_t decay){ 
603   //  Calculates the form factor for Dalitz decays A->B+l+l
604   //  Returns: |F(q^2)|^2
605   //
606   //  References: L.G. Landsberg, Physics Reports 128 No.6 (1985) 301-376. 
607  
608   Double_t ff2, mass2;
609   Double_t n2, n4, m2; 
610   // Lepton-G
611   
612   Double_t lambda2inv = 0;
613   switch (decay) { 
614   case 0:   // eta -> mu mu gamma  
615   // eta   -> l+ l- gamma: pole approximation
616     lambda2inv = 1.95; 
617     mass2 = fParticle[kEtaLMR]->GetMass() * fParticle[kEtaLMR]->GetMass(); 
618     if (q2 < mass2) ff2 = TMath::Power(1./(1.-lambda2inv*q2),2);
619     else ff2 = 0; 
620     break;
621   case 1:   // omega -> mu mu pi0 
622     // omega -> l+ l- pi0: pole approximation
623     mass2 = fParticle[kOmegaLMR]->GetMass() * fParticle[kOmegaLMR]->GetMass(); 
624     lambda2inv = 2.26; 
625     if (q2 < mass2) ff2 = TMath::Power(1./(1.-lambda2inv*q2),2);
626     else ff2 = 0; 
627     break;
628   case 2:   // etaPrime -> mu mu gamma 
629     mass2 = fParticle[kEtaPrimeLMR]->GetMass() * fParticle[kEtaPrimeLMR]->GetMass(); 
630     // eta'  -> l+ l- gamma: Breit-Wigner fitted to data
631     n2 = 0.764 * 0.764; 
632     n4 = n2 * n2; 
633     m2 = 0.1020 * 0.1020;
634     if (q2 < mass2) ff2 = n4 / (TMath::Power(n2-q2,2) + m2 * n2); 
635     else ff2 = 0; 
636     break;
637   default:
638     printf ("FormFactor: Decay not found\n"); 
639     return 0; 
640     break; 
641   }
642   return ff2; 
643 }
644
645 //____________________________________________________________
646
647 Double_t AliGenMUONLMR::RhoLineShapeNew(const Double_t *x, const Double_t* /*para*/){
648   //new parameterization implemented by Hiroyuki Sako (GSI)
649   Double_t mass = *x;
650   double r, GammaTot;
651   Double_t mRho    = TDatabasePDG::Instance()->GetParticle("rho0")->Mass();
652   Double_t mPi     = TDatabasePDG::Instance()->GetParticle("pi0")->Mass();
653   Double_t mMu     = TDatabasePDG::Instance()->GetParticle("mu-")->Mass();
654   Double_t Gamma0  = TDatabasePDG::Instance()->GetParticle("rho0")->Width();
655
656   const double Norm = 0.0744416*1.01;  
657   // 0.0744416 at m = 0.72297
658   // is the max number with Norm=1 (for rho)
659   
660   double mThreshold = 2.*mPi;
661
662   const double T = 0.170; // Assumption of pi+ temperature [GeV/c^2]
663   //const double T = 0.11; // Taken from fit to pi+ temperature [GeV/c^2]
664   // with Reference: LEBC-EHS collab., Z. Phys. C 50 (1991) 405
665
666   if (mass < mThreshold) {
667     r = 0.;
668     return r;
669   }
670
671   double k = sqrt(0.25*mass*mass-(mThreshold/2)*(mThreshold/2));
672   double k0 = sqrt(0.25*mRho*mRho-(mThreshold/2)*(mThreshold/2));
673
674   GammaTot = (k/k0)*(k/k0)*(k/k0)*(mRho/mass)*(mRho/mass)*Gamma0;
675
676   double FormFactor2 = 1/((mass*mass-mRho*mRho)*(mass*mass-mRho*mRho)+
677                           mass*mass*GammaTot*GammaTot);
678
679   r = pow(mass,1.5)*pow((1-mThreshold*mThreshold/(mass*mass)),1.5)*
680     ((mass*mass+2*mMu*mMu)/(mass*mass))*(pow((mass*mass-4*mMu*mMu),0.5)/mass)*FormFactor2
681     *exp(-mass/T)/Norm;
682
683   return r;
684 }
685
686