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