[u/mrichter/AliRoot.git] / FASTSIM / AliFastGlauber.cxx

/**************************************************************************
 * 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$
Revision 1.3  2003/04/21 09:35:53  morsch
Protection against division by 0 in Binaries().

Revision 1.2  2003/04/14 14:23:44  morsch
Correction in Binaries().

Revision 1.1  2003/04/10 08:48:13  morsch
First commit.

*/

// 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)
{
//
//  Number of binary collisions
//
    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)
{
//
//  Number of hard processes per event
//
    Double_t b     = x[0];
    Double_t y;
    if (b != 0.) {
	y = fWSbinary->Eval(b)/fWSgeo->Eval(b);
    } else {
	y = fWSbinary->Eval(1.e-4)/fWSgeo->Eval(1.e-4);
    }
    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(0.) * mu;
}

void AliFastGlauber::GetRandom(Int_t& bin, Bool_t& hard)
{
    //
    // Gives back a random impact parameter bin, and hard trigger decission
    //
	Float_t b  = fWSgeo->GetRandom();
	Float_t mu = fWSN->Eval(b) * fSigmaHard;
	Float_t p  = 1.-TMath::Exp(-mu);
	if (b < 5.) {
	    bin = 1;
	} else if (b <  8.6) {
	    bin = 2;
	} else if (b < 11.2) {
	    bin = 3;
	} else if (b < 13.2) {
	    bin = 4;
	} else {
	    bin = 5;
	}
	
	hard = kFALSE;
	
	Float_t r = gRandom->Rndm();
	
	if (r < p) hard = kTRUE;
}


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