]> git.uio.no Git - u/mrichter/AliRoot.git/blame - FASTSIM/AliFastGlauber.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Implementing SetTreeAddress
[u/mrichter/AliRoot.git] / FASTSIM / AliFastGlauber.cxx
CommitLineData
5b3a5a5d 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
88cb7938 16/* $Id$ */
5b3a5a5d 17
18// from AliRoot
19#include "AliFastGlauber.h"
20// from root
21#include <TH1F.h>
22#include <TF1.h>
23#include <TF2.h>
24#include <TCanvas.h>
25#include <TRandom.h>
26#include <TFile.h>
27
28ClassImp(AliFastGlauber)
29
f3a04204 30TF1* AliFastGlauber::fWSb = NULL;
31TF2* AliFastGlauber::fWSbz = NULL;
32TF1* AliFastGlauber::fWSz = NULL;
33TF1* AliFastGlauber::fWSta = NULL;
34TF2* AliFastGlauber::fWStarfi = NULL;
35TF2* AliFastGlauber::fWAlmond = NULL;
36TF1* AliFastGlauber::fWStaa = NULL;
37TF1* AliFastGlauber::fWSgeo = NULL;
38TF1* AliFastGlauber::fWSbinary = NULL;
39TF1* AliFastGlauber::fWSN = NULL;
40TF1* AliFastGlauber::fWPathLength0 = NULL;
41TF1* AliFastGlauber::fWPathLength = NULL;
42TF1* AliFastGlauber::fWIntRadius = NULL;
43Float_t AliFastGlauber::fbMax = 0.;
5b3a5a5d 44
45AliFastGlauber::AliFastGlauber()
46{
47// Default Constructor
48//
49// Defaults for Pb
50 SetWoodSaxonParameters(6.624, 0.549, 0.00, 7.69e-4);
51 SetHardCrossSection();
52 SetMaxImpact();
53}
54
55void AliFastGlauber::Init(Int_t mode)
56{
57// Initialisation
58//
59// Wood-Saxon
60//
61 fWSb = new TF1("WSb", WSb, 0, fbMax, 4);
62 fWSb->SetParameter(0, fWSr0);
63 fWSb->SetParameter(1, fWSd);
64 fWSb->SetParameter(2, fWSw);
65 fWSb->SetParameter(3, fWSn);
66
67 fWSbz = new TF2("WSbz", WSbz, 0, fbMax, 4);
68 fWSbz->SetParameter(0, fWSr0);
69 fWSbz->SetParameter(1, fWSd);
70 fWSbz->SetParameter(2, fWSw);
71 fWSbz->SetParameter(3, fWSn);
72
73 fWSz = new TF1("WSz", WSz, 0, fbMax, 5);
74 fWSz->SetParameter(0, fWSr0);
75 fWSz->SetParameter(1, fWSd);
76 fWSz->SetParameter(2, fWSw);
77 fWSz->SetParameter(3, fWSn);
78
79//
80// Thickness
81//
82 fWSta = new TF1("WSta", WSta, 0., fbMax, 0);
83
84//
85// Overlap Kernel
86//
87 fWStarfi = new TF2("WStarfi", WStarfi, 0., fbMax, 0., TMath::Pi(), 1);
88 fWStarfi->SetParameter(0, 0.);
89 fWStarfi->SetNpx(200);
90 fWStarfi->SetNpy(20);
91//
f3a04204 92// Almond shaped interaction region
93//
94 fWAlmond = new TF2("WAlmond", WAlmond, -fbMax, fbMax, -fbMax, fbMax, 1);
95 fWAlmond->SetParameter(0, 0.);
96 fWAlmond->SetNpx(200);
97 fWAlmond->SetNpy(200);
98//
99// Path Length as a function of Phi
100//
101 fWPathLength0 = new TF1("WPathLength0", WPathLength0, -TMath::Pi(), TMath::Pi(), 1);
102 fWPathLength0->SetParameter(0, 0.);
103 fWPathLength = new TF1("WPathLength", WPathLength, -TMath::Pi(), TMath::Pi(), 2);
104// Impact Parameter
105 fWPathLength->SetParameter(0, 0.);
106// Number of interactions used for average
107 fWPathLength->SetParameter(1, 1000.);
108
109 fWIntRadius = new TF1("WIntRadius", WIntRadius, 0., fbMax, 1);
110 fWIntRadius->SetParameter(0, 0.);
111//
5b3a5a5d 112// Overlap
113//
114 if (! mode) {
115 fWStaa = new TF1("WStaa", WStaa, 0., fbMax, 0);
116 fWStaa->SetNpx(100);
117 } else {
118 TFile* f = new TFile("$(ALICE_ROOT)/FASTSIM/data/glauberPbPb.root");
119 fWStaa = (TF1*) f->Get("WStaa");
120 }
121
122//
123// Geometrical Cross-Section
124//
125 fWSgeo = new TF1("WSgeo", WSgeo, 0., fbMax, 0);
126 fWSgeo->SetNpx(100);
127//
128// Hard cross section (~ binary collisions)
129//
130 fWSbinary = new TF1("WSbinary", WSbinary, 0., fbMax, 1);
131 fWSbinary->SetParameter(0, fSigmaHard); // mb
132 fWSbinary->SetNpx(100);
133//
134// Hard collisions per event
135//
136 fWSN = new TF1("WSN", WSN, 0., fbMax, 1);
137 fWSN->SetNpx(100);
138}
139
140void AliFastGlauber::DrawWSb()
141{
142//
143// Draw Wood-Saxon Nuclear Density Function
144//
145 TCanvas *c1 = new TCanvas("c1","Wood Saxon",400,10,600,700);
146 c1->cd();
147 fWSb->Draw();
148}
149
150void AliFastGlauber::DrawOverlap()
151{
152//
153// Draw Overlap Function
154//
155 TCanvas *c2 = new TCanvas("c2","Overlap",400,10,600,700);
156 c2->cd();
157 fWStaa->Draw();
158}
159
160void AliFastGlauber::DrawThickness()
161{
162//
163// Draw Thickness Function
164//
165 TCanvas *c3 = new TCanvas("c3","Thickness",400,10,600,700);
166 c3->cd();
167 fWSta->Draw();
168}
169
170void AliFastGlauber::DrawGeo()
171{
172//
173// Draw Geometrical Cross-Section
174//
175 TCanvas *c3 = new TCanvas("c3","Geometrical Cross-Section",400,10,600,700);
176 c3->cd();
177 fWSgeo->Draw();
178}
179
180void AliFastGlauber::DrawBinary()
181{
182//
183// Draw Wood-Saxon Nuclear Density Function
184//
185 TCanvas *c4 = new TCanvas("c4","Binary Cross-Section",400,10,600,700);
186 c4->cd();
187 fWSbinary->Draw();
188}
189
190void AliFastGlauber::DrawN()
191{
192//
193// Draw Binaries per event
194//
195 TCanvas *c5 = new TCanvas("c5","Binaries per event",400,10,600,700);
196 c5->cd();
197 fWSN->Draw();
198}
199
f3a04204 200void AliFastGlauber::DrawKernel(Double_t b)
5b3a5a5d 201{
202//
203// Draw Kernel
204//
205 TCanvas *c6 = new TCanvas("c6","Kernel",400,10,600,700);
206 c6->cd();
f3a04204 207 fWStarfi->SetParameter(0, b);
5b3a5a5d 208 fWStarfi->Draw();
209}
210
f3a04204 211void AliFastGlauber::DrawAlmond(Double_t b)
212{
213//
214// Draw Kernel
215//
216 TCanvas *c7 = new TCanvas("c7","Almond",400,10,600,700);
217 c7->cd();
218 fWAlmond->SetParameter(0, b);
219 fWAlmond->Draw();
220}
221
222void AliFastGlauber::DrawPathLength0(Double_t b)
223{
224//
225// Draw Kernel
226//
227 TCanvas *c8 = new TCanvas("c8","Path Length",400,10,600,700);
228 c8->cd();
229 fWPathLength0->SetParameter(0, b);
230 fWPathLength0->SetMinimum(0.);
231 fWPathLength0->SetMaximum(10.);
232 fWPathLength0->Draw();
233}
234
235void AliFastGlauber::DrawPathLength(Double_t b , Int_t ni)
236{
237//
238// Draw Kernel
239//
240 TCanvas *c9 = new TCanvas("c9","Path Length",400,10,600,700);
241 c9->cd();
242 fWAlmond->SetParameter(0, b);
243
244 fWPathLength->SetParameter(0, b);
245 fWPathLength->SetParameter(1, Double_t (ni));
246 fWPathLength->SetMinimum(0.);
247 fWPathLength->SetMaximum(10.);
248 fWPathLength->Draw();
249}
250
251void AliFastGlauber::DrawIntRadius(Double_t b)
252{
253//
254// Draw Kernel
255//
256 TCanvas *c10 = new TCanvas("c10","Interaction Radius",400,10,600,700);
257 c10->cd();
258 fWIntRadius->SetParameter(0, b);
259 fWIntRadius->SetMinimum(0.);
260 fWIntRadius->Draw();
261}
262
5b3a5a5d 263Double_t AliFastGlauber::WSb(Double_t* x, Double_t* par)
264{
265//
266// Woods-Saxon Parameterisation
267// as a function of radius
268//
269 Double_t xx = x[0];
270 Double_t r0 = par[0];
271 Double_t d = par[1];
272 Double_t w = par[2];
273 Double_t n = par[3];
274
275 Double_t y = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
276 return y;
277}
278
279Double_t AliFastGlauber::WSbz(Double_t* x, Double_t* par)
280{
281//
282// Wood Saxon Parameterisation
283// as a function of z and b
284//
285 Double_t bb = x[0];
286 Double_t zz = x[1];
287 Double_t r0 = par[0];
288 Double_t d = par[1];
289 Double_t w = par[2];
290 Double_t n = par[3];
291 Double_t xx = TMath::Sqrt(bb*bb+zz*zz);
292 Double_t y = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
293 return y;
294}
295
296Double_t AliFastGlauber::WSz(Double_t* x, Double_t* par)
297{
298//
299// Wood Saxon Parameterisation
300// as a function of z for fixed b
301//
302 Double_t bb = par[4];
303 Double_t zz = x[0];
304 Double_t r0 = par[0];
305 Double_t d = par[1];
306 Double_t w = par[2];
307 Double_t n = par[3];
308 Double_t xx = TMath::Sqrt(bb*bb+zz*zz);
309 Double_t y = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
310// if (y < 1.e-6) y = 0;
311 return y;
312}
313
f86dad79 314Double_t AliFastGlauber::WSta(Double_t* x, Double_t* /*par*/)
5b3a5a5d 315{
316//
317// Thickness function
318//
319 Double_t b = x[0];
320 fWSz->SetParameter(4, b);
321 Double_t y = 2. * fWSz->Integral(0., fbMax);
322 return y;
323}
324
325
326
327Double_t AliFastGlauber::WStarfi(Double_t* x, Double_t* par)
328{
329//
330// Kernel for overlap function
331//
332 Double_t b = par[0];
333 Double_t r1 = x[0];
334 Double_t phi = x[1];
335 Double_t r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
336 Double_t y = r1 * fWSta->Eval(r1) * fWSta->Eval(r2);
337 return y;
338}
339
340
f3a04204 341Double_t AliFastGlauber::WAlmond(Double_t* x, Double_t* par)
342{
343//
344// Almond shaped interaction region
345//
346 Double_t b = par[0];
347 Double_t xx = x[0] + b/2.;
348 Double_t yy = x[1];
349 Double_t r1 = TMath::Sqrt(xx * xx + yy * yy);
350 Double_t phi = TMath::ATan2(yy,xx);
351
352 Double_t r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
353//
354// Interaction probability calculated as product of thicknesses
355//
356 Double_t y = fWSta->Eval(r1) * fWSta->Eval(r2);
357 return y;
358}
359
360Double_t AliFastGlauber::WIntRadius(Double_t* x, Double_t* par)
361{
362//
363// Average radius at which interaction takes place
364//
365// Radius in the Almond
366 Double_t r = x[0];
367// Impact parameter
368 Double_t b = par[0];
369 fWAlmond->SetParameter(0, b);
370// Steps in phi
371 Double_t dphi = 2. * TMath::Pi() / 100.;
372// Average over phi
373 Double_t phi = 0.;
374 Double_t y = 0.;
375
376 for (Int_t i = 0; i < 100; i++) {
377 Double_t xx = r * TMath::Cos(phi);
378 Double_t yy = r * TMath::Sin(phi);
379 y += fWAlmond->Eval(xx,yy);
380 phi += dphi;
381 } // phi loop
382// Result multiplied by Jacobian (2 pi r)
383 return (2. * TMath::Pi() * y * r / 100.);
384}
385
386Double_t AliFastGlauber::WPathLength0(Double_t* x, Double_t* par)
387{
388//
389// Path Length as a function of phi for interaction point fixed at (0,0)
390//
391//
392// Steps in r
393 const Int_t np = 100;
394 const Double_t dr = fbMax/Double_t(np);
395// Impact parameter
396 Double_t b = par[0];
397// Phi direction in Almond
398 Double_t phi0 = x[0];
399 Double_t r = 0.;
400 Double_t rw = 0.;
401 Double_t w = 0.;
402// Step along radial direction phi
403 for (Int_t i = 0; i < np; i++) {
404//
405// Transform into target frame
406//
407 Double_t xx = r * TMath::Cos(phi0) + b / 2.;
408 Double_t yy = r * TMath::Sin(phi0);
409 Double_t phi = TMath::ATan2(yy, xx);
410
411 Double_t r1 = TMath::Sqrt(xx * xx + yy * yy);
412// Radius in projectile frame
413 Double_t r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
414 Double_t y = fWSta->Eval(r1) * fWSta->Eval(r2);
415
416 rw += y * r;
417 w += y;
418 r += dr;
419 } // radial steps
420//
421// My length definition (is exact for hard sphere)
422 return (2. * rw / w);
423}
424
425Double_t AliFastGlauber::WPathLength(Double_t* x, Double_t* par)
426{
427//
428// Path Length as a function of phi
429// Interaction point from random distribution
430//
431//
432// r-steps
433//
434 const Int_t np = 100;
435 const Double_t dr = fbMax/Double_t(np);
436// Number of interactions
437 const Int_t npi = Int_t (par[1]);
438
439//
440// Impact parameter
441 Double_t b = par[0];
442// Phi direction
443 Double_t phi0 = x[0];
444
445 printf("phi0 %f \n", phi0);
446
447// Path length
448 Double_t l = 0.;
449
450 for (Int_t in = 0; in < npi; in ++) {
451 Double_t rw = 0.;
452 Double_t w = 0.;
453
454 // Interaction point
455 Double_t x0, y0;
456 fWAlmond->GetRandom2(x0, y0);
457// Initial radius
458 Double_t r0 = TMath::Sqrt(x0 * x0 + y0 * y0);
459 Int_t nps = Int_t ((fbMax - r0)/dr) - 1;
460
461 Double_t r = 0.;
462// Radial steps
463 for (Int_t i = 0; (i < nps ); i++) {
464
465// Transform into target frame
466 Double_t xx = x0 + r * TMath::Cos(phi0) + b / 2.;
467 Double_t yy = y0 + r * TMath::Sin(phi0);
468 Double_t phi = TMath::ATan2(yy, xx);
469 Double_t r1 = TMath::Sqrt(xx * xx + yy * yy);
470// Radius in projectile frame
471 Double_t r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
472 Double_t y = fWSta->Eval(r1) * fWSta->Eval(r2);
473
474 rw += y * r;
475 w += y;
476 r += dr;
477 } // steps
478// Average over interactions
479 l += 2. * rw / w;
480 } // interactions
481 return (l / Double_t(npi));
482}
483
f86dad79 484Double_t AliFastGlauber::WStaa(Double_t* x, Double_t* /*par*/)
5b3a5a5d 485{
486//
487// Overlap function
488//
489 Double_t b = x[0];
490 fWStarfi->SetParameter(0, b);
491/*
492 Double_t al[2];
493 Double_t bl[2];
494 al[0] = 0.;
495 al[1] = 0.;
496 bl[0] = 6.6;
497 bl[1] = TMath::Pi();
498 Double_t err;
499
500 Double_t y = 2. * fWStarfi->IntegralMultiple(2, al, bl, 0.001, err);
501 printf("WStaa: %f %f %f\n", b, y, err);
502*/
503//
504// MC Integration
505//
506 Double_t y = 0;
507 for (Int_t i = 0; i < 100000; i++)
508 {
509 Double_t phi = TMath::Pi() * gRandom->Rndm();
510 Double_t b1 = fbMax * gRandom->Rndm();
511 y += fWStarfi->Eval(b1, phi);
512 }
513 y *= 2. * 0.1 * 208. * 208. * TMath::Pi() * fbMax / 100000.;
514 return y;
515}
516
f86dad79 517Double_t AliFastGlauber::WSgeo(Double_t* x, Double_t* /*par*/)
5b3a5a5d 518{
519//
520// Geometrical Cross-Section
521//
522 Double_t b = x[0];
523 Double_t taa = fWStaa->Eval(b);
524 const Double_t sigma = 55.6; // mbarn
525
526 Double_t y = 2. * TMath::Pi() * b * (1. - TMath::Exp(- sigma * taa)); // fm
527 return y;
528}
529
530
531Double_t AliFastGlauber::WSbinary(Double_t* x, Double_t* par)
532{
533//
c2140715 534// Number of binary collisions
5b3a5a5d 535//
536 Double_t b = x[0];
537 Double_t sigma = par[0];
538 Double_t taa = fWStaa->Eval(b);
539
540 Double_t y = 2. * TMath::Pi() * b * sigma * taa; // fm
541 return y;
542}
543
f86dad79 544Double_t AliFastGlauber::WSN(Double_t* x, Double_t* /*par*/)
5b3a5a5d 545{
546//
c2140715 547// Number of hard processes per event
5b3a5a5d 548//
549 Double_t b = x[0];
88cb7938 550 Double_t y = fWSbinary->Eval(b)/fWSgeo->Eval(b);
5b3a5a5d 551 return y;
552}
553
554void AliFastGlauber::SimulateTrigger(Int_t n)
555{
556 //
557 // Simulates Trigger
558 //
559 TH1F* mbtH = new TH1F("mbtH", "MB Trigger b-Distribution", 100, 0., 20.);
560 TH1F* hdtH = new TH1F("hdtH", "Hard Trigger b-Distribution", 100, 0., 20.);
561 TH1F* mbmH = new TH1F("mbmH", "MB Trigger Multiplicity Distribution", 100, 0., 8000.);
562 TH1F* hdmH = new TH1F("hdmH", "Hard Trigger Multiplicity Distribution", 100, 0., 8000.);
563
564 mbtH->SetXTitle("b [fm]");
565 hdtH->SetXTitle("b [fm]");
566 mbmH->SetXTitle("Multiplicity");
567 hdmH->SetXTitle("Multiplicity");
568
569 TCanvas *c0 = new TCanvas("c0","Trigger Simulation",400,10,600,700);
570 c0->Divide(2,1);
571 TCanvas *c1 = new TCanvas("c1","Trigger Simulation",400,10,600,700);
572 c1->Divide(1,2);
573
574 //
575 //
576 Init(1);
577 for (Int_t iev = 0; iev < n; iev++)
578 {
579 Float_t b, p, mult;
580 GetRandom(b, p, mult);
581 mbtH->Fill(b,1.);
582 hdtH->Fill(b, p);
583 mbmH->Fill(mult, 1.);
584 hdmH->Fill(mult, p);
585
586 c0->cd(1);
587 mbtH->Draw();
588 c0->cd(2);
589 hdtH->Draw();
590 c0->Update();
591
592 c1->cd(1);
593 mbmH->Draw();
594 c1->cd(2);
595 hdmH->Draw();
596 c1->Update();
597 }
598}
599
600void AliFastGlauber::GetRandom(Float_t& b, Float_t& p, Float_t& mult)
601{
602 //
603 // Gives back a random impact parameter, hard trigger probability and multiplicity
604 //
605 b = fWSgeo->GetRandom();
606 Float_t mu = fWSN->Eval(b);
607 p = 1.-TMath::Exp(-mu);
88cb7938 608 mult = 6000./fWSN->Eval(1.) * mu;
5b3a5a5d 609}
610
c2140715 611void AliFastGlauber::GetRandom(Int_t& bin, Bool_t& hard)
612{
613 //
614 // Gives back a random impact parameter bin, and hard trigger decission
615 //
616 Float_t b = fWSgeo->GetRandom();
617 Float_t mu = fWSN->Eval(b) * fSigmaHard;
618 Float_t p = 1.-TMath::Exp(-mu);
619 if (b < 5.) {
620 bin = 1;
621 } else if (b < 8.6) {
622 bin = 2;
623 } else if (b < 11.2) {
624 bin = 3;
625 } else if (b < 13.2) {
626 bin = 4;
204e011f 627 } else if (b < 15.0) {
c2140715 628 bin = 5;
204e011f 629 } else {
630 bin = 6;
c2140715 631 }
632
633 hard = kFALSE;
634
635 Float_t r = gRandom->Rndm();
636
637 if (r < p) hard = kTRUE;
638}
639
640
5b3a5a5d 641Float_t AliFastGlauber::GetRandomImpactParameter(Float_t bmin, Float_t bmax)
642{
643 //
644 // Gives back a random impact parameter in the range bmin .. bmax
645 //
646
647 Float_t b = -1.;
648 while(b < bmin || b > bmax)
649 b = fWSgeo->GetRandom();
650 return b;
651}
652
653Double_t AliFastGlauber::CrossSection(Double_t b1, Double_t b2)
654{
655 //
656 // Return cross-section integrated from b1 to b2
657 //
658
659 return fWSgeo->Integral(b1, b2)/100.;
660}
661
662Double_t AliFastGlauber::FractionOfHardCrossSection(Double_t b1, Double_t b2)
663{
664 //
665 // Return raction of hard cross-section integrated from b1 to b2
666 //
667
668 return fWSbinary->Integral(b1, b2)/fWSbinary->Integral(0., 100.);
669}
670
671
672Double_t AliFastGlauber::Binaries(Double_t b)
673{
674 //
675 // Return number of binary collisions normalized to 1 at b=0
676 //
677
270181a5 678 return fWSN->Eval(b)/fWSN->Eval(0.);
5b3a5a5d 679}