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