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