]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/FLOW/AliFlowTools/glauberMC/AliGlauberMC.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Do not add VirtualMC to the Configuration folder, it is added to Event from AliMC
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowTools / glauberMC / AliGlauberMC.cxx
CommitLineData
ec852657 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//
18// Glauber MC implementation
19//
20// origin: PHOBOS experiment
21// alification: Mikolaj Krzewicki, Nikhef, mikolaj.krzewicki@cern.ch
22//
23////////////////////////////////////////////////////////////////////////////////
24
25#include <Riostream.h>
26#include <TMath.h>
27#include <TEllipse.h>
28#include <TRandom.h>
29#include <TNamed.h>
30#include <TObjArray.h>
31#include <TNtuple.h>
32#include <TFile.h>
bb442091 33#include <TTree.h>
34#include <TF1.h>
35#include <TH1F.h>
36#include <TArray.h>
ec852657 37
38#include "AliGlauberNucleon.h"
39#include "AliGlauberNucleus.h"
40#include "AliGlauberMC.h"
41
42ClassImp(AliGlauberMC)
43
44//______________________________________________________________________________
45AliGlauberMC::AliGlauberMC(Option_t* NA, Option_t* NB, Double_t xsect) :
566fc3d0 46 TNamed(),
bb442091 47 fANucleus(NA),
48 fBNucleus(NB),
566fc3d0 49 fXSect(xsect),
bb442091 50 fNucleonsA(0),
51 fNucleonsB(0),
566fc3d0 52 fAN(0),
53 fBN(0),
bb442091 54 fnt(0),
55 fMeanX2(0),
56 fMeanY2(0),
57 fMeanXY(0),
58 fMeanXParts(0),
59 fMeanYParts(0),
60 fMeanXColl(0),
61 fMeanYColl(0),
62 fMeanX2Coll(0),
63 fMeanY2Coll(0),
64 fMeanXYColl(0),
65 fMeanXSystem(0),
66 fMeanYSystem(0),
43215add 67 fMeanXA(0),
68 fMeanYA(0),
69 fMeanXB(0),
70 fMeanYB(0),
71 fBMC(0),
bb442091 72 fEvents(0),
73 fTotalEvents(0),
566fc3d0 74 fBMin(0.),
75 fBMax(20.),
bb442091 76 fMaxNpartFound(0),
77 fNpart(0),
78 fNcoll(0),
566fc3d0 79 fSx2(0.),
80 fSy2(0.),
81 fSxy(0.),
82 fSx2Coll(0.),
83 fSy2Coll(0.),
84 fSxyColl(0.),
bb442091 85 fX(0.13),
14244e14 86 fNpp(8.),
fee95d0f 87 fDoPartProd(kFALSE)
ec852657 88{
43215add 89 //ctor
bb442091 90 fdNdEtaParam[0] = 8.0;
91 fdNdEtaParam[1] = 0.13;
92 fdNdEtaGBWParam[0] = 0.79;
93 fdNdEtaGBWParam[1] = 0.288;
94 fdNdEtaGBWParam[2] = 30.25;
95 fdNdEtaNBDParam[0] = 3.0;
96 fdNdEtaNBDParam[1] = 4.0;
97 fdNdEtaNBDParam[2] = 0.13;
566fc3d0 98 fdNdEtaTwoNBDParam[0] = 3.0;
99 fdNdEtaTwoNBDParam[1] = 4.0;
100 fdNdEtaTwoNBDParam[2] = 2.0;
101 fdNdEtaTwoNBDParam[3] = 11.0;
102 fdNdEtaTwoNBDParam[4] = 0.4;
bb442091 103 fdNdEtaTwoNBDParam[5] = 0.13;
104
566fc3d0 105 SetName(Form("Glauber_%s_%s",fANucleus.GetName(),fBNucleus.GetName()));
106 SetTitle(Form("Glauber %s+%s Version",fANucleus.GetName(),fBNucleus.GetName()));
107}
108
109//______________________________________________________________________________
110AliGlauberMC::~AliGlauberMC()
111{
112 //dtor
113 delete fnt;
114}
115
116//______________________________________________________________________________
117AliGlauberMC::AliGlauberMC(const AliGlauberMC& in):
118 TNamed(in),
119 fANucleus(in.fANucleus),
120 fBNucleus(in.fBNucleus),
121 fXSect(in.fXSect),
122 fNucleonsA(in.fNucleonsA),
123 fNucleonsB(in.fNucleonsB),
124 fAN(in.fAN),
125 fBN(in.fBN),
126 fnt(in.fnt),
127 fMeanX2(in.fMeanX2),
128 fMeanY2(in.fMeanY2),
129 fMeanXY(in.fMeanXY),
130 fMeanXParts(in.fMeanXParts),
131 fMeanYParts(in.fMeanYParts),
132 fMeanXColl(in.fMeanXColl),
133 fMeanYColl(in.fMeanYColl),
134 fMeanX2Coll(in.fMeanX2Coll),
135 fMeanY2Coll(in.fMeanY2Coll),
136 fMeanXYColl(in.fMeanXYColl),
137 fMeanXSystem(in.fMeanXSystem),
138 fMeanYSystem(in.fMeanYSystem),
43215add 139 fMeanXA(in.fMeanXA),
140 fMeanYA(in.fMeanYA),
141 fMeanXB(in.fMeanXB),
142 fMeanYB(in.fMeanYB),
143 fBMC(in.fBMC),
566fc3d0 144 fEvents(in.fEvents),
145 fTotalEvents(in.fTotalEvents),
146 fBMin(in.fBMin),
147 fBMax(in.fBMax),
148 fMaxNpartFound(in.fMaxNpartFound),
149 fNpart(in.fNpart),
150 fNcoll(in.fNcoll),
151 fSx2(in.fSx2),
152 fSy2(in.fSy2),
153 fSxy(in.fSxy),
154 fSx2Coll(in.fSx2Coll),
155 fSy2Coll(in.fSy2Coll),
156 fSxyColl(in.fSxyColl),
157 fX(in.fX),
97d98e48 158 fNpp(in.fNpp),
fee95d0f 159 fDoPartProd(kFALSE)
566fc3d0 160{
161 //copy ctor
162 memcpy(fdNdEtaParam,in.fdNdEtaParam,2*sizeof(Double_t));
163 memcpy(fdNdEtaGBWParam,in.fdNdEtaGBWParam,3*sizeof(Double_t));
164 memcpy(fdNdEtaNBDParam,in.fdNdEtaNBDParam,3*sizeof(Double_t));
165 memcpy(fdNdEtaTwoNBDParam,in.fdNdEtaTwoNBDParam,6*sizeof(Double_t));
166}
167
168//______________________________________________________________________________
169AliGlauberMC& AliGlauberMC::operator=(const AliGlauberMC& in)
170{
171 //assignment
172 fANucleus=in.fANucleus;
173 fBNucleus=in.fBNucleus;
174 fXSect=in.fXSect;
175 fNucleonsA=in.fNucleonsA;
176 fNucleonsB=in.fNucleonsB;
177 fAN=in.fAN;
178 fBN=in.fBN;
179 fnt=in.fnt;
180 fMeanX2=in.fMeanX2;
181 fMeanY2=in.fMeanY2;
182 fMeanXY=in.fMeanXY;
183 fMeanXParts=in.fMeanXParts;
184 fMeanYParts=in.fMeanYParts;
185 fMeanXColl=in.fMeanXColl;
186 fMeanYColl=in.fMeanYColl;
187 fMeanX2Coll=in.fMeanX2Coll;
188 fMeanY2Coll=in.fMeanY2Coll;
189 fMeanXYColl=in.fMeanXYColl;
190 fMeanXSystem=in.fMeanXSystem;
191 fMeanYSystem=in.fMeanYSystem;
43215add 192 fMeanXA=in.fMeanXA;
193 fMeanYA=in.fMeanYA;
194 fMeanXB=in.fMeanXB;
195 fMeanYB=in.fMeanYB;
196 fBMC=in.fBMC;
566fc3d0 197 fEvents=in.fEvents;
198 fTotalEvents=in.fTotalEvents;
199 fBMin=in.fBMin;
200 fBMax=in.fBMax;
201 memcpy(fdNdEtaParam,in.fdNdEtaParam,2*sizeof(Double_t));
202 memcpy(fdNdEtaGBWParam,in.fdNdEtaGBWParam,3*sizeof(Double_t));
203 memcpy(fdNdEtaNBDParam,in.fdNdEtaNBDParam,3*sizeof(Double_t));
204 memcpy(fdNdEtaTwoNBDParam,in.fdNdEtaTwoNBDParam,6*sizeof(Double_t));
205 fMaxNpartFound=in.fMaxNpartFound;
206 fNpart=in.fNpart;
207 fNcoll=in.fNcoll;
208 fSx2=in.fSx2;
209 fSy2=in.fSy2;
210 fSxy=in.fSxy;
211 fSx2Coll=in.fSx2Coll;
212 fSy2Coll=in.fSy2Coll;
213 fSxyColl=in.fSxyColl;
214 fX=in.fX;
215 fNpp=in.fNpp;
216 return *this;
ec852657 217}
218
219//______________________________________________________________________________
220Bool_t AliGlauberMC::CalcEvent(Double_t bgen)
221{
bb442091 222 // prepare event
223 fANucleus.ThrowNucleons(-bgen/2.);
224 fNucleonsA = fANucleus.GetNucleons();
225 fAN = fANucleus.GetN();
226 for (Int_t i = 0; i<fAN; i++)
227 {
14244e14 228 AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->UncheckedAt(i));
bb442091 229 nucleonA->SetInNucleusA();
230 }
231 fBNucleus.ThrowNucleons(bgen/2.);
232 fNucleonsB = fBNucleus.GetNucleons();
233 fBN = fBNucleus.GetN();
234 for (Int_t i = 0; i<fBN; i++)
235 {
14244e14 236 AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->UncheckedAt(i));
bb442091 237 nucleonB->SetInNucleusB();
238 }
239
240 // "ball" diameter = distance at which two balls interact
241 Double_t d2 = (Double_t)fXSect/(TMath::Pi()*10); // in fm^2
242
243 // for each of the A nucleons in nucleus B
244 for (Int_t i = 0; i<fBN; i++)
245 {
14244e14 246 AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->UncheckedAt(i));
bb442091 247 for (Int_t j = 0 ; j < fAN ; j++)
248 {
14244e14 249 AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->UncheckedAt(j));
bb442091 250 Double_t dx = nucleonB->GetX()-nucleonA->GetX();
251 Double_t dy = nucleonB->GetY()-nucleonA->GetY();
252 Double_t dij = dx*dx+dy*dy;
253 if (dij < d2)
254 {
255 nucleonB->Collide();
256 nucleonA->Collide();
ec852657 257 }
bb442091 258 }
ec852657 259 }
bb442091 260
ec852657 261 return CalcResults(bgen);
262}
263
264//______________________________________________________________________________
265Bool_t AliGlauberMC::CalcResults(Double_t bgen)
266{
bb442091 267 // calc results for the given event
268 //return true if we have participants
269
270 fNpart=0;
271 fNcoll=0;
272 fMeanX2=0;
273 fMeanY2=0;
274 fMeanXY=0;
275 fMeanXParts=0;
276 fMeanYParts=0;
a0278ee7 277 fMeanXColl=0;
278 fMeanYColl=0;
279 fMeanX2Coll=0;
280 fMeanY2Coll=0;
281 fMeanXYColl=0;
bb442091 282 fMeanXSystem=0;
283 fMeanYSystem=0;
43215add 284 fMeanXA=0;
285 fMeanYA=0;
286 fMeanXB=0;
287 fMeanYB=0;
bb442091 288
289 for (Int_t i = 0; i<fAN; i++)
290 {
14244e14 291 AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->UncheckedAt(i));
bb442091 292 Double_t xA=nucleonA->GetX();
293 Double_t yA=nucleonA->GetY();
294 fMeanXSystem += xA;
295 fMeanYSystem += yA;
43215add 296 fMeanXA += xA;
297 fMeanYA += yA;
bb442091 298
299 if(nucleonA->IsWounded())
300 {
301 fNpart++;
302 fMeanXParts += xA;
303 fMeanYParts += yA;
304 fMeanX2 += xA * xA;
305 fMeanY2 += yA * yA;
306 fMeanXY += xA * yA;
307 }
308 }
309
310 for (Int_t i = 0; i<fBN; i++)
311 {
14244e14 312 AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->UncheckedAt(i));
bb442091 313 Double_t xB=nucleonB->GetX();
314 Double_t yB=nucleonB->GetY();
315 fMeanXSystem += xB;
316 fMeanYSystem += yB;
43215add 317 fMeanXB += xB;
318 fMeanYB += yB;
bb442091 319
320 if(nucleonB->IsWounded())
321 {
322 Int_t ncoll = nucleonB->GetNColl();
323 fNpart++;
324 fMeanXParts += xB;
325 fMeanXColl += xB*ncoll;
326 fMeanYParts += yB;
327 fMeanYColl += yB*ncoll;
328 fMeanX2 += xB * xB;
329 fMeanX2Coll += xB*xB*ncoll*ncoll;
330 fMeanY2 += yB * yB;
331 fMeanY2Coll += yB*yB*ncoll*ncoll;
332 fMeanXY += xB * yB;
333 fMeanXYColl += xB*yB*ncoll*ncoll;
334 fNcoll += nucleonB->GetNColl();
335 }
336 }
337
338 if (fNpart>0)
339 {
340 fMeanXParts /= fNpart;
341 fMeanYParts /= fNpart;
342 fMeanX2 /= fNpart;
343 fMeanY2 /= fNpart;
344 fMeanXY /= fNpart;
345 }
346 else
347 {
348 fMeanXParts = 0;
349 fMeanYParts = 0;
350 fMeanX2 = 0;
351 fMeanY2 = 0;
352 fMeanXY = 0;
353 }
354
355 if (fNcoll>0)
356 {
357 fMeanXColl /= fNcoll;
358 fMeanYColl /= fNcoll;
359 fMeanX2Coll /= fNcoll;
360 fMeanY2Coll /= fNcoll;
361 fMeanXYColl /= fNcoll;
362 }
363 else
364 {
365 fMeanXColl = 0;
366 fMeanYColl = 0;
367 fMeanX2Coll = 0;
368 fMeanY2Coll = 0;
369 fMeanXYColl = 0;
370 }
371
372 if(fAN+fBN>0)
373 {
374 fMeanXSystem /= (fAN + fBN);
375 fMeanYSystem /= (fAN + fBN);
376 }
377 else
378 {
379 fMeanXSystem = 0;
380 fMeanYSystem = 0;
381 }
382 if(fAN>0)
383 {
43215add 384 fMeanXA /= fAN;
385 fMeanYA /= fAN;
bb442091 386 }
387 else
388 {
43215add 389 fMeanXA = 0;
390 fMeanYA = 0;
bb442091 391 }
392
393 if(fBN>0)
394 {
43215add 395 fMeanXB /= fBN;
396 fMeanYB /= fBN;
bb442091 397 }
398 else
399 {
43215add 400 fMeanXB = 0;
401 fMeanYB = 0;
bb442091 402 }
403
404 fSx2=fMeanX2-(fMeanXParts*fMeanXParts);
405 fSy2=fMeanY2-(fMeanYParts*fMeanYParts);
406 fSxy=fMeanXY-fMeanXParts*fMeanYParts;
407
408 fSx2Coll=fMeanX2Coll-(fMeanXColl*fMeanXColl);
409 fSy2Coll=fMeanY2Coll-(fMeanYColl*fMeanYColl);
410 fSxyColl=fMeanXYColl-fMeanXColl*fMeanYColl;
411
43215add 412 fBMC = bgen;
bb442091 413
414 fTotalEvents++;
415 if (fNpart>0) fEvents++;
416
417 if (fNpart==0) return kFALSE;
418 if (fNpart > fMaxNpartFound) fMaxNpartFound = fNpart;
419
420 return kTRUE;
ec852657 421}
422
423//______________________________________________________________________________
424void AliGlauberMC::Draw(Option_t* /*option*/)
425{
43215add 426 //draw method
bb442091 427 fANucleus.Draw(fXSect, 2);
428 fBNucleus.Draw(fXSect, 4);
429
430 TEllipse e;
431 e.SetFillColor(0);
432 e.SetLineColor(1);
433 e.SetLineStyle(2);
434 e.SetLineWidth(1);
435 e.DrawEllipse(GetB()/2,0,fBNucleus.GetR(),fBNucleus.GetR(),0,360,0);
436 e.DrawEllipse(-GetB()/2,0,fANucleus.GetR(),fANucleus.GetR(),0,360,0);
ec852657 437}
438
439//______________________________________________________________________________
440Double_t AliGlauberMC::GetTotXSect() const
441{
43215add 442 //total xsection
bb442091 443 return (1.*fEvents/fTotalEvents)*TMath::Pi()*fBMax*fBMax/100;
ec852657 444}
445
446//______________________________________________________________________________
447Double_t AliGlauberMC::GetTotXSectErr() const
448{
43215add 449 //total xsection error
bb442091 450 return GetTotXSect()/TMath::Sqrt((Double_t)fEvents) *
451 TMath::Sqrt(Double_t(1.-fEvents/fTotalEvents));
ec852657 452}
453
454//______________________________________________________________________________
bb442091 455TObjArray *AliGlauberMC::GetNucleons()
ec852657 456{
43215add 457 //get array of nucleons
bb442091 458 if(!fNucleonsA || !fNucleonsB) return 0;
459 fNucleonsA->SetOwner(0);
460 fNucleonsB->SetOwner(0);
461 TObjArray *allnucleons=new TObjArray(fAN+fBN);
462 allnucleons->SetOwner();
463 for (Int_t i = 0; i<fAN; i++)
464 {
14244e14 465 allnucleons->Add(fNucleonsA->UncheckedAt(i));
bb442091 466 }
467 for (Int_t i = 0; i<fBN; i++)
468 {
14244e14 469 allnucleons->Add(fNucleonsB->UncheckedAt(i));
bb442091 470 }
471 return allnucleons;
472}
473//______________________________________________________________________________
474Double_t AliGlauberMC::NegativeBinomialDistribution(Int_t x, Int_t k, Double_t nmean)
475{
43215add 476 //produce a number distributed acc. neg.bin.dist
bb442091 477 if(k<=0)
478 {
479 cout << "Error, zero or negative K" << endl;
480 return 0;
481 }
482 return (TMath::Binomial(x+k-1,x))
483 *TMath::Power(((nmean/Double_t(k))/(1+nmean/Double_t(k))),Double_t(x))
484 *(1/(TMath::Power((1+nmean/Double_t(k)),Double_t(k))));
485}
486//______________________________________________________________________________
43215add 487Int_t AliGlauberMC::NegativeBinomialRandom(Int_t k, Double_t nmean) const
bb442091 488{
43215add 489 //return random integer from a Negative Binomial Distribution
bb442091 490 static const Int_t fMaxPlot = 1000;
491 Double_t array[fMaxPlot];
492 array[0] = NegativeBinomialDistribution(0,k,nmean);
493 for (Int_t i=1; i<fMaxPlot; i++)
494 {
495 array[i] = NegativeBinomialDistribution(i,k,nmean) + array[i-1];
496 }
497 Double_t r = gRandom->Uniform(0,1);
498 return TMath::BinarySearch(fMaxPlot,array,r)+2;
499}
500//______________________________________________________________________________
501Int_t AliGlauberMC::DoubleNegativeBinomialRandom( Int_t k,
502 Double_t nmean,
503 Int_t k2,
504 Double_t nmean2,
43215add 505 Double_t alpha ) const
bb442091 506{
507 //return random integer from a Double Negative Binomial Distribution
508 static const Int_t fMaxPlot = 1000;
509 Double_t array[fMaxPlot];
510 array[0] = alpha*NegativeBinomialDistribution(0,k,nmean)+(1-alpha)*NegativeBinomialDistribution(0,k2,nmean2);
511 for (Int_t i=1; i<fMaxPlot; i++)
512 {
513 array[i] = alpha*NegativeBinomialDistribution(i,k,nmean)+(1-alpha)*NegativeBinomialDistribution(i,k2,nmean2) + array[i-1];
514 }
515 Double_t r = gRandom->Uniform(0,1);
516 return TMath::BinarySearch(fMaxPlot,array,r)+2;
517}
518
bb442091 519//______________________________________________________________________________
520void AliGlauberMC::SetdNdEtaParam(Double_t nnp, Double_t x)
521{
522 // set parameters used for calculating multiplicity, see GetdNdEta() for commments
523 fdNdEtaParam[0]=nnp;
524 fdNdEtaParam[1]=x;
ec852657 525}
526
1b4934d5 527//______________________________________________________________________________
bb442091 528void AliGlauberMC::SetdNdEtaGBWParam(Double_t delta, Double_t lambda, Double_t snn)
529{
530 // set parameters used for calculating multiplicity see GetdNdEtaGBW() for commments
531 fdNdEtaGBWParam[0]=delta;
532 fdNdEtaGBWParam[1]=lambda;
533 fdNdEtaGBWParam[2]=snn;
534}
14244e14 535
bb442091 536//______________________________________________________________________________
537void AliGlauberMC::SetdNdEtaNBDParam(Double_t k, Double_t nmean, Double_t beta)
538{
539 // set parameters used for calculating multiplicity see GetdNdEtaNBD() for commments
540 fdNdEtaNBDParam[0]=k;
541 fdNdEtaNBDParam[1]=nmean;
542 fdNdEtaNBDParam[2]=beta;
543}
14244e14 544
bb442091 545//______________________________________________________________________________
566fc3d0 546void AliGlauberMC::SetdNdEtaTwoNBDParam( Double_t k1,
547 Double_t nmean1,
548 Double_t k2,
549 Double_t nmean2,
550 Double_t alpha,
551 Double_t beta)
bb442091 552{
553 // set parameters used for calculating multiplicity see GetdNdEtaTwoNBD() for commments
566fc3d0 554 fdNdEtaTwoNBDParam[0]=k1;
555 fdNdEtaTwoNBDParam[1]=nmean1;
556 fdNdEtaTwoNBDParam[2]=k2;
557 fdNdEtaTwoNBDParam[3]=nmean2;
558 fdNdEtaTwoNBDParam[4]=alpha;
bb442091 559 fdNdEtaTwoNBDParam[5]=beta;
560}
14244e14 561
bb442091 562//______________________________________________________________________________
43215add 563Double_t AliGlauberMC::GetdNdEta(Double_t nnp, Double_t x) const
1b4934d5 564{
565 //Get particle density per unit of rapidity
bb442091 566 //using two component model
567 //Parameters: npp, x
568 return nnp*((1.-x)*fNpart/2.+x*fNcoll);
0448b9f3 569}
570
571//______________________________________________________________________________
bb442091 572Double_t AliGlauberMC::GetdNdEtaGBW( Double_t delta,
573 Double_t lambda,
43215add 574 Double_t snn) const
0448b9f3 575{
576 //Get particle density per unit of rapidity
577 //using the GBW model
bb442091 578 //Parameters: delta, lambda, snn
579 return fNpart*0.47*TMath::Sqrt(TMath::Power(snn,lambda))
580 * TMath::Power(fNpart,(1.-delta)/3./delta);
1b4934d5 581}
582
14244e14 583//_______________________________________________________________________________
43215add 584Double_t AliGlauberMC::GetdNdEtaNBD ( Int_t k, Double_t nmean, Double_t beta) const
bb442091 585{
586 //Get particle density per unit of rapidity
587 //using a aandomized number from a negative binomial distrubution
588 //Parameters: k = related to distribition width
589 // nmean = mean of distribution
590 // beta = set contribution of participants / binary collisions to multiplicity
591 Double_t mulnp=0.0;
592 for(Int_t i = 0; i<fNpart; i++)
593 {
594 mulnp+=NegativeBinomialRandom(k,nmean);
595 }
596 Double_t mulnb=0.0;
597 for(Int_t i = 0; i<fNcoll; i++)
598 {
599 mulnb+=NegativeBinomialRandom(k,nmean);
600 }
601 return (1-beta)*mulnp/2+beta*mulnb;
602}
bb442091 603
14244e14 604//______________________________________________________________________________
566fc3d0 605Double_t AliGlauberMC::GetdNdEtaTwoNBD ( Int_t k1,
606 Double_t nmean1,
607 Int_t k2,
608 Double_t nmean2,
609 Double_t alpha,
43215add 610 Double_t beta ) const
bb442091 611{
612 //Get particle density per unit of rapidity
613 //using random numbers from two negative binomial distributions
614 //Parameters: k1 = related to distribition width of distribution 1
615 // nmean1 = mean of distribution 1
616 // k2 = related to distribition width of distribution 2
617 // nmean2 = mean of distribution 2
618 // alpha = set contributions of distrubitin 1 / distribution 2
619 // beta = set contribution of participants / binary collisions to multiplicity
620 Double_t mulnp=0.0;
621 for(Int_t i = 0; i<fNpart; i++)
622 {
623 mulnp+=DoubleNegativeBinomialRandom(k1,nmean1,k2,nmean2,alpha);
624 }
625 Double_t mulnb=0.0;
626 for(Int_t i = 0; i<fNcoll; i++)
627 {
628 mulnb+=DoubleNegativeBinomialRandom(k1,nmean1,k2,nmean2,alpha);
629 }
630 Double_t mul=(1-beta)*mulnp/2+beta*mulnb;
631 return mul;
bb442091 632}
14244e14 633
1b4934d5 634//______________________________________________________________________________
635Double_t AliGlauberMC::GetEccentricityPart() const
636{
bb442091 637 //get participant eccentricity of participants
a0278ee7 638 if (fNpart<2) return 0.0;
1b4934d5 639 return (TMath::Sqrt((fSy2-fSx2)*(fSy2-fSx2)+4*fSxy*fSxy)/(fSy2+fSx2));
640}
641
bb442091 642//______________________________________________________________________________
643Double_t AliGlauberMC::GetEccentricityPartColl() const
644{
645 //get participant eccentricity of binary collisions
a0278ee7 646 if (fNcoll<2) return 0.0;
bb442091 647 return (TMath::Sqrt((fSy2Coll-fSx2Coll)*(fSy2Coll-fSx2Coll)+4*fSxyColl*fSxyColl)/(fSy2Coll+fSx2Coll));
648}
649
1b4934d5 650//______________________________________________________________________________
651Double_t AliGlauberMC::GetEccentricity() const
652{
bb442091 653 //get standard eccentricity of participants
74d0740e 654 if (fNpart<2) return 0.0;
1b4934d5 655 return ((fSy2-fSx2)/(fSy2+fSx2));
656}
bb442091 657
658//______________________________________________________________________________
659Double_t AliGlauberMC::GetEccentricityColl() const
660{
661 //get standard eccentricity of binary collisions
b7166644 662 if (fNcoll<2) return 0.0;
bb442091 663 return ((fSy2Coll-fSx2Coll)/(fSy2Coll+fSx2Coll));
664}
14244e14 665
ec852657 666//______________________________________________________________________________
667Bool_t AliGlauberMC::NextEvent(Double_t bgen)
668{
1b4934d5 669 //Make a new event
bb442091 670 Int_t nAttempts = 10; // set indices, max attempts and max comparisons
671 Bool_t succes = kFALSE;
672 for(Int_t j=0; j<nAttempts; j++)
673 {
674 if(bgen<0||!succes) //get impactparameter
675 {
676 bgen = TMath::Sqrt((fBMax*fBMax-fBMin*fBMin)*gRandom->Rndm()+fBMin*fBMin);
677 }
97d98e48 678 if ( (succes=CalcEvent(bgen)) ) break; //ends if we have particparts
bb442091 679 }
680 return succes;
ec852657 681}
682
683//______________________________________________________________________________
684void AliGlauberMC::Run(Int_t nevents)
685{
43215add 686 //example run
bb442091 687 cout << "Generating " << nevents << " events..." << endl;
688 TString name(Form("nt_%s_%s",fANucleus.GetName(),fBNucleus.GetName()));
689 TString title(Form("%s + %s (x-sect = %d mb)",fANucleus.GetName(),fBNucleus.GetName(),(Int_t) fXSect));
690 if (fnt == 0)
691 {
692 fnt = new TNtuple(name,title,
14244e14 693 "Npart:Ncoll:B:MeanX:MeanY:MeanX2:MeanY2:MeanXY:VarX:VarY:VarXY:MeanXSystem:MeanYSystem:MeanXA:MeanYA:MeanXB:MeanYB:VarE:VarEColl:VarEPart:VarEPartColl:dNdEta:dNdEtaGBW:dNdEtaNBD:dNdEtaTwoNBD:xsect:tAA");
bb442091 694 fnt->SetDirectory(0);
695 }
696 Int_t q = 0;
697 Int_t u = 0;
698 for (Int_t i = 0; i<nevents; i++)
699 {
700
701 if(!NextEvent())
702 {
703 u++;
704 continue;
705 }
a0278ee7 706
bb442091 707 q++;
14244e14 708 Float_t v[27];
bb442091 709 v[0] = GetNpart();
710 v[1] = GetNcoll();
43215add 711 v[2] = fBMC;
bb442091 712 v[3] = fMeanXParts;
713 v[4] = fMeanYParts;
714 v[5] = fMeanX2;
715 v[6] = fMeanY2;
716 v[7] = fMeanXY;
717 v[8] = fSx2;
718 v[9] = fSy2;
719 v[10] = fSxy;
720 v[11] = fMeanXSystem;
721 v[12] = fMeanYSystem;
43215add 722 v[13] = fMeanXA;
723 v[14] = fMeanYA;
724 v[15] = fMeanXB;
725 v[16] = fMeanYB;
bb442091 726 v[17] = GetEccentricity();
727 v[18] = GetEccentricityColl();
728 v[19] = GetEccentricityPart();
729 v[20] = GetEccentricityPartColl();
14244e14 730 if (fDoPartProd) {
731 v[21] = GetdNdEta( fdNdEtaParam[0],fdNdEtaParam[1] );
732 v[22] = GetdNdEtaGBW( fdNdEtaGBWParam[0],fdNdEtaGBWParam[1],fdNdEtaGBWParam[2] );
733 v[23] = GetdNdEtaNBD( TMath::Nint(fdNdEtaNBDParam[0]),
734 fdNdEtaNBDParam[1],
735 fdNdEtaNBDParam[2] );
736 v[24] = GetdNdEtaTwoNBD( TMath::Nint(fdNdEtaTwoNBDParam[0]),
737 fdNdEtaTwoNBDParam[1],
738 TMath::Nint(fdNdEtaTwoNBDParam[2]),
739 fdNdEtaTwoNBDParam[3],
740 fdNdEtaTwoNBDParam[4],
741 fdNdEtaTwoNBDParam[5] );
742 } else {
743 v[21] = 0;
744 v[22] = 0;
745 v[23] = 0;
746 v[24] = 0;
747 }
748 v[25]=fXSect;
749
750 Float_t mytAA=-999;
751 if (GetNcoll()>0) mytAA=GetNcoll()/fXSect;
752 v[26]=mytAA;
bb442091 753 fnt->Fill(v);
754
755 if ((i%50)==0) std::cout << "Generating Event # " << i << "... \r" << flush;
756 }
757 std::cout << "Generating Event # " << nevents << "... \r" << endl << "Done! Succesfull events: " << q << " discarded events: " << u <<"."<< endl;
ec852657 758}
759
760//---------------------------------------------------------------------------------
43215add 761void AliGlauberMC::RunAndSaveNtuple( Int_t n,
762 const Option_t *sysA,
763 const Option_t *sysB,
ec852657 764 Double_t signn,
765 Double_t mind,
43215add 766 Double_t r,
767 Double_t a,
ec852657 768 const char *fname)
769{
43215add 770 //example run
566fc3d0 771 AliGlauberMC mcg(sysA,sysB,signn);
772 mcg.SetMinDistance(mind);
14244e14 773 mcg.Setr(r);
774 mcg.Seta(a);
566fc3d0 775 mcg.Run(n);
776 TNtuple *nt=mcg.GetNtuple();
bb442091 777 TFile out(fname,"recreate",fname,9);
778 if(nt) nt->Write();
14244e14 779 printf("total cross section with a nucleon-nucleon cross section \t%f is \t%f",signn,mcg.GetTotXSect());
bb442091 780 out.Close();
ec852657 781}
782
783//---------------------------------------------------------------------------------
43215add 784void AliGlauberMC::RunAndSaveNucleons( Int_t n,
785 const Option_t *sysA,
bb442091 786 Option_t *sysB,
43215add 787 const Double_t signn,
ec852657 788 Double_t mind,
789 Bool_t verbose,
790 const char *fname)
791{
43215add 792 //example run
566fc3d0 793 AliGlauberMC mcg(sysA,sysB,signn);
794 mcg.SetMinDistance(mind);
bb442091 795 TFile *out=0;
796 if(fname)
797 out=new TFile(fname,"recreate",fname,9);
798
799 for(Int_t ievent=0; ievent<n; ievent++)
800 {
801
802 //get an event with at least one collision
566fc3d0 803 while(!mcg.NextEvent()) {}
bb442091 804
805 //access, save and (if wanted) print out nucleons
566fc3d0 806 TObjArray* nucleons=mcg.GetNucleons();
bb442091 807 if(!nucleons) continue;
808 if(out)
809 nucleons->Write(Form("nucleonarray%d",ievent),TObject::kSingleKey);
810
811 if(verbose)
812 {
813 cout<<endl<<endl<<"EVENT NO: "<<ievent<<endl;
566fc3d0 814 cout<<"B = "<<mcg.GetB()<<" Npart = "<<mcg.GetNpart()<<endl<<endl;
bb442091 815 printf("Nucleus\t X\t Y\t Z\tNcoll\n");
816 Int_t nNucls=nucleons->GetEntries();
817 for(Int_t iNucl=0; iNucl<nNucls; iNucl++)
818 {
14244e14 819 AliGlauberNucleon *nucl=(AliGlauberNucleon *)nucleons->UncheckedAt(iNucl);
bb442091 820 Char_t nucleus='A';
821 if(nucl->IsInNucleusB()) nucleus='B';
822 Double_t x=nucl->GetX();
823 Double_t y=nucl->GetY();
824 Double_t z=nucl->GetZ();
825 Int_t ncoll=nucl->GetNColl();
826 printf(" %c\t%2.2f\t%2.2f\t%2.2f\t%3d\n",nucleus,x,y,z,ncoll);
ec852657 827 }
bb442091 828 }
829 }
830 if(out) delete out;
ec852657 831}
832
566fc3d0 833//---------------------------------------------------------------------------------
834void AliGlauberMC::Reset()
835{
836 //delete the ntuple
837 delete fnt; fnt=NULL;
838}