Added taa, participants and almonds for AuAu.
[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$ */
041f7f97 17//
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).
24//
25//
26// Author: andreas.morsch@cern.ch
5b3a5a5d 27
28// from AliRoot
29#include "AliFastGlauber.h"
30// from root
a2f2f511 31#include <TStyle.h>
5b3a5a5d 32#include <TH1F.h>
a2f2f511 33#include <TH2F.h>
5b3a5a5d 34#include <TF1.h>
35#include <TF2.h>
36#include <TCanvas.h>
37#include <TRandom.h>
38#include <TFile.h>
39
40ClassImp(AliFastGlauber)
41
041f7f97 42TF1* AliFastGlauber::fgWSb = NULL;
43TF2* AliFastGlauber::fgWSbz = NULL;
44TF1* AliFastGlauber::fgWSz = NULL;
45TF1* AliFastGlauber::fgWSta = NULL;
46TF2* AliFastGlauber::fgWStarfi = NULL;
47TF2* AliFastGlauber::fgWAlmond = NULL;
a2f2f511 48//TF2* AliFastGlauber::fWAlmondCurrent = NULL;
041f7f97 49TF1* AliFastGlauber::fgWStaa = NULL;
50TF1* AliFastGlauber::fgWSgeo = NULL;
51TF1* AliFastGlauber::fgWSbinary = NULL;
52TF1* AliFastGlauber::fgWSN = NULL;
53TF1* AliFastGlauber::fgWPathLength0 = NULL;
54TF1* AliFastGlauber::fgWPathLength = NULL;
55TF1* AliFastGlauber::fgWEnergyDensity = NULL;
56TF1* AliFastGlauber::fgWIntRadius = NULL;
1bc228f5 57TF2* AliFastGlauber::fgWKParticipants = NULL;
58TF1* AliFastGlauber::fgWParticipants = NULL;
041f7f97 59Float_t AliFastGlauber::fgBMax = 0.;
5b3a5a5d 60
61AliFastGlauber::AliFastGlauber()
62{
63// Default Constructor
64//
65// Defaults for Pb
66 SetWoodSaxonParameters(6.624, 0.549, 0.00, 7.69e-4);
67 SetHardCrossSection();
68 SetMaxImpact();
a2f2f511 69 SetLengthDefinition();
5b3a5a5d 70}
71
72void AliFastGlauber::Init(Int_t mode)
73{
74// Initialisation
a2f2f511 75// mode = 0; all functions are calculated
76// mode = 1; overlap function is read from file (for Pb-Pb only)
77// mode = 2; interaction almond functions are read from file
78// (for Pb-Pb only); USE THIS FOR PATH LENGTH CALC.!
79//
80
5b3a5a5d 81//
82// Wood-Saxon
83//
041f7f97 84 fgWSb = new TF1("WSb", WSb, 0, fgBMax, 4);
85 fgWSb->SetParameter(0, fWSr0);
86 fgWSb->SetParameter(1, fWSd);
87 fgWSb->SetParameter(2, fWSw);
88 fgWSb->SetParameter(3, fWSn);
5b3a5a5d 89
041f7f97 90 fgWSbz = new TF2("WSbz", WSbz, 0, fgBMax, 4);
91 fgWSbz->SetParameter(0, fWSr0);
92 fgWSbz->SetParameter(1, fWSd);
93 fgWSbz->SetParameter(2, fWSw);
94 fgWSbz->SetParameter(3, fWSn);
5b3a5a5d 95
041f7f97 96 fgWSz = new TF1("WSz", WSz, 0, fgBMax, 5);
97 fgWSz->SetParameter(0, fWSr0);
98 fgWSz->SetParameter(1, fWSd);
99 fgWSz->SetParameter(2, fWSw);
100 fgWSz->SetParameter(3, fWSn);
5b3a5a5d 101
102//
103// Thickness
104//
041f7f97 105 fgWSta = new TF1("WSta", WSta, 0., fgBMax, 0);
5b3a5a5d 106
107//
108// Overlap Kernel
109//
041f7f97 110 fgWStarfi = new TF2("WStarfi", WStarfi, 0., fgBMax, 0., TMath::Pi(), 1);
111 fgWStarfi->SetParameter(0, 0.);
112 fgWStarfi->SetNpx(200);
1bc228f5 113 fgWStarfi->SetNpy(20);
114//
115// Participants Kernel
116//
117 fgWKParticipants = new TF2("WKParticipants", WKParticipants, 0., fgBMax, 0., TMath::Pi(), 1);
118 fgWKParticipants->SetParameter(0, 0.);
119 fgWKParticipants->SetNpx(200);
120 fgWKParticipants->SetNpy(20);
5b3a5a5d 121//
f3a04204 122// Almond shaped interaction region
123//
041f7f97 124 fgWAlmond = new TF2("WAlmond", WAlmond, -fgBMax, fgBMax, -fgBMax, fgBMax, 1);
125 fgWAlmond->SetParameter(0, 0.);
126 fgWAlmond->SetNpx(200);
a2f2f511 127 fgWAlmond->SetNpy(200);
128
129 if(mode==2) {
130 printf(" Reading interaction almonds from file\n");
131 Char_t almondName[100];
132 TFile* ff = new TFile("$(ALICE_ROOT)/FASTSIM/data/glauberPbPb.root");
133 for(Int_t k=0; k<40; k++) {
134 sprintf(almondName,"WAlmondFixedB%d",k);
135 fWAlmondCurrent = (TF2*)ff->Get(almondName);
136 new(&fWAlmondFixedB[k]) TF2(*fWAlmondCurrent);
137 }
138 }
f3a04204 139//
140// Path Length as a function of Phi
141//
041f7f97 142 fgWPathLength0 = new TF1("WPathLength0", WPathLength0, -TMath::Pi(), TMath::Pi(), 2);
143 fgWPathLength0->SetParameter(0, 0.);
8de7e046 144// Pathlength definition
041f7f97 145 fgWPathLength0->SetParameter(1, 0.);
8de7e046 146
041f7f97 147 fgWPathLength = new TF1("WPathLength", WPathLength, -TMath::Pi(), TMath::Pi(), 3);
f3a04204 148// Impact Parameter
041f7f97 149 fgWPathLength->SetParameter(0, 0.);
f3a04204 150// Number of interactions used for average
041f7f97 151 fgWPathLength->SetParameter(1, 1000.);
8de7e046 152// Pathlength definition
041f7f97 153 fgWPathLength->SetParameter(2, 0);
f3a04204 154
041f7f97 155 fgWIntRadius = new TF1("WIntRadius", WIntRadius, 0., fgBMax, 1);
156 fgWIntRadius->SetParameter(0, 0.);
2a103154 157
158
f3a04204 159//
5b3a5a5d 160// Overlap
161//
162 if (! mode) {
041f7f97 163 fgWStaa = new TF1("WStaa", WStaa, 0., fgBMax, 0);
164 fgWStaa->SetNpx(100);
1bc228f5 165 fgWParticipants = new TF1("WParticipants", WParticipants, 0., fgBMax, 0);
166 fgWParticipants->SetNpx(100);
167
5b3a5a5d 168 } else {
169 TFile* f = new TFile("$(ALICE_ROOT)/FASTSIM/data/glauberPbPb.root");
041f7f97 170 fgWStaa = (TF1*) f->Get("WStaa");
1bc228f5 171 fgWParticipants = (TF1*) f->Get("WParticipants");
5b3a5a5d 172 }
173
174//
1bc228f5 175// Participants
176//
177
178//
041f7f97 179 fgWEnergyDensity = new TF1("WEnergyDensity", WEnergyDensity, 0., 2. * fWSr0, 1);
180 fgWEnergyDensity->SetParameter(0, fWSr0 + 1.);
2a103154 181
182//
5b3a5a5d 183// Geometrical Cross-Section
184//
041f7f97 185 fgWSgeo = new TF1("WSgeo", WSgeo, 0., fgBMax, 0);
186 fgWSgeo->SetNpx(100);
5b3a5a5d 187//
188// Hard cross section (~ binary collisions)
189//
041f7f97 190 fgWSbinary = new TF1("WSbinary", WSbinary, 0., fgBMax, 1);
191 fgWSbinary->SetParameter(0, fSigmaHard); // mb
192 fgWSbinary->SetNpx(100);
5b3a5a5d 193//
194// Hard collisions per event
195//
041f7f97 196 fgWSN = new TF1("WSN", WSN, 0., fgBMax, 1);
197 fgWSN->SetNpx(100);
5b3a5a5d 198}
199
200void AliFastGlauber::DrawWSb()
201{
202//
203// Draw Wood-Saxon Nuclear Density Function
204//
205 TCanvas *c1 = new TCanvas("c1","Wood Saxon",400,10,600,700);
206 c1->cd();
041f7f97 207 fgWSb->Draw();
5b3a5a5d 208}
209
210void AliFastGlauber::DrawOverlap()
211{
212//
213// Draw Overlap Function
214//
215 TCanvas *c2 = new TCanvas("c2","Overlap",400,10,600,700);
216 c2->cd();
041f7f97 217 fgWStaa->Draw();
5b3a5a5d 218}
219
1bc228f5 220void AliFastGlauber::DrawParticipants()
221{
222//
223// Draw Number of Participants
224//
225 TCanvas *c2 = new TCanvas("c2","Participants",400,10,600,700);
226 c2->cd();
227 fgWParticipants->Draw();
228}
229
5b3a5a5d 230void AliFastGlauber::DrawThickness()
231{
232//
233// Draw Thickness Function
234//
235 TCanvas *c3 = new TCanvas("c3","Thickness",400,10,600,700);
236 c3->cd();
041f7f97 237 fgWSta->Draw();
5b3a5a5d 238}
239
240void AliFastGlauber::DrawGeo()
241{
242//
243// Draw Geometrical Cross-Section
244//
245 TCanvas *c3 = new TCanvas("c3","Geometrical Cross-Section",400,10,600,700);
246 c3->cd();
041f7f97 247 fgWSgeo->Draw();
5b3a5a5d 248}
249
250void AliFastGlauber::DrawBinary()
251{
252//
a2f2f511 253// Draw Binary Cross-Section
5b3a5a5d 254//
255 TCanvas *c4 = new TCanvas("c4","Binary Cross-Section",400,10,600,700);
256 c4->cd();
041f7f97 257 fgWSbinary->Draw();
5b3a5a5d 258}
259
260void AliFastGlauber::DrawN()
261{
262//
263// Draw Binaries per event
264//
265 TCanvas *c5 = new TCanvas("c5","Binaries per event",400,10,600,700);
266 c5->cd();
041f7f97 267 fgWSN->Draw();
5b3a5a5d 268}
269
f3a04204 270void AliFastGlauber::DrawKernel(Double_t b)
5b3a5a5d 271{
272//
273// Draw Kernel
274//
275 TCanvas *c6 = new TCanvas("c6","Kernel",400,10,600,700);
276 c6->cd();
041f7f97 277 fgWStarfi->SetParameter(0, b);
278 fgWStarfi->Draw();
5b3a5a5d 279}
280
f3a04204 281void AliFastGlauber::DrawAlmond(Double_t b)
282{
283//
2a103154 284// Draw Interaction Almond
f3a04204 285//
286 TCanvas *c7 = new TCanvas("c7","Almond",400,10,600,700);
287 c7->cd();
041f7f97 288 fgWAlmond->SetParameter(0, b);
289 fgWAlmond->Draw();
f3a04204 290}
291
8de7e046 292void AliFastGlauber::DrawPathLength0(Double_t b, Int_t iopt)
f3a04204 293{
294//
2a103154 295// Draw Path Length
f3a04204 296//
297 TCanvas *c8 = new TCanvas("c8","Path Length",400,10,600,700);
298 c8->cd();
041f7f97 299 fgWPathLength0->SetParameter(0, b);
300 fgWPathLength0->SetParameter(1, Double_t(iopt));
301 fgWPathLength0->SetMinimum(0.);
302 fgWPathLength0->SetMaximum(10.);
303 fgWPathLength0->Draw();
f3a04204 304}
305
8de7e046 306void AliFastGlauber::DrawPathLength(Double_t b , Int_t ni, Int_t iopt)
f3a04204 307{
308//
2a103154 309// Draw Path Length
f3a04204 310//
311 TCanvas *c9 = new TCanvas("c9","Path Length",400,10,600,700);
312 c9->cd();
041f7f97 313 fgWAlmond->SetParameter(0, b);
f3a04204 314
041f7f97 315 fgWPathLength->SetParameter(0, b);
316 fgWPathLength->SetParameter(1, Double_t (ni));
317 fgWPathLength->SetParameter(2, Double_t (iopt));
318 fgWPathLength->SetMinimum(0.);
319 fgWPathLength->SetMaximum(10.);
320 fgWPathLength->Draw();
f3a04204 321}
322
323void AliFastGlauber::DrawIntRadius(Double_t b)
324{
325//
2a103154 326// Draw Interaction Radius
f3a04204 327//
328 TCanvas *c10 = new TCanvas("c10","Interaction Radius",400,10,600,700);
329 c10->cd();
041f7f97 330 fgWIntRadius->SetParameter(0, b);
331 fgWIntRadius->SetMinimum(0.);
332 fgWIntRadius->Draw();
f3a04204 333}
334
2a103154 335void AliFastGlauber::DrawEnergyDensity()
336{
337//
338// Draw energy density
339//
340 TCanvas *c11 = new TCanvas("c11","Energy Density",400, 10, 600, 700);
341 c11->cd();
041f7f97 342 fgWEnergyDensity->SetMinimum(0.);
343 fgWEnergyDensity->Draw();
2a103154 344}
345
5b3a5a5d 346Double_t AliFastGlauber::WSb(Double_t* x, Double_t* par)
347{
348//
349// Woods-Saxon Parameterisation
350// as a function of radius
351//
352 Double_t xx = x[0];
353 Double_t r0 = par[0];
354 Double_t d = par[1];
355 Double_t w = par[2];
356 Double_t n = par[3];
357
358 Double_t y = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
2a103154 359
5b3a5a5d 360 return y;
361}
362
363Double_t AliFastGlauber::WSbz(Double_t* x, Double_t* par)
364{
365//
366// Wood Saxon Parameterisation
367// as a function of z and b
368//
369 Double_t bb = x[0];
370 Double_t zz = x[1];
371 Double_t r0 = par[0];
372 Double_t d = par[1];
373 Double_t w = par[2];
374 Double_t n = par[3];
375 Double_t xx = TMath::Sqrt(bb*bb+zz*zz);
376 Double_t y = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
2a103154 377
5b3a5a5d 378 return y;
379}
380
381Double_t AliFastGlauber::WSz(Double_t* x, Double_t* par)
382{
383//
384// Wood Saxon Parameterisation
385// as a function of z for fixed b
386//
387 Double_t bb = par[4];
388 Double_t zz = x[0];
389 Double_t r0 = par[0];
390 Double_t d = par[1];
391 Double_t w = par[2];
392 Double_t n = par[3];
393 Double_t xx = TMath::Sqrt(bb*bb+zz*zz);
394 Double_t y = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
2a103154 395
5b3a5a5d 396 return y;
397}
398
f86dad79 399Double_t AliFastGlauber::WSta(Double_t* x, Double_t* /*par*/)
5b3a5a5d 400{
401//
402// Thickness function
403//
404 Double_t b = x[0];
041f7f97 405 fgWSz->SetParameter(4, b);
406 Double_t y = 2. * fgWSz->Integral(0., fgBMax);
5b3a5a5d 407 return y;
408}
409
410
411
412Double_t AliFastGlauber::WStarfi(Double_t* x, Double_t* par)
413{
414//
415// Kernel for overlap function
416//
417 Double_t b = par[0];
418 Double_t r1 = x[0];
419 Double_t phi = x[1];
420 Double_t r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
041f7f97 421 Double_t y = r1 * fgWSta->Eval(r1) * fgWSta->Eval(r2);
5b3a5a5d 422 return y;
423}
424
1bc228f5 425Double_t AliFastGlauber::WKParticipants(Double_t* x, Double_t* par)
426{
427//
428// Kernel for number of participants
429//
430 Double_t b = par[0];
431 Double_t r1 = x[0];
432 Double_t phi = x[1];
433 Double_t r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
434 Double_t xsi = fgWSta->Eval(r2) * 5.6;
435 Double_t a = 208;
436 Double_t sum = a * xsi;
437 Double_t y = sum;
438 for (Int_t i = 1; i <= 208; i++)
439 {
440 a--;
441 sum *= (-xsi) * a / Float_t(i+1);
442 y += sum;
443 }
444
445 y = r1 * fgWSta->Eval(r1) * y;
446 return y;
447}
448
5b3a5a5d 449
f3a04204 450Double_t AliFastGlauber::WAlmond(Double_t* x, Double_t* par)
451{
452//
453// Almond shaped interaction region
454//
455 Double_t b = par[0];
456 Double_t xx = x[0] + b/2.;
457 Double_t yy = x[1];
458 Double_t r1 = TMath::Sqrt(xx * xx + yy * yy);
459 Double_t phi = TMath::ATan2(yy,xx);
460
461 Double_t r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
462//
463// Interaction probability calculated as product of thicknesses
464//
041f7f97 465 Double_t y = fgWSta->Eval(r1) * fgWSta->Eval(r2);
f3a04204 466 return y;
467}
468
469Double_t AliFastGlauber::WIntRadius(Double_t* x, Double_t* par)
470{
471//
472// Average radius at which interaction takes place
473//
474// Radius in the Almond
475 Double_t r = x[0];
476// Impact parameter
477 Double_t b = par[0];
041f7f97 478 fgWAlmond->SetParameter(0, b);
f3a04204 479// Steps in phi
480 Double_t dphi = 2. * TMath::Pi() / 100.;
481// Average over phi
482 Double_t phi = 0.;
483 Double_t y = 0.;
484
485 for (Int_t i = 0; i < 100; i++) {
486 Double_t xx = r * TMath::Cos(phi);
487 Double_t yy = r * TMath::Sin(phi);
041f7f97 488 y += fgWAlmond->Eval(xx,yy);
f3a04204 489 phi += dphi;
490 } // phi loop
491// Result multiplied by Jacobian (2 pi r)
492 return (2. * TMath::Pi() * y * r / 100.);
493}
494
495Double_t AliFastGlauber::WPathLength0(Double_t* x, Double_t* par)
496{
497//
498// Path Length as a function of phi for interaction point fixed at (0,0)
499//
500//
501// Steps in r
041f7f97 502 const Int_t kNp = 100;
503 const Double_t kDr = fgBMax/Double_t(kNp);
f3a04204 504// Impact parameter
505 Double_t b = par[0];
8de7e046 506// Path Length definition
507 Int_t iopt = Int_t(par[1]);
508
f3a04204 509// Phi direction in Almond
510 Double_t phi0 = x[0];
511 Double_t r = 0.;
512 Double_t rw = 0.;
513 Double_t w = 0.;
514// Step along radial direction phi
041f7f97 515 for (Int_t i = 0; i < kNp; i++) {
f3a04204 516//
517// Transform into target frame
518//
519 Double_t xx = r * TMath::Cos(phi0) + b / 2.;
520 Double_t yy = r * TMath::Sin(phi0);
521 Double_t phi = TMath::ATan2(yy, xx);
522
523 Double_t r1 = TMath::Sqrt(xx * xx + yy * yy);
524// Radius in projectile frame
525 Double_t r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
041f7f97 526 Double_t y = fgWSta->Eval(r1) * fgWSta->Eval(r2);
f3a04204 527
528 rw += y * r;
529 w += y;
041f7f97 530 r += kDr;
f3a04204 531 } // radial steps
532//
b6d0ed6b 533// My length definition (is exact for hard disk)
8de7e046 534 if (!iopt) {
535 return (2. * rw / w);
536 } else {
041f7f97 537 return TMath::Sqrt(2. * rw * kDr / fgWSta->Eval(0.01) / fgWSta->Eval(0.01));
8de7e046 538 }
f3a04204 539}
540
541Double_t AliFastGlauber::WPathLength(Double_t* x, Double_t* par)
542{
543//
544// Path Length as a function of phi
545// Interaction point from random distribution
546//
547//
548// r-steps
549//
041f7f97 550 const Int_t kNp = 100;
551 const Double_t kDr = fgBMax/Double_t(kNp);
f3a04204 552// Number of interactions
041f7f97 553 const Int_t kNpi = Int_t (par[1]);
f3a04204 554
555//
556// Impact parameter
557 Double_t b = par[0];
8de7e046 558// Path Length definition
559 Int_t iopt = Int_t(par[2]);
f3a04204 560// Phi direction
561 Double_t phi0 = x[0];
562
563 printf("phi0 %f \n", phi0);
564
565// Path length
566 Double_t l = 0.;
567
041f7f97 568 for (Int_t in = 0; in < kNpi; in ++) {
f3a04204 569 Double_t rw = 0.;
570 Double_t w = 0.;
571
572 // Interaction point
573 Double_t x0, y0;
041f7f97 574 fgWAlmond->GetRandom2(x0, y0);
f3a04204 575// Initial radius
576 Double_t r0 = TMath::Sqrt(x0 * x0 + y0 * y0);
041f7f97 577 Int_t nps = Int_t ((fgBMax - r0)/kDr) - 1;
f3a04204 578
579 Double_t r = 0.;
580// Radial steps
581 for (Int_t i = 0; (i < nps ); i++) {
582
583// Transform into target frame
584 Double_t xx = x0 + r * TMath::Cos(phi0) + b / 2.;
585 Double_t yy = y0 + r * TMath::Sin(phi0);
586 Double_t phi = TMath::ATan2(yy, xx);
587 Double_t r1 = TMath::Sqrt(xx * xx + yy * yy);
588// Radius in projectile frame
589 Double_t r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
041f7f97 590 Double_t y = fgWSta->Eval(r1) * fgWSta->Eval(r2);
f3a04204 591
592 rw += y * r;
593 w += y;
041f7f97 594 r += kDr;
f3a04204 595 } // steps
596// Average over interactions
8de7e046 597 if (!iopt) {
598 l += (2. * rw / w);
599 } else {
041f7f97 600 l+= 2. * rw * kDr / fgWSta->Eval(0.01) / fgWSta->Eval(0.01);
8de7e046 601 }
f3a04204 602 } // interactions
9db56006 603 if (!iopt)
041f7f97 604 return (l / Double_t(kNpi));
9db56006 605 else
041f7f97 606 return (TMath::Sqrt(l / Double_t(kNpi)));
f3a04204 607}
608
f86dad79 609Double_t AliFastGlauber::WStaa(Double_t* x, Double_t* /*par*/)
5b3a5a5d 610{
611//
612// Overlap function
613//
614 Double_t b = x[0];
041f7f97 615 fgWStarfi->SetParameter(0, b);
5b3a5a5d 616/*
617 Double_t al[2];
618 Double_t bl[2];
619 al[0] = 0.;
620 al[1] = 0.;
621 bl[0] = 6.6;
622 bl[1] = TMath::Pi();
623 Double_t err;
624
041f7f97 625 Double_t y = 2. * fgWStarfi->IntegralMultiple(2, al, bl, 0.001, err);
5b3a5a5d 626 printf("WStaa: %f %f %f\n", b, y, err);
627*/
628//
629// MC Integration
630//
631 Double_t y = 0;
632 for (Int_t i = 0; i < 100000; i++)
633 {
634 Double_t phi = TMath::Pi() * gRandom->Rndm();
041f7f97 635 Double_t b1 = fgBMax * gRandom->Rndm();
636 y += fgWStarfi->Eval(b1, phi);
5b3a5a5d 637 }
041f7f97 638 y *= 2. * 0.1 * 208. * 208. * TMath::Pi() * fgBMax / 100000.;
5b3a5a5d 639 return y;
640}
641
1bc228f5 642Double_t AliFastGlauber::WParticipants(Double_t* x, Double_t* /*par*/)
643{
644//
645// Overlap function
646//
647 Double_t b = x[0];
648 fgWKParticipants->SetParameter(0, b);
649//
650// MC Integration
651//
652 Double_t y = 0;
653 for (Int_t i = 0; i < 100000; i++)
654 {
655 Double_t phi = TMath::Pi() * gRandom->Rndm();
656 Double_t b1 = fgBMax * gRandom->Rndm();
657 y += fgWKParticipants->Eval(b1, phi);
658 }
659 y *= 4. * 208. * TMath::Pi() * fgBMax / 100000.;
660 return y;
661}
662
663
f86dad79 664Double_t AliFastGlauber::WSgeo(Double_t* x, Double_t* /*par*/)
5b3a5a5d 665{
666//
667// Geometrical Cross-Section
668//
669 Double_t b = x[0];
041f7f97 670 Double_t taa = fgWStaa->Eval(b);
671 const Double_t kSigma = 55.6; // mbarn
5b3a5a5d 672
041f7f97 673 Double_t y = 2. * TMath::Pi() * b * (1. - TMath::Exp(- kSigma * taa)); // fm
5b3a5a5d 674 return y;
675}
676
677
678Double_t AliFastGlauber::WSbinary(Double_t* x, Double_t* par)
679{
680//
c2140715 681// Number of binary collisions
5b3a5a5d 682//
683 Double_t b = x[0];
684 Double_t sigma = par[0];
041f7f97 685 Double_t taa = fgWStaa->Eval(b);
5b3a5a5d 686
687 Double_t y = 2. * TMath::Pi() * b * sigma * taa; // fm
688 return y;
689}
690
f86dad79 691Double_t AliFastGlauber::WSN(Double_t* x, Double_t* /*par*/)
5b3a5a5d 692{
693//
c2140715 694// Number of hard processes per event
5b3a5a5d 695//
696 Double_t b = x[0];
041f7f97 697 Double_t y = fgWSbinary->Eval(b)/fgWSgeo->Eval(b);
5b3a5a5d 698 return y;
699}
700
2a103154 701Double_t AliFastGlauber::WEnergyDensity(Double_t* x, Double_t* par)
702{
703//
704// Initial energy density as a function of the impact parameter
705//
706 Double_t b = x[0];
707 Double_t rA = par[0];
b6d0ed6b 708//
709// Attention: area of transverse reaction zone in hard-sphere approximation !
710 Double_t saa = (TMath::Pi() - 2. * TMath::ASin(b/ 2./ rA)) * rA * rA
711 - b * TMath::Sqrt(rA * rA - b * b/ 4.);
041f7f97 712 Double_t taa = fgWStaa->Eval(b);
2a103154 713
714 return (taa/saa);
715}
716
5b3a5a5d 717void AliFastGlauber::SimulateTrigger(Int_t n)
718{
719 //
720 // Simulates Trigger
721 //
722 TH1F* mbtH = new TH1F("mbtH", "MB Trigger b-Distribution", 100, 0., 20.);
723 TH1F* hdtH = new TH1F("hdtH", "Hard Trigger b-Distribution", 100, 0., 20.);
724 TH1F* mbmH = new TH1F("mbmH", "MB Trigger Multiplicity Distribution", 100, 0., 8000.);
725 TH1F* hdmH = new TH1F("hdmH", "Hard Trigger Multiplicity Distribution", 100, 0., 8000.);
726
727 mbtH->SetXTitle("b [fm]");
728 hdtH->SetXTitle("b [fm]");
729 mbmH->SetXTitle("Multiplicity");
730 hdmH->SetXTitle("Multiplicity");
731
732 TCanvas *c0 = new TCanvas("c0","Trigger Simulation",400,10,600,700);
733 c0->Divide(2,1);
734 TCanvas *c1 = new TCanvas("c1","Trigger Simulation",400,10,600,700);
735 c1->Divide(1,2);
736
737 //
738 //
739 Init(1);
740 for (Int_t iev = 0; iev < n; iev++)
741 {
742 Float_t b, p, mult;
743 GetRandom(b, p, mult);
744 mbtH->Fill(b,1.);
745 hdtH->Fill(b, p);
746 mbmH->Fill(mult, 1.);
747 hdmH->Fill(mult, p);
748
749 c0->cd(1);
750 mbtH->Draw();
751 c0->cd(2);
752 hdtH->Draw();
753 c0->Update();
754
755 c1->cd(1);
756 mbmH->Draw();
757 c1->cd(2);
758 hdmH->Draw();
759 c1->Update();
760 }
761}
762
763void AliFastGlauber::GetRandom(Float_t& b, Float_t& p, Float_t& mult)
764{
765 //
766 // Gives back a random impact parameter, hard trigger probability and multiplicity
767 //
041f7f97 768 b = fgWSgeo->GetRandom();
769 Float_t mu = fgWSN->Eval(b);
5b3a5a5d 770 p = 1.-TMath::Exp(-mu);
041f7f97 771 mult = 6000./fgWSN->Eval(1.) * mu;
5b3a5a5d 772}
773
c2140715 774void AliFastGlauber::GetRandom(Int_t& bin, Bool_t& hard)
775{
776 //
777 // Gives back a random impact parameter bin, and hard trigger decission
778 //
041f7f97 779 Float_t b = fgWSgeo->GetRandom();
780 Float_t mu = fgWSN->Eval(b) * fSigmaHard;
c2140715 781 Float_t p = 1.-TMath::Exp(-mu);
782 if (b < 5.) {
783 bin = 1;
784 } else if (b < 8.6) {
785 bin = 2;
786 } else if (b < 11.2) {
787 bin = 3;
788 } else if (b < 13.2) {
789 bin = 4;
204e011f 790 } else if (b < 15.0) {
c2140715 791 bin = 5;
204e011f 792 } else {
793 bin = 6;
c2140715 794 }
795
796 hard = kFALSE;
797
798 Float_t r = gRandom->Rndm();
799
800 if (r < p) hard = kTRUE;
801}
802
803
1bc228f5 804
5b3a5a5d 805Float_t AliFastGlauber::GetRandomImpactParameter(Float_t bmin, Float_t bmax)
806{
807 //
808 // Gives back a random impact parameter in the range bmin .. bmax
809 //
810
811 Float_t b = -1.;
812 while(b < bmin || b > bmax)
041f7f97 813 b = fgWSgeo->GetRandom();
5b3a5a5d 814 return b;
815}
816
817Double_t AliFastGlauber::CrossSection(Double_t b1, Double_t b2)
818{
819 //
820 // Return cross-section integrated from b1 to b2
821 //
822
041f7f97 823 return fgWSgeo->Integral(b1, b2)/100.;
5b3a5a5d 824}
825
1bc228f5 826Float_t AliFastGlauber::GetNumberOfParticipants(Float_t b)
827{
828 //
829 // Gives back the number of participants for impact parameter b
830 //
831
832 return (fgWParticipants->Eval(b));
833}
834
5b3a5a5d 835Double_t AliFastGlauber::FractionOfHardCrossSection(Double_t b1, Double_t b2)
836{
837 //
838 // Return raction of hard cross-section integrated from b1 to b2
839 //
840
041f7f97 841 return fgWSbinary->Integral(b1, b2)/fgWSbinary->Integral(0., 100.);
5b3a5a5d 842}
843
844
845Double_t AliFastGlauber::Binaries(Double_t b)
846{
847 //
848 // Return number of binary collisions normalized to 1 at b=0
849 //
850
041f7f97 851 return fgWSN->Eval(b)/fgWSN->Eval(0.001);
5b3a5a5d 852}
a2f2f511 853
854//=================== Added by A. Dainese 11/02/04 ===========================
855void AliFastGlauber::SetCentralityClass(Double_t xsecFrLow,Double_t xsecFrUp)
856{
857 //
858 // Set limits of centrality class as fractions of the geomtrical
859 // cross section
860 //
861 if(xsecFrLow>1. || xsecFrUp>1. || xsecFrLow>xsecFrUp) {
862 printf(" Please set 0 <= xsecFrLow <= xsecFrUp <= 1\n");
863 return;
864 }
865
866 Double_t bLow=0.,bUp=0.;
867 Double_t xsecFr=0.;
868
869 while(xsecFr<xsecFrLow) {
870 xsecFr = fgWSgeo->Integral(0.,bLow)/fgWSgeo->Integral(0.,100.);
871 bLow += 0.1;
872 }
873 bUp = bLow;
874 while(xsecFr<xsecFrUp) {
875 xsecFr = fgWSgeo->Integral(0.,bUp)/fgWSgeo->Integral(0.,100.);
876 bUp += 0.1;
877 }
878
879 printf(" Centrality class: %4.2f-%4.2f; %4.1f < b < %4.1f fm\n",
880 xsecFrLow,xsecFrUp,bLow,bUp);
881
882 fgWSbinary->SetRange(bLow,bUp);
883
884 return;
885}
886
887void AliFastGlauber::GetRandomBHard(Double_t& b)
888{
889 //
890 // Get random impact parameter according to distribution of
891 // hard (binary) cross-section, in the range defined by the centrality class
892 //
893 b = fgWSbinary->GetRandom();
894 Int_t bin = 2*(Int_t)b;
895 if( (b-(Int_t)b) > 0.5) bin++;
896 fWAlmondCurrent = &fWAlmondFixedB[bin];
897 return;
898}
899
900void AliFastGlauber::GetRandomXY(Double_t& x,Double_t& y)
901{
902 //
903 // Get random position of parton production point according to
904 // product of thickness functions
905 //
906 fWAlmondCurrent->GetRandom2(x,y);
907 return;
908}
909
910void AliFastGlauber::GetRandomPhi(Double_t& phi)
911{
912 //
913 // Get random parton azimuthal propagation direction
914 //
915 phi = 2.*TMath::Pi()*gRandom->Rndm();
916 return;
917}
918
919Double_t AliFastGlauber::CalculateLength(Double_t b,
920 Double_t x0,Double_t y0,Double_t phi0)
921{
922 //
923 // Calculate path length for a parton with production point (x0,y0)
924 // and propagation direction (ux=cos(phi0),uy=sin(phi0))
925 // in a collision with impact parameter b
926 //
927
928 // number of steps in l
929 const Int_t kNp = 100;
930 const Double_t kDl = fgBMax/Double_t(kNp);
931
932 Double_t l,integral,integral1,integral2,r0,xx,yy,phi,r1,r2,ell;
933 Double_t rhocoll,rhocollHalfMax;
934 Int_t i,nps;
935
936 if(fEllDef==1) {
937 //
938 // Definition 1:
939 //
940 // ell = 2 * \int_0^\infty dl*l*\rho_coll(x0+l*ux,y0+l*uy) /
941 // \int_0^\infty dl*\rho_coll(x0+l*ux,y0+l*uy)
942 //
943
944
945 l = 0.;
946 integral1 = 0.;
947 integral2 = 0.;
948
949 // Initial radius
950 r0 = TMath::Sqrt(x0 * x0 + y0 * y0);
951 nps = Int_t ((fgBMax - r0)/kDl) - 1;
952
953 // Radial steps
954 for (i = 0; i < nps; i++) {
955
956 // Transform into target frame
957 xx = x0 + l * TMath::Cos(phi0) + b / 2.;
958 yy = y0 + l * TMath::Sin(phi0);
959 phi = TMath::ATan2(yy, xx);
960 r1 = TMath::Sqrt(xx * xx + yy * yy);
961 // Radius in projectile frame
962 r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
963 rhocoll = fgWSta->Eval(r1) * fgWSta->Eval(r2);
964
965 integral1 += rhocoll * l * kDl;
966 integral2 += rhocoll * kDl;
967 l += kDl;
968 } // steps
969
970 ell = (2. * integral1 / integral2);
971 return ell;
972
973 } else if(fEllDef==2) {
974 //
975 // Definition 2:
976 //
977 // ell = \int_0^\infty dl*
978 // \Theta(\rho_coll(x0+l*ux,y0+l*uy)-0.5*\rho_coll(0,0))
979 //
980
981 l = 0.;
982 integral = 0.;
983
984 // Initial radius
985 r0 = TMath::Sqrt(x0 * x0 + y0 * y0);
986 nps = Int_t ((fgBMax - r0)/kDl) - 1;
987
988 rhocollHalfMax = 0.5*fWAlmondCurrent->Eval(0.,0.);
989
990 // Radial steps
991 for (i = 0; i < nps; i++) {
992
993 // Transform into target frame
994 xx = x0 + l * TMath::Cos(phi0) + b / 2.;
995 yy = y0 + l * TMath::Sin(phi0);
996 phi = TMath::ATan2(yy, xx);
997 r1 = TMath::Sqrt(xx * xx + yy * yy);
998 // Radius in projectile frame
999 r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
1000 rhocoll = fgWSta->Eval(r1) * fgWSta->Eval(r2);
1001
1002 if(rhocoll>rhocollHalfMax) integral += kDl;
1003
1004 l += kDl;
1005 } // steps
1006
1007 ell = integral;
1008 return ell;
1009
1010 } else {
1011 printf(" Wrong length definition setting!\n");
1012 return -1.;
1013 }
1014}
1015
1016void AliFastGlauber::GetLength(Double_t& ell,Double_t b)
1017{
1018 //
1019 // Return length from random b, x0, y0, phi0
1020 //
1021 Double_t x0,y0,phi0;
1022 if(b<0.) GetRandomBHard(b);
1023 GetRandomXY(x0,y0);
1024 GetRandomPhi(phi0);
1025
1026 ell = CalculateLength(b,x0,y0,phi0);
1027
1028 return;
1029}
1030
1031void AliFastGlauber::GetLengthsBackToBack(Double_t& ell1,Double_t& ell2,
1032 Double_t b)
1033{
1034 //
1035 // Return 2 lengths back to back from random b, x0, y0, phi0
1036 //
1037 Double_t x0,y0,phi0;
1038 if(b<0.) GetRandomBHard(b);
1039 GetRandomXY(x0,y0);
1040 GetRandomPhi(phi0);
1041 Double_t phi0plusPi = phi0+TMath::Pi();
1042
1043 ell1 = CalculateLength(b,x0,y0,phi0);
1044 ell2 = CalculateLength(b,x0,y0,phi0plusPi);
1045
1046 return;
1047}
1048
1bc228f5 1049void AliFastGlauber::GetLengthsForPythia(Int_t n,Double_t* phi,Double_t* ell, Double_t b)
a2f2f511 1050{
1051 //
1052 // Returns lenghts for n partons with azimuthal angles phi[n]
1053 // from random b, x0, y0
1054 //
1bc228f5 1055 Double_t x0, y0;
1056 if(b < 0.) GetRandomBHard(b);
a2f2f511 1057 GetRandomXY(x0,y0);
1058
1bc228f5 1059 for(Int_t i = 0; i< n; i++) ell[i] = CalculateLength(b,x0,y0,phi[i]);
a2f2f511 1060
1061 return;
1062}
1063
1064void AliFastGlauber::PlotBDistr(Int_t n)
1065{
1066 //
1067 // Plot distribution of n impact parameters
1068 //
1069 Double_t b;
1070 TH1F *hB = new TH1F("hB","dN/db",100,0,fgBMax);
1071 hB->SetXTitle("b [fm]");
1072 hB->SetYTitle("dN/db [a.u.]");
1073 hB->SetFillColor(3);
1074
1075 for(Int_t i=0; i<n; i++) {
1076 GetRandomBHard(b);
1077 hB->Fill(b);
1078 }
1079
1080 TCanvas *cB = new TCanvas("cB","Impact parameter distribution",0,0,500,500);
1081 cB->cd();
1082 hB->Draw();
1083
1084 return;
1085}
1086
1087void AliFastGlauber::PlotLengthDistr(Int_t n,Bool_t save,Char_t *fname)
1088{
1089 //
1090 // Plot length distribution
1091 //
1092 Double_t ell;
1093 TH1F *hEll = new TH1F("hEll","Length distribution",16,-0.5,15.5);
1094 hEll->SetXTitle("Transverse path length, L [fm]");
1095 hEll->SetYTitle("Probability");
1096 hEll->SetFillColor(2);
1097
1098 for(Int_t i=0; i<n; i++) {
1099 GetLength(ell);
1100 hEll->Fill(ell);
1101 }
1102 hEll->Scale(1/(Double_t)n);
1103
1104 TCanvas *cL = new TCanvas("cL","Length distribution",0,0,500,500);
1105 cL->cd();
1106 hEll->Draw();
1107
1108 if(save) {
1109 TFile *f = new TFile(fname,"recreate");
1110 hEll->Write();
1111 f->Close();
1112 }
1113 return;
1114}
1115
1116void AliFastGlauber::PlotLengthB2BDistr(Int_t n,Bool_t save,Char_t *fname)
1117{
1118 //
1119 // Plot lengths back-to-back distributions
1120 //
1121 Double_t ell1,ell2;
1122 TH2F *hElls = new TH2F("hElls","Lengths back-to-back",100,0,15,100,0,15);
1123 hElls->SetXTitle("Transverse path length, L [fm]");
1124 hElls->SetYTitle("Transverse path length, L [fm]");
1125
1126 for(Int_t i=0; i<n; i++) {
1127 GetLengthsBackToBack(ell1,ell2);
1128 hElls->Fill(ell1,ell2);
1129 }
1130 hElls->Scale(1/(Double_t)n);
1131
1132
1133 TCanvas *cLs = new TCanvas("cLs","Length back-to-back distribution",0,0,500,500);
1134 gStyle->SetPalette(1,0);
1135 cLs->cd();
1136 hElls->Draw("col,Z");
1137
1138 if(save) {
1139 TFile *f = new TFile(fname,"recreate");
1140 hElls->Write();
1141 f->Close();
1142 }
1143 return;
1144}
1145
1146void AliFastGlauber::StoreAlmonds()
1147{
1148 //
1149 // Store in glauberPbPb.root 40 almonds for b = (0.25+k*0.5) fm (k=0->39)
1150 //
1151
1152 Char_t almondName[100];
1153 TFile* ff = new TFile("$(ALICE_ROOT)/FASTSIM/data/glauberPbPb.root","update");
1154 for(Int_t k=0; k<40; k++) {
1155 sprintf(almondName,"WAlmondFixedB%d",k);
1156 Double_t b = 0.25+k*0.5;
1157 printf(" b = %f\n",b);
1158 fgWAlmond->SetParameter(0,b);
1159
1160 fgWAlmond->Write(almondName);
1161 }
1162 ff->Close();
1163
1164 return;
1165}
1166
1167void AliFastGlauber::PlotAlmonds()
1168{
1169 //
1170 // Plot almonds for some impact parameters
1171 //
1172
1173 TCanvas *c = new TCanvas("c","Almonds",0,0,500,500);
1174 gStyle->SetPalette(1,0);
1175 c->Divide(2,2);
1176 c->cd(1);
1177 fWAlmondFixedB[0].Draw("cont1");
1178 c->cd(2);
1179 fWAlmondFixedB[10].Draw("cont1");
1180 c->cd(3);
1181 fWAlmondFixedB[20].Draw("cont1");
1182 c->cd(4);
1183 fWAlmondFixedB[30].Draw("cont1");
1184
1185 return;
1186}