fix macro control plots npart eccentricity
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowTools / glauberMC / AliGlauberMC.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 //
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>
33 #include <TTree.h>
34 #include <TF1.h>
35 #include <TH1F.h>
36 #include <TArray.h>
37
38 #include "AliGlauberNucleon.h"
39 #include "AliGlauberNucleus.h"
40 #include "AliGlauberMC.h"
41
42 ClassImp(AliGlauberMC)
43
44 //______________________________________________________________________________
45 AliGlauberMC::AliGlauberMC(Option_t* NA, Option_t* NB, Double_t xsect) :
46   TNamed(),
47   fANucleus(NA),
48   fBNucleus(NB),
49   fXSect(xsect),
50   fNucleonsA(0),
51   fNucleonsB(0),
52   fAN(0),
53   fBN(0),
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),
74   fBMin(0.),
75   fBMax(20.),
76   fMaxNpartFound(0),
77   fNpart(0),
78   fNcoll(0),
79   fSx2(0.),
80   fSy2(0.),
81   fSxy(0.),
82   fSx2Coll(0.),
83   fSy2Coll(0.),
84   fSxyColl(0.),
85   fX(0.13),
86   fNpp(8.),
87   fDoPartProd(1)
88 {
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;
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;
102   fdNdEtaTwoNBDParam[5] = 0.13;
103
104   SetName(Form("Glauber_%s_%s",fANucleus.GetName(),fBNucleus.GetName()));
105   SetTitle(Form("Glauber %s+%s Version",fANucleus.GetName(),fBNucleus.GetName()));
106 }
107
108 //______________________________________________________________________________
109 AliGlauberMC::~AliGlauberMC()
110 {
111   //dtor
112   delete fnt;
113 }
114
115 //______________________________________________________________________________
116 AliGlauberMC::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),
157   fNpp(in.fNpp),
158   fDoPartProd(1)
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 //______________________________________________________________________________
168 AliGlauberMC& 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;
216 }
217
218 //______________________________________________________________________________
219 Bool_t AliGlauberMC::CalcEvent(Double_t bgen)
220 {
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   {
227     AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->UncheckedAt(i));
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   {
235     AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->UncheckedAt(i));
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   {
245     AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->UncheckedAt(i));
246     for (Int_t j = 0 ; j < fAN ; j++)
247     {
248       AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->UncheckedAt(j));
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();
256       }
257     }
258   }
259
260   return CalcResults(bgen);
261 }
262
263 //______________________________________________________________________________
264 Bool_t AliGlauberMC::CalcResults(Double_t bgen)
265 {
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;
276   fMeanXColl=0;
277   fMeanYColl=0;
278   fMeanX2Coll=0;
279   fMeanY2Coll=0;
280   fMeanXYColl=0;
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   {
290     AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->UncheckedAt(i));
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   {
311     AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->UncheckedAt(i));
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;
420 }
421
422 //______________________________________________________________________________
423 void AliGlauberMC::Draw(Option_t* /*option*/)
424 {
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);
435 }
436
437 //______________________________________________________________________________
438 Double_t AliGlauberMC::GetTotXSect() const
439 {
440   return (1.*fEvents/fTotalEvents)*TMath::Pi()*fBMax*fBMax/100;
441 }
442
443 //______________________________________________________________________________
444 Double_t AliGlauberMC::GetTotXSectErr() const
445 {
446   return GetTotXSect()/TMath::Sqrt((Double_t)fEvents) *
447          TMath::Sqrt(Double_t(1.-fEvents/fTotalEvents));
448 }
449
450 //______________________________________________________________________________
451 TObjArray *AliGlauberMC::GetNucleons()
452 {
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   {
460     allnucleons->Add(fNucleonsA->UncheckedAt(i));
461   }
462   for (Int_t i = 0; i<fBN; i++)
463   {
464     allnucleons->Add(fNucleonsB->UncheckedAt(i));
465   }
466   return allnucleons;
467 }
468 //______________________________________________________________________________
469 Double_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 //______________________________________________________________________________
481 Int_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 //______________________________________________________________________________
495 Int_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
513 //______________________________________________________________________________
514 void 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;
519 }
520
521 //______________________________________________________________________________
522 void 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 }
529
530 //______________________________________________________________________________
531 void 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 }
538
539 //______________________________________________________________________________
540 void 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)
546 {
547   // set parameters used for calculating multiplicity see GetdNdEtaTwoNBD() for commments
548   fdNdEtaTwoNBDParam[0]=k1;
549   fdNdEtaTwoNBDParam[1]=nmean1;
550   fdNdEtaTwoNBDParam[2]=k2;
551   fdNdEtaTwoNBDParam[3]=nmean2;
552   fdNdEtaTwoNBDParam[4]=alpha;
553   fdNdEtaTwoNBDParam[5]=beta;
554 }
555
556 //______________________________________________________________________________
557 Double_t AliGlauberMC::GetdNdEta(Double_t nnp, Double_t x)
558 {
559   //Get particle density per unit of rapidity
560   //using two component model
561   //Parameters: npp, x
562   return nnp*((1.-x)*fNpart/2.+x*fNcoll);
563 }
564
565 //______________________________________________________________________________
566 Double_t AliGlauberMC::GetdNdEtaGBW( Double_t delta,
567                                      Double_t lambda,
568                                      Double_t  snn)
569 {
570   //Get particle density per unit of rapidity
571   //using the GBW model
572   //Parameters: delta, lambda, snn
573   return fNpart*0.47*TMath::Sqrt(TMath::Power(snn,lambda))
574          * TMath::Power(fNpart,(1.-delta)/3./delta);
575 }
576
577 //_______________________________________________________________________________
578 Double_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 }
597
598 //______________________________________________________________________________
599 Double_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 )
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;
626 }
627
628 //______________________________________________________________________________
629 Double_t AliGlauberMC::GetEccentricityPart() const
630 {
631   //get participant eccentricity of participants
632   if (fNpart<2) return 0.0;
633   return (TMath::Sqrt((fSy2-fSx2)*(fSy2-fSx2)+4*fSxy*fSxy)/(fSy2+fSx2));
634 }
635
636 //______________________________________________________________________________
637 Double_t AliGlauberMC::GetEccentricityPartColl() const
638 {
639   //get participant eccentricity of binary collisions
640   if (fNcoll<2) return 0.0;
641   return (TMath::Sqrt((fSy2Coll-fSx2Coll)*(fSy2Coll-fSx2Coll)+4*fSxyColl*fSxyColl)/(fSy2Coll+fSx2Coll));
642 }
643
644 //______________________________________________________________________________
645 Double_t AliGlauberMC::GetEccentricity() const
646 {
647   //get standard eccentricity of participants
648   if (fNpart<2) return 0.0;
649   return ((fSy2-fSx2)/(fSy2+fSx2));
650 }
651
652 //______________________________________________________________________________
653 Double_t AliGlauberMC::GetEccentricityColl() const
654 {
655   //get standard eccentricity of binary collisions
656   if (fNcoll<2) return 0.0;
657   return ((fSy2Coll-fSx2Coll)/(fSy2Coll+fSx2Coll));
658 }
659
660 //______________________________________________________________________________
661 Bool_t AliGlauberMC::NextEvent(Double_t bgen)
662 {
663   //Make a new event
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     }
672     if ( (succes=CalcEvent(bgen)) ) break; //ends if we have particparts
673   }
674   return succes;
675 }
676
677 //______________________________________________________________________________
678 void AliGlauberMC::Run(Int_t nevents)
679 {
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,
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");
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     }
699
700     q++;
701     Float_t v[27];
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();
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;
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;
751 }
752
753 //---------------------------------------------------------------------------------
754 void AliGlauberMC::runAndSaveNtuple( Int_t n,
755                                      Option_t *sysA,
756                                      Option_t *sysB,
757                                      Double_t signn,
758                                      Double_t mind,
759                                      Double_t r,
760                                      Double_t a,
761                                      const char *fname)
762 {
763   AliGlauberMC mcg(sysA,sysB,signn);
764   mcg.SetMinDistance(mind);
765   mcg.Setr(r);
766   mcg.Seta(a);
767   mcg.Run(n);
768   TNtuple  *nt=mcg.GetNtuple();
769   TFile out(fname,"recreate",fname,9);
770   if(nt) nt->Write();
771   printf("total cross section with a nucleon-nucleon cross section \t%f is \t%f",signn,mcg.GetTotXSect());
772   out.Close();
773 }
774
775 //---------------------------------------------------------------------------------
776 void AliGlauberMC::runAndSaveNucleons( Int_t n,
777                                        Option_t *sysA,
778                                        Option_t *sysB,
779                                        Double_t signn,
780                                        Double_t mind,
781                                        Bool_t verbose,
782                                        const char *fname)
783 {
784   AliGlauberMC mcg(sysA,sysB,signn);
785   mcg.SetMinDistance(mind);
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
794     while(!mcg.NextEvent()) {}
795
796     //access, save and (if wanted) print out nucleons
797     TObjArray* nucleons=mcg.GetNucleons();
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;
805       cout<<"B = "<<mcg.GetB()<<"  Npart = "<<mcg.GetNpart()<<endl<<endl;
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       {
810         AliGlauberNucleon *nucl=(AliGlauberNucleon *)nucleons->UncheckedAt(iNucl);
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);
818       }
819     }
820   }
821   if(out) delete out;
822 }
823
824 //---------------------------------------------------------------------------------
825 void AliGlauberMC::Reset()
826 {
827   //delete the ntuple
828   delete fnt; fnt=NULL;
829 }