PWG/Glauber/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 //  update:      You Zhou, Nikhef, yzhou@nikhef.nl
  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
  36 #include "AliGlauberNucleon.h"
  37 #include "AliGlauberNucleus.h"
  38 #include "AliGlauberMC.h"
  39
  40 using std::flush;
  41 ClassImp(AliGlauberMC)
  42
  43 //______________________________________________________________________________
  44 AliGlauberMC::AliGlauberMC(Option_t* NA, Option_t* NB, Double_t xsect) :
  45   TNamed(),
  46   fANucleus(NA),
  47   fBNucleus(NB),
  48   fXSect(xsect),
  49   fNucleonsA(0),
  50   fNucleonsB(0),
  51   fAN(0),
  52   fQAN(0),
  53   fBN(0),
  54   fQBN(0),
  55   fnt(0),
  56   fMeanX2(0),
  57   fMeanY2(0),
  58   fMeanXY(0),
  59   fMeanX2Parts(0),
  60   fMeanY2Parts(0),
  61   fMeanXYParts(0),
  62   fMeanXParts(0),
  63   fMeanYParts(0),
  64   fMeanOXParts(0),
  65   fMeanOYParts(0),
  66   fMeanXColl(0),
  67   fMeanYColl(0),
  68   fMeanOXColl(0),
  69   fMeanOYColl(0),
  70   fMeanX2Coll(0),
  71   fMeanY2Coll(0),
  72   fMeanXYColl(0),
  73   fMeanXCom(0),
  74   fMeanYCom(0),
  75   fMeanOXCom(0),
  76   fMeanOYCom(0),
  77   fMeanX2Com(0),
  78   fMeanY2Com(0),
  79   fMeanXYCom(0),
  80   fMeanXSystem(0),
  81   fMeanYSystem(0),
  82   fMeanXA(0),
  83   fMeanYA(0),
  84   fMeanXB(0),
  85   fMeanYB(0),
  86   fMeanOXA(0),
  87   fMeanOYA(0),
  88   fMeanOXB(0),
  89   fMeanOYB(0),
  90   fBMC(0),
  91   fEvents(0),
  92   fTotalEvents(0),
  93   fBMin(0.),
  94   fBMax(20.),
  95   fMultType(kNBDSV),
  96   fMaxNpartFound(0),
  97   fONpart(0),
  98   fONcoll(0),
  99   fONcom(0),
 100   fNpart(0),
 101   fNcoll(0),
 102   fNcollw(0),
 103   fNcom(0),
 104   fMeanr2(0),
 105   fMeanr3(0),
 106   fMeanr4(0),
 107   fMeanr5(0),
 108   fMeanr2Cos2Phi(0),
 109   fMeanr2Sin2Phi(0),
 110   fMeanr2Cos3Phi(0),
 111   fMeanr2Sin3Phi(0),
 112   fMeanr2Cos4Phi(0),
 113   fMeanr2Sin4Phi(0),
 114   fMeanr2Cos5Phi(0),
 115   fMeanr2Sin5Phi(0),
 116   fMeanr3Cos3Phi(0),
 117   fMeanr3Sin3Phi(0),
 118   fMeanr4Cos4Phi(0),
 119   fMeanr4Sin4Phi(0),
 120   fMeanr5Cos5Phi(0),
 121   fMeanr5Sin5Phi(0),
 122   fMeanr2Coll(0),
 123   fMeanr3Coll(0),
 124   fMeanr4Coll(0),
 125   fMeanr5Coll(0),
 126   fMeanr2Cos2PhiColl(0),
 127   fMeanr2Sin2PhiColl(0),
 128   fMeanr2Cos3PhiColl(0),
 129   fMeanr2Sin3PhiColl(0),
 130   fMeanr2Cos4PhiColl(0),
 131   fMeanr2Sin4PhiColl(0),
 132   fMeanr2Cos5PhiColl(0),
 133   fMeanr2Sin5PhiColl(0),
 134   fMeanr3Cos3PhiColl(0),
 135   fMeanr3Sin3PhiColl(0),
 136   fMeanr4Cos4PhiColl(0),
 137   fMeanr4Sin4PhiColl(0),
 138   fMeanr5Cos5PhiColl(0),
 139   fMeanr5Sin5PhiColl(0),
 140   fMeanr2Com(0),
 141   fMeanr3Com(0),
 142   fMeanr4Com(0),
 143   fMeanr5Com(0),
 144   fMeanr2Cos2PhiCom(0),
 145   fMeanr2Sin2PhiCom(0),
 146   fMeanr2Cos3PhiCom(0),
 147   fMeanr2Sin3PhiCom(0),
 148   fMeanr2Cos4PhiCom(0),
 149   fMeanr2Sin4PhiCom(0),
 150   fMeanr2Cos5PhiCom(0),
 151   fMeanr2Sin5PhiCom(0),
 152   fMeanr3Cos3PhiCom(0),
 153   fMeanr3Sin3PhiCom(0),
 154   fMeanr4Cos4PhiCom(0),
 155   fMeanr4Sin4PhiCom(0),
 156   fMeanr5Cos5PhiCom(0),
 157   fMeanr5Sin5PhiCom(0),
 158   fSx2Parts(0.),
 159   fSy2Parts(0.),
 160   fSxyParts(0.),
 161   fSx2Coll(0.),
 162   fSy2Coll(0.),
 163   fSxyColl(0.),
 164   fSx2Com(0.),
 165   fSy2Com(0.),
 166   fSxyCom(0.),
 167   fX(0.13),
 168   fNpp(8.),
 169   fDoPartProd(kFALSE),
 170   fBNN(0.),
 171   fDoFluc(kFALSE),
 172   fOmega(0),
 173   fSig0(0),
 174   fLambda(0),
 175   fSigFluc(0)
 176 {
 177   //ctor
 178   for (UInt_t i=0; i<(sizeof(fdNdEtaParam)/sizeof(fdNdEtaParam[0])); i++)
 179   {
 180     fdNdEtaParam[i]=0.0;
 181   }
 182
 183   SetName(Form("Glauber_%s_%s",fANucleus.GetName(),fBNucleus.GetName()));
 184   SetTitle(Form("Glauber %s+%s Version",fANucleus.GetName(),fBNucleus.GetName()));
 185 }
 186
 187 //______________________________________________________________________________
 188 AliGlauberMC::~AliGlauberMC()
 189 {
 190   //dtor
 191   delete fnt;
 192 }
 193
 194 //______________________________________________________________________________
 195 AliGlauberMC::AliGlauberMC(const AliGlauberMC& in):
 196   TNamed(in),
 197   fANucleus(in.fANucleus),
 198   fBNucleus(in.fBNucleus),
 199   fXSect(in.fXSect),
 200   fNucleonsA(in.fNucleonsA),
 201   fNucleonsB(in.fNucleonsB),
 202   fAN(in.fAN),
 203   fQAN(in.fQAN),
 204   fBN(in.fBN),
 205   fQBN(in.fQBN),
 206   fnt(in.fnt),
 207   fMeanX2(in.fMeanX2),
 208   fMeanY2(in.fMeanY2),
 209   fMeanXY(in.fMeanXY),
 210   fMeanX2Parts(in.fMeanX2Parts),
 211   fMeanY2Parts(in.fMeanY2Parts),
 212   fMeanXYParts(in.fMeanXYParts),
 213   fMeanXParts(in.fMeanXParts),
 214   fMeanYParts(in.fMeanYParts),
 215   fMeanOXParts(in.fMeanOXParts),
 216   fMeanOYParts(in.fMeanOYParts),
 217   fMeanXColl(in.fMeanXColl),
 218   fMeanYColl(in.fMeanYColl),
 219   fMeanOXColl(in.fMeanOXColl),
 220   fMeanOYColl(in.fMeanOYColl),
 221   fMeanX2Coll(in.fMeanX2Coll),
 222   fMeanY2Coll(in.fMeanY2Coll),
 223   fMeanXYColl(in.fMeanXYColl),
 224   fMeanXCom(in.fMeanXCom),
 225   fMeanYCom(in.fMeanYCom),
 226   fMeanOXCom(in.fMeanOXCom),
 227   fMeanOYCom(in.fMeanOYCom),
 228   fMeanX2Com(in.fMeanX2Com),
 229   fMeanY2Com(in.fMeanY2Com),
 230   fMeanXYCom(in.fMeanXYCom),
 231   fMeanXSystem(in.fMeanXSystem),
 232   fMeanYSystem(in.fMeanYSystem),
 233   fMeanXA(in.fMeanXA),
 234   fMeanYA(in.fMeanYA),
 235   fMeanXB(in.fMeanXB),
 236   fMeanYB(in.fMeanYB),
 237   fMeanOXA(in.fMeanOXA),
 238   fMeanOYA(in.fMeanOYA),
 239   fMeanOXB(in.fMeanOXB),
 240   fMeanOYB(in.fMeanOYB),
 241   fBMC(in.fBMC),
 242   fEvents(in.fEvents),
 243   fTotalEvents(in.fTotalEvents),
 244   fBMin(in.fBMin),
 245   fBMax(in.fBMax),
 246   fMultType(in.fMultType),
 247   fMaxNpartFound(in.fMaxNpartFound),
 248   fONpart(in.fONpart),
 249   fONcoll(in.fONcoll),
 250   fONcom(in.fONcom),
 251   fNpart(in.fNpart),
 252   fNcoll(in.fNcoll),
 253   fNcom(in.fNcom),
 254   fMeanr2(in.fMeanr2),
 255   fMeanr3(in.fMeanr3),
 256   fMeanr4(in.fMeanr4),
 257   fMeanr5(in.fMeanr5),
 258   fMeanr2Cos2Phi(in.fMeanr2Cos2Phi),
 259   fMeanr2Sin2Phi(in.fMeanr2Sin2Phi),
 260   fMeanr2Cos3Phi(in.fMeanr2Cos3Phi),
 261   fMeanr2Sin3Phi(in.fMeanr2Sin3Phi),
 262   fMeanr2Cos4Phi(in.fMeanr2Cos4Phi),
 263   fMeanr2Sin4Phi(in.fMeanr2Sin4Phi),
 264   fMeanr2Cos5Phi(in.fMeanr2Cos5Phi),
 265   fMeanr2Sin5Phi(in.fMeanr2Sin5Phi),
 266   fMeanr3Cos3Phi(in.fMeanr3Cos3Phi),
 267   fMeanr3Sin3Phi(in.fMeanr3Sin3Phi),
 268   fMeanr4Cos4Phi(in.fMeanr4Cos4Phi),
 269   fMeanr4Sin4Phi(in.fMeanr4Sin4Phi),
 270   fMeanr5Cos5Phi(in.fMeanr5Cos5Phi),
 271   fMeanr5Sin5Phi(in.fMeanr5Sin5Phi),
 272   fMeanr2Coll(in.fMeanr2Coll),
 273   fMeanr3Coll(in.fMeanr3Coll),
 274   fMeanr4Coll(in.fMeanr4Coll),
 275   fMeanr5Coll(in.fMeanr5Coll),
 276   fMeanr2Cos2PhiColl(in.fMeanr2Cos2PhiColl),
 277   fMeanr2Sin2PhiColl(in.fMeanr2Sin2PhiColl),
 278   fMeanr2Cos3PhiColl(in.fMeanr2Cos3PhiColl),
 279   fMeanr2Sin3PhiColl(in.fMeanr2Sin3PhiColl),
 280   fMeanr2Cos4PhiColl(in.fMeanr2Cos4PhiColl),
 281   fMeanr2Sin4PhiColl(in.fMeanr2Sin4PhiColl),
 282   fMeanr2Cos5PhiColl(in.fMeanr2Cos5PhiColl),
 283   fMeanr2Sin5PhiColl(in.fMeanr2Sin5PhiColl),
 284   fMeanr3Cos3PhiColl(in.fMeanr3Cos3PhiColl),
 285   fMeanr3Sin3PhiColl(in.fMeanr3Sin3PhiColl),
 286   fMeanr4Cos4PhiColl(in.fMeanr4Cos4PhiColl),
 287   fMeanr4Sin4PhiColl(in.fMeanr4Sin4PhiColl),
 288   fMeanr5Cos5PhiColl(in.fMeanr5Cos5PhiColl),
 289   fMeanr5Sin5PhiColl(in.fMeanr5Sin5PhiColl),
 290   fMeanr2Com(in.fMeanr2Com),
 291   fMeanr3Com(in.fMeanr3Com),
 292   fMeanr4Com(in.fMeanr4Com),
 293   fMeanr5Com(in.fMeanr5Com),
 294   fMeanr2Cos2PhiCom(in.fMeanr2Cos2PhiCom),
 295   fMeanr2Sin2PhiCom(in.fMeanr2Sin2PhiCom),
 296   fMeanr2Cos3PhiCom(in.fMeanr2Cos3PhiCom),
 297   fMeanr2Sin3PhiCom(in.fMeanr2Sin3PhiCom),
 298   fMeanr2Cos4PhiCom(in.fMeanr2Cos4PhiCom),
 299   fMeanr2Sin4PhiCom(in.fMeanr2Sin4PhiCom),
 300   fMeanr2Cos5PhiCom(in.fMeanr2Cos5PhiCom),
 301   fMeanr2Sin5PhiCom(in.fMeanr2Sin5PhiCom),
 302   fMeanr3Cos3PhiCom(in.fMeanr3Cos3PhiCom),
 303   fMeanr3Sin3PhiCom(in.fMeanr3Sin3PhiCom),
 304   fMeanr4Cos4PhiCom(in.fMeanr4Cos4PhiCom),
 305   fMeanr4Sin4PhiCom(in.fMeanr4Sin4PhiCom),
 306   fMeanr5Cos5PhiCom(in.fMeanr5Cos5PhiCom),
 307   fMeanr5Sin5PhiCom(in.fMeanr5Sin5PhiCom),
 308   fSx2Parts(in.fSx2Parts),
 309   fSy2Parts(in.fSy2Parts),
 310   fSxyParts(in.fSxyParts),
 311   fSx2Coll(in.fSx2Coll),
 312   fSy2Coll(in.fSy2Coll),
 313   fSxyColl(in.fSxyColl),
 314   fSx2Com(in.fSx2Com),
 315   fSy2Com(in.fSy2Com),
 316   fSxyCom(in.fSxyCom),
 317   fX(in.fX),
 318   fNpp(in.fNpp),
 319   fDoPartProd(in.fDoPartProd),
 320   fBNN(in.fBNN),
 321   fDoFluc(in.fDoFluc),
 322   fOmega(in.fOmega),
 323   fSig0(in.fSig0),
 324   fLambda(in.fLambda),
 325   fSigFluc(in.fSigFluc)
 326 {
 327   //copy ctor
 328   memcpy(fdNdEtaParam,in.fdNdEtaParam,sizeof(fdNdEtaParam));
 329 }
 330
 331 //______________________________________________________________________________
 332 AliGlauberMC& AliGlauberMC::operator=(const AliGlauberMC& in)
 333 {
 334   //assignment
 335   if (&in==this) return *this;
 336   fANucleus=in.fANucleus;
 337   fBNucleus=in.fBNucleus;
 338   fXSect=in.fXSect;
 339   fNucleonsA=in.fNucleonsA;
 340   fNucleonsB=in.fNucleonsB;
 341   fAN=in.fAN;
 342   fQAN=in.fQAN;
 343   fBN=in.fBN;
 344   fQBN=in.fQBN;
 345   fnt=in.fnt;
 346   fMeanX2=in.fMeanX2;
 347   fMeanY2=in.fMeanY2;
 348   fMeanXY=in.fMeanXY;
 349   fMeanr2=in.fMeanr2;
 350   fMeanr3=in.fMeanr3;
 351   fMeanr4=in.fMeanr4;
 352   fMeanr5=in.fMeanr5;
 353   fMeanXParts=in.fMeanXParts;
 354   fMeanYParts=in.fMeanYParts;
 355   fMeanOXParts=in.fMeanOXParts;
 356   fMeanOYParts=in.fMeanOYParts;
 357   fMeanXColl=in.fMeanXColl;
 358   fMeanYColl=in.fMeanYColl;
 359   fMeanOXColl=in.fMeanOXColl;
 360   fMeanOYColl=in.fMeanOYColl;
 361   fMeanX2Coll=in.fMeanX2Coll;
 362   fMeanY2Coll=in.fMeanY2Coll;
 363   fMeanXYColl=in.fMeanXYColl;
 364   fMeanr2Coll=in.fMeanr2Coll;
 365   fMeanr3Coll=in.fMeanr3Coll;
 366   fMeanr4Coll=in.fMeanr4Coll;
 367   fMeanr5Coll=in.fMeanr5Coll;
 368   fMeanXCom=in.fMeanXColl;
 369   fMeanYCom=in.fMeanYColl;
 370   fMeanOXCom=in.fMeanOXCom;
 371   fMeanOYCom=in.fMeanOYCom;
 372   fMeanX2Com=in.fMeanX2Com;
 373   fMeanY2Com=in.fMeanY2Com;
 374   fMeanXYCom=in.fMeanXYCom;
 375   fMeanr2Com=in.fMeanr2Com;
 376   fMeanr3Com=in.fMeanr3Com;
 377   fMeanr4Com=in.fMeanr4Com;
 378   fMeanr5Com=in.fMeanr5Com;
 379   fMeanXSystem=in.fMeanXSystem;
 380   fMeanYSystem=in.fMeanYSystem;
 381   fMeanXA=in.fMeanXA;
 382   fMeanYA=in.fMeanYA;
 383   fMeanXB=in.fMeanXB;
 384   fMeanYB=in.fMeanYB;
 385   fMeanOXA=in.fMeanOXA;
 386   fMeanOYA=in.fMeanOYA;
 387   fMeanOXB=in.fMeanOXB;
 388   fMeanOYB=in.fMeanOYB;
 389   fBMC=in.fBMC;
 390   fEvents=in.fEvents;
 391   fTotalEvents=in.fTotalEvents;
 392   fBMin=in.fBMin;
 393   fBMax=in.fBMax;
 394   fMultType=in.fMultType,
 395   memcpy(fdNdEtaParam,in.fdNdEtaParam,sizeof(fdNdEtaParam));
 396   fMaxNpartFound=in.fMaxNpartFound;
 397   fNpart=in.fNpart;
 398   fNcoll=in.fNcoll;
 399   fNcom=in.fNcom;
 400   fONpart=in.fONpart;
 401   fONcoll=in.fONcoll;
 402   fONcom=in.fONcom;
 403   fMeanr2Cos2Phi=in.fMeanr2Cos2Phi;
 404   fMeanr2Sin2Phi=in.fMeanr2Sin2Phi;
 405   fMeanr2Cos3Phi=in.fMeanr2Cos3Phi;
 406   fMeanr2Sin3Phi=in.fMeanr2Sin3Phi;
 407   fMeanr2Cos4Phi=in.fMeanr2Cos4Phi;
 408   fMeanr2Sin4Phi=in.fMeanr2Sin4Phi;
 409   fMeanr2Cos5Phi=in.fMeanr2Cos5Phi;
 410   fMeanr2Sin5Phi=in.fMeanr2Sin5Phi;
 411   fMeanr3Cos3Phi=in.fMeanr3Cos3Phi;
 412   fMeanr3Sin3Phi=in.fMeanr3Sin3Phi;
 413   fMeanr4Cos4Phi=in.fMeanr4Cos4Phi;
 414   fMeanr4Sin4Phi=in.fMeanr4Sin4Phi;
 415   fMeanr5Cos5Phi=in.fMeanr5Cos5Phi;
 416   fMeanr5Sin5Phi=in.fMeanr5Sin5Phi;
 417   fMeanr2Cos2PhiColl=in.fMeanr2Cos2PhiColl;
 418   fMeanr2Sin2PhiColl=in.fMeanr2Sin2PhiColl;
 419   fMeanr2Cos3PhiColl=in.fMeanr2Cos3PhiColl;
 420   fMeanr2Sin3PhiColl=in.fMeanr2Sin3PhiColl;
 421   fMeanr2Cos4PhiColl=in.fMeanr2Cos4PhiColl;
 422   fMeanr2Sin4PhiColl=in.fMeanr2Sin4PhiColl;
 423   fMeanr2Cos5PhiColl=in.fMeanr2Cos5PhiColl;
 424   fMeanr2Sin5PhiColl=in.fMeanr2Sin5PhiColl;
 425   fMeanr3Cos3PhiColl=in.fMeanr3Cos3PhiColl;
 426   fMeanr3Sin3PhiColl=in.fMeanr3Sin3PhiColl;
 427   fMeanr4Cos4PhiColl=in.fMeanr4Cos4PhiColl;
 428   fMeanr4Sin4PhiColl=in.fMeanr4Sin4PhiColl;
 429   fMeanr5Cos5PhiColl=in.fMeanr5Cos5PhiColl;
 430   fMeanr5Sin5PhiColl=in.fMeanr5Sin5PhiColl;
 431   fMeanr2Cos2PhiCom=in.fMeanr2Cos2PhiCom;
 432   fMeanr2Sin2PhiCom=in.fMeanr2Sin2PhiCom;
 433   fMeanr2Cos3PhiCom=in.fMeanr2Cos3PhiCom;
 434   fMeanr2Sin3PhiCom=in.fMeanr2Sin3PhiCom;
 435   fMeanr2Cos4PhiCom=in.fMeanr2Cos4PhiCom;
 436   fMeanr2Sin4PhiCom=in.fMeanr2Sin4PhiCom;
 437   fMeanr2Cos5PhiCom=in.fMeanr2Cos5PhiCom;
 438   fMeanr2Sin5PhiCom=in.fMeanr2Sin5PhiCom;
 439   fMeanr3Cos3PhiCom=in.fMeanr3Cos3PhiCom;
 440   fMeanr3Sin3PhiCom=in.fMeanr3Sin3PhiCom;
 441   fMeanr4Cos4PhiCom=in.fMeanr4Cos4PhiCom;
 442   fMeanr4Sin4PhiCom=in.fMeanr4Sin4PhiCom;
 443   fMeanr5Cos5PhiCom=in.fMeanr5Cos5PhiCom;
 444   fMeanr5Sin5PhiCom=in.fMeanr5Sin5PhiCom;
 445   fSx2Parts=in.fSx2Parts;
 446   fSy2Parts=in.fSy2Parts;
 447   fSxyParts=in.fSxyParts;
 448   fSx2Coll=in.fSx2Coll;
 449   fSy2Coll=in.fSy2Coll;
 450   fSxyColl=in.fSxyColl;
 451   fSx2Com=in.fSx2Com;
 452   fSy2Com=in.fSy2Com;
 453   fSxyCom=in.fSxyCom;
 454   fX=in.fX;
 455   fNpp=in.fNpp;
 456   return *this;
 457 }
 458
 459 //______________________________________________________________________________
 460 Bool_t AliGlauberMC::CalcEvent(Double_t bgen)
 461 {
 462   // prepare event
 463
 464   if (fDoFluc) {
 465     if (!fSigFluc) {
 466       fSigFluc = new TF1("fSigFluc","[0]*x/[3]/(x/[3]+[1])*exp(-((x/[1]/[3]-1)/[2])^2)",0,250);
 467       fSigFluc->SetParameters(1,fSig0,fOmega,fLambda);
 468       cout << "Setting fluc: " << fSig0 << " " << fOmega << " " << fLambda << endl;
 469     }
 470   }
 471
 472   fANucleus.ThrowNucleons(-bgen/2.);
 473   fNucleonsA = fANucleus.GetNucleons();
 474   fAN = fANucleus.GetN();
 475   fQAN = fAN * 3;
 476   //fAN = 3 * fANucleus.GetN(); // for Pb, Number of quark = 3*208;
 477   for (Int_t i = 0; i<fAN; i++)
 478   {
 479     AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->UncheckedAt(i));
 480     nucleonA->SetInNucleusA();
 481     nucleonA->SetSigNN(fXSect);
 482     if (fDoFluc)
 483       nucleonA->SetSigNN(fSigFluc->GetRandom());
 484   }
 485   fBNucleus.ThrowNucleons(bgen/2.);
 486   fNucleonsB = fBNucleus.GetNucleons();
 487   //fBN = 3 * fBNucleus.GetN(); // Number of quark = number of nucleus*3;
 488   fBN = fBNucleus.GetN();
 489   fQBN = fBN * 3;
 490   for (Int_t i = 0; i<fBN; i++)
 491   {
 492     AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->UncheckedAt(i));
 493     nucleonB->SetInNucleusB();
 494     nucleonB->SetSigNN(fXSect);
 495     if (fDoFluc)
 496       nucleonB->SetSigNN(fSigFluc->GetRandom());
 497   }
 498
 499   if (fDoFluc) {
 500     if (!fSigFluc) {
 501       fSigFluc = new TF1("fSigFluc","[0]*x/[3]/(x/[3]+[1])*exp(-((x/[1]/[3]-1)/[2])^2)",0,250);
 502       fSigFluc->SetParameters(1,fSig0,fOmega,fLambda);
 503       cout << "Setting fluc: " << fSig0 << " " << fOmega << " " << fLambda << endl;
 504     }
 505     fXSect = fSigFluc->GetRandom();
 506   }
 507   // "ball" diameter = distance at which two balls interact
 508   Double_t d2 = (Double_t)fXSect/(TMath::Pi()*10); // in fm^2
 509
 510   Double_t bNN   = 0;
 511   Double_t Nco   = 0;
 512   Double_t Ncohc = 0; // hard core
 513
 514   // for each of the A nucleons in nucleus B
 515   for (Int_t i = 0; i<fBN; i++)
 516   {
 517     AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->UncheckedAt(i));
 518     for (Int_t j = 0 ; j < fAN ; j++)
 519     {
 520       AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->UncheckedAt(j));
 521       Double_t dx = nucleonB->GetX()-nucleonA->GetX();
 522       Double_t dy = nucleonB->GetY()-nucleonA->GetY();
 523       Double_t dij = dx*dx+dy*dy;
 524       if (fDoFluc) {
 525         //fXSect = nucleonA->GetSigNN();
 526         //fXSect = (nucleonA->GetSigNN()+nucleonB->GetSigNN())/2.;
 527         fXSect = TMath::Max(nucleonA->GetSigNN(),nucleonB->GetSigNN());
 528         d2 = (Double_t)fXSect/(TMath::Pi()*10); // in fm^2
 529       }
 530       if (dij < d2)
 531       {
 532         bNN += dij;
 533         ++Nco;
 534         nucleonB->Collide();
 535         nucleonA->Collide();
 536         if (dij<d2/4)
 537           ++Ncohc;
 538       }
 539     }
 540   }
 541
 542   if (Nco>0) {
 543     fNcollw = Ncohc;
 544     fBNN = bNN/Nco;
 545   } else {
 546     fNcollw = 0;
 547     fBNN    = 0.;
 548   }
 549
 550   if (Nco>0)
 551     fBNN = bNN/Nco;
 552   else
 553     fBNN = 0.;
 554   return CalcResults(bgen);
 555 }
 556
 557 //______________________________________________________________________________
 558 Bool_t AliGlauberMC::CalcResults(Double_t bgen)
 559 {
 560   // calc results for the given event
 561   //return true if we have participants
 562
 563   fNpart=0;
 564   fNcoll=0;
 565   fNcom=0;
 566   fONpart=0;
 567   fONcoll=0;
 568   fONcom=0;
 569   fMeanX2=0.;
 570   fMeanY2=0.;
 571   fMeanXY=0.;
 572   fMeanXParts=0.;
 573   fMeanYParts=0.;
 574   fMeanOXParts=0.;
 575   fMeanOYParts=0.;
 576   fMeanXColl=0.;
 577   fMeanYColl=0.;
 578   fMeanOXColl=0.;
 579   fMeanOYColl=0.;
 580   fMeanX2Coll=0.;
 581   fMeanY2Coll=0.;
 582   fMeanXYColl=0.;
 583   fMeanXCom=0.;
 584   fMeanYCom=0.;
 585   fMeanOXCom=0.;
 586   fMeanOYCom=0.;
 587   fMeanX2Com=0.;
 588   fMeanY2Com=0.;
 589   fMeanXYCom=0.;
 590   fMeanXSystem=0.;
 591   fMeanYSystem=0.;
 592   fMeanXA=0.;
 593   fMeanYA=0.;
 594   fMeanXB=0.;
 595   fMeanYB=0.;
 596   fMeanOXA=0.;
 597   fMeanOYA=0.;
 598   fMeanOXB=0.;
 599   fMeanOYB=0.;
 600   fMeanr2=0.;
 601   fMeanr3=0.;
 602   fMeanr4=0.;
 603   fMeanr5=0.;
 604   fMeanr2Cos2Phi=0.;
 605   fMeanr2Sin2Phi=0.;
 606   fMeanr2Cos3Phi=0.;
 607   fMeanr2Sin3Phi=0.;
 608   fMeanr2Cos4Phi=0.;
 609   fMeanr2Sin4Phi=0.;
 610   fMeanr2Cos5Phi=0.;
 611   fMeanr2Sin5Phi=0.;
 612   fMeanr3Cos3Phi=0.;
 613   fMeanr3Sin3Phi=0.;
 614   fMeanr4Cos4Phi=0.;
 615   fMeanr4Sin4Phi=0.;
 616   fMeanr5Cos5Phi=0.;
 617   fMeanr5Sin5Phi=0.;
 618   fMeanr2Coll=0.;
 619   fMeanr3Coll=0.;
 620   fMeanr4Coll=0.;
 621   fMeanr5Coll=0.;
 622   fMeanr2Cos2PhiColl=0.;
 623   fMeanr2Sin2PhiColl=0.;
 624   fMeanr2Cos3PhiColl=0.;
 625   fMeanr2Sin3PhiColl=0.;
 626   fMeanr2Cos4PhiColl=0.;
 627   fMeanr2Sin4PhiColl=0.;
 628   fMeanr2Cos5PhiColl=0.;
 629   fMeanr2Sin5PhiColl=0.;
 630   fMeanr3Cos3PhiColl=0.;
 631   fMeanr3Sin3PhiColl=0.;
 632   fMeanr4Cos4PhiColl=0.;
 633   fMeanr4Sin4PhiColl=0.;
 634   fMeanr5Cos5PhiColl=0.;
 635   fMeanr5Sin5PhiColl=0.;
 636   fMeanr2Com=0.;
 637   fMeanr3Com=0.;
 638   fMeanr4Com=0.;
 639   fMeanr5Com=0.;
 640   fMeanr2Cos2PhiCom=0.;
 641   fMeanr2Sin2PhiCom=0.;
 642   fMeanr2Cos3PhiCom=0.;
 643   fMeanr2Sin3PhiCom=0.;
 644   fMeanr2Cos4PhiCom=0.;
 645   fMeanr2Sin4PhiCom=0.;
 646   fMeanr2Cos5PhiCom=0.;
 647   fMeanr2Sin5PhiCom=0.;
 648   fMeanr3Cos3PhiCom=0.;
 649   fMeanr3Sin3PhiCom=0.;
 650   fMeanr4Cos4PhiCom=0.;
 651   fMeanr4Sin4PhiCom=0.;
 652   fMeanr5Cos5PhiCom=0.;
 653   fMeanr5Sin5PhiCom=0.;
 654
 655   for (Int_t i = 0; i<fAN; i++)
 656   {
 657     AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->UncheckedAt(i));
 658     Double_t oXA = nucleonA->GetX();
 659     Double_t oYA = nucleonA->GetY();
 660     //fMeanOXSystem  += oXA;
 661     //fMeanOYSystem  += oYA;
 662     fMeanOXA  += oXA;
 663     fMeanOYA  += oYA;
 664
 665     if(nucleonA->IsWounded())
 666     {
 667       fONpart++;
 668       fMeanOXParts  += oXA;
 669       fMeanOYParts  += oYA;
 670
 671       fONcom += (1-0.150);
 672       fMeanOXCom  += oXA*(1-0.150);
 673       fMeanOYCom += oXA*(1-0.150);
 674     }
 675   }
 676
 677   for (Int_t i = 0; i<fBN; i++)
 678   {
 679     AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->UncheckedAt(i));
 680     Double_t oXB=nucleonB->GetX();
 681     Double_t oYB=nucleonB->GetY();
 682
 683     if(nucleonB->IsWounded())
 684     {
 685       Int_t oNcoll = nucleonB->GetNColl();
 686       fONpart++;
 687       fMeanOXParts  += oXB;
 688       fMeanOXColl  += oXB*oNcoll;
 689       fMeanOXCom  += oXB*((1-0.150)+0.150*oNcoll);
 690       fMeanOYParts  += oYB;
 691       fMeanOYColl += oYB*oNcoll;
 692       fMeanOYCom += oYB*((1-0.150)+0.150*oNcoll);
 693       fONcoll += oNcoll;
 694       fONcom += ((1-0.150)+0.150*oNcoll);
 695     }
 696   }
 697
 698   if (fONpart>0)
 699   {
 700     fMeanOXParts /= fONpart;
 701     fMeanOYParts /= fONpart;
 702   }
 703   else
 704   {
 705     fMeanOXParts = 0;
 706     fMeanOYParts = 0;
 707   }
 708
 709   if (fONcoll>0)
 710   {
 711     fMeanOXColl /= fONcoll;
 712     fMeanOYColl /= fONcoll;
 713   }
 714   else
 715   {
 716     fMeanOXColl = 0;
 717     fMeanOYColl = 0;
 718   }
 719
 720  if (fONcom>0)
 721   {
 722     fMeanOXCom /= fONcom;
 723     fMeanOYCom /= fONcom;
 724   }
 725   else
 726   {
 727     fMeanOXCom = 0;
 728     fMeanOYCom = 0;
 729   }
 730
 731   //////////////////////////////////////////////////////////////////
 732   for (Int_t i = 0; i<fAN; i++)
 733   {
 734     AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->UncheckedAt(i));
 735     Double_t xAA = nucleonA->GetX(); // X
 736     Double_t yAA = nucleonA->GetY(); // Y
 737     Double_t xAPart = xAA - fMeanOXParts; // X'
 738     Double_t yAPart = yAA - fMeanOYParts; // Y'
 739     Double_t r2APart = xAPart *xAPart+yAPart*yAPart;     // r'^2
 740     Double_t rAPart = TMath::Sqrt(r2APart);  // r'
 741     Double_t r3APart = r2APart*rAPart;
 742     Double_t r4APart = r3APart*rAPart;
 743     Double_t r5APart = r4APart*rAPart;
 744     Double_t phiAPart = TMath::ATan2(yAPart,xAPart);
 745     Double_t sin2PhiAPart = TMath::Sin(2.*phiAPart);
 746     Double_t cos2PhiAPart = TMath::Cos(2.*phiAPart);
 747     Double_t sin3PhiAPart = TMath::Sin(3.*phiAPart);
 748     Double_t cos3PhiAPart = TMath::Cos(3.*phiAPart);
 749     Double_t sin4PhiAPart = TMath::Sin(4.*phiAPart);
 750     Double_t cos4PhiAPart = TMath::Cos(4.*phiAPart);
 751     Double_t sin5PhiAPart = TMath::Sin(5.*phiAPart);
 752     Double_t cos5PhiAPart = TMath::Cos(5.*phiAPart);
 753
 754     // for combine
 755     Double_t xACom = xAA - fMeanOXCom; // X'-Combine
 756     Double_t yACom = yAA - fMeanOYCom; // Y'-C
 757     Double_t r2ACom = xACom *xACom+yACom*yACom;     // r'^2-C
 758     Double_t rACom = TMath::Sqrt(r2ACom);  // r'-C
 759     Double_t r3ACom = r2ACom*rACom;
 760     Double_t r4ACom = r3ACom*rACom;
 761     Double_t r5ACom = r4ACom*rACom;
 762     Double_t phiACom = TMath::ATan2(yACom,xACom);
 763     Double_t sin2PhiACom = TMath::Sin(2.*phiACom);
 764     Double_t cos2PhiACom = TMath::Cos(2.*phiACom);
 765     Double_t sin3PhiACom = TMath::Sin(3.*phiACom);
 766     Double_t cos3PhiACom = TMath::Cos(3.*phiACom);
 767     Double_t sin4PhiACom = TMath::Sin(4.*phiACom);
 768     Double_t cos4PhiACom = TMath::Cos(4.*phiACom);
 769     Double_t sin5PhiACom = TMath::Sin(5.*phiACom);
 770     Double_t cos5PhiACom = TMath::Cos(5.*phiACom);
 771
 772     fMeanXSystem  += xAA;
 773     fMeanYSystem  += yAA;
 774     fMeanXA  += xAA;
 775     fMeanYA  += yAA;
 776     fMeanX2 += xAA * xAA;
 777     fMeanY2 += yAA * yAA;
 778     fMeanXY += xAA * yAA;
 779
 780     if(nucleonA->IsWounded())
 781      {
 782        //Wounded
 783       fNpart++;
 784       fMeanXParts  += xAPart;
 785       fMeanYParts  += yAPart;
 786       fMeanX2Parts += xAPart * xAPart;
 787       fMeanY2Parts += yAPart * yAPart;
 788       fMeanXYParts += xAPart * yAPart;
 789       fMeanr2 += r2APart;
 790       fMeanr3 += r3APart;
 791       fMeanr4 += r4APart;
 792       fMeanr5 += r5APart;
 793       fMeanr2Cos2Phi += r2APart*cos2PhiAPart;
 794       fMeanr2Sin2Phi += r2APart*sin2PhiAPart;
 795       fMeanr2Cos3Phi += r2APart*cos3PhiAPart;
 796       fMeanr2Sin3Phi += r2APart*sin3PhiAPart;
 797       fMeanr2Cos4Phi += r2APart*cos4PhiAPart;
 798       fMeanr2Sin4Phi += r2APart*sin4PhiAPart;
 799       fMeanr2Cos5Phi += r2APart*cos5PhiAPart;
 800       fMeanr2Sin5Phi += r2APart*sin5PhiAPart;
 801       fMeanr3Cos3Phi += r3APart*cos3PhiAPart;
 802       fMeanr3Sin3Phi += r3APart*sin3PhiAPart;
 803       fMeanr4Cos4Phi += r4APart*cos4PhiAPart;
 804       fMeanr4Sin4Phi += r4APart*sin4PhiAPart;
 805       fMeanr5Cos5Phi += r5APart*cos5PhiAPart;
 806       fMeanr5Sin5Phi += r5APart*sin5PhiAPart;
 807
 808       // Combined
 809       fMeanXCom  += xACom*(1-0.150);
 810       fMeanYCom += yACom*(1-0.150);
 811       fMeanX2Com += xACom*xACom*(1-0.150);
 812       fMeanY2Com += yACom*yACom*(1-0.150);
 813       fMeanXYCom += xACom*yACom*(1-0.150);
 814       fNcom += (1-0.150);
 815       fMeanr2Com += r2ACom*(1-0.150);
 816       fMeanr3Com += r3ACom*(1-0.150);
 817       fMeanr4Com += r4ACom*(1-0.150);
 818       fMeanr5Com += r5ACom*(1-0.150);
 819       fMeanr2Cos2PhiCom += r2ACom*cos2PhiACom*(1-0.150);
 820       fMeanr2Sin2PhiCom += r2ACom*sin2PhiACom*(1-0.150);
 821       fMeanr2Cos3PhiCom += r2ACom*cos3PhiACom*(1-0.150);
 822       fMeanr2Sin3PhiCom += r2ACom*sin3PhiACom*(1-0.150);
 823       fMeanr2Cos4PhiCom += r2ACom*cos4PhiACom*(1-0.150);
 824       fMeanr2Sin4PhiCom += r2ACom*sin4PhiACom*(1-0.150);
 825       fMeanr2Cos5PhiCom += r2ACom*cos5PhiACom*(1-0.150);
 826       fMeanr2Sin5PhiCom += r2ACom*sin5PhiACom*(1-0.150);
 827       fMeanr3Cos3PhiCom += r3ACom*cos3PhiACom*(1-0.150);
 828       fMeanr3Sin3PhiCom += r3ACom*sin3PhiACom*(1-0.150);
 829       fMeanr4Cos4PhiCom += r4ACom*cos4PhiACom*(1-0.150);
 830       fMeanr4Sin4PhiCom += r4ACom*sin4PhiACom*(1-0.150);
 831       fMeanr5Cos5PhiCom += r5ACom*cos5PhiACom*(1-0.150);
 832       fMeanr5Sin5PhiCom += r5ACom*sin5PhiACom*(1-0.150);
 833
 834     }
 835   }
 836
 837   for (Int_t i = 0; i<fBN; i++)
 838     {
 839       AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->UncheckedAt(i));
 840       Double_t xBB = nucleonB->GetX();
 841       Double_t yBB = nucleonB->GetY();
 842       // for Wounded
 843       Double_t xBPart = xBB - fMeanOXParts; // X'
 844       Double_t yBPart = yBB - fMeanOYParts; // Y'
 845       Double_t r2BPart = xBPart*xBPart+yBPart*yBPart;     // r'^2
 846       Double_t rBPart = TMath::Sqrt(r2BPart);  // r'
 847       Double_t r3BPart = r2BPart*rBPart;
 848       Double_t r4BPart = r3BPart*rBPart;
 849       Double_t r5BPart = r4BPart*rBPart;
 850       Double_t phiBPart = TMath::ATan2(yBPart,xBPart);
 851       Double_t sin2PhiBPart = TMath::Sin(2.*phiBPart);
 852       Double_t cos2PhiBPart = TMath::Cos(2.*phiBPart);
 853       Double_t sin3PhiBPart = TMath::Sin(3.*phiBPart);
 854       Double_t cos3PhiBPart = TMath::Cos(3.*phiBPart);
 855       Double_t sin4PhiBPart = TMath::Sin(4.*phiBPart);
 856       Double_t cos4PhiBPart = TMath::Cos(4.*phiBPart);
 857       Double_t sin5PhiBPart = TMath::Sin(5.*phiBPart);
 858       Double_t cos5PhiBPart = TMath::Cos(5.*phiBPart);
 859       // for Binary
 860       Double_t xBColl = xBB - fMeanOXColl; // X'-Binary
 861       Double_t yBColl = yBB - fMeanOYColl; // Y'-B
 862       Double_t r2BColl = xBColl*xBColl+yBColl*yBColl;     // r'^2-B
 863       Double_t rBColl = TMath::Sqrt(r2BColl);  // r'-B
 864       Double_t r3BColl = r2BColl*rBColl;
 865       Double_t r4BColl = r3BColl*rBColl;
 866       Double_t r5BColl = r4BColl*rBColl;
 867       Double_t phiBColl = TMath::ATan2(yBColl,xBColl);
 868       Double_t sin2PhiBColl = TMath::Sin(2.*phiBColl);
 869       Double_t cos2PhiBColl = TMath::Cos(2.*phiBColl);
 870       Double_t sin3PhiBColl = TMath::Sin(3.*phiBColl);
 871       Double_t cos3PhiBColl = TMath::Cos(3.*phiBColl);
 872       Double_t sin4PhiBColl = TMath::Sin(4.*phiBColl);
 873       Double_t cos4PhiBColl = TMath::Cos(4.*phiBColl);
 874       Double_t sin5PhiBColl = TMath::Sin(5.*phiBColl);
 875       Double_t cos5PhiBColl = TMath::Cos(5.*phiBColl);
 876       // for combine
 877       Double_t xBCom = xBB - fMeanOXCom; // X'-Combine
 878       Double_t yBCom = yBB - fMeanOYCom; // Y'-C
 879       Double_t r2BCom = xBCom *xBCom+yBCom*yBCom;     // r'^2-C
 880       Double_t rBCom = TMath::Sqrt(r2BCom);  // r'-C
 881       Double_t r3BCom = r2BCom*rBCom;
 882       Double_t r4BCom = r3BCom*rBCom;
 883       Double_t r5BCom = r4BCom*rBCom;
 884       Double_t phiBCom = TMath::ATan2(yBCom,xBCom);
 885       Double_t sin2PhiBCom = TMath::Sin(2.*phiBCom);
 886       Double_t cos2PhiBCom = TMath::Cos(2.*phiBCom);
 887       Double_t sin3PhiBCom = TMath::Sin(3.*phiBCom);
 888       Double_t cos3PhiBCom = TMath::Cos(3.*phiBCom);
 889       Double_t sin4PhiBCom = TMath::Sin(4.*phiBCom);
 890       Double_t cos4PhiBCom = TMath::Cos(4.*phiBCom);
 891       Double_t sin5PhiBCom = TMath::Sin(5.*phiBCom);
 892       Double_t cos5PhiBCom = TMath::Cos(5.*phiBCom);
 893
 894       fMeanXSystem  += xBB;
 895       fMeanYSystem  += yBB;
 896       fMeanXB  += xBB;
 897       fMeanYB  += yBB;
 898       fMeanX2 += xBB*xBB;
 899       fMeanY2 += yBB*yBB;
 900       fMeanXY += xBB*yBB;
 901
 902       if(nucleonB->IsWounded())
 903         {
 904           Int_t ncoll = nucleonB->GetNColl();
 905           fNpart++;
 906           fMeanXParts  += xBPart;
 907           fMeanXColl  += xBColl*ncoll;
 908           fMeanXCom  += xBCom*((1-0.150)+0.150*ncoll);
 909           fMeanYParts  += yBPart;
 910           fMeanYColl += yBColl*ncoll;
 911           fMeanYCom += yBCom*((1-0.150)+0.150*ncoll);
 912           fMeanX2Parts += xBPart * xBPart;
 913           fMeanX2Coll += xBColl*xBColl*ncoll;
 914           fMeanX2Com += xBCom*xBCom*((1-0.150)+0.150*ncoll);
 915           fMeanY2Parts += yBPart * yBPart;
 916           fMeanY2Coll += yBColl*yBColl*ncoll;
 917           fMeanY2Com += yBCom*yBCom*((1-0.150)+0.150*ncoll);
 918           fMeanXYParts += xBPart * yBPart;
 919           fMeanXYColl += xBColl*yBColl*ncoll;
 920           fMeanXYCom += xBCom*yBCom*((1-0.150)+0.150*ncoll);
 921           fNcoll += ncoll;
 922           fNcom += ((1-0.150)+0.150*ncoll);
 923           fMeanr2 += r2BPart;
 924           fMeanr3 += r3BPart;
 925           fMeanr4 += r4BPart;
 926           fMeanr5 += r5BPart;
 927           fMeanr2Cos2Phi += r2BPart*cos2PhiBPart;
 928           fMeanr2Sin2Phi += r2BPart*sin2PhiBPart;
 929           fMeanr2Cos3Phi += r2BPart*cos3PhiBPart;
 930           fMeanr2Sin3Phi += r2BPart*sin3PhiBPart;
 931           fMeanr2Cos4Phi += r2BPart*cos4PhiBPart;
 932           fMeanr2Sin4Phi += r2BPart*sin4PhiBPart;
 933           fMeanr2Cos5Phi += r2BPart*cos5PhiBPart;
 934           fMeanr2Sin5Phi += r2BPart*sin5PhiBPart;
 935           fMeanr3Cos3Phi += r3BPart*cos3PhiBPart;
 936           fMeanr3Sin3Phi += r3BPart*sin3PhiBPart;
 937           fMeanr4Cos4Phi += r4BPart*cos4PhiBPart;
 938           fMeanr4Sin4Phi += r4BPart*sin4PhiBPart;
 939           fMeanr5Cos5Phi += r5BPart*cos5PhiBPart;
 940           fMeanr5Sin5Phi += r5BPart*sin5PhiBPart;
 941           fMeanr2Coll += r2BColl*ncoll;
 942           fMeanr3Coll += r3BColl*ncoll;
 943           fMeanr4Coll += r4BColl*ncoll;
 944           fMeanr5Coll += r5BColl*ncoll;
 945           fMeanr2Cos2PhiColl += r2BColl*cos2PhiBColl*ncoll;
 946           fMeanr2Sin2PhiColl += r2BColl*sin2PhiBColl*ncoll;
 947           fMeanr2Cos3PhiColl += r2BColl*cos3PhiBColl*ncoll;
 948           fMeanr2Sin3PhiColl += r2BColl*sin3PhiBColl*ncoll;
 949           fMeanr2Cos4PhiColl += r2BColl*cos4PhiBColl*ncoll;
 950           fMeanr2Sin4PhiColl += r2BColl*sin4PhiBColl*ncoll;
 951           fMeanr2Cos5PhiColl += r2BColl*cos5PhiBColl*ncoll;
 952           fMeanr2Sin5PhiColl += r2BColl*sin5PhiBColl*ncoll;
 953           fMeanr3Cos3PhiColl += r3BColl*cos3PhiBColl*ncoll;
 954           fMeanr3Sin3PhiColl += r3BColl*sin3PhiBColl*ncoll;
 955           fMeanr4Cos4PhiColl += r4BColl*cos4PhiBColl*ncoll;
 956           fMeanr4Sin4PhiColl += r4BColl*sin4PhiBColl*ncoll;
 957           fMeanr5Cos5PhiColl += r5BColl*cos5PhiBColl*ncoll;
 958           fMeanr5Sin5PhiColl += r5BColl*sin5PhiBColl*ncoll;
 959           fMeanr2Com += r2BCom*((1-0.150)+0.150*ncoll);
 960           fMeanr3Com += r3BCom*((1-0.150)+0.150*ncoll);
 961           fMeanr4Com += r4BCom*((1-0.150)+0.150*ncoll);
 962           fMeanr5Com += r5BCom*((1-0.150)+0.150*ncoll);
 963           fMeanr2Cos2PhiCom += r2BCom*cos2PhiBCom*((1-0.150)+0.150*ncoll);
 964           fMeanr2Sin2PhiCom += r2BCom*sin2PhiBCom*((1-0.150)+0.150*ncoll);
 965           fMeanr2Cos3PhiCom += r2BCom*cos3PhiBCom*((1-0.150)+0.150*ncoll);
 966           fMeanr2Sin3PhiCom += r2BCom*sin3PhiBCom*((1-0.150)+0.150*ncoll);
 967           fMeanr2Cos4PhiCom += r2BCom*cos4PhiBCom*((1-0.150)+0.150*ncoll);
 968           fMeanr2Sin4PhiCom += r2BCom*sin4PhiBCom*((1-0.150)+0.150*ncoll);
 969           fMeanr2Cos5PhiCom += r2BCom*cos5PhiBCom*((1-0.150)+0.150*ncoll);
 970           fMeanr2Sin5PhiCom += r2BCom*sin5PhiBCom*((1-0.150)+0.150*ncoll);
 971           fMeanr3Cos3PhiCom += r3BCom*cos3PhiBCom*((1-0.150)+0.150*ncoll);
 972           fMeanr3Sin3PhiCom += r3BCom*sin3PhiBCom*((1-0.150)+0.150*ncoll);
 973           fMeanr4Cos4PhiCom += r4BCom*cos4PhiBCom*((1-0.150)+0.150*ncoll);
 974           fMeanr4Sin4PhiCom += r4BCom*sin4PhiBCom*((1-0.150)+0.150*ncoll);
 975           fMeanr5Cos5PhiCom += r5BCom*cos5PhiBCom*((1-0.150)+0.150*ncoll);
 976           fMeanr5Sin5PhiCom += r5BCom*sin5PhiBCom*((1-0.150)+0.150*ncoll);
 977         }
 978     }
 979
 980   if (fNpart>0)
 981     {
 982       fMeanXParts /= fNpart;
 983       fMeanYParts /= fNpart;
 984       fMeanX2Parts /= fNpart;
 985       fMeanY2Parts /= fNpart;
 986       fMeanXYParts /= fNpart;
 987       fMeanr2 /= fNpart;
 988       fMeanr3 /= fNpart;
 989       fMeanr4 /= fNpart;
 990       fMeanr5 /= fNpart;
 991       fMeanr2Cos2Phi /= fNpart;
 992       fMeanr2Sin2Phi /= fNpart;
 993       fMeanr2Cos3Phi /= fNpart;
 994       fMeanr2Sin3Phi /= fNpart;
 995       fMeanr2Cos4Phi /= fNpart;
 996       fMeanr2Sin4Phi /= fNpart;
 997       fMeanr2Cos5Phi /= fNpart;
 998       fMeanr2Sin5Phi /= fNpart;
 999       fMeanr3Cos3Phi /= fNpart;
1000       fMeanr3Sin3Phi /= fNpart;
1001       fMeanr4Cos4Phi /= fNpart;
1002       fMeanr4Sin4Phi /= fNpart;
1003       fMeanr5Cos5Phi /= fNpart;
1004       fMeanr5Sin5Phi /= fNpart;
1005     }
1006
1007   else
1008     {
1009       fMeanXParts = 0;
1010       fMeanYParts = 0;
1011       fMeanX2Parts = 0;
1012       fMeanY2Parts = 0;
1013       fMeanXYParts = 0;
1014       fMeanr2 = 0;
1015       fMeanr3 = 0;
1016       fMeanr4 = 0;
1017       fMeanr5 = 0;
1018       fMeanr2Cos2Phi = 0;
1019       fMeanr2Sin2Phi = 0;
1020       fMeanr2Cos3Phi = 0;
1021       fMeanr2Sin3Phi = 0;
1022       fMeanr2Cos4Phi = 0;
1023       fMeanr2Sin4Phi = 0;
1024       fMeanr2Cos5Phi = 0;
1025       fMeanr2Sin5Phi = 0;
1026       fMeanr3Cos3Phi = 0;
1027       fMeanr3Sin3Phi = 0;
1028       fMeanr4Cos4Phi = 0;
1029       fMeanr4Sin4Phi = 0;
1030       fMeanr5Cos5Phi = 0;
1031       fMeanr5Sin5Phi = 0;
1032     }
1033
1034   if (fNcoll>0)
1035     {
1036       fMeanXColl /= fNcoll;
1037       fMeanYColl /= fNcoll;
1038       fMeanX2Coll /= fNcoll;
1039       fMeanY2Coll /= fNcoll;
1040       fMeanXYColl /= fNcoll;
1041       fMeanr2Coll /= fNcoll;
1042       fMeanr3Coll /= fNcoll;
1043       fMeanr4Coll /= fNcoll;
1044       fMeanr5Coll /= fNcoll;
1045       fMeanr2Cos2PhiColl /= fNcoll;
1046       fMeanr2Sin2PhiColl /= fNcoll;
1047       fMeanr2Cos3PhiColl /= fNcoll;
1048       fMeanr2Sin3PhiColl /= fNcoll;
1049       fMeanr2Cos4PhiColl /= fNcoll;
1050       fMeanr2Sin4PhiColl /= fNcoll;
1051       fMeanr2Cos5PhiColl /= fNcoll;
1052       fMeanr2Sin5PhiColl /= fNcoll;
1053       fMeanr3Cos3PhiColl /= fNcoll;
1054       fMeanr3Sin3PhiColl /= fNcoll;
1055       fMeanr4Cos4PhiColl /= fNcoll;
1056       fMeanr4Sin4PhiColl /= fNcoll;
1057       fMeanr5Cos5PhiColl /= fNcoll;
1058       fMeanr5Sin5PhiColl /= fNcoll;
1059   }
1060
1061   else
1062   {
1063     fMeanXColl = 0;
1064     fMeanYColl = 0;
1065     fMeanX2Coll = 0;
1066     fMeanY2Coll = 0;
1067     fMeanXYColl = 0;
1068     fMeanr2Cos2PhiColl =0;
1069     fMeanr2Sin2PhiColl =0;
1070     fMeanr2Cos3PhiColl =0;
1071     fMeanr2Sin3PhiColl =0;
1072     fMeanr2Cos4PhiColl =0;
1073     fMeanr2Sin4PhiColl =0;
1074     fMeanr2Cos5PhiColl =0;
1075     fMeanr2Sin5PhiColl =0;
1076     fMeanr3Cos3PhiColl =0;
1077     fMeanr3Sin3PhiColl =0;
1078     fMeanr4Cos4PhiColl =0;
1079     fMeanr4Sin4PhiColl =0;
1080     fMeanr5Cos5PhiColl =0;
1081     fMeanr5Sin5PhiColl =0;
1082   }
1083
1084  if (fNcom>0)
1085     {
1086       fMeanXCom /= fNcom;
1087       fMeanYCom /= fNcom;
1088       fMeanX2Com /= fNcom;
1089       fMeanY2Com /= fNcom;
1090       fMeanXYCom /= fNcom;
1091       fMeanr2Com /= fNcom;
1092       fMeanr3Com /= fNcom;
1093       fMeanr4Com /= fNcom;
1094       fMeanr5Com /= fNcom;
1095       fMeanr2Cos2PhiCom /= fNcom;
1096       fMeanr2Sin2PhiCom /= fNcom;
1097       fMeanr2Cos3PhiCom /= fNcom;
1098       fMeanr2Sin3PhiCom /= fNcom;
1099       fMeanr2Cos4PhiCom /= fNcom;
1100       fMeanr2Sin4PhiCom /= fNcom;
1101       fMeanr2Cos5PhiCom /= fNcom;
1102       fMeanr2Sin5PhiCom /= fNcom;
1103       fMeanr3Cos3PhiCom /= fNcom;
1104       fMeanr3Sin3PhiCom /= fNcom;
1105       fMeanr4Cos4PhiCom /= fNcom;
1106       fMeanr4Sin4PhiCom /= fNcom;
1107       fMeanr5Cos5PhiCom /= fNcom;
1108       fMeanr5Sin5PhiCom /= fNcom;
1109   }
1110
1111  else
1112   {
1113     fMeanXCom = 0;
1114     fMeanYCom = 0;
1115     fMeanX2Com = 0;
1116     fMeanY2Com = 0;
1117     fMeanXYCom = 0;
1118     fMeanr2Cos2PhiCom =0;
1119     fMeanr2Sin2PhiCom =0;
1120     fMeanr2Cos3PhiCom =0;
1121     fMeanr2Sin3PhiCom =0;
1122     fMeanr2Cos4PhiCom =0;
1123     fMeanr2Sin4PhiCom =0;
1124     fMeanr2Cos5PhiCom =0;
1125     fMeanr2Sin5PhiCom =0;
1126     fMeanr3Cos3PhiCom =0;
1127     fMeanr3Sin3PhiCom =0;
1128     fMeanr4Cos4PhiCom =0;
1129     fMeanr4Sin4PhiCom =0;
1130     fMeanr5Cos5PhiCom =0;
1131     fMeanr5Sin5PhiCom =0;
1132   }
1133
1134   if(fAN+fBN>0)
1135   {
1136     fMeanXSystem /= (fAN + fBN);
1137     fMeanYSystem /= (fAN + fBN);
1138     fMeanX2 /= (fAN + fBN);
1139     fMeanY2 /= (fAN + fBN);
1140     fMeanXY /= (fAN + fBN);
1141   }
1142   else
1143   {
1144     fMeanXSystem = 0;
1145     fMeanYSystem = 0;
1146     fMeanX2 = 0;
1147     fMeanY2 = 0;
1148     fMeanXY = 0;
1149   }
1150   if(fAN>0)
1151   {
1152     fMeanXA /= fAN;
1153     fMeanYA /= fAN;
1154   }
1155   else
1156   {
1157     fMeanXA = 0;
1158     fMeanYA = 0;
1159   }
1160
1161   if(fBN>0)
1162   {
1163     fMeanXB /= fBN;
1164     fMeanYB /= fBN;
1165   }
1166   else
1167   {
1168     fMeanXB = 0;
1169     fMeanYB = 0;
1170   }
1171
1172
1173
1174
1175   //////////////////////////////////////////////////////////////////
1176   fSx2Parts=fMeanX2Parts-(fMeanXParts*fMeanXParts);
1177   fSy2Parts=fMeanY2Parts-(fMeanYParts*fMeanYParts);
1178   fSxyParts=fMeanXYParts-fMeanXParts*fMeanYParts;
1179   fSx2Coll=fMeanX2Coll-(fMeanXColl*fMeanXColl);
1180   fSy2Coll=fMeanY2Coll-(fMeanYColl*fMeanYColl);
1181   fSxyColl=fMeanXYColl-fMeanXColl*fMeanYColl;
1182   fSx2Com=fMeanX2Com-(fMeanXCom*fMeanXCom);
1183   fSy2Com=fMeanY2Com-(fMeanYCom*fMeanYCom);
1184   fSxyCom=fMeanXYCom-fMeanXCom*fMeanYCom;
1185   fBMC = bgen;
1186   fTotalEvents++;
1187   if (fNpart>0) fEvents++;
1188   if (fNpart==0) return kFALSE;
1189   if (fNpart > fMaxNpartFound) fMaxNpartFound = fNpart;
1190
1191   return kTRUE;
1192 }
1193
1194 //______________________________________________________________________________
1195 void AliGlauberMC::Draw(Option_t* /*option*/)
1196 {
1197   //draw method
1198   fANucleus.Draw(fXSect, 2);
1199   fBNucleus.Draw(fXSect, 4);
1200
1201   TEllipse e;
1202   e.SetFillColor(0);
1203   e.SetLineColor(1);
1204   e.SetLineStyle(2);
1205   e.SetLineWidth(1);
1206   e.DrawEllipse(GetB()/2,0,fBNucleus.GetR(),fBNucleus.GetR(),0,360,0);
1207   e.DrawEllipse(-GetB()/2,0,fANucleus.GetR(),fANucleus.GetR(),0,360,0);
1208 }
1209
1210 //______________________________________________________________________________
1211 Double_t AliGlauberMC::GetTotXSect() const
1212 {
1213   //total xsection
1214   return (1.*fEvents/fTotalEvents)*TMath::Pi()*fBMax*fBMax/100;
1215 }
1216
1217 //______________________________________________________________________________
1218 Double_t AliGlauberMC::GetTotXSectErr() const
1219 {
1220   //total xsection error
1221   return GetTotXSect()/TMath::Sqrt((Double_t)fEvents) *
1222          TMath::Sqrt(Double_t(1.-fEvents/fTotalEvents));
1223 }
1224
1225 //______________________________________________________________________________
1226 TObjArray *AliGlauberMC::GetNucleons()
1227 {
1228   //get array of nucleons
1229   if(!fNucleonsA || !fNucleonsB) return 0;
1230   fNucleonsA->SetOwner(0);
1231   fNucleonsB->SetOwner(0);
1232   TObjArray *allnucleons=new TObjArray(fAN+fBN);
1233   allnucleons->SetOwner();
1234   for (Int_t i = 0; i<fAN; i++)
1235   {
1236     allnucleons->Add(fNucleonsA->UncheckedAt(i));
1237   }
1238   for (Int_t i = 0; i<fBN; i++)
1239   {
1240     allnucleons->Add(fNucleonsB->UncheckedAt(i));
1241   }
1242   return allnucleons;
1243 }
1244
1245 //______________________________________________________________________________
1246 Double_t AliGlauberMC::NegativeBinomialDistribution(Int_t x, Int_t k, Double_t nmean)
1247 {
1248   //produce a number distributed acc. neg.bin.dist
1249   if(k<=0)
1250   {
1251     cout << "Error, zero or negative K" << endl;
1252     return 0;
1253   }
1254   return (TMath::Binomial(x+k-1,x))
1255          *TMath::Power(((nmean/Double_t(k))/(1+nmean/Double_t(k))),Double_t(x))
1256          *(1/(TMath::Power((1+nmean/Double_t(k)),Double_t(k))));
1257 }
1258
1259 //______________________________________________________________________________
1260 Int_t AliGlauberMC::NegativeBinomialRandom(Int_t k, Double_t nmean) const
1261 {
1262   //return random integer from a Negative Binomial Distribution
1263   static const Int_t fMaxPlot = 1000;
1264   Double_t array[fMaxPlot];
1265   array[0] = NegativeBinomialDistribution(0,k,nmean);
1266   for (Int_t i=1; i<fMaxPlot; i++)
1267   {
1268     array[i] = NegativeBinomialDistribution(i,k,nmean) + array[i-1];
1269   }
1270   Double_t r = gRandom->Uniform(0,1);
1271   return TMath::BinarySearch(fMaxPlot,array,r)+2;
1272
1273 }
1274
1275 //______________________________________________________________________________
1276 Int_t AliGlauberMC::NegativeBinomialRandomSV(Double_t k, Double_t nbar) const
1277 {
1278   // negative binomial distribution generator, S. Voloshin, 09-May-2007
1279   Double_t sum=0.;
1280   Int_t i=0;
1281   Double_t ran=gRandom->Rndm();
1282   Double_t trm=1./pow(1.+nbar/k,k);
1283   if (trm==0.)
1284   {
1285     cout<<"NBD overflow"<<"  nbar="<<nbar<<"   k="<<k<<endl;
1286     return -1;
1287   }
1288   for(i=0; i<2000 && sum<ran ; i++)
1289   {
1290     sum += trm;
1291     trm *= (k+i)/(i+1.)*(nbar/(nbar+k));
1292   }
1293   return i-1;
1294 }
1295
1296 //______________________________________________________________________________
1297 Int_t AliGlauberMC::DoubleNegativeBinomialRandom( Int_t k,
1298     Double_t nmean,
1299     Int_t k2,
1300     Double_t nmean2,
1301     Double_t alpha ) const
1302 {
1303   //return random integer from a Double Negative Binomial Distribution
1304   static const Int_t fMaxPlot = 1000;
1305   Double_t array[fMaxPlot];
1306   array[0] = alpha*NegativeBinomialDistribution(0,k,nmean)+(1-alpha)*NegativeBinomialDistribution(0,k2,nmean2);
1307   for (Int_t i=1; i<fMaxPlot; i++)
1308   {
1309     array[i] = alpha*NegativeBinomialDistribution(i,k,nmean)+(1-alpha)*NegativeBinomialDistribution(i,k2,nmean2) + array[i-1];
1310   }
1311   Double_t r = gRandom->Uniform(0,1);
1312   return TMath::BinarySearch(fMaxPlot,array,r)+2;
1313 }
1314
1315 //______________________________________________________________________________
1316 Double_t AliGlauberMC::GetdNdEta() const
1317 {
1318   switch (fMultType)
1319   {
1320     case kSimple:
1321       return GetdNdEtaSimple(fdNdEtaParam);
1322     case kNBD:
1323       return GetdNdEtaNBD(fdNdEtaParam);
1324     case kNBDSV:
1325       return GetdNdEtaNBDSV(fdNdEtaParam);
1326     case kTwoNBD:
1327       return GetdNdEtaTwoNBD(fdNdEtaParam);
1328     case kGBW:
1329       return GetdNdEtaGBW(fdNdEtaParam);
1330     default:
1331       return 0.0;
1332   }
1333 }
1334
1335 //______________________________________________________________________________
1336 Double_t AliGlauberMC::GetdNdEtaSimple(const Double_t* p) const
1337 {
1338   //Get particle density per unit of rapidity
1339   //using two component model
1340   //Parameters: npp, x
1341   Double_t nnp = p[0]; //=8.0
1342   Double_t x = p[1]; //=0.13
1343
1344   return nnp*((1.-x)*fNpart/2.+x*fNcoll);
1345 }
1346
1347 //______________________________________________________________________________
1348 Double_t AliGlauberMC::GetdNdEtaGBW( const Double_t* p ) const
1349 {
1350   //Get particle density per unit of rapidity
1351   //using the GBW model
1352   //Parameters: delta, lambda, snn
1353   Double_t delta = p[0]; //=0.79
1354   Double_t lambda = p[1]; //=0.288
1355   Double_t  snn = p[2]; //=30.25
1356
1357   return fNpart*0.47*TMath::Sqrt(TMath::Power(snn,lambda))
1358          * TMath::Power(fNpart,(1.-delta)/3./delta);
1359 }
1360
1361 //_______________________________________________________________________________
1362 Double_t AliGlauberMC::GetdNdEtaNBDSV ( const Double_t* p ) const
1363 {
1364   //Get particle density per unit of rapidity (from Sergei)
1365   Double_t npp = p[0];          //=2.49
1366   Double_t ratioSgm2Mu = p[1];  //=1.7
1367   Double_t xhard = p[2];        //=0.13
1368
1369   double knb=npp/(ratioSgm2Mu-1.);  // sgm^2/mu = 1+ mu/k
1370   double scale = (1.-xhard)*fNpart/2.+xhard*fNcoll; // effectively get number of pp collisions
1371   float nch1=-99.;
1372   if (knb*scale <1000.)
1373   {
1374     nch1=(float) NegativeBinomialRandomSV( npp*scale, knb*scale );
1375   }
1376   else
1377   {
1378     nch1=(float) NegativeBinomialRandomSV( npp*scale/2., knb*scale/2. ) +
1379          (float) NegativeBinomialRandomSV( npp*scale/2., knb*scale/2. );
1380   }
1381   return nch1;
1382 }
1383
1384 //_______________________________________________________________________________
1385 Double_t AliGlauberMC::GetdNdEtaNBD ( const Double_t* p ) const
1386 {
1387   //Get particle density per unit of rapidity
1388   //using a aandomized number from a negative binomial distrubution
1389   //Parameters:   k  = related to distribition width=3
1390   //              nmean = mean of distribution=4
1391   //              beta = set contribution of participants / binary collisions to multiplicity=0.13
1392   Int_t k = TMath::Nint(p[0]);
1393   Double_t nmean = p[1];
1394   Double_t beta = p[2];
1395
1396   Double_t mulnp=0.0;
1397   for(Int_t i = 0; i<fNpart; i++)
1398   {
1399     mulnp+=NegativeBinomialRandom(k,nmean);
1400   }
1401   Double_t mulnb=0.0;
1402   for(Int_t i = 0; i<fNcoll; i++)
1403   {
1404     mulnb+=NegativeBinomialRandom(k,nmean);
1405   }
1406   return (1-beta)*mulnp/2+beta*mulnb;
1407 }
1408
1409 //______________________________________________________________________________
1410 Double_t AliGlauberMC::GetdNdEtaTwoNBD ( const Double_t* p ) const
1411 {
1412   //Get particle density per unit of rapidity
1413   //using random numbers from two negative binomial distributions
1414   //Parameters:   k1 = related to distribition width of distribution 1=3
1415   //              nmean1 = mean of distribution 1=4
1416   //              k2 = related to distribition width of distribution 2=2
1417   //              nmean2 = mean of distribution 2=11
1418   //              alpha = set contributions of distrubitin 1 / distribution 2=0.4
1419   //              beta = set contribution of participants / binary collisions to multiplicity =0.13
1420   Int_t k1 = TMath::Nint(p[0]);
1421   Double_t nmean1 = p[1];
1422   Int_t k2 = TMath::Nint(p[2]);
1423   Double_t nmean2 = p[3];
1424   Double_t alpha = p[4];
1425   Double_t beta = p[6];
1426
1427   Double_t mulnp=0.0;
1428   for(Int_t i = 0; i<fNpart; i++)
1429   {
1430     mulnp+=DoubleNegativeBinomialRandom(k1,nmean1,k2,nmean2,alpha);
1431   }
1432   Double_t mulnb=0.0;
1433   for(Int_t i = 0; i<fNcoll; i++)
1434   {
1435     mulnb+=DoubleNegativeBinomialRandom(k1,nmean1,k2,nmean2,alpha);
1436   }
1437   Double_t mul=(1-beta)*mulnp/2+beta*mulnb;
1438   return mul;
1439 }
1440
1441 //______________________________________________________________________________
1442 Double_t AliGlauberMC::GetEccentricityPart() const
1443 {
1444   //get participant eccentricity of participants
1445   if (fNpart<2) return 0.0;
1446   return (TMath::Sqrt((fSy2Parts-fSx2Parts)*(fSy2Parts-fSx2Parts)+4*fSxyParts*fSxyParts)/(fSy2Parts+fSx2Parts));
1447 }
1448
1449 //_____________________________________________________________________________
1450 Double_t AliGlauberMC::GetEccentricityPartColl() const
1451 {
1452   //get participant eccentricity of binary collisions
1453   if (fNcoll<0) return 0.0;
1454   if (fSy2Coll==0.0) return 0.0;
1455   return (TMath::Sqrt((fSy2Coll-fSx2Coll)*(fSy2Coll-fSx2Coll)+4*fSxyColl*fSxyColl)/(fSy2Coll+fSx2Coll));
1456 }
1457
1458 //_____________________________________________________________________________
1459 Double_t AliGlauberMC::GetEccentricityPartCom() const
1460 {
1461   //get participant eccentricity of binary collisions
1462   if (fNcom<0) return 0.0;
1463   return (TMath::Sqrt((fSy2Com-fSx2Com)*(fSy2Com-fSx2Com)+4*fSxyCom*fSxyCom)/(fSy2Com+fSx2Com));
1464 }
1465 //______________________________________________________________________________
1466 Double_t AliGlauberMC::GetEccentricity() const
1467 {
1468   //get standard eccentricity of participants
1469   if (fNpart<2) return 0.0;
1470   return ((fSy2Parts-fSx2Parts)/(fSy2Parts+fSx2Parts));
1471 }
1472
1473 //______________________________________________________________________________
1474 Double_t AliGlauberMC::GetEccentricityColl() const
1475 {
1476   //get standard eccentricity of binary collisions
1477   if (fNcoll<0) return 0.0;
1478   if (fSy2Coll==0.0) return 0.0;
1479   return ((fSy2Coll-fSx2Coll)/(fSy2Coll+fSx2Coll));
1480 }
1481
1482 //______________________________________________________________________________
1483 Double_t AliGlauberMC::GetEccentricityCom() const
1484 {
1485   //get standard eccentricity of combined weight
1486   if (fNcom<0) return 0.0;
1487   return ((fSy2Com-fSx2Com)/(fSy2Com+fSx2Com));
1488   //return (((fMeanY2Com-(fMeanYCom*fMeanYCom)) - (fMeanX2Com-(fMeanXCom*fMeanXCom))) / ((fMeanY2Com-(fMeanYCom*fMeanYCom)) + (fMeanX2Com-(fMeanXCom*fMeanXCom))) );
1489 }
1490
1491 //______________________________________________________________________________
1492 Double_t AliGlauberMC::GetStoa() const
1493 {
1494   //get standard Transverse Overlap Area
1495   if (fNpart<2) return 0.0;
1496   return ( TMath::Pi()*(TMath::Sqrt(fSx2Parts))*(TMath::Sqrt(fSy2Parts)));
1497 }
1498
1499 //______________________________________________________________________________
1500 Bool_t AliGlauberMC::NextEvent(Double_t bgen)
1501 {
1502   //Make a new event
1503   Int_t nAttempts = 10; // set indices, max attempts and max comparisons
1504   Bool_t succes = kFALSE;
1505   for(Int_t j=0; j<nAttempts; j++)
1506   {
1507     if(bgen<0||!succes) //get impactparameter
1508     {
1509       bgen = TMath::Sqrt((fBMax*fBMax-fBMin*fBMin)*gRandom->Rndm()+fBMin*fBMin);
1510     }
1511     if ( (succes=CalcEvent(bgen)) ) break; //ends if we have particparts
1512   }
1513   return succes;
1514 }
1515 //______________________________________________________________________________
1516 Double_t AliGlauberMC::GetEpsilon2Part() const
1517 {
1518   //get participant eccentricity of participants
1519   if (fNpart<2) return 0.0;
1520   return (TMath::Sqrt(fMeanr2Cos2Phi*fMeanr2Cos2Phi+fMeanr2Sin2Phi*fMeanr2Sin2Phi)/fMeanr2);
1521 }
1522 //______________________________________________________________________________
1523 Double_t AliGlauberMC::GetEpsilon3Part() const
1524 {
1525   //get participant eccentricity of participants
1526   if (fNpart<2) return 0.0;
1527   return (TMath::Sqrt(fMeanr2Cos3Phi*fMeanr2Cos3Phi+fMeanr2Sin3Phi*fMeanr2Sin3Phi)/fMeanr2);
1528   //return (TMath::Sqrt(fMeanr3Cos3Phi*fMeanr3Cos3Phi+fMeanr3Sin3Phi*fMeanr3Sin3Phi)/fMeanr3);
1529 }
1530 //______________________________________________________________________________
1531 Double_t AliGlauberMC::GetEpsilon4Part() const
1532 {
1533   //get participant eccentricity of participants
1534   if (fNpart<2) return 0.0;
1535   return (TMath::Sqrt(fMeanr2Cos4Phi*fMeanr2Cos4Phi+fMeanr2Sin4Phi*fMeanr2Sin4Phi)/fMeanr2);
1536   //return (TMath::Sqrt(fMeanr4Cos4Phi*fMeanr4Cos4Phi+fMeanr4Sin4Phi*fMeanr4Sin4Phi)/fMeanr4);
1537 }
1538
1539 //______________________________________________________________________________
1540 Double_t AliGlauberMC::GetEpsilon5Part() const
1541 {
1542   //get participant eccentricity of participants
1543   if (fNpart<2) return 0.0;
1544   return (TMath::Sqrt(fMeanr2Cos5Phi*fMeanr2Cos5Phi+fMeanr2Sin5Phi*fMeanr2Sin5Phi)/fMeanr2);
1545   //return (TMath::Sqrt(fMeanr5Cos5Phi*fMeanr5Cos5Phi+fMeanr5Sin5Phi*fMeanr5Sin5Phi)/fMeanr5);
1546 }
1547
1548 //______________________________________________________________________________
1549 Double_t AliGlauberMC::GetEpsilon2Coll() const
1550 {
1551   //get epsilon2 of binary collisions
1552   if (fNcoll<0) return 0.0;
1553   if (fMeanr2Coll==0.0) return 0.0;
1554   return (TMath::Sqrt(fMeanr2Cos2PhiColl*fMeanr2Cos2PhiColl+fMeanr2Sin2PhiColl*fMeanr2Sin2PhiColl)/fMeanr2Coll);
1555 }
1556
1557 //______________________________________________________________________________
1558 Double_t AliGlauberMC::GetEpsilon3Coll() const
1559 {
1560   //get epsilon3 of binary collisions
1561   if (fNcoll<0) return 0.0;
1562   if (fMeanr2Coll==0.0) return 0.0;
1563   return (TMath::Sqrt(fMeanr2Cos3PhiColl*fMeanr2Cos3PhiColl+fMeanr2Sin3PhiColl*fMeanr2Sin3PhiColl)/fMeanr2Coll);
1564   //return (TMath::Sqrt(fMeanr3Cos3PhiColl*fMeanr3Cos3PhiColl+fMeanr3Sin3PhiColl*fMeanr3Sin3PhiColl)/fMeanr3Coll);
1565 }
1566
1567 //______________________________________________________________________________
1568 Double_t AliGlauberMC::GetEpsilon4Coll() const
1569 {
1570   //get epsilon4 of binary collisions
1571   if (fNcoll<0) return 0.0;
1572   if (fMeanr2Coll==0.0) return 0.0;
1573   return (TMath::Sqrt(fMeanr2Cos4PhiColl*fMeanr2Cos4PhiColl+fMeanr2Sin4PhiColl*fMeanr2Sin4PhiColl)/fMeanr2Coll);
1574   //return (TMath::Sqrt(fMeanr4Cos4PhiColl*fMeanr4Cos4PhiColl+fMeanr4Sin4PhiColl*fMeanr4Sin4PhiColl)/fMeanr4Coll);
1575 }
1576
1577 //______________________________________________________________________________
1578 Double_t AliGlauberMC::GetEpsilon5Coll() const
1579 {
1580   //get epsilon5 of binary collisions
1581   if (fNcoll<0) return 0.0;
1582   if (fMeanr2Coll==0.0) return 0.0;
1583   return (TMath::Sqrt(fMeanr2Cos5PhiColl*fMeanr2Cos5PhiColl+fMeanr2Sin5PhiColl*fMeanr2Sin5PhiColl)/fMeanr2Coll);
1584   //return (TMath::Sqrt(fMeanr5Cos5PhiColl*fMeanr5Cos5PhiColl+fMeanr5Sin5PhiColl*fMeanr5Sin5PhiColl)/fMeanr5Coll);
1585 }
1586
1587 //______________________________________________________________________________
1588 Double_t AliGlauberMC::GetEpsilon2Com() const
1589 {
1590   //get epsilon2 of binary collisions
1591   if (fNcom<0) return 0.0;
1592   return (TMath::Sqrt(fMeanr2Cos2PhiCom*fMeanr2Cos2PhiCom+fMeanr2Sin2PhiCom*fMeanr2Sin2PhiCom)/fMeanr2Com);
1593 }
1594
1595 //______________________________________________________________________________
1596 Double_t AliGlauberMC::GetEpsilon3Com() const
1597 {
1598   //get epsilon3 of binary collisions
1599   if (fNcom<0) return 0.0;
1600   return (TMath::Sqrt(fMeanr2Cos3PhiCom*fMeanr2Cos3PhiCom+fMeanr2Sin3PhiCom*fMeanr2Sin3PhiCom)/fMeanr2Com);
1601   //return (TMath::Sqrt(fMeanr3Cos3PhiCom*fMeanr3Cos3PhiCom+fMeanr3Sin3PhiCom*fMeanr3Sin3PhiCom)/fMeanr3Com);
1602 }
1603
1604 //______________________________________________________________________________
1605 Double_t AliGlauberMC::GetEpsilon4Com() const
1606 {
1607   //get epsilon4 of binary collisions
1608   if (fNcom<0) return 0.0;
1609   return (TMath::Sqrt(fMeanr2Cos4PhiCom*fMeanr2Cos4PhiCom+fMeanr2Sin4PhiCom*fMeanr2Sin4PhiCom)/fMeanr2Com);
1610   //return (TMath::Sqrt(fMeanr4Cos4PhiCom*fMeanr4Cos4PhiCom+fMeanr4Sin4PhiCom*fMeanr4Sin4PhiCom)/fMeanr4Com);
1611 }
1612
1613 //______________________________________________________________________________
1614 Double_t AliGlauberMC::GetEpsilon5Com() const
1615 {
1616   //get epsilon5 of binary collisions
1617   if (fNcom<0) return 0.0;
1618   return (TMath::Sqrt(fMeanr2Cos5PhiCom*fMeanr2Cos5PhiCom+fMeanr2Sin5PhiCom*fMeanr2Sin5PhiCom)/fMeanr2Com);
1619   //return (TMath::Sqrt(fMeanr5Cos5PhiCom*fMeanr5Cos5PhiCom+fMeanr5Sin5PhiCom*fMeanr5Sin5PhiCom)/fMeanr5Com);
1620 }
1621
1622 //______________________________________________________________________________
1623 Double_t AliGlauberMC::GetPsi2() const
1624 {
1625   return ((TMath::ATan2(fMeanr2Sin2Phi,fMeanr2Cos2Phi)+TMath::Pi())/2);
1626 }
1627
1628 //______________________________________________________________________________
1629 Double_t AliGlauberMC::GetPsi3() const
1630 {
1631   return ((TMath::ATan2(fMeanr2Sin3Phi,fMeanr2Cos3Phi)+TMath::Pi())/3);
1632   //return ((TMath::ATan2(fMeanr3Sin3Phi,fMeanr3Cos3Phi)+TMath::Pi())/3);
1633 }
1634
1635 //______________________________________________________________________________
1636 Double_t AliGlauberMC::GetPsi4() const
1637 {
1638   return ((TMath::ATan2(fMeanr2Sin4Phi,fMeanr2Cos4Phi)+TMath::Pi())/4);
1639   //return ((TMath::ATan2(fMeanr4Sin4Phi,fMeanr4Cos4Phi)+TMath::Pi())/4);
1640 }
1641
1642 //______________________________________________________________________________
1643 Double_t AliGlauberMC::GetPsi5() const
1644 {
1645   return ((TMath::ATan2(fMeanr2Sin5Phi,fMeanr2Cos5Phi)+TMath::Pi())/5);
1646   //return ((TMath::ATan2(fMeanr5Sin5Phi,fMeanr5Cos5Phi)+TMath::Pi())/5);
1647 }
1648
1649 /*_____________________________________________________________________________
1650 Double_t AliGlauberMC::GetE43Part() const
1651 {
1652   //get participant eccentricity of participants
1653   if (fNpart<2) return 0.0;
1654   //return ((TMath::Sqrt(fMeanr2Cos2Phi*fMeanr2Cos2Phi+fMeanr2Sin2Phi*fMeanr2Sin2Phi)/fMeanr2) * (TMath::Sqrt(fMeanr2Cos2Phi*fMeanr2Cos2Phi+fMeanr2Sin2Phi*fMeanr2Sin2Phi)/fMeanr2) * (TMath::Sqrt(fMeanr2Cos4Phi*fMeanr2Cos4Phi+fMeanr2Sin4Phi*fMeanr2Sin4Phi)/fMeanr2) * TMath::Cos (((TMath::ATan2(fMeanr2Sin2Phi,fMeanr2Cos2Phi)+TMath::Pi())/2) - ((TMath::ATan2(fMeanr2Sin4Phi,fMeanr2Cos4Phi)+TMath::Pi())/4)));
1655   return ((TMath::Sqrt(fMeanr2Cos2Phi*fMeanr2Cos2Phi+fMeanr2Sin2Phi*fMeanr2Sin2Phi)/fMeanr2) * (TMath::Sqrt(fMeanr2Cos2Phi*fMeanr2Cos2Phi+fMeanr2Sin2Phi*fMeanr2Sin2Phi)/fMeanr2) * (TMath::Sqrt(fMeanr4Cos4Phi*fMeanr4Cos4Phi+fMeanr4Sin4Phi*fMeanr4Sin4Phi)/fMeanr4) * TMath::Cos(((TMath::ATan2(fMeanr2Sin2Phi,fMeanr2Cos2Phi)+TMath::Pi())/2) - ((TMath::ATan2(fMeanr4Sin4Phi,fMeanr4Cos4Phi)+TMath::Pi())/4)));
1656 }
1657
1658
1659 //______________________________________________________________________________
1660 Double_t AliGlauberMC::GetE43Coll() const
1661 {
1662   //get epsilon5 of binary collisions
1663   if (fNcoll<2) return 0.0;
1664   //return ((TMath::Sqrt(fMeanr2Cos2PhiColl*fMeanr2Cos2PhiColl+fMeanr2Sin2PhiColl*fMeanr2Sin2PhiColl)/fMeanr2Coll) * (TMath::Sqrt(fMeanr2Cos2PhiColl*fMeanr2Cos2PhiColl+fMeanr2Sin2PhiColl*fMeanr2Sin2PhiColl)/fMeanr2Coll) * (TMath::Sqrt(fMeanr2Cos4PhiColl*fMeanr2Cos4PhiColl+fMeanr2Sin4PhiColl*fMeanr2Sin4PhiColl)/fMeanr2Coll) * TMath::Cos(((TMath::ATan2(fMeanr2Sin2Phi,fMeanr2Cos2Phi)+TMath::Pi())/2) - ((TMath::ATan2(fMeanr2Sin4Phi,fMeanr2Cos4Phi)+TMath::Pi())/4)));
1665   return ((TMath::Sqrt(fMeanr2Cos2PhiColl*fMeanr2Cos2PhiColl+fMeanr2Sin2PhiColl*fMeanr2Sin2PhiColl)/fMeanr2Coll) * (TMath::Sqrt(fMeanr2Cos2PhiColl*fMeanr2Cos2PhiColl+fMeanr2Sin2PhiColl*fMeanr2Sin2PhiColl)/fMeanr2Coll) * (TMath::Sqrt(fMeanr4Cos4PhiColl*fMeanr4Cos4PhiColl+fMeanr4Sin4PhiColl*fMeanr4Sin4PhiColl)/fMeanr4Coll) * TMath::Cos(((TMath::ATan2(fMeanr2Sin2Phi,fMeanr2Cos2Phi)+TMath::Pi())/2) - ((TMath::ATan2(fMeanr4Sin4Phi,fMeanr4Cos4Phi)+TMath::Pi())/4)));
1666 }
1667
1668 //______________________________________________________________________________
1669 Double_t AliGlauberMC::GetE43Com() const
1670 {
1671   //get epsilon5 of binary collisions
1672   if (fNcom<2) return 0.0;
1673   //return ((TMath::Sqrt(fMeanr2Cos2PhiCom*fMeanr2Cos2PhiCom+fMeanr2Sin2PhiCom*fMeanr2Sin2PhiCom)/fMeanr2Com) * (TMath::Sqrt(fMeanr2Cos2PhiCom*fMeanr2Cos2PhiCom+fMeanr2Sin2PhiCom*fMeanr2Sin2PhiCom)/fMeanr2Com) * (TMath::Sqrt(fMeanr2Cos4PhiCom*fMeanr2Cos4PhiCom+fMeanr2Sin4PhiCom*fMeanr2Sin4PhiCom)/fMeanr2Com) * TMath::Cos(((TMath::ATan2(fMeanr2Sin2Phi,fMeanr2Cos2Phi)+TMath::Pi())/2) - ((TMath::ATan2(fMeanr2Sin4Phi,fMeanr2Cos4Phi)+TMath::Pi())/4)));
1674   return ((TMath::Sqrt(fMeanr2Cos2PhiCom*fMeanr2Cos2PhiCom+fMeanr2Sin2PhiCom*fMeanr2Sin2PhiCom)/fMeanr2Com) * (TMath::Sqrt(fMeanr2Cos2PhiCom*fMeanr2Cos2PhiCom+fMeanr2Sin2PhiCom*fMeanr2Sin2PhiCom)/fMeanr2Com) * (TMath::Sqrt(fMeanr4Cos4PhiCom*fMeanr4Cos4PhiCom+fMeanr4Sin4PhiCom*fMeanr4Sin4PhiCom)/fMeanr4Com) * TMath::Cos(((TMath::ATan2(fMeanr2Sin2Phi,fMeanr2Cos2Phi)+TMath::Pi())/2) - ((TMath::ATan2(fMeanr4Sin4Phi,fMeanr4Cos4Phi)+TMath::Pi())/4)));
1675 }
1676 //___________________________________________________________________________
1677
1678 Double_t AliGlauberMC::GetPsi4m2() const
1679 {
1680   //return (TMath::Cos(4*((TMath::ATan2(fMeanr2Sin4Phi,fMeanr2Cos4Phi)+TMath::Pi())/4) - ((TMath::ATan2(fMeanr2Sin2Phi,fMeanr2Cos2Phi)+TMath::Pi())/2)));
1681   return (TMath::Cos(4*(((TMath::ATan2(fMeanr4Sin4Phi,fMeanr4Cos4Phi)+TMath::Pi())/4)-((TMath::ATan2(fMeanr2Sin2Phi,fMeanr2Cos2Phi)+TMath::Pi())/2))));
1682 }
1683 */
1684 //______________________________________________________________________________
1685 void AliGlauberMC::Run(Int_t nevents)
1686 {
1687   //example run
1688   cout << "Generating " << nevents << " events..." << endl;
1689   TString name(Form("nt_%s_%s",fANucleus.GetName(),fBNucleus.GetName()));
1690   TString title(Form("%s + %s (x-sect = %d mb)",fANucleus.GetName(),fBNucleus.GetName(),(Int_t) fXSect));
1691   if (fnt == 0)
1692   {
1693     fnt = new TNtuple(name,title,
1694                       "Npart:Ncoll:B:MeanX:MeanY:MeanX2:MeanY2:MeanXY:VarX:VarY:VarXY:MeanXSystem:MeanYSystem:MeanXA:MeanYA:MeanXB:MeanYB:VarE:Stoa:VarEColl:VarECom:VarEPart:VarEPartColl:VarEPartCom:dNdEta:dNdEtaGBW:dNdEtaTwoNBD:xsect:tAA:Epsl2:Epsl3:Epsl4:Epsl5:E2Coll:E3Coll:E4Coll:E5Coll:E2Com:E3Com:E4Com:E5Com:Psi2:Psi3:Psi4:Psi5:BNN:signn:Ncollw");
1695     fnt->SetDirectory(0);
1696   }
1697   Int_t q = 0;
1698   Int_t u = 0;
1699   for (Int_t i = 0; i<nevents; i++)
1700   {
1701
1702     if(!NextEvent())
1703     {
1704       u++;
1705       continue;
1706     }
1707
1708     q++;
1709     Float_t v[48];
1710     v[0]  = GetNpart();
1711     v[1]  = GetNcoll();
1712     v[2]  = fBMC;
1713     v[3]  = fMeanXParts;
1714     v[4]  = fMeanYParts;
1715     v[5]  = fMeanX2Parts;
1716     v[6]  = fMeanY2Parts;
1717     v[7]  = fMeanXYParts;
1718     v[8]  = fSx2Parts;
1719     v[9]  = fSy2Parts;
1720     v[10] = fSxyParts;
1721     v[11] = fMeanXSystem;
1722     v[12] = fMeanYSystem;
1723     v[13] = fMeanXA;
1724     v[14] = fMeanYA;
1725     v[15] = fMeanXB;
1726     v[16] = fMeanYB;
1727     v[17] = GetEccentricity();
1728     v[18] = GetStoa();
1729     v[19] = GetEccentricityColl();
1730     v[20] = GetEccentricityCom();
1731     v[21] = GetEccentricityPart();
1732     v[22] = GetEccentricityPartColl();
1733     v[23] = GetEccentricityPartCom();
1734     if (fDoPartProd)
1735     {
1736       v[24] = GetdNdEta();
1737       v[25] = GetdNdEta();
1738       v[26] = v[24]+v[25];
1739     }
1740     else
1741     {
1742       v[24] = 0;
1743       v[25] = 0;
1744       v[26] = 0;
1745     }
1746     v[27]=fXSect;
1747
1748     Float_t mytAA=-999;
1749     if (GetNcoll()>0) mytAA=GetNcoll()/fXSect;
1750     v[28]=mytAA;
1751     //_____________epsilon2,3,4,4_______
1752     v[29] = GetEpsilon2Part();
1753     v[30] = GetEpsilon3Part();
1754     v[31] = GetEpsilon4Part();
1755     v[32] = GetEpsilon5Part();
1756     v[33] = GetEpsilon2Coll();
1757     v[34] = GetEpsilon3Coll();
1758     v[35] = GetEpsilon4Coll();
1759     v[36] = GetEpsilon5Coll();
1760     v[37] = GetEpsilon2Com();
1761     v[38] = GetEpsilon3Com();
1762     v[39] = GetEpsilon4Com();
1763     v[40] = GetEpsilon5Com();
1764     v[41] = GetPsi2();
1765     v[42] = GetPsi3();
1766     v[43] = GetPsi4();
1767     v[44] = GetPsi5();
1768     v[45] = fBNN;
1769     v[46] = fXSect;
1770     v[47] = fNcollw;
1771
1772     //always at the end
1773     fnt->Fill(v);
1774
1775     if ((i%100)==0) std::cout << "Generating Event # " << i << "... \r" << flush;
1776   }
1777   std::cout << "Generating Event # " << nevents << "... \r" << endl << "Done! Succesfull events:  " << q << "  discarded events:  " << u <<"."<< endl;
1778 }
1779
1780 //---------------------------------------------------------------------------------
1781 void AliGlauberMC::RunAndSaveNtuple( Int_t n,
1782                                      const Option_t *sysA,
1783                                      const Option_t *sysB,
1784                                      Double_t signn,
1785                                      Double_t mind,
1786                                      Double_t r,
1787                                      Double_t a,
1788                                      const char *fname)
1789 {
1790   //example run
1791   AliGlauberMC mcg(sysA,sysB,signn);
1792   mcg.SetMinDistance(mind);
1793   mcg.Setr(r);
1794   mcg.Seta(a);
1795   mcg.Run(n);
1796   TNtuple  *nt=mcg.GetNtuple();
1797   TFile out(fname,"recreate",fname,9);
1798   if(nt) nt->Write();
1799   printf("total cross section with a nucleon-nucleon cross section \t%f is \t%f",signn,mcg.GetTotXSect());
1800   out.Close();
1801 }
1802
1803 //---------------------------------------------------------------------------------
1804 void AliGlauberMC::RunAndSaveNucleons( Int_t n,
1805                                        const Option_t *sysA,
1806                                        Option_t *sysB,
1807                                        const Double_t signn,
1808                                        Double_t mind,
1809                                        Bool_t verbose,
1810                                        const char *fname)
1811 {
1812   //example run
1813   AliGlauberMC mcg(sysA,sysB,signn);
1814   mcg.SetMinDistance(mind);
1815   TFile *out=0;
1816   if(fname)
1817     out=new TFile(fname,"recreate",fname,9);
1818
1819   for(Int_t ievent=0; ievent<n; ievent++)
1820   {
1821
1822     //get an event with at least one collision
1823     while(!mcg.NextEvent()) {}
1824
1825     //access, save and (if wanted) print out nucleons
1826     TObjArray* nucleons=mcg.GetNucleons();
1827     if(!nucleons) continue;
1828     if(out)
1829       nucleons->Write(Form("nucleonarray%d",ievent),TObject::kSingleKey);
1830
1831     if(verbose)
1832     {
1833       cout<<endl<<endl<<"EVENT NO: "<<ievent<<endl;
1834       cout<<"B = "<<mcg.GetB()<<"  Npart = "<<mcg.GetNpart()<<endl<<endl;
1835       printf("Nucleus\t X\t Y\t Z\tNcoll\n");
1836       Int_t nNucls=nucleons->GetEntries();
1837       for(Int_t iNucl=0; iNucl<nNucls; iNucl++)
1838       {
1839         AliGlauberNucleon *nucl=(AliGlauberNucleon *)nucleons->UncheckedAt(iNucl);
1840         Char_t nucleus='A';
1841         if(nucl->IsInNucleusB()) nucleus='B';
1842         Double_t x=nucl->GetX();
1843         Double_t y=nucl->GetY();
1844         Double_t z=nucl->GetZ();
1845         Int_t ncoll=nucl->GetNColl();
1846         printf("   %c\t%2.2f\t%2.2f\t%2.2f\t%3d\n",nucleus,x,y,z,ncoll);
1847       }
1848     }
1849   }
1850   if(out) delete out;
1851 }
1852
1853 //---------------------------------------------------------------------------------
1854 void AliGlauberMC::Reset()
1855 {
1856   //delete the ntuple
1857   delete fnt;
1858   fnt=NULL;
1859 }