One mult for Glauber
authorsnelling <snelling@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 21 Mar 2011 16:49:07 +0000 (16:49 +0000)
committersnelling <snelling@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 21 Mar 2011 16:49:07 +0000 (16:49 +0000)
PWG2/FLOW/AliFlowTasks/AliAnalysisTaskPIDflowQA.cxx
PWG2/FLOW/AliFlowTasks/AliAnalysisTaskPIDflowQA.h
PWG2/FLOW/AliFlowTasks/AliFlowTrackCuts.cxx
PWG2/FLOW/AliFlowTasks/AliFlowTrackCuts.h
PWG2/FLOW/AliFlowTools/glauberMC/AliGlauberMC.cxx
PWG2/FLOW/AliFlowTools/glauberMC/AliGlauberMC.h

index 4b6d280..be8beb3 100644 (file)
@@ -61,6 +61,11 @@ AliAnalysisTaskPIDflowQA:: AliAnalysisTaskPIDflowQA():
   fTPCsignalPimc(NULL),
   fTPCsignalKmc(NULL),
   fTPCsignalPmc(NULL),
+  fTOFtime(NULL),
+  fTOFtimeE(NULL),
+  fTOFtimePi(NULL),
+  fTOFtimeK(NULL),
+  fTOFtimeP(NULL),
   fTOFbeta(NULL),
   fTOFbetaE(NULL),
   fTOFbetaPi(NULL),
@@ -170,6 +175,11 @@ AliAnalysisTaskPIDflowQA:: AliAnalysisTaskPIDflowQA(const char *name):
   fTPCsignalPimc(NULL),
   fTPCsignalKmc(NULL),
   fTPCsignalPmc(NULL),
+  fTOFtime(NULL),
+  fTOFtimeE(NULL),
+  fTOFtimePi(NULL),
+  fTOFtimeK(NULL),
+  fTOFtimeP(NULL),
   fTOFbeta(NULL),
   fTOFbetaE(NULL),
   fTOFbetaPi(NULL),
@@ -366,6 +376,17 @@ void  AliAnalysisTaskPIDflowQA::UserCreateOutputObjects()
     fOutputList->Add(fTPCsignalPmc);
   }
 
+  fTOFtime=new TH2F("fTOFtime",";p[GeV/c];#time",kPBins,binsPDummy,1000, 0.4, 1.1);//
+  fOutputList->Add(fTOFtime);
+  fTOFtimeE=new TH2F("fTOFtimeE",";p [GeV/c];#time-#time_{#pi}",kPBins,binsPDummy,500, -0.25, 0.25);//
+  fTOFtimePi=new TH2F("fTOFtimePi",";p [GeV/c];#time-#time_{#pi}",kPBins,binsPDummy,500, -0.25, 0.25);//
+  fTOFtimeK=new TH2F("fTOFtimeK",";p [GeV/c];#time-#time_{K}",kPBins,binsPDummy,500, -0.25, 0.25);//
+  fTOFtimeP=new TH2F("fTOFtimeP",";p [GeV/c];#time-#time_{p}",kPBins,binsPDummy,500, -0.25, 0.25);//
+  //fOutputList->Add(fTOFtimeE);
+  //fOutputList->Add(fTOFtimePi);
+  //fOutputList->Add(fTOFtimeK);
+  //fOutputList->Add(fTOFtimeP);
+
   fTOFbeta=new TH2F("fTOFbeta",";p[GeV/c];#beta",kPBins,binsPDummy,1000, 0.4, 1.1);//
   fOutputList->Add(fTOFbeta);
   fTOFbetaE=new TH2F("fTOFbetaE",";p [GeV/c];#beta-#beta_{#pi}",kPBins,binsPDummy,500, -0.25, 0.25);//
index 344e45e..8c6b3c9 100644 (file)
@@ -52,6 +52,12 @@ private:
   TH2F* fTPCsignalKmc;//!TPC PID signal as function of p for K+
   TH2F* fTPCsignalPmc;//!TPC PID signal as function of p for p
 
+  TH2F* fTOFtime;//!vs time
+  TH2F* fTOFtimeE;//!vs time
+  TH2F* fTOFtimePi;//!vs time
+  TH2F* fTOFtimeK;//!vs time
+  TH2F* fTOFtimeP;//!vs time
+
   TH2F* fTOFbeta;//!vs beta
   TH2F* fTOFbetaE;//!vs beta
   TH2F* fTOFbetaPi;//!vs beta
index 92695a7..4463ad7 100644 (file)
@@ -1227,16 +1227,16 @@ Bool_t AliFlowTrackCuts::PassesTOFbetaSimpleCut(const AliESDtrack* track )
   {
     case AliPID::kPion:
       return ( (s[2]<0.015) && (s[2]>-0.015) &&
-               (s[3]>0.025) && (s[3]<-0.025) &&
-               (s[4]>0.03) && (s[4]<-0.03) );
+               (s[3]>0.025) &&
+               (s[4]>0.03) );
     case AliPID::kKaon:
       return ( (s[3]<0.015) && (s[3]>-0.015) &&
-               (s[2]>0.03) && (s[2]<-0.03) &&
-               (s[4]<-0.03) && (s[4]>0.03) );
+               (s[2]<-0.03) &&
+               (s[4]>0.03) );
     case AliPID::kProton:
       return ( (s[4]<0.015) && (s[2]>-0.015) &&
-               (s[3]<-0.025) && (s[3]>0.025) &&
-               (s[2]<-0.025) && (s[2]>0.025) );
+               (s[3]<-0.025) &&
+               (s[2]<-0.025) );
     default:
       return kFALSE;
   }
index 9d1153e..5308b8a 100644 (file)
@@ -203,7 +203,7 @@ class AliFlowTrackCuts : public AliFlowTrackSimpleCuts {
   void HandleVParticle(AliVParticle* track);
   void DefineHistograms();
   void InitPIDcuts();
-  void InitESDcuts() {if (!fAliESDtrackCuts) fAliESDtrackCuts=new AliESDtrackCuts();}
+  void InitESDcuts() {if (!fAliESDtrackCuts) {fAliESDtrackCuts=new AliESDtrackCuts();}}
   // part added by F. Noferini
   Bool_t PassesTOFbayesianCut(AliESDtrack* track); 
   void SetPriors(); // set my favourite priors
index 1432e92..49c2ba7 100644 (file)
@@ -19,7 +19,7 @@
 //
 //  origin: PHOBOS experiment
 //  alification: Mikolaj Krzewicki, Nikhef, mikolaj.krzewicki@cern.ch
-//  update:      You Zhou, Nikhef, yzhou@nikhef.nl 
+//  update:      You Zhou, Nikhef, yzhou@nikhef.nl
 ////////////////////////////////////////////////////////////////////////////////
 
 #include <Riostream.h>
@@ -31,9 +31,6 @@
 #include <TNtuple.h>
 #include <TFile.h>
 #include <TTree.h>
-#include <TF1.h>
-#include <TH1F.h>
-#include <TArray.h>
 
 #include "AliGlauberNucleon.h"
 #include "AliGlauberNucleus.h"
@@ -73,6 +70,7 @@ AliGlauberMC::AliGlauberMC(Option_t* NA, Option_t* NB, Double_t xsect) :
   fTotalEvents(0),
   fBMin(0.),
   fBMax(20.),
+  fMultType(kNBDSV),
   fMaxNpartFound(0),
   fNpart(0),
   fNcoll(0),
@@ -123,20 +121,10 @@ AliGlauberMC::AliGlauberMC(Option_t* NA, Option_t* NB, Double_t xsect) :
   fDoPartProd(kFALSE)
 {
   //ctor
-  fdNdEtaParam[0] = 8.0;
-  fdNdEtaParam[1] = 0.13;
-  fdNdEtaGBWParam[0] = 0.79;
-  fdNdEtaGBWParam[1] = 0.288;
-  fdNdEtaGBWParam[2] = 30.25;
-  fdNdEtaNBDParam[0] = 3.0;
-  fdNdEtaNBDParam[1] = 4.0;
-  fdNdEtaNBDParam[2] = 0.13;
-  fdNdEtaTwoNBDParam[0] = 3.0;
-  fdNdEtaTwoNBDParam[1] = 4.0;
-  fdNdEtaTwoNBDParam[2] = 2.0;
-  fdNdEtaTwoNBDParam[3] = 11.0;
-  fdNdEtaTwoNBDParam[4] = 0.4;
-  fdNdEtaTwoNBDParam[5] = 0.13;
+  for (UInt_t i=0; i<(sizeof(fdNdEtaParam)/sizeof(fdNdEtaParam[0])); i++)
+  {
+    fdNdEtaParam[i]=0.0;
+  }
 
   SetName(Form("Glauber_%s_%s",fANucleus.GetName(),fBNucleus.GetName()));
   SetTitle(Form("Glauber %s+%s Version",fANucleus.GetName(),fBNucleus.GetName()));
@@ -181,6 +169,7 @@ AliGlauberMC::AliGlauberMC(const AliGlauberMC& in):
   fTotalEvents(in.fTotalEvents),
   fBMin(in.fBMin),
   fBMax(in.fBMax),
+  fMultType(in.fMultType),
   fMaxNpartFound(in.fMaxNpartFound),
   fNpart(in.fNpart),
   fNcoll(in.fNcoll),
@@ -231,27 +220,24 @@ AliGlauberMC::AliGlauberMC(const AliGlauberMC& in):
   fDoPartProd(kFALSE)
 {
   //copy ctor
-  memcpy(fdNdEtaParam,in.fdNdEtaParam,2*sizeof(Double_t));
-  memcpy(fdNdEtaGBWParam,in.fdNdEtaGBWParam,3*sizeof(Double_t));
-  memcpy(fdNdEtaNBDParam,in.fdNdEtaNBDParam,3*sizeof(Double_t));
-  memcpy(fdNdEtaTwoNBDParam,in.fdNdEtaTwoNBDParam,6*sizeof(Double_t));
+  memcpy(fdNdEtaParam,in.fdNdEtaParam,sizeof(fdNdEtaParam));
 }
 
 //______________________________________________________________________________
 AliGlauberMC& AliGlauberMC::operator=(const AliGlauberMC& in)
 {
   //assignment
-  fANucleus=in.fANucleus; 
-  fBNucleus=in.fBNucleus; 
-  fXSect=in.fXSect;    
+  fANucleus=in.fANucleus;
+  fBNucleus=in.fBNucleus;
+  fXSect=in.fXSect;
   fNucleonsA=in.fNucleonsA;
   fNucleonsB=in.fNucleonsB;
-  fAN=in.fAN;       
-  fBN=in.fBN;       
-  fnt=in.fnt;       
-  fMeanX2=in.fMeanX2;   
-  fMeanY2=in.fMeanY2;   
-  fMeanXY=in.fMeanXY; 
+  fAN=in.fAN;
+  fBN=in.fBN;
+  fnt=in.fnt;
+  fMeanX2=in.fMeanX2;
+  fMeanY2=in.fMeanY2;
+  fMeanXY=in.fMeanXY;
   fMeanr2=in.fMeanr2;
   fMeanr3=in.fMeanr3;
   fMeanr4=in.fMeanr4;
@@ -269,22 +255,20 @@ AliGlauberMC& AliGlauberMC::operator=(const AliGlauberMC& in)
   fMeanr5Coll=in.fMeanr5Coll;
   fMeanXSystem=in.fMeanXSystem;
   fMeanYSystem=in.fMeanYSystem;
-  fMeanXA=in.fMeanXA;  
-  fMeanYA=in.fMeanYA;  
-  fMeanXB=in.fMeanXB;  
-  fMeanYB=in.fMeanYB;  
-  fBMC=in.fBMC;     
-  fEvents=in.fEvents;   
+  fMeanXA=in.fMeanXA;
+  fMeanYA=in.fMeanYA;
+  fMeanXB=in.fMeanXB;
+  fMeanYB=in.fMeanYB;
+  fBMC=in.fBMC;
+  fEvents=in.fEvents;
   fTotalEvents=in.fTotalEvents;
-  fBMin=in.fBMin;     
-  fBMax=in.fBMax;     
-  memcpy(fdNdEtaParam,in.fdNdEtaParam,2*sizeof(Double_t));
-  memcpy(fdNdEtaGBWParam,in.fdNdEtaGBWParam,3*sizeof(Double_t));
-  memcpy(fdNdEtaNBDParam,in.fdNdEtaNBDParam,3*sizeof(Double_t));
-  memcpy(fdNdEtaTwoNBDParam,in.fdNdEtaTwoNBDParam,6*sizeof(Double_t));
+  fBMin=in.fBMin;
+  fBMax=in.fBMax;
+  fMultType=in.fMultType,
+  memcpy(fdNdEtaParam,in.fdNdEtaParam,sizeof(fdNdEtaParam));
   fMaxNpartFound=in.fMaxNpartFound;
-  fNpart=in.fNpart;   
-  fNcoll=in.fNcoll; 
+  fNpart=in.fNpart;
+  fNcoll=in.fNcoll;
   fMeanr2Cos2Phi=in.fMeanr2Cos2Phi;
   fMeanr2Sin2Phi=in.fMeanr2Sin2Phi;
   fMeanr2Cos3Phi=in.fMeanr2Cos3Phi;
@@ -313,14 +297,14 @@ AliGlauberMC& AliGlauberMC::operator=(const AliGlauberMC& in)
   fMeanr4Sin4PhiColl=in.fMeanr4Sin4PhiColl;
   fMeanr5Cos5PhiColl=in.fMeanr5Cos5PhiColl;
   fMeanr5Sin5PhiColl=in.fMeanr5Sin5PhiColl;
-  fSx2=in.fSx2;      
-  fSy2=in.fSy2;      
-  fSxy=in.fSxy;      
-  fSx2Coll=in.fSx2Coll;  
-  fSy2Coll=in.fSy2Coll;  
-  fSxyColl=in.fSxyColl;  
-  fX=in.fX;      
-  fNpp=in.fNpp;       
+  fSx2=in.fSx2;
+  fSy2=in.fSy2;
+  fSxy=in.fSxy;
+  fSx2Coll=in.fSx2Coll;
+  fSy2Coll=in.fSy2Coll;
+  fSxyColl=in.fSxyColl;
+  fX=in.fX;
+  fNpp=in.fNpp;
   return *this;
 }
 
@@ -332,40 +316,40 @@ Bool_t AliGlauberMC::CalcEvent(Double_t bgen)
   fNucleonsA = fANucleus.GetNucleons();
   fAN = fANucleus.GetN();
   for (Int_t i = 0; i<fAN; i++)
-    {
-      AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->UncheckedAt(i));
-      nucleonA->SetInNucleusA();
-    }
+  {
+    AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->UncheckedAt(i));
+    nucleonA->SetInNucleusA();
+  }
   fBNucleus.ThrowNucleons(bgen/2.);
   fNucleonsB = fBNucleus.GetNucleons();
   fBN = fBNucleus.GetN();
   for (Int_t i = 0; i<fBN; i++)
-    {
-      AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->UncheckedAt(i));
-      nucleonB->SetInNucleusB();
-    }
-  
+  {
+    AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->UncheckedAt(i));
+    nucleonB->SetInNucleusB();
+  }
+
   // "ball" diameter = distance at which two balls interact
   Double_t d2 = (Double_t)fXSect/(TMath::Pi()*10); // in fm^2
-  
+
   // for each of the A nucleons in nucleus B
   for (Int_t i = 0; i<fBN; i++)
+  {
+    AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->UncheckedAt(i));
+    for (Int_t j = 0 ; j < fAN ; j++)
     {
-      AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->UncheckedAt(i));
-      for (Int_t j = 0 ; j < fAN ; j++)
-       {
-         AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->UncheckedAt(j));
-         Double_t dx = nucleonB->GetX()-nucleonA->GetX();
-         Double_t dy = nucleonB->GetY()-nucleonA->GetY();
-         Double_t dij = dx*dx+dy*dy;
-         if (dij < d2)
-           {
-             nucleonB->Collide();
-             nucleonA->Collide();
-           }
-       }
+      AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->UncheckedAt(j));
+      Double_t dx = nucleonB->GetX()-nucleonA->GetX();
+      Double_t dy = nucleonB->GetY()-nucleonA->GetY();
+      Double_t dij = dx*dx+dy*dy;
+      if (dij < d2)
+      {
+        nucleonB->Collide();
+        nucleonA->Collide();
+      }
     }
-  
+  }
+
   return CalcResults(bgen);
 }
 
@@ -374,7 +358,7 @@ Bool_t AliGlauberMC::CalcResults(Double_t bgen)
 {
   // calc results for the given event
   //return true if we have participants
-  
+
   fNpart=0;
   fNcoll=0;
   fMeanX2=0.;
@@ -429,272 +413,272 @@ Bool_t AliGlauberMC::CalcResults(Double_t bgen)
   fMeanr4Sin4PhiColl=0.;
   fMeanr5Cos5PhiColl=0.;
   fMeanr5Sin5PhiColl=0.;
-  
+
   for (Int_t i = 0; i<fAN; i++)
+  {
+    AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->UncheckedAt(i));
+    Double_t xA = nucleonA->GetX();
+    Double_t yA = nucleonA->GetY();
+    Double_t r2A = xA*xA+yA*yA;
+    Double_t rA = TMath::Sqrt(r2A);
+    Double_t r3A = r2A*rA;
+    Double_t r4A = r3A*rA;
+    Double_t r5A = r4A*rA;
+    Double_t phiA = TMath::ATan2(yA,xA);
+    Double_t sin2PhiA = TMath::Sin(2.*phiA);
+    Double_t cos2PhiA = TMath::Cos(2.*phiA);
+    Double_t sin3PhiA = TMath::Sin(3.*phiA);
+    Double_t cos3PhiA = TMath::Cos(3.*phiA);
+    Double_t sin4PhiA = TMath::Sin(4.*phiA);
+    Double_t cos4PhiA = TMath::Cos(4.*phiA);
+    Double_t sin5PhiA = TMath::Sin(5.*phiA);
+    Double_t cos5PhiA = TMath::Cos(5.*phiA);
+    //printf("x: %.3f, y:%.3f, r:%.3f, phi:%.3f, sp:%.6f,cp:%.6f\n",xA,yA,TMath::Sqrt(r2A),PhiA,Sin2PhiA, Cos2PhiA);
+    fMeanXSystem  += xA;
+    fMeanYSystem  += yA;
+    fMeanXA  += xA;
+    fMeanYA  += yA;
+
+    if(nucleonA->IsWounded())
     {
-      AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->UncheckedAt(i));
-      Double_t xA = nucleonA->GetX();
-      Double_t yA = nucleonA->GetY();
-      Double_t r2A = xA*xA+yA*yA;
-      Double_t rA = TMath::Sqrt(r2A);
-      Double_t r3A = r2A*rA;
-      Double_t r4A = r3A*rA;
-      Double_t r5A = r4A*rA;
-      Double_t phiA = TMath::ATan2(yA,xA);
-      Double_t sin2PhiA = TMath::Sin(2.*phiA);
-      Double_t cos2PhiA = TMath::Cos(2.*phiA);
-      Double_t sin3PhiA = TMath::Sin(3.*phiA);
-      Double_t cos3PhiA = TMath::Cos(3.*phiA);
-      Double_t sin4PhiA = TMath::Sin(4.*phiA);
-      Double_t cos4PhiA = TMath::Cos(4.*phiA);
-      Double_t sin5PhiA = TMath::Sin(5.*phiA);
-      Double_t cos5PhiA = TMath::Cos(5.*phiA);
-      //printf("x: %.3f, y:%.3f, r:%.3f, phi:%.3f, sp:%.6f,cp:%.6f\n",xA,yA,TMath::Sqrt(r2A),PhiA,Sin2PhiA, Cos2PhiA);
-      fMeanXSystem  += xA;
-      fMeanYSystem  += yA;
-      fMeanXA  += xA;
-      fMeanYA  += yA;
-      
-      if(nucleonA->IsWounded())
-       {
-         fNpart++;
-         fMeanXParts  += xA;
-         fMeanYParts  += yA;
-         fMeanX2 += xA * xA;
-         fMeanY2 += yA * yA;
-         fMeanXY += xA * yA;
-         fMeanr2 += r2A;
-         fMeanr3 += r3A;
-         fMeanr4 += r4A;
-         fMeanr5 += r5A;
-         fMeanr2Cos2Phi += r2A*cos2PhiA;
-         fMeanr2Sin2Phi += r2A*sin2PhiA;
-         fMeanr2Cos3Phi += r2A*cos3PhiA;
-         fMeanr2Sin3Phi += r2A*sin3PhiA;
-         fMeanr2Cos4Phi += r2A*cos4PhiA;
-         fMeanr2Sin4Phi += r2A*sin4PhiA;
-         fMeanr2Cos5Phi += r2A*cos5PhiA;
-         fMeanr2Sin5Phi += r2A*sin5PhiA;
-         fMeanr3Cos3Phi += r3A*cos3PhiA;
-         fMeanr3Sin3Phi += r3A*sin3PhiA;
-         fMeanr4Cos4Phi += r4A*cos4PhiA;
-         fMeanr4Sin4Phi += r4A*sin4PhiA;
-         fMeanr5Cos5Phi += r5A*cos5PhiA;
-         fMeanr5Sin5Phi += r5A*sin5PhiA;
-       }
+      fNpart++;
+      fMeanXParts  += xA;
+      fMeanYParts  += yA;
+      fMeanX2 += xA * xA;
+      fMeanY2 += yA * yA;
+      fMeanXY += xA * yA;
+      fMeanr2 += r2A;
+      fMeanr3 += r3A;
+      fMeanr4 += r4A;
+      fMeanr5 += r5A;
+      fMeanr2Cos2Phi += r2A*cos2PhiA;
+      fMeanr2Sin2Phi += r2A*sin2PhiA;
+      fMeanr2Cos3Phi += r2A*cos3PhiA;
+      fMeanr2Sin3Phi += r2A*sin3PhiA;
+      fMeanr2Cos4Phi += r2A*cos4PhiA;
+      fMeanr2Sin4Phi += r2A*sin4PhiA;
+      fMeanr2Cos5Phi += r2A*cos5PhiA;
+      fMeanr2Sin5Phi += r2A*sin5PhiA;
+      fMeanr3Cos3Phi += r3A*cos3PhiA;
+      fMeanr3Sin3Phi += r3A*sin3PhiA;
+      fMeanr4Cos4Phi += r4A*cos4PhiA;
+      fMeanr4Sin4Phi += r4A*sin4PhiA;
+      fMeanr5Cos5Phi += r5A*cos5PhiA;
+      fMeanr5Sin5Phi += r5A*sin5PhiA;
     }
-  
+  }
+
   for (Int_t i = 0; i<fBN; i++)
+  {
+    AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->UncheckedAt(i));
+    Double_t xB=nucleonB->GetX();
+    Double_t yB=nucleonB->GetY();
+    Double_t r2B=xB*xB+yB*yB;
+    Double_t rB = TMath::Sqrt(r2B);
+    Double_t r3B = r2B*rB;
+    Double_t r4B = r3B*rB;
+    Double_t r5B = r4B*rB;
+    Double_t phiB = TMath::ATan2(yB,xB);
+    Double_t sin2PhiB = TMath::Sin(2*phiB);
+    Double_t cos2PhiB = TMath::Cos(2*phiB);
+    Double_t sin3PhiB = TMath::Sin(3*phiB);
+    Double_t cos3PhiB = TMath::Cos(3*phiB);
+    Double_t sin4PhiB = TMath::Sin(4*phiB);
+    Double_t cos4PhiB = TMath::Cos(4*phiB);
+    Double_t sin5PhiB = TMath::Sin(5*phiB);
+    Double_t cos5PhiB = TMath::Cos(5*phiB);
+    fMeanXSystem  += xB;
+    fMeanYSystem  += yB;
+    fMeanXB  += xB;
+    fMeanYB  += yB;
+
+    if(nucleonB->IsWounded())
     {
-      AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->UncheckedAt(i));
-      Double_t xB=nucleonB->GetX();
-      Double_t yB=nucleonB->GetY();
-      Double_t r2B=xB*xB+yB*yB;
-      Double_t rB = TMath::Sqrt(r2B);
-      Double_t r3B = r2B*rB;
-      Double_t r4B = r3B*rB;
-      Double_t r5B = r4B*rB;
-      Double_t phiB = TMath::ATan2(yB,xB);
-      Double_t sin2PhiB = TMath::Sin(2*phiB);
-      Double_t cos2PhiB = TMath::Cos(2*phiB);
-      Double_t sin3PhiB = TMath::Sin(3*phiB);
-      Double_t cos3PhiB = TMath::Cos(3*phiB);
-      Double_t sin4PhiB = TMath::Sin(4*phiB);
-      Double_t cos4PhiB = TMath::Cos(4*phiB);
-      Double_t sin5PhiB = TMath::Sin(5*phiB);
-      Double_t cos5PhiB = TMath::Cos(5*phiB);
-      fMeanXSystem  += xB;
-      fMeanYSystem  += yB;
-      fMeanXB  += xB;
-      fMeanYB  += yB;
-      
-      if(nucleonB->IsWounded())
-       {
-         Int_t ncoll = nucleonB->GetNColl();
-         fNpart++;
-         fMeanXParts  += xB;
-         fMeanXColl  += xB*ncoll;
-         fMeanYParts  += yB;
-         fMeanYColl += yB*ncoll;
-         fMeanX2 += xB * xB;
-         fMeanX2Coll += xB*xB*ncoll*ncoll;
-         fMeanY2 += yB * yB;
-         fMeanY2Coll += yB*yB*ncoll*ncoll;
-         fMeanXY += xB * yB;
-         fMeanXYColl += xB*yB*ncoll*ncoll;
-         fNcoll += nucleonB->GetNColl();
-         fMeanr2 += r2B;
-         fMeanr3 += r3B;
-         fMeanr4 += r4B;
-         fMeanr5 += r5B;
-         fMeanr2Cos2Phi += r2B*cos2PhiB;
-         fMeanr2Sin2Phi += r2B*sin2PhiB;
-         fMeanr2Cos3Phi += r2B*cos3PhiB;
-         fMeanr2Sin3Phi += r2B*sin3PhiB;
-         fMeanr2Cos4Phi += r2B*cos4PhiB;
-         fMeanr2Sin4Phi += r2B*sin4PhiB;
-         fMeanr2Cos5Phi += r2B*cos5PhiB;
-         fMeanr2Sin5Phi += r2B*sin5PhiB;
-         fMeanr3Cos3Phi += r3B*cos3PhiB;
-         fMeanr3Sin3Phi += r3B*sin3PhiB;
-         fMeanr4Cos4Phi += r4B*cos4PhiB;
-         fMeanr4Sin4Phi += r4B*sin4PhiB;
-         fMeanr5Cos5Phi += r5B*cos5PhiB;
-         fMeanr5Sin5Phi += r5B*sin5PhiB;
-         fMeanr2Coll += r2B*ncoll*ncoll;
-         fMeanr3Coll += r3B*ncoll*ncoll;
-         fMeanr4Coll += r4B*ncoll*ncoll;
-         fMeanr5Coll += r5B*ncoll*ncoll;
-         fMeanr2Cos2PhiColl += r2B*cos2PhiB*ncoll*ncoll;
-         fMeanr2Sin2PhiColl += r2B*sin2PhiB*ncoll*ncoll;
-         fMeanr2Cos3PhiColl += r2B*cos3PhiB*ncoll*ncoll;
-         fMeanr2Sin3PhiColl += r2B*sin3PhiB*ncoll*ncoll;
-         fMeanr2Cos4PhiColl += r2B*cos4PhiB*ncoll*ncoll;
-         fMeanr2Sin4PhiColl += r2B*sin4PhiB*ncoll*ncoll;
-         fMeanr2Cos5PhiColl += r2B*cos5PhiB*ncoll*ncoll;
-         fMeanr2Sin5PhiColl += r2B*sin5PhiB*ncoll*ncoll;
-         fMeanr3Cos3PhiColl += r3B*cos3PhiB*ncoll*ncoll;
-         fMeanr3Sin3PhiColl += r3B*sin3PhiB*ncoll*ncoll;
-         fMeanr4Cos4PhiColl += r4B*cos4PhiB*ncoll*ncoll;
-         fMeanr4Sin4PhiColl += r4B*sin4PhiB*ncoll*ncoll;
-         fMeanr5Cos5PhiColl += r5B*cos5PhiB*ncoll*ncoll;
-         fMeanr5Sin5PhiColl += r5B*sin5PhiB*ncoll*ncoll;
-       }
+      Int_t ncoll = nucleonB->GetNColl();
+      fNpart++;
+      fMeanXParts  += xB;
+      fMeanXColl  += xB*ncoll;
+      fMeanYParts  += yB;
+      fMeanYColl += yB*ncoll;
+      fMeanX2 += xB * xB;
+      fMeanX2Coll += xB*xB*ncoll*ncoll;
+      fMeanY2 += yB * yB;
+      fMeanY2Coll += yB*yB*ncoll*ncoll;
+      fMeanXY += xB * yB;
+      fMeanXYColl += xB*yB*ncoll*ncoll;
+      fNcoll += nucleonB->GetNColl();
+      fMeanr2 += r2B;
+      fMeanr3 += r3B;
+      fMeanr4 += r4B;
+      fMeanr5 += r5B;
+      fMeanr2Cos2Phi += r2B*cos2PhiB;
+      fMeanr2Sin2Phi += r2B*sin2PhiB;
+      fMeanr2Cos3Phi += r2B*cos3PhiB;
+      fMeanr2Sin3Phi += r2B*sin3PhiB;
+      fMeanr2Cos4Phi += r2B*cos4PhiB;
+      fMeanr2Sin4Phi += r2B*sin4PhiB;
+      fMeanr2Cos5Phi += r2B*cos5PhiB;
+      fMeanr2Sin5Phi += r2B*sin5PhiB;
+      fMeanr3Cos3Phi += r3B*cos3PhiB;
+      fMeanr3Sin3Phi += r3B*sin3PhiB;
+      fMeanr4Cos4Phi += r4B*cos4PhiB;
+      fMeanr4Sin4Phi += r4B*sin4PhiB;
+      fMeanr5Cos5Phi += r5B*cos5PhiB;
+      fMeanr5Sin5Phi += r5B*sin5PhiB;
+      fMeanr2Coll += r2B*ncoll*ncoll;
+      fMeanr3Coll += r3B*ncoll*ncoll;
+      fMeanr4Coll += r4B*ncoll*ncoll;
+      fMeanr5Coll += r5B*ncoll*ncoll;
+      fMeanr2Cos2PhiColl += r2B*cos2PhiB*ncoll*ncoll;
+      fMeanr2Sin2PhiColl += r2B*sin2PhiB*ncoll*ncoll;
+      fMeanr2Cos3PhiColl += r2B*cos3PhiB*ncoll*ncoll;
+      fMeanr2Sin3PhiColl += r2B*sin3PhiB*ncoll*ncoll;
+      fMeanr2Cos4PhiColl += r2B*cos4PhiB*ncoll*ncoll;
+      fMeanr2Sin4PhiColl += r2B*sin4PhiB*ncoll*ncoll;
+      fMeanr2Cos5PhiColl += r2B*cos5PhiB*ncoll*ncoll;
+      fMeanr2Sin5PhiColl += r2B*sin5PhiB*ncoll*ncoll;
+      fMeanr3Cos3PhiColl += r3B*cos3PhiB*ncoll*ncoll;
+      fMeanr3Sin3PhiColl += r3B*sin3PhiB*ncoll*ncoll;
+      fMeanr4Cos4PhiColl += r4B*cos4PhiB*ncoll*ncoll;
+      fMeanr4Sin4PhiColl += r4B*sin4PhiB*ncoll*ncoll;
+      fMeanr5Cos5PhiColl += r5B*cos5PhiB*ncoll*ncoll;
+      fMeanr5Sin5PhiColl += r5B*sin5PhiB*ncoll*ncoll;
     }
-  
+  }
+
   if (fNpart>0)
-    {
-      fMeanXParts /= fNpart;
-      fMeanYParts /= fNpart;
-      fMeanX2 /= fNpart;
-      fMeanY2 /= fNpart;
-      fMeanXY /= fNpart;
-      fMeanr2 /= fNpart;
-      fMeanr3 /= fNpart;
-      fMeanr4 /= fNpart;
-      fMeanr5 /= fNpart;
-      fMeanr2Cos2Phi /= fNpart;
-      fMeanr2Sin2Phi /= fNpart;
-      fMeanr2Cos3Phi /= fNpart;
-      fMeanr2Sin3Phi /= fNpart;
-      fMeanr2Cos4Phi /= fNpart;
-      fMeanr2Sin4Phi /= fNpart;
-      fMeanr2Cos5Phi /= fNpart;
-      fMeanr2Sin5Phi /= fNpart;
-      fMeanr3Cos3Phi /= fNpart;
-      fMeanr3Sin3Phi /= fNpart;
-      fMeanr4Cos4Phi /= fNpart;
-      fMeanr4Sin4Phi /= fNpart;
-      fMeanr5Cos5Phi /= fNpart;
-      fMeanr5Sin5Phi /= fNpart;
-    }
+  {
+    fMeanXParts /= fNpart;
+    fMeanYParts /= fNpart;
+    fMeanX2 /= fNpart;
+    fMeanY2 /= fNpart;
+    fMeanXY /= fNpart;
+    fMeanr2 /= fNpart;
+    fMeanr3 /= fNpart;
+    fMeanr4 /= fNpart;
+    fMeanr5 /= fNpart;
+    fMeanr2Cos2Phi /= fNpart;
+    fMeanr2Sin2Phi /= fNpart;
+    fMeanr2Cos3Phi /= fNpart;
+    fMeanr2Sin3Phi /= fNpart;
+    fMeanr2Cos4Phi /= fNpart;
+    fMeanr2Sin4Phi /= fNpart;
+    fMeanr2Cos5Phi /= fNpart;
+    fMeanr2Sin5Phi /= fNpart;
+    fMeanr3Cos3Phi /= fNpart;
+    fMeanr3Sin3Phi /= fNpart;
+    fMeanr4Cos4Phi /= fNpart;
+    fMeanr4Sin4Phi /= fNpart;
+    fMeanr5Cos5Phi /= fNpart;
+    fMeanr5Sin5Phi /= fNpart;
+  }
   else
-    {
-      fMeanXParts = 0;
-      fMeanYParts = 0;
-      fMeanX2 = 0;
-      fMeanY2 = 0;
-      fMeanXY = 0;
-      fMeanr2 = 0;
-      fMeanr3 = 0;
-      fMeanr4 = 0;
-      fMeanr5 = 0;
-      fMeanr2Cos2Phi = 0;
-      fMeanr2Sin2Phi = 0;
-      fMeanr2Cos3Phi = 0;
-      fMeanr2Sin3Phi = 0;
-      fMeanr2Cos4Phi = 0;
-      fMeanr2Sin4Phi = 0;
-      fMeanr2Cos5Phi = 0;
-      fMeanr2Sin5Phi = 0;
-      fMeanr3Cos3Phi = 0;
-      fMeanr3Sin3Phi = 0;
-      fMeanr4Cos4Phi = 0;
-      fMeanr4Sin4Phi = 0;
-      fMeanr5Cos5Phi = 0;
-      fMeanr5Sin5Phi = 0;
-    }
-  
+  {
+    fMeanXParts = 0;
+    fMeanYParts = 0;
+    fMeanX2 = 0;
+    fMeanY2 = 0;
+    fMeanXY = 0;
+    fMeanr2 = 0;
+    fMeanr3 = 0;
+    fMeanr4 = 0;
+    fMeanr5 = 0;
+    fMeanr2Cos2Phi = 0;
+    fMeanr2Sin2Phi = 0;
+    fMeanr2Cos3Phi = 0;
+    fMeanr2Sin3Phi = 0;
+    fMeanr2Cos4Phi = 0;
+    fMeanr2Sin4Phi = 0;
+    fMeanr2Cos5Phi = 0;
+    fMeanr2Sin5Phi = 0;
+    fMeanr3Cos3Phi = 0;
+    fMeanr3Sin3Phi = 0;
+    fMeanr4Cos4Phi = 0;
+    fMeanr4Sin4Phi = 0;
+    fMeanr5Cos5Phi = 0;
+    fMeanr5Sin5Phi = 0;
+  }
+
   if (fNcoll>0)
-    {
-      fMeanXColl /= fNcoll;
-      fMeanYColl /= fNcoll;
-      fMeanX2Coll /= fNcoll;
-      fMeanY2Coll /= fNcoll;
-      fMeanXYColl /= fNcoll;
-      fMeanr2Coll /= fNcoll;
-      fMeanr3Coll /= fNcoll;
-      fMeanr4Coll /= fNcoll;
-      fMeanr5Coll /= fNcoll;
-      fMeanr2Cos2PhiColl /= fNcoll;
-      fMeanr2Sin2PhiColl /= fNcoll;
-      fMeanr2Cos3PhiColl /= fNcoll;
-      fMeanr2Sin3PhiColl /= fNcoll;
-      fMeanr2Cos4PhiColl /= fNcoll;
-      fMeanr2Sin4PhiColl /= fNcoll;
-      fMeanr2Cos5PhiColl /= fNcoll;
-      fMeanr2Sin5PhiColl /= fNcoll;
-      fMeanr3Cos3PhiColl /= fNcoll;
-      fMeanr3Sin3PhiColl /= fNcoll;
-      fMeanr4Cos4PhiColl /= fNcoll;
-      fMeanr4Sin4PhiColl /= fNcoll;
-      fMeanr5Cos5PhiColl /= fNcoll;
-      fMeanr5Sin5PhiColl /= fNcoll;
-    }
+  {
+    fMeanXColl /= fNcoll;
+    fMeanYColl /= fNcoll;
+    fMeanX2Coll /= fNcoll;
+    fMeanY2Coll /= fNcoll;
+    fMeanXYColl /= fNcoll;
+    fMeanr2Coll /= fNcoll;
+    fMeanr3Coll /= fNcoll;
+    fMeanr4Coll /= fNcoll;
+    fMeanr5Coll /= fNcoll;
+    fMeanr2Cos2PhiColl /= fNcoll;
+    fMeanr2Sin2PhiColl /= fNcoll;
+    fMeanr2Cos3PhiColl /= fNcoll;
+    fMeanr2Sin3PhiColl /= fNcoll;
+    fMeanr2Cos4PhiColl /= fNcoll;
+    fMeanr2Sin4PhiColl /= fNcoll;
+    fMeanr2Cos5PhiColl /= fNcoll;
+    fMeanr2Sin5PhiColl /= fNcoll;
+    fMeanr3Cos3PhiColl /= fNcoll;
+    fMeanr3Sin3PhiColl /= fNcoll;
+    fMeanr4Cos4PhiColl /= fNcoll;
+    fMeanr4Sin4PhiColl /= fNcoll;
+    fMeanr5Cos5PhiColl /= fNcoll;
+    fMeanr5Sin5PhiColl /= fNcoll;
+  }
   else
-    {
-      fMeanXColl = 0;
-      fMeanYColl = 0;
-      fMeanX2Coll = 0;
-      fMeanY2Coll = 0;
-      fMeanXYColl = 0;
-      fMeanr2Cos2PhiColl =0;
-      fMeanr2Sin2PhiColl =0;
-      fMeanr2Cos3PhiColl =0;
-      fMeanr2Sin3PhiColl =0;
-      fMeanr2Cos4PhiColl =0;
-      fMeanr2Sin4PhiColl =0;
-      fMeanr2Cos5PhiColl =0;
-      fMeanr2Sin5PhiColl =0;
-      fMeanr3Cos3PhiColl =0;
-      fMeanr3Sin3PhiColl =0;
-      fMeanr4Cos4PhiColl =0;
-      fMeanr4Sin4PhiColl =0;
-      fMeanr5Cos5PhiColl =0;
-      fMeanr5Sin5PhiColl =0;
-    }
+  {
+    fMeanXColl = 0;
+    fMeanYColl = 0;
+    fMeanX2Coll = 0;
+    fMeanY2Coll = 0;
+    fMeanXYColl = 0;
+    fMeanr2Cos2PhiColl =0;
+    fMeanr2Sin2PhiColl =0;
+    fMeanr2Cos3PhiColl =0;
+    fMeanr2Sin3PhiColl =0;
+    fMeanr2Cos4PhiColl =0;
+    fMeanr2Sin4PhiColl =0;
+    fMeanr2Cos5PhiColl =0;
+    fMeanr2Sin5PhiColl =0;
+    fMeanr3Cos3PhiColl =0;
+    fMeanr3Sin3PhiColl =0;
+    fMeanr4Cos4PhiColl =0;
+    fMeanr4Sin4PhiColl =0;
+    fMeanr5Cos5PhiColl =0;
+    fMeanr5Sin5PhiColl =0;
+  }
   if(fAN+fBN>0)
-    {
-      fMeanXSystem /= (fAN + fBN);
-      fMeanYSystem /= (fAN + fBN);
-    }
+  {
+    fMeanXSystem /= (fAN + fBN);
+    fMeanYSystem /= (fAN + fBN);
+  }
   else
-    {
-      fMeanXSystem = 0;
-      fMeanYSystem = 0;
-    }
+  {
+    fMeanXSystem = 0;
+    fMeanYSystem = 0;
+  }
   if(fAN>0)
-    {
-      fMeanXA /= fAN;
-      fMeanYA /= fAN;
-    }
+  {
+    fMeanXA /= fAN;
+    fMeanYA /= fAN;
+  }
   else
-    {
-      fMeanXA = 0;
-      fMeanYA = 0;
-    }
-  
+  {
+    fMeanXA = 0;
+    fMeanYA = 0;
+  }
+
   if(fBN>0)
-    {
-      fMeanXB /= fBN;
-      fMeanYB /= fBN;
-    }
+  {
+    fMeanXB /= fBN;
+    fMeanYB /= fBN;
+  }
   else
-    {
-      fMeanXB = 0;
-      fMeanYB = 0;
-    }
-  
+  {
+    fMeanXB = 0;
+    fMeanYB = 0;
+  }
+
   fSx2=fMeanX2-(fMeanXParts*fMeanXParts);
   fSy2=fMeanY2-(fMeanYParts*fMeanYParts);
   fSxy=fMeanXY-fMeanXParts*fMeanYParts;
@@ -738,7 +722,7 @@ Double_t AliGlauberMC::GetTotXSectErr() const
 {
   //total xsection error
   return GetTotXSect()/TMath::Sqrt((Double_t)fEvents) *
-    TMath::Sqrt(Double_t(1.-fEvents/fTotalEvents));
+         TMath::Sqrt(Double_t(1.-fEvents/fTotalEvents));
 }
 
 //______________________________________________________________________________
@@ -751,13 +735,13 @@ TObjArray *AliGlauberMC::GetNucleons()
   TObjArray *allnucleons=new TObjArray(fAN+fBN);
   allnucleons->SetOwner();
   for (Int_t i = 0; i<fAN; i++)
-    {
-      allnucleons->Add(fNucleonsA->UncheckedAt(i));
-    }
+  {
+    allnucleons->Add(fNucleonsA->UncheckedAt(i));
+  }
   for (Int_t i = 0; i<fBN; i++)
-    {
-      allnucleons->Add(fNucleonsB->UncheckedAt(i));
-    }
+  {
+    allnucleons->Add(fNucleonsB->UncheckedAt(i));
+  }
   return allnucleons;
 }
 
@@ -766,13 +750,13 @@ Double_t AliGlauberMC::NegativeBinomialDistribution(Int_t x, Int_t k, Double_t n
 {
   //produce a number distributed acc. neg.bin.dist
   if(k<=0)
-    {
-      cout << "Error, zero or negative K" << endl;
-      return 0;
-    }
+  {
+    cout << "Error, zero or negative K" << endl;
+    return 0;
+  }
   return (TMath::Binomial(x+k-1,x))
-    *TMath::Power(((nmean/Double_t(k))/(1+nmean/Double_t(k))),Double_t(x))
-    *(1/(TMath::Power((1+nmean/Double_t(k)),Double_t(k))));
+         *TMath::Power(((nmean/Double_t(k))/(1+nmean/Double_t(k))),Double_t(x))
+         *(1/(TMath::Power((1+nmean/Double_t(k)),Double_t(k))));
 }
 
 //______________________________________________________________________________
@@ -783,133 +767,166 @@ Int_t AliGlauberMC::NegativeBinomialRandom(Int_t k, Double_t nmean) const
   Double_t array[fMaxPlot];
   array[0] = NegativeBinomialDistribution(0,k,nmean);
   for (Int_t i=1; i<fMaxPlot; i++)
-    {
-      array[i] = NegativeBinomialDistribution(i,k,nmean) + array[i-1];
-    }
+  {
+    array[i] = NegativeBinomialDistribution(i,k,nmean) + array[i-1];
+  }
   Double_t r = gRandom->Uniform(0,1);
   return TMath::BinarySearch(fMaxPlot,array,r)+2;
 
 }
+
+//______________________________________________________________________________
+Int_t AliGlauberMC::NegativeBinomialRandomSV(Double_t k, Double_t nbar) const
+{
+  // negative binomial distribution generator, S. Voloshin, 09-May-2007
+  Double_t sum=0.;
+  Int_t i=0;
+  Double_t ran=gRandom->Rndm();
+  Double_t trm=1./pow(1.+nbar/k,k);
+  if (trm==0.)
+  {
+    cout<<"NBD overflow"<<"  nbar="<<nbar<<"   k="<<k<<endl;
+    return -1;
+  }
+  for(i=0; i<2000 && sum<ran ; i++)
+  {
+    sum += trm;
+    trm *= (k+i)/(i+1.)*(nbar/(nbar+k));
+  }
+  return i-1;
+}
+
 //______________________________________________________________________________
 Int_t AliGlauberMC::DoubleNegativeBinomialRandom( Int_t k,
-                                                  Double_t nmean,
-                                                  Int_t k2,
-                                                  Double_t nmean2,
-                                                  Double_t alpha ) const
+    Double_t nmean,
+    Int_t k2,
+    Double_t nmean2,
+    Double_t alpha ) const
 {
   //return random integer from a Double Negative Binomial Distribution
   static const Int_t fMaxPlot = 1000;
   Double_t array[fMaxPlot];
   array[0] = alpha*NegativeBinomialDistribution(0,k,nmean)+(1-alpha)*NegativeBinomialDistribution(0,k2,nmean2);
   for (Int_t i=1; i<fMaxPlot; i++)
-    {
-      array[i] = alpha*NegativeBinomialDistribution(i,k,nmean)+(1-alpha)*NegativeBinomialDistribution(i,k2,nmean2) + array[i-1];
-    }
+  {
+    array[i] = alpha*NegativeBinomialDistribution(i,k,nmean)+(1-alpha)*NegativeBinomialDistribution(i,k2,nmean2) + array[i-1];
+  }
   Double_t r = gRandom->Uniform(0,1);
   return TMath::BinarySearch(fMaxPlot,array,r)+2;
 }
 
 //______________________________________________________________________________
-void AliGlauberMC::SetdNdEtaParam(Double_t nnp, Double_t x)
-{
-  // set parameters used for calculating multiplicity, see GetdNdEta() for commments
-  fdNdEtaParam[0]=nnp;
-  fdNdEtaParam[1]=x;
-}
-
-//______________________________________________________________________________
-void AliGlauberMC::SetdNdEtaGBWParam(Double_t delta, Double_t lambda, Double_t snn)
-{
-  // set parameters used for calculating multiplicity see GetdNdEtaGBW() for commments
-  fdNdEtaGBWParam[0]=delta;
-  fdNdEtaGBWParam[1]=lambda;
-  fdNdEtaGBWParam[2]=snn;
-}
-
-//______________________________________________________________________________
-void AliGlauberMC::SetdNdEtaNBDParam(Double_t k, Double_t nmean, Double_t beta)
+Double_t AliGlauberMC::GetdNdEta() const
 {
-  // set parameters used for calculating multiplicity see GetdNdEtaNBD() for commments
-  fdNdEtaNBDParam[0]=k;
-  fdNdEtaNBDParam[1]=nmean;
-  fdNdEtaNBDParam[2]=beta;
-}
-
-//______________________________________________________________________________
-void AliGlauberMC::SetdNdEtaTwoNBDParam( Double_t k1, 
-                                         Double_t nmean1, 
-                                         Double_t k2, 
-                                         Double_t nmean2, 
-                                         Double_t alpha,
-                                         Double_t beta)
-{
-  // set parameters used for calculating multiplicity see GetdNdEtaTwoNBD() for commments
-  fdNdEtaTwoNBDParam[0]=k1;
-  fdNdEtaTwoNBDParam[1]=nmean1;
-  fdNdEtaTwoNBDParam[2]=k2;
-  fdNdEtaTwoNBDParam[3]=nmean2;
-  fdNdEtaTwoNBDParam[4]=alpha;
-  fdNdEtaTwoNBDParam[5]=beta;
+  switch (fMultType)
+  {
+    case kSimple:
+      return GetdNdEtaSimple(fdNdEtaParam);
+    case kNBD:
+      return GetdNdEtaNBD(fdNdEtaParam);
+    case kNBDSV:
+      return GetdNdEtaNBDSV(fdNdEtaParam);
+    case kTwoNBD:
+      return GetdNdEtaTwoNBD(fdNdEtaParam);
+    case kGBW:
+      return GetdNdEtaGBW(fdNdEtaParam);
+    default:
+      return 0.0;
+  }
 }
 
 //______________________________________________________________________________
-Double_t AliGlauberMC::GetdNdEta(Double_t nnp, Double_t x) const
+Double_t AliGlauberMC::GetdNdEtaSimple(const Double_t* p) const 
 {
   //Get particle density per unit of rapidity
   //using two component model
   //Parameters: npp, x
+  Double_t nnp = p[0]; //=8.0
+  Double_t x = p[1]; //=0.13
+
   return nnp*((1.-x)*fNpart/2.+x*fNcoll);
 }
 
 //______________________________________________________________________________
-Double_t AliGlauberMC::GetdNdEtaGBW( Double_t delta,
-                                     Double_t lambda,
-                                     Double_t  snn) const
+Double_t AliGlauberMC::GetdNdEtaGBW( const Double_t* p ) const
 {
   //Get particle density per unit of rapidity
   //using the GBW model
   //Parameters: delta, lambda, snn
+  Double_t delta = p[0]; //=0.79
+  Double_t lambda = p[1]; //=0.288
+  Double_t  snn = p[2]; //=30.25
+
   return fNpart*0.47*TMath::Sqrt(TMath::Power(snn,lambda))
-    * TMath::Power(fNpart,(1.-delta)/3./delta);
+         * TMath::Power(fNpart,(1.-delta)/3./delta);
 }
 
 //_______________________________________________________________________________
-Double_t AliGlauberMC::GetdNdEtaNBD ( Int_t k, Double_t nmean, Double_t beta) const
+Double_t AliGlauberMC::GetdNdEtaNBDSV ( const Double_t* p ) const 
+{
+  //Get particle density per unit of rapidity (from Sergei)
+  Double_t npp = p[0];          //=2.49
+  Double_t ratioSgm2Mu = p[1];  //=1.7
+  Double_t xhard = p[2];        //=0.13
+
+  double knb=npp/(ratioSgm2Mu-1.);  // sgm^2/mu = 1+ mu/k
+  double scale = (1.-xhard)*fNpart/2.+xhard*fNcoll; // effectively get number of pp collisions
+  float nch1=-99.;
+  if (knb*scale <1000.)
+  {
+    nch1=(float) NegativeBinomialRandomSV( npp*scale, knb*scale );
+  }
+  else
+  {
+    nch1=(float) NegativeBinomialRandomSV( npp*scale/2., knb*scale/2. ) +
+         (float) NegativeBinomialRandomSV( npp*scale/2., knb*scale/2. );
+  }
+  return nch1;
+}
+
+//_______________________________________________________________________________
+Double_t AliGlauberMC::GetdNdEtaNBD ( const Double_t* p ) const
 {
   //Get particle density per unit of rapidity
   //using a aandomized number from a negative binomial distrubution
-  //Parameters:   k  = related to distribition width
-  //              nmean = mean of distribution
-  //              beta = set contribution of participants / binary collisions to multiplicity
+  //Parameters:   k  = related to distribition width=3
+  //              nmean = mean of distribution=4
+  //              beta = set contribution of participants / binary collisions to multiplicity=0.13
+  Int_t k = TMath::Nint(p[0]);
+  Double_t nmean = p[1];
+  Double_t beta = p[2];
+
   Double_t mulnp=0.0;
   for(Int_t i = 0; i<fNpart; i++)
-    {
-      mulnp+=NegativeBinomialRandom(k,nmean);
-    }
+  {
+    mulnp+=NegativeBinomialRandom(k,nmean);
+  }
   Double_t mulnb=0.0;
   for(Int_t i = 0; i<fNcoll; i++)
-    {
-      mulnb+=NegativeBinomialRandom(k,nmean);
-    }
+  {
+    mulnb+=NegativeBinomialRandom(k,nmean);
+  }
   return (1-beta)*mulnp/2+beta*mulnb;
 }
 
 //______________________________________________________________________________
-Double_t AliGlauberMC::GetdNdEtaTwoNBD ( Int_t k1,
-                                         Double_t nmean1,
-                                         Int_t k2,
-                                         Double_t nmean2,
-                                         Double_t alpha,
-                                         Double_t beta ) const
+Double_t AliGlauberMC::GetdNdEtaTwoNBD ( const Double_t* p ) const
 {
   //Get particle density per unit of rapidity
   //using random numbers from two negative binomial distributions
-  //Parameters:   k1 = related to distribition width of distribution 1
-  //              nmean1 = mean of distribution 1
-  //              k2 = related to distribition width of distribution 2
-  //              nmean2 = mean of distribution 2
-  //              alpha = set contributions of distrubitin 1 / distribution 2
-  //              beta = set contribution of participants / binary collisions to multiplicity
+  //Parameters:   k1 = related to distribition width of distribution 1=3
+  //              nmean1 = mean of distribution 1=4
+  //              k2 = related to distribition width of distribution 2=2
+  //              nmean2 = mean of distribution 2=11
+  //              alpha = set contributions of distrubitin 1 / distribution 2=0.4
+  //              beta = set contribution of participants / binary collisions to multiplicity =0.13
+  Int_t k1 = TMath::Nint(p[0]);
+  Double_t nmean1 = p[1];
+  Int_t k2 = TMath::Nint(p[2]);
+  Double_t nmean2 = p[3];
+  Double_t alpha = p[4];
+  Double_t beta = p[6];
+
   Double_t mulnp=0.0;
   for(Int_t i = 0; i<fNpart; i++)
   {
@@ -1019,7 +1036,7 @@ Double_t AliGlauberMC::GetEpsilon2Coll() const
 //______________________________________________________________________________
 Double_t AliGlauberMC::GetEpsilon3Coll() const
 {
-  //get epsilon3 of binary collisions 
+  //get epsilon3 of binary collisions
   if (fNcoll<2) return 0.0;
   return (TMath::Sqrt(fMeanr2Cos3PhiColl*fMeanr2Cos3PhiColl+fMeanr2Sin3PhiColl*fMeanr2Sin3PhiColl)/fMeanr2Coll);
   //return (TMath::Sqrt(fMeanr3Cos3PhiColl*fMeanr3Cos3PhiColl+fMeanr3Sin3PhiColl*fMeanr3Sin3PhiColl)/fMeanr3Coll);
@@ -1028,7 +1045,7 @@ Double_t AliGlauberMC::GetEpsilon3Coll() const
 //______________________________________________________________________________
 Double_t AliGlauberMC::GetEpsilon4Coll() const
 {
-   //get epsilon4 of binary collisions
+  //get epsilon4 of binary collisions
   if (fNcoll<2) return 0.0;
   return (TMath::Sqrt(fMeanr2Cos4PhiColl*fMeanr2Cos4PhiColl+fMeanr2Sin4PhiColl*fMeanr2Sin4PhiColl)/fMeanr2Coll);
   //return (TMath::Sqrt(fMeanr4Cos4PhiColl*fMeanr4Cos4PhiColl+fMeanr4Sin4PhiColl*fMeanr4Sin4PhiColl)/fMeanr4Coll);
@@ -1053,84 +1070,78 @@ void AliGlauberMC::Run(Int_t nevents)
   if (fnt == 0)
   {
     fnt = new TNtuple(name,title,
-                      "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:Epsl2:Epsl3:Epsl4:Epsl5:E2Coll:E3Coll:E4Coll:E5Coll");
+                      "Npart:Ncoll:B:MeanX:MeanY:MeanX2:MeanY2:MeanXY:VarX:VarY:"
+                      "VarXY:MeanXSystem:MeanYSystem:MeanXA:MeanYA:MeanXB:MeanYB:VarE:VarEColl:VarEPart:"
+                      "VarEPartColl:dNdEta1:dNdEta2:dNdEta:xsect:tAA:Epsl2:Epsl3:Epsl4:Epsl5:"
+                      "E2Coll:E3Coll:E4Coll:E5Coll" );
     fnt->SetDirectory(0);
   }
   Int_t q = 0;
   Int_t u = 0;
   for (Int_t i = 0; i<nevents; i++)
+  {
+
+    if(!NextEvent())
+    {
+      u++;
+      continue;
+    }
+
+    q++;
+    //Float_t v[27];
+    Float_t v[34];
+    v[0]  = GetNpart();
+    v[1]  = GetNcoll();
+    v[2]  = fBMC;
+    v[3]  = fMeanXParts;
+    v[4]  = fMeanYParts;
+    v[5]  = fMeanX2;
+    v[6]  = fMeanY2;
+    v[7]  = fMeanXY;
+    v[8]  = fSx2;
+    v[9]  = fSy2;
+    v[10] = fSxy;
+    v[11] = fMeanXSystem;
+    v[12] = fMeanYSystem;
+    v[13] = fMeanXA;
+    v[14] = fMeanYA;
+    v[15] = fMeanXB;
+    v[16] = fMeanYB;
+    v[17] = GetEccentricity();
+    v[18] = GetEccentricityColl();
+    v[19] = GetEccentricityPart();
+    v[20] = GetEccentricityPartColl();
+    if (fDoPartProd)
+    {
+      v[21] = GetdNdEta();
+      v[22] = GetdNdEta();
+      v[23] = v[21]+v[22];
+    }
+    else
     {
-      
-      if(!NextEvent())
-       {
-         u++;
-         continue;
-       }
-      
-      q++;
-      //Float_t v[27];
-      Float_t v[35];
-      v[0]  = GetNpart();
-      v[1]  = GetNcoll();
-      v[2]  = fBMC;
-      v[3]  = fMeanXParts;
-      v[4]  = fMeanYParts;
-      v[5]  = fMeanX2;
-      v[6]  = fMeanY2;
-      v[7]  = fMeanXY;
-      v[8]  = fSx2;
-      v[9]  = fSy2;
-      v[10] = fSxy;
-      v[11] = fMeanXSystem;
-      v[12] = fMeanYSystem;
-      v[13] = fMeanXA;
-      v[14] = fMeanYA;
-      v[15] = fMeanXB;
-      v[16] = fMeanYB;
-      v[17] = GetEccentricity();
-      v[18] = GetEccentricityColl();
-      v[19] = GetEccentricityPart();
-      v[20] = GetEccentricityPartColl();
-      if (fDoPartProd) 
-       {
-         v[21] = GetdNdEta( fdNdEtaParam[0],fdNdEtaParam[1] );
-         v[22] = GetdNdEtaGBW( fdNdEtaGBWParam[0],fdNdEtaGBWParam[1],fdNdEtaGBWParam[2] );
-         v[23] = GetdNdEtaNBD( TMath::Nint(fdNdEtaNBDParam[0]),
-                               fdNdEtaNBDParam[1],
-                               fdNdEtaNBDParam[2] );
-         v[24] = GetdNdEtaTwoNBD( TMath::Nint(fdNdEtaTwoNBDParam[0]),
-                                  fdNdEtaTwoNBDParam[1],
-                                  TMath::Nint(fdNdEtaTwoNBDParam[2]),
-                                  fdNdEtaTwoNBDParam[3],
-                                  fdNdEtaTwoNBDParam[4],
-                                  fdNdEtaTwoNBDParam[5] );
-       }
-      else 
-       {
-         v[21] = 0;
-         v[22] = 0;
-         v[23] = 0;
-         v[24] = 0;
-       }
-      v[25]=fXSect;
-      
-      Float_t mytAA=-999;
-      if (GetNcoll()>0) mytAA=GetNcoll()/fXSect;
-      v[26]=mytAA;    
-      //_____________epsilon2,3,4,4_______
-      v[27] = GetEpsilon2Part();
-      v[28] = GetEpsilon3Part();
-      v[29] = GetEpsilon4Part();
-      v[30] = GetEpsilon5Part();
-      v[31] = GetEpsilon2Coll();
-      v[32] = GetEpsilon3Coll();
-      v[33] = GetEpsilon4Coll();
-      v[34] = GetEpsilon5Coll();
-      //always at the end
-      fnt->Fill(v);  
-
-      if ((i%100)==0) std::cout << "Generating Event # " << i << "... \r" << flush;
+      v[21] = 0;
+      v[22] = 0;
+      v[23] = 0;
     }
+    v[24]=fXSect;
+
+    Float_t mytAA=-999;
+    if (GetNcoll()>0) mytAA=GetNcoll()/fXSect;
+    v[25]=mytAA;
+    //_____________epsilon2,3,4,4_______
+    v[26] = GetEpsilon2Part();
+    v[27] = GetEpsilon3Part();
+    v[28] = GetEpsilon4Part();
+    v[29] = GetEpsilon5Part();
+    v[30] = GetEpsilon2Coll();
+    v[31] = GetEpsilon3Coll();
+    v[32] = GetEpsilon4Coll();
+    v[33] = GetEpsilon5Coll();
+    //always at the end
+    fnt->Fill(v);
+
+    if ((i%100)==0) std::cout << "Generating Event # " << i << "... \r" << flush;
+  }
   std::cout << "Generating Event # " << nevents << "... \r" << endl << "Done! Succesfull events:  " << q << "  discarded events:  " << u <<"."<< endl;
 }
 
@@ -1140,8 +1151,8 @@ void AliGlauberMC::RunAndSaveNtuple( Int_t n,
                                      const Option_t *sysB,
                                      Double_t signn,
                                      Double_t mind,
-                                    Double_t r,
-                                    Double_t a,
+                                     Double_t r,
+                                     Double_t a,
                                      const char *fname)
 {
   //example run
@@ -1174,36 +1185,36 @@ void AliGlauberMC::RunAndSaveNucleons( Int_t n,
     out=new TFile(fname,"recreate",fname,9);
 
   for(Int_t ievent=0; ievent<n; ievent++)
+  {
+
+    //get an event with at least one collision
+    while(!mcg.NextEvent()) {}
+
+    //access, save and (if wanted) print out nucleons
+    TObjArray* nucleons=mcg.GetNucleons();
+    if(!nucleons) continue;
+    if(out)
+      nucleons->Write(Form("nucleonarray%d",ievent),TObject::kSingleKey);
+
+    if(verbose)
     {
-      
-      //get an event with at least one collision
-      while(!mcg.NextEvent()) {}
-      
-      //access, save and (if wanted) print out nucleons
-      TObjArray* nucleons=mcg.GetNucleons();
-      if(!nucleons) continue;
-      if(out)
-       nucleons->Write(Form("nucleonarray%d",ievent),TObject::kSingleKey);
-      
-      if(verbose)
-       {
-         cout<<endl<<endl<<"EVENT NO: "<<ievent<<endl;
-         cout<<"B = "<<mcg.GetB()<<"  Npart = "<<mcg.GetNpart()<<endl<<endl;
-         printf("Nucleus\t X\t Y\t Z\tNcoll\n");
-         Int_t nNucls=nucleons->GetEntries();
-         for(Int_t iNucl=0; iNucl<nNucls; iNucl++)
-           {
-             AliGlauberNucleon *nucl=(AliGlauberNucleon *)nucleons->UncheckedAt(iNucl);
-             Char_t nucleus='A';
-             if(nucl->IsInNucleusB()) nucleus='B';
-             Double_t x=nucl->GetX();
-             Double_t y=nucl->GetY();
-             Double_t z=nucl->GetZ();
-             Int_t ncoll=nucl->GetNColl();
-             printf("   %c\t%2.2f\t%2.2f\t%2.2f\t%3d\n",nucleus,x,y,z,ncoll);
-           }
-       }
+      cout<<endl<<endl<<"EVENT NO: "<<ievent<<endl;
+      cout<<"B = "<<mcg.GetB()<<"  Npart = "<<mcg.GetNpart()<<endl<<endl;
+      printf("Nucleus\t X\t Y\t Z\tNcoll\n");
+      Int_t nNucls=nucleons->GetEntries();
+      for(Int_t iNucl=0; iNucl<nNucls; iNucl++)
+      {
+        AliGlauberNucleon *nucl=(AliGlauberNucleon *)nucleons->UncheckedAt(iNucl);
+        Char_t nucleus='A';
+        if(nucl->IsInNucleusB()) nucleus='B';
+        Double_t x=nucl->GetX();
+        Double_t y=nucl->GetY();
+        Double_t z=nucl->GetZ();
+        Int_t ncoll=nucl->GetNColl();
+        printf("   %c\t%2.2f\t%2.2f\t%2.2f\t%3d\n",nucleus,x,y,z,ncoll);
+      }
     }
+  }
   if(out) delete out;
 }
 
@@ -1211,5 +1222,6 @@ void AliGlauberMC::RunAndSaveNucleons( Int_t n,
 void AliGlauberMC::Reset()
 {
   //delete the ntuple
-  delete fnt; fnt=NULL;
+  delete fnt;
+  fnt=NULL;
 }
index a372bbe..f5f2522 100644 (file)
 ////////////////////////////////////////////////////////////////////////////////
 
 #include "AliGlauberNucleus.h"
-
 #include <Riostream.h>
-class TNamed;
+#include <TNamed.h>
+
 class TObjArray;
 class TNtuple;
-class TArray;
 
 class AliGlauberMC : public TNamed {
 public:
+   enum EdNdEtaType { kSimple,
+                      kNBD,
+                      kNBDSV,
+                      kTwoNBD,
+                      kGBW,
+                      kNone };
+
    AliGlauberMC(Option_t* NA = "Pb", Option_t* NB = "Pb", Double_t xsect = 64);
    virtual     ~AliGlauberMC();
    AliGlauberMC(const AliGlauberMC& in);
@@ -34,20 +40,12 @@ public:
    Bool_t       CalcEvent(Double_t bgen);
 
    //various ways to calculate multiplicity
-   Double_t     GetdNdEta(    Double_t nnp=8.0,
-                              Double_t x=0.13 ) const;
-   Double_t     GetdNdEtaGBW( Double_t delta=0.79,
-                              Double_t lambda=0.288,
-                              Double_t snn=30.25 ) const;
-   Double_t    GetdNdEtaNBD(     Int_t k=3,
-                              Double_t nmean = 4,
-                              Double_t beta = 0.13 ) const;
-   Double_t    GetdNdEtaTwoNBD(  Int_t k1=3,
-                              Double_t nmean1=4,
-                              Int_t k2=2,
-                              Double_t nmean2=11,
-                              Double_t alpha=0.4,
-                              Double_t beta=0.13 ) const;
+   Double_t GetdNdEta() const;
+   Double_t GetdNdEtaSimple( const Double_t* param ) const;
+   Double_t GetdNdEtaGBW(    const Double_t* param ) const;
+   Double_t    GetdNdEtaNBD(    const Double_t* param ) const;
+   Double_t    GetdNdEtaNBDSV(  const Double_t* param ) const;
+   Double_t    GetdNdEtaTwoNBD( const Double_t* param ) const;
 
    Double_t     GetEccentricity()    const;
    Double_t     GetEccentricityColl()      const;
@@ -74,13 +72,12 @@ public:
    void         Reset();
    static Double_t     NegativeBinomialDistribution(Int_t x, Int_t k, Double_t nmean);
    Int_t  NegativeBinomialRandom(Int_t k, Double_t nmean) const;
+   Int_t  NegativeBinomialRandomSV(Double_t nbar, Double_t k) const;
    Int_t  DoubleNegativeBinomialRandom(Int_t k1, Double_t nmean1, Int_t k2, Double_t nmean2, Double_t alpha) const;
+   Double_t* GetdNdEtaParam() {return fdNdEtaParam;}
+   void      SetdNdEtaType(EdNdEtaType method) {fMultType=method;}
    void   SetBmin(Double_t bmin)      {fBMin = bmin;}
    void   SetBmax(Double_t bmax)      {fBMax = bmax;}
-   void   SetdNdEtaParam( Double_t nnp = 8., Double_t x = 0.13);
-   void   SetdNdEtaGBWParam( Double_t delta = 0.79, Double_t lambda = 0.288, Double_t snn = 30.25);
-   void   SetdNdEtaNBDParam(Double_t k=3, Double_t nmean=4, Double_t beta=0.13);
-   void   SetdNdEtaTwoNBDParam(Double_t alpha = 0.4, Double_t k1 = 3, Double_t nmean1 = 4., Double_t k2 = 2., Double_t nmean2 = 11., Double_t beta=0.13);
    void   SetMinDistance(Double_t d)  {fANucleus.SetMinDist(d); fBNucleus.SetMinDist(d);}
    void   SetDoPartProduction(Bool_t b) { fDoPartProd = b; }
    void   Setr(Double_t r)  {fANucleus.SetR(r); fBNucleus.SetR(r);}
@@ -133,10 +130,8 @@ private:
    Int_t        fTotalEvents;    //All events within selected impact parameter range
    Double_t     fBMin;           //Minimum impact parameter to be generated
    Double_t     fBMax;           //Maximum impact parameter to be generated
-   Double_t    fdNdEtaParam[2];           //Parameters: nnp, x
-   Double_t     fdNdEtaGBWParam[3];  //Parameters: delta, lambda, snn
-   Double_t     fdNdEtaNBDParam[3];       //Parameters:  k, nmean, beta
-   Double_t     fdNdEtaTwoNBDParam[6];    //Parameters: k1, nmean1, k2, nmean2, alpha, beta
+   Double_t        fdNdEtaParam[10];//Parameters for multiplicity calculation: meaning depends on method selection
+   EdNdEtaType fMultType;//mutliplicity method selection  
    Int_t        fMaxNpartFound;  //Largest value of Npart obtained
    Int_t        fNpart;          //Number of wounded (participating) nucleons in current event
    Int_t        fNcoll;          //Number of binary collisions in current event