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