1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 // Utility class to make simple Glauber type calculations for collision geometries:
19 // Impact parameter, production points, reaction plane dependence
20 // The SimulateTrigger method can be used for simple MB and hard-process
21 // (binary scaling) trigger studies.
22 // Some basic quantities can be visualized directly.
23 // The default set-up for PbPb collisions can be read from a file calling Init(1).
26 // Author: andreas.morsch@cern.ch
29 #include "AliFastGlauber.h"
38 ClassImp(AliFastGlauber)
40 TF1* AliFastGlauber::fgWSb = NULL;
41 TF2* AliFastGlauber::fgWSbz = NULL;
42 TF1* AliFastGlauber::fgWSz = NULL;
43 TF1* AliFastGlauber::fgWSta = NULL;
44 TF2* AliFastGlauber::fgWStarfi = NULL;
45 TF2* AliFastGlauber::fgWAlmond = NULL;
46 TF1* AliFastGlauber::fgWStaa = NULL;
47 TF1* AliFastGlauber::fgWSgeo = NULL;
48 TF1* AliFastGlauber::fgWSbinary = NULL;
49 TF1* AliFastGlauber::fgWSN = NULL;
50 TF1* AliFastGlauber::fgWPathLength0 = NULL;
51 TF1* AliFastGlauber::fgWPathLength = NULL;
52 TF1* AliFastGlauber::fgWEnergyDensity = NULL;
53 TF1* AliFastGlauber::fgWIntRadius = NULL;
54 Float_t AliFastGlauber::fgBMax = 0.;
56 AliFastGlauber::AliFastGlauber()
58 // Default Constructor
61 SetWoodSaxonParameters(6.624, 0.549, 0.00, 7.69e-4);
62 SetHardCrossSection();
66 void AliFastGlauber::Init(Int_t mode)
72 fgWSb = new TF1("WSb", WSb, 0, fgBMax, 4);
73 fgWSb->SetParameter(0, fWSr0);
74 fgWSb->SetParameter(1, fWSd);
75 fgWSb->SetParameter(2, fWSw);
76 fgWSb->SetParameter(3, fWSn);
78 fgWSbz = new TF2("WSbz", WSbz, 0, fgBMax, 4);
79 fgWSbz->SetParameter(0, fWSr0);
80 fgWSbz->SetParameter(1, fWSd);
81 fgWSbz->SetParameter(2, fWSw);
82 fgWSbz->SetParameter(3, fWSn);
84 fgWSz = new TF1("WSz", WSz, 0, fgBMax, 5);
85 fgWSz->SetParameter(0, fWSr0);
86 fgWSz->SetParameter(1, fWSd);
87 fgWSz->SetParameter(2, fWSw);
88 fgWSz->SetParameter(3, fWSn);
93 fgWSta = new TF1("WSta", WSta, 0., fgBMax, 0);
98 fgWStarfi = new TF2("WStarfi", WStarfi, 0., fgBMax, 0., TMath::Pi(), 1);
99 fgWStarfi->SetParameter(0, 0.);
100 fgWStarfi->SetNpx(200);
101 fgWStarfi->SetNpy(20);
103 // Almond shaped interaction region
105 fgWAlmond = new TF2("WAlmond", WAlmond, -fgBMax, fgBMax, -fgBMax, fgBMax, 1);
106 fgWAlmond->SetParameter(0, 0.);
107 fgWAlmond->SetNpx(200);
108 fgWAlmond->SetNpy(200);
110 // Path Length as a function of Phi
112 fgWPathLength0 = new TF1("WPathLength0", WPathLength0, -TMath::Pi(), TMath::Pi(), 2);
113 fgWPathLength0->SetParameter(0, 0.);
114 // Pathlength definition
115 fgWPathLength0->SetParameter(1, 0.);
117 fgWPathLength = new TF1("WPathLength", WPathLength, -TMath::Pi(), TMath::Pi(), 3);
119 fgWPathLength->SetParameter(0, 0.);
120 // Number of interactions used for average
121 fgWPathLength->SetParameter(1, 1000.);
122 // Pathlength definition
123 fgWPathLength->SetParameter(2, 0);
125 fgWIntRadius = new TF1("WIntRadius", WIntRadius, 0., fgBMax, 1);
126 fgWIntRadius->SetParameter(0, 0.);
133 fgWStaa = new TF1("WStaa", WStaa, 0., fgBMax, 0);
134 fgWStaa->SetNpx(100);
136 TFile* f = new TFile("$(ALICE_ROOT)/FASTSIM/data/glauberPbPb.root");
137 fgWStaa = (TF1*) f->Get("WStaa");
141 fgWEnergyDensity = new TF1("WEnergyDensity", WEnergyDensity, 0., 2. * fWSr0, 1);
142 fgWEnergyDensity->SetParameter(0, fWSr0 + 1.);
145 // Geometrical Cross-Section
147 fgWSgeo = new TF1("WSgeo", WSgeo, 0., fgBMax, 0);
148 fgWSgeo->SetNpx(100);
150 // Hard cross section (~ binary collisions)
152 fgWSbinary = new TF1("WSbinary", WSbinary, 0., fgBMax, 1);
153 fgWSbinary->SetParameter(0, fSigmaHard); // mb
154 fgWSbinary->SetNpx(100);
156 // Hard collisions per event
158 fgWSN = new TF1("WSN", WSN, 0., fgBMax, 1);
162 void AliFastGlauber::DrawWSb()
165 // Draw Wood-Saxon Nuclear Density Function
167 TCanvas *c1 = new TCanvas("c1","Wood Saxon",400,10,600,700);
172 void AliFastGlauber::DrawOverlap()
175 // Draw Overlap Function
177 TCanvas *c2 = new TCanvas("c2","Overlap",400,10,600,700);
182 void AliFastGlauber::DrawThickness()
185 // Draw Thickness Function
187 TCanvas *c3 = new TCanvas("c3","Thickness",400,10,600,700);
192 void AliFastGlauber::DrawGeo()
195 // Draw Geometrical Cross-Section
197 TCanvas *c3 = new TCanvas("c3","Geometrical Cross-Section",400,10,600,700);
202 void AliFastGlauber::DrawBinary()
205 // Draw Wood-Saxon Nuclear Density Function
207 TCanvas *c4 = new TCanvas("c4","Binary Cross-Section",400,10,600,700);
212 void AliFastGlauber::DrawN()
215 // Draw Binaries per event
217 TCanvas *c5 = new TCanvas("c5","Binaries per event",400,10,600,700);
222 void AliFastGlauber::DrawKernel(Double_t b)
227 TCanvas *c6 = new TCanvas("c6","Kernel",400,10,600,700);
229 fgWStarfi->SetParameter(0, b);
233 void AliFastGlauber::DrawAlmond(Double_t b)
236 // Draw Interaction Almond
238 TCanvas *c7 = new TCanvas("c7","Almond",400,10,600,700);
240 fgWAlmond->SetParameter(0, b);
244 void AliFastGlauber::DrawPathLength0(Double_t b, Int_t iopt)
249 TCanvas *c8 = new TCanvas("c8","Path Length",400,10,600,700);
251 fgWPathLength0->SetParameter(0, b);
252 fgWPathLength0->SetParameter(1, Double_t(iopt));
253 fgWPathLength0->SetMinimum(0.);
254 fgWPathLength0->SetMaximum(10.);
255 fgWPathLength0->Draw();
258 void AliFastGlauber::DrawPathLength(Double_t b , Int_t ni, Int_t iopt)
263 TCanvas *c9 = new TCanvas("c9","Path Length",400,10,600,700);
265 fgWAlmond->SetParameter(0, b);
267 fgWPathLength->SetParameter(0, b);
268 fgWPathLength->SetParameter(1, Double_t (ni));
269 fgWPathLength->SetParameter(2, Double_t (iopt));
270 fgWPathLength->SetMinimum(0.);
271 fgWPathLength->SetMaximum(10.);
272 fgWPathLength->Draw();
275 void AliFastGlauber::DrawIntRadius(Double_t b)
278 // Draw Interaction Radius
280 TCanvas *c10 = new TCanvas("c10","Interaction Radius",400,10,600,700);
282 fgWIntRadius->SetParameter(0, b);
283 fgWIntRadius->SetMinimum(0.);
284 fgWIntRadius->Draw();
287 void AliFastGlauber::DrawEnergyDensity()
290 // Draw energy density
292 TCanvas *c11 = new TCanvas("c11","Energy Density",400, 10, 600, 700);
294 fgWEnergyDensity->SetMinimum(0.);
295 fgWEnergyDensity->Draw();
298 Double_t AliFastGlauber::WSb(Double_t* x, Double_t* par)
301 // Woods-Saxon Parameterisation
302 // as a function of radius
305 Double_t r0 = par[0];
310 Double_t y = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
315 Double_t AliFastGlauber::WSbz(Double_t* x, Double_t* par)
318 // Wood Saxon Parameterisation
319 // as a function of z and b
323 Double_t r0 = par[0];
327 Double_t xx = TMath::Sqrt(bb*bb+zz*zz);
328 Double_t y = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
333 Double_t AliFastGlauber::WSz(Double_t* x, Double_t* par)
336 // Wood Saxon Parameterisation
337 // as a function of z for fixed b
339 Double_t bb = par[4];
341 Double_t r0 = par[0];
345 Double_t xx = TMath::Sqrt(bb*bb+zz*zz);
346 Double_t y = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
351 Double_t AliFastGlauber::WSta(Double_t* x, Double_t* /*par*/)
354 // Thickness function
357 fgWSz->SetParameter(4, b);
358 Double_t y = 2. * fgWSz->Integral(0., fgBMax);
364 Double_t AliFastGlauber::WStarfi(Double_t* x, Double_t* par)
367 // Kernel for overlap function
372 Double_t r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
373 Double_t y = r1 * fgWSta->Eval(r1) * fgWSta->Eval(r2);
378 Double_t AliFastGlauber::WAlmond(Double_t* x, Double_t* par)
381 // Almond shaped interaction region
384 Double_t xx = x[0] + b/2.;
386 Double_t r1 = TMath::Sqrt(xx * xx + yy * yy);
387 Double_t phi = TMath::ATan2(yy,xx);
389 Double_t r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
391 // Interaction probability calculated as product of thicknesses
393 Double_t y = fgWSta->Eval(r1) * fgWSta->Eval(r2);
397 Double_t AliFastGlauber::WIntRadius(Double_t* x, Double_t* par)
400 // Average radius at which interaction takes place
402 // Radius in the Almond
406 fgWAlmond->SetParameter(0, b);
408 Double_t dphi = 2. * TMath::Pi() / 100.;
413 for (Int_t i = 0; i < 100; i++) {
414 Double_t xx = r * TMath::Cos(phi);
415 Double_t yy = r * TMath::Sin(phi);
416 y += fgWAlmond->Eval(xx,yy);
419 // Result multiplied by Jacobian (2 pi r)
420 return (2. * TMath::Pi() * y * r / 100.);
423 Double_t AliFastGlauber::WPathLength0(Double_t* x, Double_t* par)
426 // Path Length as a function of phi for interaction point fixed at (0,0)
430 const Int_t kNp = 100;
431 const Double_t kDr = fgBMax/Double_t(kNp);
434 // Path Length definition
435 Int_t iopt = Int_t(par[1]);
437 // Phi direction in Almond
438 Double_t phi0 = x[0];
442 // Step along radial direction phi
443 for (Int_t i = 0; i < kNp; i++) {
445 // Transform into target frame
447 Double_t xx = r * TMath::Cos(phi0) + b / 2.;
448 Double_t yy = r * TMath::Sin(phi0);
449 Double_t phi = TMath::ATan2(yy, xx);
451 Double_t r1 = TMath::Sqrt(xx * xx + yy * yy);
452 // Radius in projectile frame
453 Double_t r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
454 Double_t y = fgWSta->Eval(r1) * fgWSta->Eval(r2);
461 // My length definition (is exact for hard disk)
463 return (2. * rw / w);
465 return TMath::Sqrt(2. * rw * kDr / fgWSta->Eval(0.01) / fgWSta->Eval(0.01));
469 Double_t AliFastGlauber::WPathLength(Double_t* x, Double_t* par)
472 // Path Length as a function of phi
473 // Interaction point from random distribution
478 const Int_t kNp = 100;
479 const Double_t kDr = fgBMax/Double_t(kNp);
480 // Number of interactions
481 const Int_t kNpi = Int_t (par[1]);
486 // Path Length definition
487 Int_t iopt = Int_t(par[2]);
489 Double_t phi0 = x[0];
491 printf("phi0 %f \n", phi0);
496 for (Int_t in = 0; in < kNpi; in ++) {
502 fgWAlmond->GetRandom2(x0, y0);
504 Double_t r0 = TMath::Sqrt(x0 * x0 + y0 * y0);
505 Int_t nps = Int_t ((fgBMax - r0)/kDr) - 1;
509 for (Int_t i = 0; (i < nps ); i++) {
511 // Transform into target frame
512 Double_t xx = x0 + r * TMath::Cos(phi0) + b / 2.;
513 Double_t yy = y0 + r * TMath::Sin(phi0);
514 Double_t phi = TMath::ATan2(yy, xx);
515 Double_t r1 = TMath::Sqrt(xx * xx + yy * yy);
516 // Radius in projectile frame
517 Double_t r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
518 Double_t y = fgWSta->Eval(r1) * fgWSta->Eval(r2);
524 // Average over interactions
528 l+= 2. * rw * kDr / fgWSta->Eval(0.01) / fgWSta->Eval(0.01);
532 return (l / Double_t(kNpi));
534 return (TMath::Sqrt(l / Double_t(kNpi)));
537 Double_t AliFastGlauber::WStaa(Double_t* x, Double_t* /*par*/)
543 fgWStarfi->SetParameter(0, b);
553 Double_t y = 2. * fgWStarfi->IntegralMultiple(2, al, bl, 0.001, err);
554 printf("WStaa: %f %f %f\n", b, y, err);
560 for (Int_t i = 0; i < 100000; i++)
562 Double_t phi = TMath::Pi() * gRandom->Rndm();
563 Double_t b1 = fgBMax * gRandom->Rndm();
564 y += fgWStarfi->Eval(b1, phi);
566 y *= 2. * 0.1 * 208. * 208. * TMath::Pi() * fgBMax / 100000.;
570 Double_t AliFastGlauber::WSgeo(Double_t* x, Double_t* /*par*/)
573 // Geometrical Cross-Section
576 Double_t taa = fgWStaa->Eval(b);
577 const Double_t kSigma = 55.6; // mbarn
579 Double_t y = 2. * TMath::Pi() * b * (1. - TMath::Exp(- kSigma * taa)); // fm
584 Double_t AliFastGlauber::WSbinary(Double_t* x, Double_t* par)
587 // Number of binary collisions
590 Double_t sigma = par[0];
591 Double_t taa = fgWStaa->Eval(b);
593 Double_t y = 2. * TMath::Pi() * b * sigma * taa; // fm
597 Double_t AliFastGlauber::WSN(Double_t* x, Double_t* /*par*/)
600 // Number of hard processes per event
603 Double_t y = fgWSbinary->Eval(b)/fgWSgeo->Eval(b);
607 Double_t AliFastGlauber::WEnergyDensity(Double_t* x, Double_t* par)
610 // Initial energy density as a function of the impact parameter
613 Double_t rA = par[0];
615 // Attention: area of transverse reaction zone in hard-sphere approximation !
616 Double_t saa = (TMath::Pi() - 2. * TMath::ASin(b/ 2./ rA)) * rA * rA
617 - b * TMath::Sqrt(rA * rA - b * b/ 4.);
618 Double_t taa = fgWStaa->Eval(b);
623 void AliFastGlauber::SimulateTrigger(Int_t n)
628 TH1F* mbtH = new TH1F("mbtH", "MB Trigger b-Distribution", 100, 0., 20.);
629 TH1F* hdtH = new TH1F("hdtH", "Hard Trigger b-Distribution", 100, 0., 20.);
630 TH1F* mbmH = new TH1F("mbmH", "MB Trigger Multiplicity Distribution", 100, 0., 8000.);
631 TH1F* hdmH = new TH1F("hdmH", "Hard Trigger Multiplicity Distribution", 100, 0., 8000.);
633 mbtH->SetXTitle("b [fm]");
634 hdtH->SetXTitle("b [fm]");
635 mbmH->SetXTitle("Multiplicity");
636 hdmH->SetXTitle("Multiplicity");
638 TCanvas *c0 = new TCanvas("c0","Trigger Simulation",400,10,600,700);
640 TCanvas *c1 = new TCanvas("c1","Trigger Simulation",400,10,600,700);
646 for (Int_t iev = 0; iev < n; iev++)
649 GetRandom(b, p, mult);
652 mbmH->Fill(mult, 1.);
669 void AliFastGlauber::GetRandom(Float_t& b, Float_t& p, Float_t& mult)
672 // Gives back a random impact parameter, hard trigger probability and multiplicity
674 b = fgWSgeo->GetRandom();
675 Float_t mu = fgWSN->Eval(b);
676 p = 1.-TMath::Exp(-mu);
677 mult = 6000./fgWSN->Eval(1.) * mu;
680 void AliFastGlauber::GetRandom(Int_t& bin, Bool_t& hard)
683 // Gives back a random impact parameter bin, and hard trigger decission
685 Float_t b = fgWSgeo->GetRandom();
686 Float_t mu = fgWSN->Eval(b) * fSigmaHard;
687 Float_t p = 1.-TMath::Exp(-mu);
690 } else if (b < 8.6) {
692 } else if (b < 11.2) {
694 } else if (b < 13.2) {
696 } else if (b < 15.0) {
704 Float_t r = gRandom->Rndm();
706 if (r < p) hard = kTRUE;
710 Float_t AliFastGlauber::GetRandomImpactParameter(Float_t bmin, Float_t bmax)
713 // Gives back a random impact parameter in the range bmin .. bmax
717 while(b < bmin || b > bmax)
718 b = fgWSgeo->GetRandom();
722 Double_t AliFastGlauber::CrossSection(Double_t b1, Double_t b2)
725 // Return cross-section integrated from b1 to b2
728 return fgWSgeo->Integral(b1, b2)/100.;
731 Double_t AliFastGlauber::FractionOfHardCrossSection(Double_t b1, Double_t b2)
734 // Return raction of hard cross-section integrated from b1 to b2
737 return fgWSbinary->Integral(b1, b2)/fgWSbinary->Integral(0., 100.);
741 Double_t AliFastGlauber::Binaries(Double_t b)
744 // Return number of binary collisions normalized to 1 at b=0
747 return fgWSN->Eval(b)/fgWSN->Eval(0.001);