git://git.uio.no / u/mrichter/AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Update master to aliroot
[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//
65aa45f2 18// Utility class to make simple Glauber type calculations
a42548b0 19// for SYMMETRIC collision geometries (AA):
041f7f97 20// Impact parameter, production points, reaction plane dependence
65aa45f2 21//
041f7f97 22// The SimulateTrigger method can be used for simple MB and hard-process
23// (binary scaling) trigger studies.
65aa45f2 24//
041f7f97 25// Some basic quantities can be visualized directly.
041f7f97 26//
65aa45f2 27// The default set-up for PbPb or AUAu collisions can be read from a file
28// calling Init(1) or Init(2) if you want to read Almonds too.
29//
30// ***** If you change settings dont forget to call init afterwards, *****
31// ***** in order to update the formulas with the new parameters. *****
041f7f97 32//
33// Author: andreas.morsch@cern.ch
65aa45f2 34//=================== Added by A. Dainese 11/02/04 ===========================
35// Calculate path length for a parton with production point (x0,y0)
36// and propagation direction (ux=cos(phi0),uy=sin(phi0))
37// in a collision with impact parameter b and functions that make use
38// of it.
39//=================== Added by A. Dainese 05/03/04 ===========================
40// Calculation of line integrals I0 and I1
41// integral0 = \int_0^ellCut dl*(T_A*T_B)(x0+l*ux,y0+l*uy)
42// integral1 = \int_0^ellCut dl*l*(T_A*T_B)(x0+l*ux,y0+l*uy)
43// mostly for use in the Quenching class
44//=================== Added by C. Loizdes 27/03/04 ===========================
45// Handling of AuAu collisions
46// More get/set functions
47// Comments, units and clearing of code
48//
5b3a5a5d 49
50// from AliRoot
51#include "AliFastGlauber.h"
52// from root
9edefa04 53#include <TCanvas.h>
5b3a5a5d 54#include <TF1.h>
55#include <TF2.h>
5b3a5a5d 56#include <TFile.h>
9edefa04 57#include <TH1F.h>
58#include <TH2F.h>
59#include <TLegend.h>
60#include <TMath.h>
9edefa04 61#include <TRandom.h>
62#include <TStyle.h>
5b3a5a5d 63
64ClassImp(AliFastGlauber)
65
65aa45f2 66Float_t AliFastGlauber::fgBMax = 0.;
041f7f97 67TF1* AliFastGlauber::fgWSb = NULL;
882b3108 68TF1* AliFastGlauber::fgRWSb = NULL;
041f7f97 69TF2* AliFastGlauber::fgWSbz = NULL;
70TF1* AliFastGlauber::fgWSz = NULL;
71TF1* AliFastGlauber::fgWSta = NULL;
72TF2* AliFastGlauber::fgWStarfi = NULL;
73TF2* AliFastGlauber::fgWAlmond = NULL;
74TF1* AliFastGlauber::fgWStaa = NULL;
75TF1* AliFastGlauber::fgWSgeo = NULL;
76TF1* AliFastGlauber::fgWSbinary = NULL;
77TF1* AliFastGlauber::fgWSN = NULL;
78TF1* AliFastGlauber::fgWPathLength0 = NULL;
79TF1* AliFastGlauber::fgWPathLength = NULL;
80TF1* AliFastGlauber::fgWEnergyDensity = NULL;
81TF1* AliFastGlauber::fgWIntRadius = NULL;
1bc228f5 82TF2* AliFastGlauber::fgWKParticipants = NULL;
83TF1* AliFastGlauber::fgWParticipants = NULL;
65aa45f2 84TF2* AliFastGlauber::fgWAlmondCurrent = NULL;
7f2f270b 85TF2* AliFastGlauber::fgWAlmondFixedB[40];
65aa45f2 86const Int_t AliFastGlauber::fgkMCInts = 100000;
18b7a4a1 87AliFastGlauber* AliFastGlauber::fgGlauber = NULL;
88
5b3a5a5d 89
e6e76983 90AliFastGlauber::AliFastGlauber():
91 fWSr0(0.),
92 fWSd(0.),
93 fWSw(0.),
94 fWSn(0.),
95 fSigmaHard(0.),
96 fSigmaNN(0.),
97 fA(0),
98 fBmin(0.),
99 fBmax(0.),
100 fEllDef(0),
101 fName()
5b3a5a5d 102{
a42548b0 103 // Default Constructor
65aa45f2 104 // Defaults for Pb
105 SetMaxImpact();
106 SetLengthDefinition();
107 SetPbPbLHC();
5373d3f7 108 fXY[0] = fXY[1] = 0;
462421a4 109 fI0I1[0] = fI0I1[1] = 0;
65aa45f2 110}
111
a42548b0 112AliFastGlauber::AliFastGlauber(const AliFastGlauber & gl)
e6e76983 113 :TObject(gl),
114 fWSr0(0.),
115 fWSd(0.),
116 fWSw(0.),
117 fWSn(0.),
118 fSigmaHard(0.),
119 fSigmaNN(0.),
120 fA(0),
121 fBmin(0.),
122 fBmax(0.),
123 fEllDef(0),
124 fName()
a42548b0 125{
126// Copy constructor
127 gl.Copy(*this);
5373d3f7 128 fXY[0] = fXY[1] = 0;
462421a4 129 fI0I1[0] = fI0I1[1] = 0;
a42548b0 130}
131
18b7a4a1 132AliFastGlauber* AliFastGlauber::Instance()
133{
134// Set random number generator
135 if (fgGlauber) {
136 return fgGlauber;
137 } else {
138 fgGlauber = new AliFastGlauber();
139 return fgGlauber;
140 }
141}
142
65aa45f2 143AliFastGlauber::~AliFastGlauber()
144{
a42548b0 145// Destructor
7f2f270b 146 for(Int_t k=0; k<40; k++) delete fgWAlmondFixedB[k];
65aa45f2 147}
148
149void AliFastGlauber::SetAuAuRhic()
150{
151 //Set all parameters for RHIC
152 SetWoodSaxonParametersAu();
153 SetHardCrossSection();
710a8d90 154 SetNNCrossSection(42);
65aa45f2 155 SetNucleus(197);
156 SetFileName("$(ALICE_ROOT)/FASTSIM/data/glauberAuAu.root");
157}
158
159void AliFastGlauber::SetPbPbLHC()
160{
161 //Set all parameters for LHC
162 SetWoodSaxonParametersPb();
163 SetHardCrossSection();
164 SetNNCrossSection();
165 SetNucleus();
166 SetFileName();
5b3a5a5d 167}
168
169void AliFastGlauber::Init(Int_t mode)
170{
65aa45f2 171 // Initialisation
172 // mode = 0; all functions are calculated
173 // mode = 1; overlap function is read from file (for Pb-Pb only)
174 // mode = 2; interaction almond functions are read from file
175 // USE THIS FOR PATH LENGTH CALC.!
176 //
a2f2f511 177
65aa45f2 178 //
179 Reset();
180 //
5b3a5a5d 181
65aa45f2 182 //
183 // Wood-Saxon
184 //
185 fgWSb = new TF1("WSb", WSb, 0, fgBMax, 4);
186 fgWSb->SetParameter(0, fWSr0);
187 fgWSb->SetParameter(1, fWSd);
188 fgWSb->SetParameter(2, fWSw);
189 fgWSb->SetParameter(3, fWSn);
190
882b3108 191 fgRWSb = new TF1("RWSb", RWSb, 0, fgBMax, 4);
192 fgRWSb->SetParameter(0, fWSr0);
193 fgRWSb->SetParameter(1, fWSd);
194 fgRWSb->SetParameter(2, fWSw);
195 fgRWSb->SetParameter(3, fWSn);
196
fac5662b 197 fgWSbz = new TF2("WSbz", WSbz, 0, fgBMax, 0, fgBMax, 4);
65aa45f2 198 fgWSbz->SetParameter(0, fWSr0);
199 fgWSbz->SetParameter(1, fWSd);
200 fgWSbz->SetParameter(2, fWSw);
201 fgWSbz->SetParameter(3, fWSn);
202
203 fgWSz = new TF1("WSz", WSz, 0, fgBMax, 5);
204 fgWSz->SetParameter(0, fWSr0);
205 fgWSz->SetParameter(1, fWSd);
206 fgWSz->SetParameter(2, fWSw);
207 fgWSz->SetParameter(3, fWSn);
208
209 //
210 // Thickness
211 //
212 fgWSta = new TF1("WSta", WSta, 0., fgBMax, 0);
5b3a5a5d 213
65aa45f2 214 //
215 // Overlap Kernel
216 //
217 fgWStarfi = new TF2("WStarfi", WStarfi, 0., fgBMax, 0., TMath::Pi(), 1);
218 fgWStarfi->SetParameter(0, 0.);
219 fgWStarfi->SetNpx(200);
220 fgWStarfi->SetNpy(20);
8de7e046 221
65aa45f2 222 //
223 // Participants Kernel
224 //
225 fgWKParticipants = new TF2("WKParticipants", WKParticipants, 0., fgBMax, 0., TMath::Pi(), 3);
226 fgWKParticipants->SetParameter(0, 0.);
227 fgWKParticipants->SetParameter(1, fSigmaNN);
228 fgWKParticipants->SetParameter(2, fA);
229 fgWKParticipants->SetNpx(200);
230 fgWKParticipants->SetNpy(20);
f3a04204 231
65aa45f2 232 //
233 // Overlap and Participants
234 //
235 if (!mode) {
236 fgWStaa = new TF1("WStaa", WStaa, 0., fgBMax, 1);
237 fgWStaa->SetNpx(100);
238 fgWStaa->SetParameter(0,fA);
239 fgWStaa->SetNpx(100);
240 fgWParticipants = new TF1("WParticipants", WParticipants, 0., fgBMax, 2);
241 fgWParticipants->SetParameter(0, fSigmaNN);
242 fgWParticipants->SetParameter(1, fA);
243 fgWParticipants->SetNpx(100);
244 } else {
245 Info("Init","Reading overlap function from file %s",fName.Data());
246 TFile* f = new TFile(fName.Data());
247 if(!f->IsOpen()){
248 Fatal("Init", "Could not open file %s",fName.Data());
249 }
250 fgWStaa = (TF1*) f->Get("WStaa");
251 fgWParticipants = (TF1*) f->Get("WParticipants");
252 delete f;
253 }
2a103154 254
65aa45f2 255 //
256 // Energy Density
257 //
258 fgWEnergyDensity = new TF1("WEnergyDensity", WEnergyDensity, 0., 2. * fWSr0, 1);
259 fgWEnergyDensity->SetParameter(0, fWSr0 + 1.);
260
261 //
262 // Geometrical Cross-Section
263 //
264 fgWSgeo = new TF1("WSgeo", WSgeo, 0., fgBMax, 1);
265 fgWSgeo->SetParameter(0,fSigmaNN); //mbarn
266 fgWSgeo->SetNpx(100);
2a103154 267
65aa45f2 268 //
269 // Hard cross section (binary collisions)
270 //
271 fgWSbinary = new TF1("WSbinary", WSbinary, 0., fgBMax, 1);
272 fgWSbinary->SetParameter(0, fSigmaHard); //mbarn
273 fgWSbinary->SetNpx(100);
1bc228f5 274
65aa45f2 275 //
276 // Hard collisions per event
277 //
882b3108 278 fgWSN = new TF1("WSN", WSN, 0.01, fgBMax, 1);
65aa45f2 279 fgWSN->SetNpx(100);
280
281 //
282 // Almond shaped interaction region
283 //
284 fgWAlmond = new TF2("WAlmond", WAlmond, -fgBMax, fgBMax, -fgBMax, fgBMax, 1);
285 fgWAlmond->SetParameter(0, 0.);
286 fgWAlmond->SetNpx(200);
287 fgWAlmond->SetNpy(200);
288
289 if(mode==2) {
290 Info("Init","Reading interaction almonds from file: %s",fName.Data());
291 Char_t almondName[100];
292 TFile* ff = new TFile(fName.Data());
293 for(Int_t k=0; k<40; k++) {
bae6da3c 294 snprintf(almondName,100, "WAlmondFixedB%d",k);
65aa45f2 295 fgWAlmondCurrent = (TF2*)ff->Get(almondName);
7f2f270b 296 fgWAlmondFixedB[k] = fgWAlmondCurrent;
5b3a5a5d 297 }
65aa45f2 298 delete ff;
299 }
300
301 fgWIntRadius = new TF1("WIntRadius", WIntRadius, 0., fgBMax, 1);
302 fgWIntRadius->SetParameter(0, 0.);
303
304 //
305 // Path Length as a function of Phi
306 //
307 fgWPathLength0 = new TF1("WPathLength0", WPathLength0, -TMath::Pi(), TMath::Pi(), 2);
308 fgWPathLength0->SetParameter(0, 0.);
309 fgWPathLength0->SetParameter(1, 0.); //Pathlength definition
310
311 fgWPathLength = new TF1("WPathLength", WPathLength, -TMath::Pi(), TMath::Pi(), 3);
312 fgWPathLength->SetParameter(0, 0.); //Impact Parameter
313 fgWPathLength->SetParameter(1, 1000.); //Number of interactions used for average
314 fgWPathLength->SetParameter(2, 0); //Pathlength definition
315}
316
a42548b0 317void AliFastGlauber::Reset() const
65aa45f2 318{
319 //
320 // Reset dynamic allocated formulas
321 // in case init is called twice
322
323 if(fgWSb) delete fgWSb;
882b3108 324 if(fgRWSb) delete fgRWSb;
65aa45f2 325 if(fgWSbz) delete fgWSbz;
326 if(fgWSz) delete fgWSz;
327 if(fgWSta) delete fgWSta;
328 if(fgWStarfi) delete fgWStarfi;
329 if(fgWAlmond) delete fgWAlmond;
330 if(fgWStaa) delete fgWStaa;
331 if(fgWSgeo) delete fgWSgeo;
332 if(fgWSbinary) delete fgWSbinary;
333 if(fgWSN) delete fgWSN;
334 if(fgWPathLength0) delete fgWPathLength0;
335 if(fgWPathLength) delete fgWPathLength;
336 if(fgWEnergyDensity) delete fgWEnergyDensity;
337 if(fgWIntRadius) delete fgWIntRadius;
338 if(fgWKParticipants) delete fgWKParticipants;
339 if(fgWParticipants) delete fgWParticipants;
5b3a5a5d 340}
341
710a8d90 342void AliFastGlauber::DrawWSb() const
5b3a5a5d 343{
65aa45f2 344 //
345 // Draw Wood-Saxon Nuclear Density Function
346 //
347 TCanvas *c1 = new TCanvas("c1","Wood Saxon",400,10,600,700);
348 c1->cd();
882b3108 349 Double_t max=fgWSb->GetMaximum(0.,fgBMax)*1.01;
65aa45f2 350 TH2F *h2f=new TH2F("h2fwsb","Wood Saxon: #rho(r) = n (1-#omega(r/r_{0})^2)/(1+exp((r-r_{0})/d)) [fm^{-3}]",2,0,fgBMax,2,0,max);
351 h2f->SetStats(0);
352 h2f->GetXaxis()->SetTitle("r [fm]");
353 h2f->GetYaxis()->SetNoExponent(kTRUE);
354 h2f->GetYaxis()->SetTitle("#rho [fm^{-3}]");
355 h2f->Draw();
356 fgWSb->Draw("same");
357 TLegend *l1a = new TLegend(0.45,0.6,.90,0.8);
358 l1a->SetFillStyle(0);
359 l1a->SetBorderSize(0);
360 Char_t label[100];
bae6da3c 361 snprintf(label,100, "r_{0} = %.2f fm",fWSr0);
65aa45f2 362 l1a->AddEntry(fgWSb,label,"");
bae6da3c 363 snprintf(label,100, "d = %.2f fm",fWSd);
65aa45f2 364 l1a->AddEntry(fgWSb,label,"");
bae6da3c 365 snprintf(label,100, "n = %.2e fm^{-3}",fWSn);
65aa45f2 366 l1a->AddEntry(fgWSb,label,"");
bae6da3c 367 snprintf(label,100, "#omega = %.2f",fWSw);
65aa45f2 368 l1a->AddEntry(fgWSb,label,"");
369 l1a->Draw();
370 c1->Update();
5b3a5a5d 371}
372
710a8d90 373void AliFastGlauber::DrawOverlap() const
5b3a5a5d 374{
65aa45f2 375 //
376 // Draw Overlap Function
377 //
378 TCanvas *c2 = new TCanvas("c2","Overlap",400,10,600,700);
379 c2->cd();
380 Double_t max=fgWStaa->GetMaximum(0,fgBMax)*1.01;
a42548b0 381 TH2F *h2f=new TH2F("h2ftaa","Overlap function: T_{AB} [mbarn^{-1}]",2,0,fgBMax,2,0, max);
65aa45f2 382 h2f->SetStats(0);
383 h2f->GetXaxis()->SetTitle("b [fm]");
384 h2f->GetYaxis()->SetTitle("T_{AB} [mbarn^{-1}]");
385 h2f->Draw();
386 fgWStaa->Draw("same");
5b3a5a5d 387}
388
710a8d90 389void AliFastGlauber::DrawParticipants() const
1bc228f5 390{
65aa45f2 391 //
392 // Draw Number of Participants Npart
393 //
394 TCanvas *c3 = new TCanvas("c3","Participants",400,10,600,700);
395 c3->cd();
396 Double_t max=fgWParticipants->GetMaximum(0,fgBMax)*1.01;
397 TH2F *h2f=new TH2F("h2fpart","Number of Participants",2,0,fgBMax,2,0,max);
398 h2f->SetStats(0);
399 h2f->GetXaxis()->SetTitle("b [fm]");
400 h2f->GetYaxis()->SetTitle("N_{part}");
401 h2f->Draw();
402 fgWParticipants->Draw("same");
403 TLegend *l1a = new TLegend(0.50,0.75,.90,0.9);
404 l1a->SetFillStyle(0);
405 l1a->SetBorderSize(0);
406 Char_t label[100];
bae6da3c 407 snprintf(label,100, "#sigma^{inel.}_{NN} = %.1f mbarn",fSigmaNN);
65aa45f2 408 l1a->AddEntry(fgWParticipants,label,"");
409 l1a->Draw();
410 c3->Update();
1bc228f5 411}
412
710a8d90 413void AliFastGlauber::DrawThickness() const
5b3a5a5d 414{
65aa45f2 415 //
416 // Draw Thickness Function
417 //
418 TCanvas *c4 = new TCanvas("c4","Thickness",400,10,600,700);
419 c4->cd();
420 Double_t max=fgWSta->GetMaximum(0,fgBMax)*1.01;
421 TH2F *h2f=new TH2F("h2fta","Thickness function: T_{A} [fm^{-2}]",2,0,fgBMax,2,0,max);
422 h2f->SetStats(0);
423 h2f->GetXaxis()->SetTitle("b [fm]");
424 h2f->GetYaxis()->SetTitle("T_{A} [fm^{-2}]");
425 h2f->Draw();
426 fgWSta->Draw("same");
5b3a5a5d 427}
428
710a8d90 429void AliFastGlauber::DrawGeo() const
5b3a5a5d 430{
65aa45f2 431 //
432 // Draw Geometrical Cross-Section
433 //
434 TCanvas *c5 = new TCanvas("c5","Geometrical Cross-Section",400,10,600,700);
435 c5->cd();
436 Double_t max=fgWSgeo->GetMaximum(0,fgBMax)*1.01;
437 TH2F *h2f=new TH2F("h2fgeo","Differential Geometrical Cross-Section: d#sigma^{geo}_{AB}/db [fm]",2,0,fgBMax,2,0,max);
438 h2f->SetStats(0);
439 h2f->GetXaxis()->SetTitle("b [fm]");
440 h2f->GetYaxis()->SetTitle("d#sigma^{geo}_{AB}/db [fm]");
441 h2f->Draw();
442 fgWSgeo->Draw("same");
443 TLegend *l1a = new TLegend(0.10,0.8,.40,0.9);
444 l1a->SetFillStyle(0);
445 l1a->SetBorderSize(0);
446 Char_t label[100];
bae6da3c 447 snprintf(label,100, "#sigma_{NN}^{inel.} = %.1f mbarn",fSigmaNN);
65aa45f2 448 l1a->AddEntry(fgWSgeo,label,"");
449 l1a->Draw();
450 c5->Update();
5b3a5a5d 451}
452
710a8d90 453void AliFastGlauber::DrawBinary() const
5b3a5a5d 454{
65aa45f2 455 //
456 // Draw Binary Cross-Section
457 //
458 TCanvas *c6 = new TCanvas("c6","Binary Cross-Section",400,10,600,700);
459 c6->cd();
460 Double_t max=fgWSbinary->GetMaximum(0,fgBMax)*1.01;
461 TH2F *h2f=new TH2F("h2fbinary","Differential Binary Cross-Section: #sigma^{hard}_{NN} dT_{AB}/db [fm]",2,0,fgBMax,2,0,max);
462 h2f->SetStats(0);
463 h2f->GetXaxis()->SetTitle("b [fm]");
464 h2f->GetYaxis()->SetTitle("d#sigma^{hard}_{AB}/db [fm]");
465 h2f->Draw();
466 fgWSbinary->Draw("same");
467 TLegend *l1a = new TLegend(0.50,0.8,.90,0.9);
468 l1a->SetFillStyle(0);
469 l1a->SetBorderSize(0);
470 Char_t label[100];
bae6da3c 471 snprintf(label,100, "#sigma_{NN}^{hard} = %.1f mbarn",fSigmaHard);
65aa45f2 472 l1a->AddEntry(fgWSb,label,"");
473 l1a->Draw();
474 c6->Update();
5b3a5a5d 475}
476
710a8d90 477void AliFastGlauber::DrawN() const
5b3a5a5d 478{
65aa45f2 479 //
480 // Draw Binaries per event (Ncoll)
481 //
482 TCanvas *c7 = new TCanvas("c7","Binaries per event",400,10,600,700);
483 c7->cd();
882b3108 484 Double_t max=fgWSN->GetMaximum(0.01,fgBMax)*1.01;
65aa45f2 485 TH2F *h2f=new TH2F("h2fhardcols","Number of hard collisions: T_{AB} #sigma^{hard}_{NN}/#sigma_{AB}^{geo}",2,0,fgBMax,2,0,max);
486 h2f->SetStats(0);
487 h2f->GetXaxis()->SetTitle("b [fm]");
488 h2f->GetYaxis()->SetTitle("N_{coll}");
489 h2f->Draw();
490 fgWSN->Draw("same");
491 TLegend *l1a = new TLegend(0.50,0.75,.90,0.9);
492 l1a->SetFillStyle(0);
493 l1a->SetBorderSize(0);
494 Char_t label[100];
bae6da3c 495 snprintf(label,100, "#sigma^{hard}_{NN} = %.1f mbarn",fSigmaHard);
65aa45f2 496 l1a->AddEntry(fgWSN,label,"");
bae6da3c 497 snprintf(label,100, "#sigma^{inel.}_{NN} = %.1f mbarn",fSigmaNN);
65aa45f2 498 l1a->AddEntry(fgWSN,label,"");
499 l1a->Draw();
500 c7->Update();
5b3a5a5d 501}
502
710a8d90 503void AliFastGlauber::DrawKernel(Double_t b) const
5b3a5a5d 504{
65aa45f2 505 //
506 // Draw Kernel
507 //
508 TCanvas *c8 = new TCanvas("c8","Kernel",400,10,600,700);
509 c8->cd();
510 fgWStarfi->SetParameter(0, b);
511 TH2F *h2f=new TH2F("h2fkernel","Kernel of Overlap function: d^{2}T_{AB}/dr/d#phi [fm^{-3}]",2,0,fgBMax,2,0,TMath::Pi());
512 h2f->SetStats(0);
513 h2f->GetXaxis()->SetTitle("r [fm]");
514 h2f->GetYaxis()->SetTitle("#phi [rad]");
515 h2f->Draw();
516 fgWStarfi->Draw("same");
517 TLegend *l1a = new TLegend(0.65,0.8,.90,0.9);
518 l1a->SetFillStyle(0);
519 l1a->SetBorderSize(0);
520 Char_t label[100];
bae6da3c 521 snprintf(label, 100, "b = %.1f fm",b);
65aa45f2 522 l1a->AddEntry(fgWStarfi,label,"");
523 l1a->Draw();
524 c8->Update();
5b3a5a5d 525}
526
710a8d90 527void AliFastGlauber::DrawAlmond(Double_t b) const
f3a04204 528{
65aa45f2 529 //
530 // Draw Interaction Almond
531 //
532 TCanvas *c9 = new TCanvas("c9","Almond",400,10,600,700);
533 c9->cd();
534 fgWAlmond->SetParameter(0, b);
882b3108 535 TH2F *h2f=new TH2F("h2falmond","Interaction Almond [fm^{-4}]",2,-fgBMax, fgBMax, 2, -fgBMax, fgBMax);
65aa45f2 536 h2f->SetStats(0);
537 h2f->GetXaxis()->SetTitle("x [fm]");
538 h2f->GetYaxis()->SetTitle("y [fm]");
882b3108 539 h2f->Draw("");
540 gStyle->SetPalette(1);
541 fgWAlmond->Draw("colzsame");
65aa45f2 542 TLegend *l1a = new TLegend(0.65,0.8,.90,0.9);
543 l1a->SetFillStyle(0);
544 l1a->SetBorderSize(0);
545 Char_t label[100];
bae6da3c 546 snprintf(label, 100, "b = %.1f fm",b);
65aa45f2 547 l1a->AddEntry(fgWAlmond,label,"");
548 l1a->Draw();
549 c9->Update();
f3a04204 550}
551
710a8d90 552void AliFastGlauber::DrawEnergyDensity() const
f3a04204 553{
65aa45f2 554 //
555 // Draw energy density
556 //
557 TCanvas *c10 = new TCanvas("c10","Energy Density",400, 10, 600, 700);
558 c10->cd();
559 fgWEnergyDensity->SetMinimum(0.);
560 Double_t max=fgWEnergyDensity->GetMaximum(0,fgWEnergyDensity->GetParameter(0))*1.01;
561 TH2F *h2f=new TH2F("h2fenergydens","Energy density",2,0,fgBMax,2,0,max);
562 h2f->SetStats(0);
563 h2f->GetXaxis()->SetTitle("b [fm]");
564 h2f->GetYaxis()->SetTitle("fm^{-4}");
565 h2f->Draw();
566 fgWEnergyDensity->Draw("same");
567 c10->Update();
f3a04204 568}
569
710a8d90 570void AliFastGlauber::DrawPathLength0(Double_t b, Int_t iopt) const
f3a04204 571{
65aa45f2 572 //
573 // Draw Path Length
574 //
575 TCanvas *c11 = new TCanvas("c11","Path Length",400,10,600,700);
576 c11->cd();
577 fgWPathLength0->SetParameter(0, b);
578 fgWPathLength0->SetParameter(1, Double_t(iopt));
579 fgWPathLength0->SetMinimum(0.);
580 fgWPathLength0->SetMaximum(10.);
581 TH2F *h2f=new TH2F("h2fpathlength0","Path length",2,-TMath::Pi(), TMath::Pi(),2,0,10.);
582 h2f->SetStats(0);
583 h2f->GetXaxis()->SetTitle("#phi [rad]");
584 h2f->GetYaxis()->SetTitle("l [fm]");
585 h2f->Draw();
586 fgWPathLength0->Draw("same");
f3a04204 587}
588
710a8d90 589void AliFastGlauber::DrawPathLength(Double_t b , Int_t ni, Int_t iopt) const
f3a04204 590{
65aa45f2 591 //
592 // Draw Path Length
593 //
594 TCanvas *c12 = new TCanvas("c12","Path Length",400,10,600,700);
595 c12->cd();
596 fgWAlmond->SetParameter(0, b);
597 fgWPathLength->SetParameter(0, b);
598 fgWPathLength->SetParameter(1, Double_t (ni));
599 fgWPathLength->SetParameter(2, Double_t (iopt));
600 fgWPathLength->SetMinimum(0.);
601 fgWPathLength->SetMaximum(10.);
602 TH2F *h2f=new TH2F("h2fpathlength","Path length",2,-TMath::Pi(), TMath::Pi(),2,0,10.);
603 h2f->SetStats(0);
604 h2f->GetXaxis()->SetTitle("#phi [rad]");
605 h2f->GetYaxis()->SetTitle("l [fm]");
606 h2f->Draw();
607 fgWPathLength->Draw("same");
f3a04204 608}
609
710a8d90 610void AliFastGlauber::DrawIntRadius(Double_t b) const
2a103154 611{
65aa45f2 612 //
613 // Draw Interaction Radius
614 //
615 TCanvas *c13 = new TCanvas("c13","Interaction Radius",400,10,600,700);
616 c13->cd();
617 fgWIntRadius->SetParameter(0, b);
618 fgWIntRadius->SetMinimum(0);
619 Double_t max=fgWIntRadius->GetMaximum(0,fgBMax)*1.01;
620 TH2F *h2f=new TH2F("h2fintradius","Interaction Density",2,0.,fgBMax,2,0,max);
621 h2f->SetStats(0);
622 h2f->GetXaxis()->SetTitle("r [fm]");
623 h2f->GetYaxis()->SetTitle("[fm^{-3}]");
624 h2f->Draw();
625 fgWIntRadius->Draw("same");
2a103154 626}
627
0c13ae32 628Double_t AliFastGlauber::WSb(const Double_t* x, const Double_t* par)
5b3a5a5d 629{
65aa45f2 630 //
631 // Woods-Saxon Parameterisation
632 // as a function of radius (xx)
633 //
710a8d90 634 const Double_t kxx = x[0]; //fm
635 const Double_t kr0 = par[0]; //fm
636 const Double_t kd = par[1]; //fm
637 const Double_t kw = par[2]; //no units
638 const Double_t kn = par[3]; //fm^-3 (used to normalize integral to one)
639 Double_t y = kn * (1.+kw*(kxx/kr0)*(kxx/kr0))/(1.+TMath::Exp((kxx-kr0)/kd));
65aa45f2 640 return y; //fm^-3
5b3a5a5d 641}
642
882b3108 643Double_t AliFastGlauber::RWSb(const Double_t* x, const Double_t* par)
644{
645 //
646 // Woods-Saxon Parameterisation
647 // as a function of radius (xx)
648 // times r**2
649 const Double_t kxx = x[0]; //fm
650 const Double_t kr0 = par[0]; //fm
651 const Double_t kd = par[1]; //fm
652 const Double_t kw = par[2]; //no units
653 const Double_t kn = par[3]; //fm^-3 (used to normalize integral to one)
654 Double_t y = kxx * kxx * kn * (1.+kw*(kxx/kr0)*(kxx/kr0))/(1.+TMath::Exp((kxx-kr0)/kd));
655
656 return y; //fm^-1
657}
658
0c13ae32 659Double_t AliFastGlauber::WSbz(const Double_t* x, const Double_t* par)
5b3a5a5d 660{
65aa45f2 661 //
662 // Wood Saxon Parameterisation
663 // as a function of z and b
664 //
710a8d90 665 const Double_t kbb = x[0]; //fm
666 const Double_t kzz = x[1]; //fm
667 const Double_t kr0 = par[0]; //fm
668 const Double_t kd = par[1]; //fm
669 const Double_t kw = par[2]; //no units
670 const Double_t kn = par[3]; //fm^-3 (used to normalize integral to one)
671 const Double_t kxx = TMath::Sqrt(kbb*kbb+kzz*kzz);
672 Double_t y = kn * (1.+kw*(kxx/kr0)*(kxx/kr0))/(1.+TMath::Exp((kxx-kr0)/kd));
65aa45f2 673 return y; //fm^-3
5b3a5a5d 674}
675
0c13ae32 676Double_t AliFastGlauber::WSz(const Double_t* x, const Double_t* par)
5b3a5a5d 677{
65aa45f2 678 //
679 // Wood Saxon Parameterisation
680 // as a function of z for fixed b
681 //
710a8d90 682 const Double_t kzz = x[0]; //fm
683 const Double_t kr0 = par[0]; //fm
684 const Double_t kd = par[1]; //fm
685 const Double_t kw = par[2]; //no units
686 const Double_t kn = par[3]; //fm^-3 (used to normalize integral to one)
687 const Double_t kbb = par[4]; //fm
688 const Double_t kxx = TMath::Sqrt(kbb*kbb+kzz*kzz);
689 Double_t y = kn * (1.+kw*(kxx/kr0)*(kxx/kr0))/(1.+TMath::Exp((kxx-kr0)/kd));
65aa45f2 690 return y; //fm^-3
5b3a5a5d 691}
692
0c13ae32 693Double_t AliFastGlauber::WSta(const Double_t* x, const Double_t* /*par*/)
5b3a5a5d 694{
65aa45f2 695 //
696 // Thickness function T_A
697 // as a function of b
698 //
710a8d90 699 const Double_t kb = x[0];
700 fgWSz->SetParameter(4, kb);
701 Double_t y = 2. * fgWSz->Integral(0., fgBMax);
65aa45f2 702 return y; //fm^-2
5b3a5a5d 703}
704
0c13ae32 705Double_t AliFastGlauber::WStarfi(const Double_t* x, const Double_t* par)
5b3a5a5d 706{
65aa45f2 707 //
708 // Kernel for overlap function: T_A(s)*T_A(s-b)
709 // as a function of r and phi
710a8d90 710 const Double_t kr1 = x[0];
711 const Double_t kphi = x[1];
712 const Double_t kb = par[0];
713 const Double_t kr2 = TMath::Sqrt(kr1*kr1 + kb*kb - 2.*kr1*kb*TMath::Cos(kphi));
714 Double_t y = kr1 * fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
65aa45f2 715 return y; //fm^-3
716}
717
0c13ae32 718Double_t AliFastGlauber::WStaa(const Double_t* x, const Double_t* par)
65aa45f2 719{
720 //
721 // Overlap function
722 // T_{AB}=Int d2s T_A(s)*T_B(s-b)
723 // as a function of b
724 // (normalized to fA*fB)
725 //
710a8d90 726 const Double_t kb = x[0];
727 const Double_t ka = par[0];
728 fgWStarfi->SetParameter(0, kb);
65aa45f2 729
730 // root integration seems to fail
731 /*
732 Double_t al[2];
733 Double_t bl[2];
734 al[0] = 1e-6;
735 al[1] = fgBMax;
736 bl[0] = 0.;
737 bl[1] = TMath::Pi();
738 Double_t err;
739
740 Double_t y = 2. * 208. * 208. * fgWStarfi->IntegralMultiple(2, al, bl, 0.001, err);
741 printf("WStaa: %.5e %.5e %.5e\n", b, y, err);
742 */
743
744 //
745 // MC Integration
746 //
747 Double_t y = 0;
a42548b0 748
749
65aa45f2 750 for (Int_t i = 0; i < fgkMCInts; i++)
751 {
a42548b0 752
710a8d90 753 const Double_t kphi = TMath::Pi() * gRandom->Rndm();
754 const Double_t kb1 = fgBMax * gRandom->Rndm();
755 y += fgWStarfi->Eval(kb1, kphi);
65aa45f2 756 }
757 y *= 2. * TMath::Pi() * fgBMax / fgkMCInts; //fm^-2
710a8d90 758 y *= ka * ka * 0.1; //mbarn^-1
65aa45f2 759 return y;
5b3a5a5d 760}
761
0c13ae32 762Double_t AliFastGlauber::WKParticipants(const Double_t* x, const Double_t* par)
1bc228f5 763{
65aa45f2 764 //
765 // Kernel for number of participants
766 // as a function of r and phi
767 //
710a8d90 768 const Double_t kr1 = x[0];
769 const Double_t kphi = x[1];
770 const Double_t kb = par[0]; //fm
771 const Double_t ksig = par[1]; //mbarn
772 const Double_t ka = par[2]; //mass number
773 const Double_t kr2 = TMath::Sqrt(kr1*kr1 +kb*kb - 2.*kr1*kb*TMath::Cos(kphi));
774 const Double_t kxsi = fgWSta->Eval(kr2) * ksig * 0.1; //no units
65aa45f2 775 /*
710a8d90 776 Double_t y=(1-TMath::Power((1-xsi),aa))
65aa45f2 777 */
710a8d90 778 Double_t a = ka;
779 Double_t sum = ka * kxsi;
65aa45f2 780 Double_t y = sum;
710a8d90 781 for (Int_t i = 1; i <= ka; i++)
1bc228f5 782 {
65aa45f2 783 a--;
710a8d90 784 sum *= (-kxsi) * a / Float_t(i+1);
65aa45f2 785 y += sum;
1bc228f5 786 }
710a8d90 787 y *= kr1 * fgWSta->Eval(kr1);
65aa45f2 788 return y; //fm^-1
1bc228f5 789}
790
0c13ae32 791Double_t AliFastGlauber::WParticipants(const Double_t* x, const Double_t* par)
65aa45f2 792{
793 //
794 // Number of Participants as
795 // a function of b
796 //
710a8d90 797 const Double_t kb = x[0];
798 const Double_t ksig = par[0]; //mbarn
799 const Double_t ka = par[1]; //mass number
800 fgWKParticipants->SetParameter(0, kb);
801 fgWKParticipants->SetParameter(1, ksig);
802 fgWKParticipants->SetParameter(2, ka);
65aa45f2 803
804 //
805 // MC Integration
806 //
807 Double_t y = 0;
808 for (Int_t i = 0; i < fgkMCInts; i++)
809 {
710a8d90 810 const Double_t kphi = TMath::Pi() * gRandom->Rndm();
811 const Double_t kb1 = fgBMax * gRandom->Rndm();
812 y += fgWKParticipants->Eval(kb1, kphi);
65aa45f2 813 }
710a8d90 814 y *= 2. * ka * 2. * TMath::Pi() * fgBMax / fgkMCInts;
65aa45f2 815 return y; //no units
816}
817
0c13ae32 818Double_t AliFastGlauber::WSgeo(const Double_t* x, const Double_t* par)
65aa45f2 819{
820 //
821 // Geometrical Cross-Section
822 // as a function of b
823 //
710a8d90 824 const Double_t kb = x[0]; //fm
825 const Double_t ksigNN = par[0]; //mbarn
826 const Double_t ktaa = fgWStaa->Eval(kb); //mbarn^-1
827 Double_t y = 2. * TMath::Pi() * kb * (1. - TMath::Exp(- ksigNN * ktaa));
65aa45f2 828 return y; //fm
829}
830
0c13ae32 831Double_t AliFastGlauber::WSbinary(const Double_t* x, const Double_t* par)
65aa45f2 832{
833 //
710a8d90 834 // Number of binary hard collisions
65aa45f2 835 // as a function of b
836 //
710a8d90 837 const Double_t kb = x[0]; //fm
838 const Double_t ksig = par[0]; //mbarn
839 const Double_t ktaa = fgWStaa->Eval(kb); //mbarn^-1
840 Double_t y = 2. * TMath::Pi() * kb * ksig * ktaa;
65aa45f2 841 return y; //fm
842}
843
0c13ae32 844Double_t AliFastGlauber::WSN(const Double_t* x, const Double_t* /*par*/)
65aa45f2 845{
846 //
847 // Number of hard processes per event
848 // as a function of b
710a8d90 849 const Double_t kb = x[0];
850 Double_t y = fgWSbinary->Eval(kb)/fgWSgeo->Eval(kb);
65aa45f2 851 return y; //no units
852}
853
0c13ae32 854Double_t AliFastGlauber::WEnergyDensity(const Double_t* x, const Double_t* par)
65aa45f2 855{
856 //
857 // Initial energy density
858 // as a function of the impact parameter
859 //
710a8d90 860 const Double_t kb = x[0];
861 const Double_t krA = par[0];
65aa45f2 862 //
863 // Attention: area of transverse reaction zone in hard-sphere approximation !
710a8d90 864 const Double_t krA2=krA*krA;
865 const Double_t kb2=kb*kb;
866 const Double_t ksaa = (TMath::Pi() - 2. * TMath::ASin(kb/ 2./ krA)) * krA2
867 - kb * TMath::Sqrt(krA2 - kb2/ 4.); //fm^2
868 const Double_t ktaa = fgWStaa->Eval(kb); //mbarn^-1
869 Double_t y=ktaa/ksaa*10;
65aa45f2 870 return y; //fm^-4
871}
5b3a5a5d 872
0c13ae32 873Double_t AliFastGlauber::WAlmond(const Double_t* x, const Double_t* par)
f3a04204 874{
65aa45f2 875 //
876 // Almond shaped interaction region
877 // as a function of cartesian x,y.
878 //
710a8d90 879 const Double_t kb = par[0];
880 const Double_t kxx = x[0] + kb/2.;
881 const Double_t kyy = x[1];
882 const Double_t kr1 = TMath::Sqrt(kxx*kxx + kyy*kyy);
883 const Double_t kphi = TMath::ATan2(kyy,kxx);
884 const Double_t kr2 = TMath::Sqrt(kr1*kr1 + kb*kb - 2.*kr1*kb*TMath::Cos(kphi));
65aa45f2 885 //
886 // Interaction probability calculated as product of thicknesses
887 //
710a8d90 888 Double_t y = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
65aa45f2 889 return y; //fm^-4
f3a04204 890}
891
0c13ae32 892Double_t AliFastGlauber::WIntRadius(const Double_t* x, const Double_t* par)
f3a04204 893{
65aa45f2 894 //
895 // Average interaction density over radius
896 // at which interaction takes place
897 // as a function of radius
898 //
710a8d90 899 const Double_t kr = x[0];
900 const Double_t kb = par[0];
901 fgWAlmond->SetParameter(0, kb);
65aa45f2 902 // Average over phi in small steps
710a8d90 903 const Double_t kdphi = 2. * TMath::Pi() / 100.;
65aa45f2 904 Double_t phi = 0.;
905 Double_t y = 0.;
906 for (Int_t i = 0; i < 100; i++) {
710a8d90 907 const Double_t kxx = kr * TMath::Cos(phi);
908 const Double_t kyy = kr * TMath::Sin(phi);
909 y += fgWAlmond->Eval(kxx,kyy);
910 phi += kdphi;
65aa45f2 911 } // phi loop
912 // Result multiplied by Jacobian (2 pi r)
710a8d90 913 y *= 2. * TMath::Pi() * kr / 100.;
65aa45f2 914 return y; //fm^-3
f3a04204 915}
916
0c13ae32 917Double_t AliFastGlauber::WPathLength0(const Double_t* x, const Double_t* par)
f3a04204 918{
65aa45f2 919 //
920 // Path Length as a function of phi
921 // for interaction point fixed at (0,0)
922 // as a function of phi-direction
923 //
924 // Phi direction in Almond
710a8d90 925 const Double_t kphi0 = x[0];
926 const Double_t kb = par[0];
65aa45f2 927 // Path Length definition
710a8d90 928 const Int_t kiopt = Int_t(par[1]);
65aa45f2 929
930 // Step along radial direction phi
931 const Int_t kNp = 100; // Steps in r
932 const Double_t kDr = fgBMax/kNp;
933 Double_t r = 0.;
934 Double_t rw = 0.;
935 Double_t w = 0.;
936 for (Int_t i = 0; i < kNp; i++) {
937 //
938 // Transform into target frame
939 //
710a8d90 940 const Double_t kxx = r * TMath::Cos(kphi0) + kb / 2.;
941 const Double_t kyy = r * TMath::Sin(kphi0);
942 const Double_t kphi = TMath::ATan2(kyy, kxx);
943 const Double_t kr1 = TMath::Sqrt(kxx*kxx + kyy*kyy);
65aa45f2 944 // Radius in projectile frame
710a8d90 945 const Double_t kr2 = TMath::Sqrt(kr1*kr1 + kb*kb - 2.*kr1*kb*TMath::Cos(kphi));
946 const Double_t ky = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
65aa45f2 947
710a8d90 948 rw += ky * r;
949 w += ky;
65aa45f2 950 r += kDr;
951 } // radial steps
952
953 Double_t y=0.;
ed79c723 954 if (!kiopt) { // My length definition (is exact for hard disk)
955 if(w) y= 2. * rw / w;
956 } else {
957 const Double_t knorm=fgWSta->Eval(1e-4);
958 if(knorm) y = TMath::Sqrt(2. * rw * kDr / knorm / knorm);
65aa45f2 959 }
960 return y; //fm
961}
962
0c13ae32 963Double_t AliFastGlauber::WPathLength(const Double_t* x, const Double_t* par)
65aa45f2 964{
965 //
966 // Path Length as a function of phi
967 // Interaction point from random distribution
968 // as a function of the phi-direction
710a8d90 969 const Double_t kphi0 = x[0];
970 const Double_t kb = par[0];
971 fgWAlmond->SetParameter(0, kb);
972 const Int_t kNpi = Int_t (par[1]); //Number of interactions
973 const Int_t kiopt = Int_t(par[2]); //Path Length definition
65aa45f2 974
975 //
976 // r-steps
977 //
978 const Int_t kNp = 100;
979 const Double_t kDr = fgBMax/Double_t(kNp);
980 Double_t l = 0.; // Path length
981 for (Int_t in = 0; in < kNpi; in ++) {
f3a04204 982 Double_t rw = 0.;
983 Double_t w = 0.;
65aa45f2 984 // Interaction point
985 Double_t x0, y0;
986 fgWAlmond->GetRandom2(x0, y0);
987 // Initial radius
710a8d90 988 const Double_t kr0 = TMath::Sqrt(x0*x0 + y0*y0);
989 const Int_t knps = Int_t ((fgBMax - kr0)/kDr) - 1;
f3a04204 990
65aa45f2 991 // Radial steps
992 Double_t r = 0.;
710a8d90 993 for (Int_t i = 0; (i < knps ); i++) {
65aa45f2 994 // Transform into target frame
710a8d90 995 const Double_t kxx = x0 + r * TMath::Cos(kphi0) + kb / 2.;
996 const Double_t kyy = y0 + r * TMath::Sin(kphi0);
997 const Double_t kphi = TMath::ATan2(kyy, kxx);
998 const Double_t kr1 = TMath::Sqrt(kxx*kxx + kyy*kyy);
65aa45f2 999 // Radius in projectile frame
710a8d90 1000 const Double_t kr2 = TMath::Sqrt(kr1*kr1 + kb*kb - 2.*kr1*kb*TMath::Cos(kphi));
1001 const Double_t ky = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
65aa45f2 1002
710a8d90 1003 rw += ky * r;
1004 w += ky;
65aa45f2 1005 r += kDr;
1006 } // steps
1007 // Average over interactions
710a8d90 1008 if (!kiopt) {
65aa45f2 1009 if(w) l += (2. * rw / w);
8de7e046 1010 } else {
710a8d90 1011 const Double_t knorm=fgWSta->Eval(1e-4);
1012 if(knorm) l+= 2. * rw * kDr / knorm / knorm;
8de7e046 1013 }
65aa45f2 1014 } // interactions
1015 Double_t ret=0;
710a8d90 1016 if (!kiopt)
65aa45f2 1017 ret= l / kNpi;
1018 else
1019 ret=TMath::Sqrt( l / kNpi);
1020 return ret; //fm
f3a04204 1021}
1022
710a8d90 1023Double_t AliFastGlauber::CrossSection(Double_t b1, Double_t b2) const
f3a04204 1024{
65aa45f2 1025 //
1026 // Return the geometrical cross-section integrated from b1 to b2
1027 //
1028 return fgWSgeo->Integral(b1, b2)*10.; //mbarn
1029}
f3a04204 1030
710a8d90 1031Double_t AliFastGlauber::HardCrossSection(Double_t b1, Double_t b2) const
5b3a5a5d 1032{
65aa45f2 1033 //
1034 // Return the hard cross-section integrated from b1 to b2
1035 //
1036 return fgWSbinary->Integral(b1, b2)*10.; //mbarn
5b3a5a5d 1037}
1038
710a8d90 1039Double_t AliFastGlauber::FractionOfHardCrossSection(Double_t b1, Double_t b2) const
1bc228f5 1040{
65aa45f2 1041 //
148c5ce5 1042 // Return fraction of hard cross-section integrated from b1 to b2
65aa45f2 1043 //
1044 return fgWSbinary->Integral(b1, b2)/fgWSbinary->Integral(0., 100.);
1bc228f5 1045}
1046
0c13ae32 1047Double_t AliFastGlauber::NHard(const Double_t b1, const Double_t b2) const
f762082f 1048{
1049 //
1050 // Number of binary hard collisions
1051 // as a function of b (nucl/ex/0302016 eq. 19)
1052 //
1053 const Double_t kshard=HardCrossSection(b1,b2);
1054 const Double_t ksgeo=CrossSection(b1,b2);
1055 if(ksgeo>0)
1056 return kshard/ksgeo;
1057 else return -1;
1058}
1059
710a8d90 1060Double_t AliFastGlauber::Binaries(Double_t b) const
5b3a5a5d 1061{
65aa45f2 1062 //
710a8d90 1063 // Return number of binary hard collisions normalized to 1 at b=0
65aa45f2 1064 //
d206dc95 1065 if(b < 1.e-4) b = 1e-4;
65aa45f2 1066 return fgWSN->Eval(b)/fgWSN->Eval(1e-4);
5b3a5a5d 1067}
1068
a42548b0 1069Double_t AliFastGlauber::MeanOverlap(Double_t b1, Double_t b2)
1070{
1071//
1072// Calculate the mean overlap for impact parameter range b1 .. b2
1073//
1074 Double_t sum = 0.;
1075 Double_t sumc = 0.;
1076 Double_t b = b1;
1077
1078 while (b < b2-0.005) {
1079 Double_t nc = GetNumberOfCollisions(b);
1080 sum += 10. * fgWStaa->Eval(b) * fgWSgeo->Eval(b) * 0.01 / (1. - TMath::Exp(-nc));
1081 sumc += 10. * fgWSgeo->Eval(b) * 0.01;
1082 b += 0.01;
1083 }
1084 return (sum / CrossSection(b1, b2));
1085}
1086
1087
1088Double_t AliFastGlauber::MeanNumberOfCollisionsPerEvent(Double_t b1, Double_t b2)
1089{
1090//
1091// Calculate the mean number of collisions per event for impact parameter range b1 .. b2
1092//
1093 Double_t sum = 0.;
1094 Double_t sumc = 0.;
1095 Double_t b = b1;
1096
1097 while (b < b2-0.005) {
1098 Double_t nc = GetNumberOfCollisions(b);
1099 sum += nc / (1. - TMath::Exp(-nc)) * 10. * fgWSgeo->Eval(b) * 0.01;
1100 sumc += 10. * fgWSgeo->Eval(b) * 0.01;
1101 b += 0.01;
1102 }
1103 return (sum / CrossSection(b1, b2));
1104}
1105
1106
710a8d90 1107Double_t AliFastGlauber::GetNumberOfBinaries(Double_t b) const
5b3a5a5d 1108{
65aa45f2 1109 //
710a8d90 1110 // Return number of binary hard collisions at b
65aa45f2 1111 //
d206dc95 1112 if(b<1.e-4) b=1e-4;
65aa45f2 1113 return fgWSN->Eval(b);
5b3a5a5d 1114}
1115
710a8d90 1116Double_t AliFastGlauber::Participants(Double_t b) const
5b3a5a5d 1117{
65aa45f2 1118 //
1119 // Return the number of participants normalized to 1 at b=0
1120 //
d206dc95 1121 if(b<1.e-4) b=1e-4;
65aa45f2 1122 return (fgWParticipants->Eval(b)/fgWParticipants->Eval(1e-4));
5b3a5a5d 1123}
1124
710a8d90 1125Double_t AliFastGlauber::GetNumberOfParticipants(Double_t b) const
2a103154 1126{
65aa45f2 1127 //
1128 // Return the number of participants for impact parameter b
1129 //
d206dc95 1130 if(b<1.e-4) b=1e-4;
65aa45f2 1131 return (fgWParticipants->Eval(b));
2a103154 1132}
1133
710a8d90 1134Double_t AliFastGlauber::GetNumberOfCollisions(Double_t b) const
1135{
1136 //
1137 // Return the number of collisions for impact parameter b
1138 //
d206dc95 1139 if(b<1.e-4) b=1e-4;
710a8d90 1140 return (fgWStaa->Eval(b)*fSigmaNN);
1141}
1142
148c5ce5 1143Double_t AliFastGlauber::GetNumberOfCollisionsPerEvent(Double_t b) const
1144{
1145 //
1146 // Return the number of collisions per event (at least one collision)
1147 // for impact parameter b
1148 //
1149 Double_t n = GetNumberOfCollisions(b);
1150 if (n > 0.) {
1151 return (n / (1. - TMath::Exp(- n)));
1152 } else {
1153 return (0.);
1154 }
1155}
1156
5b3a5a5d 1157void AliFastGlauber::SimulateTrigger(Int_t n)
1158{
65aa45f2 1159 //
1160 // Simulates Trigger
1161 //
1162 TH1F* mbtH = new TH1F("mbtH", "MB Trigger b-Distribution", 100, 0., 20.);
1163 TH1F* hdtH = new TH1F("hdtH", "Hard Trigger b-Distribution", 100, 0., 20.);
1164 TH1F* mbmH = new TH1F("mbmH", "MB Trigger Multiplicity Distribution", 100, 0., 8000.);
1165 TH1F* hdmH = new TH1F("hdmH", "Hard Trigger Multiplicity Distribution", 100, 0., 8000.);
5b3a5a5d 1166
65aa45f2 1167 mbtH->SetXTitle("b [fm]");
1168 hdtH->SetXTitle("b [fm]");
1169 mbmH->SetXTitle("Multiplicity");
1170 hdmH->SetXTitle("Multiplicity");
5b3a5a5d 1171
65aa45f2 1172 TCanvas *c0 = new TCanvas("c0","Trigger Simulation",400,10,600,700);
1173 c0->Divide(2,1);
1174 TCanvas *c1 = new TCanvas("c1","Trigger Simulation",400,10,600,700);
1175 c1->Divide(1,2);
5b3a5a5d 1176
65aa45f2 1177 //
1178 //
1179 Init(1);
1180 for (Int_t iev = 0; iev < n; iev++)
5b3a5a5d 1181 {
65aa45f2 1182 Float_t b, p, mult;
1183 GetRandom(b, p, mult);
1184 mbtH->Fill(b,1.);
1185 hdtH->Fill(b, p);
1186 mbmH->Fill(mult, 1.);
1187 hdmH->Fill(mult, p);
1188
1189 c0->cd(1);
1190 mbtH->Draw();
1191 c0->cd(2);
1192 hdtH->Draw();
1193 c0->Update();
1194
1195 c1->cd(1);
1196 mbmH->Draw();
1197 c1->cd(2);
1198 hdmH->Draw();
1199 c1->Update();
5b3a5a5d 1200 }
1201}
1202
1203void AliFastGlauber::GetRandom(Float_t& b, Float_t& p, Float_t& mult)
1204{
65aa45f2 1205 //
1206 // Gives back a random impact parameter, hard trigger probability and multiplicity
1207 //
1208 b = fgWSgeo->GetRandom();
710a8d90 1209 const Float_t kmu = fgWSN->Eval(b);
1210 p = 1.-TMath::Exp(-kmu);
1211 mult = 6000./fgWSN->Eval(1.) * kmu;
5b3a5a5d 1212}
1213
c2140715 1214void AliFastGlauber::GetRandom(Int_t& bin, Bool_t& hard)
1215{
65aa45f2 1216 //
1217 // Gives back a random impact parameter bin, and hard trigger decission
1218 //
710a8d90 1219 const Float_t kb = fgWSgeo->GetRandom();
1220 const Float_t kmu = fgWSN->Eval(kb) * fSigmaHard;
1221 const Float_t kp = 1.-TMath::Exp(-kmu);
1222 if (kb < 5.) {
65aa45f2 1223 bin = 1;
710a8d90 1224 } else if (kb < 8.6) {
65aa45f2 1225 bin = 2;
710a8d90 1226 } else if (kb < 11.2) {
65aa45f2 1227 bin = 3;
710a8d90 1228 } else if (kb < 13.2) {
65aa45f2 1229 bin = 4;
710a8d90 1230 } else if (kb < 15.0) {
65aa45f2 1231 bin = 5;
1232 } else {
1233 bin = 6;
1234 }
1235 hard = kFALSE;
710a8d90 1236 const Float_t kr = gRandom->Rndm();
1237 if (kr < kp) hard = kTRUE;
c2140715 1238}
1239
65aa45f2 1240Double_t AliFastGlauber::GetRandomImpactParameter(Double_t bmin, Double_t bmax)
5b3a5a5d 1241{
65aa45f2 1242 //
1243 // Gives back a random impact parameter in the range bmin .. bmax
1244 //
1245 Float_t b = -1.;
1246 while(b < bmin || b > bmax)
1247 b = fgWSgeo->GetRandom();
1248 return b;
5b3a5a5d 1249}
1250
710a8d90 1251void AliFastGlauber::StoreFunctions() const
5b3a5a5d 1252{
65aa45f2 1253 //
1254 // Store in file functions
1255 //
1256 TFile* ff = new TFile(fName.Data(),"recreate");
1257 fgWStaa->Write("WStaa");
1258 fgWParticipants->Write("WParticipants");
1259 ff->Close();
1260 return;
5b3a5a5d 1261}
1262
65aa45f2 1263//=================== Added by A. Dainese 11/02/04 ===========================
1bc228f5 1264
710a8d90 1265void AliFastGlauber::StoreAlmonds() const
5b3a5a5d 1266{
65aa45f2 1267 //
1268 // Store in file
1269 // 40 almonds for b = (0.25+k*0.5) fm (k=0->39)
1270 //
1271 Char_t almondName[100];
1272 TFile* ff = new TFile(fName.Data(),"update");
1273 for(Int_t k=0; k<40; k++) {
bae6da3c 1274 snprintf(almondName, 100, "WAlmondFixedB%d",k);
65aa45f2 1275 Double_t b = 0.25+k*0.5;
1276 Info("StoreAlmonds"," b = %f\n",b);
1277 fgWAlmond->SetParameter(0,b);
1278 fgWAlmond->Write(almondName);
1279 }
1280 ff->Close();
1281 return;
5b3a5a5d 1282}
a2f2f511 1283
a2f2f511 1284void AliFastGlauber::SetCentralityClass(Double_t xsecFrLow,Double_t xsecFrUp)
1285{
1286 //
65aa45f2 1287 // Set limits of centrality class as fractions
1288 // of the geomtrical cross section
a2f2f511 1289 //
1290 if(xsecFrLow>1. || xsecFrUp>1. || xsecFrLow>xsecFrUp) {
65aa45f2 1291 Error("SetCentralityClass", "Please set 0 <= xsecFrLow <= xsecFrUp <= 1\n");
a2f2f511 1292 return;
1293 }
1294
1295 Double_t bLow=0.,bUp=0.;
1296 Double_t xsecFr=0.;
710a8d90 1297 const Double_t knorm=fgWSgeo->Integral(0.,100.);
a2f2f511 1298 while(xsecFr<xsecFrLow) {
710a8d90 1299 xsecFr = fgWSgeo->Integral(0.,bLow)/knorm;
a2f2f511 1300 bLow += 0.1;
1301 }
1302 bUp = bLow;
1303 while(xsecFr<xsecFrUp) {
710a8d90 1304 xsecFr = fgWSgeo->Integral(0.,bUp)/knorm;
a2f2f511 1305 bUp += 0.1;
1306 }
1307
710a8d90 1308 Info("SetCentralityClass", "Centrality class: %4.2f-%4.2f; %4.1f < b < %4.1f fm",
a2f2f511 1309 xsecFrLow,xsecFrUp,bLow,bUp);
a2f2f511 1310 fgWSbinary->SetRange(bLow,bUp);
710a8d90 1311 fBmin=bLow;
1312 fBmax=bUp;
a2f2f511 1313 return;
1314}
1315
1316void AliFastGlauber::GetRandomBHard(Double_t& b)
1317{
1318 //
1319 // Get random impact parameter according to distribution of
1320 // hard (binary) cross-section, in the range defined by the centrality class
1321 //
1322 b = fgWSbinary->GetRandom();
1323 Int_t bin = 2*(Int_t)b;
1324 if( (b-(Int_t)b) > 0.5) bin++;
7f2f270b 1325 fgWAlmondCurrent = fgWAlmondFixedB[bin];
a2f2f511 1326 return;
1327}
1328
1329void AliFastGlauber::GetRandomXY(Double_t& x,Double_t& y)
1330{
1331 //
1332 // Get random position of parton production point according to
1333 // product of thickness functions
1334 //
65aa45f2 1335 fgWAlmondCurrent->GetRandom2(x,y);
a2f2f511 1336 return;
1337}
1338
1339void AliFastGlauber::GetRandomPhi(Double_t& phi)
1340{
1341 //
1342 // Get random parton azimuthal propagation direction
1343 //
1344 phi = 2.*TMath::Pi()*gRandom->Rndm();
1345 return;
1346}
1347
65aa45f2 1348Double_t AliFastGlauber::CalculateLength(Double_t b,Double_t x0,Double_t y0,Double_t phi0)
a2f2f511 1349{
1350 //
1351 // Calculate path length for a parton with production point (x0,y0)
1352 // and propagation direction (ux=cos(phi0),uy=sin(phi0))
1353 // in a collision with impact parameter b
1354 //
1355
1356 // number of steps in l
1357 const Int_t kNp = 100;
1358 const Double_t kDl = fgBMax/Double_t(kNp);
1359
a2f2f511 1360 if(fEllDef==1) {
1361 //
1362 // Definition 1:
1363 //
65aa45f2 1364 // ell = 2 * \int_0^\infty dl*l*(T_A*T_B)(x0+l*ux,y0+l*uy) /
1365 // \int_0^\infty dl*(T_A*T_B)(x0+l*ux,y0+l*uy)
a2f2f511 1366 //
a2f2f511 1367
a2f2f511 1368 // Initial radius
710a8d90 1369 const Double_t kr0 = TMath::Sqrt(x0*x0 + y0*y0);
1370 const Int_t knps = Int_t ((fgBMax - kr0)/kDl) - 1;
65aa45f2 1371 Double_t l = 0.;
1372 Double_t integral1 = 0.;
1373 Double_t integral2 = 0.;
a2f2f511 1374 // Radial steps
710a8d90 1375 for (Int_t i = 0; i < knps; i++) {
a2f2f511 1376
1377 // Transform into target frame
710a8d90 1378 const Double_t kxx = x0 + l * TMath::Cos(phi0) + b / 2.;
1379 const Double_t kyy = y0 + l * TMath::Sin(phi0);
1380 const Double_t kphi = TMath::ATan2(kyy, kxx);
1381 const Double_t kr1 = TMath::Sqrt(kxx*kxx + kyy*kyy);
a2f2f511 1382 // Radius in projectile frame
710a8d90 1383 const Double_t kr2 = TMath::Sqrt(kr1*kr1 + b*b - 2.*kr1*b*TMath::Cos(kphi));
1384 const Double_t kprodTATB = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
a2f2f511 1385
710a8d90 1386 integral1 += kprodTATB * l * kDl;
1387 integral2 += kprodTATB * kDl;
a2f2f511 1388 l += kDl;
1389 } // steps
1390
65aa45f2 1391 Double_t ell=0.;
1392 if(integral2)
1393 ell = (2. * integral1 / integral2);
a2f2f511 1394 return ell;
a2f2f511 1395 } else if(fEllDef==2) {
1396 //
1397 // Definition 2:
1398 //
1399 // ell = \int_0^\infty dl*
65aa45f2 1400 // \Theta((T_A*T_B)(x0+l*ux,y0+l*uy)-0.5*(T_A*T_B)(0,0))
a2f2f511 1401 //
1402
a2f2f511 1403 // Initial radius
710a8d90 1404 const Double_t kr0 = TMath::Sqrt(x0*x0 + y0*y0);
1405 const Int_t knps = Int_t ((fgBMax - kr0)/kDl) - 1;
1406 const Double_t kprodTATBHalfMax = 0.5*fgWAlmondCurrent->Eval(0.,0.);
a2f2f511 1407 // Radial steps
65aa45f2 1408 Double_t l = 0.;
1409 Double_t integral = 0.;
710a8d90 1410 for (Int_t i = 0; i < knps; i++) {
a2f2f511 1411 // Transform into target frame
710a8d90 1412 const Double_t kxx = x0 + l * TMath::Cos(phi0) + b / 2.;
1413 const Double_t kyy = y0 + l * TMath::Sin(phi0);
1414 const Double_t kphi = TMath::ATan2(kyy, kxx);
1415 const Double_t kr1 = TMath::Sqrt(kxx*kxx + kyy*kyy);
a2f2f511 1416 // Radius in projectile frame
710a8d90 1417 const Double_t kr2 = TMath::Sqrt(kr1*kr1 + b*b - 2.*kr1*b*TMath::Cos(kphi));
1418 const Double_t kprodTATB = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
1419 if(kprodTATB>kprodTATBHalfMax) integral += kDl;
a2f2f511 1420 l += kDl;
1421 } // steps
65aa45f2 1422 Double_t ell = integral;
a2f2f511 1423 return ell;
a2f2f511 1424 } else {
65aa45f2 1425 Error("CalculateLength","Wrong length definition setting: %d !\n",fEllDef);
a2f2f511 1426 return -1.;
1427 }
1428}
1429
83f67d08 1430void AliFastGlauber::GetLengthAndPhi(Double_t& ell,Double_t& phi,Double_t b)
a2f2f511 1431{
1432 //
1433 // Return length from random b, x0, y0, phi0
83f67d08 1434 // Return also phi0
a2f2f511 1435 //
1436 Double_t x0,y0,phi0;
1437 if(b<0.) GetRandomBHard(b);
1438 GetRandomXY(x0,y0);
1439 GetRandomPhi(phi0);
83f67d08 1440 phi = phi0;
a2f2f511 1441 ell = CalculateLength(b,x0,y0,phi0);
a2f2f511 1442 return;
1443}
1444
83f67d08 1445void AliFastGlauber::GetLength(Double_t& ell,Double_t b)
a2f2f511 1446{
1447 //
83f67d08 1448 // Return length from random b, x0, y0, phi0
a2f2f511 1449 //
83f67d08 1450 Double_t phi;
1451 GetLengthAndPhi(ell,phi,b);
1452 return;
1453}
1454
1455void AliFastGlauber::GetLengthsBackToBackAndPhi(Double_t& ell1,Double_t& ell2,Double_t &phi,Double_t b)
1456{
1457 //
1458 // Return 2 lengths back to back from random b, x0, y0, phi0
1459 // Return also phi0
1460 //
a2f2f511 1461 Double_t x0,y0,phi0;
1462 if(b<0.) GetRandomBHard(b);
1463 GetRandomXY(x0,y0);
1464 GetRandomPhi(phi0);
710a8d90 1465 const Double_t kphi0plusPi = phi0+TMath::Pi();
83f67d08 1466 phi = phi0;
a2f2f511 1467 ell1 = CalculateLength(b,x0,y0,phi0);
710a8d90 1468 ell2 = CalculateLength(b,x0,y0,kphi0plusPi);
a2f2f511 1469 return;
1470}
1471
83f67d08 1472void AliFastGlauber::GetLengthsBackToBack(Double_t& ell1,Double_t& ell2,
1473 Double_t b)
1474{
1475 //
1476 // Return 2 lengths back to back from random b, x0, y0, phi0
1477 //
1478 Double_t phi;
1479 GetLengthsBackToBackAndPhi(ell1,ell2,phi,b);
1480 return;
1481}
1482
620294e8 1483void AliFastGlauber::GetLengthsForPythia(Int_t n,Double_t* const phi,Double_t* ell, Double_t b)
a2f2f511 1484{
1485 //
1486 // Returns lenghts for n partons with azimuthal angles phi[n]
1487 // from random b, x0, y0
1488 //
1bc228f5 1489 Double_t x0, y0;
1490 if(b < 0.) GetRandomBHard(b);
a2f2f511 1491 GetRandomXY(x0,y0);
1bc228f5 1492 for(Int_t i = 0; i< n; i++) ell[i] = CalculateLength(b,x0,y0,phi[i]);
a2f2f511 1493 return;
1494}
1495
1496void AliFastGlauber::PlotBDistr(Int_t n)
65aa45f2 1497{
a2f2f511 1498 //
1499 // Plot distribution of n impact parameters
1500 //
1501 Double_t b;
1502 TH1F *hB = new TH1F("hB","dN/db",100,0,fgBMax);
1503 hB->SetXTitle("b [fm]");
1504 hB->SetYTitle("dN/db [a.u.]");
1505 hB->SetFillColor(3);
a2f2f511 1506 for(Int_t i=0; i<n; i++) {
1507 GetRandomBHard(b);
1508 hB->Fill(b);
1509 }
a2f2f511 1510 TCanvas *cB = new TCanvas("cB","Impact parameter distribution",0,0,500,500);
1511 cB->cd();
1512 hB->Draw();
a2f2f511 1513 return;
1514}
1515
d3d4a92f 1516void AliFastGlauber::PlotLengthDistr(Int_t n,Bool_t save,const char *fname)
a2f2f511 1517{
1518 //
1519 // Plot length distribution
1520 //
1521 Double_t ell;
710a8d90 1522 TH1F *hEll = new TH1F("hEll","Length distribution",64,-0.5,15);
a2f2f511 1523 hEll->SetXTitle("Transverse path length, L [fm]");
1524 hEll->SetYTitle("Probability");
1525 hEll->SetFillColor(2);
a2f2f511 1526 for(Int_t i=0; i<n; i++) {
1527 GetLength(ell);
1528 hEll->Fill(ell);
1529 }
1530 hEll->Scale(1/(Double_t)n);
a2f2f511 1531 TCanvas *cL = new TCanvas("cL","Length distribution",0,0,500,500);
1532 cL->cd();
1533 hEll->Draw();
1534
1535 if(save) {
1536 TFile *f = new TFile(fname,"recreate");
1537 hEll->Write();
1538 f->Close();
1539 }
1540 return;
1541}
1542
d3d4a92f 1543void AliFastGlauber::PlotLengthB2BDistr(Int_t n,Bool_t save,const char *fname)
a2f2f511 1544{
1545 //
1546 // Plot lengths back-to-back distributions
1547 //
1548 Double_t ell1,ell2;
1549 TH2F *hElls = new TH2F("hElls","Lengths back-to-back",100,0,15,100,0,15);
1550 hElls->SetXTitle("Transverse path length, L [fm]");
1551 hElls->SetYTitle("Transverse path length, L [fm]");
a2f2f511 1552 for(Int_t i=0; i<n; i++) {
1553 GetLengthsBackToBack(ell1,ell2);
1554 hElls->Fill(ell1,ell2);
1555 }
1556 hElls->Scale(1/(Double_t)n);
a2f2f511 1557 TCanvas *cLs = new TCanvas("cLs","Length back-to-back distribution",0,0,500,500);
1558 gStyle->SetPalette(1,0);
1559 cLs->cd();
1560 hElls->Draw("col,Z");
a2f2f511 1561 if(save) {
1562 TFile *f = new TFile(fname,"recreate");
1563 hElls->Write();
1564 f->Close();
1565 }
1566 return;
1567}
1568
710a8d90 1569void AliFastGlauber::PlotAlmonds() const
1570{
1571 //
1572 // Plot almonds for some impact parameters
1573 //
1574 TCanvas *c = new TCanvas("c","Almonds",0,0,500,500);
1575 gStyle->SetPalette(1,0);
1576 c->Divide(2,2);
1577 c->cd(1);
7f2f270b 1578 fgWAlmondFixedB[0]->Draw("cont1");
710a8d90 1579 c->cd(2);
7f2f270b 1580 fgWAlmondFixedB[10]->Draw("cont1");
710a8d90 1581 c->cd(3);
7f2f270b 1582 fgWAlmondFixedB[20]->Draw("cont1");
710a8d90 1583 c->cd(4);
7f2f270b 1584 fgWAlmondFixedB[30]->Draw("cont1");
710a8d90 1585 return;
1586}
1587
65aa45f2 1588//=================== Added by A. Dainese 05/03/04 ===========================
1589
1590void AliFastGlauber::CalculateI0I1(Double_t& integral0,Double_t& integral1,
1591 Double_t b,Double_t x0,Double_t y0,
710a8d90 1592 Double_t phi0,Double_t ellCut) const
a2f2f511 1593{
65aa45f2 1594 //
1595 // Calculate integrals:
e9663638 1596 // integral0 = \int_0^ellCut dl*(T_A*T_B)(x0+l*ux,y0+l*uy)
1597 // integral1 = \int_0^ellCut dl*l*(T_A*T_B)(x0+l*ux,y0+l*uy)
a2f2f511 1598 //
65aa45f2 1599 // for a parton with production point (x0,y0)
1600 // and propagation direction (ux=cos(phi0),uy=sin(phi0))
1601 // in a collision with impact parameter b
1602 //
1603
1604 // number of steps in l
1605 const Int_t kNp = 100;
1606 const Double_t kDl = fgBMax/Double_t(kNp);
1607
1608 // Initial radius
710a8d90 1609 const Double_t kr0 = TMath::Sqrt(x0 * x0 + y0 * y0);
1610 const Int_t knps = Int_t ((fgBMax - kr0)/kDl) - 1;
65aa45f2 1611
1612 // Radial steps
1613 Double_t l = 0.;
1614 integral0 = 0.;
1615 integral1 = 0.;
1616 Int_t i = 0;
710a8d90 1617 while((i < knps) && (l < ellCut)) {
65aa45f2 1618 // Transform into target frame
710a8d90 1619 const Double_t kxx = x0 + l * TMath::Cos(phi0) + b / 2.;
1620 const Double_t kyy = y0 + l * TMath::Sin(phi0);
1621 const Double_t kphi = TMath::ATan2(kyy, kxx);
1622 const Double_t kr1 = TMath::Sqrt(kxx*kxx + kyy*kyy);
65aa45f2 1623 // Radius in projectile frame
710a8d90 1624 const Double_t kr2 = TMath::Sqrt(kr1*kr1 + b*b - 2.*kr1*b*TMath::Cos(kphi));
1625 const Double_t kprodTATB = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
1626 integral0 += kprodTATB * kDl;
1627 integral1 += kprodTATB * l * kDl;
65aa45f2 1628 l += kDl;
1629 i++;
1630 } // steps
1631 return;
1632}
1633
83f67d08 1634void AliFastGlauber::GetI0I1AndPhi(Double_t& integral0,Double_t& integral1,
1635 Double_t& phi,
1636 Double_t ellCut,Double_t b)
65aa45f2 1637{
a2f2f511 1638 //
65aa45f2 1639 // Return I0 and I1 from random b, x0, y0, phi0
83f67d08 1640 // Return also phi
65aa45f2 1641 //
1642 Double_t x0,y0,phi0;
1643 if(b<0.) GetRandomBHard(b);
1644 GetRandomXY(x0,y0);
1645 GetRandomPhi(phi0);
83f67d08 1646 phi = phi0;
65aa45f2 1647 CalculateI0I1(integral0,integral1,b,x0,y0,phi0,ellCut);
1648 return;
1649}
a2f2f511 1650
83f67d08 1651void AliFastGlauber::GetI0I1(Double_t& integral0,Double_t& integral1,
1652 Double_t ellCut,Double_t b)
1653{
1654 //
1655 // Return I0 and I1 from random b, x0, y0, phi0
1656 //
1657 Double_t phi;
1658 GetI0I1AndPhi(integral0,integral1,phi,ellCut,b);
1659 return;
1660}
1661
1662void AliFastGlauber::GetI0I1BackToBackAndPhi(Double_t& integral01,Double_t& integral11,
1663 Double_t& integral02,Double_t& integral12,
1664 Double_t& phi,
1665 Double_t ellCut,Double_t b)
65aa45f2 1666{
1667 //
1668 // Return 2 pairs of I0 and I1 back to back from random b, x0, y0, phi0
83f67d08 1669 // Return also phi0
65aa45f2 1670 //
1671 Double_t x0,y0,phi0;
1672 if(b<0.) GetRandomBHard(b);
1673 GetRandomXY(x0,y0);
1674 GetRandomPhi(phi0);
83f67d08 1675 phi = phi0;
710a8d90 1676 const Double_t kphi0plusPi = phi0+TMath::Pi();
c54404bf 1677 CalculateI0I1(integral01,integral11,b,x0,y0,phi0,ellCut);
1678 CalculateI0I1(integral02,integral12,b,x0,y0,kphi0plusPi,ellCut);
1679 return;
1680}
1681
1682void AliFastGlauber::GetI0I1BackToBackAndPhiAndXY(Double_t& integral01,Double_t& integral11,
1683 Double_t& integral02,Double_t& integral12,
1684 Double_t& phi,Double_t &x,Double_t &y,
1685 Double_t ellCut,Double_t b)
1686{
1687 //
1688 // Return 2 pairs of I0 and I1 back to back from random b, x0, y0, phi0
1689 // Return also phi0
1690 //
1691 Double_t x0,y0,phi0;
1692 if(b<0.) GetRandomBHard(b);
1693 GetRandomXY(x0,y0);
1694 GetRandomPhi(phi0);
1695 phi = phi0; x=x0; y=y0;
1696 const Double_t kphi0plusPi = phi0+TMath::Pi();
65aa45f2 1697 CalculateI0I1(integral01,integral11,b,x0,y0,phi0,ellCut);
710a8d90 1698 CalculateI0I1(integral02,integral12,b,x0,y0,kphi0plusPi,ellCut);
65aa45f2 1699 return;
1700}
a2f2f511 1701
83f67d08 1702void AliFastGlauber::GetI0I1BackToBack(Double_t& integral01,Double_t& integral11,
1703 Double_t& integral02,Double_t& integral12,
1704 Double_t ellCut,Double_t b)
1705{
1706 //
1707 // Return 2 pairs of I0 and I1 back to back from random b, x0, y0, phi0
1708 //
1709 Double_t phi;
1710 GetI0I1BackToBackAndPhi(integral01,integral11,integral02,integral12,
1711 phi,ellCut,b);
1712 return;
1713}
1714
65aa45f2 1715void AliFastGlauber::GetI0I1ForPythia(Int_t n,Double_t* phi,
1716 Double_t* integral0,Double_t* integral1,
1717 Double_t ellCut,Double_t b)
1718{
1719 //
1720 // Returns I0 and I1 pairs for n partons with azimuthal angles phi[n]
1721 // from random b, x0, y0
1722 //
1723 Double_t x0,y0;
1724 if(b<0.) GetRandomBHard(b);
1725 GetRandomXY(x0,y0);
1726 for(Int_t i=0; i<n; i++)
1727 CalculateI0I1(integral0[i],integral1[i],b,x0,y0,phi[i],ellCut);
1728 return;
1729}
1730
2e3b5c95 1731void AliFastGlauber::GetI0I1ForPythiaAndXY(Int_t n,Double_t* phi,
1732 Double_t* integral0,Double_t* integral1,
1733 Double_t &x,Double_t& y,
1734 Double_t ellCut,Double_t b)
1735{
1736 //
1737 // Returns I0 and I1 pairs for n partons with azimuthal angles phi[n]
1738 // from random b, x0, y0 and return x0,y0
1739 //
1740 Double_t x0,y0;
1741 if(b<0.) GetRandomBHard(b);
1742 GetRandomXY(x0,y0);
1743 for(Int_t i=0; i<n; i++)
1744 CalculateI0I1(integral0[i],integral1[i],b,x0,y0,phi[i],ellCut);
1745 x=x0;
1746 y=y0;
1747 return;
1748}
1749
65aa45f2 1750void AliFastGlauber::PlotI0I1Distr(Int_t n,Double_t ellCut,
d3d4a92f 1751 Bool_t save,const char *fname)
65aa45f2 1752{
1753 //
bbf8513d 1754 // Plot I0-I1 distribution
65aa45f2 1755 //
1756 Double_t i0,i1;
bbf8513d 1757 TH2F *hI0I1s = new TH2F("hI0I1s","I_{0} versus I_{1}",1000,0,0.001,1000,0,0.01);
1758 hI0I1s->SetXTitle("I_{0} [fm^{-3}]");
1759 hI0I1s->SetYTitle("I_{1} [fm^{-2}]");
1760
65aa45f2 1761 TH1F *hI0 = new TH1F("hI0","I_{0} = #hat{q}L / k",
bbf8513d 1762 1000,0,0.001);
65aa45f2 1763 hI0->SetXTitle("I_{0} [fm^{-3}]");
1764 hI0->SetYTitle("Probability");
1765 hI0->SetFillColor(3);
1766 TH1F *hI1 = new TH1F("hI1","I_{1} = #omega_{c} / k",
bbf8513d 1767 1000,0,0.01);
65aa45f2 1768 hI1->SetXTitle("I_{1} [fm^{-2}]");
1769 hI1->SetYTitle("Probability");
1770 hI1->SetFillColor(4);
1771 TH1F *h2 = new TH1F("h2","2 I_{1}^{2}/I_{0} = R / k",
1772 100,0,0.02);
1773 h2->SetXTitle("2 I_{1}^{2}/I_{0} [fm^{-1}]");
1774 h2->SetYTitle("Probability");
1775 h2->SetFillColor(2);
1776 TH1F *h3 = new TH1F("h3","2 I_{1}/I_{0} = L",
1777 100,0,15);
1778 h3->SetXTitle("2 I_{1}/I_{0} [fm]");
1779 h3->SetYTitle("Probability");
1780 h3->SetFillColor(5);
1781 TH1F *h4 = new TH1F("h4","I_{0}^{2}/(2 I_{1}) = #hat{q} / k",
1782 100,0,0.00015);
1783 h4->SetXTitle("I_{0}^{2}/(2 I_{1}) [fm^{-4}]");
1784 h4->SetYTitle("Probability");
1785 h4->SetFillColor(7);
1786
1787 for(Int_t i=0; i<n; i++) {
1788 GetI0I1(i0,i1,ellCut);
bbf8513d 1789 hI0I1s->Fill(i0,i1);
65aa45f2 1790 hI0->Fill(i0);
1791 hI1->Fill(i1);
1792 h2->Fill(2.*i1*i1/i0);
1793 h3->Fill(2.*i1/i0);
1794 h4->Fill(i0*i0/2./i1);
a2f2f511 1795 }
65aa45f2 1796 hI0->Scale(1/(Double_t)n);
1797 hI1->Scale(1/(Double_t)n);
bbf8513d 1798 h2->Scale(1/(Double_t)n);
1799 h3->Scale(1/(Double_t)n);
1800 h4->Scale(1/(Double_t)n);
1801 hI0I1s->Scale(1/(Double_t)n);
65aa45f2 1802
1803 TCanvas *cI0I1 = new TCanvas("cI0I1","I0 and I1",0,0,900,700);
1804 cI0I1->Divide(3,2);
1805 cI0I1->cd(1);
1806 hI0->Draw();
1807 cI0I1->cd(2);
1808 hI1->Draw();
1809 cI0I1->cd(3);
1810 h2->Draw();
1811 cI0I1->cd(4);
1812 h3->Draw();
1813 cI0I1->cd(5);
1814 h4->Draw();
bbf8513d 1815 cI0I1->cd(6);
1816 gStyle->SetPalette(1,0);
1817 hI0I1s->Draw("col,Z");
a2f2f511 1818
65aa45f2 1819 if(save) {
1820 TFile *f = new TFile(fname,"recreate");
bbf8513d 1821 hI0I1s->Write();
65aa45f2 1822 hI0->Write();
1823 hI1->Write();
1824 h2->Write();
1825 h3->Write();
1826 h4->Write();
1827 f->Close();
1828 }
a2f2f511 1829 return;
1830}
1831
65aa45f2 1832void AliFastGlauber::PlotI0I1B2BDistr(Int_t n,Double_t ellCut,
d3d4a92f 1833 Bool_t save,const char *fname)
a2f2f511 1834{
1835 //
bbf8513d 1836 // Plot I0-I1 back-to-back distributions
a2f2f511 1837 //
65aa45f2 1838 Double_t i01,i11,i02,i12;
1839 TH2F *hI0s = new TH2F("hI0s","I_{0}'s back-to-back",100,0,100,100,0,100);
1840 hI0s->SetXTitle("I_{0} [fm^{-3}]");
1841 hI0s->SetYTitle("I_{0} [fm^{-3}]");
1842 TH2F *hI1s = new TH2F("hI1s","I_{1}'s back-to-back",100,0,100,100,0,100);
1843 hI1s->SetXTitle("I_{1} [fm^{-2}]");
1844 hI1s->SetYTitle("I_{1} [fm^{-2}]");
a2f2f511 1845
65aa45f2 1846 for(Int_t i=0; i<n; i++) {
1847 GetI0I1BackToBack(i01,i11,i02,i12,ellCut);
1848 hI0s->Fill(i01,i02);
1849 hI1s->Fill(i11,i12);
1850 }
1851 hI0s->Scale(1/(Double_t)n);
1852 hI1s->Scale(1/(Double_t)n);
1853
1854 TCanvas *cI0I1s = new TCanvas("cI0I1s","I0 and I1 back-to-back distributions",0,0,800,400);
a2f2f511 1855 gStyle->SetPalette(1,0);
65aa45f2 1856 cI0I1s->Divide(2,1);
1857 cI0I1s->cd(1);
1858 hI0s->Draw("col,Z");
1859 cI0I1s->cd(2);
1860 hI1s->Draw("col,Z");
a2f2f511 1861
65aa45f2 1862 if(save) {
1863 TFile *f = new TFile(fname,"recreate");
1864 hI0s->Write();
1865 hI1s->Write();
1866 f->Close();
1867 }
a2f2f511 1868 return;
1869}
710a8d90 1870
a42548b0 1871AliFastGlauber& AliFastGlauber::operator=(const AliFastGlauber& rhs)
1872{
1873// Assignment operator
1874 rhs.Copy(*this);
1875 return *this;
1876}
1877
1878void AliFastGlauber::Copy(TObject&) const
1879{
1880 //
1881 // Copy
1882 //
1883 Fatal("Copy","Not implemented!\n");
1884}
1885