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