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 /*
  17 $Log$
  18 Revision 1.3  2003/04/21 09:35:53  morsch
  19 Protection against division by 0 in Binaries().
  20
  21 Revision 1.2  2003/04/14 14:23:44  morsch
  22 Correction in Binaries().
  23
  24 Revision 1.1  2003/04/10 08:48:13  morsch
  25 First commit.
  26
  27 */
  28
  29 // from AliRoot
  30 #include "AliFastGlauber.h"
  31 // from root
  32 #include <TH1F.h>
  33 #include <TF1.h>
  34 #include <TF2.h>
  35 #include <TCanvas.h>
  36 #include <TRandom.h>
  37 #include <TFile.h>
  38
  39 ClassImp(AliFastGlauber)
  40
  41 TF1*    AliFastGlauber::fWSb      = NULL;
  42 TF2*    AliFastGlauber::fWSbz     = NULL;
  43 TF1*    AliFastGlauber::fWSz      = NULL;
  44 TF1*    AliFastGlauber::fWSta     = NULL;
  45 TF2*    AliFastGlauber::fWStarfi  = NULL;
  46 TF1*    AliFastGlauber::fWStaa    = NULL;
  47 TF1*    AliFastGlauber::fWSgeo    = NULL;
  48 TF1*    AliFastGlauber::fWSbinary = NULL;
  49 TF1*    AliFastGlauber::fWSN      = NULL;
  50 Float_t AliFastGlauber::fbMax     = 0.;
  51
  52 AliFastGlauber::AliFastGlauber()
  53 {
  54 //  Default Constructor
  55 //
  56 //  Defaults for Pb
  57     SetWoodSaxonParameters(6.624, 0.549, 0.00, 7.69e-4);
  58     SetHardCrossSection();
  59     SetMaxImpact();
  60 }
  61
  62 void AliFastGlauber::Init(Int_t mode)
  63 {
  64 // Initialisation
  65 //
  66 //  Wood-Saxon
  67 //
  68     fWSb = new TF1("WSb", WSb, 0, fbMax, 4);
  69     fWSb->SetParameter(0, fWSr0);
  70     fWSb->SetParameter(1, fWSd);
  71     fWSb->SetParameter(2, fWSw);
  72     fWSb->SetParameter(3, fWSn);
  73
  74     fWSbz = new TF2("WSbz", WSbz, 0, fbMax, 4);
  75     fWSbz->SetParameter(0, fWSr0);
  76     fWSbz->SetParameter(1, fWSd);
  77     fWSbz->SetParameter(2, fWSw);
  78     fWSbz->SetParameter(3, fWSn);
  79
  80     fWSz = new TF1("WSz", WSz, 0, fbMax, 5);
  81     fWSz->SetParameter(0, fWSr0);
  82     fWSz->SetParameter(1, fWSd);
  83     fWSz->SetParameter(2, fWSw);
  84     fWSz->SetParameter(3, fWSn);
  85
  86 //
  87 //  Thickness
  88 //
  89     fWSta = new TF1("WSta", WSta, 0., fbMax, 0);
  90
  91 //
  92 //  Overlap Kernel
  93 //
  94     fWStarfi = new TF2("WStarfi", WStarfi, 0., fbMax, 0., TMath::Pi(), 1);
  95     fWStarfi->SetParameter(0, 0.);
  96     fWStarfi->SetNpx(200);
  97     fWStarfi->SetNpy(20);
  98 //
  99 //  Overlap
 100 //
 101     if (! mode) {
 102         fWStaa = new TF1("WStaa", WStaa, 0., fbMax, 0);
 103         fWStaa->SetNpx(100);
 104     } else {
 105         TFile* f = new TFile("$(ALICE_ROOT)/FASTSIM/data/glauberPbPb.root");
 106         fWStaa = (TF1*) f->Get("WStaa");
 107     }
 108
 109 //
 110 //  Geometrical Cross-Section
 111 //
 112     fWSgeo = new TF1("WSgeo", WSgeo, 0., fbMax, 0);
 113     fWSgeo->SetNpx(100);
 114 //
 115 //  Hard cross section (~ binary collisions)
 116 //
 117     fWSbinary = new TF1("WSbinary", WSbinary, 0., fbMax, 1);
 118     fWSbinary->SetParameter(0, fSigmaHard); // mb
 119     fWSbinary->SetNpx(100);
 120 //
 121 // Hard collisions per event
 122 //
 123     fWSN = new TF1("WSN", WSN, 0., fbMax, 1);
 124     fWSN->SetNpx(100);
 125 }
 126
 127 void AliFastGlauber::DrawWSb()
 128 {
 129 //
 130 //  Draw Wood-Saxon Nuclear Density Function
 131 //
 132     TCanvas *c1 = new TCanvas("c1","Wood Saxon",400,10,600,700);
 133     c1->cd();
 134     fWSb->Draw();
 135 }
 136
 137 void AliFastGlauber::DrawOverlap()
 138 {
 139 //
 140 //  Draw Overlap Function
 141 //
 142     TCanvas *c2 = new TCanvas("c2","Overlap",400,10,600,700);
 143     c2->cd();
 144     fWStaa->Draw();
 145 }
 146
 147 void AliFastGlauber::DrawThickness()
 148 {
 149 //
 150 //  Draw Thickness Function
 151 //
 152     TCanvas *c3 = new TCanvas("c3","Thickness",400,10,600,700);
 153     c3->cd();
 154     fWSta->Draw();
 155 }
 156
 157 void AliFastGlauber::DrawGeo()
 158 {
 159 //
 160 //  Draw Geometrical Cross-Section
 161 //
 162     TCanvas *c3 = new TCanvas("c3","Geometrical Cross-Section",400,10,600,700);
 163     c3->cd();
 164     fWSgeo->Draw();
 165 }
 166
 167 void AliFastGlauber::DrawBinary()
 168 {
 169 //
 170 //  Draw Wood-Saxon Nuclear Density Function
 171 //
 172     TCanvas *c4 = new TCanvas("c4","Binary Cross-Section",400,10,600,700);
 173     c4->cd();
 174     fWSbinary->Draw();
 175 }
 176
 177 void AliFastGlauber::DrawN()
 178 {
 179 //
 180 //  Draw Binaries per event
 181 //
 182     TCanvas *c5 = new TCanvas("c5","Binaries per event",400,10,600,700);
 183     c5->cd();
 184     fWSN->Draw();
 185 }
 186
 187 void AliFastGlauber::DrawKernel()
 188 {
 189 //
 190 //  Draw Kernel
 191 //
 192     TCanvas *c6 = new TCanvas("c6","Kernel",400,10,600,700);
 193     c6->cd();
 194     fWStarfi->Draw();
 195 }
 196
 197 Double_t AliFastGlauber::WSb(Double_t* x, Double_t* par)
 198 {
 199 //
 200 //  Woods-Saxon Parameterisation
 201 //  as a function of radius
 202 //
 203     Double_t xx  = x[0];
 204     Double_t r0  = par[0];
 205     Double_t d   = par[1];
 206     Double_t w   = par[2];
 207     Double_t n   = par[3];
 208
 209     Double_t y  = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
 210     return y;
 211 }
 212
 213 Double_t AliFastGlauber::WSbz(Double_t* x, Double_t* par)
 214 {
 215 //
 216 //  Wood Saxon Parameterisation
 217 //  as a function of z and  b
 218 //
 219     Double_t bb  = x[0];
 220     Double_t zz  = x[1];
 221     Double_t r0  = par[0];
 222     Double_t d   = par[1];
 223     Double_t w   = par[2];
 224     Double_t n   = par[3];
 225     Double_t xx  = TMath::Sqrt(bb*bb+zz*zz);
 226     Double_t y  = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
 227     return y;
 228 }
 229
 230 Double_t AliFastGlauber::WSz(Double_t* x, Double_t* par)
 231 {
 232 //
 233 //  Wood Saxon Parameterisation
 234 //  as a function of z for fixed b
 235 //
 236     Double_t bb  = par[4];
 237     Double_t zz  = x[0];
 238     Double_t r0  = par[0];
 239     Double_t d   = par[1];
 240     Double_t w   = par[2];
 241     Double_t n   = par[3];
 242     Double_t xx  = TMath::Sqrt(bb*bb+zz*zz);
 243     Double_t y  = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
 244 //    if (y < 1.e-6) y = 0;
 245     return y;
 246 }
 247
 248 Double_t AliFastGlauber::WSta(Double_t* x, Double_t* par)
 249 {
 250 //
 251 //  Thickness function
 252 //
 253     Double_t b  = x[0];
 254     fWSz->SetParameter(4, b);
 255     Double_t y  = 2. * fWSz->Integral(0., fbMax);
 256     return y;
 257 }
 258
 259
 260
 261 Double_t AliFastGlauber::WStarfi(Double_t* x, Double_t* par)
 262 {
 263 //
 264 //  Kernel for overlap function
 265 //
 266     Double_t b    = par[0];
 267     Double_t r1   = x[0];
 268     Double_t phi  = x[1];
 269     Double_t r2   = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
 270     Double_t y    = r1 * fWSta->Eval(r1) * fWSta->Eval(r2);
 271     return y;
 272 }
 273
 274
 275 Double_t AliFastGlauber::WStaa(Double_t* x, Double_t* par)
 276 {
 277 //
 278 //  Overlap function
 279 //
 280     Double_t b    = x[0];
 281     fWStarfi->SetParameter(0, b);
 282 /*
 283     Double_t al[2];
 284     Double_t bl[2];
 285     al[0] = 0.;
 286     al[1] = 0.;
 287     bl[0] = 6.6;
 288     bl[1] = TMath::Pi();
 289     Double_t err;
 290
 291     Double_t y =  2. * fWStarfi->IntegralMultiple(2, al, bl, 0.001, err);
 292     printf("WStaa: %f %f %f\n", b, y, err);
 293 */
 294 //
 295 //  MC Integration
 296 //
 297     Double_t y = 0;
 298     for (Int_t i = 0; i < 100000; i++)
 299     {
 300         Double_t phi = TMath::Pi() * gRandom->Rndm();
 301         Double_t b1  = fbMax       * gRandom->Rndm();
 302         y += fWStarfi->Eval(b1, phi);
 303     }
 304     y *= 2. * 0.1 *  208. * 208. * TMath::Pi() * fbMax / 100000.;
 305     return y;
 306 }
 307
 308 Double_t AliFastGlauber::WSgeo(Double_t* x, Double_t* par)
 309 {
 310 //
 311 //  Geometrical Cross-Section
 312 //
 313     Double_t b    = x[0];
 314     Double_t taa  = fWStaa->Eval(b);
 315     const Double_t sigma = 55.6; // mbarn
 316
 317     Double_t y    = 2. * TMath::Pi() * b * (1. - TMath::Exp(- sigma * taa)); // fm
 318     return y;
 319 }
 320
 321
 322 Double_t AliFastGlauber::WSbinary(Double_t* x, Double_t* par)
 323 {
 324 //
 325 //  Number of binary collisions
 326 //
 327     Double_t b     = x[0];
 328     Double_t sigma = par[0];
 329     Double_t taa   = fWStaa->Eval(b);
 330
 331     Double_t y    = 2. * TMath::Pi() * b * sigma * taa; // fm
 332     return y;
 333 }
 334
 335 Double_t AliFastGlauber::WSN(Double_t* x, Double_t* par)
 336 {
 337 //
 338 //  Number of hard processes per event
 339 //
 340     Double_t b     = x[0];
 341     Double_t y;
 342     if (b != 0.) {
 343         y = fWSbinary->Eval(b)/fWSgeo->Eval(b);
 344     } else {
 345         y = fWSbinary->Eval(1.e-4)/fWSgeo->Eval(1.e-4);
 346     }
 347     return y;
 348 }
 349
 350 void AliFastGlauber::SimulateTrigger(Int_t n)
 351 {
 352     //
 353     //  Simulates Trigger
 354     //
 355     TH1F* mbtH = new TH1F("mbtH", "MB Trigger b-Distribution",   100, 0., 20.);
 356     TH1F* hdtH = new TH1F("hdtH", "Hard Trigger b-Distribution", 100, 0., 20.);
 357     TH1F* mbmH = new TH1F("mbmH", "MB Trigger Multiplicity Distribution",   100, 0., 8000.);
 358     TH1F* hdmH = new TH1F("hdmH", "Hard Trigger Multiplicity Distribution", 100, 0., 8000.);
 359
 360     mbtH->SetXTitle("b [fm]");
 361     hdtH->SetXTitle("b [fm]");
 362     mbmH->SetXTitle("Multiplicity");
 363     hdmH->SetXTitle("Multiplicity");
 364
 365     TCanvas *c0 = new TCanvas("c0","Trigger Simulation",400,10,600,700);
 366     c0->Divide(2,1);
 367     TCanvas *c1 = new TCanvas("c1","Trigger Simulation",400,10,600,700);
 368     c1->Divide(1,2);
 369
 370     //
 371     //
 372     Init(1);
 373     for (Int_t iev = 0; iev < n; iev++)
 374     {
 375         Float_t b, p, mult;
 376         GetRandom(b, p, mult);
 377         mbtH->Fill(b,1.);
 378         hdtH->Fill(b, p);
 379         mbmH->Fill(mult, 1.);
 380         hdmH->Fill(mult, p);
 381
 382         c0->cd(1);
 383         mbtH->Draw();
 384         c0->cd(2);
 385         hdtH->Draw();
 386         c0->Update();
 387
 388         c1->cd(1);
 389         mbmH->Draw();
 390         c1->cd(2);
 391         hdmH->Draw();
 392         c1->Update();
 393     }
 394 }
 395
 396 void AliFastGlauber::GetRandom(Float_t& b, Float_t& p, Float_t& mult)
 397 {
 398     //
 399     // Gives back a random impact parameter, hard trigger probability and multiplicity
 400     //
 401         b = fWSgeo->GetRandom();
 402         Float_t mu = fWSN->Eval(b);
 403         p = 1.-TMath::Exp(-mu);
 404         mult = 6000./fWSN->Eval(0.) * mu;
 405 }
 406
 407 void AliFastGlauber::GetRandom(Int_t& bin, Bool_t& hard)
 408 {
 409     //
 410     // Gives back a random impact parameter bin, and hard trigger decission
 411     //
 412         Float_t b  = fWSgeo->GetRandom();
 413         Float_t mu = fWSN->Eval(b) * fSigmaHard;
 414         Float_t p  = 1.-TMath::Exp(-mu);
 415         if (b < 5.) {
 416             bin = 1;
 417         } else if (b <  8.6) {
 418             bin = 2;
 419         } else if (b < 11.2) {
 420             bin = 3;
 421         } else if (b < 13.2) {
 422             bin = 4;
 423         } else {
 424             bin = 5;
 425         }
 426
 427         hard = kFALSE;
 428
 429         Float_t r = gRandom->Rndm();
 430
 431         if (r < p) hard = kTRUE;
 432 }
 433
 434
 435 Float_t  AliFastGlauber::GetRandomImpactParameter(Float_t bmin, Float_t bmax)
 436 {
 437     //
 438     // Gives back a random impact parameter in the range bmin .. bmax
 439     //
 440
 441     Float_t b = -1.;
 442     while(b < bmin || b > bmax)
 443         b = fWSgeo->GetRandom();
 444     return b;
 445 }
 446
 447 Double_t AliFastGlauber::CrossSection(Double_t b1, Double_t b2)
 448 {
 449     //
 450     // Return cross-section integrated from b1 to b2
 451     //
 452
 453     return fWSgeo->Integral(b1, b2)/100.;
 454 }
 455
 456 Double_t AliFastGlauber::FractionOfHardCrossSection(Double_t b1, Double_t b2)
 457 {
 458     //
 459     // Return raction of hard cross-section integrated from b1 to b2
 460     //
 461
 462     return fWSbinary->Integral(b1, b2)/fWSbinary->Integral(0., 100.);
 463 }
 464
 465
 466 Double_t AliFastGlauber::Binaries(Double_t b)
 467 {
 468     //
 469     // Return number of binary collisions normalized to 1 at b=0
 470     //
 471
 472     return fWSN->Eval(b)/fWSN->Eval(0.);
 473 }