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