Including basic QA histograms + steering macro to run the code
authorpchrist <pchrist@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 26 Dec 2011 00:38:57 +0000 (00:38 +0000)
committerpchrist <pchrist@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 26 Dec 2011 00:38:57 +0000 (00:38 +0000)
PWG2/EBYE/AliAnalysisTaskToyModel.cxx
PWG2/EBYE/AliAnalysisTaskToyModel.h
PWG2/EBYE/macros/runBalanceFunctionToyModel.C [new file with mode: 0755]

index 70cf480..b566b15 100755 (executable)
@@ -9,8 +9,8 @@
 #include "TArrayF.h"\r
 #include "TF1.h"\r
 #include "TRandom.h"\r
+#include "TFile.h"\r
 \r
-#include "AliAnalysisTaskSE.h"\r
 #include "AliAnalysisManager.h"\r
 #include "AliLog.h"\r
 \r
 ClassImp(AliAnalysisTaskToyModel)\r
 \r
 //________________________________________________________________________\r
-AliAnalysisTaskToyModel::AliAnalysisTaskToyModel(const char *name) \r
-: AliAnalysisTaskSE(name), \r
-  fBalance(0),\r
-  fRunShuffling(kFALSE), fShuffledBalance(0),\r
-  fList(0), fListBF(0), fListBFS(0),\r
-  fHistEventStats(0),\r
-  fTotalMultiplicityMean(100.), fTotalMultiplicitySigma(0.0),\r
-  fNetChargeMean(0.0), fNetChargeSigma(0.0),\r
-  fPtMin(0.0), fPtMax(100.0),\r
-  fEtaMin(-1.0), fEtaMax(1.0),\r
-  fUseAcceptanceParameterization(kFALSE), fAcceptanceParameterization(0),\r
-  fUseAllCharges(kTRUE), fParticleMass(0.0),\r
-  fPtSpectraAllCharges(0), fTemperatureAllCharges(100.),\r
-  fReactionPlane(0.0),\r
-  fAzimuthalAngleAllCharges(0), fDirectedFlowAllCharges(0.0), \r
-  fEllipticFlowAllCharges(0.0), fTriangularFlowAllCharges(0.0),\r
-  fQuandrangularFlowAllCharges(0.0), fPentangularFlowAllCharges(0.0),\r
-  fPionPercentage(0.8), fPionMass(0.0),\r
-  fPtSpectraPions(0), fTemperaturePions(100.),\r
-  fAzimuthalAnglePions(0), fDirectedFlowPions(0.0), \r
-  fEllipticFlowPions(0.0), fTriangularFlowPions(0.0), \r
-  fQuandrangularFlowPions(0.0), fPentangularFlowPions(0.0),\r
-  fKaonPercentage(0.8), fKaonMass(0.0),\r
-  fPtSpectraKaons(0), fTemperatureKaons(100.),\r
-  fAzimuthalAngleKaons(0), fDirectedFlowKaons(0.0), \r
-  fEllipticFlowKaons(0.0), fTriangularFlowKaons(0.0),\r
-  fQuandrangularFlowKaons(0.0), fPentangularFlowKaons(0.0),\r
-  fProtonPercentage(0.8), fProtonMass(0.0),\r
-  fPtSpectraProtons(0), fTemperatureProtons(100.),\r
-  fAzimuthalAngleProtons(0), fDirectedFlowProtons(0.0), \r
-  fEllipticFlowProtons(0.0), fTriangularFlowProtons(0.0),\r
-  fQuandrangularFlowProtons(0.0), fPentangularFlowProtons(0.0) {\r
+AliAnalysisTaskToyModel::AliAnalysisTaskToyModel() \r
+: TObject(),\r
+    fUseDebug(kFALSE),\r
+    fBalance(0),\r
+    fRunShuffling(kFALSE), fShuffledBalance(0),\r
+    fList(0), fListBF(0), fListBFS(0),\r
+    fHistEventStats(0),\r
+    fHistNumberOfAcceptedParticles(0),\r
+    fHistReactionPlane(0),\r
+    fHistEtaTotal(0), fHistEta(0),\r
+    fHistRapidityTotal(0), fHistRapidity(0),\r
+    fHistRapidityPions(0), fHistRapidityKaons(0), fHistRapidityProtons(0),\r
+    fHistPhiTotal(0), fHistPhi(0),\r
+    fHistPhiPions(0), fHistPhiKaons(0), fHistPhiProtons(0),\r
+    fHistPtTotal(0), fHistPt(0),\r
+    fHistPtPions(0), fHistPtKaons(0), fHistPtProtons(0),\r
+    fTotalMultiplicityMean(100.), fTotalMultiplicitySigma(0.0),\r
+    fNetChargeMean(0.0), fNetChargeSigma(0.0),\r
+    fPtMin(0.0), fPtMax(100.0),\r
+    fEtaMin(-1.0), fEtaMax(1.0),\r
+    fUseAcceptanceParameterization(kFALSE), fAcceptanceParameterization(0),\r
+    fUseAllCharges(kFALSE), fParticleMass(0.0),\r
+    fPtSpectraAllCharges(0), fTemperatureAllCharges(100.),\r
+    fReactionPlane(0.0),\r
+    fAzimuthalAngleAllCharges(0), fDirectedFlowAllCharges(0.0), \r
+    fEllipticFlowAllCharges(0.0), fTriangularFlowAllCharges(0.0),\r
+    fQuandrangularFlowAllCharges(0.0), fPentangularFlowAllCharges(0.0),\r
+    fPionPercentage(0.8), fPionMass(0.0),\r
+    fPtSpectraPions(0), fTemperaturePions(100.),\r
+    fAzimuthalAnglePions(0), fDirectedFlowPions(0.0), \r
+    fEllipticFlowPions(0.0), fTriangularFlowPions(0.0), \r
+    fQuandrangularFlowPions(0.0), fPentangularFlowPions(0.0),\r
+    fKaonPercentage(0.8), fKaonMass(0.0),\r
+    fPtSpectraKaons(0), fTemperatureKaons(100.),\r
+    fAzimuthalAngleKaons(0), fDirectedFlowKaons(0.0), \r
+    fEllipticFlowKaons(0.0), fTriangularFlowKaons(0.0),\r
+    fQuandrangularFlowKaons(0.0), fPentangularFlowKaons(0.0),\r
+    fProtonPercentage(0.8), fProtonMass(0.0),\r
+    fPtSpectraProtons(0), fTemperatureProtons(100.),\r
+    fAzimuthalAngleProtons(0), fDirectedFlowProtons(0.0), \r
+    fEllipticFlowProtons(0.0), fTriangularFlowProtons(0.0),\r
+    fQuandrangularFlowProtons(0.0), fPentangularFlowProtons(0.0) {\r
   // Constructor\r
-\r
-  // Define input and output slots here\r
-  // Input slot #0 works with a TChain\r
-  DefineInput(0, TChain::Class());\r
-  // Output slot #0 writes into a TH1 container\r
-  DefineOutput(1, TList::Class());\r
-  DefineOutput(2, TList::Class());\r
-  DefineOutput(3, TList::Class());\r
 }\r
 \r
 //________________________________________________________________________\r
@@ -98,93 +100,73 @@ void AliAnalysisTaskToyModel::Init() {
 \r
   if(fUseAllCharges) {\r
     fParticleMass = fPionMass;\r
-    fPtSpectraAllCharges = new TF1("fPtSpectraAllCharges","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,100.);\r
+    fPtSpectraAllCharges = new TF1("fPtSpectraAllCharges","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);\r
     fPtSpectraAllCharges->SetParName(0,"Mass");\r
-    fPtSpectraAllCharges->SetParameter(0,fParticleMass);\r
     fPtSpectraAllCharges->SetParName(1,"Temperature");\r
-    fPtSpectraAllCharges->SetParameter(1,fTemperatureAllCharges);\r
   }\r
   else {\r
-    fPtSpectraPions = new TF1("fPtSpectraPions","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,100.);\r
+    fPtSpectraPions = new TF1("fPtSpectraPions","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);\r
     fPtSpectraPions->SetParName(0,"Mass");\r
-    fPtSpectraPions->SetParameter(0,fPionMass);\r
     fPtSpectraPions->SetParName(1,"Temperature");\r
-    fPtSpectraPions->SetParameter(1,fTemperaturePions);\r
     \r
-    fPtSpectraKaons = new TF1("fPtSpectraKaons","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,100.);\r
+    fPtSpectraKaons = new TF1("fPtSpectraKaons","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);\r
     fPtSpectraKaons->SetParName(0,"Mass");\r
-    fPtSpectraKaons->SetParameter(0,fKaonMass);\r
     fPtSpectraKaons->SetParName(1,"Temperature");\r
-    fPtSpectraKaons->SetParameter(1,fTemperatureKaons);\r
     \r
     fPtSpectraProtons = new TF1("fPtSpectraProtons","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);\r
     fPtSpectraProtons->SetParName(0,"Mass");\r
-    fPtSpectraProtons->SetParameter(0,fProtonMass);\r
     fPtSpectraProtons->SetParName(1,"Temperature");\r
-    fPtSpectraProtons->SetParameter(1,fTemperatureProtons);\r
   }\r
   //==============Particles and spectra==============//\r
 \r
   //==============Flow values==============//\r
   if(fUseAllCharges) {\r
+    if(fUseDebug) {\r
+      Printf("Directed flow: %lf",fDirectedFlowAllCharges);\r
+      Printf("Elliptic flow: %lf",fEllipticFlowAllCharges);\r
+      Printf("Triangular flow: %lf",fTriangularFlowAllCharges);\r
+      Printf("Quandrangular flow: %lf",fQuandrangularFlowAllCharges);\r
+      Printf("Pentangular flow: %lf",fPentangularFlowAllCharges);\r
+    }\r
+\r
     fAzimuthalAngleAllCharges = new TF1("fAzimuthalAngleAllCharges","1+2.*[1]*TMath::Cos(x-[0])+2.*[2]*TMath::Cos(2*(x-[0]))+2.*[3]*TMath::Cos(3*(x-[0]))+2.*[4]*TMath::Cos(4*(x-[0]))+2.*[5]*TMath::Cos(5*(x-[0]))",0.,2.*TMath::Pi());\r
     fAzimuthalAngleAllCharges->SetParName(0,"Reaction Plane");\r
     fAzimuthalAngleAllCharges->SetParName(1,"Directed flow");\r
-    fAzimuthalAngleAllCharges->SetParameter(1,fDirectedFlowAllCharges);\r
     fAzimuthalAngleAllCharges->SetParName(2,"Elliptic flow"); \r
-    fAzimuthalAngleAllCharges->SetParameter(2,fEllipticFlowAllCharges);\r
     fAzimuthalAngleAllCharges->SetParName(3,"Triangular flow");\r
-    fAzimuthalAngleAllCharges->SetParameter(3,fTriangularFlowAllCharges);\r
     fAzimuthalAngleAllCharges->SetParName(4,"Quandrangular flow");\r
-    fAzimuthalAngleAllCharges->SetParameter(4,fQuandrangularFlowAllCharges);\r
     fAzimuthalAngleAllCharges->SetParName(5,"Pentangular flow");\r
-    fAzimuthalAngleAllCharges->SetParameter(5,fPentangularFlowAllCharges);\r
   }\r
   else {\r
     fAzimuthalAnglePions = new TF1("fAzimuthalAnglePions","1+2.*[1]*TMath::Cos(x-[0])+2.*[2]*TMath::Cos(2*(x-[0]))+2.*[3]*TMath::Cos(3*(x-[0]))+2.*[4]*TMath::Cos(4*(x-[0]))+2.*[5]*TMath::Cos(5*(x-[0]))",0.,2.*TMath::Pi());\r
     fAzimuthalAnglePions->SetParName(0,"Reaction Plane");\r
     fAzimuthalAnglePions->SetParName(1,"Directed flow");\r
-    fAzimuthalAnglePions->SetParameter(1,fDirectedFlowPions);\r
     fAzimuthalAnglePions->SetParName(2,"Elliptic flow"); \r
-    fAzimuthalAnglePions->SetParameter(2,fEllipticFlowPions);\r
     fAzimuthalAnglePions->SetParName(3,"Triangular flow");\r
-    fAzimuthalAnglePions->SetParameter(3,fTriangularFlowPions);\r
     fAzimuthalAnglePions->SetParName(4,"Quandrangular flow");\r
-    fAzimuthalAnglePions->SetParameter(4,fQuandrangularFlowPions);\r
     fAzimuthalAnglePions->SetParName(5,"Pentangular flow");\r
-    fAzimuthalAnglePions->SetParameter(5,fPentangularFlowPions);\r
     \r
     fAzimuthalAngleKaons = new TF1("fAzimuthalAngleKaons","1+2.*[1]*TMath::Cos(x-[0])+2.*[2]*TMath::Cos(2*(x-[0]))+2.*[3]*TMath::Cos(3*(x-[0]))+2.*[4]*TMath::Cos(4*(x-[0]))+2.*[5]*TMath::Cos(5*(x-[0]))",0.,2.*TMath::Pi());\r
     fAzimuthalAngleKaons->SetParName(0,"Reaction Plane");\r
     fAzimuthalAngleKaons->SetParName(1,"Directed flow");\r
-    fAzimuthalAngleKaons->SetParameter(1,fDirectedFlowKaons);\r
     fAzimuthalAngleKaons->SetParName(2,"Elliptic flow"); \r
-    fAzimuthalAngleKaons->SetParameter(2,fEllipticFlowKaons);\r
     fAzimuthalAngleKaons->SetParName(3,"Triangular flow");\r
-    fAzimuthalAngleKaons->SetParameter(3,fTriangularFlowKaons);\r
     fAzimuthalAngleKaons->SetParName(4,"Quandrangular flow");\r
-    fAzimuthalAngleKaons->SetParameter(4,fQuandrangularFlowKaons);\r
     fAzimuthalAngleKaons->SetParName(5,"Pentangular flow");\r
-    fAzimuthalAngleKaons->SetParameter(5,fPentangularFlowKaons);\r
     \r
     fAzimuthalAngleProtons = new TF1("fAzimuthalAngleProtons","1+2.*[1]*TMath::Cos(x-[0])+2.*[2]*TMath::Cos(2*(x-[0]))+2.*[3]*TMath::Cos(3*(x-[0]))+2.*[4]*TMath::Cos(4*(x-[0]))+2.*[5]*TMath::Cos(5*(x-[0]))",0.,2.*TMath::Pi());\r
     fAzimuthalAngleProtons->SetParName(0,"Reaction Plane");\r
     fAzimuthalAngleProtons->SetParName(1,"Directed flow");\r
-    fAzimuthalAngleProtons->SetParameter(1,fDirectedFlowProtons);\r
     fAzimuthalAngleProtons->SetParName(2,"Elliptic flow"); \r
-    fAzimuthalAngleProtons->SetParameter(2,fEllipticFlowProtons);\r
     fAzimuthalAngleProtons->SetParName(3,"Triangular flow");\r
-    fAzimuthalAngleProtons->SetParameter(3,fTriangularFlowProtons);\r
     fAzimuthalAngleProtons->SetParName(4,"Quandrangular flow");\r
-    fAzimuthalAngleProtons->SetParameter(4,fQuandrangularFlowProtons);\r
     fAzimuthalAngleProtons->SetParName(5,"Pentangular flow");\r
-    fAzimuthalAngleProtons->SetParameter(5,fPentangularFlowProtons);\r
   }\r
   //==============Flow values==============//\r
 }\r
 \r
 //________________________________________________________________________\r
-void AliAnalysisTaskToyModel::UserCreateOutputObjects() {\r
+void AliAnalysisTaskToyModel::CreateOutputObjects() {\r
   // Create histograms\r
   // Called once\r
   if(!fBalance) {\r
@@ -214,6 +196,7 @@ void AliAnalysisTaskToyModel::UserCreateOutputObjects() {
     fListBFS->SetOwner();\r
   }\r
 \r
+  //==============QA================//\r
   //Event stats.\r
   TString gCutName[4] = {"Total","Offline trigger",\r
                          "Vertex","Analyzed"};\r
@@ -224,7 +207,62 @@ void AliAnalysisTaskToyModel::UserCreateOutputObjects() {
     fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());\r
   fList->Add(fHistEventStats);\r
 \r
-  // Balance function histograms\r
+  fHistNumberOfAcceptedParticles = new TH1F("fHistNumberOfAcceptedParticles",";N_{acc.};Entries",10000,0,10000);\r
+  fList->Add(fHistNumberOfAcceptedParticles);\r
+  \r
+  fHistReactionPlane = new TH1F("fHistReactionPlane","Reaction plane angle;#Psi [rad];Entries",1000,0.,2.*TMath::Pi());\r
+  fList->Add(fHistReactionPlane);\r
+\r
+  //Pseudo-rapidity\r
+  fHistEtaTotal = new TH1F("fHistEtaTotal","Pseudo-rapidity (full phase space);#eta;Entries",1000,-15.,15.); \r
+  fList->Add(fHistEtaTotal);\r
+\r
+  fHistEta = new TH1F("fHistEta","Pseudo-rapidity (acceptance);#eta;Entries",1000,-1.5,1.5); \r
+  fList->Add(fHistEta);\r
+\r
+  //Rapidity\r
+  fHistRapidityTotal = new TH1F("fHistRapidityTotal","Rapidity (full phase space);y;Entries",1000,-15.,15.); \r
+  fList->Add(fHistRapidityTotal);\r
+\r
+  fHistRapidity = new TH1F("fHistRapidity","Rapidity (acceptance);y;Entries",1000,-1.5,1.5); \r
+  fList->Add(fHistRapidity);\r
+  fHistRapidityPions = new TH1F("fHistRapidityPions","Rapidity (acceptance - pions);y;Entries",1000,-1.5,1.5); \r
+  fList->Add(fHistRapidityPions);\r
+  fHistRapidityKaons = new TH1F("fHistRapidityKaons","Rapidity (acceptance - kaons);y;Entries",1000,-1.5,1.5); \r
+  fList->Add(fHistRapidityKaons);\r
+  fHistRapidityProtons = new TH1F("fHistRapidityProtons","Rapidity (acceptance - protons);y;Entries",1000,-1.5,1.5); \r
+  fList->Add(fHistRapidityProtons);\r
+\r
+  //Phi\r
+  fHistPhiTotal = new TH1F("fHistPhiTotal","Phi (full phase space);#phi (rad);Entries",1000,0.,2*TMath::Pi());\r
+  fList->Add(fHistPhiTotal);\r
+\r
+  fHistPhi = new TH1F("fHistPhi","Phi (acceptance);#phi (rad);Entries",1000,0.,2*TMath::Pi());\r
+  fList->Add(fHistPhi);\r
+\r
+  fHistPhiPions = new TH1F("fHistPhiPions","Phi (acceptance - pions);#phi (rad);Entries",1000,0.,2*TMath::Pi());\r
+  fList->Add(fHistPhiPions);\r
+  fHistPhiKaons = new TH1F("fHistPhiKaons","Phi (acceptance - kaons);#phi (rad);Entries",1000,0.,2*TMath::Pi());\r
+  fList->Add(fHistPhiKaons);\r
+  fHistPhiProtons = new TH1F("fHistPhiProtons","Phi (acceptance - protons);#phi (rad);Entries",1000,0.,2*TMath::Pi());\r
+  fList->Add(fHistPhiProtons);\r
+\r
+\r
+  //Pt\r
+  fHistPtTotal = new TH1F("fHistPtTotal","Pt (full phase space);p_{t} (GeV/c);Entries",1000,0.,10.);\r
+  fList->Add(fHistPtTotal);\r
+\r
+  fHistPt = new TH1F("fHistPt","Pt (acceptance);p_{t} (GeV/c);Entries",1000,0.,10.);\r
+  fList->Add(fHistPt);\r
+\r
+  fHistPtPions = new TH1F("fHistPtPions","Pt (acceptance - pions);p_{t} (GeV/c);Entries",1000,0.,10.);\r
+  fList->Add(fHistPtPions);\r
+  fHistPtKaons = new TH1F("fHistPtKaons","Pt (acceptance - kaons);p_{t} (GeV/c);Entries",1000,0.,10.);\r
+  fList->Add(fHistPtKaons);\r
+  fHistPtProtons = new TH1F("fHistPtProtons","Pt (acceptance - protons);p_{t} (GeV/c);Entries",1000,0.,10.);\r
+  fList->Add(fHistPtProtons);\r
+\r
+  //==============Balance function histograms================//\r
   // Initialize histograms if not done yet\r
   if(!fBalance->GetHistNp(0)){\r
     AliWarning("Histograms not yet initialized! --> Will be done now");\r
@@ -259,24 +297,15 @@ void AliAnalysisTaskToyModel::UserCreateOutputObjects() {
   }\r
 \r
   // Post output data.\r
-  PostData(1, fList);\r
-  PostData(2, fListBF);\r
-  if(fRunShuffling) PostData(3, fListBFS);\r
+  //PostData(1, fList);\r
+  //PostData(2, fListBF);\r
+  //if(fRunShuffling) PostData(3, fListBFS);\r
 }\r
 \r
 //________________________________________________________________________\r
-void AliAnalysisTaskToyModel::UserExec(Option_t *) {\r
+void AliAnalysisTaskToyModel::Run(Int_t nEvents) {\r
   // Main loop\r
   // Called for each event\r
-\r
-  // vector holding the charges/kinematics of all tracks (charge,y,eta,phi,p0,p1,p2,pt,E)\r
-  vector<Double_t> *chargeVectorShuffle[9];   // this will be shuffled\r
-  vector<Double_t> *chargeVector[9];          // original charge\r
-  for(Int_t i = 0; i < 9; i++){\r
-    chargeVectorShuffle[i] = new vector<Double_t>;\r
-    chargeVector[i]        = new vector<Double_t>;\r
-  }\r
-\r
   Double_t v_charge;\r
   Double_t v_y;\r
   Double_t v_eta;\r
@@ -284,140 +313,236 @@ void AliAnalysisTaskToyModel::UserExec(Option_t *) {
   Double_t v_p[3];\r
   Double_t v_pt;\r
   Double_t v_E;\r
+  Bool_t isPion = kFALSE, isKaon = kFALSE, isProton = kFALSE;\r
 \r
-  //Multiplicities\r
-  Int_t nMultiplicity = (Int_t)(gRandom->Gaus(fTotalMultiplicityMean,fTotalMultiplicitySigma));\r
-  Int_t nNetCharge = (Int_t)(gRandom->Gaus(fNetChargeMean,fNetChargeSigma));\r
-  Int_t nGeneratedPositive = 0, nGeneratedNegative = 0;\r
-  Int_t nGeneratedPions = 0, nGeneratedKaons = 0, nGeneratedProtons = 0;\r
+  if(fUseAllCharges) {\r
+    fPtSpectraAllCharges->SetParameter(0,fParticleMass);\r
+    fPtSpectraAllCharges->SetParameter(1,fTemperatureAllCharges);\r
 \r
-  //Randomization of the reaction plane\r
-  fReactionPlane = 2.0*TMath::Pi()*gRandom->Rndm();\r
-  if(fUseAllCharges) \r
-    fAzimuthalAngleAllCharges->SetParameter(0,fReactionPlane);\r
+    fAzimuthalAngleAllCharges->SetParameter(1,fDirectedFlowAllCharges);\r
+    fAzimuthalAngleAllCharges->SetParameter(2,fEllipticFlowAllCharges);\r
+    fAzimuthalAngleAllCharges->SetParameter(3,fTriangularFlowAllCharges);\r
+    fAzimuthalAngleAllCharges->SetParameter(4,fQuandrangularFlowAllCharges);\r
+    fAzimuthalAngleAllCharges->SetParameter(5,fPentangularFlowAllCharges);\r
+  }\r
   else {\r
-    fAzimuthalAnglePions->SetParameter(0,fReactionPlane);\r
-    fAzimuthalAngleKaons->SetParameter(0,fReactionPlane);\r
-    fAzimuthalAngleProtons->SetParameter(0,fReactionPlane);\r
+    fPtSpectraPions->SetParameter(0,fPionMass);\r
+    fPtSpectraPions->SetParameter(1,fTemperaturePions);\r
+    fPtSpectraKaons->SetParameter(0,fKaonMass);\r
+    fPtSpectraKaons->SetParameter(1,fTemperatureKaons);\r
+    fPtSpectraProtons->SetParameter(0,fProtonMass);\r
+    fPtSpectraProtons->SetParameter(1,fTemperatureProtons);\r
+\r
+    fAzimuthalAnglePions->SetParameter(1,fDirectedFlowPions);\r
+    fAzimuthalAnglePions->SetParameter(2,fEllipticFlowPions);\r
+    fAzimuthalAnglePions->SetParameter(3,fTriangularFlowPions);\r
+    fAzimuthalAnglePions->SetParameter(4,fQuandrangularFlowPions);\r
+    fAzimuthalAnglePions->SetParameter(5,fPentangularFlowPions);\r
+\r
+    fAzimuthalAngleKaons->SetParameter(1,fDirectedFlowKaons);\r
+    fAzimuthalAngleKaons->SetParameter(2,fEllipticFlowKaons);\r
+    fAzimuthalAngleKaons->SetParameter(3,fTriangularFlowKaons);\r
+    fAzimuthalAngleKaons->SetParameter(4,fQuandrangularFlowKaons);\r
+    fAzimuthalAngleKaons->SetParameter(5,fPentangularFlowKaons);\r
+\r
+    fAzimuthalAngleProtons->SetParameter(1,fDirectedFlowProtons);\r
+    fAzimuthalAngleProtons->SetParameter(2,fEllipticFlowProtons);\r
+    fAzimuthalAngleProtons->SetParameter(3,fTriangularFlowProtons);\r
+    fAzimuthalAngleProtons->SetParameter(4,fQuandrangularFlowProtons);\r
+    fAzimuthalAngleProtons->SetParameter(5,fPentangularFlowProtons);\r
   }\r
 \r
-  Int_t gNumberOfAcceptedParticles = 0;\r
-  //Generate particles\r
-  for(Int_t iParticleCount = 0; iParticleCount < nMultiplicity; iParticleCount++) {\r
-    //Decide the charge\r
-    Double_t randomNumberCharge = gRandom->Rndm();\r
-    if(randomNumberCharge <= 0.5*(nMultiplicity + nNetCharge)/nMultiplicity) {\r
-      v_charge = 1.0;\r
-      nGeneratedPositive += 1;\r
+  for(Int_t iEvent = 0; iEvent < nEvents; iEvent++) {\r
+    // vector holding the charges/kinematics of all tracks (charge,y,eta,phi,p0,p1,p2,pt,E)\r
+    vector<Double_t> *chargeVectorShuffle[9];   // this will be shuffled\r
+    vector<Double_t> *chargeVector[9];          // original charge\r
+    for(Int_t i = 0; i < 9; i++){\r
+      chargeVectorShuffle[i] = new vector<Double_t>;\r
+      chargeVector[i]        = new vector<Double_t>;\r
     }\r
+\r
+    fHistEventStats->Fill(1);\r
+    fHistEventStats->Fill(2);\r
+    fHistEventStats->Fill(3);\r
+\r
+    if((iEvent%1000) == 0) \r
+      cout<<"Event: "<<iEvent<<"/"<<nEvents<<endl;\r
+    \r
+    //Multiplicities\r
+    Int_t nMultiplicity = (Int_t)(gRandom->Gaus(fTotalMultiplicityMean,fTotalMultiplicitySigma));\r
+    Int_t nNetCharge = (Int_t)(gRandom->Gaus(fNetChargeMean,fNetChargeSigma));\r
+    Int_t nGeneratedPositive = 0, nGeneratedNegative = 0;\r
+    Int_t nGeneratedPions = 0, nGeneratedKaons = 0, nGeneratedProtons = 0;\r
+    \r
+    //Randomization of the reaction plane\r
+    //fReactionPlane = 2.0*TMath::Pi()*gRandom->Rndm();\r
+    fReactionPlane = 0.0;\r
+    if(fUseAllCharges) \r
+      fAzimuthalAngleAllCharges->SetParameter(0,fReactionPlane);\r
     else {\r
-      nGeneratedNegative += 1;\r
-      v_charge = -1.0;\r
+      fAzimuthalAnglePions->SetParameter(0,fReactionPlane);\r
+      fAzimuthalAngleKaons->SetParameter(0,fReactionPlane);\r
+      fAzimuthalAngleProtons->SetParameter(0,fReactionPlane);\r
     }\r
-\r
-    if(!fUseAllCharges) {\r
-      //Decide the specie\r
-      Double_t randomNumberSpecies = gRandom->Rndm();\r
-      if((randomNumberSpecies >= 0.0)&&(randomNumberSpecies < fPionPercentage)) {\r
-       nGeneratedPions += 1;\r
-       v_pt = fPtSpectraPions->GetRandom();\r
-       v_phi = fAzimuthalAnglePions->GetRandom();\r
-       fParticleMass = fPionMass;\r
+    \r
+    Int_t gNumberOfAcceptedParticles = 0;\r
+    //Generate particles\r
+    for(Int_t iParticleCount = 0; iParticleCount < nMultiplicity; iParticleCount++) {\r
+      isPion = kFALSE; isKaon = kFALSE; isProton = kFALSE;\r
+\r
+      //Decide the charge\r
+      Double_t randomNumberCharge = gRandom->Rndm();\r
+      if(randomNumberCharge <= 0.5*(nMultiplicity + nNetCharge)/nMultiplicity) {\r
+       v_charge = 1.0;\r
+       nGeneratedPositive += 1;\r
       }\r
-      else if((randomNumberSpecies >= fPionPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage)) {\r
-       nGeneratedKaons += 1;\r
-       v_pt = fPtSpectraKaons->GetRandom();\r
-       v_phi = fAzimuthalAngleKaons->GetRandom();\r
-       fParticleMass = fKaonMass;\r
+      else {\r
+       nGeneratedNegative += 1;\r
+       v_charge = -1.0;\r
       }\r
-      else if((randomNumberSpecies >= fPionPercentage + fKaonPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage + fProtonPercentage)) {\r
-       nGeneratedProtons += 1;\r
-       v_pt = fPtSpectraProtons->GetRandom();\r
-       v_phi = fAzimuthalAngleProtons->GetRandom();\r
-       fParticleMass = fProtonMass;\r
+      \r
+      if(!fUseAllCharges) {\r
+       //Decide the specie\r
+       Double_t randomNumberSpecies = gRandom->Rndm();\r
+       if((randomNumberSpecies >= 0.0)&&(randomNumberSpecies < fPionPercentage)) {\r
+         nGeneratedPions += 1;\r
+         v_pt = fPtSpectraPions->GetRandom();\r
+         v_phi = fAzimuthalAnglePions->GetRandom();\r
+         fParticleMass = fPionMass;\r
+         isPion = kTRUE;\r
+       }\r
+       else if((randomNumberSpecies >= fPionPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage)) {\r
+         nGeneratedKaons += 1;\r
+         v_pt = fPtSpectraKaons->GetRandom();\r
+         v_phi = fAzimuthalAngleKaons->GetRandom();\r
+         fParticleMass = fKaonMass;\r
+         isKaon = kTRUE;\r
+       }\r
+       else if((randomNumberSpecies >= fPionPercentage + fKaonPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage + fProtonPercentage)) {\r
+         nGeneratedProtons += 1;\r
+         v_pt = fPtSpectraProtons->GetRandom();\r
+         v_phi = fAzimuthalAngleProtons->GetRandom();\r
+         fParticleMass = fProtonMass;\r
+         isProton = kTRUE;\r
+       }\r
       }\r
-    }\r
-    else {\r
-      v_pt = fPtSpectraAllCharges->GetRandom();\r
-      v_phi = fAzimuthalAngleAllCharges->GetRandom();\r
-    }\r
-    v_eta = gRandom->Gaus(0.0,4.0);\r
-\r
-    v_p[0] = v_pt*TMath::Cos(v_phi);\r
-    v_p[1] = v_pt*TMath::Sin(v_phi);\r
-    v_p[2] = v_pt*TMath::SinH(v_eta);\r
-    v_E = TMath::Sqrt(TMath::Power(fParticleMass,2) +\r
-                     TMath::Power(v_p[0],2) +\r
-                     TMath::Power(v_p[1],2) +\r
-                     TMath::Power(v_p[2],2));\r
-    \r
-    v_y = 0.5*TMath::Log((v_E + v_p[2])/(v_E - v_p[2]));\r
-\r
-    //Acceptance\r
-    if((v_eta < fEtaMin) || (v_eta > fEtaMax)) continue;\r
-    if((v_pt < fPtMin) || (v_pt > fPtMax)) continue;\r
-\r
-    // fill charge vector\r
-    chargeVector[0]->push_back(v_charge);\r
-    chargeVector[1]->push_back(v_y);\r
-    chargeVector[2]->push_back(v_eta);\r
-    chargeVector[3]->push_back(v_phi);\r
-    chargeVector[4]->push_back(v_p[0]);\r
-    chargeVector[5]->push_back(v_p[1]);\r
-    chargeVector[6]->push_back(v_p[2]);\r
-    chargeVector[7]->push_back(v_pt);\r
-    chargeVector[8]->push_back(v_E);\r
+      else {\r
+       v_pt = fPtSpectraAllCharges->GetRandom();\r
+       v_phi = fAzimuthalAngleAllCharges->GetRandom();\r
+      }\r
+      v_eta = gRandom->Gaus(0.0,4.0);\r
+      \r
+      v_p[0] = v_pt*TMath::Cos(v_phi);\r
+      v_p[1] = v_pt*TMath::Sin(v_phi);\r
+      v_p[2] = v_pt*TMath::SinH(v_eta);\r
+      v_E = TMath::Sqrt(TMath::Power(fParticleMass,2) +\r
+                       TMath::Power(v_p[0],2) +\r
+                       TMath::Power(v_p[1],2) +\r
+                       TMath::Power(v_p[2],2));\r
+      \r
+      v_y = 0.5*TMath::Log((v_E + v_p[2])/(v_E - v_p[2]));\r
+      \r
+      //Fill QA histograms (full phase space)\r
+      fHistEtaTotal->Fill(v_eta);\r
+      fHistRapidityTotal->Fill(v_y);\r
+      fHistPhiTotal->Fill(v_phi);\r
+      fHistPtTotal->Fill(v_pt);\r
+\r
+      //Acceptance\r
+      if((v_eta < fEtaMin) || (v_eta > fEtaMax)) continue;\r
+      if((v_pt < fPtMin) || (v_pt > fPtMax)) continue;\r
+      //Printf("pt: %lf - mins: %lf - max: %lf",v_pt,fPtMin,fPtMax);\r
+\r
+      //Fill QA histograms (acceptance);\r
+      fHistEta->Fill(v_eta);\r
+      fHistRapidity->Fill(v_y);\r
+      fHistPhi->Fill(v_phi);\r
+      fHistPt->Fill(v_pt);\r
+      if(isPion) {\r
+       fHistRapidityPions->Fill(v_y);\r
+       fHistPhiPions->Fill(v_phi);\r
+       fHistPtPions->Fill(v_pt);\r
+      }\r
+      else if(isKaon) {\r
+       fHistRapidityKaons->Fill(v_y);\r
+       fHistPhiKaons->Fill(v_phi);\r
+       fHistPtKaons->Fill(v_pt);\r
+      }\r
+      else if(isProton) {\r
+       fHistRapidityProtons->Fill(v_y);\r
+       fHistPhiProtons->Fill(v_phi);\r
+       fHistPtProtons->Fill(v_pt);\r
+      }\r
+\r
+      // fill charge vector\r
+      chargeVector[0]->push_back(v_charge);\r
+      chargeVector[1]->push_back(v_y);\r
+      chargeVector[2]->push_back(v_eta);\r
+      chargeVector[3]->push_back(v_phi);\r
+      chargeVector[4]->push_back(v_p[0]);\r
+      chargeVector[5]->push_back(v_p[1]);\r
+      chargeVector[6]->push_back(v_p[2]);\r
+      chargeVector[7]->push_back(v_pt);\r
+      chargeVector[8]->push_back(v_E);\r
+      \r
+      if(fRunShuffling) {\r
+       chargeVectorShuffle[0]->push_back(v_charge);\r
+       chargeVectorShuffle[1]->push_back(v_y);\r
+       chargeVectorShuffle[2]->push_back(v_eta);\r
+       chargeVectorShuffle[3]->push_back(v_phi);\r
+       chargeVectorShuffle[4]->push_back(v_p[0]);\r
+       chargeVectorShuffle[5]->push_back(v_p[1]);\r
+       chargeVectorShuffle[6]->push_back(v_p[2]);\r
+       chargeVectorShuffle[7]->push_back(v_pt);\r
+       chargeVectorShuffle[8]->push_back(v_E);\r
+      }\r
+      gNumberOfAcceptedParticles += 1;\r
+      \r
+    }//generated particle loop\r
     \r
-    if(fRunShuffling) {\r
-      chargeVectorShuffle[0]->push_back(v_charge);\r
-      chargeVectorShuffle[1]->push_back(v_y);\r
-      chargeVectorShuffle[2]->push_back(v_eta);\r
-      chargeVectorShuffle[3]->push_back(v_phi);\r
-      chargeVectorShuffle[4]->push_back(v_p[0]);\r
-      chargeVectorShuffle[5]->push_back(v_p[1]);\r
-      chargeVectorShuffle[6]->push_back(v_p[2]);\r
-      chargeVectorShuffle[7]->push_back(v_pt);\r
-      chargeVectorShuffle[8]->push_back(v_E);\r
+    if(fUseDebug) {\r
+      Printf("=======================================================");\r
+      Printf("Total: %d - Total positive: %d - Total negative: %d",nMultiplicity,nGeneratedPositive,nGeneratedNegative);\r
+      if(!fUseAllCharges)\r
+       Printf("Pions: %lf - Kaons: %lf - Protons: %lf",1.*nGeneratedPions/nMultiplicity,1.*nGeneratedKaons/nMultiplicity,1.*nGeneratedProtons/nMultiplicity);\r
+      Printf("Calculating the balance function for %d particles",chargeVector[0]->size());\r
     }\r
-    gNumberOfAcceptedParticles += 1;\r
-           \r
-  }//generated particle loop\r
-\r
-  Printf("=======================================================");\r
-  Printf("Total: %d - Total positive: %d - Total negative: %d",nMultiplicity,nGeneratedPositive,nGeneratedNegative);\r
-  Printf("Pions: %lf - Kaons: %lf - Protons: %lf",1.*nGeneratedPions/nMultiplicity,1.*nGeneratedKaons/nMultiplicity,1.*nGeneratedProtons/nMultiplicity);\r
 \r
-  //Calculate the balance function\r
-  fBalance->CalculateBalance(gNumberOfAcceptedParticles,chargeVector);\r
+    fHistEventStats->Fill(4);\r
+    fHistNumberOfAcceptedParticles->Fill(gNumberOfAcceptedParticles);\r
+    fHistReactionPlane->Fill(fReactionPlane);\r
 \r
-  if(fRunShuffling) {\r
-    // shuffle charges\r
-    random_shuffle( chargeVectorShuffle[0]->begin(), chargeVectorShuffle[0]->end() );\r
-    fShuffledBalance->CalculateBalance(gNumberOfAcceptedParticles,chargeVectorShuffle);\r
+    //Calculate the balance function\r
+    fBalance->CalculateBalance(gNumberOfAcceptedParticles,chargeVector);\r
+    \r
+    if(fRunShuffling) {\r
+      // shuffle charges\r
+      random_shuffle( chargeVectorShuffle[0]->begin(), chargeVectorShuffle[0]->end() );\r
+      fShuffledBalance->CalculateBalance(gNumberOfAcceptedParticles,chargeVectorShuffle);\r
+    }\r
   }\r
 }      \r
 \r
 //________________________________________________________________________\r
-void  AliAnalysisTaskToyModel::FinishTaskOutput() {\r
+void  AliAnalysisTaskToyModel::FinishOutput() {\r
   //Printf("END BF");\r
 \r
+  TFile *gOutput = TFile::Open("outputToyModel.root","recreate");\r
+  fList->Write();\r
+\r
   if (!fBalance) {\r
     Printf("ERROR: fBalance not available");\r
     return;\r
   }  \r
+  fListBF->Write();\r
+\r
   if(fRunShuffling) {\r
     if (!fShuffledBalance) {\r
       Printf("ERROR: fShuffledBalance not available");\r
       return;\r
     }\r
+    fListBFS->Write();\r
   }\r
+  gOutput->Close();\r
 }\r
 \r
-//________________________________________________________________________\r
-void AliAnalysisTaskToyModel::Terminate(Option_t *) {\r
-  // Draw result to the screen\r
-  // Called once at the end of the query\r
-\r
-  // not implemented ...\r
-\r
-}\r
index 31a2aa4..0348de2 100755 (executable)
@@ -12,21 +12,18 @@ class TF1;
 \r
 class AliBalance;\r
 \r
-#include "AliAnalysisTaskSE.h"\r
-#include "AliBalance.h"\r
-\r
-\r
-class AliAnalysisTaskToyModel : public AliAnalysisTaskSE {\r
+class AliAnalysisTaskToyModel : public TObject {\r
  public:\r
-  AliAnalysisTaskToyModel(const char *name = "AliAnalysisTaskToyModel");\r
+  AliAnalysisTaskToyModel();\r
   virtual ~AliAnalysisTaskToyModel(); \r
   \r
   virtual void   Init();\r
-  virtual void   UserCreateOutputObjects();\r
-  virtual void   UserExec(Option_t *option);\r
-  virtual void   FinishTaskOutput();\r
-  virtual void   Terminate(Option_t *);\r
+  virtual void   CreateOutputObjects();\r
+  virtual void   Run(Int_t nEvents);\r
+  virtual void   FinishOutput();\r
   \r
+  void SetDebugFlag() {fUseDebug = kTRUE;}\r
+\r
   void SetAnalysisObject(AliBalance *const analysis) {\r
     fBalance         = analysis;\r
   }\r
@@ -67,7 +64,7 @@ class AliAnalysisTaskToyModel : public AliAnalysisTaskSE {
   void SetTriangularFlowForAllCharges(Double_t v3) {\r
     fUseAllCharges = kTRUE;\r
     fTriangularFlowAllCharges = v3;}\r
-  void SetQuadrangularFlowForAllCharges(Double_t v4) {\r
+  void SetQuandrangularFlowForAllCharges(Double_t v4) {\r
     fUseAllCharges = kTRUE;\r
     fQuandrangularFlowAllCharges = v4;}\r
   void SetPentangularFlowForAllCharges(Double_t v5) {\r
@@ -85,7 +82,7 @@ class AliAnalysisTaskToyModel : public AliAnalysisTaskSE {
     fEllipticFlowPions = v2;}\r
   void SetTriangularFlowForPions(Double_t v3) {\r
     fTriangularFlowPions = v3;}\r
-  void SetQuadrangularFlowForPions(Double_t v4) {\r
+  void SetQuandrangularFlowForPions(Double_t v4) {\r
     fQuandrangularFlowPions = v4;}\r
   void SetPentangularFlowForPions(Double_t v5) {\r
     fPentangularFlowPions = v5;}\r
@@ -101,7 +98,7 @@ class AliAnalysisTaskToyModel : public AliAnalysisTaskSE {
     fEllipticFlowKaons = v2;}\r
   void SetTriangularFlowForKaons(Double_t v3) {\r
     fTriangularFlowKaons = v3;}\r
-  void SetQuadrangularFlowForKaons(Double_t v4) {\r
+  void SetQuandrangularFlowForKaons(Double_t v4) {\r
     fQuandrangularFlowKaons = v4;}\r
   void SetPentangularFlowForKaons(Double_t v5) {\r
     fPentangularFlowKaons = v5;}\r
@@ -117,13 +114,15 @@ class AliAnalysisTaskToyModel : public AliAnalysisTaskSE {
     fEllipticFlowProtons = v2;}\r
   void SetTriangularFlowForProtons(Double_t v3) {\r
     fTriangularFlowProtons = v3;}\r
-  void SetQuadrangularFlowForProtons(Double_t v4) {\r
+  void SetQuandrangularFlowForProtons(Double_t v4) {\r
     fQuandrangularFlowProtons = v4;}\r
   void SetPentangularFlowForProtons(Double_t v5) {\r
     fPentangularFlowProtons = v5;}\r
   //============Toy model: List of setters============//\r
 \r
  private:\r
+  Bool_t fUseDebug; //Debug flag\r
+\r
   AliBalance *fBalance; //BF object\r
   Bool_t fRunShuffling;//run shuffling or not\r
   AliBalance *fShuffledBalance; //BF object (shuffled)\r
@@ -132,6 +131,25 @@ class AliAnalysisTaskToyModel : public AliAnalysisTaskSE {
   TList *fListBFS; //fList object\r
 \r
   TH1F *fHistEventStats; //event stats\r
+  TH1F *fHistNumberOfAcceptedParticles; //number of accepted particles\r
+  TH1F *fHistReactionPlane; //reaction plane angle\r
+  TH1F *fHistEtaTotal; //pseudo-rapidity (full phase space)\r
+  TH1F *fHistEta; //pseudo-rapidity (acceptance)\r
+  TH1F *fHistRapidityTotal; //rapidity (full phase space)\r
+  TH1F *fHistRapidity; //rapidity (acceptance)\r
+  TH1F *fHistRapidityPions; //rapidity (acceptance)\r
+  TH1F *fHistRapidityKaons; //rapidity (acceptance)\r
+  TH1F *fHistRapidityProtons; //rapidity (acceptance)\r
+  TH1F *fHistPhiTotal; //phi (full phase space)\r
+  TH1F *fHistPhi; //phi (acceptance)\r
+  TH1F *fHistPhiPions; //phi (acceptance)\r
+  TH1F *fHistPhiKaons; //phi (acceptance)\r
+  TH1F *fHistPhiProtons; //phi (acceptance)\r
+  TH1F *fHistPtTotal; //pt (full phase space)\r
+  TH1F *fHistPt; //pt (acceptance)\r
+  TH1F *fHistPtPions; //pt (acceptance)\r
+  TH1F *fHistPtKaons; //pt (acceptance)\r
+  TH1F *fHistPtProtons; //pt (acceptance)\r
 \r
   //Toy model input\r
   Double_t fTotalMultiplicityMean; //mean for the total multiplicity\r
diff --git a/PWG2/EBYE/macros/runBalanceFunctionToyModel.C b/PWG2/EBYE/macros/runBalanceFunctionToyModel.C
new file mode 100755 (executable)
index 0000000..f9e5c5c
--- /dev/null
@@ -0,0 +1,118 @@
+//=========Total multiplicity=========//\r
+Double_t nTotalMultiplicityMean = 100.;\r
+Double_t nTotalMultiplicitySigma = 10.;\r
+\r
+//=========Net charge=========//\r
+Double_t nNetChargeMean = 0.0;\r
+Double_t nNetChargeSigma = 5.0;\r
+\r
+//==============Particles and spectra==============//\r
+Double_t gAllChargesTemperature = 0.11; //in GeV\r
+Double_t gPionPercentage = 0.8;\r
+Double_t gPionTemperature = 0.1; //in GeV\r
+Double_t gKaonPercentage = 0.12;\r
+Double_t gKaonTemperature = 0.12; //in GeV\r
+Double_t gProtonPercentage = 0.08;\r
+Double_t gProtonTemperature = 0.2; //in GeV\r
+//==============Particles and spectra==============//\r
+\r
+//==============Flow values==============//\r
+Double_t gDirectedFlow = 0.0;\r
+Double_t gEllipticFlow = 0.07;\r
+Double_t gTriangularFlow = 0.0;\r
+Double_t gQuandrangularFlow = 0.0;\r
+Double_t gPentangularFlow = 0.0;\r
+//==============Flow values==============//\r
+\r
+//=========Acceptance definition=========//\r
+Double_t gEtaMin = -1.0;\r
+Double_t gEtaMax = 1.0;\r
+Double_t gPtMin = 0.1;\r
+Double_t gPtMax = 2.0;\r
+//=========Acceptance definition=========//\r
+\r
+Bool_t kUseDebug = kFALSE;\r
+\r
+// Run macro used for the toy model analysis\r
+// Author: Panos.Christakoglou@nikhef.nl\r
+\r
+//______________________________________________________________________________\r
+void runBalanceFunctionToyModel(Int_t nEvents = 10000,\r
+                               Bool_t kUseAllCharges = kTRUE) {\r
+  TStopwatch timer;\r
+  timer.Start();\r
+\r
+  // load libraries\r
+  gSystem->Load("libCore.so");        \r
+  gSystem->Load("libGeom.so");\r
+  gSystem->Load("libVMC.so");\r
+  gSystem->Load("libPhysics.so");\r
+  gSystem->Load("libTree.so");\r
+  gSystem->Load("libSTEERBase.so");\r
+  gSystem->Load("libESD.so");\r
+  gSystem->Load("libAOD.so");\r
+  gSystem->Load("libANALYSIS.so");\r
+  gSystem->Load("libANALYSISalice.so");\r
+  gSystem->Load("libPWG2ebye.so");\r
+  \r
+  //configure the bf objects\r
+  gROOT->LoadMacro("configBalanceFunctionAnalysis.C");\r
+  AliBalance *bf  = GetBalanceFunctionObject("MC");\r
+  AliBalance *bfs = GetBalanceFunctionObject("MC",kTRUE);\r
+  \r
+  //Configure the toy model object\r
+  AliAnalysisTaskToyModel *toyModelAnalysis = new AliAnalysisTaskToyModel();\r
+  if(kUseDebug) toyModelAnalysis->SetDebugFlag();\r
+  toyModelAnalysis->SetAnalysisObject(bf);\r
+  toyModelAnalysis->SetShufflingObject(bfs);\r
+  toyModelAnalysis->SetTotalMultiplicity(nTotalMultiplicityMean,nTotalMultiplicitySigma);\r
+  toyModelAnalysis->SetNetCharge(nNetChargeMean,nNetChargeSigma);\r
+  toyModelAnalysis->SetKinematicsCutsMC(gPtMin,gPtMax,gEtaMin,gEtaMax);\r
+\r
+  if(kUseAllCharges) {\r
+    toyModelAnalysis->SetSpectraTemperatureForAllCharges(gAllChargesTemperature);\r
+    toyModelAnalysis->SetDirectedFlowForAllCharges(gDirectedFlow);\r
+    toyModelAnalysis->SetEllipticFlowForAllCharges(gEllipticFlow);\r
+    toyModelAnalysis->SetTriangularFlowForAllCharges(gTriangularFlow);\r
+    toyModelAnalysis->SetQuandrangularFlowForAllCharges(gQuandrangularFlow);\r
+    toyModelAnalysis->SetPentangularFlowForAllCharges(gPentangularFlow);\r
+  }\r
+  else {\r
+    //Pions\r
+    toyModelAnalysis->SetPionPercentage(gPionPercentage);\r
+    toyModelAnalysis->SetSpectraTemperatureForPions(gPionTemperature);\r
+    toyModelAnalysis->SetDirectedFlowForPions(gDirectedFlow);\r
+    toyModelAnalysis->SetEllipticFlowForPions(gEllipticFlow);\r
+    toyModelAnalysis->SetTriangularFlowForPions(gTriangularFlow);\r
+    toyModelAnalysis->SetQuandrangularFlowForPions(gQuandrangularFlow);\r
+    toyModelAnalysis->SetPentangularFlowForPions(gPentangularFlow);\r
+\r
+    //Kaons\r
+    toyModelAnalysis->SetKaonPercentage(gKaonPercentage);\r
+    toyModelAnalysis->SetSpectraTemperatureForKaons(gKaonTemperature);\r
+    toyModelAnalysis->SetDirectedFlowForKaons(gDirectedFlow);\r
+    toyModelAnalysis->SetEllipticFlowForKaons(gEllipticFlow);\r
+    toyModelAnalysis->SetTriangularFlowForKaons(gTriangularFlow);\r
+    toyModelAnalysis->SetQuandrangularFlowForKaons(gQuandrangularFlow);\r
+    toyModelAnalysis->SetPentangularFlowForKaons(gPentangularFlow);\r
+\r
+    //Protons\r
+    toyModelAnalysis->SetProtonPercentage(gProtonPercentage);\r
+    toyModelAnalysis->SetSpectraTemperatureForProtons(gProtonTemperature);\r
+    toyModelAnalysis->SetDirectedFlowForProtons(gDirectedFlow);\r
+    toyModelAnalysis->SetEllipticFlowForProtons(gEllipticFlow);\r
+    toyModelAnalysis->SetTriangularFlowForProtons(gTriangularFlow);\r
+    toyModelAnalysis->SetQuandrangularFlowForProtons(gQuandrangularFlow);\r
+    toyModelAnalysis->SetPentangularFlowForProtons(gPentangularFlow);\r
+  }\r
+\r
+  toyModelAnalysis->Init();\r
+  toyModelAnalysis->CreateOutputObjects();\r
+  toyModelAnalysis->Run(nEvents);\r
+  toyModelAnalysis->FinishOutput();\r
+\r
+  // Print real and CPU time used for analysis:  \r
+  timer.Stop();\r
+  timer.Print();\r
+}\r
+\r