3efc6d2f014f8376f8f859cadfa0fcb2390a3951
[u/mrichter/AliRoot.git] / FASTSIM / AliFastGlauber.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /*
17 $Log$
18 Revision 1.2  2003/04/14 14:23:44  morsch
19 Correction in Binaries().
20
21 Revision 1.1  2003/04/10 08:48:13  morsch
22 First commit.
23
24 */
25
26 // from AliRoot
27 #include "AliFastGlauber.h"
28 // from root
29 #include <TH1F.h>
30 #include <TF1.h>
31 #include <TF2.h>
32 #include <TCanvas.h>
33 #include <TRandom.h>
34 #include <TFile.h>
35
36 ClassImp(AliFastGlauber)
37
38 TF1*    AliFastGlauber::fWSb      = NULL;     
39 TF2*    AliFastGlauber::fWSbz     = NULL;    
40 TF1*    AliFastGlauber::fWSz      = NULL;     
41 TF1*    AliFastGlauber::fWSta     = NULL;    
42 TF2*    AliFastGlauber::fWStarfi  = NULL; 
43 TF1*    AliFastGlauber::fWStaa    = NULL;   
44 TF1*    AliFastGlauber::fWSgeo    = NULL;   
45 TF1*    AliFastGlauber::fWSbinary = NULL;   
46 TF1*    AliFastGlauber::fWSN      = NULL;   
47 Float_t AliFastGlauber::fbMax     = 0.;
48
49 AliFastGlauber::AliFastGlauber()
50 {
51 //  Default Constructor
52 //
53 //  Defaults for Pb
54     SetWoodSaxonParameters(6.624, 0.549, 0.00, 7.69e-4);
55     SetHardCrossSection();
56     SetMaxImpact();
57 }
58
59 void AliFastGlauber::Init(Int_t mode)
60 {
61 // Initialisation
62 //
63 //  Wood-Saxon
64 //
65     fWSb = new TF1("WSb", WSb, 0, fbMax, 4);
66     fWSb->SetParameter(0, fWSr0);
67     fWSb->SetParameter(1, fWSd);
68     fWSb->SetParameter(2, fWSw);
69     fWSb->SetParameter(3, fWSn);
70
71     fWSbz = new TF2("WSbz", WSbz, 0, fbMax, 4);
72     fWSbz->SetParameter(0, fWSr0);
73     fWSbz->SetParameter(1, fWSd);
74     fWSbz->SetParameter(2, fWSw);
75     fWSbz->SetParameter(3, fWSn);
76
77     fWSz = new TF1("WSz", WSz, 0, fbMax, 5);
78     fWSz->SetParameter(0, fWSr0);
79     fWSz->SetParameter(1, fWSd);
80     fWSz->SetParameter(2, fWSw);
81     fWSz->SetParameter(3, fWSn);
82
83 //
84 //  Thickness
85 //
86     fWSta = new TF1("WSta", WSta, 0., fbMax, 0);
87     
88 //
89 //  Overlap Kernel
90 //
91     fWStarfi = new TF2("WStarfi", WStarfi, 0., fbMax, 0., TMath::Pi(), 1);
92     fWStarfi->SetParameter(0, 0.);     
93     fWStarfi->SetNpx(200);     
94     fWStarfi->SetNpy(20);     
95 //
96 //  Overlap
97 //
98     if (! mode) {
99         fWStaa = new TF1("WStaa", WStaa, 0., fbMax, 0);
100         fWStaa->SetNpx(100);
101     } else {
102         TFile* f = new TFile("$(ALICE_ROOT)/FASTSIM/data/glauberPbPb.root");
103         fWStaa = (TF1*) f->Get("WStaa");
104     }
105     
106 //
107 //  Geometrical Cross-Section
108 //
109     fWSgeo = new TF1("WSgeo", WSgeo, 0., fbMax, 0);
110     fWSgeo->SetNpx(100);
111 //
112 //  Hard cross section (~ binary collisions)
113 //
114     fWSbinary = new TF1("WSbinary", WSbinary, 0., fbMax, 1);
115     fWSbinary->SetParameter(0, fSigmaHard); // mb
116     fWSbinary->SetNpx(100);
117 //
118 // Hard collisions per event
119 //
120     fWSN = new TF1("WSN", WSN, 0., fbMax, 1);
121     fWSN->SetNpx(100);
122 }
123
124 void AliFastGlauber::DrawWSb()
125 {
126 //
127 //  Draw Wood-Saxon Nuclear Density Function
128 //
129     TCanvas *c1 = new TCanvas("c1","Wood Saxon",400,10,600,700);
130     c1->cd();
131     fWSb->Draw();
132 }
133
134 void AliFastGlauber::DrawOverlap()
135 {
136 //
137 //  Draw Overlap Function
138 //
139     TCanvas *c2 = new TCanvas("c2","Overlap",400,10,600,700);
140     c2->cd();
141     fWStaa->Draw();
142 }
143
144 void AliFastGlauber::DrawThickness()
145 {
146 //
147 //  Draw Thickness Function
148 //
149     TCanvas *c3 = new TCanvas("c3","Thickness",400,10,600,700);
150     c3->cd();
151     fWSta->Draw();
152 }
153
154 void AliFastGlauber::DrawGeo()
155 {
156 //
157 //  Draw Geometrical Cross-Section
158 //
159     TCanvas *c3 = new TCanvas("c3","Geometrical Cross-Section",400,10,600,700);
160     c3->cd();
161     fWSgeo->Draw();
162 }
163
164 void AliFastGlauber::DrawBinary()
165 {
166 //
167 //  Draw Wood-Saxon Nuclear Density Function
168 //
169     TCanvas *c4 = new TCanvas("c4","Binary Cross-Section",400,10,600,700);
170     c4->cd();
171     fWSbinary->Draw();
172 }
173
174 void AliFastGlauber::DrawN()
175 {
176 //
177 //  Draw Binaries per event
178 //
179     TCanvas *c5 = new TCanvas("c5","Binaries per event",400,10,600,700);
180     c5->cd();
181     fWSN->Draw();
182 }
183
184 void AliFastGlauber::DrawKernel()
185 {
186 //
187 //  Draw Kernel
188 //
189     TCanvas *c6 = new TCanvas("c6","Kernel",400,10,600,700);
190     c6->cd();
191     fWStarfi->Draw();
192 }
193
194 Double_t AliFastGlauber::WSb(Double_t* x, Double_t* par)
195 {
196 //
197 //  Woods-Saxon Parameterisation
198 //  as a function of radius
199 //
200     Double_t xx  = x[0];
201     Double_t r0  = par[0];
202     Double_t d   = par[1];
203     Double_t w   = par[2];
204     Double_t n   = par[3];
205     
206     Double_t y  = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
207     return y;
208 }
209
210 Double_t AliFastGlauber::WSbz(Double_t* x, Double_t* par)
211 {
212 //
213 //  Wood Saxon Parameterisation
214 //  as a function of z and  b
215 //
216     Double_t bb  = x[0];
217     Double_t zz  = x[1];
218     Double_t r0  = par[0];
219     Double_t d   = par[1];
220     Double_t w   = par[2];
221     Double_t n   = par[3];
222     Double_t xx  = TMath::Sqrt(bb*bb+zz*zz);
223     Double_t y  = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
224     return y;
225 }
226
227 Double_t AliFastGlauber::WSz(Double_t* x, Double_t* par)
228 {
229 //
230 //  Wood Saxon Parameterisation
231 //  as a function of z for fixed b
232 //
233     Double_t bb  = par[4];
234     Double_t zz  = x[0];
235     Double_t r0  = par[0];
236     Double_t d   = par[1];
237     Double_t w   = par[2];
238     Double_t n   = par[3];
239     Double_t xx  = TMath::Sqrt(bb*bb+zz*zz);
240     Double_t y  = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
241 //    if (y < 1.e-6) y = 0;
242     return y;
243 }
244
245 Double_t AliFastGlauber::WSta(Double_t* x, Double_t* par)
246 {
247 //
248 //  Thickness function 
249 //
250     Double_t b  = x[0];
251     fWSz->SetParameter(4, b);
252     Double_t y  = 2. * fWSz->Integral(0., fbMax);
253     return y;
254 }
255
256
257
258 Double_t AliFastGlauber::WStarfi(Double_t* x, Double_t* par)
259 {
260 //
261 //  Kernel for overlap function
262 //
263     Double_t b    = par[0];
264     Double_t r1   = x[0];
265     Double_t phi  = x[1];
266     Double_t r2   = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi)); 
267     Double_t y    = r1 * fWSta->Eval(r1) * fWSta->Eval(r2);
268     return y;
269 }
270
271
272 Double_t AliFastGlauber::WStaa(Double_t* x, Double_t* par)
273 {
274 //
275 //  Overlap function
276 //
277     Double_t b    = x[0];
278     fWStarfi->SetParameter(0, b);
279 /*
280     Double_t al[2];
281     Double_t bl[2];
282     al[0] = 0.;
283     al[1] = 0.;
284     bl[0] = 6.6;
285     bl[1] = TMath::Pi();
286     Double_t err;
287     
288     Double_t y =  2. * fWStarfi->IntegralMultiple(2, al, bl, 0.001, err);
289     printf("WStaa: %f %f %f\n", b, y, err);
290 */
291 //
292 //  MC Integration
293 //
294     Double_t y = 0;
295     for (Int_t i = 0; i < 100000; i++)
296     {
297         Double_t phi = TMath::Pi() * gRandom->Rndm();
298         Double_t b1  = fbMax       * gRandom->Rndm();   
299         y += fWStarfi->Eval(b1, phi);
300     }
301     y *= 2. * 0.1 *  208. * 208. * TMath::Pi() * fbMax / 100000.;
302     return y;
303 }
304
305 Double_t AliFastGlauber::WSgeo(Double_t* x, Double_t* par)
306 {
307 //
308 //  Geometrical Cross-Section
309 //
310     Double_t b    = x[0];
311     Double_t taa  = fWStaa->Eval(b);
312     const Double_t sigma = 55.6; // mbarn
313     
314     Double_t y    = 2. * TMath::Pi() * b * (1. - TMath::Exp(- sigma * taa)); // fm
315     return y;
316 }
317
318
319 Double_t AliFastGlauber::WSbinary(Double_t* x, Double_t* par)
320 {
321 //
322 //  Geometrical Cross-Section
323 //
324     Double_t b     = x[0];
325     Double_t sigma = par[0];
326     Double_t taa   = fWStaa->Eval(b);
327     
328     Double_t y    = 2. * TMath::Pi() * b * sigma * taa; // fm
329     return y;
330 }
331
332 Double_t AliFastGlauber::WSN(Double_t* x, Double_t* par)
333 {
334 //
335 //  Geometrical Cross-Section
336 //
337     Double_t b     = x[0];
338     Double_t y;
339     if (b != 0.) {
340         y = fWSbinary->Eval(b)/fWSgeo->Eval(b);
341     } else {
342         y = fWSbinary->Eval(1.e-4)/fWSgeo->Eval(1.e-4);
343     }
344     return y;
345 }
346
347 void AliFastGlauber::SimulateTrigger(Int_t n)
348 {
349     //
350     //  Simulates Trigger
351     //
352     TH1F* mbtH = new TH1F("mbtH", "MB Trigger b-Distribution",   100, 0., 20.);
353     TH1F* hdtH = new TH1F("hdtH", "Hard Trigger b-Distribution", 100, 0., 20.);   
354     TH1F* mbmH = new TH1F("mbmH", "MB Trigger Multiplicity Distribution",   100, 0., 8000.);
355     TH1F* hdmH = new TH1F("hdmH", "Hard Trigger Multiplicity Distribution", 100, 0., 8000.);   
356
357     mbtH->SetXTitle("b [fm]");
358     hdtH->SetXTitle("b [fm]");    
359     mbmH->SetXTitle("Multiplicity");
360     hdmH->SetXTitle("Multiplicity");    
361
362     TCanvas *c0 = new TCanvas("c0","Trigger Simulation",400,10,600,700);    
363     c0->Divide(2,1);
364     TCanvas *c1 = new TCanvas("c1","Trigger Simulation",400,10,600,700);    
365     c1->Divide(1,2);
366
367     //
368     //
369     Init(1);
370     for (Int_t iev = 0; iev < n; iev++)
371     {
372         Float_t b, p, mult;
373         GetRandom(b, p, mult);
374         mbtH->Fill(b,1.);
375         hdtH->Fill(b, p);
376         mbmH->Fill(mult, 1.);
377         hdmH->Fill(mult, p);
378
379         c0->cd(1);
380         mbtH->Draw();
381         c0->cd(2);
382         hdtH->Draw();   
383         c0->Update();
384
385         c1->cd(1);
386         mbmH->Draw();
387         c1->cd(2);
388         hdmH->Draw();   
389         c1->Update();
390     }
391 }
392
393 void AliFastGlauber::GetRandom(Float_t& b, Float_t& p, Float_t& mult)
394 {
395     //
396     // Gives back a random impact parameter, hard trigger probability and multiplicity
397     //
398         b = fWSgeo->GetRandom();
399         Float_t mu = fWSN->Eval(b);
400         p = 1.-TMath::Exp(-mu);
401         mult = 6000./fWSN->Eval(0.) * mu;
402 }
403
404 Float_t  AliFastGlauber::GetRandomImpactParameter(Float_t bmin, Float_t bmax)
405 {
406     //
407     // Gives back a random impact parameter in the range bmin .. bmax
408     //
409
410     Float_t b = -1.;
411     while(b < bmin || b > bmax)
412         b = fWSgeo->GetRandom();
413     return b;
414 }
415
416 Double_t AliFastGlauber::CrossSection(Double_t b1, Double_t b2)
417 {
418     //
419     // Return cross-section integrated from b1 to b2 
420     //
421     
422     return fWSgeo->Integral(b1, b2)/100.;
423 }
424
425 Double_t AliFastGlauber::FractionOfHardCrossSection(Double_t b1, Double_t b2)
426 {
427     //
428     // Return raction of hard cross-section integrated from b1 to b2 
429     //
430     
431     return fWSbinary->Integral(b1, b2)/fWSbinary->Integral(0., 100.);
432 }
433
434
435 Double_t AliFastGlauber::Binaries(Double_t b)
436 {
437     //
438     // Return number of binary collisions normalized to 1 at b=0
439     //
440     
441     return fWSN->Eval(b)/fWSN->Eval(0.);
442 }