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 }