--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+/*
+$Log$
+*/
+
+// from AliRoot
+#include "AliFastGlauber.h"
+// from root
+#include <TH1F.h>
+#include <TF1.h>
+#include <TF2.h>
+#include <TCanvas.h>
+#include <TRandom.h>
+#include <TFile.h>
+
+ClassImp(AliFastGlauber)
+
+TF1* AliFastGlauber::fWSb = NULL;
+TF2* AliFastGlauber::fWSbz = NULL;
+TF1* AliFastGlauber::fWSz = NULL;
+TF1* AliFastGlauber::fWSta = NULL;
+TF2* AliFastGlauber::fWStarfi = NULL;
+TF1* AliFastGlauber::fWStaa = NULL;
+TF1* AliFastGlauber::fWSgeo = NULL;
+TF1* AliFastGlauber::fWSbinary = NULL;
+TF1* AliFastGlauber::fWSN = NULL;
+Float_t AliFastGlauber::fbMax = 0.;
+
+AliFastGlauber::AliFastGlauber()
+{
+// Default Constructor
+//
+// Defaults for Pb
+ SetWoodSaxonParameters(6.624, 0.549, 0.00, 7.69e-4);
+ SetHardCrossSection();
+ SetMaxImpact();
+}
+
+void AliFastGlauber::Init(Int_t mode)
+{
+// Initialisation
+//
+// Wood-Saxon
+//
+ fWSb = new TF1("WSb", WSb, 0, fbMax, 4);
+ fWSb->SetParameter(0, fWSr0);
+ fWSb->SetParameter(1, fWSd);
+ fWSb->SetParameter(2, fWSw);
+ fWSb->SetParameter(3, fWSn);
+
+ fWSbz = new TF2("WSbz", WSbz, 0, fbMax, 4);
+ fWSbz->SetParameter(0, fWSr0);
+ fWSbz->SetParameter(1, fWSd);
+ fWSbz->SetParameter(2, fWSw);
+ fWSbz->SetParameter(3, fWSn);
+
+ fWSz = new TF1("WSz", WSz, 0, fbMax, 5);
+ fWSz->SetParameter(0, fWSr0);
+ fWSz->SetParameter(1, fWSd);
+ fWSz->SetParameter(2, fWSw);
+ fWSz->SetParameter(3, fWSn);
+
+//
+// Thickness
+//
+ fWSta = new TF1("WSta", WSta, 0., fbMax, 0);
+
+//
+// Overlap Kernel
+//
+ fWStarfi = new TF2("WStarfi", WStarfi, 0., fbMax, 0., TMath::Pi(), 1);
+ fWStarfi->SetParameter(0, 0.);
+ fWStarfi->SetNpx(200);
+ fWStarfi->SetNpy(20);
+//
+// Overlap
+//
+ if (! mode) {
+ fWStaa = new TF1("WStaa", WStaa, 0., fbMax, 0);
+ fWStaa->SetNpx(100);
+ } else {
+ TFile* f = new TFile("$(ALICE_ROOT)/FASTSIM/data/glauberPbPb.root");
+ fWStaa = (TF1*) f->Get("WStaa");
+ }
+
+//
+// Geometrical Cross-Section
+//
+ fWSgeo = new TF1("WSgeo", WSgeo, 0., fbMax, 0);
+ fWSgeo->SetNpx(100);
+//
+// Hard cross section (~ binary collisions)
+//
+ fWSbinary = new TF1("WSbinary", WSbinary, 0., fbMax, 1);
+ fWSbinary->SetParameter(0, fSigmaHard); // mb
+ fWSbinary->SetNpx(100);
+//
+// Hard collisions per event
+//
+ fWSN = new TF1("WSN", WSN, 0., fbMax, 1);
+ fWSN->SetNpx(100);
+}
+
+void AliFastGlauber::DrawWSb()
+{
+//
+// Draw Wood-Saxon Nuclear Density Function
+//
+ TCanvas *c1 = new TCanvas("c1","Wood Saxon",400,10,600,700);
+ c1->cd();
+ fWSb->Draw();
+}
+
+void AliFastGlauber::DrawOverlap()
+{
+//
+// Draw Overlap Function
+//
+ TCanvas *c2 = new TCanvas("c2","Overlap",400,10,600,700);
+ c2->cd();
+ fWStaa->Draw();
+}
+
+void AliFastGlauber::DrawThickness()
+{
+//
+// Draw Thickness Function
+//
+ TCanvas *c3 = new TCanvas("c3","Thickness",400,10,600,700);
+ c3->cd();
+ fWSta->Draw();
+}
+
+void AliFastGlauber::DrawGeo()
+{
+//
+// Draw Geometrical Cross-Section
+//
+ TCanvas *c3 = new TCanvas("c3","Geometrical Cross-Section",400,10,600,700);
+ c3->cd();
+ fWSgeo->Draw();
+}
+
+void AliFastGlauber::DrawBinary()
+{
+//
+// Draw Wood-Saxon Nuclear Density Function
+//
+ TCanvas *c4 = new TCanvas("c4","Binary Cross-Section",400,10,600,700);
+ c4->cd();
+ fWSbinary->Draw();
+}
+
+void AliFastGlauber::DrawN()
+{
+//
+// Draw Binaries per event
+//
+ TCanvas *c5 = new TCanvas("c5","Binaries per event",400,10,600,700);
+ c5->cd();
+ fWSN->Draw();
+}
+
+void AliFastGlauber::DrawKernel()
+{
+//
+// Draw Kernel
+//
+ TCanvas *c6 = new TCanvas("c6","Kernel",400,10,600,700);
+ c6->cd();
+ fWStarfi->Draw();
+}
+
+Double_t AliFastGlauber::WSb(Double_t* x, Double_t* par)
+{
+//
+// Woods-Saxon Parameterisation
+// as a function of radius
+//
+ Double_t xx = x[0];
+ Double_t r0 = par[0];
+ Double_t d = par[1];
+ Double_t w = par[2];
+ Double_t n = par[3];
+
+ Double_t y = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
+ return y;
+}
+
+Double_t AliFastGlauber::WSbz(Double_t* x, Double_t* par)
+{
+//
+// Wood Saxon Parameterisation
+// as a function of z and b
+//
+ Double_t bb = x[0];
+ Double_t zz = x[1];
+ Double_t r0 = par[0];
+ Double_t d = par[1];
+ Double_t w = par[2];
+ Double_t n = par[3];
+ Double_t xx = TMath::Sqrt(bb*bb+zz*zz);
+ Double_t y = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
+ return y;
+}
+
+Double_t AliFastGlauber::WSz(Double_t* x, Double_t* par)
+{
+//
+// Wood Saxon Parameterisation
+// as a function of z for fixed b
+//
+ Double_t bb = par[4];
+ Double_t zz = x[0];
+ Double_t r0 = par[0];
+ Double_t d = par[1];
+ Double_t w = par[2];
+ Double_t n = par[3];
+ Double_t xx = TMath::Sqrt(bb*bb+zz*zz);
+ Double_t y = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
+// if (y < 1.e-6) y = 0;
+ return y;
+}
+
+Double_t AliFastGlauber::WSta(Double_t* x, Double_t* par)
+{
+//
+// Thickness function
+//
+ Double_t b = x[0];
+ fWSz->SetParameter(4, b);
+ Double_t y = 2. * fWSz->Integral(0., fbMax);
+ return y;
+}
+
+
+
+Double_t AliFastGlauber::WStarfi(Double_t* x, Double_t* par)
+{
+//
+// Kernel for overlap function
+//
+ Double_t b = par[0];
+ Double_t r1 = x[0];
+ Double_t phi = x[1];
+ Double_t r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
+ Double_t y = r1 * fWSta->Eval(r1) * fWSta->Eval(r2);
+ return y;
+}
+
+
+Double_t AliFastGlauber::WStaa(Double_t* x, Double_t* par)
+{
+//
+// Overlap function
+//
+ Double_t b = x[0];
+ fWStarfi->SetParameter(0, b);
+/*
+ Double_t al[2];
+ Double_t bl[2];
+ al[0] = 0.;
+ al[1] = 0.;
+ bl[0] = 6.6;
+ bl[1] = TMath::Pi();
+ Double_t err;
+
+ Double_t y = 2. * fWStarfi->IntegralMultiple(2, al, bl, 0.001, err);
+ printf("WStaa: %f %f %f\n", b, y, err);
+*/
+//
+// MC Integration
+//
+ Double_t y = 0;
+ for (Int_t i = 0; i < 100000; i++)
+ {
+ Double_t phi = TMath::Pi() * gRandom->Rndm();
+ Double_t b1 = fbMax * gRandom->Rndm();
+ y += fWStarfi->Eval(b1, phi);
+ }
+ y *= 2. * 0.1 * 208. * 208. * TMath::Pi() * fbMax / 100000.;
+ return y;
+}
+
+Double_t AliFastGlauber::WSgeo(Double_t* x, Double_t* par)
+{
+//
+// Geometrical Cross-Section
+//
+ Double_t b = x[0];
+ Double_t taa = fWStaa->Eval(b);
+ const Double_t sigma = 55.6; // mbarn
+
+ Double_t y = 2. * TMath::Pi() * b * (1. - TMath::Exp(- sigma * taa)); // fm
+ return y;
+}
+
+
+Double_t AliFastGlauber::WSbinary(Double_t* x, Double_t* par)
+{
+//
+// Geometrical Cross-Section
+//
+ Double_t b = x[0];
+ Double_t sigma = par[0];
+ Double_t taa = fWStaa->Eval(b);
+
+ Double_t y = 2. * TMath::Pi() * b * sigma * taa; // fm
+ return y;
+}
+
+Double_t AliFastGlauber::WSN(Double_t* x, Double_t* par)
+{
+//
+// Geometrical Cross-Section
+//
+ Double_t b = x[0];
+ Double_t y = fWSbinary->Eval(b)/fWSgeo->Eval(b);
+ return y;
+}
+
+void AliFastGlauber::SimulateTrigger(Int_t n)
+{
+ //
+ // Simulates Trigger
+ //
+ TH1F* mbtH = new TH1F("mbtH", "MB Trigger b-Distribution", 100, 0., 20.);
+ TH1F* hdtH = new TH1F("hdtH", "Hard Trigger b-Distribution", 100, 0., 20.);
+ TH1F* mbmH = new TH1F("mbmH", "MB Trigger Multiplicity Distribution", 100, 0., 8000.);
+ TH1F* hdmH = new TH1F("hdmH", "Hard Trigger Multiplicity Distribution", 100, 0., 8000.);
+
+ mbtH->SetXTitle("b [fm]");
+ hdtH->SetXTitle("b [fm]");
+ mbmH->SetXTitle("Multiplicity");
+ hdmH->SetXTitle("Multiplicity");
+
+ TCanvas *c0 = new TCanvas("c0","Trigger Simulation",400,10,600,700);
+ c0->Divide(2,1);
+ TCanvas *c1 = new TCanvas("c1","Trigger Simulation",400,10,600,700);
+ c1->Divide(1,2);
+
+ //
+ //
+ Init(1);
+ for (Int_t iev = 0; iev < n; iev++)
+ {
+ Float_t b, p, mult;
+ GetRandom(b, p, mult);
+ mbtH->Fill(b,1.);
+ hdtH->Fill(b, p);
+ mbmH->Fill(mult, 1.);
+ hdmH->Fill(mult, p);
+
+ c0->cd(1);
+ mbtH->Draw();
+ c0->cd(2);
+ hdtH->Draw();
+ c0->Update();
+
+ c1->cd(1);
+ mbmH->Draw();
+ c1->cd(2);
+ hdmH->Draw();
+ c1->Update();
+ }
+}
+
+void AliFastGlauber::GetRandom(Float_t& b, Float_t& p, Float_t& mult)
+{
+ //
+ // Gives back a random impact parameter, hard trigger probability and multiplicity
+ //
+ b = fWSgeo->GetRandom();
+ Float_t mu = fWSN->Eval(b);
+ p = 1.-TMath::Exp(-mu);
+ mult = 6000./fWSN->Eval(1.) * mu;
+}
+
+Float_t AliFastGlauber::GetRandomImpactParameter(Float_t bmin, Float_t bmax)
+{
+ //
+ // Gives back a random impact parameter in the range bmin .. bmax
+ //
+
+ Float_t b = -1.;
+ while(b < bmin || b > bmax)
+ b = fWSgeo->GetRandom();
+ return b;
+}
+
+Double_t AliFastGlauber::CrossSection(Double_t b1, Double_t b2)
+{
+ //
+ // Return cross-section integrated from b1 to b2
+ //
+
+ return fWSgeo->Integral(b1, b2)/100.;
+}
+
+Double_t AliFastGlauber::FractionOfHardCrossSection(Double_t b1, Double_t b2)
+{
+ //
+ // Return raction of hard cross-section integrated from b1 to b2
+ //
+
+ return fWSbinary->Integral(b1, b2)/fWSbinary->Integral(0., 100.);
+}
+
+
+Double_t AliFastGlauber::Binaries(Double_t b)
+{
+ //
+ // Return number of binary collisions normalized to 1 at b=0
+ //
+
+ return fWSN->Eval(b)/fWSN->Eval(1.);
+}
--- /dev/null
+#ifndef ALIFASTGLAUBER_H
+#define ALIFASTGLAUBER_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+/* $Id$ */
+
+#include <TObject.h>
+class TF1;
+class TF2;
+
+class AliFastGlauber : public TObject {
+ public:
+ AliFastGlauber();
+ virtual ~AliFastGlauber(){;}
+ void SetWoodSaxonParameters(Double_t r0, Double_t d, Double_t w, Double_t n)
+ {fWSr0 = r0; fWSd = d; fWSw = w; fWSn = n;}
+ void SetMaxImpact(Float_t bmax = 20.) {fbMax = bmax;};
+ void SetHardCrossSection(Float_t xs = 6.6) {fSigmaHard = xs;}
+
+ static Double_t WSb (Double_t *xx, Double_t *par);
+ static Double_t WSbz (Double_t *xx, Double_t *par);
+ static Double_t WSz (Double_t *xx, Double_t *par);
+ static Double_t WSta (Double_t *xx, Double_t *par);
+ static Double_t WStarfi (Double_t *xx, Double_t *par);
+ static Double_t WStaa (Double_t *xx, Double_t *par);
+ static Double_t WSgeo (Double_t *xx, Double_t *par);
+ static Double_t WSbinary (Double_t *xx, Double_t *par);
+ static Double_t WSN (Double_t *xx, Double_t *par);
+ void Init(Int_t mode = 0);
+ void DrawWSb();
+ void DrawThickness();
+ void DrawOverlap();
+ void DrawGeo();
+ void DrawBinary();
+ void DrawN();
+ void DrawKernel();
+ Double_t CrossSection(Double_t b1, Double_t b2);
+ Double_t FractionOfHardCrossSection(Double_t b1, Double_t b2);
+ Double_t Binaries(Double_t b);
+ TF2* Kernel() {return fWStarfi;}
+ TF1* Overlap() {return fWStaa;}
+ void SimulateTrigger(Int_t n);
+ void GetRandom(Float_t& b, Float_t& p, Float_t& mult);
+ Float_t GetRandomImpactParameter(Float_t bmin, Float_t bmax);
+ protected:
+ static TF1* fWSb; // Wood-Saxon Function (b)
+ static TF2* fWSbz; // Wood-Saxon Function (b, z)
+ static TF1* fWSz; // Wood-Saxon Function (b = b0, z)
+ static TF1* fWSta; // Thickness Function
+ static TF2* fWStarfi; // Kernel for Overlap Function
+ static TF1* fWStaa; // Overlap Function
+ static TF1* fWSgeo; // dSigma/db geometric
+ static TF1* fWSbinary; // dSigma/db binary
+ static TF1* fWSN; // dN/db binary
+
+ Float_t fWSr0; // Wood-Saxon Parameter r0
+ Float_t fWSd; // Wood-Saxon Parameter d
+ Float_t fWSw; // Wood-Saxon Parameter w
+ Float_t fWSn; // Wood-Saxon Parameter n
+ Float_t fSigmaHard; // Hard Cross Section
+ static Float_t fbMax; // Maximum Impact Parameter
+
+ ClassDef(AliFastGlauber,1) // Event geometry simulation in the Glauber Model
+};
+
+#endif
+
+
+
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff