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