FASTSIM/AliFastGlauber.cxx

   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
  16 /* $Id$ */
  17
  18 // from AliRoot
  19 #include "AliFastGlauber.h"
  20 // from root
  21 #include <TH1F.h>
  22 #include <TF1.h>
  23 #include <TF2.h>
  24 #include <TCanvas.h>
  25 #include <TRandom.h>
  26 #include <TFile.h>
  27
  28 ClassImp(AliFastGlauber)
  29
  30 TF1*    AliFastGlauber::fWSb            = NULL;
  31 TF2*    AliFastGlauber::fWSbz           = NULL;
  32 TF1*    AliFastGlauber::fWSz            = NULL;
  33 TF1*    AliFastGlauber::fWSta           = NULL;
  34 TF2*    AliFastGlauber::fWStarfi        = NULL;
  35 TF2*    AliFastGlauber::fWAlmond        = NULL;
  36 TF1*    AliFastGlauber::fWStaa          = NULL;
  37 TF1*    AliFastGlauber::fWSgeo          = NULL;
  38 TF1*    AliFastGlauber::fWSbinary       = NULL;
  39 TF1*    AliFastGlauber::fWSN            = NULL;
  40 TF1*    AliFastGlauber::fWPathLength0   = NULL;
  41 TF1*    AliFastGlauber::fWPathLength    = NULL;
  42 TF1*    AliFastGlauber::fWEnergyDensity = NULL;
  43 TF1*    AliFastGlauber::fWIntRadius     = NULL;
  44 Float_t AliFastGlauber::fbMax           = 0.;
  45
  46 AliFastGlauber::AliFastGlauber()
  47 {
  48 //  Default Constructor
  49 //
  50 //  Defaults for Pb
  51     SetWoodSaxonParameters(6.624, 0.549, 0.00, 7.69e-4);
  52     SetHardCrossSection();
  53     SetMaxImpact();
  54 }
  55
  56 void AliFastGlauber::Init(Int_t mode)
  57 {
  58 // Initialisation
  59 //
  60 //  Wood-Saxon
  61 //
  62     fWSb = new TF1("WSb", WSb, 0, fbMax, 4);
  63     fWSb->SetParameter(0, fWSr0);
  64     fWSb->SetParameter(1, fWSd);
  65     fWSb->SetParameter(2, fWSw);
  66     fWSb->SetParameter(3, fWSn);
  67
  68     fWSbz = new TF2("WSbz", WSbz, 0, fbMax, 4);
  69     fWSbz->SetParameter(0, fWSr0);
  70     fWSbz->SetParameter(1, fWSd);
  71     fWSbz->SetParameter(2, fWSw);
  72     fWSbz->SetParameter(3, fWSn);
  73
  74     fWSz = new TF1("WSz", WSz, 0, fbMax, 5);
  75     fWSz->SetParameter(0, fWSr0);
  76     fWSz->SetParameter(1, fWSd);
  77     fWSz->SetParameter(2, fWSw);
  78     fWSz->SetParameter(3, fWSn);
  79
  80 //
  81 //  Thickness
  82 //
  83     fWSta = new TF1("WSta", WSta, 0., fbMax, 0);
  84
  85 //
  86 //  Overlap Kernel
  87 //
  88     fWStarfi = new TF2("WStarfi", WStarfi, 0., fbMax, 0., TMath::Pi(), 1);
  89     fWStarfi->SetParameter(0, 0.);
  90     fWStarfi->SetNpx(200);
  91     fWStarfi->SetNpy(20);
  92 //
  93 //  Almond shaped interaction region
  94 //
  95     fWAlmond = new TF2("WAlmond", WAlmond, -fbMax, fbMax, -fbMax, fbMax, 1);
  96     fWAlmond->SetParameter(0, 0.);
  97     fWAlmond->SetNpx(200);
  98     fWAlmond->SetNpy(200);
  99 //
 100 //  Path Length as a function of Phi
 101 //
 102     fWPathLength0 = new TF1("WPathLength0", WPathLength0, -TMath::Pi(), TMath::Pi(), 1);
 103     fWPathLength0->SetParameter(0, 0.);
 104     fWPathLength = new TF1("WPathLength", WPathLength, -TMath::Pi(), TMath::Pi(), 2);
 105 //  Impact Parameter
 106     fWPathLength->SetParameter(0, 0.);
 107 //  Number of interactions used for average
 108     fWPathLength->SetParameter(1, 1000.);
 109
 110     fWIntRadius = new TF1("WIntRadius", WIntRadius, 0., fbMax, 1);
 111     fWIntRadius->SetParameter(0, 0.);
 112
 113
 114 //
 115 //  Overlap
 116 //
 117     if (! mode) {
 118         fWStaa = new TF1("WStaa", WStaa, 0., fbMax, 0);
 119         fWStaa->SetNpx(100);
 120     } else {
 121         TFile* f = new TFile("$(ALICE_ROOT)/FASTSIM/data/glauberPbPb.root");
 122         fWStaa = (TF1*) f->Get("WStaa");
 123     }
 124
 125 //
 126     fWEnergyDensity = new TF1("WEnergyDensity", WEnergyDensity, 0., 2. * fWSr0, 1);
 127     fWEnergyDensity->SetParameter(0, fWSr0 + 1.);
 128
 129 //
 130 //  Geometrical Cross-Section
 131 //
 132     fWSgeo = new TF1("WSgeo", WSgeo, 0., fbMax, 0);
 133     fWSgeo->SetNpx(100);
 134 //
 135 //  Hard cross section (~ binary collisions)
 136 //
 137     fWSbinary = new TF1("WSbinary", WSbinary, 0., fbMax, 1);
 138     fWSbinary->SetParameter(0, fSigmaHard); // mb
 139     fWSbinary->SetNpx(100);
 140 //
 141 // Hard collisions per event
 142 //
 143     fWSN = new TF1("WSN", WSN, 0., fbMax, 1);
 144     fWSN->SetNpx(100);
 145 }
 146
 147 void AliFastGlauber::DrawWSb()
 148 {
 149 //
 150 //  Draw Wood-Saxon Nuclear Density Function
 151 //
 152     TCanvas *c1 = new TCanvas("c1","Wood Saxon",400,10,600,700);
 153     c1->cd();
 154     fWSb->Draw();
 155 }
 156
 157 void AliFastGlauber::DrawOverlap()
 158 {
 159 //
 160 //  Draw Overlap Function
 161 //
 162     TCanvas *c2 = new TCanvas("c2","Overlap",400,10,600,700);
 163     c2->cd();
 164     fWStaa->Draw();
 165 }
 166
 167 void AliFastGlauber::DrawThickness()
 168 {
 169 //
 170 //  Draw Thickness Function
 171 //
 172     TCanvas *c3 = new TCanvas("c3","Thickness",400,10,600,700);
 173     c3->cd();
 174     fWSta->Draw();
 175 }
 176
 177 void AliFastGlauber::DrawGeo()
 178 {
 179 //
 180 //  Draw Geometrical Cross-Section
 181 //
 182     TCanvas *c3 = new TCanvas("c3","Geometrical Cross-Section",400,10,600,700);
 183     c3->cd();
 184     fWSgeo->Draw();
 185 }
 186
 187 void AliFastGlauber::DrawBinary()
 188 {
 189 //
 190 //  Draw Wood-Saxon Nuclear Density Function
 191 //
 192     TCanvas *c4 = new TCanvas("c4","Binary Cross-Section",400,10,600,700);
 193     c4->cd();
 194     fWSbinary->Draw();
 195 }
 196
 197 void AliFastGlauber::DrawN()
 198 {
 199 //
 200 //  Draw Binaries per event
 201 //
 202     TCanvas *c5 = new TCanvas("c5","Binaries per event",400,10,600,700);
 203     c5->cd();
 204     fWSN->Draw();
 205 }
 206
 207 void AliFastGlauber::DrawKernel(Double_t b)
 208 {
 209 //
 210 //  Draw Kernel
 211 //
 212     TCanvas *c6 = new TCanvas("c6","Kernel",400,10,600,700);
 213     c6->cd();
 214     fWStarfi->SetParameter(0, b);
 215     fWStarfi->Draw();
 216 }
 217
 218 void AliFastGlauber::DrawAlmond(Double_t b)
 219 {
 220 //
 221 //  Draw Interaction Almond
 222 //
 223     TCanvas *c7 = new TCanvas("c7","Almond",400,10,600,700);
 224     c7->cd();
 225     fWAlmond->SetParameter(0, b);
 226     fWAlmond->Draw();
 227 }
 228
 229 void AliFastGlauber::DrawPathLength0(Double_t b)
 230 {
 231 //
 232 //  Draw Path Length
 233 //
 234     TCanvas *c8 = new TCanvas("c8","Path Length",400,10,600,700);
 235     c8->cd();
 236     fWPathLength0->SetParameter(0, b);
 237     fWPathLength0->SetMinimum(0.);
 238     fWPathLength0->SetMaximum(10.);
 239     fWPathLength0->Draw();
 240 }
 241
 242 void AliFastGlauber::DrawPathLength(Double_t b , Int_t ni)
 243 {
 244 //
 245 //  Draw Path Length
 246 //
 247     TCanvas *c9 = new TCanvas("c9","Path Length",400,10,600,700);
 248     c9->cd();
 249     fWAlmond->SetParameter(0, b);
 250
 251     fWPathLength->SetParameter(0, b);
 252     fWPathLength->SetParameter(1, Double_t (ni));
 253     fWPathLength->SetMinimum(0.);
 254     fWPathLength->SetMaximum(10.);
 255     fWPathLength->Draw();
 256 }
 257
 258 void AliFastGlauber::DrawIntRadius(Double_t b)
 259 {
 260 //
 261 //  Draw Interaction Radius
 262 //
 263     TCanvas *c10 = new TCanvas("c10","Interaction Radius",400,10,600,700);
 264     c10->cd();
 265     fWIntRadius->SetParameter(0, b);
 266     fWIntRadius->SetMinimum(0.);
 267     fWIntRadius->Draw();
 268 }
 269
 270 void AliFastGlauber::DrawEnergyDensity()
 271 {
 272 //
 273 //  Draw energy density
 274 //
 275     TCanvas *c11 = new TCanvas("c11","Energy Density",400, 10, 600, 700);
 276     c11->cd();
 277     fWEnergyDensity->SetMinimum(0.);
 278     fWEnergyDensity->Draw();
 279 }
 280
 281 Double_t AliFastGlauber::WSb(Double_t* x, Double_t* par)
 282 {
 283 //
 284 //  Woods-Saxon Parameterisation
 285 //  as a function of radius
 286 //
 287     Double_t xx  = x[0];
 288     Double_t r0  = par[0];
 289     Double_t d   = par[1];
 290     Double_t w   = par[2];
 291     Double_t n   = par[3];
 292
 293     Double_t y  = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
 294
 295     return y;
 296 }
 297
 298 Double_t AliFastGlauber::WSbz(Double_t* x, Double_t* par)
 299 {
 300 //
 301 //  Wood Saxon Parameterisation
 302 //  as a function of z and  b
 303 //
 304     Double_t bb  = x[0];
 305     Double_t zz  = x[1];
 306     Double_t r0  = par[0];
 307     Double_t d   = par[1];
 308     Double_t w   = par[2];
 309     Double_t n   = par[3];
 310     Double_t xx  = TMath::Sqrt(bb*bb+zz*zz);
 311     Double_t y  = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
 312
 313     return y;
 314 }
 315
 316 Double_t AliFastGlauber::WSz(Double_t* x, Double_t* par)
 317 {
 318 //
 319 //  Wood Saxon Parameterisation
 320 //  as a function of z for fixed b
 321 //
 322     Double_t bb  = par[4];
 323     Double_t zz  = x[0];
 324     Double_t r0  = par[0];
 325     Double_t d   = par[1];
 326     Double_t w   = par[2];
 327     Double_t n   = par[3];
 328     Double_t xx  = TMath::Sqrt(bb*bb+zz*zz);
 329     Double_t y  = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
 330
 331     return y;
 332 }
 333
 334 Double_t AliFastGlauber::WSta(Double_t* x, Double_t* /*par*/)
 335 {
 336 //
 337 //  Thickness function
 338 //
 339     Double_t b  = x[0];
 340     fWSz->SetParameter(4, b);
 341     Double_t y  = 2. * fWSz->Integral(0., fbMax);
 342     return y;
 343 }
 344
 345
 346
 347 Double_t AliFastGlauber::WStarfi(Double_t* x, Double_t* par)
 348 {
 349 //
 350 //  Kernel for overlap function
 351 //
 352     Double_t b    = par[0];
 353     Double_t r1   = x[0];
 354     Double_t phi  = x[1];
 355     Double_t r2   = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
 356     Double_t y    = r1 * fWSta->Eval(r1) * fWSta->Eval(r2);
 357     return y;
 358 }
 359
 360
 361 Double_t AliFastGlauber::WAlmond(Double_t* x, Double_t* par)
 362 {
 363 //
 364 //  Almond shaped interaction region
 365 //
 366     Double_t b    = par[0];
 367     Double_t xx   = x[0] + b/2.;
 368     Double_t yy   = x[1];
 369     Double_t r1   = TMath::Sqrt(xx * xx + yy * yy);
 370     Double_t phi  = TMath::ATan2(yy,xx);
 371
 372     Double_t r2   = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
 373 //
 374 //  Interaction probability calculated as product of thicknesses
 375 //
 376     Double_t y    = fWSta->Eval(r1) * fWSta->Eval(r2);
 377     return y;
 378 }
 379
 380 Double_t AliFastGlauber::WIntRadius(Double_t* x, Double_t* par)
 381 {
 382 //
 383 //  Average radius at which interaction takes place
 384 //
 385 //  Radius in the Almond
 386     Double_t r    = x[0];
 387 //  Impact parameter
 388     Double_t b    = par[0];
 389     fWAlmond->SetParameter(0, b);
 390 //  Steps in phi
 391     Double_t dphi = 2. * TMath::Pi() / 100.;
 392 //  Average over phi
 393     Double_t phi  = 0.;
 394     Double_t y    = 0.;
 395
 396     for (Int_t i = 0; i < 100; i++) {
 397         Double_t xx = r * TMath::Cos(phi);
 398         Double_t yy = r * TMath::Sin(phi);
 399         y   += fWAlmond->Eval(xx,yy);
 400         phi += dphi;
 401     } // phi loop
 402 // Result multiplied by Jacobian (2 pi r)
 403     return (2. * TMath::Pi() * y * r / 100.);
 404 }
 405
 406 Double_t AliFastGlauber::WPathLength0(Double_t* x, Double_t* par)
 407 {
 408 //
 409 //  Path Length as a function of phi for interaction point fixed at (0,0)
 410 //
 411 //
 412 //  Steps in r
 413     const Int_t    np  = 100;
 414     const Double_t dr  = fbMax/Double_t(np);
 415 //  Impact parameter
 416     Double_t b      = par[0];
 417 //  Phi direction in Almond
 418     Double_t phi0   = x[0];
 419     Double_t r  = 0.;
 420     Double_t rw = 0.;
 421     Double_t w  = 0.;
 422 //  Step along radial direction phi
 423     for (Int_t i = 0; i < np; i++) {
 424 //
 425 //  Transform into target frame
 426 //
 427         Double_t xx   = r * TMath::Cos(phi0) + b / 2.;
 428         Double_t yy   = r * TMath::Sin(phi0);
 429         Double_t phi  = TMath::ATan2(yy, xx);
 430
 431         Double_t r1   = TMath::Sqrt(xx * xx + yy * yy);
 432 // Radius in projectile frame
 433         Double_t r2   = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
 434         Double_t y    = fWSta->Eval(r1) * fWSta->Eval(r2);
 435
 436         rw += y * r;
 437         w  += y;
 438         r  += dr;
 439     } // radial steps
 440 //
 441 //  My length definition (is exact for hard sphere)
 442     return (2. * rw / w);
 443 }
 444
 445 Double_t AliFastGlauber::WPathLength(Double_t* x, Double_t* par)
 446 {
 447 //
 448 //  Path Length as a function of phi
 449 //  Interaction point from random distribution
 450 //
 451 //
 452 //  r-steps
 453 //
 454     const Int_t    np   = 100;
 455     const Double_t dr  = fbMax/Double_t(np);
 456 //  Number of interactions
 457     const Int_t    npi  = Int_t (par[1]);
 458
 459 //
 460 //  Impact parameter
 461     Double_t b      = par[0];
 462 //  Phi direction
 463     Double_t phi0   = x[0];
 464
 465     printf("phi0 %f \n", phi0);
 466
 467 //  Path length
 468     Double_t l = 0.;
 469
 470     for (Int_t in = 0; in < npi; in ++) {
 471         Double_t rw = 0.;
 472         Double_t w  = 0.;
 473
 474         // Interaction point
 475         Double_t x0, y0;
 476         fWAlmond->GetRandom2(x0, y0);
 477 // Initial radius
 478         Double_t r0  = TMath::Sqrt(x0 * x0 + y0 * y0);
 479         Int_t    nps = Int_t ((fbMax - r0)/dr) - 1;
 480
 481         Double_t r  = 0.;
 482 // Radial steps
 483         for (Int_t i = 0; (i < nps ); i++) {
 484
 485 // Transform into target frame
 486             Double_t xx   = x0 + r * TMath::Cos(phi0) + b / 2.;
 487             Double_t yy   = y0 + r * TMath::Sin(phi0);
 488             Double_t phi  = TMath::ATan2(yy, xx);
 489             Double_t r1   = TMath::Sqrt(xx * xx + yy * yy);
 490 // Radius in projectile frame
 491             Double_t r2   = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
 492             Double_t y    = fWSta->Eval(r1) * fWSta->Eval(r2);
 493
 494             rw += y * r;
 495             w  += y;
 496             r  += dr;
 497         } // steps
 498 // Average over interactions
 499         l += 2. * rw / w;
 500     } // interactions
 501     return (l / Double_t(npi));
 502 }
 503
 504 Double_t AliFastGlauber::WStaa(Double_t* x, Double_t* /*par*/)
 505 {
 506 //
 507 //  Overlap function
 508 //
 509     Double_t b    = x[0];
 510     fWStarfi->SetParameter(0, b);
 511 /*
 512     Double_t al[2];
 513     Double_t bl[2];
 514     al[0] = 0.;
 515     al[1] = 0.;
 516     bl[0] = 6.6;
 517     bl[1] = TMath::Pi();
 518     Double_t err;
 519
 520     Double_t y =  2. * fWStarfi->IntegralMultiple(2, al, bl, 0.001, err);
 521     printf("WStaa: %f %f %f\n", b, y, err);
 522 */
 523 //
 524 //  MC Integration
 525 //
 526     Double_t y = 0;
 527     for (Int_t i = 0; i < 100000; i++)
 528     {
 529         Double_t phi = TMath::Pi() * gRandom->Rndm();
 530         Double_t b1  = fbMax       * gRandom->Rndm();
 531         y += fWStarfi->Eval(b1, phi);
 532     }
 533     y *= 2. * 0.1 *  208. * 208. * TMath::Pi() * fbMax / 100000.;
 534     return y;
 535 }
 536
 537 Double_t AliFastGlauber::WSgeo(Double_t* x, Double_t* /*par*/)
 538 {
 539 //
 540 //  Geometrical Cross-Section
 541 //
 542     Double_t b    = x[0];
 543     Double_t taa  = fWStaa->Eval(b);
 544     const Double_t sigma = 55.6; // mbarn
 545
 546     Double_t y    = 2. * TMath::Pi() * b * (1. - TMath::Exp(- sigma * taa)); // fm
 547     return y;
 548 }
 549
 550
 551 Double_t AliFastGlauber::WSbinary(Double_t* x, Double_t* par)
 552 {
 553 //
 554 //  Number of binary collisions
 555 //
 556     Double_t b     = x[0];
 557     Double_t sigma = par[0];
 558     Double_t taa   = fWStaa->Eval(b);
 559
 560     Double_t y    = 2. * TMath::Pi() * b * sigma * taa; // fm
 561     return y;
 562 }
 563
 564 Double_t AliFastGlauber::WSN(Double_t* x, Double_t* /*par*/)
 565 {
 566 //
 567 //  Number of hard processes per event
 568 //
 569     Double_t b     = x[0];
 570     Double_t y     = fWSbinary->Eval(b)/fWSgeo->Eval(b);
 571     return y;
 572 }
 573
 574 Double_t AliFastGlauber::WEnergyDensity(Double_t* x, Double_t* par)
 575 {
 576 //
 577 //  Initial energy density as a function of the impact parameter
 578 //
 579     Double_t b     = x[0];
 580     Double_t rA    = par[0];
 581
 582     Double_t saa   = (TMath::Pi() - 2. * TMath::ASin(b/ 2./ rA)) * rA * rA - b * TMath::Sqrt(rA * rA - b * b/ 4.);
 583     Double_t taa   = fWStaa->Eval(b);
 584
 585     return (taa/saa);
 586 }
 587
 588 void AliFastGlauber::SimulateTrigger(Int_t n)
 589 {
 590     //
 591     //  Simulates Trigger
 592     //
 593     TH1F* mbtH = new TH1F("mbtH", "MB Trigger b-Distribution",   100, 0., 20.);
 594     TH1F* hdtH = new TH1F("hdtH", "Hard Trigger b-Distribution", 100, 0., 20.);
 595     TH1F* mbmH = new TH1F("mbmH", "MB Trigger Multiplicity Distribution",   100, 0., 8000.);
 596     TH1F* hdmH = new TH1F("hdmH", "Hard Trigger Multiplicity Distribution", 100, 0., 8000.);
 597
 598     mbtH->SetXTitle("b [fm]");
 599     hdtH->SetXTitle("b [fm]");
 600     mbmH->SetXTitle("Multiplicity");
 601     hdmH->SetXTitle("Multiplicity");
 602
 603     TCanvas *c0 = new TCanvas("c0","Trigger Simulation",400,10,600,700);
 604     c0->Divide(2,1);
 605     TCanvas *c1 = new TCanvas("c1","Trigger Simulation",400,10,600,700);
 606     c1->Divide(1,2);
 607
 608     //
 609     //
 610     Init(1);
 611     for (Int_t iev = 0; iev < n; iev++)
 612     {
 613         Float_t b, p, mult;
 614         GetRandom(b, p, mult);
 615         mbtH->Fill(b,1.);
 616         hdtH->Fill(b, p);
 617         mbmH->Fill(mult, 1.);
 618         hdmH->Fill(mult, p);
 619
 620         c0->cd(1);
 621         mbtH->Draw();
 622         c0->cd(2);
 623         hdtH->Draw();
 624         c0->Update();
 625
 626         c1->cd(1);
 627         mbmH->Draw();
 628         c1->cd(2);
 629         hdmH->Draw();
 630         c1->Update();
 631     }
 632 }
 633
 634 void AliFastGlauber::GetRandom(Float_t& b, Float_t& p, Float_t& mult)
 635 {
 636     //
 637     // Gives back a random impact parameter, hard trigger probability and multiplicity
 638     //
 639         b = fWSgeo->GetRandom();
 640         Float_t mu = fWSN->Eval(b);
 641         p = 1.-TMath::Exp(-mu);
 642         mult = 6000./fWSN->Eval(1.) * mu;
 643 }
 644
 645 void AliFastGlauber::GetRandom(Int_t& bin, Bool_t& hard)
 646 {
 647     //
 648     // Gives back a random impact parameter bin, and hard trigger decission
 649     //
 650         Float_t b  = fWSgeo->GetRandom();
 651         Float_t mu = fWSN->Eval(b) * fSigmaHard;
 652         Float_t p  = 1.-TMath::Exp(-mu);
 653         if (b < 5.) {
 654             bin = 1;
 655         } else if (b <  8.6) {
 656             bin = 2;
 657         } else if (b < 11.2) {
 658             bin = 3;
 659         } else if (b < 13.2) {
 660             bin = 4;
 661         } else if (b < 15.0) {
 662             bin = 5;
 663         } else {
 664             bin = 6;
 665         }
 666
 667         hard = kFALSE;
 668
 669         Float_t r = gRandom->Rndm();
 670
 671         if (r < p) hard = kTRUE;
 672 }
 673
 674
 675 Float_t  AliFastGlauber::GetRandomImpactParameter(Float_t bmin, Float_t bmax)
 676 {
 677     //
 678     // Gives back a random impact parameter in the range bmin .. bmax
 679     //
 680
 681     Float_t b = -1.;
 682     while(b < bmin || b > bmax)
 683         b = fWSgeo->GetRandom();
 684     return b;
 685 }
 686
 687 Double_t AliFastGlauber::CrossSection(Double_t b1, Double_t b2)
 688 {
 689     //
 690     // Return cross-section integrated from b1 to b2
 691     //
 692
 693     return fWSgeo->Integral(b1, b2)/100.;
 694 }
 695
 696 Double_t AliFastGlauber::FractionOfHardCrossSection(Double_t b1, Double_t b2)
 697 {
 698     //
 699     // Return raction of hard cross-section integrated from b1 to b2
 700     //
 701
 702     return fWSbinary->Integral(b1, b2)/fWSbinary->Integral(0., 100.);
 703 }
 704
 705
 706 Double_t AliFastGlauber::Binaries(Double_t b)
 707 {
 708     //
 709     // Return number of binary collisions normalized to 1 at b=0
 710     //
 711
 712     return fWSN->Eval(b)/fWSN->Eval(0.);
 713 }