end-of-line normalization
[u/mrichter/AliRoot.git] / PWGCF / EBYE / BalanceFunctions / AliAnalysisTaskToyModel.cxx
index 99ba00c..8985f57 100755 (executable)
-#include "TChain.h"\r
-#include "TList.h"\r
-#include "TCanvas.h"\r
-#include "TParticle.h"\r
-#include "TLorentzVector.h"\r
-#include "TGraphErrors.h"\r
-#include "TH1.h"\r
-#include "TH2.h"\r
-#include "TH3.h"\r
-#include "TArrayF.h"\r
-#include "TF1.h"\r
-#include "TRandom.h"\r
-#include "TFile.h"\r
-\r
-#include "AliAnalysisManager.h"\r
-#include "AliLog.h"\r
-\r
-#include "AliEventPoolManager.h"           \r
-\r
-#include "AliAnalysisTaskToyModel.h"\r
-#include "AliBalance.h"\r
-#include "AliBalancePsi.h"\r
-#include "AliAnalysisTaskTriggeredBF.h"\r
-\r
-\r
-// Analysis task for the toy model analysis\r
-// Authors: Panos.Christakoglou@nikhef.nl\r
-\r
-using std::cout;\r
-using std::endl;\r
-\r
-ClassImp(AliAnalysisTaskToyModel)\r
-\r
-//________________________________________________________________________\r
-AliAnalysisTaskToyModel::AliAnalysisTaskToyModel() \r
-: TObject(),\r
-  fUseDebug(kFALSE),\r
-  fBalance(0),\r
-  fRunShuffling(kFALSE), fShuffledBalance(0),\r
-  fRunMixing(kFALSE), fMixedBalance(0), fPoolMgr(0),\r
-  fList(0), fListBF(0), fListBFS(0), fListBFM(0),\r
-  fHistEventStats(0),\r
-  fHistNumberOfAcceptedParticles(0),\r
-  fHistReactionPlane(0),\r
-  fHistEtaTotal(0), fHistEta(0),\r
-  fHistEtaPhiPos(0), fHistEtaPhiNeg(0),\r
-  fHistRapidity(0),\r
-  fHistRapidityPions(0), fHistRapidityKaons(0), fHistRapidityProtons(0),\r
-  fHistPhi(0),\r
-  fHistPhiPions(0), fHistPhiKaons(0), fHistPhiProtons(0),\r
-  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
-  fSimulateDetectorEffects(kFALSE),\r
-  fNumberOfInefficientSectors(0),\r
-  fInefficiencyFactorInPhi(1.0),\r
-  fNumberOfDeadSectors(0),\r
-  fEfficiencyDropNearEtaEdges(kFALSE),\r
-  fEfficiencyMatrix(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
-  fUseDynamicalCorrelations(kFALSE), fDynamicalCorrelationsPercentage(0.1),\r
-  fUseJets(kFALSE), fPtAssoc(0) {\r
-  // Constructor\r
-}\r
-\r
-//________________________________________________________________________\r
-AliAnalysisTaskToyModel::~AliAnalysisTaskToyModel() {\r
-  //Destructor\r
-  if(fUseAllCharges) {\r
-    delete fPtSpectraAllCharges;\r
-    delete fAzimuthalAngleAllCharges;\r
-  }\r
-  else {\r
-    delete fPtSpectraPions;\r
-    delete fAzimuthalAnglePions;\r
-    delete fPtSpectraKaons;\r
-    delete fAzimuthalAngleKaons;\r
-    delete fPtSpectraProtons;\r
-    delete fAzimuthalAngleProtons;\r
-  }\r
-  if(fUseJets) delete fPtAssoc;\r
-}\r
-\r
-//________________________________________________________________________\r
-void AliAnalysisTaskToyModel::Init() {\r
-  //Initialize objects\r
-  //==============gRandom Seed=======================//\r
-  gRandom->SetSeed(0);  //seed is set to a random value (depending on machine clock)\r
-  //==============gRandom Seed=======================//\r
-\r
-  //==============Particles and spectra==============//\r
-  TParticle *pion = new TParticle();\r
-  pion->SetPdgCode(211);\r
-  fPionMass = pion->GetMass();\r
-\r
-  TParticle *kaon = new TParticle();\r
-  kaon->SetPdgCode(321);\r
-  fKaonMass = kaon->GetMass();\r
-  \r
-  TParticle *proton = new TParticle();\r
-  proton->SetPdgCode(2212);\r
-  fProtonMass = proton->GetMass();\r
-\r
-  if(fUseAllCharges) {\r
-    fParticleMass = fPionMass;\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->SetParName(1,"Temperature");\r
-    //fPtSpectraAllCharges = new TF1("fPtSpectraAllCharges","(x^2/TMath::Sqrt(TMath::Power(x,2) + TMath::Power(0.139,2)))*TMath::Power((1. + x/[0]),-[1])",0.,20.);\r
-    //fPtSpectraAllCharges->SetParName(0,"pt0");\r
-    //fPtSpectraAllCharges->SetParName(1,"b");\r
-  }\r
-  else {\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->SetParName(1,"Temperature");\r
-    \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->SetParName(1,"Temperature");\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->SetParName(1,"Temperature");\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->SetParName(2,"Elliptic flow"); \r
-    fAzimuthalAngleAllCharges->SetParName(3,"Triangular flow");\r
-    fAzimuthalAngleAllCharges->SetParName(4,"Quandrangular flow");\r
-    fAzimuthalAngleAllCharges->SetParName(5,"Pentangular flow");\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->SetParName(2,"Elliptic flow"); \r
-    fAzimuthalAnglePions->SetParName(3,"Triangular flow");\r
-    fAzimuthalAnglePions->SetParName(4,"Quandrangular flow");\r
-    fAzimuthalAnglePions->SetParName(5,"Pentangular flow");\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->SetParName(2,"Elliptic flow"); \r
-    fAzimuthalAngleKaons->SetParName(3,"Triangular flow");\r
-    fAzimuthalAngleKaons->SetParName(4,"Quandrangular flow");\r
-    fAzimuthalAngleKaons->SetParName(5,"Pentangular flow");\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->SetParName(2,"Elliptic flow"); \r
-    fAzimuthalAngleProtons->SetParName(3,"Triangular flow");\r
-    fAzimuthalAngleProtons->SetParName(4,"Quandrangular flow");\r
-    fAzimuthalAngleProtons->SetParName(5,"Pentangular flow");\r
-  }\r
-  //==============Flow values==============//\r
-\r
-  //===================Jets===================//\r
-  if(fUseJets) {\r
-    fPtAssoc = new TF1("fPtAssoc","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,20.);\r
-    fPtAssoc->SetParName(0,"pt0");\r
-    fPtAssoc->SetParName(1,"b");\r
-    fPtAssoc->SetParameter(0,0.139);\r
-    fPtAssoc->SetParameter(1,0.5);\r
-    fPtAssoc->SetLineColor(1);\r
-  }\r
-\r
-  //===================Jets===================//\r
-\r
-  //==============Efficiency matrix==============//\r
-  if(fSimulateDetectorEffects) SetupEfficiencyMatrix();\r
-  //==============Efficiency matrix==============//\r
-}\r
-\r
-//________________________________________________________________________\r
-void AliAnalysisTaskToyModel::SetupEfficiencyMatrix() {\r
-  //Setup the efficiency matrix\r
-  TH1F *hPt = new TH1F("hPt","",200,0.1,20.1);\r
-  TH1F *hEta = new TH1F("hEta","",20,-0.95,0.95);\r
-  TH1F *hPhi = new TH1F("hPhi","",72,0.,2.*TMath::Pi());\r
-  fEfficiencyMatrix = new TH3F("fEfficiencyMatrix","",\r
-                              hEta->GetNbinsX(),\r
-                              hEta->GetXaxis()->GetXmin(),\r
-                              hEta->GetXaxis()->GetXmax(),\r
-                              hPt->GetNbinsX(),\r
-                              hPt->GetXaxis()->GetXmin(),\r
-                              hPt->GetXaxis()->GetXmax(),\r
-                              hPhi->GetNbinsX(),\r
-                              hPhi->GetXaxis()->GetXmin(),\r
-                              hPhi->GetXaxis()->GetXmax());\r
-\r
-  //Efficiency in pt\r
-  Double_t epsilon[20] = {0.3,0.6,0.77,0.79,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80};\r
-  for(Int_t i=1;i<=20;i++) {\r
-    hPt->SetBinContent(i,epsilon[i-1]);\r
-    hPt->SetBinError(i,0.01);\r
-  }\r
-  for(Int_t i=21;i<=200;i++) {\r
-    hPt->SetBinContent(i,epsilon[19]);\r
-    hPt->SetBinError(i,0.01);\r
-  }\r
-\r
-  //Efficiency in eta\r
-  for(Int_t i=1;i<=hEta->GetNbinsX();i++) {\r
-    hEta->SetBinContent(i,1.0);\r
-    hEta->SetBinError(i,0.01);\r
-  }\r
-  if(fEfficiencyDropNearEtaEdges) {\r
-    hEta->SetBinContent(1,0.7); hEta->SetBinContent(2,0.8);\r
-    hEta->SetBinContent(3,0.9);\r
-    hEta->SetBinContent(18,0.9); hEta->SetBinContent(19,0.8);\r
-    hEta->SetBinContent(20,0.7);\r
-  }\r
-\r
-  //Efficiency in phi\r
-  for(Int_t i=1;i<=hPhi->GetNbinsX();i++) {\r
-    hPhi->SetBinContent(i,1.0);\r
-    hPhi->SetBinError(i,0.01);\r
-  }\r
-  for(Int_t i=1;i<=fNumberOfInefficientSectors;i++)\r
-    hPhi->SetBinContent(hPhi->FindBin(hPhi->GetRandom()),fInefficiencyFactorInPhi);\r
-  for(Int_t i=1;i<=fNumberOfDeadSectors;i++)\r
-    hPhi->SetBinContent(hPhi->FindBin(hPhi->GetRandom()),0.0);\r
-  \r
-  //Fill the 3D efficiency map\r
-  for(Int_t iBinX = 1; iBinX<=fEfficiencyMatrix->GetXaxis()->GetNbins();iBinX++) {\r
-    //cout<<"==================================="<<endl;\r
-    for(Int_t iBinY = 1; iBinY<=fEfficiencyMatrix->GetYaxis()->GetNbins();iBinY++) {\r
-      //cout<<"==================================="<<endl;\r
-      for(Int_t iBinZ = 1; iBinZ<=fEfficiencyMatrix->GetZaxis()->GetNbins();iBinZ++) {\r
-       fEfficiencyMatrix->SetBinContent(iBinX,iBinY,iBinZ,hEta->GetBinContent(iBinX)*hPt->GetBinContent(iBinY)*hPhi->GetBinContent(iBinZ));\r
-       //cout<<"Eta: "<<hEta->GetBinCenter(iBinX)<<" - Pt: "<<hPt->GetBinCenter(iBinY)<<" - Phi: "<<hPhi->GetBinCenter(iBinZ)<<" - "<<hEta->GetBinContent(iBinX)<<" , "<<hPt->GetBinContent(iBinY)<<" , "<<hPhi->GetBinContent(iBinZ)<<" - Efficiency: "<<hEta->GetBinContent(iBinX)*hPt->GetBinContent(iBinY)*hPhi->GetBinContent(iBinZ)<<endl;\r
-      }\r
-    }\r
-  }\r
-}\r
-\r
-//________________________________________________________________________\r
-void AliAnalysisTaskToyModel::CreateOutputObjects() {\r
-  // Create histograms\r
-  // Called once\r
-\r
-  // global switch disabling the reference \r
-  // (to avoid "Replacing existing TH1" if several wagons are created in train)\r
-  Bool_t oldStatus = TH1::AddDirectoryStatus();\r
-  TH1::AddDirectory(kFALSE);\r
-\r
-  if(!fBalance) {\r
-    fBalance = new AliBalancePsi();\r
-    fBalance->SetDeltaEtaMax(2.0);\r
-  }\r
-  if(fRunShuffling) {\r
-    if(!fShuffledBalance) {\r
-      fShuffledBalance = new AliBalancePsi();\r
-      fShuffledBalance->SetDeltaEtaMax(2.0);\r
-    }\r
-  }\r
-  if(fRunMixing) {\r
-    if(!fMixedBalance) {\r
-      fMixedBalance = new AliBalancePsi();\r
-      fMixedBalance->SetDeltaEtaMax(2.0);\r
-    }\r
-  }\r
-\r
-  //QA list\r
-  fList = new TList();\r
-  fList->SetName("listQA");\r
-  fList->SetOwner();\r
-\r
-  //Balance Function list\r
-  fListBF = new TList();\r
-  fListBF->SetName("listBF");\r
-  fListBF->SetOwner();\r
-\r
-  if(fRunShuffling) {\r
-    fListBFS = new TList();\r
-    fListBFS->SetName("listBFShuffled");\r
-    fListBFS->SetOwner();\r
-  }\r
-\r
-  if(fRunMixing) {\r
-    fListBFM = new TList();\r
-    fListBFM->SetName("listBFMixed");\r
-    fListBFM->SetOwner();\r
-  }\r
-\r
-  //==============QA================//\r
-  //Event stats.\r
-  TString gCutName[4] = {"Total","Offline trigger",\r
-                         "Vertex","Analyzed"};\r
-  fHistEventStats = new TH1F("fHistEventStats",\r
-                             "Event statistics;;N_{events}",\r
-                             4,0.5,4.5);\r
-  for(Int_t i = 1; i <= 4; i++)\r
-    fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());\r
-  fList->Add(fHistEventStats);\r
-\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
-  fHistEtaPhiPos = new TH2F("fHistEtaPhiPos","#eta-#phi distribution (+);#eta;#varphi (rad)",80,-2.,2.,72,-TMath::Pi()/2.,3.*TMath::Pi()/2.);\r
-  fList->Add(fHistEtaPhiPos);\r
-  fHistEtaPhiNeg = new TH2F("fHistEtaPhiNeg","#eta-#phi distribution (-);#eta;#varphi (rad)",80,-2.,2.,72,-TMath::Pi()/2.,3.*TMath::Pi()/2.);\r
-  fList->Add(fHistEtaPhiNeg);\r
-\r
-  //Rapidity\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
-  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
-  //Pt\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
-  if(fEfficiencyMatrix) fList->Add(fEfficiencyMatrix);\r
-\r
-  //==============Balance function histograms================//\r
-  // Initialize histograms if not done yet\r
-  if(!fBalance->GetHistNp()){\r
-    AliWarning("Histograms not yet initialized! --> Will be done now");\r
-    AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
-    fBalance->InitHistograms();\r
-  }\r
-\r
-  if(fRunShuffling) {\r
-    if(!fShuffledBalance->GetHistNp()) {\r
-      AliWarning("Histograms (shuffling) not yet initialized! --> Will be done now");\r
-      AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
-      fShuffledBalance->InitHistograms();\r
-    }\r
-  }\r
-\r
-  fListBF->Add(fBalance->GetHistNp());\r
-  fListBF->Add(fBalance->GetHistNn());\r
-  fListBF->Add(fBalance->GetHistNpn());\r
-  fListBF->Add(fBalance->GetHistNnn());\r
-  fListBF->Add(fBalance->GetHistNpp());\r
-  fListBF->Add(fBalance->GetHistNnp());\r
-\r
-  if(fRunShuffling) {\r
-    fListBFS->Add(fShuffledBalance->GetHistNp());\r
-    fListBFS->Add(fShuffledBalance->GetHistNn());\r
-    fListBFS->Add(fShuffledBalance->GetHistNpn());\r
-    fListBFS->Add(fShuffledBalance->GetHistNnn());\r
-    fListBFS->Add(fShuffledBalance->GetHistNpp());\r
-    fListBFS->Add(fShuffledBalance->GetHistNnp());\r
-  }\r
-\r
-  if(fRunMixing) {\r
-    fListBFM->Add(fMixedBalance->GetHistNp());\r
-    fListBFM->Add(fMixedBalance->GetHistNn());\r
-    fListBFM->Add(fMixedBalance->GetHistNpn());\r
-    fListBFM->Add(fMixedBalance->GetHistNnn());\r
-    fListBFM->Add(fMixedBalance->GetHistNpp());\r
-    fListBFM->Add(fMixedBalance->GetHistNnp());\r
-  }\r
-\r
-  // Event Mixing\r
-  if(fRunMixing){\r
-    Int_t fMixingTracks = 2000;\r
-    Int_t trackDepth = fMixingTracks; \r
-    Int_t poolsize   = 1000;  // Maximum number of events, ignored in the present implemented of AliEventPoolManager\r
-    \r
-    // centrality bins\r
-    // Double_t centralityBins[] = {0.,1.,2.,3.,4.,5.,7.,10.,20.,30.,40.,50.,60.,70.,80.,100.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
-    // Double_t* centbins        = centralityBins;\r
-    // Int_t nCentralityBins     = sizeof(centralityBins) / sizeof(Double_t) - 1;\r
-\r
-    // multiplicity bins\r
-    Double_t multiplicityBins[] = {0,10,20,30,40,50,60,70,80,100,100000}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
-    Double_t* multbins        = multiplicityBins;\r
-    Int_t nMultiplicityBins     = sizeof(multiplicityBins) / sizeof(Double_t) - 1;\r
-    \r
-    // Zvtx bins\r
-    Double_t vertexBins[] = {-10., -7., -5., -3., -1., 1., 3., 5., 7., 10.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
-    Double_t* vtxbins     = vertexBins;\r
-    Int_t nVertexBins     = sizeof(vertexBins) / sizeof(Double_t) - 1;\r
-    \r
-    // Event plane angle (Psi) bins\r
-    // Double_t psiBins[] = {0.,45.,135.,215.,305.,360.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
-    // Double_t* psibins     = psiBins;\r
-    // Int_t nPsiBins     = sizeof(psiBins) / sizeof(Double_t) - 1;\r
-    \r
-    // // run the event mixing also in bins of event plane (statistics!)\r
-    // if(fRunMixingEventPlane){\r
-    //   if(fEventClass=="Multiplicity"){\r
-    //         fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins, nPsiBins, psibins);\r
-    //   }\r
-    //   else{\r
-    //         fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins, nPsiBins, psibins);\r
-    //   }\r
-    // }\r
-    // else{\r
-    //if(fEventClass=="Multiplicity"){\r
-    fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins);\r
-    //}\r
-    //else{\r
-    //fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);\r
-    //}\r
-  }\r
-\r
-  TH1::AddDirectory(oldStatus);\r
-}\r
-\r
-//________________________________________________________________________\r
-void AliAnalysisTaskToyModel::Run(Int_t nEvents) {\r
-  // Main loop\r
-  // Called for each event\r
-  Short_t vCharge = 0;\r
-  Float_t vY = 0.0;\r
-  Float_t vEta = 0.0;\r
-  Float_t vPhi = 0.0;\r
-  Float_t vP[3] = {0.,0.,0.};\r
-  Float_t vPt = 0.0;\r
-  Float_t vE = 0.0;\r
-  Bool_t isPion = kFALSE, isKaon = kFALSE, isProton = kFALSE;\r
-\r
-  Double_t gDecideCharge = 0.;\r
-\r
-  if(fUseAllCharges) {\r
-    //fPtSpectraAllCharges->SetParameter(0,fParticleMass);\r
-    fPtSpectraAllCharges->SetParameter(0,1.05);\r
-    fPtSpectraAllCharges->SetParameter(1,fTemperatureAllCharges);\r
-\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
-    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
-\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
-    // TObjArray for the accepted particles \r
-    // (has to be done here, otherwise mxing with event pool does not work, overwriting pointers!)\r
-    TObjArray *tracksMain = new TObjArray();\r
-    tracksMain->SetOwner(kTRUE);\r
-    TObjArray *tracksMixing = 0x0;\r
-    if(fRunMixing) {\r
-      tracksMixing = new TObjArray();\r
-      tracksMixing->SetOwner(kTRUE);\r
-    }\r
-\r
-    tracksMain->Clear();\r
-    if(fRunMixing) tracksMixing->Clear();\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
-    \r
-    Int_t nGeneratedPositive = (Int_t)((nMultiplicity/2) + nNetCharge);\r
-    Int_t nGeneratedNegative = nMultiplicity - nGeneratedPositive;\r
-    if(fUseDebug) \r
-      Printf("Total multiplicity: %d - Generated positive: %d - Generated negative: %d",nMultiplicity,nGeneratedPositive,nGeneratedNegative);\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
-      fAzimuthalAnglePions->SetParameter(0,fReactionPlane);\r
-      fAzimuthalAngleKaons->SetParameter(0,fReactionPlane);\r
-      fAzimuthalAngleProtons->SetParameter(0,fReactionPlane);\r
-    }\r
-    \r
-    Int_t gNumberOfAcceptedParticles = 0;\r
-    Int_t gNumberOfAcceptedPositiveParticles = 0;\r
-    Int_t gNumberOfAcceptedNegativeParticles = 0;\r
-    Int_t nGeneratedPions = 0, nGeneratedKaons = 0, nGeneratedProtons = 0;\r
\r
-    //Generate particles\r
-    for(Int_t iParticleCount = 0; iParticleCount < nMultiplicity; iParticleCount++) {\r
-      isPion = kFALSE; isKaon = kFALSE; isProton = kFALSE;\r
-      if(fUseDebug) \r
-       Printf("Generating positive: %d(%d)",iParticleCount+1,nGeneratedPositive);\r
-\r
-      //Pseudo-rapidity sampled from a Gaussian centered @ 0\r
-      vEta = gRandom->Gaus(0.0,4.0);\r
-\r
-      //Fill QA histograms (full phase space)\r
-      fHistEtaTotal->Fill(vEta);\r
-\r
-      //Decide the charge\r
-      gDecideCharge = gRandom->Rndm();\r
-      if(gDecideCharge <= 1.*nGeneratedPositive/nMultiplicity)       \r
-       vCharge = 1;\r
-      else \r
-       vCharge = -1;\r
-      \r
-      //Acceptance\r
-      if((vEta < fEtaMin) || (vEta > fEtaMax)) continue;\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
-         vPt = fPtSpectraPions->GetRandom();\r
-         vPhi = fAzimuthalAnglePions->GetRandom();\r
-         fParticleMass = fPionMass;\r
-         isPion = kTRUE;\r
-       }\r
-       else if((randomNumberSpecies >= fPionPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage)) {\r
-         nGeneratedKaons += 1;\r
-         vPt = fPtSpectraKaons->GetRandom();\r
-         vPhi = fAzimuthalAngleKaons->GetRandom();\r
-         fParticleMass = fKaonMass;\r
-         isKaon = kTRUE;\r
-       }\r
-       else if((randomNumberSpecies >= fPionPercentage + fKaonPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage + fProtonPercentage)) {\r
-         nGeneratedProtons += 1;\r
-         vPt = fPtSpectraProtons->GetRandom();\r
-         vPhi = fAzimuthalAngleProtons->GetRandom();\r
-         fParticleMass = fProtonMass;\r
-         isProton = kTRUE;\r
-       }\r
-      }\r
-      else {\r
-       vPt = fPtSpectraAllCharges->GetRandom();\r
-       vPhi = fAzimuthalAngleAllCharges->GetRandom();\r
-      }\r
-      \r
-      vP[0] = vPt*TMath::Cos(vPhi);\r
-      vP[1] = vPt*TMath::Sin(vPhi);\r
-      vP[2] = vPt*TMath::SinH(vEta);\r
-      vE = TMath::Sqrt(TMath::Power(fParticleMass,2) +\r
-                       TMath::Power(vP[0],2) +\r
-                       TMath::Power(vP[1],2) +\r
-                       TMath::Power(vP[2],2));\r
-      \r
-      vY = 0.5*TMath::Log((vE + vP[2])/(vE - vP[2]));\r
-      \r
-      //pt coverage\r
-      if((vPt < fPtMin) || (vPt > fPtMax)) continue;\r
-      //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);\r
-\r
-      //acceptance filter\r
-      if(fUseAcceptanceParameterization) {\r
-       Double_t gRandomNumberForAcceptance = gRandom->Rndm();\r
-       if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(vPt)) \r
-         continue;\r
-      }\r
-\r
-      //Detector effects\r
-      if(fSimulateDetectorEffects) {\r
-       Double_t randomNumber = gRandom->Rndm();\r
-       if(randomNumber > fEfficiencyMatrix->GetBinContent(fEfficiencyMatrix->FindBin(vEta,vPt,vPhi)))\r
-         continue;\r
-      }\r
-\r
-      //Fill QA histograms (acceptance);\r
-      if(vCharge > 0) {\r
-       gNumberOfAcceptedPositiveParticles += 1;\r
-       if(vPhi > 3.*TMath::Pi()/2.)\r
-         fHistEtaPhiPos->Fill(vEta,vPhi-2.*TMath::Pi());\r
-       else\r
-         fHistEtaPhiPos->Fill(vEta,vPhi);\r
-      }\r
-      else {\r
-       gNumberOfAcceptedNegativeParticles += 1;\r
-       if(vPhi > 3.*TMath::Pi()/2.)\r
-         fHistEtaPhiNeg->Fill(vEta,vPhi-2.*TMath::Pi());\r
-       else\r
-         fHistEtaPhiNeg->Fill(vEta,vPhi);\r
-      }\r
-\r
-      fHistEta->Fill(vEta);\r
-      fHistRapidity->Fill(vY);\r
-      fHistPhi->Fill(vPhi);\r
-      fHistPt->Fill(vPt);\r
-      if(isPion) {\r
-       fHistRapidityPions->Fill(vY);\r
-       fHistPhiPions->Fill(vPhi);\r
-       fHistPtPions->Fill(vPt);\r
-      }\r
-      else if(isKaon) {\r
-       fHistRapidityKaons->Fill(vY);\r
-       fHistPhiKaons->Fill(vPhi);\r
-       fHistPtKaons->Fill(vPt);\r
-      }\r
-      else if(isProton) {\r
-       fHistRapidityProtons->Fill(vY);\r
-       fHistPhiProtons->Fill(vPhi);\r
-       fHistPtProtons->Fill(vPt);\r
-      }\r
-      gNumberOfAcceptedParticles += 1;\r
-\r
-      // add the track to the TObjArray\r
-      tracksMain->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, 1.0));  \r
-      if(fRunMixing) \r
-       tracksMixing->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, 1.0));  \r
-    }//generated positive particle loop\r
-\r
-    //Jets\r
-    if(fUseJets) {\r
-      const Int_t nAssociated = 3;\r
-\r
-      Double_t gPtTrig1 = 0., gPtTrig2 = 0., gPtAssoc = 0.;\r
-      Double_t gPhiTrig1 = 0., gPhiTrig2 = 0., gPhiAssoc = 0.;\r
-      Double_t gEtaTrig1 = 0., gEtaTrig2 = 0., gEtaAssoc = 0.;\r
-      Short_t gChargeTrig1 = 0, gChargeTrig2 = 0, gChargeAssoc = 0;\r
\r
-      Double_t gJetCone = 0.2;\r
-\r
-      //First leading particle\r
-      gPtTrig1 = gRandom->Uniform(3.,5.);\r
-      gEtaTrig1 = gRandom->Uniform(-0.8,0.8);\r
-      gPhiTrig1 = gRandom->Uniform(0.,TMath::TwoPi());\r
-\r
-      //Decide the charge\r
-      gDecideCharge = gRandom->Rndm();\r
-      if(gDecideCharge <= 0.5)\r
-       gChargeTrig1 = 1;\r
-      else \r
-       gChargeTrig1 = -1;\r
-      \r
-      //Acceptance\r
-      if((gEtaTrig1 < fEtaMin) || (gEtaTrig1 > fEtaMax)) continue;\r
-      //pt coverage\r
-      if((gPtTrig1 < fPtMin) || (gPtTrig1 > fPtMax)) continue;\r
-      //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);\r
-\r
-      //acceptance filter\r
-      if(fUseAcceptanceParameterization) {\r
-       Double_t gRandomNumberForAcceptance = gRandom->Rndm();\r
-       if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(gPtTrig1)) \r
-         continue;\r
-      }\r
-\r
-      //Detector effects\r
-      if(fSimulateDetectorEffects) {\r
-       Double_t randomNumber = gRandom->Rndm();\r
-       if(randomNumber > fEfficiencyMatrix->GetBinContent(fEfficiencyMatrix->FindBin(gEtaTrig1,gPtTrig1,gPhiTrig1)))\r
-         continue;\r
-      }\r
-\r
-      gNumberOfAcceptedParticles += 1;\r
-\r
-      // add the track to the TObjArray\r
-      tracksMain->Add(new AliBFBasicParticle(gEtaTrig1, gPhiTrig1, gPtTrig1, gChargeTrig1, 1.0));  \r
-      if(fRunMixing) \r
-       tracksMixing->Add(new AliBFBasicParticle(gEtaTrig1, gPhiTrig1, gPtTrig1, gChargeTrig1, 1.0));  \r
-\r
-      Int_t iAssociated = 0; \r
-      while(iAssociated < nAssociated) {\r
-       gPtAssoc = fPtAssoc->GetRandom();\r
-       if(gPtAssoc < gPtTrig1) {\r
-         gEtaAssoc = gRandom->Uniform(gEtaTrig1 - gJetCone/2.,gEtaTrig1 + gJetCone/2.);\r
-         gPhiAssoc = gRandom->Uniform(gPhiTrig1 - gJetCone/2.,gPhiTrig1 + gJetCone/2.);\r
-         if(gPhiAssoc < 0.) gPhiAssoc += TMath::TwoPi();\r
-         else if(gPhiAssoc > TMath::TwoPi()) gPhiAssoc -= TMath::TwoPi();\r
-         \r
-         iAssociated += 1;\r
-\r
-         gDecideCharge = gRandom->Rndm();\r
-         if(gDecideCharge <= 0.5)\r
-           gChargeAssoc = 1;\r
-         else \r
-           gChargeAssoc = -1;\r
-         \r
-         //Acceptance\r
-         if((gEtaAssoc < fEtaMin) || (gEtaAssoc > fEtaMax)) continue;\r
-         //pt coverage\r
-         if((gPtAssoc < fPtMin) || (gPtAssoc > fPtMax)) continue;\r
-         //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);\r
-\r
-         //acceptance filter\r
-         if(fUseAcceptanceParameterization) {\r
-           Double_t gRandomNumberForAcceptance = gRandom->Rndm();\r
-           if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(gPtAssoc)) \r
-             continue;\r
-         }\r
-\r
-         //Detector effects\r
-         if(fSimulateDetectorEffects) {\r
-           Double_t randomNumber = gRandom->Rndm();\r
-           if(randomNumber > fEfficiencyMatrix->GetBinContent(fEfficiencyMatrix->FindBin(gEtaAssoc,gPtAssoc,gPhiAssoc)))\r
-             continue;\r
-         }\r
-\r
-         gNumberOfAcceptedParticles += 1;\r
-\r
-         // add the track to the TObjArray\r
-         tracksMain->Add(new AliBFBasicParticle(gEtaAssoc, gPhiAssoc, gPtAssoc, gChargeAssoc, 1.0));  \r
-         if(fRunMixing) \r
-           tracksMixing->Add(new AliBFBasicParticle(gEtaAssoc, gPhiAssoc, gPtAssoc, gChargeAssoc, 1.0));  \r
-       }//pt,assoc < pt,trig\r
-      }//associated\r
-      \r
-      //back2back\r
-      gPtTrig2 = gPtTrig1;\r
-      gEtaTrig2 = -gEtaTrig1;\r
-      gPhiTrig2 = TMath::Pi() + gPhiTrig1;\r
-      if(gPhiTrig2 < 0.) gPhiTrig2 += TMath::TwoPi();\r
-      else if(gPhiTrig2 > TMath::TwoPi()) gPhiTrig2 -= TMath::TwoPi();\r
-\r
-      //Decide the charge\r
-      gDecideCharge = gRandom->Rndm();\r
-      if(gDecideCharge <= 0.5)\r
-       gChargeTrig2 = 1;\r
-      else \r
-       gChargeTrig2 = -1;\r
-      \r
-      //Acceptance\r
-      if((gEtaTrig2 < fEtaMin) || (gEtaTrig2 > fEtaMax)) continue;\r
-      //pt coverage\r
-      if((gPtTrig2 < fPtMin) || (gPtTrig2 > fPtMax)) continue;\r
-      //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);\r
-\r
-      //acceptance filter\r
-      if(fUseAcceptanceParameterization) {\r
-       Double_t gRandomNumberForAcceptance = gRandom->Rndm();\r
-       if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(gPtTrig2)) \r
-         continue;\r
-      }\r
-\r
-      //Detector effects\r
-      if(fSimulateDetectorEffects) {\r
-       Double_t randomNumber = gRandom->Rndm();\r
-       if(randomNumber > fEfficiencyMatrix->GetBinContent(fEfficiencyMatrix->FindBin(gEtaTrig2,gPtTrig2,gPhiTrig2)))\r
-         continue;\r
-      }\r
-\r
-      gNumberOfAcceptedParticles += 1;\r
-\r
-      // add the track to the TObjArray\r
-      tracksMain->Add(new AliBFBasicParticle(gEtaTrig2, gPhiTrig2, gPtTrig2, gChargeTrig2, 1.0));  \r
-      if(fRunMixing) \r
-       tracksMixing->Add(new AliBFBasicParticle(gEtaTrig2, gPhiTrig2, gPtTrig2, gChargeTrig2, 1.0));  \r
-\r
-      iAssociated = 0; \r
-      while(iAssociated < nAssociated) {\r
-       gPtAssoc = fPtAssoc->GetRandom();\r
-       if(gPtAssoc < gPtTrig2) {\r
-         gEtaAssoc = gRandom->Uniform(gEtaTrig2 - gJetCone/2.,gEtaTrig2 + gJetCone/2.);\r
-         gPhiAssoc = gRandom->Uniform(gPhiTrig2 - gJetCone/2.,gPhiTrig2 + gJetCone/2.);\r
-         if(gPhiAssoc < 0.) gPhiAssoc += TMath::TwoPi();\r
-         else if(gPhiAssoc > TMath::TwoPi()) gPhiAssoc -= TMath::TwoPi();\r
-         \r
-         iAssociated += 1;\r
-\r
-         gDecideCharge = gRandom->Rndm();\r
-         if(gDecideCharge <= 0.5)\r
-           gChargeAssoc = 1;\r
-         else \r
-           gChargeAssoc = -1;\r
-         \r
-         //Acceptance\r
-         if((gEtaAssoc < fEtaMin) || (gEtaAssoc > fEtaMax)) continue;\r
-         //pt coverage\r
-         if((gPtAssoc < fPtMin) || (gPtAssoc > fPtMax)) continue;\r
-         //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);\r
-\r
-         //acceptance filter\r
-         if(fUseAcceptanceParameterization) {\r
-           Double_t gRandomNumberForAcceptance = gRandom->Rndm();\r
-           if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(gPtAssoc)) \r
-             continue;\r
-         }\r
-\r
-         //Detector effects\r
-         if(fSimulateDetectorEffects) {\r
-           Double_t randomNumber = gRandom->Rndm();\r
-           if(randomNumber > fEfficiencyMatrix->GetBinContent(fEfficiencyMatrix->FindBin(gEtaAssoc,gPtAssoc,gPhiAssoc)))\r
-             continue;\r
-         }\r
-\r
-         gNumberOfAcceptedParticles += 1;\r
-\r
-         // add the track to the TObjArray\r
-         tracksMain->Add(new AliBFBasicParticle(gEtaAssoc, gPhiAssoc, gPtAssoc, gChargeAssoc, 1.0));  \r
-         if(fRunMixing) \r
-           tracksMixing->Add(new AliBFBasicParticle(gEtaAssoc, gPhiAssoc, gPtAssoc, gChargeAssoc, 1.0));  \r
-       }//pt,assoc < pt,trig\r
-      }//associated\r
-    }//Jet usage\r
-\r
-\r
-\r
-\r
-\r
-\r
-\r
-\r
-    \r
-    //Dynamical correlations\r
-    Int_t nGeneratedPositiveDynamicalCorrelations = 0;\r
-    Int_t nGeneratedNegativeDynamicalCorrelations = 0;\r
-    /*Double_t vChargePrime = 0;\r
-    Double_t vYPrime = 0.0;\r
-    Double_t vEtaPrime = 0.0;\r
-    Double_t vPhiPrime = 0.0;\r
-    Double_t vPPrime[3] = {0.,0.,0.};\r
-    Double_t vPtPrime = 0.0;\r
-    Double_t vEPrime = 0.0;\r
-    //Generate "correlated" particles \r
-    if(fUseDynamicalCorrelations) {\r
-      Int_t gNumberOfDynamicalCorrelations = (Int_t)(0.5*gNumberOfAcceptedParticles*fDynamicalCorrelationsPercentage);\r
-      for(Int_t iDynamicalCorrelations = 0; iDynamicalCorrelations < gNumberOfDynamicalCorrelations; iDynamicalCorrelations++) {\r
-       isPion = kFALSE; isKaon = kFALSE; isProton = kFALSE;\r
-       \r
-       //Pseudo-rapidity sampled from a Gaussian centered @ 0\r
-       vEta = gRandom->Gaus(0.0,0.1);\r
-       vCharge = 1.0;\r
-       nGeneratedPositiveDynamicalCorrelations += 1;\r
-       \r
-       vEtaPrime = -vEta;\r
-       vChargePrime = -1.0;\r
-       nGeneratedNegativeDynamicalCorrelations += 1;\r
-         \r
-       //Acceptance\r
-       if((vEta < fEtaMin) || (vEta > fEtaMax)) continue;\r
-       if((vEtaPrime < fEtaMin) || (vEtaPrime > fEtaMax)) continue;\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
-           vPt = fPtSpectraPions->GetRandom();\r
-           vPhi = fAzimuthalAnglePions->GetRandom();\r
-           fParticleMass = fPionMass;\r
-           isPion = kTRUE;\r
-         }\r
-         else if((randomNumberSpecies >= fPionPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage)) {\r
-           nGeneratedKaons += 1;\r
-           vPt = fPtSpectraKaons->GetRandom();\r
-           vPhi = fAzimuthalAngleKaons->GetRandom();\r
-           fParticleMass = fKaonMass;\r
-           isKaon = kTRUE;\r
-         }\r
-         else if((randomNumberSpecies >= fPionPercentage + fKaonPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage + fProtonPercentage)) {\r
-           nGeneratedProtons += 1;\r
-           vPt = fPtSpectraProtons->GetRandom();\r
-           vPtPrime = vPt;\r
-           vPhi = fAzimuthalAngleProtons->GetRandom();\r
-           fParticleMass = fProtonMass;\r
-           isProton = kTRUE;\r
-         }\r
-       }\r
-       else {\r
-         vPt = fPtSpectraAllCharges->GetRandom();\r
-         vPhi = fAzimuthalAngleAllCharges->GetRandom();\r
-       }\r
-       vPtPrime = vPt;\r
-       vPhiPrime = vPhi;\r
-\r
-       vP[0] = vPt*TMath::Cos(vPhi);\r
-       vP[1] = vPt*TMath::Sin(vPhi);\r
-       vP[2] = vPt*TMath::SinH(vEta);\r
-       vE = TMath::Sqrt(TMath::Power(fParticleMass,2) +\r
-                         TMath::Power(vP[0],2) +\r
-                         TMath::Power(vP[1],2) +\r
-                         TMath::Power(vP[2],2));\r
-       \r
-       vY = 0.5*TMath::Log((vE + vP[2])/(vE - vP[2]));\r
-\r
-       vPPrime[0] = vPtPrime*TMath::Cos(vPhiPrime);\r
-       vPPrime[1] = vPtPrime*TMath::Sin(vPhiPrime);\r
-       vPPrime[2] = vPtPrime*TMath::SinH(vEtaPrime);\r
-       vEPrime = TMath::Sqrt(TMath::Power(fParticleMass,2) +\r
-                         TMath::Power(vPPrime[0],2) +\r
-                         TMath::Power(vPPrime[1],2) +\r
-                         TMath::Power(vPPrime[2],2));\r
-       \r
-       vYPrime = 0.5*TMath::Log((vEPrime + vPPrime[2])/(vEPrime - vPPrime[2]));\r
-      \r
-       //pt coverage\r
-       if((vPt < fPtMin) || (vPt > fPtMax)) continue;\r
-       if((vPtPrime < fPtMin) || (vPtPrime > fPtMax)) continue;\r
-\r
-       //acceptance filter\r
-       if(fUseAcceptanceParameterization) {\r
-         Double_t gRandomNumberForAcceptance = gRandom->Rndm();\r
-         if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(vPt)) \r
-           continue;\r
-         \r
-         Double_t gRandomNumberForAcceptancePrime = gRandom->Rndm();\r
-         if(gRandomNumberForAcceptancePrime > fAcceptanceParameterization->Eval(vPtPrime)) \r
-           continue;\r
-       }\r
-      \r
-       // fill charge vector (positive)\r
-       chargeVector[0]->push_back(vCharge);\r
-       chargeVector[1]->push_back(vY);\r
-       chargeVector[2]->push_back(vEta);\r
-       chargeVector[3]->push_back(vPhi);\r
-       chargeVector[4]->push_back(vP[0]);\r
-       chargeVector[5]->push_back(vP[1]);\r
-       chargeVector[6]->push_back(vP[2]);\r
-       chargeVector[7]->push_back(vPt);\r
-       chargeVector[8]->push_back(vE);\r
-       \r
-       if(fRunShuffling) {\r
-         chargeVectorShuffle[0]->push_back(vCharge);\r
-         chargeVectorShuffle[1]->push_back(vY);\r
-         chargeVectorShuffle[2]->push_back(vEta);\r
-         chargeVectorShuffle[3]->push_back(vPhi);\r
-         chargeVectorShuffle[4]->push_back(vP[0]);\r
-         chargeVectorShuffle[5]->push_back(vP[1]);\r
-         chargeVectorShuffle[6]->push_back(vP[2]);\r
-         chargeVectorShuffle[7]->push_back(vPt);\r
-         chargeVectorShuffle[8]->push_back(vE);\r
-       }\r
-\r
-       // fill charge vector (negative)\r
-       chargeVector[0]->push_back(vChargePrime);\r
-       chargeVector[1]->push_back(vYPrime);\r
-       chargeVector[2]->push_back(vEtaPrime);\r
-       chargeVector[3]->push_back(vPhiPrime);\r
-       chargeVector[4]->push_back(vPPrime[0]);\r
-       chargeVector[5]->push_back(vPPrime[1]);\r
-       chargeVector[6]->push_back(vPPrime[2]);\r
-       chargeVector[7]->push_back(vPtPrime);\r
-       chargeVector[8]->push_back(vEPrime);\r
-       \r
-       if(fRunShuffling) {\r
-         chargeVectorShuffle[0]->push_back(vChargePrime);\r
-         chargeVectorShuffle[1]->push_back(vYPrime);\r
-         chargeVectorShuffle[2]->push_back(vEtaPrime);\r
-         chargeVectorShuffle[3]->push_back(vPhiPrime);\r
-         chargeVectorShuffle[4]->push_back(vPPrime[0]);\r
-         chargeVectorShuffle[5]->push_back(vPPrime[1]);\r
-         chargeVectorShuffle[6]->push_back(vPPrime[2]);\r
-         chargeVectorShuffle[7]->push_back(vPtPrime);\r
-         chargeVectorShuffle[8]->push_back(vEPrime);\r
-       }\r
-\r
-       gNumberOfAcceptedParticles += 2;\r
-       }//loop over the dynamical correlations\r
-       }*/// usage of dynamical correlations\r
-\r
-    if(fUseDebug) {\r
-      Printf("=======================================================");\r
-      Printf("Total: %d - Total positive: %d - Total negative: %d",nMultiplicity,nGeneratedPositive,nGeneratedNegative);\r
-      Printf("Accepted positive: %d - Accepted negative: %d",gNumberOfAcceptedPositiveParticles,gNumberOfAcceptedNegativeParticles);\r
-      Printf("Correlations: %d - Correlations positive: %d - Correlations negative: %d",nGeneratedPositiveDynamicalCorrelations+nGeneratedNegativeDynamicalCorrelations,nGeneratedPositiveDynamicalCorrelations,nGeneratedNegativeDynamicalCorrelations);\r
-      Printf("Number of accepted particles: %d",gNumberOfAcceptedParticles);\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
-\r
-    fHistEventStats->Fill(4);\r
-    fHistNumberOfAcceptedParticles->Fill(gNumberOfAcceptedParticles);\r
-    fHistReactionPlane->Fill(fReactionPlane);\r
-\r
-    // Event mixing \r
-    if (fRunMixing)\r
-      {\r
-       // 1. First get an event pool corresponding in mult (cent) and\r
-       //    zvertex to the current event. Once initialized, the pool\r
-       //    should contain nMix (reduced) events. This routine does not\r
-       //    pre-scan the chain. The first several events of every chain\r
-       //    will be skipped until the needed pools are filled to the\r
-       //    specified depth. If the pool categories are not too rare, this\r
-       //    should not be a problem. If they are rare, you could lose`\r
-       //    statistics.\r
-       \r
-       // 2. Collect the whole pool's content of tracks into one TObjArray\r
-       //    (bgTracks), which is effectively a single background super-event.\r
-       \r
-       // 3. The reduced and bgTracks arrays must both be passed into\r
-       //    FillCorrelations(). Also nMix should be passed in, so a weight\r
-       //    of 1./nMix can be applied.\r
-       \r
-       AliEventPool* pool = fPoolMgr->GetEventPool(1., 0.,fReactionPlane);\r
-      \r
-       if (!pool){\r
-         AliFatal(Form("No pool found for centrality = %f, zVtx = %f, psi = %f", 1., 0.,fReactionPlane));\r
-       }\r
-       else{\r
-         \r
-         //pool->SetDebug(1);\r
-         \r
-         if (pool->IsReady() || pool->NTracksInPool() > 50000 / 10 || pool->GetCurrentNEvents() >= 5){ \r
-           \r
-           \r
-           Int_t nMix = pool->GetCurrentNEvents();\r
-           //cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << ", tracks in this event = "<<gNumberOfAcceptedParticles <<", event Plane = "<<fReactionPlane<<endl;\r
-           \r
-           //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(2);\r
-           //((TH2F*) fListOfHistos->FindObject("mixedDist"))->Fill(centrality, pool->NTracksInPool());\r
-           //if (pool->IsReady())\r
-           //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(3);\r
-           \r
-           // Fill mixed-event histos here  \r
-           for (Int_t jMix=0; jMix<nMix; jMix++) \r
-             {\r
-               TObjArray* tracksMixed = pool->GetEvent(jMix);\r
-               fMixedBalance->CalculateBalance(fReactionPlane,tracksMain,tracksMixed,1,1.,0.);\r
-             }\r
-         }\r
-         \r
-         // Update the Event pool\r
-         pool->UpdatePool(tracksMain);\r
-         //pool->PrintInfo();\r
-         \r
-       }//pool NULL check  \r
-      }//run mixing\r
-    \r
-    //Calculate the balance function\r
-    fBalance->CalculateBalance(fReactionPlane,tracksMain,NULL,1,1.,0.);\r
-    //if(fRunMixing)\r
-    //  fMixedBalance->CalculateBalance(fReactionPlane,tracksMixing,NULL,1,1.,0.);          \r
-\r
-  }//event loop\r
-}      \r
-\r
-//________________________________________________________________________\r
-void  AliAnalysisTaskToyModel::FinishOutput() {\r
-  //Printf("END BF");\r
-\r
-  TFile *gOutput = TFile::Open("outputToyModel.root","recreate");\r
-  TDirectoryFile *dir = new TDirectoryFile("PWGCFEbyE.outputBalanceFunctionPsiAnalysis","PWGCFEbyE.outputBalanceFunctionPsiAnalysis");\r
-  \r
-  fList->SetName("listQA");\r
-  fList->SetOwner(kTRUE);\r
-  dir->Add(fList);\r
-  //fList->Write("listQA",TObject::kSingleKey);\r
-\r
-  if (!fBalance) {\r
-    AliError("ERROR: fBalance not available");\r
-    return;\r
-  }  \r
-  fListBF->SetName("listBF");\r
-  fListBF->SetOwner(kTRUE);\r
-  dir->Add(fListBF);\r
-  //fListBF->Write("listBF",TObject::kSingleKey);\r
-\r
-  if(fRunShuffling) {\r
-    if (!fShuffledBalance) {\r
-      AliError("ERROR: fShuffledBalance not available");\r
-      return;\r
-    }\r
-    fListBFS->SetName("listBFShuffled");\r
-    fListBFS->SetOwner(kTRUE);\r
-    dir->Add(fListBFS);\r
-    //fListBFS->Write("listBFShuffled",TObject::kSingleKey);\r
-  }\r
-\r
-  if(fRunMixing) {\r
-    if (!fMixedBalance) {\r
-      AliError("ERROR: fMixedBalance not available");\r
-      return;\r
-    }\r
-    fListBFM->SetName("listBFMixed");\r
-    fListBFM->SetOwner(kTRUE);\r
-    dir->Add(fListBFM);\r
-    //fListBFM->Write("listBFMixed",TObject::kSingleKey);\r
-  }\r
-\r
-  dir->Write(dir->GetName(),TObject::kSingleKey);\r
-  gOutput->Close();\r
-}\r
-\r
+#include "TChain.h"
+#include "TList.h"
+#include "TCanvas.h"
+#include "TParticle.h"
+#include "TLorentzVector.h"
+#include "TGraphErrors.h"
+#include "TH1.h"
+#include "TH2.h"
+#include "TH3.h"
+#include "TArrayF.h"
+#include "TF1.h"
+#include "TRandom.h"
+#include "TFile.h"
+
+#include "AliAnalysisManager.h"
+#include "AliLog.h"
+
+#include "AliEventPoolManager.h"           
+
+#include "AliAnalysisTaskToyModel.h"
+#include "AliBalance.h"
+#include "AliBalancePsi.h"
+#include "AliAnalysisTaskTriggeredBF.h"
+
+
+// Analysis task for the toy model analysis
+// Authors: Panos.Christakoglou@nikhef.nl
+
+using std::cout;
+using std::endl;
+
+ClassImp(AliAnalysisTaskToyModel)
+
+//________________________________________________________________________
+AliAnalysisTaskToyModel::AliAnalysisTaskToyModel() 
+: TObject(),
+  fUseDebug(kFALSE),
+  fBalance(0),
+  fRunShuffling(kFALSE), fShuffledBalance(0),
+  fRunMixing(kFALSE), fMixedBalance(0), fPoolMgr(0),
+  fList(0), fListBF(0), fListBFS(0), fListBFM(0),
+  fHistEventStats(0),
+  fHistNumberOfAcceptedParticles(0),
+  fHistReactionPlane(0),
+  fHistEtaTotal(0), fHistEta(0),
+  fHistEtaPhiPos(0), fHistEtaPhiNeg(0),
+  fHistRapidity(0),
+  fHistRapidityPions(0), fHistRapidityKaons(0), fHistRapidityProtons(0),
+  fHistPhi(0),
+  fHistPhiPions(0), fHistPhiKaons(0), fHistPhiProtons(0),
+  fHistPt(0),
+  fHistPtPions(0), fHistPtKaons(0), fHistPtProtons(0),
+  fTotalMultiplicityMean(100.), fTotalMultiplicitySigma(0.0),
+  fNetChargeMean(0.0), fNetChargeSigma(0.0),
+  fPtMin(0.0), fPtMax(100.0),
+  fEtaMin(-1.0), fEtaMax(1.0),
+  fUseAcceptanceParameterization(kFALSE), fAcceptanceParameterization(0),
+  fSimulateDetectorEffects(kFALSE),
+  fNumberOfInefficientSectors(0),
+  fInefficiencyFactorInPhi(1.0),
+  fNumberOfDeadSectors(0),
+  fEfficiencyDropNearEtaEdges(kFALSE),
+  fEfficiencyMatrix(0),
+  fUseAllCharges(kFALSE), fParticleMass(0.0),
+  fPtSpectraAllCharges(0), fTemperatureAllCharges(100.),
+  fReactionPlane(0.0),
+  fAzimuthalAngleAllCharges(0), fDirectedFlowAllCharges(0.0), 
+  fEllipticFlowAllCharges(0.0), fTriangularFlowAllCharges(0.0),
+  fQuandrangularFlowAllCharges(0.0), fPentangularFlowAllCharges(0.0),
+  fPionPercentage(0.8), fPionMass(0.0),
+  fPtSpectraPions(0), fTemperaturePions(100.),
+  fAzimuthalAnglePions(0), fDirectedFlowPions(0.0), 
+  fEllipticFlowPions(0.0), fTriangularFlowPions(0.0), 
+  fQuandrangularFlowPions(0.0), fPentangularFlowPions(0.0),
+  fKaonPercentage(0.8), fKaonMass(0.0),
+  fPtSpectraKaons(0), fTemperatureKaons(100.),
+  fAzimuthalAngleKaons(0), fDirectedFlowKaons(0.0), 
+  fEllipticFlowKaons(0.0), fTriangularFlowKaons(0.0),
+  fQuandrangularFlowKaons(0.0), fPentangularFlowKaons(0.0),
+  fProtonPercentage(0.8), fProtonMass(0.0),
+  fPtSpectraProtons(0), fTemperatureProtons(100.),
+  fAzimuthalAngleProtons(0), fDirectedFlowProtons(0.0), 
+  fEllipticFlowProtons(0.0), fTriangularFlowProtons(0.0),
+  fQuandrangularFlowProtons(0.0), fPentangularFlowProtons(0.0),
+  fUseDynamicalCorrelations(kFALSE), fDynamicalCorrelationsPercentage(0.1),
+  fUseJets(kFALSE), fPtAssoc(0) {
+  // Constructor
+}
+
+//________________________________________________________________________
+AliAnalysisTaskToyModel::~AliAnalysisTaskToyModel() {
+  //Destructor
+  if(fUseAllCharges) {
+    delete fPtSpectraAllCharges;
+    delete fAzimuthalAngleAllCharges;
+  }
+  else {
+    delete fPtSpectraPions;
+    delete fAzimuthalAnglePions;
+    delete fPtSpectraKaons;
+    delete fAzimuthalAngleKaons;
+    delete fPtSpectraProtons;
+    delete fAzimuthalAngleProtons;
+  }
+  if(fUseJets) delete fPtAssoc;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskToyModel::Init() {
+  //Initialize objects
+  //==============gRandom Seed=======================//
+  gRandom->SetSeed(0);  //seed is set to a random value (depending on machine clock)
+  //==============gRandom Seed=======================//
+
+  //==============Particles and spectra==============//
+  TParticle *pion = new TParticle();
+  pion->SetPdgCode(211);
+  fPionMass = pion->GetMass();
+
+  TParticle *kaon = new TParticle();
+  kaon->SetPdgCode(321);
+  fKaonMass = kaon->GetMass();
+  
+  TParticle *proton = new TParticle();
+  proton->SetPdgCode(2212);
+  fProtonMass = proton->GetMass();
+
+  if(fUseAllCharges) {
+    fParticleMass = fPionMass;
+    fPtSpectraAllCharges = new TF1("fPtSpectraAllCharges","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);
+    fPtSpectraAllCharges->SetParName(0,"Mass");
+    fPtSpectraAllCharges->SetParName(1,"Temperature");
+    //fPtSpectraAllCharges = new TF1("fPtSpectraAllCharges","(x^2/TMath::Sqrt(TMath::Power(x,2) + TMath::Power(0.139,2)))*TMath::Power((1. + x/[0]),-[1])",0.,20.);
+    //fPtSpectraAllCharges->SetParName(0,"pt0");
+    //fPtSpectraAllCharges->SetParName(1,"b");
+  }
+  else {
+    fPtSpectraPions = new TF1("fPtSpectraPions","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);
+    fPtSpectraPions->SetParName(0,"Mass");
+    fPtSpectraPions->SetParName(1,"Temperature");
+    
+    fPtSpectraKaons = new TF1("fPtSpectraKaons","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);
+    fPtSpectraKaons->SetParName(0,"Mass");
+    fPtSpectraKaons->SetParName(1,"Temperature");
+    
+    fPtSpectraProtons = new TF1("fPtSpectraProtons","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);
+    fPtSpectraProtons->SetParName(0,"Mass");
+    fPtSpectraProtons->SetParName(1,"Temperature");
+  }
+  //==============Particles and spectra==============//
+
+  //==============Flow values==============//
+  if(fUseAllCharges) {
+    if(fUseDebug) {
+      Printf("Directed flow: %lf",fDirectedFlowAllCharges);
+      Printf("Elliptic flow: %lf",fEllipticFlowAllCharges);
+      Printf("Triangular flow: %lf",fTriangularFlowAllCharges);
+      Printf("Quandrangular flow: %lf",fQuandrangularFlowAllCharges);
+      Printf("Pentangular flow: %lf",fPentangularFlowAllCharges);
+    }
+
+    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());
+    fAzimuthalAngleAllCharges->SetParName(0,"Reaction Plane");
+    fAzimuthalAngleAllCharges->SetParName(1,"Directed flow");
+    fAzimuthalAngleAllCharges->SetParName(2,"Elliptic flow"); 
+    fAzimuthalAngleAllCharges->SetParName(3,"Triangular flow");
+    fAzimuthalAngleAllCharges->SetParName(4,"Quandrangular flow");
+    fAzimuthalAngleAllCharges->SetParName(5,"Pentangular flow");
+  }
+  else {
+    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());
+    fAzimuthalAnglePions->SetParName(0,"Reaction Plane");
+    fAzimuthalAnglePions->SetParName(1,"Directed flow");
+    fAzimuthalAnglePions->SetParName(2,"Elliptic flow"); 
+    fAzimuthalAnglePions->SetParName(3,"Triangular flow");
+    fAzimuthalAnglePions->SetParName(4,"Quandrangular flow");
+    fAzimuthalAnglePions->SetParName(5,"Pentangular flow");
+    
+    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());
+    fAzimuthalAngleKaons->SetParName(0,"Reaction Plane");
+    fAzimuthalAngleKaons->SetParName(1,"Directed flow");
+    fAzimuthalAngleKaons->SetParName(2,"Elliptic flow"); 
+    fAzimuthalAngleKaons->SetParName(3,"Triangular flow");
+    fAzimuthalAngleKaons->SetParName(4,"Quandrangular flow");
+    fAzimuthalAngleKaons->SetParName(5,"Pentangular flow");
+    
+    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());
+    fAzimuthalAngleProtons->SetParName(0,"Reaction Plane");
+    fAzimuthalAngleProtons->SetParName(1,"Directed flow");
+    fAzimuthalAngleProtons->SetParName(2,"Elliptic flow"); 
+    fAzimuthalAngleProtons->SetParName(3,"Triangular flow");
+    fAzimuthalAngleProtons->SetParName(4,"Quandrangular flow");
+    fAzimuthalAngleProtons->SetParName(5,"Pentangular flow");
+  }
+  //==============Flow values==============//
+
+  //===================Jets===================//
+  if(fUseJets) {
+    fPtAssoc = new TF1("fPtAssoc","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,20.);
+    fPtAssoc->SetParName(0,"pt0");
+    fPtAssoc->SetParName(1,"b");
+    fPtAssoc->SetParameter(0,0.139);
+    fPtAssoc->SetParameter(1,0.5);
+    fPtAssoc->SetLineColor(1);
+  }
+
+  //===================Jets===================//
+
+  //==============Efficiency matrix==============//
+  if(fSimulateDetectorEffects) SetupEfficiencyMatrix();
+  //==============Efficiency matrix==============//
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskToyModel::SetupEfficiencyMatrix() {
+  //Setup the efficiency matrix
+  TH1F *hPt = new TH1F("hPt","",200,0.1,20.1);
+  TH1F *hEta = new TH1F("hEta","",20,-0.95,0.95);
+  TH1F *hPhi = new TH1F("hPhi","",72,0.,2.*TMath::Pi());
+  fEfficiencyMatrix = new TH3F("fEfficiencyMatrix","",
+                              hEta->GetNbinsX(),
+                              hEta->GetXaxis()->GetXmin(),
+                              hEta->GetXaxis()->GetXmax(),
+                              hPt->GetNbinsX(),
+                              hPt->GetXaxis()->GetXmin(),
+                              hPt->GetXaxis()->GetXmax(),
+                              hPhi->GetNbinsX(),
+                              hPhi->GetXaxis()->GetXmin(),
+                              hPhi->GetXaxis()->GetXmax());
+
+  //Efficiency in pt
+  Double_t epsilon[20] = {0.3,0.6,0.77,0.79,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80};
+  for(Int_t i=1;i<=20;i++) {
+    hPt->SetBinContent(i,epsilon[i-1]);
+    hPt->SetBinError(i,0.01);
+  }
+  for(Int_t i=21;i<=200;i++) {
+    hPt->SetBinContent(i,epsilon[19]);
+    hPt->SetBinError(i,0.01);
+  }
+
+  //Efficiency in eta
+  for(Int_t i=1;i<=hEta->GetNbinsX();i++) {
+    hEta->SetBinContent(i,1.0);
+    hEta->SetBinError(i,0.01);
+  }
+  if(fEfficiencyDropNearEtaEdges) {
+    hEta->SetBinContent(1,0.7); hEta->SetBinContent(2,0.8);
+    hEta->SetBinContent(3,0.9);
+    hEta->SetBinContent(18,0.9); hEta->SetBinContent(19,0.8);
+    hEta->SetBinContent(20,0.7);
+  }
+
+  //Efficiency in phi
+  for(Int_t i=1;i<=hPhi->GetNbinsX();i++) {
+    hPhi->SetBinContent(i,1.0);
+    hPhi->SetBinError(i,0.01);
+  }
+  for(Int_t i=1;i<=fNumberOfInefficientSectors;i++)
+    hPhi->SetBinContent(hPhi->FindBin(hPhi->GetRandom()),fInefficiencyFactorInPhi);
+  for(Int_t i=1;i<=fNumberOfDeadSectors;i++)
+    hPhi->SetBinContent(hPhi->FindBin(hPhi->GetRandom()),0.0);
+  
+  //Fill the 3D efficiency map
+  for(Int_t iBinX = 1; iBinX<=fEfficiencyMatrix->GetXaxis()->GetNbins();iBinX++) {
+    //cout<<"==================================="<<endl;
+    for(Int_t iBinY = 1; iBinY<=fEfficiencyMatrix->GetYaxis()->GetNbins();iBinY++) {
+      //cout<<"==================================="<<endl;
+      for(Int_t iBinZ = 1; iBinZ<=fEfficiencyMatrix->GetZaxis()->GetNbins();iBinZ++) {
+       fEfficiencyMatrix->SetBinContent(iBinX,iBinY,iBinZ,hEta->GetBinContent(iBinX)*hPt->GetBinContent(iBinY)*hPhi->GetBinContent(iBinZ));
+       //cout<<"Eta: "<<hEta->GetBinCenter(iBinX)<<" - Pt: "<<hPt->GetBinCenter(iBinY)<<" - Phi: "<<hPhi->GetBinCenter(iBinZ)<<" - "<<hEta->GetBinContent(iBinX)<<" , "<<hPt->GetBinContent(iBinY)<<" , "<<hPhi->GetBinContent(iBinZ)<<" - Efficiency: "<<hEta->GetBinContent(iBinX)*hPt->GetBinContent(iBinY)*hPhi->GetBinContent(iBinZ)<<endl;
+      }
+    }
+  }
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskToyModel::CreateOutputObjects() {
+  // Create histograms
+  // Called once
+
+  // global switch disabling the reference 
+  // (to avoid "Replacing existing TH1" if several wagons are created in train)
+  Bool_t oldStatus = TH1::AddDirectoryStatus();
+  TH1::AddDirectory(kFALSE);
+
+  if(!fBalance) {
+    fBalance = new AliBalancePsi();
+    fBalance->SetDeltaEtaMax(2.0);
+  }
+  if(fRunShuffling) {
+    if(!fShuffledBalance) {
+      fShuffledBalance = new AliBalancePsi();
+      fShuffledBalance->SetDeltaEtaMax(2.0);
+    }
+  }
+  if(fRunMixing) {
+    if(!fMixedBalance) {
+      fMixedBalance = new AliBalancePsi();
+      fMixedBalance->SetDeltaEtaMax(2.0);
+    }
+  }
+
+  //QA list
+  fList = new TList();
+  fList->SetName("listQA");
+  fList->SetOwner();
+
+  //Balance Function list
+  fListBF = new TList();
+  fListBF->SetName("listBF");
+  fListBF->SetOwner();
+
+  if(fRunShuffling) {
+    fListBFS = new TList();
+    fListBFS->SetName("listBFShuffled");
+    fListBFS->SetOwner();
+  }
+
+  if(fRunMixing) {
+    fListBFM = new TList();
+    fListBFM->SetName("listBFMixed");
+    fListBFM->SetOwner();
+  }
+
+  //==============QA================//
+  //Event stats.
+  TString gCutName[4] = {"Total","Offline trigger",
+                         "Vertex","Analyzed"};
+  fHistEventStats = new TH1F("fHistEventStats",
+                             "Event statistics;;N_{events}",
+                             4,0.5,4.5);
+  for(Int_t i = 1; i <= 4; i++)
+    fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
+  fList->Add(fHistEventStats);
+
+  fHistNumberOfAcceptedParticles = new TH1F("fHistNumberOfAcceptedParticles",";N_{acc.};Entries",10000,0,10000);
+  fList->Add(fHistNumberOfAcceptedParticles);
+  
+  fHistReactionPlane = new TH1F("fHistReactionPlane","Reaction plane angle;#Psi [rad];Entries",1000,0.,2.*TMath::Pi());
+  fList->Add(fHistReactionPlane);
+
+  //Pseudo-rapidity
+  fHistEtaTotal = new TH1F("fHistEtaTotal","Pseudo-rapidity (full phase space);#eta;Entries",1000,-15.,15.); 
+  fList->Add(fHistEtaTotal);
+
+  fHistEta = new TH1F("fHistEta","Pseudo-rapidity (acceptance);#eta;Entries",1000,-1.5,1.5); 
+  fList->Add(fHistEta);
+
+  fHistEtaPhiPos = new TH2F("fHistEtaPhiPos","#eta-#phi distribution (+);#eta;#varphi (rad)",80,-2.,2.,72,-TMath::Pi()/2.,3.*TMath::Pi()/2.);
+  fList->Add(fHistEtaPhiPos);
+  fHistEtaPhiNeg = new TH2F("fHistEtaPhiNeg","#eta-#phi distribution (-);#eta;#varphi (rad)",80,-2.,2.,72,-TMath::Pi()/2.,3.*TMath::Pi()/2.);
+  fList->Add(fHistEtaPhiNeg);
+
+  //Rapidity
+  fHistRapidity = new TH1F("fHistRapidity","Rapidity (acceptance);y;Entries",1000,-1.5,1.5); 
+  fList->Add(fHistRapidity);
+  fHistRapidityPions = new TH1F("fHistRapidityPions","Rapidity (acceptance - pions);y;Entries",1000,-1.5,1.5); 
+  fList->Add(fHistRapidityPions);
+  fHistRapidityKaons = new TH1F("fHistRapidityKaons","Rapidity (acceptance - kaons);y;Entries",1000,-1.5,1.5); 
+  fList->Add(fHistRapidityKaons);
+  fHistRapidityProtons = new TH1F("fHistRapidityProtons","Rapidity (acceptance - protons);y;Entries",1000,-1.5,1.5); 
+  fList->Add(fHistRapidityProtons);
+
+  //Phi
+  fHistPhi = new TH1F("fHistPhi","Phi (acceptance);#phi (rad);Entries",1000,0.,2*TMath::Pi());
+  fList->Add(fHistPhi);
+
+  fHistPhiPions = new TH1F("fHistPhiPions","Phi (acceptance - pions);#phi (rad);Entries",1000,0.,2*TMath::Pi());
+  fList->Add(fHistPhiPions);
+  fHistPhiKaons = new TH1F("fHistPhiKaons","Phi (acceptance - kaons);#phi (rad);Entries",1000,0.,2*TMath::Pi());
+  fList->Add(fHistPhiKaons);
+  fHistPhiProtons = new TH1F("fHistPhiProtons","Phi (acceptance - protons);#phi (rad);Entries",1000,0.,2*TMath::Pi());
+  fList->Add(fHistPhiProtons);
+
+  //Pt
+  fHistPt = new TH1F("fHistPt","Pt (acceptance);p_{t} (GeV/c);Entries",1000,0.,10.);
+  fList->Add(fHistPt);
+
+  fHistPtPions = new TH1F("fHistPtPions","Pt (acceptance - pions);p_{t} (GeV/c);Entries",1000,0.,10.);
+  fList->Add(fHistPtPions);
+  fHistPtKaons = new TH1F("fHistPtKaons","Pt (acceptance - kaons);p_{t} (GeV/c);Entries",1000,0.,10.);
+  fList->Add(fHistPtKaons);
+  fHistPtProtons = new TH1F("fHistPtProtons","Pt (acceptance - protons);p_{t} (GeV/c);Entries",1000,0.,10.);
+  fList->Add(fHistPtProtons);
+
+  if(fEfficiencyMatrix) fList->Add(fEfficiencyMatrix);
+
+  //==============Balance function histograms================//
+  // Initialize histograms if not done yet
+  if(!fBalance->GetHistNp()){
+    AliWarning("Histograms not yet initialized! --> Will be done now");
+    AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
+    fBalance->InitHistograms();
+  }
+
+  if(fRunShuffling) {
+    if(!fShuffledBalance->GetHistNp()) {
+      AliWarning("Histograms (shuffling) not yet initialized! --> Will be done now");
+      AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
+      fShuffledBalance->InitHistograms();
+    }
+  }
+
+  fListBF->Add(fBalance->GetHistNp());
+  fListBF->Add(fBalance->GetHistNn());
+  fListBF->Add(fBalance->GetHistNpn());
+  fListBF->Add(fBalance->GetHistNnn());
+  fListBF->Add(fBalance->GetHistNpp());
+  fListBF->Add(fBalance->GetHistNnp());
+
+  if(fRunShuffling) {
+    fListBFS->Add(fShuffledBalance->GetHistNp());
+    fListBFS->Add(fShuffledBalance->GetHistNn());
+    fListBFS->Add(fShuffledBalance->GetHistNpn());
+    fListBFS->Add(fShuffledBalance->GetHistNnn());
+    fListBFS->Add(fShuffledBalance->GetHistNpp());
+    fListBFS->Add(fShuffledBalance->GetHistNnp());
+  }
+
+  if(fRunMixing) {
+    fListBFM->Add(fMixedBalance->GetHistNp());
+    fListBFM->Add(fMixedBalance->GetHistNn());
+    fListBFM->Add(fMixedBalance->GetHistNpn());
+    fListBFM->Add(fMixedBalance->GetHistNnn());
+    fListBFM->Add(fMixedBalance->GetHistNpp());
+    fListBFM->Add(fMixedBalance->GetHistNnp());
+  }
+
+  // Event Mixing
+  if(fRunMixing){
+    Int_t fMixingTracks = 2000;
+    Int_t trackDepth = fMixingTracks; 
+    Int_t poolsize   = 1000;  // Maximum number of events, ignored in the present implemented of AliEventPoolManager
+    
+    // centrality bins
+    // Double_t centralityBins[] = {0.,1.,2.,3.,4.,5.,7.,10.,20.,30.,40.,50.,60.,70.,80.,100.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!
+    // Double_t* centbins        = centralityBins;
+    // Int_t nCentralityBins     = sizeof(centralityBins) / sizeof(Double_t) - 1;
+
+    // multiplicity bins
+    Double_t multiplicityBins[] = {0,10,20,30,40,50,60,70,80,100,100000}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!
+    Double_t* multbins        = multiplicityBins;
+    Int_t nMultiplicityBins     = sizeof(multiplicityBins) / sizeof(Double_t) - 1;
+    
+    // Zvtx bins
+    Double_t vertexBins[] = {-10., -7., -5., -3., -1., 1., 3., 5., 7., 10.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!
+    Double_t* vtxbins     = vertexBins;
+    Int_t nVertexBins     = sizeof(vertexBins) / sizeof(Double_t) - 1;
+    
+    // Event plane angle (Psi) bins
+    // Double_t psiBins[] = {0.,45.,135.,215.,305.,360.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!
+    // Double_t* psibins     = psiBins;
+    // Int_t nPsiBins     = sizeof(psiBins) / sizeof(Double_t) - 1;
+    
+    // // run the event mixing also in bins of event plane (statistics!)
+    // if(fRunMixingEventPlane){
+    //   if(fEventClass=="Multiplicity"){
+    //         fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins, nPsiBins, psibins);
+    //   }
+    //   else{
+    //         fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins, nPsiBins, psibins);
+    //   }
+    // }
+    // else{
+    //if(fEventClass=="Multiplicity"){
+    fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins);
+    //}
+    //else{
+    //fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);
+    //}
+  }
+
+  TH1::AddDirectory(oldStatus);
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskToyModel::Run(Int_t nEvents) {
+  // Main loop
+  // Called for each event
+  Short_t vCharge = 0;
+  Float_t vY = 0.0;
+  Float_t vEta = 0.0;
+  Float_t vPhi = 0.0;
+  Float_t vP[3] = {0.,0.,0.};
+  Float_t vPt = 0.0;
+  Float_t vE = 0.0;
+  Bool_t isPion = kFALSE, isKaon = kFALSE, isProton = kFALSE;
+
+  Double_t gDecideCharge = 0.;
+
+  if(fUseAllCharges) {
+    //fPtSpectraAllCharges->SetParameter(0,fParticleMass);
+    fPtSpectraAllCharges->SetParameter(0,1.05);
+    fPtSpectraAllCharges->SetParameter(1,fTemperatureAllCharges);
+
+    fAzimuthalAngleAllCharges->SetParameter(1,fDirectedFlowAllCharges);
+    fAzimuthalAngleAllCharges->SetParameter(2,fEllipticFlowAllCharges);
+    fAzimuthalAngleAllCharges->SetParameter(3,fTriangularFlowAllCharges);
+    fAzimuthalAngleAllCharges->SetParameter(4,fQuandrangularFlowAllCharges);
+    fAzimuthalAngleAllCharges->SetParameter(5,fPentangularFlowAllCharges);
+  }
+  else {
+    fPtSpectraPions->SetParameter(0,fPionMass);
+    fPtSpectraPions->SetParameter(1,fTemperaturePions);
+    fPtSpectraKaons->SetParameter(0,fKaonMass);
+    fPtSpectraKaons->SetParameter(1,fTemperatureKaons);
+    fPtSpectraProtons->SetParameter(0,fProtonMass);
+    fPtSpectraProtons->SetParameter(1,fTemperatureProtons);
+
+    fAzimuthalAnglePions->SetParameter(1,fDirectedFlowPions);
+    fAzimuthalAnglePions->SetParameter(2,fEllipticFlowPions);
+    fAzimuthalAnglePions->SetParameter(3,fTriangularFlowPions);
+    fAzimuthalAnglePions->SetParameter(4,fQuandrangularFlowPions);
+    fAzimuthalAnglePions->SetParameter(5,fPentangularFlowPions);
+
+    fAzimuthalAngleKaons->SetParameter(1,fDirectedFlowKaons);
+    fAzimuthalAngleKaons->SetParameter(2,fEllipticFlowKaons);
+    fAzimuthalAngleKaons->SetParameter(3,fTriangularFlowKaons);
+    fAzimuthalAngleKaons->SetParameter(4,fQuandrangularFlowKaons);
+    fAzimuthalAngleKaons->SetParameter(5,fPentangularFlowKaons);
+
+    fAzimuthalAngleProtons->SetParameter(1,fDirectedFlowProtons);
+    fAzimuthalAngleProtons->SetParameter(2,fEllipticFlowProtons);
+    fAzimuthalAngleProtons->SetParameter(3,fTriangularFlowProtons);
+    fAzimuthalAngleProtons->SetParameter(4,fQuandrangularFlowProtons);
+    fAzimuthalAngleProtons->SetParameter(5,fPentangularFlowProtons);
+  }
+
+
+  for(Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
+    // vector holding the charges/kinematics of all tracks (charge,y,eta,phi,p0,p1,p2,pt,E)
+    //vector<Double_t> *chargeVectorShuffle[9];   // this will be shuffled
+    //vector<Double_t> *chargeVector[9];          // original charge
+    //for(Int_t i = 0; i < 9; i++){
+    //chargeVectorShuffle[i] = new vector<Double_t>;
+    //chargeVector[i]        = new vector<Double_t>;
+    //}
+
+    // TObjArray for the accepted particles 
+    // (has to be done here, otherwise mxing with event pool does not work, overwriting pointers!)
+    TObjArray *tracksMain = new TObjArray();
+    tracksMain->SetOwner(kTRUE);
+    TObjArray *tracksMixing = 0x0;
+    if(fRunMixing) {
+      tracksMixing = new TObjArray();
+      tracksMixing->SetOwner(kTRUE);
+    }
+
+    tracksMain->Clear();
+    if(fRunMixing) tracksMixing->Clear();
+    
+    fHistEventStats->Fill(1);
+    fHistEventStats->Fill(2);
+    fHistEventStats->Fill(3);
+
+    if((iEvent%1000) == 0) 
+      cout<<"Event: "<<iEvent<<"/"<<nEvents<<endl;
+    
+    //Multiplicities
+    Int_t nMultiplicity = (Int_t)(gRandom->Gaus(fTotalMultiplicityMean,fTotalMultiplicitySigma));
+    Int_t nNetCharge = (Int_t)(gRandom->Gaus(fNetChargeMean,fNetChargeSigma));
+    
+    Int_t nGeneratedPositive = (Int_t)((nMultiplicity/2) + nNetCharge);
+    Int_t nGeneratedNegative = nMultiplicity - nGeneratedPositive;
+    if(fUseDebug) 
+      Printf("Total multiplicity: %d - Generated positive: %d - Generated negative: %d",nMultiplicity,nGeneratedPositive,nGeneratedNegative);
+    
+    //Randomization of the reaction plane
+    fReactionPlane = 2.0*TMath::Pi()*gRandom->Rndm();
+    //fReactionPlane = 0.0;
+    if(fUseAllCharges) 
+      fAzimuthalAngleAllCharges->SetParameter(0,fReactionPlane);
+    else {
+      fAzimuthalAnglePions->SetParameter(0,fReactionPlane);
+      fAzimuthalAngleKaons->SetParameter(0,fReactionPlane);
+      fAzimuthalAngleProtons->SetParameter(0,fReactionPlane);
+    }
+    
+    Int_t gNumberOfAcceptedParticles = 0;
+    Int_t gNumberOfAcceptedPositiveParticles = 0;
+    Int_t gNumberOfAcceptedNegativeParticles = 0;
+    Int_t nGeneratedPions = 0, nGeneratedKaons = 0, nGeneratedProtons = 0;
+    //Generate particles
+    for(Int_t iParticleCount = 0; iParticleCount < nMultiplicity; iParticleCount++) {
+      isPion = kFALSE; isKaon = kFALSE; isProton = kFALSE;
+      if(fUseDebug) 
+       Printf("Generating positive: %d(%d)",iParticleCount+1,nGeneratedPositive);
+
+      //Pseudo-rapidity sampled from a Gaussian centered @ 0
+      vEta = gRandom->Gaus(0.0,4.0);
+
+      //Fill QA histograms (full phase space)
+      fHistEtaTotal->Fill(vEta);
+
+      //Decide the charge
+      gDecideCharge = gRandom->Rndm();
+      if(gDecideCharge <= 1.*nGeneratedPositive/nMultiplicity)       
+       vCharge = 1;
+      else 
+       vCharge = -1;
+      
+      //Acceptance
+      if((vEta < fEtaMin) || (vEta > fEtaMax)) continue;
+
+      if(!fUseAllCharges) {
+       //Decide the specie
+       Double_t randomNumberSpecies = gRandom->Rndm();
+       if((randomNumberSpecies >= 0.0)&&(randomNumberSpecies < fPionPercentage)) {
+         nGeneratedPions += 1;
+         vPt = fPtSpectraPions->GetRandom();
+         vPhi = fAzimuthalAnglePions->GetRandom();
+         fParticleMass = fPionMass;
+         isPion = kTRUE;
+       }
+       else if((randomNumberSpecies >= fPionPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage)) {
+         nGeneratedKaons += 1;
+         vPt = fPtSpectraKaons->GetRandom();
+         vPhi = fAzimuthalAngleKaons->GetRandom();
+         fParticleMass = fKaonMass;
+         isKaon = kTRUE;
+       }
+       else if((randomNumberSpecies >= fPionPercentage + fKaonPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage + fProtonPercentage)) {
+         nGeneratedProtons += 1;
+         vPt = fPtSpectraProtons->GetRandom();
+         vPhi = fAzimuthalAngleProtons->GetRandom();
+         fParticleMass = fProtonMass;
+         isProton = kTRUE;
+       }
+      }
+      else {
+       vPt = fPtSpectraAllCharges->GetRandom();
+       vPhi = fAzimuthalAngleAllCharges->GetRandom();
+      }
+      
+      vP[0] = vPt*TMath::Cos(vPhi);
+      vP[1] = vPt*TMath::Sin(vPhi);
+      vP[2] = vPt*TMath::SinH(vEta);
+      vE = TMath::Sqrt(TMath::Power(fParticleMass,2) +
+                       TMath::Power(vP[0],2) +
+                       TMath::Power(vP[1],2) +
+                       TMath::Power(vP[2],2));
+      
+      vY = 0.5*TMath::Log((vE + vP[2])/(vE - vP[2]));
+      
+      //pt coverage
+      if((vPt < fPtMin) || (vPt > fPtMax)) continue;
+      //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);
+
+      //acceptance filter
+      if(fUseAcceptanceParameterization) {
+       Double_t gRandomNumberForAcceptance = gRandom->Rndm();
+       if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(vPt)) 
+         continue;
+      }
+
+      //Detector effects
+      if(fSimulateDetectorEffects) {
+       Double_t randomNumber = gRandom->Rndm();
+       if(randomNumber > fEfficiencyMatrix->GetBinContent(fEfficiencyMatrix->FindBin(vEta,vPt,vPhi)))
+         continue;
+      }
+
+      //Fill QA histograms (acceptance);
+      if(vCharge > 0) {
+       gNumberOfAcceptedPositiveParticles += 1;
+       if(vPhi > 3.*TMath::Pi()/2.)
+         fHistEtaPhiPos->Fill(vEta,vPhi-2.*TMath::Pi());
+       else
+         fHistEtaPhiPos->Fill(vEta,vPhi);
+      }
+      else {
+       gNumberOfAcceptedNegativeParticles += 1;
+       if(vPhi > 3.*TMath::Pi()/2.)
+         fHistEtaPhiNeg->Fill(vEta,vPhi-2.*TMath::Pi());
+       else
+         fHistEtaPhiNeg->Fill(vEta,vPhi);
+      }
+
+      fHistEta->Fill(vEta);
+      fHistRapidity->Fill(vY);
+      fHistPhi->Fill(vPhi);
+      fHistPt->Fill(vPt);
+      if(isPion) {
+       fHistRapidityPions->Fill(vY);
+       fHistPhiPions->Fill(vPhi);
+       fHistPtPions->Fill(vPt);
+      }
+      else if(isKaon) {
+       fHistRapidityKaons->Fill(vY);
+       fHistPhiKaons->Fill(vPhi);
+       fHistPtKaons->Fill(vPt);
+      }
+      else if(isProton) {
+       fHistRapidityProtons->Fill(vY);
+       fHistPhiProtons->Fill(vPhi);
+       fHistPtProtons->Fill(vPt);
+      }
+      gNumberOfAcceptedParticles += 1;
+
+      // add the track to the TObjArray
+      tracksMain->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, 1.0));  
+      if(fRunMixing) 
+       tracksMixing->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, 1.0));  
+    }//generated positive particle loop
+
+    //Jets
+    if(fUseJets) {
+      const Int_t nAssociated = 3;
+
+      Double_t gPtTrig1 = 0., gPtTrig2 = 0., gPtAssoc = 0.;
+      Double_t gPhiTrig1 = 0., gPhiTrig2 = 0., gPhiAssoc = 0.;
+      Double_t gEtaTrig1 = 0., gEtaTrig2 = 0., gEtaAssoc = 0.;
+      Short_t gChargeTrig1 = 0, gChargeTrig2 = 0, gChargeAssoc = 0;
+      Double_t gJetCone = 0.2;
+
+      //First leading particle
+      gPtTrig1 = gRandom->Uniform(3.,5.);
+      gEtaTrig1 = gRandom->Uniform(-0.8,0.8);
+      gPhiTrig1 = gRandom->Uniform(0.,TMath::TwoPi());
+
+      //Decide the charge
+      gDecideCharge = gRandom->Rndm();
+      if(gDecideCharge <= 0.5)
+       gChargeTrig1 = 1;
+      else 
+       gChargeTrig1 = -1;
+      
+      //Acceptance
+      if((gEtaTrig1 < fEtaMin) || (gEtaTrig1 > fEtaMax)) continue;
+      //pt coverage
+      if((gPtTrig1 < fPtMin) || (gPtTrig1 > fPtMax)) continue;
+      //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);
+
+      //acceptance filter
+      if(fUseAcceptanceParameterization) {
+       Double_t gRandomNumberForAcceptance = gRandom->Rndm();
+       if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(gPtTrig1)) 
+         continue;
+      }
+
+      //Detector effects
+      if(fSimulateDetectorEffects) {
+       Double_t randomNumber = gRandom->Rndm();
+       if(randomNumber > fEfficiencyMatrix->GetBinContent(fEfficiencyMatrix->FindBin(gEtaTrig1,gPtTrig1,gPhiTrig1)))
+         continue;
+      }
+
+      gNumberOfAcceptedParticles += 1;
+
+      // add the track to the TObjArray
+      tracksMain->Add(new AliBFBasicParticle(gEtaTrig1, gPhiTrig1, gPtTrig1, gChargeTrig1, 1.0));  
+      if(fRunMixing) 
+       tracksMixing->Add(new AliBFBasicParticle(gEtaTrig1, gPhiTrig1, gPtTrig1, gChargeTrig1, 1.0));  
+
+      Int_t iAssociated = 0; 
+      while(iAssociated < nAssociated) {
+       gPtAssoc = fPtAssoc->GetRandom();
+       if(gPtAssoc < gPtTrig1) {
+         gEtaAssoc = gRandom->Uniform(gEtaTrig1 - gJetCone/2.,gEtaTrig1 + gJetCone/2.);
+         gPhiAssoc = gRandom->Uniform(gPhiTrig1 - gJetCone/2.,gPhiTrig1 + gJetCone/2.);
+         if(gPhiAssoc < 0.) gPhiAssoc += TMath::TwoPi();
+         else if(gPhiAssoc > TMath::TwoPi()) gPhiAssoc -= TMath::TwoPi();
+         
+         iAssociated += 1;
+
+         gDecideCharge = gRandom->Rndm();
+         if(gDecideCharge <= 0.5)
+           gChargeAssoc = 1;
+         else 
+           gChargeAssoc = -1;
+         
+         //Acceptance
+         if((gEtaAssoc < fEtaMin) || (gEtaAssoc > fEtaMax)) continue;
+         //pt coverage
+         if((gPtAssoc < fPtMin) || (gPtAssoc > fPtMax)) continue;
+         //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);
+
+         //acceptance filter
+         if(fUseAcceptanceParameterization) {
+           Double_t gRandomNumberForAcceptance = gRandom->Rndm();
+           if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(gPtAssoc)) 
+             continue;
+         }
+
+         //Detector effects
+         if(fSimulateDetectorEffects) {
+           Double_t randomNumber = gRandom->Rndm();
+           if(randomNumber > fEfficiencyMatrix->GetBinContent(fEfficiencyMatrix->FindBin(gEtaAssoc,gPtAssoc,gPhiAssoc)))
+             continue;
+         }
+
+         gNumberOfAcceptedParticles += 1;
+
+         // add the track to the TObjArray
+         tracksMain->Add(new AliBFBasicParticle(gEtaAssoc, gPhiAssoc, gPtAssoc, gChargeAssoc, 1.0));  
+         if(fRunMixing) 
+           tracksMixing->Add(new AliBFBasicParticle(gEtaAssoc, gPhiAssoc, gPtAssoc, gChargeAssoc, 1.0));  
+       }//pt,assoc < pt,trig
+      }//associated
+      
+      //back2back
+      gPtTrig2 = gPtTrig1;
+      gEtaTrig2 = -gEtaTrig1;
+      gPhiTrig2 = TMath::Pi() + gPhiTrig1;
+      if(gPhiTrig2 < 0.) gPhiTrig2 += TMath::TwoPi();
+      else if(gPhiTrig2 > TMath::TwoPi()) gPhiTrig2 -= TMath::TwoPi();
+
+      //Decide the charge
+      gDecideCharge = gRandom->Rndm();
+      if(gDecideCharge <= 0.5)
+       gChargeTrig2 = 1;
+      else 
+       gChargeTrig2 = -1;
+      
+      //Acceptance
+      if((gEtaTrig2 < fEtaMin) || (gEtaTrig2 > fEtaMax)) continue;
+      //pt coverage
+      if((gPtTrig2 < fPtMin) || (gPtTrig2 > fPtMax)) continue;
+      //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);
+
+      //acceptance filter
+      if(fUseAcceptanceParameterization) {
+       Double_t gRandomNumberForAcceptance = gRandom->Rndm();
+       if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(gPtTrig2)) 
+         continue;
+      }
+
+      //Detector effects
+      if(fSimulateDetectorEffects) {
+       Double_t randomNumber = gRandom->Rndm();
+       if(randomNumber > fEfficiencyMatrix->GetBinContent(fEfficiencyMatrix->FindBin(gEtaTrig2,gPtTrig2,gPhiTrig2)))
+         continue;
+      }
+
+      gNumberOfAcceptedParticles += 1;
+
+      // add the track to the TObjArray
+      tracksMain->Add(new AliBFBasicParticle(gEtaTrig2, gPhiTrig2, gPtTrig2, gChargeTrig2, 1.0));  
+      if(fRunMixing) 
+       tracksMixing->Add(new AliBFBasicParticle(gEtaTrig2, gPhiTrig2, gPtTrig2, gChargeTrig2, 1.0));  
+
+      iAssociated = 0; 
+      while(iAssociated < nAssociated) {
+       gPtAssoc = fPtAssoc->GetRandom();
+       if(gPtAssoc < gPtTrig2) {
+         gEtaAssoc = gRandom->Uniform(gEtaTrig2 - gJetCone/2.,gEtaTrig2 + gJetCone/2.);
+         gPhiAssoc = gRandom->Uniform(gPhiTrig2 - gJetCone/2.,gPhiTrig2 + gJetCone/2.);
+         if(gPhiAssoc < 0.) gPhiAssoc += TMath::TwoPi();
+         else if(gPhiAssoc > TMath::TwoPi()) gPhiAssoc -= TMath::TwoPi();
+         
+         iAssociated += 1;
+
+         gDecideCharge = gRandom->Rndm();
+         if(gDecideCharge <= 0.5)
+           gChargeAssoc = 1;
+         else 
+           gChargeAssoc = -1;
+         
+         //Acceptance
+         if((gEtaAssoc < fEtaMin) || (gEtaAssoc > fEtaMax)) continue;
+         //pt coverage
+         if((gPtAssoc < fPtMin) || (gPtAssoc > fPtMax)) continue;
+         //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);
+
+         //acceptance filter
+         if(fUseAcceptanceParameterization) {
+           Double_t gRandomNumberForAcceptance = gRandom->Rndm();
+           if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(gPtAssoc)) 
+             continue;
+         }
+
+         //Detector effects
+         if(fSimulateDetectorEffects) {
+           Double_t randomNumber = gRandom->Rndm();
+           if(randomNumber > fEfficiencyMatrix->GetBinContent(fEfficiencyMatrix->FindBin(gEtaAssoc,gPtAssoc,gPhiAssoc)))
+             continue;
+         }
+
+         gNumberOfAcceptedParticles += 1;
+
+         // add the track to the TObjArray
+         tracksMain->Add(new AliBFBasicParticle(gEtaAssoc, gPhiAssoc, gPtAssoc, gChargeAssoc, 1.0));  
+         if(fRunMixing) 
+           tracksMixing->Add(new AliBFBasicParticle(gEtaAssoc, gPhiAssoc, gPtAssoc, gChargeAssoc, 1.0));  
+       }//pt,assoc < pt,trig
+      }//associated
+    }//Jet usage
+
+
+
+
+
+
+
+
+    
+    //Dynamical correlations
+    Int_t nGeneratedPositiveDynamicalCorrelations = 0;
+    Int_t nGeneratedNegativeDynamicalCorrelations = 0;
+    /*Double_t vChargePrime = 0;
+    Double_t vYPrime = 0.0;
+    Double_t vEtaPrime = 0.0;
+    Double_t vPhiPrime = 0.0;
+    Double_t vPPrime[3] = {0.,0.,0.};
+    Double_t vPtPrime = 0.0;
+    Double_t vEPrime = 0.0;
+    //Generate "correlated" particles 
+    if(fUseDynamicalCorrelations) {
+      Int_t gNumberOfDynamicalCorrelations = (Int_t)(0.5*gNumberOfAcceptedParticles*fDynamicalCorrelationsPercentage);
+      for(Int_t iDynamicalCorrelations = 0; iDynamicalCorrelations < gNumberOfDynamicalCorrelations; iDynamicalCorrelations++) {
+       isPion = kFALSE; isKaon = kFALSE; isProton = kFALSE;
+       
+       //Pseudo-rapidity sampled from a Gaussian centered @ 0
+       vEta = gRandom->Gaus(0.0,0.1);
+       vCharge = 1.0;
+       nGeneratedPositiveDynamicalCorrelations += 1;
+       
+       vEtaPrime = -vEta;
+       vChargePrime = -1.0;
+       nGeneratedNegativeDynamicalCorrelations += 1;
+         
+       //Acceptance
+       if((vEta < fEtaMin) || (vEta > fEtaMax)) continue;
+       if((vEtaPrime < fEtaMin) || (vEtaPrime > fEtaMax)) continue;
+      
+       if(!fUseAllCharges) {
+         //Decide the specie
+         Double_t randomNumberSpecies = gRandom->Rndm();
+         if((randomNumberSpecies >= 0.0)&&(randomNumberSpecies < fPionPercentage)) {
+           nGeneratedPions += 1;
+           vPt = fPtSpectraPions->GetRandom();
+           vPhi = fAzimuthalAnglePions->GetRandom();
+           fParticleMass = fPionMass;
+           isPion = kTRUE;
+         }
+         else if((randomNumberSpecies >= fPionPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage)) {
+           nGeneratedKaons += 1;
+           vPt = fPtSpectraKaons->GetRandom();
+           vPhi = fAzimuthalAngleKaons->GetRandom();
+           fParticleMass = fKaonMass;
+           isKaon = kTRUE;
+         }
+         else if((randomNumberSpecies >= fPionPercentage + fKaonPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage + fProtonPercentage)) {
+           nGeneratedProtons += 1;
+           vPt = fPtSpectraProtons->GetRandom();
+           vPtPrime = vPt;
+           vPhi = fAzimuthalAngleProtons->GetRandom();
+           fParticleMass = fProtonMass;
+           isProton = kTRUE;
+         }
+       }
+       else {
+         vPt = fPtSpectraAllCharges->GetRandom();
+         vPhi = fAzimuthalAngleAllCharges->GetRandom();
+       }
+       vPtPrime = vPt;
+       vPhiPrime = vPhi;
+
+       vP[0] = vPt*TMath::Cos(vPhi);
+       vP[1] = vPt*TMath::Sin(vPhi);
+       vP[2] = vPt*TMath::SinH(vEta);
+       vE = TMath::Sqrt(TMath::Power(fParticleMass,2) +
+                         TMath::Power(vP[0],2) +
+                         TMath::Power(vP[1],2) +
+                         TMath::Power(vP[2],2));
+       
+       vY = 0.5*TMath::Log((vE + vP[2])/(vE - vP[2]));
+
+       vPPrime[0] = vPtPrime*TMath::Cos(vPhiPrime);
+       vPPrime[1] = vPtPrime*TMath::Sin(vPhiPrime);
+       vPPrime[2] = vPtPrime*TMath::SinH(vEtaPrime);
+       vEPrime = TMath::Sqrt(TMath::Power(fParticleMass,2) +
+                         TMath::Power(vPPrime[0],2) +
+                         TMath::Power(vPPrime[1],2) +
+                         TMath::Power(vPPrime[2],2));
+       
+       vYPrime = 0.5*TMath::Log((vEPrime + vPPrime[2])/(vEPrime - vPPrime[2]));
+      
+       //pt coverage
+       if((vPt < fPtMin) || (vPt > fPtMax)) continue;
+       if((vPtPrime < fPtMin) || (vPtPrime > fPtMax)) continue;
+
+       //acceptance filter
+       if(fUseAcceptanceParameterization) {
+         Double_t gRandomNumberForAcceptance = gRandom->Rndm();
+         if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(vPt)) 
+           continue;
+         
+         Double_t gRandomNumberForAcceptancePrime = gRandom->Rndm();
+         if(gRandomNumberForAcceptancePrime > fAcceptanceParameterization->Eval(vPtPrime)) 
+           continue;
+       }
+      
+       // fill charge vector (positive)
+       chargeVector[0]->push_back(vCharge);
+       chargeVector[1]->push_back(vY);
+       chargeVector[2]->push_back(vEta);
+       chargeVector[3]->push_back(vPhi);
+       chargeVector[4]->push_back(vP[0]);
+       chargeVector[5]->push_back(vP[1]);
+       chargeVector[6]->push_back(vP[2]);
+       chargeVector[7]->push_back(vPt);
+       chargeVector[8]->push_back(vE);
+       
+       if(fRunShuffling) {
+         chargeVectorShuffle[0]->push_back(vCharge);
+         chargeVectorShuffle[1]->push_back(vY);
+         chargeVectorShuffle[2]->push_back(vEta);
+         chargeVectorShuffle[3]->push_back(vPhi);
+         chargeVectorShuffle[4]->push_back(vP[0]);
+         chargeVectorShuffle[5]->push_back(vP[1]);
+         chargeVectorShuffle[6]->push_back(vP[2]);
+         chargeVectorShuffle[7]->push_back(vPt);
+         chargeVectorShuffle[8]->push_back(vE);
+       }
+
+       // fill charge vector (negative)
+       chargeVector[0]->push_back(vChargePrime);
+       chargeVector[1]->push_back(vYPrime);
+       chargeVector[2]->push_back(vEtaPrime);
+       chargeVector[3]->push_back(vPhiPrime);
+       chargeVector[4]->push_back(vPPrime[0]);
+       chargeVector[5]->push_back(vPPrime[1]);
+       chargeVector[6]->push_back(vPPrime[2]);
+       chargeVector[7]->push_back(vPtPrime);
+       chargeVector[8]->push_back(vEPrime);
+       
+       if(fRunShuffling) {
+         chargeVectorShuffle[0]->push_back(vChargePrime);
+         chargeVectorShuffle[1]->push_back(vYPrime);
+         chargeVectorShuffle[2]->push_back(vEtaPrime);
+         chargeVectorShuffle[3]->push_back(vPhiPrime);
+         chargeVectorShuffle[4]->push_back(vPPrime[0]);
+         chargeVectorShuffle[5]->push_back(vPPrime[1]);
+         chargeVectorShuffle[6]->push_back(vPPrime[2]);
+         chargeVectorShuffle[7]->push_back(vPtPrime);
+         chargeVectorShuffle[8]->push_back(vEPrime);
+       }
+
+       gNumberOfAcceptedParticles += 2;
+       }//loop over the dynamical correlations
+       }*/// usage of dynamical correlations
+
+    if(fUseDebug) {
+      Printf("=======================================================");
+      Printf("Total: %d - Total positive: %d - Total negative: %d",nMultiplicity,nGeneratedPositive,nGeneratedNegative);
+      Printf("Accepted positive: %d - Accepted negative: %d",gNumberOfAcceptedPositiveParticles,gNumberOfAcceptedNegativeParticles);
+      Printf("Correlations: %d - Correlations positive: %d - Correlations negative: %d",nGeneratedPositiveDynamicalCorrelations+nGeneratedNegativeDynamicalCorrelations,nGeneratedPositiveDynamicalCorrelations,nGeneratedNegativeDynamicalCorrelations);
+      Printf("Number of accepted particles: %d",gNumberOfAcceptedParticles);
+      if(!fUseAllCharges)
+       Printf("Pions: %lf - Kaons: %lf - Protons: %lf",1.*nGeneratedPions/nMultiplicity,1.*nGeneratedKaons/nMultiplicity,1.*nGeneratedProtons/nMultiplicity);
+      //Printf("Calculating the balance function for %d particles",chargeVector[0]->size());
+    }
+
+    fHistEventStats->Fill(4);
+    fHistNumberOfAcceptedParticles->Fill(gNumberOfAcceptedParticles);
+    fHistReactionPlane->Fill(fReactionPlane);
+
+    // Event mixing 
+    if (fRunMixing)
+      {
+       // 1. First get an event pool corresponding in mult (cent) and
+       //    zvertex to the current event. Once initialized, the pool
+       //    should contain nMix (reduced) events. This routine does not
+       //    pre-scan the chain. The first several events of every chain
+       //    will be skipped until the needed pools are filled to the
+       //    specified depth. If the pool categories are not too rare, this
+       //    should not be a problem. If they are rare, you could lose`
+       //    statistics.
+       
+       // 2. Collect the whole pool's content of tracks into one TObjArray
+       //    (bgTracks), which is effectively a single background super-event.
+       
+       // 3. The reduced and bgTracks arrays must both be passed into
+       //    FillCorrelations(). Also nMix should be passed in, so a weight
+       //    of 1./nMix can be applied.
+       
+       AliEventPool* pool = fPoolMgr->GetEventPool(1., 0.,fReactionPlane);
+      
+       if (!pool){
+         AliFatal(Form("No pool found for centrality = %f, zVtx = %f, psi = %f", 1., 0.,fReactionPlane));
+       }
+       else{
+         
+         //pool->SetDebug(1);
+         
+         if (pool->IsReady() || pool->NTracksInPool() > 50000 / 10 || pool->GetCurrentNEvents() >= 5){ 
+           
+           
+           Int_t nMix = pool->GetCurrentNEvents();
+           //cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << ", tracks in this event = "<<gNumberOfAcceptedParticles <<", event Plane = "<<fReactionPlane<<endl;
+           
+           //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(2);
+           //((TH2F*) fListOfHistos->FindObject("mixedDist"))->Fill(centrality, pool->NTracksInPool());
+           //if (pool->IsReady())
+           //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(3);
+           
+           // Fill mixed-event histos here  
+           for (Int_t jMix=0; jMix<nMix; jMix++) 
+             {
+               TObjArray* tracksMixed = pool->GetEvent(jMix);
+               fMixedBalance->CalculateBalance(fReactionPlane,tracksMain,tracksMixed,1,1.,0.);
+             }
+         }
+         
+         // Update the Event pool
+         pool->UpdatePool(tracksMain);
+         //pool->PrintInfo();
+         
+       }//pool NULL check  
+      }//run mixing
+    
+    //Calculate the balance function
+    fBalance->CalculateBalance(fReactionPlane,tracksMain,NULL,1,1.,0.);
+    //if(fRunMixing)
+    //  fMixedBalance->CalculateBalance(fReactionPlane,tracksMixing,NULL,1,1.,0.);          
+
+  }//event loop
+}      
+
+//________________________________________________________________________
+void  AliAnalysisTaskToyModel::FinishOutput() {
+  //Printf("END BF");
+
+  TFile *gOutput = TFile::Open("outputToyModel.root","recreate");
+  TDirectoryFile *dir = new TDirectoryFile("PWGCFEbyE.outputBalanceFunctionPsiAnalysis","PWGCFEbyE.outputBalanceFunctionPsiAnalysis");
+  
+  fList->SetName("listQA");
+  fList->SetOwner(kTRUE);
+  dir->Add(fList);
+  //fList->Write("listQA",TObject::kSingleKey);
+
+  if (!fBalance) {
+    AliError("ERROR: fBalance not available");
+    return;
+  }  
+  fListBF->SetName("listBF");
+  fListBF->SetOwner(kTRUE);
+  dir->Add(fListBF);
+  //fListBF->Write("listBF",TObject::kSingleKey);
+
+  if(fRunShuffling) {
+    if (!fShuffledBalance) {
+      AliError("ERROR: fShuffledBalance not available");
+      return;
+    }
+    fListBFS->SetName("listBFShuffled");
+    fListBFS->SetOwner(kTRUE);
+    dir->Add(fListBFS);
+    //fListBFS->Write("listBFShuffled",TObject::kSingleKey);
+  }
+
+  if(fRunMixing) {
+    if (!fMixedBalance) {
+      AliError("ERROR: fMixedBalance not available");
+      return;
+    }
+    fListBFM->SetName("listBFMixed");
+    fListBFM->SetOwner(kTRUE);
+    dir->Add(fListBFM);
+    //fListBFM->Write("listBFMixed",TObject::kSingleKey);
+  }
+
+  dir->Write(dir->GetName(),TObject::kSingleKey);
+  gOutput->Close();
+}
+