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