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