1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 // Utility class to make simple Glauber type calculations
19 // for SYMMTRIC collision geometries (AA):
20 // Impact parameter, production points, reaction plane dependence
22 // The SimulateTrigger method can be used for simple MB and hard-process
23 // (binary scaling) trigger studies.
25 // Some basic quantities can be visualized directly.
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.
30 // ***** If you change settings dont forget to call init afterwards, *****
31 // ***** in order to update the formulas with the new parameters. *****
33 // Author: andreas.morsch@cern.ch
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
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
51 #include "AliFastGlauber.h"
63 #include <Riostream.h>
65 ClassImp(AliFastGlauber)
67 Float_t AliFastGlauber::fgBMax = 0.;
68 TF1* AliFastGlauber::fgWSb = NULL;
69 TF2* AliFastGlauber::fgWSbz = NULL;
70 TF1* AliFastGlauber::fgWSz = NULL;
71 TF1* AliFastGlauber::fgWSta = NULL;
72 TF2* AliFastGlauber::fgWStarfi = NULL;
73 TF2* AliFastGlauber::fgWAlmond = NULL;
74 TF1* AliFastGlauber::fgWStaa = NULL;
75 TF1* AliFastGlauber::fgWSgeo = NULL;
76 TF1* AliFastGlauber::fgWSbinary = NULL;
77 TF1* AliFastGlauber::fgWSN = NULL;
78 TF1* AliFastGlauber::fgWPathLength0 = NULL;
79 TF1* AliFastGlauber::fgWPathLength = NULL;
80 TF1* AliFastGlauber::fgWEnergyDensity = NULL;
81 TF1* AliFastGlauber::fgWIntRadius = NULL;
82 TF2* AliFastGlauber::fgWKParticipants = NULL;
83 TF1* AliFastGlauber::fgWParticipants = NULL;
84 TF2* AliFastGlauber::fgWAlmondCurrent = NULL;
85 TF2 AliFastGlauber::fgWAlmondFixedB[40];
86 const Int_t AliFastGlauber::fgkMCInts = 100000;
87 Int_t AliFastGlauber::fgCounter = 0;
89 AliFastGlauber::AliFastGlauber() : fName()
91 // Default Constructor
95 Error("AliFastGlauber","More than more instance (%d) is not supported, check your code!",fgCounter);
99 SetLengthDefinition();
103 AliFastGlauber::~AliFastGlauber()
106 //if(fgCounter==0) Reset();
109 void AliFastGlauber::SetAuAuRhic()
111 //Set all parameters for RHIC
112 SetWoodSaxonParametersAu();
113 SetHardCrossSection();
114 SetNNCrossSection(42);
116 SetFileName("$(ALICE_ROOT)/FASTSIM/data/glauberAuAu.root");
119 void AliFastGlauber::SetPbPbLHC()
121 //Set all parameters for LHC
122 SetWoodSaxonParametersPb();
123 SetHardCrossSection();
129 void AliFastGlauber::Init(Int_t mode)
132 // mode = 0; all functions are calculated
133 // mode = 1; overlap function is read from file (for Pb-Pb only)
134 // mode = 2; interaction almond functions are read from file
135 // USE THIS FOR PATH LENGTH CALC.!
145 fgWSb = new TF1("WSb", WSb, 0, fgBMax, 4);
146 fgWSb->SetParameter(0, fWSr0);
147 fgWSb->SetParameter(1, fWSd);
148 fgWSb->SetParameter(2, fWSw);
149 fgWSb->SetParameter(3, fWSn);
151 fgWSbz = new TF2("WSbz", WSbz, 0, fgBMax, 4);
152 fgWSbz->SetParameter(0, fWSr0);
153 fgWSbz->SetParameter(1, fWSd);
154 fgWSbz->SetParameter(2, fWSw);
155 fgWSbz->SetParameter(3, fWSn);
157 fgWSz = new TF1("WSz", WSz, 0, fgBMax, 5);
158 fgWSz->SetParameter(0, fWSr0);
159 fgWSz->SetParameter(1, fWSd);
160 fgWSz->SetParameter(2, fWSw);
161 fgWSz->SetParameter(3, fWSn);
166 fgWSta = new TF1("WSta", WSta, 0., fgBMax, 0);
171 fgWStarfi = new TF2("WStarfi", WStarfi, 0., fgBMax, 0., TMath::Pi(), 1);
172 fgWStarfi->SetParameter(0, 0.);
173 fgWStarfi->SetNpx(200);
174 fgWStarfi->SetNpy(20);
177 // Participants Kernel
179 fgWKParticipants = new TF2("WKParticipants", WKParticipants, 0., fgBMax, 0., TMath::Pi(), 3);
180 fgWKParticipants->SetParameter(0, 0.);
181 fgWKParticipants->SetParameter(1, fSigmaNN);
182 fgWKParticipants->SetParameter(2, fA);
183 fgWKParticipants->SetNpx(200);
184 fgWKParticipants->SetNpy(20);
187 // Overlap and Participants
190 fgWStaa = new TF1("WStaa", WStaa, 0., fgBMax, 1);
191 fgWStaa->SetNpx(100);
192 fgWStaa->SetParameter(0,fA);
193 fgWStaa->SetNpx(100);
194 fgWParticipants = new TF1("WParticipants", WParticipants, 0., fgBMax, 2);
195 fgWParticipants->SetParameter(0, fSigmaNN);
196 fgWParticipants->SetParameter(1, fA);
197 fgWParticipants->SetNpx(100);
199 Info("Init","Reading overlap function from file %s",fName.Data());
200 TFile* f = new TFile(fName.Data());
202 Fatal("Init", "Could not open file %s",fName.Data());
204 fgWStaa = (TF1*) f->Get("WStaa");
205 fgWParticipants = (TF1*) f->Get("WParticipants");
212 fgWEnergyDensity = new TF1("WEnergyDensity", WEnergyDensity, 0., 2. * fWSr0, 1);
213 fgWEnergyDensity->SetParameter(0, fWSr0 + 1.);
216 // Geometrical Cross-Section
218 fgWSgeo = new TF1("WSgeo", WSgeo, 0., fgBMax, 1);
219 fgWSgeo->SetParameter(0,fSigmaNN); //mbarn
220 fgWSgeo->SetNpx(100);
223 // Hard cross section (binary collisions)
225 fgWSbinary = new TF1("WSbinary", WSbinary, 0., fgBMax, 1);
226 fgWSbinary->SetParameter(0, fSigmaHard); //mbarn
227 fgWSbinary->SetNpx(100);
230 // Hard collisions per event
232 fgWSN = new TF1("WSN", WSN, 0., fgBMax, 1);
236 // Almond shaped interaction region
238 fgWAlmond = new TF2("WAlmond", WAlmond, -fgBMax, fgBMax, -fgBMax, fgBMax, 1);
239 fgWAlmond->SetParameter(0, 0.);
240 fgWAlmond->SetNpx(200);
241 fgWAlmond->SetNpy(200);
244 Info("Init","Reading interaction almonds from file: %s",fName.Data());
245 Char_t almondName[100];
246 TFile* ff = new TFile(fName.Data());
247 for(Int_t k=0; k<40; k++) {
248 sprintf(almondName,"WAlmondFixedB%d",k);
249 fgWAlmondCurrent = (TF2*)ff->Get(almondName);
250 new(&fgWAlmondFixedB[k]) TF2(*fgWAlmondCurrent);
255 fgWIntRadius = new TF1("WIntRadius", WIntRadius, 0., fgBMax, 1);
256 fgWIntRadius->SetParameter(0, 0.);
259 // Path Length as a function of Phi
261 fgWPathLength0 = new TF1("WPathLength0", WPathLength0, -TMath::Pi(), TMath::Pi(), 2);
262 fgWPathLength0->SetParameter(0, 0.);
263 fgWPathLength0->SetParameter(1, 0.); //Pathlength definition
265 fgWPathLength = new TF1("WPathLength", WPathLength, -TMath::Pi(), TMath::Pi(), 3);
266 fgWPathLength->SetParameter(0, 0.); //Impact Parameter
267 fgWPathLength->SetParameter(1, 1000.); //Number of interactions used for average
268 fgWPathLength->SetParameter(2, 0); //Pathlength definition
271 void AliFastGlauber::Reset()
274 // Reset dynamic allocated formulas
275 // in case init is called twice
277 if(fgWSb) delete fgWSb;
278 if(fgWSbz) delete fgWSbz;
279 if(fgWSz) delete fgWSz;
280 if(fgWSta) delete fgWSta;
281 if(fgWStarfi) delete fgWStarfi;
282 if(fgWAlmond) delete fgWAlmond;
283 if(fgWStaa) delete fgWStaa;
284 if(fgWSgeo) delete fgWSgeo;
285 if(fgWSbinary) delete fgWSbinary;
286 if(fgWSN) delete fgWSN;
287 if(fgWPathLength0) delete fgWPathLength0;
288 if(fgWPathLength) delete fgWPathLength;
289 if(fgWEnergyDensity) delete fgWEnergyDensity;
290 if(fgWIntRadius) delete fgWIntRadius;
291 if(fgWKParticipants) delete fgWKParticipants;
292 if(fgWParticipants) delete fgWParticipants;
295 void AliFastGlauber::DrawWSb() const
298 // Draw Wood-Saxon Nuclear Density Function
300 TCanvas *c1 = new TCanvas("c1","Wood Saxon",400,10,600,700);
302 Double_t max=fgWSb->GetMaximum(0,fgBMax)*1.01;
303 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->GetXaxis()->SetTitle("r [fm]");
306 h2f->GetYaxis()->SetNoExponent(kTRUE);
307 h2f->GetYaxis()->SetTitle("#rho [fm^{-3}]");
310 TLegend *l1a = new TLegend(0.45,0.6,.90,0.8);
311 l1a->SetFillStyle(0);
312 l1a->SetBorderSize(0);
314 sprintf(label,"r_{0} = %.2f fm",fWSr0);
315 l1a->AddEntry(fgWSb,label,"");
316 sprintf(label,"d = %.2f fm",fWSd);
317 l1a->AddEntry(fgWSb,label,"");
318 sprintf(label,"n = %.2e fm^{-3}",fWSn);
319 l1a->AddEntry(fgWSb,label,"");
320 sprintf(label,"#omega = %.2f",fWSw);
321 l1a->AddEntry(fgWSb,label,"");
326 void AliFastGlauber::DrawOverlap() const
329 // Draw Overlap Function
331 TCanvas *c2 = new TCanvas("c2","Overlap",400,10,600,700);
333 Double_t max=fgWStaa->GetMaximum(0,fgBMax)*1.01;
334 TH2F *h2f=new TH2F("h2ftaa","Overlap function: T_{AB} [mbarn^{-1}]",2,0,fgBMax,2,0,max);
336 h2f->GetXaxis()->SetTitle("b [fm]");
337 h2f->GetYaxis()->SetTitle("T_{AB} [mbarn^{-1}]");
339 fgWStaa->Draw("same");
342 void AliFastGlauber::DrawParticipants() const
345 // Draw Number of Participants Npart
347 TCanvas *c3 = new TCanvas("c3","Participants",400,10,600,700);
349 Double_t max=fgWParticipants->GetMaximum(0,fgBMax)*1.01;
350 TH2F *h2f=new TH2F("h2fpart","Number of Participants",2,0,fgBMax,2,0,max);
352 h2f->GetXaxis()->SetTitle("b [fm]");
353 h2f->GetYaxis()->SetTitle("N_{part}");
355 fgWParticipants->Draw("same");
356 TLegend *l1a = new TLegend(0.50,0.75,.90,0.9);
357 l1a->SetFillStyle(0);
358 l1a->SetBorderSize(0);
360 sprintf(label,"#sigma^{inel.}_{NN} = %.1f mbarn",fSigmaNN);
361 l1a->AddEntry(fgWParticipants,label,"");
366 void AliFastGlauber::DrawThickness() const
369 // Draw Thickness Function
371 TCanvas *c4 = new TCanvas("c4","Thickness",400,10,600,700);
373 Double_t max=fgWSta->GetMaximum(0,fgBMax)*1.01;
374 TH2F *h2f=new TH2F("h2fta","Thickness function: T_{A} [fm^{-2}]",2,0,fgBMax,2,0,max);
376 h2f->GetXaxis()->SetTitle("b [fm]");
377 h2f->GetYaxis()->SetTitle("T_{A} [fm^{-2}]");
379 fgWSta->Draw("same");
382 void AliFastGlauber::DrawGeo() const
385 // Draw Geometrical Cross-Section
387 TCanvas *c5 = new TCanvas("c5","Geometrical Cross-Section",400,10,600,700);
389 Double_t max=fgWSgeo->GetMaximum(0,fgBMax)*1.01;
390 TH2F *h2f=new TH2F("h2fgeo","Differential Geometrical Cross-Section: d#sigma^{geo}_{AB}/db [fm]",2,0,fgBMax,2,0,max);
392 h2f->GetXaxis()->SetTitle("b [fm]");
393 h2f->GetYaxis()->SetTitle("d#sigma^{geo}_{AB}/db [fm]");
395 fgWSgeo->Draw("same");
396 TLegend *l1a = new TLegend(0.10,0.8,.40,0.9);
397 l1a->SetFillStyle(0);
398 l1a->SetBorderSize(0);
400 sprintf(label,"#sigma_{NN}^{inel.} = %.1f mbarn",fSigmaNN);
401 l1a->AddEntry(fgWSgeo,label,"");
406 void AliFastGlauber::DrawBinary() const
409 // Draw Binary Cross-Section
411 TCanvas *c6 = new TCanvas("c6","Binary Cross-Section",400,10,600,700);
413 Double_t max=fgWSbinary->GetMaximum(0,fgBMax)*1.01;
414 TH2F *h2f=new TH2F("h2fbinary","Differential Binary Cross-Section: #sigma^{hard}_{NN} dT_{AB}/db [fm]",2,0,fgBMax,2,0,max);
416 h2f->GetXaxis()->SetTitle("b [fm]");
417 h2f->GetYaxis()->SetTitle("d#sigma^{hard}_{AB}/db [fm]");
419 fgWSbinary->Draw("same");
420 TLegend *l1a = new TLegend(0.50,0.8,.90,0.9);
421 l1a->SetFillStyle(0);
422 l1a->SetBorderSize(0);
424 sprintf(label,"#sigma_{NN}^{hard} = %.1f mbarn",fSigmaHard);
425 l1a->AddEntry(fgWSb,label,"");
430 void AliFastGlauber::DrawN() const
433 // Draw Binaries per event (Ncoll)
435 TCanvas *c7 = new TCanvas("c7","Binaries per event",400,10,600,700);
437 Double_t max=fgWSN->GetMaximum(0,fgBMax)*1.01;
438 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->GetXaxis()->SetTitle("b [fm]");
441 h2f->GetYaxis()->SetTitle("N_{coll}");
444 TLegend *l1a = new TLegend(0.50,0.75,.90,0.9);
445 l1a->SetFillStyle(0);
446 l1a->SetBorderSize(0);
448 sprintf(label,"#sigma^{hard}_{NN} = %.1f mbarn",fSigmaHard);
449 l1a->AddEntry(fgWSN,label,"");
450 sprintf(label,"#sigma^{inel.}_{NN} = %.1f mbarn",fSigmaNN);
451 l1a->AddEntry(fgWSN,label,"");
456 void AliFastGlauber::DrawKernel(Double_t b) const
461 TCanvas *c8 = new TCanvas("c8","Kernel",400,10,600,700);
463 fgWStarfi->SetParameter(0, b);
464 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->GetXaxis()->SetTitle("r [fm]");
467 h2f->GetYaxis()->SetTitle("#phi [rad]");
469 fgWStarfi->Draw("same");
470 TLegend *l1a = new TLegend(0.65,0.8,.90,0.9);
471 l1a->SetFillStyle(0);
472 l1a->SetBorderSize(0);
474 sprintf(label,"b = %.1f fm",b);
475 l1a->AddEntry(fgWStarfi,label,"");
480 void AliFastGlauber::DrawAlmond(Double_t b) const
483 // Draw Interaction Almond
485 TCanvas *c9 = new TCanvas("c9","Almond",400,10,600,700);
487 fgWAlmond->SetParameter(0, b);
488 TH2F *h2f=new TH2F("h2falmond","Interaction Almond [fm^{-4}]",2,0,fgBMax,2,0,fgBMax);
490 h2f->GetXaxis()->SetTitle("x [fm]");
491 h2f->GetYaxis()->SetTitle("y [fm]");
493 fgWAlmond->Draw("same");
494 TLegend *l1a = new TLegend(0.65,0.8,.90,0.9);
495 l1a->SetFillStyle(0);
496 l1a->SetBorderSize(0);
498 sprintf(label,"b = %.1f fm",b);
499 l1a->AddEntry(fgWAlmond,label,"");
504 void AliFastGlauber::DrawEnergyDensity() const
507 // Draw energy density
509 TCanvas *c10 = new TCanvas("c10","Energy Density",400, 10, 600, 700);
511 fgWEnergyDensity->SetMinimum(0.);
512 Double_t max=fgWEnergyDensity->GetMaximum(0,fgWEnergyDensity->GetParameter(0))*1.01;
513 TH2F *h2f=new TH2F("h2fenergydens","Energy density",2,0,fgBMax,2,0,max);
515 h2f->GetXaxis()->SetTitle("b [fm]");
516 h2f->GetYaxis()->SetTitle("fm^{-4}");
518 fgWEnergyDensity->Draw("same");
522 void AliFastGlauber::DrawPathLength0(Double_t b, Int_t iopt) const
527 TCanvas *c11 = new TCanvas("c11","Path Length",400,10,600,700);
529 fgWPathLength0->SetParameter(0, b);
530 fgWPathLength0->SetParameter(1, Double_t(iopt));
531 fgWPathLength0->SetMinimum(0.);
532 fgWPathLength0->SetMaximum(10.);
533 TH2F *h2f=new TH2F("h2fpathlength0","Path length",2,-TMath::Pi(), TMath::Pi(),2,0,10.);
535 h2f->GetXaxis()->SetTitle("#phi [rad]");
536 h2f->GetYaxis()->SetTitle("l [fm]");
538 fgWPathLength0->Draw("same");
541 void AliFastGlauber::DrawPathLength(Double_t b , Int_t ni, Int_t iopt) const
546 TCanvas *c12 = new TCanvas("c12","Path Length",400,10,600,700);
548 fgWAlmond->SetParameter(0, b);
549 fgWPathLength->SetParameter(0, b);
550 fgWPathLength->SetParameter(1, Double_t (ni));
551 fgWPathLength->SetParameter(2, Double_t (iopt));
552 fgWPathLength->SetMinimum(0.);
553 fgWPathLength->SetMaximum(10.);
554 TH2F *h2f=new TH2F("h2fpathlength","Path length",2,-TMath::Pi(), TMath::Pi(),2,0,10.);
556 h2f->GetXaxis()->SetTitle("#phi [rad]");
557 h2f->GetYaxis()->SetTitle("l [fm]");
559 fgWPathLength->Draw("same");
562 void AliFastGlauber::DrawIntRadius(Double_t b) const
565 // Draw Interaction Radius
567 TCanvas *c13 = new TCanvas("c13","Interaction Radius",400,10,600,700);
569 fgWIntRadius->SetParameter(0, b);
570 fgWIntRadius->SetMinimum(0);
571 Double_t max=fgWIntRadius->GetMaximum(0,fgBMax)*1.01;
572 TH2F *h2f=new TH2F("h2fintradius","Interaction Density",2,0.,fgBMax,2,0,max);
574 h2f->GetXaxis()->SetTitle("r [fm]");
575 h2f->GetYaxis()->SetTitle("[fm^{-3}]");
577 fgWIntRadius->Draw("same");
580 Double_t AliFastGlauber::WSb(Double_t* x, Double_t* par)
583 // Woods-Saxon Parameterisation
584 // as a function of radius (xx)
586 const Double_t kxx = x[0]; //fm
587 const Double_t kr0 = par[0]; //fm
588 const Double_t kd = par[1]; //fm
589 const Double_t kw = par[2]; //no units
590 const Double_t kn = par[3]; //fm^-3 (used to normalize integral to one)
591 Double_t y = kn * (1.+kw*(kxx/kr0)*(kxx/kr0))/(1.+TMath::Exp((kxx-kr0)/kd));
595 Double_t AliFastGlauber::WSbz(Double_t* x, Double_t* par)
598 // Wood Saxon Parameterisation
599 // as a function of z and b
601 const Double_t kbb = x[0]; //fm
602 const Double_t kzz = x[1]; //fm
603 const Double_t kr0 = par[0]; //fm
604 const Double_t kd = par[1]; //fm
605 const Double_t kw = par[2]; //no units
606 const Double_t kn = par[3]; //fm^-3 (used to normalize integral to one)
607 const Double_t kxx = TMath::Sqrt(kbb*kbb+kzz*kzz);
608 Double_t y = kn * (1.+kw*(kxx/kr0)*(kxx/kr0))/(1.+TMath::Exp((kxx-kr0)/kd));
612 Double_t AliFastGlauber::WSz(Double_t* x, Double_t* par)
615 // Wood Saxon Parameterisation
616 // as a function of z for fixed b
618 const Double_t kzz = x[0]; //fm
619 const Double_t kr0 = par[0]; //fm
620 const Double_t kd = par[1]; //fm
621 const Double_t kw = par[2]; //no units
622 const Double_t kn = par[3]; //fm^-3 (used to normalize integral to one)
623 const Double_t kbb = par[4]; //fm
624 const Double_t kxx = TMath::Sqrt(kbb*kbb+kzz*kzz);
625 Double_t y = kn * (1.+kw*(kxx/kr0)*(kxx/kr0))/(1.+TMath::Exp((kxx-kr0)/kd));
629 Double_t AliFastGlauber::WSta(Double_t* x, Double_t* /*par*/)
632 // Thickness function T_A
633 // as a function of b
635 const Double_t kb = x[0];
636 fgWSz->SetParameter(4, kb);
637 Double_t y = 2. * fgWSz->Integral(0., fgBMax);
641 Double_t AliFastGlauber::WStarfi(Double_t* x, Double_t* par)
644 // Kernel for overlap function: T_A(s)*T_A(s-b)
645 // as a function of r and phi
646 const Double_t kr1 = x[0];
647 const Double_t kphi = x[1];
648 const Double_t kb = par[0];
649 const Double_t kr2 = TMath::Sqrt(kr1*kr1 + kb*kb - 2.*kr1*kb*TMath::Cos(kphi));
650 Double_t y = kr1 * fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
654 Double_t AliFastGlauber::WStaa(Double_t* x, Double_t* par)
658 // T_{AB}=Int d2s T_A(s)*T_B(s-b)
659 // as a function of b
660 // (normalized to fA*fB)
662 const Double_t kb = x[0];
663 const Double_t ka = par[0];
664 fgWStarfi->SetParameter(0, kb);
666 // root integration seems to fail
676 Double_t y = 2. * 208. * 208. * fgWStarfi->IntegralMultiple(2, al, bl, 0.001, err);
677 printf("WStaa: %.5e %.5e %.5e\n", b, y, err);
684 for (Int_t i = 0; i < fgkMCInts; i++)
686 const Double_t kphi = TMath::Pi() * gRandom->Rndm();
687 const Double_t kb1 = fgBMax * gRandom->Rndm();
688 y += fgWStarfi->Eval(kb1, kphi);
690 y *= 2. * TMath::Pi() * fgBMax / fgkMCInts; //fm^-2
691 y *= ka * ka * 0.1; //mbarn^-1
695 Double_t AliFastGlauber::WKParticipants(Double_t* x, Double_t* par)
698 // Kernel for number of participants
699 // as a function of r and phi
701 const Double_t kr1 = x[0];
702 const Double_t kphi = x[1];
703 const Double_t kb = par[0]; //fm
704 const Double_t ksig = par[1]; //mbarn
705 const Double_t ka = par[2]; //mass number
706 const Double_t kr2 = TMath::Sqrt(kr1*kr1 +kb*kb - 2.*kr1*kb*TMath::Cos(kphi));
707 const Double_t kxsi = fgWSta->Eval(kr2) * ksig * 0.1; //no units
709 Double_t y=(1-TMath::Power((1-xsi),aa))
712 Double_t sum = ka * kxsi;
714 for (Int_t i = 1; i <= ka; i++)
717 sum *= (-kxsi) * a / Float_t(i+1);
720 y *= kr1 * fgWSta->Eval(kr1);
724 Double_t AliFastGlauber::WParticipants(Double_t* x, Double_t* par)
727 // Number of Participants as
730 const Double_t kb = x[0];
731 const Double_t ksig = par[0]; //mbarn
732 const Double_t ka = par[1]; //mass number
733 fgWKParticipants->SetParameter(0, kb);
734 fgWKParticipants->SetParameter(1, ksig);
735 fgWKParticipants->SetParameter(2, ka);
741 for (Int_t i = 0; i < fgkMCInts; i++)
743 const Double_t kphi = TMath::Pi() * gRandom->Rndm();
744 const Double_t kb1 = fgBMax * gRandom->Rndm();
745 y += fgWKParticipants->Eval(kb1, kphi);
747 y *= 2. * ka * 2. * TMath::Pi() * fgBMax / fgkMCInts;
751 Double_t AliFastGlauber::WSgeo(Double_t* x, Double_t* par)
754 // Geometrical Cross-Section
755 // as a function of b
757 const Double_t kb = x[0]; //fm
758 const Double_t ksigNN = par[0]; //mbarn
759 const Double_t ktaa = fgWStaa->Eval(kb); //mbarn^-1
760 Double_t y = 2. * TMath::Pi() * kb * (1. - TMath::Exp(- ksigNN * ktaa));
764 Double_t AliFastGlauber::WSbinary(Double_t* x, Double_t* par)
767 // Number of binary hard collisions
768 // as a function of b
770 const Double_t kb = x[0]; //fm
771 const Double_t ksig = par[0]; //mbarn
772 const Double_t ktaa = fgWStaa->Eval(kb); //mbarn^-1
773 Double_t y = 2. * TMath::Pi() * kb * ksig * ktaa;
777 Double_t AliFastGlauber::WSN(Double_t* x, Double_t* /*par*/)
780 // Number of hard processes per event
781 // as a function of b
782 const Double_t kb = x[0];
783 Double_t y = fgWSbinary->Eval(kb)/fgWSgeo->Eval(kb);
787 Double_t AliFastGlauber::WEnergyDensity(Double_t* x, Double_t* par)
790 // Initial energy density
791 // as a function of the impact parameter
793 const Double_t kb = x[0];
794 const Double_t krA = par[0];
796 // Attention: area of transverse reaction zone in hard-sphere approximation !
797 const Double_t krA2=krA*krA;
798 const Double_t kb2=kb*kb;
799 const Double_t ksaa = (TMath::Pi() - 2. * TMath::ASin(kb/ 2./ krA)) * krA2
800 - kb * TMath::Sqrt(krA2 - kb2/ 4.); //fm^2
801 const Double_t ktaa = fgWStaa->Eval(kb); //mbarn^-1
802 Double_t y=ktaa/ksaa*10;
806 Double_t AliFastGlauber::WAlmond(Double_t* x, Double_t* par)
809 // Almond shaped interaction region
810 // as a function of cartesian x,y.
812 const Double_t kb = par[0];
813 const Double_t kxx = x[0] + kb/2.;
814 const Double_t kyy = x[1];
815 const Double_t kr1 = TMath::Sqrt(kxx*kxx + kyy*kyy);
816 const Double_t kphi = TMath::ATan2(kyy,kxx);
817 const Double_t kr2 = TMath::Sqrt(kr1*kr1 + kb*kb - 2.*kr1*kb*TMath::Cos(kphi));
819 // Interaction probability calculated as product of thicknesses
821 Double_t y = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
825 Double_t AliFastGlauber::WIntRadius(Double_t* x, Double_t* par)
828 // Average interaction density over radius
829 // at which interaction takes place
830 // as a function of radius
832 const Double_t kr = x[0];
833 const Double_t kb = par[0];
834 fgWAlmond->SetParameter(0, kb);
835 // Average over phi in small steps
836 const Double_t kdphi = 2. * TMath::Pi() / 100.;
839 for (Int_t i = 0; i < 100; i++) {
840 const Double_t kxx = kr * TMath::Cos(phi);
841 const Double_t kyy = kr * TMath::Sin(phi);
842 y += fgWAlmond->Eval(kxx,kyy);
845 // Result multiplied by Jacobian (2 pi r)
846 y *= 2. * TMath::Pi() * kr / 100.;
850 Double_t AliFastGlauber::WPathLength0(Double_t* x, Double_t* par)
853 // Path Length as a function of phi
854 // for interaction point fixed at (0,0)
855 // as a function of phi-direction
857 // Phi direction in Almond
858 const Double_t kphi0 = x[0];
859 const Double_t kb = par[0];
860 // Path Length definition
861 const Int_t kiopt = Int_t(par[1]);
863 // Step along radial direction phi
864 const Int_t kNp = 100; // Steps in r
865 const Double_t kDr = fgBMax/kNp;
869 for (Int_t i = 0; i < kNp; i++) {
871 // Transform into target frame
873 const Double_t kxx = r * TMath::Cos(kphi0) + kb / 2.;
874 const Double_t kyy = r * TMath::Sin(kphi0);
875 const Double_t kphi = TMath::ATan2(kyy, kxx);
876 const Double_t kr1 = TMath::Sqrt(kxx*kxx + kyy*kyy);
877 // Radius in projectile frame
878 const Double_t kr2 = TMath::Sqrt(kr1*kr1 + kb*kb - 2.*kr1*kb*TMath::Cos(kphi));
879 const Double_t ky = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
887 if (!kiopt) // My length definition (is exact for hard disk)
888 if(w) y= 2. * rw / w;
890 const Double_t knorm=fgWSta->Eval(1e-4);
891 if(knorm) y = TMath::Sqrt(2. * rw * kDr / knorm / knorm);
896 Double_t AliFastGlauber::WPathLength(Double_t* x, Double_t* par)
899 // Path Length as a function of phi
900 // Interaction point from random distribution
901 // as a function of the phi-direction
902 const Double_t kphi0 = x[0];
903 const Double_t kb = par[0];
904 fgWAlmond->SetParameter(0, kb);
905 const Int_t kNpi = Int_t (par[1]); //Number of interactions
906 const Int_t kiopt = Int_t(par[2]); //Path Length definition
911 const Int_t kNp = 100;
912 const Double_t kDr = fgBMax/Double_t(kNp);
913 Double_t l = 0.; // Path length
914 for (Int_t in = 0; in < kNpi; in ++) {
919 fgWAlmond->GetRandom2(x0, y0);
921 const Double_t kr0 = TMath::Sqrt(x0*x0 + y0*y0);
922 const Int_t knps = Int_t ((fgBMax - kr0)/kDr) - 1;
926 for (Int_t i = 0; (i < knps ); i++) {
927 // Transform into target frame
928 const Double_t kxx = x0 + r * TMath::Cos(kphi0) + kb / 2.;
929 const Double_t kyy = y0 + r * TMath::Sin(kphi0);
930 const Double_t kphi = TMath::ATan2(kyy, kxx);
931 const Double_t kr1 = TMath::Sqrt(kxx*kxx + kyy*kyy);
932 // Radius in projectile frame
933 const Double_t kr2 = TMath::Sqrt(kr1*kr1 + kb*kb - 2.*kr1*kb*TMath::Cos(kphi));
934 const Double_t ky = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
940 // Average over interactions
942 if(w) l += (2. * rw / w);
944 const Double_t knorm=fgWSta->Eval(1e-4);
945 if(knorm) l+= 2. * rw * kDr / knorm / knorm;
952 ret=TMath::Sqrt( l / kNpi);
956 Double_t AliFastGlauber::CrossSection(Double_t b1, Double_t b2) const
959 // Return the geometrical cross-section integrated from b1 to b2
961 return fgWSgeo->Integral(b1, b2)*10.; //mbarn
964 Double_t AliFastGlauber::HardCrossSection(Double_t b1, Double_t b2) const
967 // Return the hard cross-section integrated from b1 to b2
969 return fgWSbinary->Integral(b1, b2)*10.; //mbarn
972 Double_t AliFastGlauber::FractionOfHardCrossSection(Double_t b1, Double_t b2) const
975 // Return raction of hard cross-section integrated from b1 to b2
977 return fgWSbinary->Integral(b1, b2)/fgWSbinary->Integral(0., 100.);
980 Double_t AliFastGlauber::Binaries(Double_t b) const
983 // Return number of binary hard collisions normalized to 1 at b=0
986 return fgWSN->Eval(b)/fgWSN->Eval(1e-4);
989 Double_t AliFastGlauber::GetNumberOfBinaries(Double_t b) const
992 // Return number of binary hard collisions at b
995 return fgWSN->Eval(b);
998 Double_t AliFastGlauber::Participants(Double_t b) const
1001 // Return the number of participants normalized to 1 at b=0
1004 return (fgWParticipants->Eval(b)/fgWParticipants->Eval(1e-4));
1007 Double_t AliFastGlauber::GetNumberOfParticipants(Double_t b) const
1010 // Return the number of participants for impact parameter b
1013 return (fgWParticipants->Eval(b));
1016 Double_t AliFastGlauber::GetNumberOfCollisions(Double_t b) const
1019 // Return the number of collisions for impact parameter b
1022 return (fgWStaa->Eval(b)*fSigmaNN);
1025 void AliFastGlauber::SimulateTrigger(Int_t n)
1028 // Simulates Trigger
1030 TH1F* mbtH = new TH1F("mbtH", "MB Trigger b-Distribution", 100, 0., 20.);
1031 TH1F* hdtH = new TH1F("hdtH", "Hard Trigger b-Distribution", 100, 0., 20.);
1032 TH1F* mbmH = new TH1F("mbmH", "MB Trigger Multiplicity Distribution", 100, 0., 8000.);
1033 TH1F* hdmH = new TH1F("hdmH", "Hard Trigger Multiplicity Distribution", 100, 0., 8000.);
1035 mbtH->SetXTitle("b [fm]");
1036 hdtH->SetXTitle("b [fm]");
1037 mbmH->SetXTitle("Multiplicity");
1038 hdmH->SetXTitle("Multiplicity");
1040 TCanvas *c0 = new TCanvas("c0","Trigger Simulation",400,10,600,700);
1042 TCanvas *c1 = new TCanvas("c1","Trigger Simulation",400,10,600,700);
1048 for (Int_t iev = 0; iev < n; iev++)
1051 GetRandom(b, p, mult);
1054 mbmH->Fill(mult, 1.);
1055 hdmH->Fill(mult, p);
1071 void AliFastGlauber::GetRandom(Float_t& b, Float_t& p, Float_t& mult)
1074 // Gives back a random impact parameter, hard trigger probability and multiplicity
1076 b = fgWSgeo->GetRandom();
1077 const Float_t kmu = fgWSN->Eval(b);
1078 p = 1.-TMath::Exp(-kmu);
1079 mult = 6000./fgWSN->Eval(1.) * kmu;
1082 void AliFastGlauber::GetRandom(Int_t& bin, Bool_t& hard)
1085 // Gives back a random impact parameter bin, and hard trigger decission
1087 const Float_t kb = fgWSgeo->GetRandom();
1088 const Float_t kmu = fgWSN->Eval(kb) * fSigmaHard;
1089 const Float_t kp = 1.-TMath::Exp(-kmu);
1092 } else if (kb < 8.6) {
1094 } else if (kb < 11.2) {
1096 } else if (kb < 13.2) {
1098 } else if (kb < 15.0) {
1104 const Float_t kr = gRandom->Rndm();
1105 if (kr < kp) hard = kTRUE;
1108 Double_t AliFastGlauber::GetRandomImpactParameter(Double_t bmin, Double_t bmax)
1111 // Gives back a random impact parameter in the range bmin .. bmax
1114 while(b < bmin || b > bmax)
1115 b = fgWSgeo->GetRandom();
1119 void AliFastGlauber::StoreFunctions() const
1122 // Store in file functions
1124 TFile* ff = new TFile(fName.Data(),"recreate");
1125 fgWStaa->Write("WStaa");
1126 fgWParticipants->Write("WParticipants");
1131 //=================== Added by A. Dainese 11/02/04 ===========================
1133 void AliFastGlauber::StoreAlmonds() const
1137 // 40 almonds for b = (0.25+k*0.5) fm (k=0->39)
1139 Char_t almondName[100];
1140 TFile* ff = new TFile(fName.Data(),"update");
1141 for(Int_t k=0; k<40; k++) {
1142 sprintf(almondName,"WAlmondFixedB%d",k);
1143 Double_t b = 0.25+k*0.5;
1144 Info("StoreAlmonds"," b = %f\n",b);
1145 fgWAlmond->SetParameter(0,b);
1146 fgWAlmond->Write(almondName);
1152 void AliFastGlauber::SetCentralityClass(Double_t xsecFrLow,Double_t xsecFrUp)
1155 // Set limits of centrality class as fractions
1156 // of the geomtrical cross section
1158 if(xsecFrLow>1. || xsecFrUp>1. || xsecFrLow>xsecFrUp) {
1159 Error("SetCentralityClass", "Please set 0 <= xsecFrLow <= xsecFrUp <= 1\n");
1163 Double_t bLow=0.,bUp=0.;
1165 const Double_t knorm=fgWSgeo->Integral(0.,100.);
1166 while(xsecFr<xsecFrLow) {
1167 xsecFr = fgWSgeo->Integral(0.,bLow)/knorm;
1171 while(xsecFr<xsecFrUp) {
1172 xsecFr = fgWSgeo->Integral(0.,bUp)/knorm;
1176 Info("SetCentralityClass", "Centrality class: %4.2f-%4.2f; %4.1f < b < %4.1f fm",
1177 xsecFrLow,xsecFrUp,bLow,bUp);
1178 fgWSbinary->SetRange(bLow,bUp);
1184 void AliFastGlauber::GetRandomBHard(Double_t& b)
1187 // Get random impact parameter according to distribution of
1188 // hard (binary) cross-section, in the range defined by the centrality class
1190 b = fgWSbinary->GetRandom();
1191 Int_t bin = 2*(Int_t)b;
1192 if( (b-(Int_t)b) > 0.5) bin++;
1193 fgWAlmondCurrent = &fgWAlmondFixedB[bin];
1197 void AliFastGlauber::GetRandomXY(Double_t& x,Double_t& y)
1200 // Get random position of parton production point according to
1201 // product of thickness functions
1203 fgWAlmondCurrent->GetRandom2(x,y);
1207 void AliFastGlauber::GetRandomPhi(Double_t& phi)
1210 // Get random parton azimuthal propagation direction
1212 phi = 2.*TMath::Pi()*gRandom->Rndm();
1216 Double_t AliFastGlauber::CalculateLength(Double_t b,Double_t x0,Double_t y0,Double_t phi0)
1219 // Calculate path length for a parton with production point (x0,y0)
1220 // and propagation direction (ux=cos(phi0),uy=sin(phi0))
1221 // in a collision with impact parameter b
1224 // number of steps in l
1225 const Int_t kNp = 100;
1226 const Double_t kDl = fgBMax/Double_t(kNp);
1232 // ell = 2 * \int_0^\infty dl*l*(T_A*T_B)(x0+l*ux,y0+l*uy) /
1233 // \int_0^\infty dl*(T_A*T_B)(x0+l*ux,y0+l*uy)
1237 const Double_t kr0 = TMath::Sqrt(x0*x0 + y0*y0);
1238 const Int_t knps = Int_t ((fgBMax - kr0)/kDl) - 1;
1240 Double_t integral1 = 0.;
1241 Double_t integral2 = 0.;
1243 for (Int_t i = 0; i < knps; i++) {
1245 // Transform into target frame
1246 const Double_t kxx = x0 + l * TMath::Cos(phi0) + b / 2.;
1247 const Double_t kyy = y0 + l * TMath::Sin(phi0);
1248 const Double_t kphi = TMath::ATan2(kyy, kxx);
1249 const Double_t kr1 = TMath::Sqrt(kxx*kxx + kyy*kyy);
1250 // Radius in projectile frame
1251 const Double_t kr2 = TMath::Sqrt(kr1*kr1 + b*b - 2.*kr1*b*TMath::Cos(kphi));
1252 const Double_t kprodTATB = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
1254 integral1 += kprodTATB * l * kDl;
1255 integral2 += kprodTATB * kDl;
1261 ell = (2. * integral1 / integral2);
1263 } else if(fEllDef==2) {
1267 // ell = \int_0^\infty dl*
1268 // \Theta((T_A*T_B)(x0+l*ux,y0+l*uy)-0.5*(T_A*T_B)(0,0))
1272 const Double_t kr0 = TMath::Sqrt(x0*x0 + y0*y0);
1273 const Int_t knps = Int_t ((fgBMax - kr0)/kDl) - 1;
1274 const Double_t kprodTATBHalfMax = 0.5*fgWAlmondCurrent->Eval(0.,0.);
1277 Double_t integral = 0.;
1278 for (Int_t i = 0; i < knps; i++) {
1279 // Transform into target frame
1280 const Double_t kxx = x0 + l * TMath::Cos(phi0) + b / 2.;
1281 const Double_t kyy = y0 + l * TMath::Sin(phi0);
1282 const Double_t kphi = TMath::ATan2(kyy, kxx);
1283 const Double_t kr1 = TMath::Sqrt(kxx*kxx + kyy*kyy);
1284 // Radius in projectile frame
1285 const Double_t kr2 = TMath::Sqrt(kr1*kr1 + b*b - 2.*kr1*b*TMath::Cos(kphi));
1286 const Double_t kprodTATB = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
1287 if(kprodTATB>kprodTATBHalfMax) integral += kDl;
1290 Double_t ell = integral;
1293 Error("CalculateLength","Wrong length definition setting: %d !\n",fEllDef);
1298 void AliFastGlauber::GetLengthAndPhi(Double_t& ell,Double_t& phi,Double_t b)
1301 // Return length from random b, x0, y0, phi0
1304 Double_t x0,y0,phi0;
1305 if(b<0.) GetRandomBHard(b);
1309 ell = CalculateLength(b,x0,y0,phi0);
1313 void AliFastGlauber::GetLength(Double_t& ell,Double_t b)
1316 // Return length from random b, x0, y0, phi0
1319 GetLengthAndPhi(ell,phi,b);
1323 void AliFastGlauber::GetLengthsBackToBackAndPhi(Double_t& ell1,Double_t& ell2,Double_t &phi,Double_t b)
1326 // Return 2 lengths back to back from random b, x0, y0, phi0
1329 Double_t x0,y0,phi0;
1330 if(b<0.) GetRandomBHard(b);
1333 const Double_t kphi0plusPi = phi0+TMath::Pi();
1335 ell1 = CalculateLength(b,x0,y0,phi0);
1336 ell2 = CalculateLength(b,x0,y0,kphi0plusPi);
1340 void AliFastGlauber::GetLengthsBackToBack(Double_t& ell1,Double_t& ell2,
1344 // Return 2 lengths back to back from random b, x0, y0, phi0
1347 GetLengthsBackToBackAndPhi(ell1,ell2,phi,b);
1351 void AliFastGlauber::GetLengthsForPythia(Int_t n,Double_t* phi,Double_t* ell, Double_t b)
1354 // Returns lenghts for n partons with azimuthal angles phi[n]
1355 // from random b, x0, y0
1358 if(b < 0.) GetRandomBHard(b);
1360 for(Int_t i = 0; i< n; i++) ell[i] = CalculateLength(b,x0,y0,phi[i]);
1364 void AliFastGlauber::PlotBDistr(Int_t n)
1367 // Plot distribution of n impact parameters
1370 TH1F *hB = new TH1F("hB","dN/db",100,0,fgBMax);
1371 hB->SetXTitle("b [fm]");
1372 hB->SetYTitle("dN/db [a.u.]");
1373 hB->SetFillColor(3);
1374 for(Int_t i=0; i<n; i++) {
1378 TCanvas *cB = new TCanvas("cB","Impact parameter distribution",0,0,500,500);
1384 void AliFastGlauber::PlotLengthDistr(Int_t n,Bool_t save,Char_t *fname)
1387 // Plot length distribution
1390 TH1F *hEll = new TH1F("hEll","Length distribution",64,-0.5,15);
1391 hEll->SetXTitle("Transverse path length, L [fm]");
1392 hEll->SetYTitle("Probability");
1393 hEll->SetFillColor(2);
1394 for(Int_t i=0; i<n; i++) {
1398 hEll->Scale(1/(Double_t)n);
1399 TCanvas *cL = new TCanvas("cL","Length distribution",0,0,500,500);
1404 TFile *f = new TFile(fname,"recreate");
1411 void AliFastGlauber::PlotLengthB2BDistr(Int_t n,Bool_t save,Char_t *fname)
1414 // Plot lengths back-to-back distributions
1417 TH2F *hElls = new TH2F("hElls","Lengths back-to-back",100,0,15,100,0,15);
1418 hElls->SetXTitle("Transverse path length, L [fm]");
1419 hElls->SetYTitle("Transverse path length, L [fm]");
1420 for(Int_t i=0; i<n; i++) {
1421 GetLengthsBackToBack(ell1,ell2);
1422 hElls->Fill(ell1,ell2);
1424 hElls->Scale(1/(Double_t)n);
1425 TCanvas *cLs = new TCanvas("cLs","Length back-to-back distribution",0,0,500,500);
1426 gStyle->SetPalette(1,0);
1428 hElls->Draw("col,Z");
1430 TFile *f = new TFile(fname,"recreate");
1437 void AliFastGlauber::PlotAlmonds() const
1440 // Plot almonds for some impact parameters
1442 TCanvas *c = new TCanvas("c","Almonds",0,0,500,500);
1443 gStyle->SetPalette(1,0);
1446 fgWAlmondFixedB[0].Draw("cont1");
1448 fgWAlmondFixedB[10].Draw("cont1");
1450 fgWAlmondFixedB[20].Draw("cont1");
1452 fgWAlmondFixedB[30].Draw("cont1");
1456 //=================== Added by A. Dainese 05/03/04 ===========================
1458 void AliFastGlauber::CalculateI0I1(Double_t& integral0,Double_t& integral1,
1459 Double_t b,Double_t x0,Double_t y0,
1460 Double_t phi0,Double_t ellCut) const
1463 // Calculate integrals:
1464 // integral0 = \int_0^ellCut dl*(T_A*T_B)(x0+l*ux,y0+l*uy)
1465 // integral1 = \int_0^ellCut dl*l*(T_A*T_B)(x0+l*ux,y0+l*uy)
1467 // for a parton with production point (x0,y0)
1468 // and propagation direction (ux=cos(phi0),uy=sin(phi0))
1469 // in a collision with impact parameter b
1472 // number of steps in l
1473 const Int_t kNp = 100;
1474 const Double_t kDl = fgBMax/Double_t(kNp);
1477 const Double_t kr0 = TMath::Sqrt(x0 * x0 + y0 * y0);
1478 const Int_t knps = Int_t ((fgBMax - kr0)/kDl) - 1;
1485 while((i < knps) && (l < ellCut)) {
1486 // Transform into target frame
1487 const Double_t kxx = x0 + l * TMath::Cos(phi0) + b / 2.;
1488 const Double_t kyy = y0 + l * TMath::Sin(phi0);
1489 const Double_t kphi = TMath::ATan2(kyy, kxx);
1490 const Double_t kr1 = TMath::Sqrt(kxx*kxx + kyy*kyy);
1491 // Radius in projectile frame
1492 const Double_t kr2 = TMath::Sqrt(kr1*kr1 + b*b - 2.*kr1*b*TMath::Cos(kphi));
1493 const Double_t kprodTATB = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
1494 integral0 += kprodTATB * kDl;
1495 integral1 += kprodTATB * l * kDl;
1502 void AliFastGlauber::GetI0I1AndPhi(Double_t& integral0,Double_t& integral1,
1504 Double_t ellCut,Double_t b)
1507 // Return I0 and I1 from random b, x0, y0, phi0
1510 Double_t x0,y0,phi0;
1511 if(b<0.) GetRandomBHard(b);
1515 CalculateI0I1(integral0,integral1,b,x0,y0,phi0,ellCut);
1519 void AliFastGlauber::GetI0I1(Double_t& integral0,Double_t& integral1,
1520 Double_t ellCut,Double_t b)
1523 // Return I0 and I1 from random b, x0, y0, phi0
1526 GetI0I1AndPhi(integral0,integral1,phi,ellCut,b);
1530 void AliFastGlauber::GetI0I1BackToBackAndPhi(Double_t& integral01,Double_t& integral11,
1531 Double_t& integral02,Double_t& integral12,
1533 Double_t ellCut,Double_t b)
1536 // Return 2 pairs of I0 and I1 back to back from random b, x0, y0, phi0
1539 Double_t x0,y0,phi0;
1540 if(b<0.) GetRandomBHard(b);
1544 const Double_t kphi0plusPi = phi0+TMath::Pi();
1545 CalculateI0I1(integral01,integral11,b,x0,y0,phi0,ellCut);
1546 CalculateI0I1(integral02,integral12,b,x0,y0,kphi0plusPi,ellCut);
1550 void AliFastGlauber::GetI0I1BackToBack(Double_t& integral01,Double_t& integral11,
1551 Double_t& integral02,Double_t& integral12,
1552 Double_t ellCut,Double_t b)
1555 // Return 2 pairs of I0 and I1 back to back from random b, x0, y0, phi0
1558 GetI0I1BackToBackAndPhi(integral01,integral11,integral02,integral12,
1563 void AliFastGlauber::GetI0I1ForPythia(Int_t n,Double_t* phi,
1564 Double_t* integral0,Double_t* integral1,
1565 Double_t ellCut,Double_t b)
1568 // Returns I0 and I1 pairs for n partons with azimuthal angles phi[n]
1569 // from random b, x0, y0
1572 if(b<0.) GetRandomBHard(b);
1574 for(Int_t i=0; i<n; i++)
1575 CalculateI0I1(integral0[i],integral1[i],b,x0,y0,phi[i],ellCut);
1579 void AliFastGlauber::PlotI0I1Distr(Int_t n,Double_t ellCut,
1580 Bool_t save,Char_t *fname)
1583 // Plot I0-I1 distribution
1586 TH2F *hI0I1s = new TH2F("hI0I1s","I_{0} versus I_{1}",1000,0,0.001,1000,0,0.01);
1587 hI0I1s->SetXTitle("I_{0} [fm^{-3}]");
1588 hI0I1s->SetYTitle("I_{1} [fm^{-2}]");
1590 TH1F *hI0 = new TH1F("hI0","I_{0} = #hat{q}L / k",
1592 hI0->SetXTitle("I_{0} [fm^{-3}]");
1593 hI0->SetYTitle("Probability");
1594 hI0->SetFillColor(3);
1595 TH1F *hI1 = new TH1F("hI1","I_{1} = #omega_{c} / k",
1597 hI1->SetXTitle("I_{1} [fm^{-2}]");
1598 hI1->SetYTitle("Probability");
1599 hI1->SetFillColor(4);
1600 TH1F *h2 = new TH1F("h2","2 I_{1}^{2}/I_{0} = R / k",
1602 h2->SetXTitle("2 I_{1}^{2}/I_{0} [fm^{-1}]");
1603 h2->SetYTitle("Probability");
1604 h2->SetFillColor(2);
1605 TH1F *h3 = new TH1F("h3","2 I_{1}/I_{0} = L",
1607 h3->SetXTitle("2 I_{1}/I_{0} [fm]");
1608 h3->SetYTitle("Probability");
1609 h3->SetFillColor(5);
1610 TH1F *h4 = new TH1F("h4","I_{0}^{2}/(2 I_{1}) = #hat{q} / k",
1612 h4->SetXTitle("I_{0}^{2}/(2 I_{1}) [fm^{-4}]");
1613 h4->SetYTitle("Probability");
1614 h4->SetFillColor(7);
1616 for(Int_t i=0; i<n; i++) {
1617 GetI0I1(i0,i1,ellCut);
1618 hI0I1s->Fill(i0,i1);
1621 h2->Fill(2.*i1*i1/i0);
1623 h4->Fill(i0*i0/2./i1);
1625 hI0->Scale(1/(Double_t)n);
1626 hI1->Scale(1/(Double_t)n);
1627 h2->Scale(1/(Double_t)n);
1628 h3->Scale(1/(Double_t)n);
1629 h4->Scale(1/(Double_t)n);
1630 hI0I1s->Scale(1/(Double_t)n);
1632 TCanvas *cI0I1 = new TCanvas("cI0I1","I0 and I1",0,0,900,700);
1645 gStyle->SetPalette(1,0);
1646 hI0I1s->Draw("col,Z");
1649 TFile *f = new TFile(fname,"recreate");
1661 void AliFastGlauber::PlotI0I1B2BDistr(Int_t n,Double_t ellCut,
1662 Bool_t save,Char_t *fname)
1665 // Plot I0-I1 back-to-back distributions
1667 Double_t i01,i11,i02,i12;
1668 TH2F *hI0s = new TH2F("hI0s","I_{0}'s back-to-back",100,0,100,100,0,100);
1669 hI0s->SetXTitle("I_{0} [fm^{-3}]");
1670 hI0s->SetYTitle("I_{0} [fm^{-3}]");
1671 TH2F *hI1s = new TH2F("hI1s","I_{1}'s back-to-back",100,0,100,100,0,100);
1672 hI1s->SetXTitle("I_{1} [fm^{-2}]");
1673 hI1s->SetYTitle("I_{1} [fm^{-2}]");
1675 for(Int_t i=0; i<n; i++) {
1676 GetI0I1BackToBack(i01,i11,i02,i12,ellCut);
1677 hI0s->Fill(i01,i02);
1678 hI1s->Fill(i11,i12);
1680 hI0s->Scale(1/(Double_t)n);
1681 hI1s->Scale(1/(Double_t)n);
1683 TCanvas *cI0I1s = new TCanvas("cI0I1s","I0 and I1 back-to-back distributions",0,0,800,400);
1684 gStyle->SetPalette(1,0);
1685 cI0I1s->Divide(2,1);
1687 hI0s->Draw("col,Z");
1689 hI1s->Draw("col,Z");
1692 TFile *f = new TFile(fname,"recreate");