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   fMeanXA(0),
  68   fMeanYA(0),
  69   fMeanXB(0),
  70   fMeanYB(0),
  71   fBMC(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   //ctor
  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;
  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;
 103   fdNdEtaTwoNBDParam[5] = 0.13;
 104
 105   SetName(Form("Glauber_%s_%s",fANucleus.GetName(),fBNucleus.GetName()));
 106   SetTitle(Form("Glauber %s+%s Version",fANucleus.GetName(),fBNucleus.GetName()));
 107 }
 108
 109 //______________________________________________________________________________
 110 AliGlauberMC::~AliGlauberMC()
 111 {
 112   //dtor
 113   delete fnt;
 114 }
 115
 116 //______________________________________________________________________________
 117 AliGlauberMC::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),
 139   fMeanXA(in.fMeanXA),
 140   fMeanYA(in.fMeanYA),
 141   fMeanXB(in.fMeanXB),
 142   fMeanYB(in.fMeanYB),
 143   fBMC(in.fBMC),
 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),
 158   fNpp(in.fNpp),
 159   fDoPartProd(1)
 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 //______________________________________________________________________________
 169 AliGlauberMC& 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;
 192   fMeanXA=in.fMeanXA;
 193   fMeanYA=in.fMeanYA;
 194   fMeanXB=in.fMeanXB;
 195   fMeanYB=in.fMeanYB;
 196   fBMC=in.fBMC;
 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;
 217 }
 218
 219 //______________________________________________________________________________
 220 Bool_t AliGlauberMC::CalcEvent(Double_t bgen)
 221 {
 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   {
 228     AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->UncheckedAt(i));
 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   {
 236     AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->UncheckedAt(i));
 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   {
 246     AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->UncheckedAt(i));
 247     for (Int_t j = 0 ; j < fAN ; j++)
 248     {
 249       AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->UncheckedAt(j));
 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();
 257       }
 258     }
 259   }
 260
 261   return CalcResults(bgen);
 262 }
 263
 264 //______________________________________________________________________________
 265 Bool_t AliGlauberMC::CalcResults(Double_t bgen)
 266 {
 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;
 277   fMeanXColl=0;
 278   fMeanYColl=0;
 279   fMeanX2Coll=0;
 280   fMeanY2Coll=0;
 281   fMeanXYColl=0;
 282   fMeanXSystem=0;
 283   fMeanYSystem=0;
 284   fMeanXA=0;
 285   fMeanYA=0;
 286   fMeanXB=0;
 287   fMeanYB=0;
 288
 289   for (Int_t i = 0; i<fAN; i++)
 290   {
 291     AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->UncheckedAt(i));
 292     Double_t xA=nucleonA->GetX();
 293     Double_t yA=nucleonA->GetY();
 294     fMeanXSystem  += xA;
 295     fMeanYSystem  += yA;
 296     fMeanXA  += xA;
 297     fMeanYA  += yA;
 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   {
 312     AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->UncheckedAt(i));
 313     Double_t xB=nucleonB->GetX();
 314     Double_t yB=nucleonB->GetY();
 315     fMeanXSystem  += xB;
 316     fMeanYSystem  += yB;
 317     fMeanXB  += xB;
 318     fMeanYB  += yB;
 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   {
 384     fMeanXA /= fAN;
 385     fMeanYA /= fAN;
 386   }
 387   else
 388   {
 389     fMeanXA = 0;
 390     fMeanYA = 0;
 391   }
 392
 393   if(fBN>0)
 394   {
 395     fMeanXB /= fBN;
 396     fMeanYB /= fBN;
 397   }
 398   else
 399   {
 400     fMeanXB = 0;
 401     fMeanYB = 0;
 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
 412   fBMC = bgen;
 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;
 421 }
 422
 423 //______________________________________________________________________________
 424 void AliGlauberMC::Draw(Option_t* /*option*/)
 425 {
 426   //draw method
 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);
 437 }
 438
 439 //______________________________________________________________________________
 440 Double_t AliGlauberMC::GetTotXSect() const
 441 {
 442   //total xsection
 443   return (1.*fEvents/fTotalEvents)*TMath::Pi()*fBMax*fBMax/100;
 444 }
 445
 446 //______________________________________________________________________________
 447 Double_t AliGlauberMC::GetTotXSectErr() const
 448 {
 449   //total xsection error
 450   return GetTotXSect()/TMath::Sqrt((Double_t)fEvents) *
 451          TMath::Sqrt(Double_t(1.-fEvents/fTotalEvents));
 452 }
 453
 454 //______________________________________________________________________________
 455 TObjArray *AliGlauberMC::GetNucleons()
 456 {
 457   //get array of nucleons
 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   {
 465     allnucleons->Add(fNucleonsA->UncheckedAt(i));
 466   }
 467   for (Int_t i = 0; i<fBN; i++)
 468   {
 469     allnucleons->Add(fNucleonsB->UncheckedAt(i));
 470   }
 471   return allnucleons;
 472 }
 473 //______________________________________________________________________________
 474 Double_t AliGlauberMC::NegativeBinomialDistribution(Int_t x, Int_t k, Double_t nmean)
 475 {
 476   //produce a number distributed acc. neg.bin.dist
 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 //______________________________________________________________________________
 487 Int_t AliGlauberMC::NegativeBinomialRandom(Int_t k, Double_t nmean) const
 488 {
 489   //return random integer from a Negative Binomial Distribution
 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 //______________________________________________________________________________
 501 Int_t AliGlauberMC::DoubleNegativeBinomialRandom( Int_t k,
 502                                                   Double_t nmean,
 503                                                   Int_t k2,
 504                                                   Double_t nmean2,
 505                                                   Double_t alpha ) const
 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
 519 //______________________________________________________________________________
 520 void 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;
 525 }
 526
 527 //______________________________________________________________________________
 528 void 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 }
 535
 536 //______________________________________________________________________________
 537 void 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 }
 544
 545 //______________________________________________________________________________
 546 void 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)
 552 {
 553   // set parameters used for calculating multiplicity see GetdNdEtaTwoNBD() for commments
 554   fdNdEtaTwoNBDParam[0]=k1;
 555   fdNdEtaTwoNBDParam[1]=nmean1;
 556   fdNdEtaTwoNBDParam[2]=k2;
 557   fdNdEtaTwoNBDParam[3]=nmean2;
 558   fdNdEtaTwoNBDParam[4]=alpha;
 559   fdNdEtaTwoNBDParam[5]=beta;
 560 }
 561
 562 //______________________________________________________________________________
 563 Double_t AliGlauberMC::GetdNdEta(Double_t nnp, Double_t x) const
 564 {
 565   //Get particle density per unit of rapidity
 566   //using two component model
 567   //Parameters: npp, x
 568   return nnp*((1.-x)*fNpart/2.+x*fNcoll);
 569 }
 570
 571 //______________________________________________________________________________
 572 Double_t AliGlauberMC::GetdNdEtaGBW( Double_t delta,
 573                                      Double_t lambda,
 574                                      Double_t  snn) const
 575 {
 576   //Get particle density per unit of rapidity
 577   //using the GBW model
 578   //Parameters: delta, lambda, snn
 579   return fNpart*0.47*TMath::Sqrt(TMath::Power(snn,lambda))
 580          * TMath::Power(fNpart,(1.-delta)/3./delta);
 581 }
 582
 583 //_______________________________________________________________________________
 584 Double_t AliGlauberMC::GetdNdEtaNBD ( Int_t k, Double_t nmean, Double_t beta) const
 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 }
 603
 604 //______________________________________________________________________________
 605 Double_t AliGlauberMC::GetdNdEtaTwoNBD ( Int_t k1,
 606                                          Double_t nmean1,
 607                                          Int_t k2,
 608                                          Double_t nmean2,
 609                                          Double_t alpha,
 610                                          Double_t beta ) const
 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;
 632 }
 633
 634 //______________________________________________________________________________
 635 Double_t AliGlauberMC::GetEccentricityPart() const
 636 {
 637   //get participant eccentricity of participants
 638   if (fNpart<2) return 0.0;
 639   return (TMath::Sqrt((fSy2-fSx2)*(fSy2-fSx2)+4*fSxy*fSxy)/(fSy2+fSx2));
 640 }
 641
 642 //______________________________________________________________________________
 643 Double_t AliGlauberMC::GetEccentricityPartColl() const
 644 {
 645   //get participant eccentricity of binary collisions
 646   if (fNcoll<2) return 0.0;
 647   return (TMath::Sqrt((fSy2Coll-fSx2Coll)*(fSy2Coll-fSx2Coll)+4*fSxyColl*fSxyColl)/(fSy2Coll+fSx2Coll));
 648 }
 649
 650 //______________________________________________________________________________
 651 Double_t AliGlauberMC::GetEccentricity() const
 652 {
 653   //get standard eccentricity of participants
 654   if (fNpart<2) return 0.0;
 655   return ((fSy2-fSx2)/(fSy2+fSx2));
 656 }
 657
 658 //______________________________________________________________________________
 659 Double_t AliGlauberMC::GetEccentricityColl() const
 660 {
 661   //get standard eccentricity of binary collisions
 662   if (fNcoll<2) return 0.0;
 663   return ((fSy2Coll-fSx2Coll)/(fSy2Coll+fSx2Coll));
 664 }
 665
 666 //______________________________________________________________________________
 667 Bool_t AliGlauberMC::NextEvent(Double_t bgen)
 668 {
 669   //Make a new event
 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     }
 678     if ( (succes=CalcEvent(bgen)) ) break; //ends if we have particparts
 679   }
 680   return succes;
 681 }
 682
 683 //______________________________________________________________________________
 684 void AliGlauberMC::Run(Int_t nevents)
 685 {
 686   //example run
 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,
 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");
 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     }
 706
 707     q++;
 708     Float_t v[27];
 709     v[0]  = GetNpart();
 710     v[1]  = GetNcoll();
 711     v[2]  = fBMC;
 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;
 722     v[13] = fMeanXA;
 723     v[14] = fMeanYA;
 724     v[15] = fMeanXB;
 725     v[16] = fMeanYB;
 726     v[17] = GetEccentricity();
 727     v[18] = GetEccentricityColl();
 728     v[19] = GetEccentricityPart();
 729     v[20] = GetEccentricityPartColl();
 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;
 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;
 758 }
 759
 760 //---------------------------------------------------------------------------------
 761 void AliGlauberMC::RunAndSaveNtuple( Int_t n,
 762                                      const Option_t *sysA,
 763                                      const Option_t *sysB,
 764                                      Double_t signn,
 765                                      Double_t mind,
 766                                                              Double_t r,
 767                                                              Double_t a,
 768                                      const char *fname)
 769 {
 770   //example run
 771   AliGlauberMC mcg(sysA,sysB,signn);
 772   mcg.SetMinDistance(mind);
 773   mcg.Setr(r);
 774   mcg.Seta(a);
 775   mcg.Run(n);
 776   TNtuple  *nt=mcg.GetNtuple();
 777   TFile out(fname,"recreate",fname,9);
 778   if(nt) nt->Write();
 779   printf("total cross section with a nucleon-nucleon cross section \t%f is \t%f",signn,mcg.GetTotXSect());
 780   out.Close();
 781 }
 782
 783 //---------------------------------------------------------------------------------
 784 void AliGlauberMC::RunAndSaveNucleons( Int_t n,
 785                                        const Option_t *sysA,
 786                                        Option_t *sysB,
 787                                        const Double_t signn,
 788                                        Double_t mind,
 789                                        Bool_t verbose,
 790                                        const char *fname)
 791 {
 792   //example run
 793   AliGlauberMC mcg(sysA,sysB,signn);
 794   mcg.SetMinDistance(mind);
 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
 803     while(!mcg.NextEvent()) {}
 804
 805     //access, save and (if wanted) print out nucleons
 806     TObjArray* nucleons=mcg.GetNucleons();
 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;
 814       cout<<"B = "<<mcg.GetB()<<"  Npart = "<<mcg.GetNpart()<<endl<<endl;
 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       {
 819         AliGlauberNucleon *nucl=(AliGlauberNucleon *)nucleons->UncheckedAt(iNucl);
 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);
 827       }
 828     }
 829   }
 830   if(out) delete out;
 831 }
 832
 833 //---------------------------------------------------------------------------------
 834 void AliGlauberMC::Reset()
 835 {
 836   //delete the ntuple
 837   delete fnt; fnt=NULL;
 838 }