//
// origin: PHOBOS experiment
// alification: Mikolaj Krzewicki, Nikhef, mikolaj.krzewicki@cern.ch
-//
+// update: You Zhou, Nikhef, yzhou@nikhef.nl
////////////////////////////////////////////////////////////////////////////////
#include <Riostream.h>
#include <TNtuple.h>
#include <TFile.h>
#include <TTree.h>
-#include <TF1.h>
-#include <TH1F.h>
-#include <TArray.h>
#include "AliGlauberNucleon.h"
#include "AliGlauberNucleus.h"
//______________________________________________________________________________
AliGlauberMC::AliGlauberMC(Option_t* NA, Option_t* NB, Double_t xsect) :
+ TNamed(),
fANucleus(NA),
fBNucleus(NB),
- fXSect(0),
+ fXSect(xsect),
fNucleonsA(0),
fNucleonsB(0),
+ fAN(0),
+ fQAN(0),
+ fBN(0),
+ fQBN(0),
fnt(0),
fMeanX2(0),
fMeanY2(0),
fMeanXY(0),
+ fMeanX2Parts(0),
+ fMeanY2Parts(0),
+ fMeanXYParts(0),
fMeanXParts(0),
fMeanYParts(0),
+ fMeanOXParts(0),
+ fMeanOYParts(0),
fMeanXColl(0),
fMeanYColl(0),
+ fMeanOXColl(0),
+ fMeanOYColl(0),
fMeanX2Coll(0),
- fMeanY2Coll(0),
+ fMeanY2Coll(0),
fMeanXYColl(0),
+ fMeanXCom(0),
+ fMeanYCom(0),
+ fMeanOXCom(0),
+ fMeanOYCom(0),
+ fMeanX2Com(0),
+ fMeanY2Com(0),
+ fMeanXYCom(0),
fMeanXSystem(0),
fMeanYSystem(0),
- fMeanX_A(0),
- fMeanY_A(0),
- fMeanX_B(0),
- fMeanY_B(0),
- fB_MC(0),
+ fMeanXA(0),
+ fMeanYA(0),
+ fMeanXB(0),
+ fMeanYB(0),
+ fMeanOXA(0),
+ fMeanOYA(0),
+ fMeanOXB(0),
+ fMeanOYB(0),
+ fBMC(0),
fEvents(0),
fTotalEvents(0),
- fBMin(0),
- fBMax(0),
+ fBMin(0.),
+ fBMax(20.),
+ fMultType(kNBDSV),
fMaxNpartFound(0),
+ fONpart(0),
+ fONcoll(0),
+ fONcom(0),
fNpart(0),
fNcoll(0),
- fSx2(0),
- fSy2(0),
- fSxy(0),
+ fNcom(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),
+ fMeanr2Coll(0),
+ fMeanr3Coll(0),
+ fMeanr4Coll(0),
+ fMeanr5Coll(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),
+ fMeanr2Com(0),
+ fMeanr3Com(0),
+ fMeanr4Com(0),
+ fMeanr5Com(0),
+ fMeanr2Cos2PhiCom(0),
+ fMeanr2Sin2PhiCom(0),
+ fMeanr2Cos3PhiCom(0),
+ fMeanr2Sin3PhiCom(0),
+ fMeanr2Cos4PhiCom(0),
+ fMeanr2Sin4PhiCom(0),
+ fMeanr2Cos5PhiCom(0),
+ fMeanr2Sin5PhiCom(0),
+ fMeanr3Cos3PhiCom(0),
+ fMeanr3Sin3PhiCom(0),
+ fMeanr4Cos4PhiCom(0),
+ fMeanr4Sin4PhiCom(0),
+ fMeanr5Cos5PhiCom(0),
+ fMeanr5Sin5PhiCom(0),
+ fSx2Parts(0.),
+ fSy2Parts(0.),
+ fSxyParts(0.),
+ fSx2Coll(0.),
+ fSy2Coll(0.),
+ fSxyColl(0.),
+ fSx2Com(0.),
+ fSy2Com(0.),
+ fSxyCom(0.),
fX(0.13),
- fNpp(8.)
-{
- fBMin = 0;
- fBMax = 20;
- fXSect = xsect;
- 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] = 0.4;
- fdNdEtaTwoNBDParam[1] = 3.0;
- fdNdEtaTwoNBDParam[2] = 4.0;
- fdNdEtaTwoNBDParam[3] = 2.0;
- fdNdEtaTwoNBDParam[4] = 11.0;
- fdNdEtaTwoNBDParam[5] = 0.13;
-
- TString name(Form("Glauber_%s_%s",fANucleus.GetName(),fBNucleus.GetName()));
- TString title(Form("Glauber %s+%s Version",fANucleus.GetName(),fBNucleus.GetName()));
- SetName(name);
- SetTitle(title);
+ fNpp(8.),
+ fDoPartProd(kFALSE)
+{
+ //ctor
+ 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()));
+}
+
+//______________________________________________________________________________
+AliGlauberMC::~AliGlauberMC()
+{
+ //dtor
+ delete fnt;
+}
+
+//______________________________________________________________________________
+AliGlauberMC::AliGlauberMC(const AliGlauberMC& in):
+ TNamed(in),
+ fANucleus(in.fANucleus),
+ fBNucleus(in.fBNucleus),
+ fXSect(in.fXSect),
+ fNucleonsA(in.fNucleonsA),
+ fNucleonsB(in.fNucleonsB),
+ fAN(in.fAN),
+ fQAN(in.fQAN),
+ fBN(in.fBN),
+ fQBN(in.fQBN),
+ fnt(in.fnt),
+ fMeanX2(in.fMeanX2),
+ fMeanY2(in.fMeanY2),
+ fMeanXY(in.fMeanXY),
+ fMeanX2Parts(in.fMeanX2Parts),
+ fMeanY2Parts(in.fMeanY2Parts),
+ fMeanXYParts(in.fMeanXYParts),
+ fMeanXParts(in.fMeanXParts),
+ fMeanYParts(in.fMeanYParts),
+ fMeanOXParts(in.fMeanOXParts),
+ fMeanOYParts(in.fMeanOYParts),
+ fMeanXColl(in.fMeanXColl),
+ fMeanYColl(in.fMeanYColl),
+ fMeanOXColl(in.fMeanOXColl),
+ fMeanOYColl(in.fMeanOYColl),
+ fMeanX2Coll(in.fMeanX2Coll),
+ fMeanY2Coll(in.fMeanY2Coll),
+ fMeanXYColl(in.fMeanXYColl),
+ fMeanXCom(in.fMeanXCom),
+ fMeanYCom(in.fMeanYCom),
+ fMeanOXCom(in.fMeanOXCom),
+ fMeanOYCom(in.fMeanOYCom),
+ fMeanX2Com(in.fMeanX2Com),
+ fMeanY2Com(in.fMeanY2Com),
+ fMeanXYCom(in.fMeanXYCom),
+ fMeanXSystem(in.fMeanXSystem),
+ fMeanYSystem(in.fMeanYSystem),
+ fMeanXA(in.fMeanXA),
+ fMeanYA(in.fMeanYA),
+ fMeanXB(in.fMeanXB),
+ fMeanYB(in.fMeanYB),
+ fMeanOXA(in.fMeanOXA),
+ fMeanOYA(in.fMeanOYA),
+ fMeanOXB(in.fMeanOXB),
+ fMeanOYB(in.fMeanOYB),
+ fBMC(in.fBMC),
+ fEvents(in.fEvents),
+ fTotalEvents(in.fTotalEvents),
+ fBMin(in.fBMin),
+ fBMax(in.fBMax),
+ fMultType(in.fMultType),
+ fMaxNpartFound(in.fMaxNpartFound),
+ fONpart(in.fONpart),
+ fONcoll(in.fONcoll),
+ fONcom(in.fONcom),
+ fNpart(in.fNpart),
+ fNcoll(in.fNcoll),
+ fNcom(in.fNcom),
+ fMeanr2(in.fMeanr2),
+ fMeanr3(in.fMeanr3),
+ fMeanr4(in.fMeanr4),
+ fMeanr5(in.fMeanr5),
+ fMeanr2Cos2Phi(in.fMeanr2Cos2Phi),
+ fMeanr2Sin2Phi(in.fMeanr2Sin2Phi),
+ fMeanr2Cos3Phi(in.fMeanr2Cos3Phi),
+ fMeanr2Sin3Phi(in.fMeanr2Sin3Phi),
+ fMeanr2Cos4Phi(in.fMeanr2Cos4Phi),
+ fMeanr2Sin4Phi(in.fMeanr2Sin4Phi),
+ fMeanr2Cos5Phi(in.fMeanr2Cos5Phi),
+ fMeanr2Sin5Phi(in.fMeanr2Sin5Phi),
+ fMeanr3Cos3Phi(in.fMeanr3Cos3Phi),
+ fMeanr3Sin3Phi(in.fMeanr3Sin3Phi),
+ fMeanr4Cos4Phi(in.fMeanr4Cos4Phi),
+ fMeanr4Sin4Phi(in.fMeanr4Sin4Phi),
+ fMeanr5Cos5Phi(in.fMeanr5Cos5Phi),
+ fMeanr5Sin5Phi(in.fMeanr5Sin5Phi),
+ fMeanr2Coll(in.fMeanr2Coll),
+ fMeanr3Coll(in.fMeanr3Coll),
+ fMeanr4Coll(in.fMeanr4Coll),
+ fMeanr5Coll(in.fMeanr5Coll),
+ fMeanr2Cos2PhiColl(in.fMeanr2Cos2PhiColl),
+ fMeanr2Sin2PhiColl(in.fMeanr2Sin2PhiColl),
+ fMeanr2Cos3PhiColl(in.fMeanr2Cos3PhiColl),
+ fMeanr2Sin3PhiColl(in.fMeanr2Sin3PhiColl),
+ fMeanr2Cos4PhiColl(in.fMeanr2Cos4PhiColl),
+ fMeanr2Sin4PhiColl(in.fMeanr2Sin4PhiColl),
+ fMeanr2Cos5PhiColl(in.fMeanr2Cos5PhiColl),
+ fMeanr2Sin5PhiColl(in.fMeanr2Sin5PhiColl),
+ fMeanr3Cos3PhiColl(in.fMeanr3Cos3PhiColl),
+ fMeanr3Sin3PhiColl(in.fMeanr3Sin3PhiColl),
+ fMeanr4Cos4PhiColl(in.fMeanr4Cos4PhiColl),
+ fMeanr4Sin4PhiColl(in.fMeanr4Sin4PhiColl),
+ fMeanr5Cos5PhiColl(in.fMeanr5Cos5PhiColl),
+ fMeanr5Sin5PhiColl(in.fMeanr5Sin5PhiColl),
+ fMeanr2Com(in.fMeanr2Com),
+ fMeanr3Com(in.fMeanr3Com),
+ fMeanr4Com(in.fMeanr4Com),
+ fMeanr5Com(in.fMeanr5Com),
+ fMeanr2Cos2PhiCom(in.fMeanr2Cos2PhiCom),
+ fMeanr2Sin2PhiCom(in.fMeanr2Sin2PhiCom),
+ fMeanr2Cos3PhiCom(in.fMeanr2Cos3PhiCom),
+ fMeanr2Sin3PhiCom(in.fMeanr2Sin3PhiCom),
+ fMeanr2Cos4PhiCom(in.fMeanr2Cos4PhiCom),
+ fMeanr2Sin4PhiCom(in.fMeanr2Sin4PhiCom),
+ fMeanr2Cos5PhiCom(in.fMeanr2Cos5PhiCom),
+ fMeanr2Sin5PhiCom(in.fMeanr2Sin5PhiCom),
+ fMeanr3Cos3PhiCom(in.fMeanr3Cos3PhiCom),
+ fMeanr3Sin3PhiCom(in.fMeanr3Sin3PhiCom),
+ fMeanr4Cos4PhiCom(in.fMeanr4Cos4PhiCom),
+ fMeanr4Sin4PhiCom(in.fMeanr4Sin4PhiCom),
+ fMeanr5Cos5PhiCom(in.fMeanr5Cos5PhiCom),
+ fMeanr5Sin5PhiCom(in.fMeanr5Sin5PhiCom),
+ fSx2Parts(in.fSx2Parts),
+ fSy2Parts(in.fSy2Parts),
+ fSxyParts(in.fSxyParts),
+ fSx2Coll(in.fSx2Coll),
+ fSy2Coll(in.fSy2Coll),
+ fSxyColl(in.fSxyColl),
+ fSx2Com(in.fSx2Com),
+ fSy2Com(in.fSy2Com),
+ fSxyCom(in.fSxyCom),
+ fX(in.fX),
+ fNpp(in.fNpp),
+ fDoPartProd(kFALSE)
+{
+ //copy ctor
+ memcpy(fdNdEtaParam,in.fdNdEtaParam,sizeof(fdNdEtaParam));
+}
+
+//______________________________________________________________________________
+AliGlauberMC& AliGlauberMC::operator=(const AliGlauberMC& in)
+{
+ //assignment
+ if (&in==this) return *this;
+ fANucleus=in.fANucleus;
+ fBNucleus=in.fBNucleus;
+ fXSect=in.fXSect;
+ fNucleonsA=in.fNucleonsA;
+ fNucleonsB=in.fNucleonsB;
+ fAN=in.fAN;
+ fQAN=in.fQAN;
+ fBN=in.fBN;
+ fQBN=in.fQBN;
+ fnt=in.fnt;
+ fMeanX2=in.fMeanX2;
+ fMeanY2=in.fMeanY2;
+ fMeanXY=in.fMeanXY;
+ fMeanr2=in.fMeanr2;
+ fMeanr3=in.fMeanr3;
+ fMeanr4=in.fMeanr4;
+ fMeanr5=in.fMeanr5;
+ fMeanXParts=in.fMeanXParts;
+ fMeanYParts=in.fMeanYParts;
+ fMeanOXParts=in.fMeanOXParts;
+ fMeanOYParts=in.fMeanOYParts;
+ fMeanXColl=in.fMeanXColl;
+ fMeanYColl=in.fMeanYColl;
+ fMeanOXColl=in.fMeanOXColl;
+ fMeanOYColl=in.fMeanOYColl;
+ fMeanX2Coll=in.fMeanX2Coll;
+ fMeanY2Coll=in.fMeanY2Coll;
+ fMeanXYColl=in.fMeanXYColl;
+ fMeanr2Coll=in.fMeanr2Coll;
+ fMeanr3Coll=in.fMeanr3Coll;
+ fMeanr4Coll=in.fMeanr4Coll;
+ fMeanr5Coll=in.fMeanr5Coll;
+ fMeanXCom=in.fMeanXColl;
+ fMeanYCom=in.fMeanYColl;
+ fMeanOXCom=in.fMeanOXCom;
+ fMeanOYCom=in.fMeanOYCom;
+ fMeanX2Com=in.fMeanX2Com;
+ fMeanY2Com=in.fMeanY2Com;
+ fMeanXYCom=in.fMeanXYCom;
+ fMeanr2Com=in.fMeanr2Com;
+ fMeanr3Com=in.fMeanr3Com;
+ fMeanr4Com=in.fMeanr4Com;
+ fMeanr5Com=in.fMeanr5Com;
+ fMeanXSystem=in.fMeanXSystem;
+ fMeanYSystem=in.fMeanYSystem;
+ fMeanXA=in.fMeanXA;
+ fMeanYA=in.fMeanYA;
+ fMeanXB=in.fMeanXB;
+ fMeanYB=in.fMeanYB;
+ fMeanOXA=in.fMeanOXA;
+ fMeanOYA=in.fMeanOYA;
+ fMeanOXB=in.fMeanOXB;
+ fMeanOYB=in.fMeanOYB;
+ fBMC=in.fBMC;
+ fEvents=in.fEvents;
+ fTotalEvents=in.fTotalEvents;
+ fBMin=in.fBMin;
+ fBMax=in.fBMax;
+ fMultType=in.fMultType,
+ memcpy(fdNdEtaParam,in.fdNdEtaParam,sizeof(fdNdEtaParam));
+ fMaxNpartFound=in.fMaxNpartFound;
+ fNpart=in.fNpart;
+ fNcoll=in.fNcoll;
+ fNcom=in.fNcom;
+ fONpart=in.fONpart;
+ fONcoll=in.fONcoll;
+ fONcom=in.fONcom;
+ fMeanr2Cos2Phi=in.fMeanr2Cos2Phi;
+ fMeanr2Sin2Phi=in.fMeanr2Sin2Phi;
+ fMeanr2Cos3Phi=in.fMeanr2Cos3Phi;
+ fMeanr2Sin3Phi=in.fMeanr2Sin3Phi;
+ fMeanr2Cos4Phi=in.fMeanr2Cos4Phi;
+ fMeanr2Sin4Phi=in.fMeanr2Sin4Phi;
+ fMeanr2Cos5Phi=in.fMeanr2Cos5Phi;
+ fMeanr2Sin5Phi=in.fMeanr2Sin5Phi;
+ fMeanr3Cos3Phi=in.fMeanr3Cos3Phi;
+ fMeanr3Sin3Phi=in.fMeanr3Sin3Phi;
+ fMeanr4Cos4Phi=in.fMeanr4Cos4Phi;
+ fMeanr4Sin4Phi=in.fMeanr4Sin4Phi;
+ fMeanr5Cos5Phi=in.fMeanr5Cos5Phi;
+ fMeanr5Sin5Phi=in.fMeanr5Sin5Phi;
+ fMeanr2Cos2PhiColl=in.fMeanr2Cos2PhiColl;
+ fMeanr2Sin2PhiColl=in.fMeanr2Sin2PhiColl;
+ fMeanr2Cos3PhiColl=in.fMeanr2Cos3PhiColl;
+ fMeanr2Sin3PhiColl=in.fMeanr2Sin3PhiColl;
+ fMeanr2Cos4PhiColl=in.fMeanr2Cos4PhiColl;
+ fMeanr2Sin4PhiColl=in.fMeanr2Sin4PhiColl;
+ fMeanr2Cos5PhiColl=in.fMeanr2Cos5PhiColl;
+ fMeanr2Sin5PhiColl=in.fMeanr2Sin5PhiColl;
+ fMeanr3Cos3PhiColl=in.fMeanr3Cos3PhiColl;
+ fMeanr3Sin3PhiColl=in.fMeanr3Sin3PhiColl;
+ fMeanr4Cos4PhiColl=in.fMeanr4Cos4PhiColl;
+ fMeanr4Sin4PhiColl=in.fMeanr4Sin4PhiColl;
+ fMeanr5Cos5PhiColl=in.fMeanr5Cos5PhiColl;
+ fMeanr5Sin5PhiColl=in.fMeanr5Sin5PhiColl;
+ fMeanr2Cos2PhiCom=in.fMeanr2Cos2PhiCom;
+ fMeanr2Sin2PhiCom=in.fMeanr2Sin2PhiCom;
+ fMeanr2Cos3PhiCom=in.fMeanr2Cos3PhiCom;
+ fMeanr2Sin3PhiCom=in.fMeanr2Sin3PhiCom;
+ fMeanr2Cos4PhiCom=in.fMeanr2Cos4PhiCom;
+ fMeanr2Sin4PhiCom=in.fMeanr2Sin4PhiCom;
+ fMeanr2Cos5PhiCom=in.fMeanr2Cos5PhiCom;
+ fMeanr2Sin5PhiCom=in.fMeanr2Sin5PhiCom;
+ fMeanr3Cos3PhiCom=in.fMeanr3Cos3PhiCom;
+ fMeanr3Sin3PhiCom=in.fMeanr3Sin3PhiCom;
+ fMeanr4Cos4PhiCom=in.fMeanr4Cos4PhiCom;
+ fMeanr4Sin4PhiCom=in.fMeanr4Sin4PhiCom;
+ fMeanr5Cos5PhiCom=in.fMeanr5Cos5PhiCom;
+ fMeanr5Sin5PhiCom=in.fMeanr5Sin5PhiCom;
+ fSx2Parts=in.fSx2Parts;
+ fSy2Parts=in.fSy2Parts;
+ fSxyParts=in.fSxyParts;
+ fSx2Coll=in.fSx2Coll;
+ fSy2Coll=in.fSy2Coll;
+ fSxyColl=in.fSxyColl;
+ fSx2Com=in.fSx2Com;
+ fSy2Com=in.fSy2Com;
+ fSxyCom=in.fSxyCom;
+ fX=in.fX;
+ fNpp=in.fNpp;
+ return *this;
}
//______________________________________________________________________________
fANucleus.ThrowNucleons(-bgen/2.);
fNucleonsA = fANucleus.GetNucleons();
fAN = fANucleus.GetN();
+ fQAN = fAN * 3;
+ //fAN = 3 * fANucleus.GetN(); // for Pb, Number of quark = 3*208;
for (Int_t i = 0; i<fAN; i++)
{
- AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->At(i));
+ AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->UncheckedAt(i));
nucleonA->SetInNucleusA();
}
fBNucleus.ThrowNucleons(bgen/2.);
fNucleonsB = fBNucleus.GetNucleons();
+ //fBN = 3 * fBNucleus.GetN(); // Number of quark = number of nucleus*3;
fBN = fBNucleus.GetN();
+ fQBN = fBN * 3;
for (Int_t i = 0; i<fBN; i++)
{
- AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->At(i));
+ AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->UncheckedAt(i));
nucleonB->SetInNucleusB();
}
// for each of the A nucleons in nucleus B
for (Int_t i = 0; i<fBN; i++)
{
- AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->At(i));
+ AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->UncheckedAt(i));
for (Int_t j = 0 ; j < fAN ; j++)
{
- AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->At(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;
fNpart=0;
fNcoll=0;
- fMeanX2=0;
- fMeanY2=0;
- fMeanXY=0;
- fMeanXParts=0;
- fMeanYParts=0;
- fMeanXSystem=0;
- fMeanYSystem=0;
- fMeanX_A=0;
- fMeanY_A=0;
- fMeanX_B=0;
- fMeanY_B=0;
+ fNcom=0;
+ fONpart=0;
+ fONcoll=0;
+ fONcom=0;
+ fMeanX2=0.;
+ fMeanY2=0.;
+ fMeanXY=0.;
+ fMeanXParts=0.;
+ fMeanYParts=0.;
+ fMeanOXParts=0.;
+ fMeanOYParts=0.;
+ fMeanXColl=0.;
+ fMeanYColl=0.;
+ fMeanOXColl=0.;
+ fMeanOYColl=0.;
+ fMeanX2Coll=0.;
+ fMeanY2Coll=0.;
+ fMeanXYColl=0.;
+ fMeanXCom=0.;
+ fMeanYCom=0.;
+ fMeanOXCom=0.;
+ fMeanOYCom=0.;
+ fMeanX2Com=0.;
+ fMeanY2Com=0.;
+ fMeanXYCom=0.;
+ fMeanXSystem=0.;
+ fMeanYSystem=0.;
+ fMeanXA=0.;
+ fMeanYA=0.;
+ fMeanXB=0.;
+ fMeanYB=0.;
+ fMeanOXA=0.;
+ fMeanOYA=0.;
+ fMeanOXB=0.;
+ fMeanOYB=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.;
+ fMeanr2Coll=0.;
+ fMeanr3Coll=0.;
+ fMeanr4Coll=0.;
+ fMeanr5Coll=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.;
+ fMeanr2Com=0.;
+ fMeanr3Com=0.;
+ fMeanr4Com=0.;
+ fMeanr5Com=0.;
+ fMeanr2Cos2PhiCom=0.;
+ fMeanr2Sin2PhiCom=0.;
+ fMeanr2Cos3PhiCom=0.;
+ fMeanr2Sin3PhiCom=0.;
+ fMeanr2Cos4PhiCom=0.;
+ fMeanr2Sin4PhiCom=0.;
+ fMeanr2Cos5PhiCom=0.;
+ fMeanr2Sin5PhiCom=0.;
+ fMeanr3Cos3PhiCom=0.;
+ fMeanr3Sin3PhiCom=0.;
+ fMeanr4Cos4PhiCom=0.;
+ fMeanr4Sin4PhiCom=0.;
+ fMeanr5Cos5PhiCom=0.;
+ fMeanr5Sin5PhiCom=0.;
for (Int_t i = 0; i<fAN; i++)
{
- AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->At(i));
- Double_t xA=nucleonA->GetX();
- Double_t yA=nucleonA->GetY();
- fMeanXSystem += xA;
- fMeanYSystem += yA;
- fMeanX_A += xA;
- fMeanY_A += yA;
+ AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->UncheckedAt(i));
+ Double_t oXA = nucleonA->GetX();
+ Double_t oYA = nucleonA->GetY();
+ //fMeanOXSystem += oXA;
+ //fMeanOYSystem += oYA;
+ fMeanOXA += oXA;
+ fMeanOYA += oYA;
if(nucleonA->IsWounded())
{
- fNpart++;
- fMeanXParts += xA;
- fMeanYParts += yA;
- fMeanX2 += xA * xA;
- fMeanY2 += yA * yA;
- fMeanXY += xA * yA;
+ fONpart++;
+ fMeanOXParts += oXA;
+ fMeanOYParts += oYA;
+
+ fONcom += (1-0.150);
+ fMeanOXCom += oXA*(1-0.150);
+ fMeanOYCom += oXA*(1-0.150);
}
}
for (Int_t i = 0; i<fBN; i++)
{
- AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->At(i));
- Double_t xB=nucleonB->GetX();
- Double_t yB=nucleonB->GetY();
- fMeanXSystem += xB;
- fMeanYSystem += yB;
- fMeanX_B += xB;
- fMeanY_B += yB;
-
+ AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->UncheckedAt(i));
+ Double_t oXB=nucleonB->GetX();
+ Double_t oYB=nucleonB->GetY();
+
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();
+ Int_t oNcoll = nucleonB->GetNColl();
+ fONpart++;
+ fMeanOXParts += oXB;
+ fMeanOXColl += oXB*oNcoll;
+ fMeanOXCom += oXB*((1-0.150)+0.150*oNcoll);
+ fMeanOYParts += oYB;
+ fMeanOYColl += oYB*oNcoll;
+ fMeanOYCom += oYB*((1-0.150)+0.150*oNcoll);
+ fONcoll += oNcoll;
+ fONcom += ((1-0.150)+0.150*oNcoll);
}
}
- if (fNpart>0)
+ if (fONpart>0)
{
- fMeanXParts /= fNpart;
- fMeanYParts /= fNpart;
- fMeanX2 /= fNpart;
- fMeanY2 /= fNpart;
- fMeanXY /= fNpart;
+ fMeanOXParts /= fONpart;
+ fMeanOYParts /= fONpart;
}
else
{
- fMeanXParts = 0;
- fMeanYParts = 0;
- fMeanX2 = 0;
- fMeanY2 = 0;
- fMeanXY = 0;
+ fMeanOXParts = 0;
+ fMeanOYParts = 0;
}
- if (fNcoll>0)
+ if (fONcoll>0)
+ {
+ fMeanOXColl /= fONcoll;
+ fMeanOYColl /= fONcoll;
+ }
+ else
+ {
+ fMeanOXColl = 0;
+ fMeanOYColl = 0;
+ }
+
+ if (fONcom>0)
+ {
+ fMeanOXCom /= fONcom;
+ fMeanOYCom /= fONcom;
+ }
+ else
{
- fMeanXColl /= fNcoll;
- fMeanYColl /= fNcoll;
- fMeanX2Coll /= fNcoll;
- fMeanY2Coll /= fNcoll;
- fMeanXYColl /= fNcoll;
+ fMeanOXCom = 0;
+ fMeanOYCom = 0;
+ }
+
+ //////////////////////////////////////////////////////////////////
+ for (Int_t i = 0; i<fAN; i++)
+ {
+ AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->UncheckedAt(i));
+ Double_t xAA = nucleonA->GetX(); // X
+ Double_t yAA = nucleonA->GetY(); // Y
+ Double_t xAPart = xAA - fMeanOXParts; // X'
+ Double_t yAPart = yAA - fMeanOYParts; // Y'
+ Double_t r2APart = xAPart *xAPart+yAPart*yAPart; // r'^2
+ Double_t rAPart = TMath::Sqrt(r2APart); // r'
+ Double_t r3APart = r2APart*rAPart;
+ Double_t r4APart = r3APart*rAPart;
+ Double_t r5APart = r4APart*rAPart;
+ Double_t phiAPart = TMath::ATan2(yAPart,xAPart);
+ Double_t sin2PhiAPart = TMath::Sin(2.*phiAPart);
+ Double_t cos2PhiAPart = TMath::Cos(2.*phiAPart);
+ Double_t sin3PhiAPart = TMath::Sin(3.*phiAPart);
+ Double_t cos3PhiAPart = TMath::Cos(3.*phiAPart);
+ Double_t sin4PhiAPart = TMath::Sin(4.*phiAPart);
+ Double_t cos4PhiAPart = TMath::Cos(4.*phiAPart);
+ Double_t sin5PhiAPart = TMath::Sin(5.*phiAPart);
+ Double_t cos5PhiAPart = TMath::Cos(5.*phiAPart);
+
+ // for combine
+ Double_t xACom = xAA - fMeanOXCom; // X'-Combine
+ Double_t yACom = yAA - fMeanOYCom; // Y'-C
+ Double_t r2ACom = xACom *xACom+yACom*yACom; // r'^2-C
+ Double_t rACom = TMath::Sqrt(r2ACom); // r'-C
+ Double_t r3ACom = r2ACom*rACom;
+ Double_t r4ACom = r3ACom*rACom;
+ Double_t r5ACom = r4ACom*rACom;
+ Double_t phiACom = TMath::ATan2(yACom,xACom);
+ Double_t sin2PhiACom = TMath::Sin(2.*phiACom);
+ Double_t cos2PhiACom = TMath::Cos(2.*phiACom);
+ Double_t sin3PhiACom = TMath::Sin(3.*phiACom);
+ Double_t cos3PhiACom = TMath::Cos(3.*phiACom);
+ Double_t sin4PhiACom = TMath::Sin(4.*phiACom);
+ Double_t cos4PhiACom = TMath::Cos(4.*phiACom);
+ Double_t sin5PhiACom = TMath::Sin(5.*phiACom);
+ Double_t cos5PhiACom = TMath::Cos(5.*phiACom);
+
+ fMeanXSystem += xAA;
+ fMeanYSystem += yAA;
+ fMeanXA += xAA;
+ fMeanYA += yAA;
+ fMeanX2 += xAA * xAA;
+ fMeanY2 += yAA * yAA;
+ fMeanXY += xAA * yAA;
+
+ if(nucleonA->IsWounded())
+ {
+ //Wounded
+ fNpart++;
+ fMeanXParts += xAPart;
+ fMeanYParts += yAPart;
+ fMeanX2Parts += xAPart * xAPart;
+ fMeanY2Parts += yAPart * yAPart;
+ fMeanXYParts += xAPart * yAPart;
+ fMeanr2 += r2APart;
+ fMeanr3 += r3APart;
+ fMeanr4 += r4APart;
+ fMeanr5 += r5APart;
+ fMeanr2Cos2Phi += r2APart*cos2PhiAPart;
+ fMeanr2Sin2Phi += r2APart*sin2PhiAPart;
+ fMeanr2Cos3Phi += r2APart*cos3PhiAPart;
+ fMeanr2Sin3Phi += r2APart*sin3PhiAPart;
+ fMeanr2Cos4Phi += r2APart*cos4PhiAPart;
+ fMeanr2Sin4Phi += r2APart*sin4PhiAPart;
+ fMeanr2Cos5Phi += r2APart*cos5PhiAPart;
+ fMeanr2Sin5Phi += r2APart*sin5PhiAPart;
+ fMeanr3Cos3Phi += r3APart*cos3PhiAPart;
+ fMeanr3Sin3Phi += r3APart*sin3PhiAPart;
+ fMeanr4Cos4Phi += r4APart*cos4PhiAPart;
+ fMeanr4Sin4Phi += r4APart*sin4PhiAPart;
+ fMeanr5Cos5Phi += r5APart*cos5PhiAPart;
+ fMeanr5Sin5Phi += r5APart*sin5PhiAPart;
+
+ // Combined
+ fMeanXCom += xACom*(1-0.150);
+ fMeanYCom += yACom*(1-0.150);
+ fMeanX2Com += xACom*xACom*(1-0.150);
+ fMeanY2Com += yACom*yACom*(1-0.150);
+ fMeanXYCom += xACom*yACom*(1-0.150);
+ fNcom += (1-0.150);
+ fMeanr2Com += r2ACom*(1-0.150);
+ fMeanr3Com += r3ACom*(1-0.150);
+ fMeanr4Com += r4ACom*(1-0.150);
+ fMeanr5Com += r5ACom*(1-0.150);
+ fMeanr2Cos2PhiCom += r2ACom*cos2PhiACom*(1-0.150);
+ fMeanr2Sin2PhiCom += r2ACom*sin2PhiACom*(1-0.150);
+ fMeanr2Cos3PhiCom += r2ACom*cos3PhiACom*(1-0.150);
+ fMeanr2Sin3PhiCom += r2ACom*sin3PhiACom*(1-0.150);
+ fMeanr2Cos4PhiCom += r2ACom*cos4PhiACom*(1-0.150);
+ fMeanr2Sin4PhiCom += r2ACom*sin4PhiACom*(1-0.150);
+ fMeanr2Cos5PhiCom += r2ACom*cos5PhiACom*(1-0.150);
+ fMeanr2Sin5PhiCom += r2ACom*sin5PhiACom*(1-0.150);
+ fMeanr3Cos3PhiCom += r3ACom*cos3PhiACom*(1-0.150);
+ fMeanr3Sin3PhiCom += r3ACom*sin3PhiACom*(1-0.150);
+ fMeanr4Cos4PhiCom += r4ACom*cos4PhiACom*(1-0.150);
+ fMeanr4Sin4PhiCom += r4ACom*sin4PhiACom*(1-0.150);
+ fMeanr5Cos5PhiCom += r5ACom*cos5PhiACom*(1-0.150);
+ fMeanr5Sin5PhiCom += r5ACom*sin5PhiACom*(1-0.150);
+
+ }
+ }
+
+ for (Int_t i = 0; i<fBN; i++)
+ {
+ AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->UncheckedAt(i));
+ Double_t xBB = nucleonB->GetX();
+ Double_t yBB = nucleonB->GetY();
+ // for Wounded
+ Double_t xBPart = xBB - fMeanOXParts; // X'
+ Double_t yBPart = yBB - fMeanOYParts; // Y'
+ Double_t r2BPart = xBPart*xBPart+yBPart*yBPart; // r'^2
+ Double_t rBPart = TMath::Sqrt(r2BPart); // r'
+ Double_t r3BPart = r2BPart*rBPart;
+ Double_t r4BPart = r3BPart*rBPart;
+ Double_t r5BPart = r4BPart*rBPart;
+ Double_t phiBPart = TMath::ATan2(yBPart,xBPart);
+ Double_t sin2PhiBPart = TMath::Sin(2.*phiBPart);
+ Double_t cos2PhiBPart = TMath::Cos(2.*phiBPart);
+ Double_t sin3PhiBPart = TMath::Sin(3.*phiBPart);
+ Double_t cos3PhiBPart = TMath::Cos(3.*phiBPart);
+ Double_t sin4PhiBPart = TMath::Sin(4.*phiBPart);
+ Double_t cos4PhiBPart = TMath::Cos(4.*phiBPart);
+ Double_t sin5PhiBPart = TMath::Sin(5.*phiBPart);
+ Double_t cos5PhiBPart = TMath::Cos(5.*phiBPart);
+ // for Binary
+ Double_t xBColl = xBB - fMeanOXColl; // X'-Binary
+ Double_t yBColl = yBB - fMeanOYColl; // Y'-B
+ Double_t r2BColl = xBColl*xBColl+yBColl*yBColl; // r'^2-B
+ Double_t rBColl = TMath::Sqrt(r2BColl); // r'-B
+ Double_t r3BColl = r2BColl*rBColl;
+ Double_t r4BColl = r3BColl*rBColl;
+ Double_t r5BColl = r4BColl*rBColl;
+ Double_t phiBColl = TMath::ATan2(yBColl,xBColl);
+ Double_t sin2PhiBColl = TMath::Sin(2.*phiBColl);
+ Double_t cos2PhiBColl = TMath::Cos(2.*phiBColl);
+ Double_t sin3PhiBColl = TMath::Sin(3.*phiBColl);
+ Double_t cos3PhiBColl = TMath::Cos(3.*phiBColl);
+ Double_t sin4PhiBColl = TMath::Sin(4.*phiBColl);
+ Double_t cos4PhiBColl = TMath::Cos(4.*phiBColl);
+ Double_t sin5PhiBColl = TMath::Sin(5.*phiBColl);
+ Double_t cos5PhiBColl = TMath::Cos(5.*phiBColl);
+ // for combine
+ Double_t xBCom = xBB - fMeanOXCom; // X'-Combine
+ Double_t yBCom = yBB - fMeanOYCom; // Y'-C
+ Double_t r2BCom = xBCom *xBCom+yBCom*yBCom; // r'^2-C
+ Double_t rBCom = TMath::Sqrt(r2BCom); // r'-C
+ Double_t r3BCom = r2BCom*rBCom;
+ Double_t r4BCom = r3BCom*rBCom;
+ Double_t r5BCom = r4BCom*rBCom;
+ Double_t phiBCom = TMath::ATan2(yBCom,xBCom);
+ Double_t sin2PhiBCom = TMath::Sin(2.*phiBCom);
+ Double_t cos2PhiBCom = TMath::Cos(2.*phiBCom);
+ Double_t sin3PhiBCom = TMath::Sin(3.*phiBCom);
+ Double_t cos3PhiBCom = TMath::Cos(3.*phiBCom);
+ Double_t sin4PhiBCom = TMath::Sin(4.*phiBCom);
+ Double_t cos4PhiBCom = TMath::Cos(4.*phiBCom);
+ Double_t sin5PhiBCom = TMath::Sin(5.*phiBCom);
+ Double_t cos5PhiBCom = TMath::Cos(5.*phiBCom);
+
+ fMeanXSystem += xBB;
+ fMeanYSystem += yBB;
+ fMeanXB += xBB;
+ fMeanYB += yBB;
+ fMeanX2 += xBB*xBB;
+ fMeanY2 += yBB*yBB;
+ fMeanXY += xBB*yBB;
+
+ if(nucleonB->IsWounded())
+ {
+ Int_t ncoll = nucleonB->GetNColl();
+ fNpart++;
+ fMeanXParts += xBPart;
+ fMeanXColl += xBColl*ncoll;
+ fMeanXCom += xBCom*((1-0.150)+0.150*ncoll);
+ fMeanYParts += yBPart;
+ fMeanYColl += yBColl*ncoll;
+ fMeanYCom += yBCom*((1-0.150)+0.150*ncoll);
+ fMeanX2Parts += xBPart * xBPart;
+ fMeanX2Coll += xBColl*xBColl*ncoll;
+ fMeanX2Com += xBCom*xBCom*((1-0.150)+0.150*ncoll);
+ fMeanY2Parts += yBPart * yBPart;
+ fMeanY2Coll += yBColl*yBColl*ncoll;
+ fMeanY2Com += yBCom*yBCom*((1-0.150)+0.150*ncoll);
+ fMeanXYParts += xBPart * yBPart;
+ fMeanXYColl += xBColl*yBColl*ncoll;
+ fMeanXYCom += xBCom*yBCom*((1-0.150)+0.150*ncoll);
+ fNcoll += ncoll;
+ fNcom += ((1-0.150)+0.150*ncoll);
+ fMeanr2 += r2BPart;
+ fMeanr3 += r3BPart;
+ fMeanr4 += r4BPart;
+ fMeanr5 += r5BPart;
+ fMeanr2Cos2Phi += r2BPart*cos2PhiBPart;
+ fMeanr2Sin2Phi += r2BPart*sin2PhiBPart;
+ fMeanr2Cos3Phi += r2BPart*cos3PhiBPart;
+ fMeanr2Sin3Phi += r2BPart*sin3PhiBPart;
+ fMeanr2Cos4Phi += r2BPart*cos4PhiBPart;
+ fMeanr2Sin4Phi += r2BPart*sin4PhiBPart;
+ fMeanr2Cos5Phi += r2BPart*cos5PhiBPart;
+ fMeanr2Sin5Phi += r2BPart*sin5PhiBPart;
+ fMeanr3Cos3Phi += r3BPart*cos3PhiBPart;
+ fMeanr3Sin3Phi += r3BPart*sin3PhiBPart;
+ fMeanr4Cos4Phi += r4BPart*cos4PhiBPart;
+ fMeanr4Sin4Phi += r4BPart*sin4PhiBPart;
+ fMeanr5Cos5Phi += r5BPart*cos5PhiBPart;
+ fMeanr5Sin5Phi += r5BPart*sin5PhiBPart;
+ fMeanr2Coll += r2BColl*ncoll;
+ fMeanr3Coll += r3BColl*ncoll;
+ fMeanr4Coll += r4BColl*ncoll;
+ fMeanr5Coll += r5BColl*ncoll;
+ fMeanr2Cos2PhiColl += r2BColl*cos2PhiBColl*ncoll;
+ fMeanr2Sin2PhiColl += r2BColl*sin2PhiBColl*ncoll;
+ fMeanr2Cos3PhiColl += r2BColl*cos3PhiBColl*ncoll;
+ fMeanr2Sin3PhiColl += r2BColl*sin3PhiBColl*ncoll;
+ fMeanr2Cos4PhiColl += r2BColl*cos4PhiBColl*ncoll;
+ fMeanr2Sin4PhiColl += r2BColl*sin4PhiBColl*ncoll;
+ fMeanr2Cos5PhiColl += r2BColl*cos5PhiBColl*ncoll;
+ fMeanr2Sin5PhiColl += r2BColl*sin5PhiBColl*ncoll;
+ fMeanr3Cos3PhiColl += r3BColl*cos3PhiBColl*ncoll;
+ fMeanr3Sin3PhiColl += r3BColl*sin3PhiBColl*ncoll;
+ fMeanr4Cos4PhiColl += r4BColl*cos4PhiBColl*ncoll;
+ fMeanr4Sin4PhiColl += r4BColl*sin4PhiBColl*ncoll;
+ fMeanr5Cos5PhiColl += r5BColl*cos5PhiBColl*ncoll;
+ fMeanr5Sin5PhiColl += r5BColl*sin5PhiBColl*ncoll;
+ fMeanr2Com += r2BCom*((1-0.150)+0.150*ncoll);
+ fMeanr3Com += r3BCom*((1-0.150)+0.150*ncoll);
+ fMeanr4Com += r4BCom*((1-0.150)+0.150*ncoll);
+ fMeanr5Com += r5BCom*((1-0.150)+0.150*ncoll);
+ fMeanr2Cos2PhiCom += r2BCom*cos2PhiBCom*((1-0.150)+0.150*ncoll);
+ fMeanr2Sin2PhiCom += r2BCom*sin2PhiBCom*((1-0.150)+0.150*ncoll);
+ fMeanr2Cos3PhiCom += r2BCom*cos3PhiBCom*((1-0.150)+0.150*ncoll);
+ fMeanr2Sin3PhiCom += r2BCom*sin3PhiBCom*((1-0.150)+0.150*ncoll);
+ fMeanr2Cos4PhiCom += r2BCom*cos4PhiBCom*((1-0.150)+0.150*ncoll);
+ fMeanr2Sin4PhiCom += r2BCom*sin4PhiBCom*((1-0.150)+0.150*ncoll);
+ fMeanr2Cos5PhiCom += r2BCom*cos5PhiBCom*((1-0.150)+0.150*ncoll);
+ fMeanr2Sin5PhiCom += r2BCom*sin5PhiBCom*((1-0.150)+0.150*ncoll);
+ fMeanr3Cos3PhiCom += r3BCom*cos3PhiBCom*((1-0.150)+0.150*ncoll);
+ fMeanr3Sin3PhiCom += r3BCom*sin3PhiBCom*((1-0.150)+0.150*ncoll);
+ fMeanr4Cos4PhiCom += r4BCom*cos4PhiBCom*((1-0.150)+0.150*ncoll);
+ fMeanr4Sin4PhiCom += r4BCom*sin4PhiBCom*((1-0.150)+0.150*ncoll);
+ fMeanr5Cos5PhiCom += r5BCom*cos5PhiBCom*((1-0.150)+0.150*ncoll);
+ fMeanr5Sin5PhiCom += r5BCom*sin5PhiBCom*((1-0.150)+0.150*ncoll);
+ }
+ }
+
+ if (fNpart>0)
+ {
+ fMeanXParts /= fNpart;
+ fMeanYParts /= fNpart;
+ fMeanX2Parts /= fNpart;
+ fMeanY2Parts /= fNpart;
+ fMeanXYParts /= 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;
+ fMeanX2Parts = 0;
+ fMeanY2Parts = 0;
+ fMeanXYParts = 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;
}
+
else
{
fMeanXColl = 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 (fNcom>0)
+ {
+ fMeanXCom /= fNcom;
+ fMeanYCom /= fNcom;
+ fMeanX2Com /= fNcom;
+ fMeanY2Com /= fNcom;
+ fMeanXYCom /= fNcom;
+ fMeanr2Com /= fNcom;
+ fMeanr3Com /= fNcom;
+ fMeanr4Com /= fNcom;
+ fMeanr5Com /= fNcom;
+ fMeanr2Cos2PhiCom /= fNcom;
+ fMeanr2Sin2PhiCom /= fNcom;
+ fMeanr2Cos3PhiCom /= fNcom;
+ fMeanr2Sin3PhiCom /= fNcom;
+ fMeanr2Cos4PhiCom /= fNcom;
+ fMeanr2Sin4PhiCom /= fNcom;
+ fMeanr2Cos5PhiCom /= fNcom;
+ fMeanr2Sin5PhiCom /= fNcom;
+ fMeanr3Cos3PhiCom /= fNcom;
+ fMeanr3Sin3PhiCom /= fNcom;
+ fMeanr4Cos4PhiCom /= fNcom;
+ fMeanr4Sin4PhiCom /= fNcom;
+ fMeanr5Cos5PhiCom /= fNcom;
+ fMeanr5Sin5PhiCom /= fNcom;
+ }
+
+ else
+ {
+ fMeanXCom = 0;
+ fMeanYCom = 0;
+ fMeanX2Com = 0;
+ fMeanY2Com = 0;
+ fMeanXYCom = 0;
+ fMeanr2Cos2PhiCom =0;
+ fMeanr2Sin2PhiCom =0;
+ fMeanr2Cos3PhiCom =0;
+ fMeanr2Sin3PhiCom =0;
+ fMeanr2Cos4PhiCom =0;
+ fMeanr2Sin4PhiCom =0;
+ fMeanr2Cos5PhiCom =0;
+ fMeanr2Sin5PhiCom =0;
+ fMeanr3Cos3PhiCom =0;
+ fMeanr3Sin3PhiCom =0;
+ fMeanr4Cos4PhiCom =0;
+ fMeanr4Sin4PhiCom =0;
+ fMeanr5Cos5PhiCom =0;
+ fMeanr5Sin5PhiCom =0;
}
if(fAN+fBN>0)
{
fMeanXSystem /= (fAN + fBN);
fMeanYSystem /= (fAN + fBN);
+ fMeanX2 /= (fAN + fBN);
+ fMeanY2 /= (fAN + fBN);
+ fMeanXY /= (fAN + fBN);
}
else
{
fMeanXSystem = 0;
fMeanYSystem = 0;
+ fMeanX2 = 0;
+ fMeanY2 = 0;
+ fMeanXY = 0;
}
if(fAN>0)
{
- fMeanX_A /= fAN;
- fMeanY_A /= fAN;
+ fMeanXA /= fAN;
+ fMeanYA /= fAN;
}
else
{
- fMeanX_A = 0;
- fMeanY_A = 0;
+ fMeanXA = 0;
+ fMeanYA = 0;
}
if(fBN>0)
{
- fMeanX_B /= fBN;
- fMeanY_B /= fBN;
+ fMeanXB /= fBN;
+ fMeanYB /= fBN;
}
else
{
- fMeanX_B = 0;
- fMeanY_B = 0;
+ fMeanXB = 0;
+ fMeanYB = 0;
}
+
- fSx2=fMeanX2-(fMeanXParts*fMeanXParts);
- fSy2=fMeanY2-(fMeanYParts*fMeanYParts);
- fSxy=fMeanXY-fMeanXParts*fMeanYParts;
+
+ //////////////////////////////////////////////////////////////////
+ fSx2Parts=fMeanX2Parts-(fMeanXParts*fMeanXParts);
+ fSy2Parts=fMeanY2Parts-(fMeanYParts*fMeanYParts);
+ fSxyParts=fMeanXYParts-fMeanXParts*fMeanYParts;
fSx2Coll=fMeanX2Coll-(fMeanXColl*fMeanXColl);
fSy2Coll=fMeanY2Coll-(fMeanYColl*fMeanYColl);
fSxyColl=fMeanXYColl-fMeanXColl*fMeanYColl;
-
- fB_MC = bgen;
-
+ fSx2Com=fMeanX2Com-(fMeanXCom*fMeanXCom);
+ fSy2Com=fMeanY2Com-(fMeanYCom*fMeanYCom);
+ fSxyCom=fMeanXYCom-fMeanXCom*fMeanYCom;
+ fBMC = bgen;
fTotalEvents++;
if (fNpart>0) fEvents++;
-
if (fNpart==0) return kFALSE;
if (fNpart > fMaxNpartFound) fMaxNpartFound = fNpart;
//______________________________________________________________________________
void AliGlauberMC::Draw(Option_t* /*option*/)
{
+ //draw method
fANucleus.Draw(fXSect, 2);
fBNucleus.Draw(fXSect, 4);
//______________________________________________________________________________
Double_t AliGlauberMC::GetTotXSect() const
{
+ //total xsection
return (1.*fEvents/fTotalEvents)*TMath::Pi()*fBMax*fBMax/100;
}
//______________________________________________________________________________
Double_t AliGlauberMC::GetTotXSectErr() const
{
+ //total xsection error
return GetTotXSect()/TMath::Sqrt((Double_t)fEvents) *
TMath::Sqrt(Double_t(1.-fEvents/fTotalEvents));
}
//______________________________________________________________________________
TObjArray *AliGlauberMC::GetNucleons()
{
+ //get array of nucleons
if(!fNucleonsA || !fNucleonsB) return 0;
fNucleonsA->SetOwner(0);
fNucleonsB->SetOwner(0);
allnucleons->SetOwner();
for (Int_t i = 0; i<fAN; i++)
{
- allnucleons->Add(fNucleonsA->At(i));
+ allnucleons->Add(fNucleonsA->UncheckedAt(i));
}
for (Int_t i = 0; i<fBN; i++)
{
- allnucleons->Add(fNucleonsB->At(i));
+ allnucleons->Add(fNucleonsB->UncheckedAt(i));
}
return allnucleons;
}
+
//______________________________________________________________________________
Double_t AliGlauberMC::NegativeBinomialDistribution(Int_t x, Int_t k, Double_t nmean)
{
+ //produce a number distributed acc. neg.bin.dist
if(k<=0)
{
cout << "Error, zero or negative K" << endl;
*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))));
}
+
//______________________________________________________________________________
-Int_t AliGlauberMC::NegativeBinomialRandom(Int_t k, Double_t nmean)
-//return random integer from a Negative Binomial Distribution
+Int_t AliGlauberMC::NegativeBinomialRandom(Int_t k, Double_t nmean) const
{
+ //return random integer from a Negative Binomial Distribution
static const Int_t fMaxPlot = 1000;
Double_t array[fMaxPlot];
array[0] = NegativeBinomialDistribution(0,k,nmean);
}
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 )
+ 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;
return TMath::BinarySearch(fMaxPlot,array,r)+2;
}
-
//______________________________________________________________________________
-void AliGlauberMC::SetdNdEtaParam(Double_t nnp, Double_t x)
+Double_t AliGlauberMC::GetdNdEta() const
{
- // set parameters used for calculating multiplicity, see GetdNdEta() for commments
- fdNdEtaParam[0]=nnp;
- fdNdEtaParam[1]=x;
+ 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;
+ }
}
//______________________________________________________________________________
-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)
-{
- // set parameters used for calculating multiplicity see GetdNdEtaNBD() for commments
- fdNdEtaNBDParam[0]=k;
- fdNdEtaNBDParam[1]=nmean;
- fdNdEtaNBDParam[2]=beta;
-}
-//______________________________________________________________________________
-void AliGlauberMC::SetdNdEtaTwoNBDParam(Double_t alpha, Double_t k1, Double_t nmean1, Double_t k2, Double_t nmean2, Double_t beta)
-{
- // set parameters used for calculating multiplicity see GetdNdEtaTwoNBD() for commments
- fdNdEtaTwoNBDParam[0]=alpha;
- fdNdEtaTwoNBDParam[1]=k1;
- fdNdEtaTwoNBDParam[2]=nmean1;
- fdNdEtaTwoNBDParam[3]=k2;
- fdNdEtaTwoNBDParam[4]=nmean2;
- fdNdEtaTwoNBDParam[5]=beta;
-}
-//______________________________________________________________________________
-Double_t AliGlauberMC::GetdNdEta(Double_t nnp, Double_t x)
+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)
+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);
}
+
//_______________________________________________________________________________
+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_t AliGlauberMC::GetdNdEtaNBD ( Int_t k, Double_t nmean, Double_t beta)
+ 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++)
{
}
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 )
+//______________________________________________________________________________
+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++)
{
}
Double_t mul=(1-beta)*mulnp/2+beta*mulnb;
return mul;
-
}
+
//______________________________________________________________________________
Double_t AliGlauberMC::GetEccentricityPart() const
{
//get participant eccentricity of participants
- return (TMath::Sqrt((fSy2-fSx2)*(fSy2-fSx2)+4*fSxy*fSxy)/(fSy2+fSx2));
+ if (fNpart<2) return 0.0;
+ return (TMath::Sqrt((fSy2Parts-fSx2Parts)*(fSy2Parts-fSx2Parts)+4*fSxyParts*fSxyParts)/(fSy2Parts+fSx2Parts));
}
-//______________________________________________________________________________
+//_____________________________________________________________________________
Double_t AliGlauberMC::GetEccentricityPartColl() const
{
//get participant eccentricity of binary collisions
+ if (fNcoll<0) return 0.0;
return (TMath::Sqrt((fSy2Coll-fSx2Coll)*(fSy2Coll-fSx2Coll)+4*fSxyColl*fSxyColl)/(fSy2Coll+fSx2Coll));
}
+//_____________________________________________________________________________
+Double_t AliGlauberMC::GetEccentricityPartCom() const
+{
+ //get participant eccentricity of binary collisions
+ if (fNcom<0) return 0.0;
+ return (TMath::Sqrt((fSy2Com-fSx2Com)*(fSy2Com-fSx2Com)+4*fSxyCom*fSxyCom)/(fSy2Com+fSx2Com));
+}
//______________________________________________________________________________
Double_t AliGlauberMC::GetEccentricity() const
{
//get standard eccentricity of participants
- return ((fSy2-fSx2)/(fSy2+fSx2));
+ if (fNpart<2) return 0.0;
+ return ((fSy2Parts-fSx2Parts)/(fSy2Parts+fSx2Parts));
}
//______________________________________________________________________________
Double_t AliGlauberMC::GetEccentricityColl() const
{
//get standard eccentricity of binary collisions
+ if (fNcoll<0) return 0.0;
return ((fSy2Coll-fSx2Coll)/(fSy2Coll+fSx2Coll));
}
+
+//______________________________________________________________________________
+Double_t AliGlauberMC::GetEccentricityCom() const
+{
+ //get standard eccentricity of combined weight
+ if (fNcom<0) return 0.0;
+ return ((fSy2Com-fSx2Com)/(fSy2Com+fSx2Com));
+ //return (((fMeanY2Com-(fMeanYCom*fMeanYCom)) - (fMeanX2Com-(fMeanXCom*fMeanXCom))) / ((fMeanY2Com-(fMeanYCom*fMeanYCom)) + (fMeanX2Com-(fMeanXCom*fMeanXCom))) );
+}
+
+//______________________________________________________________________________
+Double_t AliGlauberMC::GetStoa() const
+{
+ //get standard Transverse Overlap Area
+ if (fNpart<2) return 0.0;
+ return ( TMath::Pi()*(TMath::Sqrt(fSx2Parts))*(TMath::Sqrt(fSy2Parts)));
+}
+
//______________________________________________________________________________
Bool_t AliGlauberMC::NextEvent(Double_t bgen)
{
{
bgen = TMath::Sqrt((fBMax*fBMax-fBMin*fBMin)*gRandom->Rndm()+fBMin*fBMin);
}
- if ( succes=CalcEvent(bgen) ) break; //ends if we have particparts
+ if ( (succes=CalcEvent(bgen)) ) break; //ends if we have particparts
}
return succes;
}
+//______________________________________________________________________________
+Double_t AliGlauberMC::GetEpsilon2Part() const
+{
+ //get participant eccentricity of participants
+ if (fNpart<2) return 0.0;
+ return (TMath::Sqrt(fMeanr2Cos2Phi*fMeanr2Cos2Phi+fMeanr2Sin2Phi*fMeanr2Sin2Phi)/fMeanr2);
+}
+//______________________________________________________________________________
+Double_t AliGlauberMC::GetEpsilon3Part() const
+{
+ //get participant eccentricity of participants
+ if (fNpart<2) return 0.0;
+ return (TMath::Sqrt(fMeanr2Cos3Phi*fMeanr2Cos3Phi+fMeanr2Sin3Phi*fMeanr2Sin3Phi)/fMeanr2);
+ //return (TMath::Sqrt(fMeanr3Cos3Phi*fMeanr3Cos3Phi+fMeanr3Sin3Phi*fMeanr3Sin3Phi)/fMeanr3);
+}
+//______________________________________________________________________________
+Double_t AliGlauberMC::GetEpsilon4Part() const
+{
+ //get participant eccentricity of participants
+ if (fNpart<2) return 0.0;
+ return (TMath::Sqrt(fMeanr2Cos4Phi*fMeanr2Cos4Phi+fMeanr2Sin4Phi*fMeanr2Sin4Phi)/fMeanr2);
+ //return (TMath::Sqrt(fMeanr4Cos4Phi*fMeanr4Cos4Phi+fMeanr4Sin4Phi*fMeanr4Sin4Phi)/fMeanr4);
+}
+//______________________________________________________________________________
+Double_t AliGlauberMC::GetEpsilon5Part() const
+{
+ //get participant eccentricity of participants
+ if (fNpart<2) return 0.0;
+ return (TMath::Sqrt(fMeanr2Cos5Phi*fMeanr2Cos5Phi+fMeanr2Sin5Phi*fMeanr2Sin5Phi)/fMeanr2);
+ //return (TMath::Sqrt(fMeanr5Cos5Phi*fMeanr5Cos5Phi+fMeanr5Sin5Phi*fMeanr5Sin5Phi)/fMeanr5);
+}
+
+//______________________________________________________________________________
+Double_t AliGlauberMC::GetEpsilon2Coll() const
+{
+ //get epsilon2 of binary collisions
+ if (fNcoll<0) return 0.0;
+ return (TMath::Sqrt(fMeanr2Cos2PhiColl*fMeanr2Cos2PhiColl+fMeanr2Sin2PhiColl*fMeanr2Sin2PhiColl)/fMeanr2Coll);
+}
+
+//______________________________________________________________________________
+Double_t AliGlauberMC::GetEpsilon3Coll() const
+{
+ //get epsilon3 of binary collisions
+ if (fNcoll<0) return 0.0;
+ return (TMath::Sqrt(fMeanr2Cos3PhiColl*fMeanr2Cos3PhiColl+fMeanr2Sin3PhiColl*fMeanr2Sin3PhiColl)/fMeanr2Coll);
+ //return (TMath::Sqrt(fMeanr3Cos3PhiColl*fMeanr3Cos3PhiColl+fMeanr3Sin3PhiColl*fMeanr3Sin3PhiColl)/fMeanr3Coll);
+}
+
+//______________________________________________________________________________
+Double_t AliGlauberMC::GetEpsilon4Coll() const
+{
+ //get epsilon4 of binary collisions
+ if (fNcoll<0) return 0.0;
+ return (TMath::Sqrt(fMeanr2Cos4PhiColl*fMeanr2Cos4PhiColl+fMeanr2Sin4PhiColl*fMeanr2Sin4PhiColl)/fMeanr2Coll);
+ //return (TMath::Sqrt(fMeanr4Cos4PhiColl*fMeanr4Cos4PhiColl+fMeanr4Sin4PhiColl*fMeanr4Sin4PhiColl)/fMeanr4Coll);
+}
+
+//______________________________________________________________________________
+Double_t AliGlauberMC::GetEpsilon5Coll() const
+{
+ //get epsilon5 of binary collisions
+ if (fNcoll<0) return 0.0;
+ return (TMath::Sqrt(fMeanr2Cos5PhiColl*fMeanr2Cos5PhiColl+fMeanr2Sin5PhiColl*fMeanr2Sin5PhiColl)/fMeanr2Coll);
+ //return (TMath::Sqrt(fMeanr5Cos5PhiColl*fMeanr5Cos5PhiColl+fMeanr5Sin5PhiColl*fMeanr5Sin5PhiColl)/fMeanr5Coll);
+}
+
+//______________________________________________________________________________
+Double_t AliGlauberMC::GetEpsilon2Com() const
+{
+ //get epsilon2 of binary collisions
+ if (fNcom<0) return 0.0;
+ return (TMath::Sqrt(fMeanr2Cos2PhiCom*fMeanr2Cos2PhiCom+fMeanr2Sin2PhiCom*fMeanr2Sin2PhiCom)/fMeanr2Com);
+}
+
+//______________________________________________________________________________
+Double_t AliGlauberMC::GetEpsilon3Com() const
+{
+ //get epsilon3 of binary collisions
+ if (fNcom<0) return 0.0;
+ return (TMath::Sqrt(fMeanr2Cos3PhiCom*fMeanr2Cos3PhiCom+fMeanr2Sin3PhiCom*fMeanr2Sin3PhiCom)/fMeanr2Com);
+ //return (TMath::Sqrt(fMeanr3Cos3PhiCom*fMeanr3Cos3PhiCom+fMeanr3Sin3PhiCom*fMeanr3Sin3PhiCom)/fMeanr3Com);
+}
+
+//______________________________________________________________________________
+Double_t AliGlauberMC::GetEpsilon4Com() const
+{
+ //get epsilon4 of binary collisions
+ if (fNcom<0) return 0.0;
+ return (TMath::Sqrt(fMeanr2Cos4PhiCom*fMeanr2Cos4PhiCom+fMeanr2Sin4PhiCom*fMeanr2Sin4PhiCom)/fMeanr2Com);
+ //return (TMath::Sqrt(fMeanr4Cos4PhiCom*fMeanr4Cos4PhiCom+fMeanr4Sin4PhiCom*fMeanr4Sin4PhiCom)/fMeanr4Com);
+}
+
+//______________________________________________________________________________
+Double_t AliGlauberMC::GetEpsilon5Com() const
+{
+ //get epsilon5 of binary collisions
+ if (fNcom<0) return 0.0;
+ return (TMath::Sqrt(fMeanr2Cos5PhiCom*fMeanr2Cos5PhiCom+fMeanr2Sin5PhiCom*fMeanr2Sin5PhiCom)/fMeanr2Com);
+ //return (TMath::Sqrt(fMeanr5Cos5PhiCom*fMeanr5Cos5PhiCom+fMeanr5Sin5PhiCom*fMeanr5Sin5PhiCom)/fMeanr5Com);
+}
+
+//______________________________________________________________________________
+Double_t AliGlauberMC::GetPsi2() const
+{
+ return ((TMath::ATan2(fMeanr2Sin2Phi,fMeanr2Cos2Phi)+TMath::Pi())/2);
+}
+
+//______________________________________________________________________________
+Double_t AliGlauberMC::GetPsi3() const
+{
+ return ((TMath::ATan2(fMeanr2Sin3Phi,fMeanr2Cos3Phi)+TMath::Pi())/3);
+ //return ((TMath::ATan2(fMeanr3Sin3Phi,fMeanr3Cos3Phi)+TMath::Pi())/3);
+}
+
+//______________________________________________________________________________
+Double_t AliGlauberMC::GetPsi4() const
+{
+ return ((TMath::ATan2(fMeanr2Sin4Phi,fMeanr2Cos4Phi)+TMath::Pi())/4);
+ //return ((TMath::ATan2(fMeanr4Sin4Phi,fMeanr4Cos4Phi)+TMath::Pi())/4);
+}
+
+//______________________________________________________________________________
+Double_t AliGlauberMC::GetPsi5() const
+{
+ return ((TMath::ATan2(fMeanr2Sin5Phi,fMeanr2Cos5Phi)+TMath::Pi())/5);
+ //return ((TMath::ATan2(fMeanr5Sin5Phi,fMeanr5Cos5Phi)+TMath::Pi())/5);
+}
+
+/*_____________________________________________________________________________
+Double_t AliGlauberMC::GetE43Part() const
+{
+ //get participant eccentricity of participants
+ if (fNpart<2) return 0.0;
+ //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)));
+ 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)));
+}
+
+
+//______________________________________________________________________________
+Double_t AliGlauberMC::GetE43Coll() const
+{
+ //get epsilon5 of binary collisions
+ if (fNcoll<2) return 0.0;
+ //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)));
+ 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)));
+}
+
+//______________________________________________________________________________
+Double_t AliGlauberMC::GetE43Com() const
+{
+ //get epsilon5 of binary collisions
+ if (fNcom<2) return 0.0;
+ //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)));
+ 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)));
+}
+//___________________________________________________________________________
+
+Double_t AliGlauberMC::GetPsi4m2() const
+{
+ //return (TMath::Cos(4*((TMath::ATan2(fMeanr2Sin4Phi,fMeanr2Cos4Phi)+TMath::Pi())/4) - ((TMath::ATan2(fMeanr2Sin2Phi,fMeanr2Cos2Phi)+TMath::Pi())/2)));
+ return (TMath::Cos(4*(((TMath::ATan2(fMeanr4Sin4Phi,fMeanr4Cos4Phi)+TMath::Pi())/4)-((TMath::ATan2(fMeanr2Sin2Phi,fMeanr2Cos2Phi)+TMath::Pi())/2))));
+}
+*/
//______________________________________________________________________________
void AliGlauberMC::Run(Int_t nevents)
{
+ //example run
cout << "Generating " << nevents << " events..." << endl;
TString name(Form("nt_%s_%s",fANucleus.GetName(),fBNucleus.GetName()));
TString title(Form("%s + %s (x-sect = %d mb)",fANucleus.GetName(),fBNucleus.GetName(),(Int_t) fXSect));
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");
+ "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");
fnt->SetDirectory(0);
}
Int_t q = 0;
u++;
continue;
}
+
q++;
- Float_t v[25];
+ //Float_t v[27];
+ Float_t v[45];
v[0] = GetNpart();
v[1] = GetNcoll();
- v[2] = fB_MC;
+ 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[5] = fMeanX2Parts;
+ v[6] = fMeanY2Parts;
+ v[7] = fMeanXYParts;
+ v[8] = fSx2Parts;
+ v[9] = fSy2Parts;
+ v[10] = fSxyParts;
v[11] = fMeanXSystem;
v[12] = fMeanYSystem;
- v[13] = fMeanX_A;
- v[14] = fMeanY_A;
- v[15] = fMeanX_B;
- v[16] = fMeanY_B;
+ v[13] = fMeanXA;
+ v[14] = fMeanYA;
+ v[15] = fMeanXB;
+ v[16] = fMeanYB;
v[17] = GetEccentricity();
- v[18] = GetEccentricityColl();
- v[19] = GetEccentricityPart();
- v[20] = GetEccentricityPartColl();
- 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] );
+ v[18] = GetStoa();
+ v[19] = GetEccentricityColl();
+ v[20] = GetEccentricityCom();
+ v[21] = GetEccentricityPart();
+ v[22] = GetEccentricityPartColl();
+ v[23] = GetEccentricityPartCom();
+ if (fDoPartProd)
+ {
+ v[24] = GetdNdEta();
+ v[25] = GetdNdEta();
+ v[26] = v[24]+v[25];
+ }
+ else
+ {
+ v[24] = 0;
+ v[25] = 0;
+ v[26] = 0;
+ }
+ v[27]=fXSect;
+
+ Float_t mytAA=-999;
+ if (GetNcoll()>0) mytAA=GetNcoll()/fXSect;
+ v[28]=mytAA;
+ //_____________epsilon2,3,4,4_______
+ v[29] = GetEpsilon2Part();
+ v[30] = GetEpsilon3Part();
+ v[31] = GetEpsilon4Part();
+ v[32] = GetEpsilon5Part();
+ v[33] = GetEpsilon2Coll();
+ v[34] = GetEpsilon3Coll();
+ v[35] = GetEpsilon4Coll();
+ v[36] = GetEpsilon5Coll();
+ v[37] = GetEpsilon2Com();
+ v[38] = GetEpsilon3Com();
+ v[39] = GetEpsilon4Com();
+ v[40] = GetEpsilon5Com();
+ v[41] = GetPsi2();
+ v[42] = GetPsi3();
+ v[43] = GetPsi4();
+ v[44] = GetPsi5();
+ //always at the end
fnt->Fill(v);
- if ((i%50)==0) std::cout << "Generating Event # " << i << "... \r" << flush;
+ 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;
}
//---------------------------------------------------------------------------------
-void AliGlauberMC::runAndSaveNtuple( Int_t n,
- Option_t *sysA,
- Option_t *sysB,
+void AliGlauberMC::RunAndSaveNtuple( Int_t n,
+ const Option_t *sysA,
+ const Option_t *sysB,
Double_t signn,
Double_t mind,
+ Double_t r,
+ Double_t a,
const char *fname)
{
- AliGlauberMC *mcg=new AliGlauberMC(sysA,sysB,signn);
- mcg->SetMinDistance(mind);
- mcg->Run(n);
- TNtuple *nt=mcg->GetNtuple();
+ //example run
+ AliGlauberMC mcg(sysA,sysB,signn);
+ mcg.SetMinDistance(mind);
+ mcg.Setr(r);
+ mcg.Seta(a);
+ mcg.Run(n);
+ TNtuple *nt=mcg.GetNtuple();
TFile out(fname,"recreate",fname,9);
if(nt) nt->Write();
+ printf("total cross section with a nucleon-nucleon cross section \t%f is \t%f",signn,mcg.GetTotXSect());
out.Close();
}
//---------------------------------------------------------------------------------
-void AliGlauberMC::runAndSaveNucleons( Int_t n,
- Option_t *sysA,
+void AliGlauberMC::RunAndSaveNucleons( Int_t n,
+ const Option_t *sysA,
Option_t *sysB,
- Double_t signn,
+ const Double_t signn,
Double_t mind,
Bool_t verbose,
const char *fname)
{
- AliGlauberMC *mcg=new AliGlauberMC(sysA,sysB,signn);
- mcg->SetMinDistance(mind);
+ //example run
+ AliGlauberMC mcg(sysA,sysB,signn);
+ mcg.SetMinDistance(mind);
TFile *out=0;
if(fname)
out=new TFile(fname,"recreate",fname,9);
{
//get an event with at least one collision
- while(!mcg->NextEvent()) {}
+ while(!mcg.NextEvent()) {}
//access, save and (if wanted) print out nucleons
- TObjArray* nucleons=mcg->GetNucleons();
+ 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;
+ 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->At(iNucl);
+ AliGlauberNucleon *nucl=(AliGlauberNucleon *)nucleons->UncheckedAt(iNucl);
Char_t nucleus='A';
if(nucl->IsInNucleusB()) nucleus='B';
Double_t x=nucl->GetX();
if(out) delete out;
}
+//---------------------------------------------------------------------------------
+void AliGlauberMC::Reset()
+{
+ //delete the ntuple
+ delete fnt;
+ fnt=NULL;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff