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