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