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