PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskToyModel.cxx

   1 #include "TChain.h"\r
   2 #include "TList.h"\r
   3 #include "TCanvas.h"\r
   4 #include "TParticle.h"\r
   5 #include "TLorentzVector.h"\r
   6 #include "TGraphErrors.h"\r
   7 #include "TH1.h"\r
   8 #include "TH2.h"\r
   9 #include "TH3.h"\r
  10 #include "TArrayF.h"\r
  11 #include "TF1.h"\r
  12 #include "TRandom.h"\r
  13 #include "TFile.h"\r
  14 \r
  15 #include "AliAnalysisManager.h"\r
  16 #include "AliLog.h"\r
  17 \r
  18 #include "AliEventPoolManager.h"           \r
  19 \r
  20 #include "AliAnalysisTaskToyModel.h"\r
  21 #include "AliBalance.h"\r
  22 #include "AliBalancePsi.h"\r
  23 #include "AliAnalysisTaskTriggeredBF.h"\r
  24 \r
  25 \r
  26 // Analysis task for the toy model analysis\r
  27 // Authors: Panos.Christakoglou@nikhef.nl\r
  28 \r
  29 using std::cout;\r
  30 using std::endl;\r
  31 \r
  32 ClassImp(AliAnalysisTaskToyModel)\r
  33 \r
  34 //________________________________________________________________________\r
  35 AliAnalysisTaskToyModel::AliAnalysisTaskToyModel() \r
  36 : TObject(),\r
  37   fUseDebug(kFALSE),\r
  38   fBalance(0),\r
  39   fRunShuffling(kFALSE), fShuffledBalance(0),\r
  40   fRunMixing(kFALSE), fMixedBalance(0), fPoolMgr(0),\r
  41   fList(0), fListBF(0), fListBFS(0), fListBFM(0),\r
  42   fHistEventStats(0),\r
  43   fHistNumberOfAcceptedParticles(0),\r
  44   fHistReactionPlane(0),\r
  45   fHistEtaTotal(0), fHistEta(0),\r
  46   fHistEtaPhiPos(0), fHistEtaPhiNeg(0),\r
  47   fHistRapidity(0),\r
  48   fHistRapidityPions(0), fHistRapidityKaons(0), fHistRapidityProtons(0),\r
  49   fHistPhi(0),\r
  50   fHistPhiPions(0), fHistPhiKaons(0), fHistPhiProtons(0),\r
  51   fHistPt(0),\r
  52   fHistPtPions(0), fHistPtKaons(0), fHistPtProtons(0),\r
  53   fTotalMultiplicityMean(100.), fTotalMultiplicitySigma(0.0),\r
  54   fNetChargeMean(0.0), fNetChargeSigma(0.0),\r
  55   fPtMin(0.0), fPtMax(100.0),\r
  56   fEtaMin(-1.0), fEtaMax(1.0),\r
  57   fUseAcceptanceParameterization(kFALSE), fAcceptanceParameterization(0),\r
  58   fSimulateDetectorEffects(kFALSE),\r
  59   fNumberOfInefficientSectors(0),\r
  60   fInefficiencyFactorInPhi(1.0),\r
  61   fNumberOfDeadSectors(0),\r
  62   fEfficiencyDropNearEtaEdges(kFALSE),\r
  63   fEfficiencyMatrix(0),\r
  64   fUseAllCharges(kFALSE), fParticleMass(0.0),\r
  65   fPtSpectraAllCharges(0), fTemperatureAllCharges(100.),\r
  66   fReactionPlane(0.0),\r
  67   fAzimuthalAngleAllCharges(0), fDirectedFlowAllCharges(0.0), \r
  68   fEllipticFlowAllCharges(0.0), fTriangularFlowAllCharges(0.0),\r
  69   fQuandrangularFlowAllCharges(0.0), fPentangularFlowAllCharges(0.0),\r
  70   fPionPercentage(0.8), fPionMass(0.0),\r
  71   fPtSpectraPions(0), fTemperaturePions(100.),\r
  72   fAzimuthalAnglePions(0), fDirectedFlowPions(0.0), \r
  73   fEllipticFlowPions(0.0), fTriangularFlowPions(0.0), \r
  74   fQuandrangularFlowPions(0.0), fPentangularFlowPions(0.0),\r
  75   fKaonPercentage(0.8), fKaonMass(0.0),\r
  76   fPtSpectraKaons(0), fTemperatureKaons(100.),\r
  77   fAzimuthalAngleKaons(0), fDirectedFlowKaons(0.0), \r
  78   fEllipticFlowKaons(0.0), fTriangularFlowKaons(0.0),\r
  79   fQuandrangularFlowKaons(0.0), fPentangularFlowKaons(0.0),\r
  80   fProtonPercentage(0.8), fProtonMass(0.0),\r
  81   fPtSpectraProtons(0), fTemperatureProtons(100.),\r
  82   fAzimuthalAngleProtons(0), fDirectedFlowProtons(0.0), \r
  83   fEllipticFlowProtons(0.0), fTriangularFlowProtons(0.0),\r
  84   fQuandrangularFlowProtons(0.0), fPentangularFlowProtons(0.0),\r
  85   fUseDynamicalCorrelations(kFALSE), fDynamicalCorrelationsPercentage(0.1),\r
  86   fUseJets(kFALSE), fPtAssoc(0) {\r
  87   // Constructor\r
  88 }\r
  89 \r
  90 //________________________________________________________________________\r
  91 AliAnalysisTaskToyModel::~AliAnalysisTaskToyModel() {\r
  92   //Destructor\r
  93   if(fUseAllCharges) {\r
  94     delete fPtSpectraAllCharges;\r
  95     delete fAzimuthalAngleAllCharges;\r
  96   }\r
  97   else {\r
  98     delete fPtSpectraPions;\r
  99     delete fAzimuthalAnglePions;\r
 100     delete fPtSpectraKaons;\r
 101     delete fAzimuthalAngleKaons;\r
 102     delete fPtSpectraProtons;\r
 103     delete fAzimuthalAngleProtons;\r
 104   }\r
 105   if(fUseJets) delete fPtAssoc;\r
 106 }\r
 107 \r
 108 //________________________________________________________________________\r
 109 void AliAnalysisTaskToyModel::Init() {\r
 110   //Initialize objects\r
 111   //==============gRandom Seed=======================//\r
 112   gRandom->SetSeed(0);  //seed is set to a random value (depending on machine clock)\r
 113   //==============gRandom Seed=======================//\r
 114 \r
 115   //==============Particles and spectra==============//\r
 116   TParticle *pion = new TParticle();\r
 117   pion->SetPdgCode(211);\r
 118   fPionMass = pion->GetMass();\r
 119 \r
 120   TParticle *kaon = new TParticle();\r
 121   kaon->SetPdgCode(321);\r
 122   fKaonMass = kaon->GetMass();\r
 123   \r
 124   TParticle *proton = new TParticle();\r
 125   proton->SetPdgCode(2212);\r
 126   fProtonMass = proton->GetMass();\r
 127 \r
 128   if(fUseAllCharges) {\r
 129     fParticleMass = fPionMass;\r
 130     fPtSpectraAllCharges = new TF1("fPtSpectraAllCharges","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);\r
 131     fPtSpectraAllCharges->SetParName(0,"Mass");\r
 132     fPtSpectraAllCharges->SetParName(1,"Temperature");\r
 133     //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
 134     //fPtSpectraAllCharges->SetParName(0,"pt0");\r
 135     //fPtSpectraAllCharges->SetParName(1,"b");\r
 136   }\r
 137   else {\r
 138     fPtSpectraPions = new TF1("fPtSpectraPions","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);\r
 139     fPtSpectraPions->SetParName(0,"Mass");\r
 140     fPtSpectraPions->SetParName(1,"Temperature");\r
 141     \r
 142     fPtSpectraKaons = new TF1("fPtSpectraKaons","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);\r
 143     fPtSpectraKaons->SetParName(0,"Mass");\r
 144     fPtSpectraKaons->SetParName(1,"Temperature");\r
 145     \r
 146     fPtSpectraProtons = new TF1("fPtSpectraProtons","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);\r
 147     fPtSpectraProtons->SetParName(0,"Mass");\r
 148     fPtSpectraProtons->SetParName(1,"Temperature");\r
 149   }\r
 150   //==============Particles and spectra==============//\r
 151 \r
 152   //==============Flow values==============//\r
 153   if(fUseAllCharges) {\r
 154     if(fUseDebug) {\r
 155       Printf("Directed flow: %lf",fDirectedFlowAllCharges);\r
 156       Printf("Elliptic flow: %lf",fEllipticFlowAllCharges);\r
 157       Printf("Triangular flow: %lf",fTriangularFlowAllCharges);\r
 158       Printf("Quandrangular flow: %lf",fQuandrangularFlowAllCharges);\r
 159       Printf("Pentangular flow: %lf",fPentangularFlowAllCharges);\r
 160     }\r
 161 \r
 162     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
 163     fAzimuthalAngleAllCharges->SetParName(0,"Reaction Plane");\r
 164     fAzimuthalAngleAllCharges->SetParName(1,"Directed flow");\r
 165     fAzimuthalAngleAllCharges->SetParName(2,"Elliptic flow"); \r
 166     fAzimuthalAngleAllCharges->SetParName(3,"Triangular flow");\r
 167     fAzimuthalAngleAllCharges->SetParName(4,"Quandrangular flow");\r
 168     fAzimuthalAngleAllCharges->SetParName(5,"Pentangular flow");\r
 169   }\r
 170   else {\r
 171     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
 172     fAzimuthalAnglePions->SetParName(0,"Reaction Plane");\r
 173     fAzimuthalAnglePions->SetParName(1,"Directed flow");\r
 174     fAzimuthalAnglePions->SetParName(2,"Elliptic flow"); \r
 175     fAzimuthalAnglePions->SetParName(3,"Triangular flow");\r
 176     fAzimuthalAnglePions->SetParName(4,"Quandrangular flow");\r
 177     fAzimuthalAnglePions->SetParName(5,"Pentangular flow");\r
 178     \r
 179     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
 180     fAzimuthalAngleKaons->SetParName(0,"Reaction Plane");\r
 181     fAzimuthalAngleKaons->SetParName(1,"Directed flow");\r
 182     fAzimuthalAngleKaons->SetParName(2,"Elliptic flow"); \r
 183     fAzimuthalAngleKaons->SetParName(3,"Triangular flow");\r
 184     fAzimuthalAngleKaons->SetParName(4,"Quandrangular flow");\r
 185     fAzimuthalAngleKaons->SetParName(5,"Pentangular flow");\r
 186     \r
 187     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
 188     fAzimuthalAngleProtons->SetParName(0,"Reaction Plane");\r
 189     fAzimuthalAngleProtons->SetParName(1,"Directed flow");\r
 190     fAzimuthalAngleProtons->SetParName(2,"Elliptic flow"); \r
 191     fAzimuthalAngleProtons->SetParName(3,"Triangular flow");\r
 192     fAzimuthalAngleProtons->SetParName(4,"Quandrangular flow");\r
 193     fAzimuthalAngleProtons->SetParName(5,"Pentangular flow");\r
 194   }\r
 195   //==============Flow values==============//\r
 196 \r
 197   //===================Jets===================//\r
 198   if(fUseJets) {\r
 199     fPtAssoc = new TF1("fPtAssoc","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,20.);\r
 200     fPtAssoc->SetParName(0,"pt0");\r
 201     fPtAssoc->SetParName(1,"b");\r
 202     fPtAssoc->SetParameter(0,0.139);\r
 203     fPtAssoc->SetParameter(1,0.5);\r
 204     fPtAssoc->SetLineColor(1);\r
 205   }\r
 206 \r
 207   //===================Jets===================//\r
 208 \r
 209   //==============Efficiency matrix==============//\r
 210   if(fSimulateDetectorEffects) SetupEfficiencyMatrix();\r
 211   //==============Efficiency matrix==============//\r
 212 }\r
 213 \r
 214 //________________________________________________________________________\r
 215 void AliAnalysisTaskToyModel::SetupEfficiencyMatrix() {\r
 216   //Setup the efficiency matrix\r
 217   TH1F *hPt = new TH1F("hPt","",200,0.1,20.1);\r
 218   TH1F *hEta = new TH1F("hEta","",20,-0.95,0.95);\r
 219   TH1F *hPhi = new TH1F("hPhi","",72,0.,2.*TMath::Pi());\r
 220   fEfficiencyMatrix = new TH3F("fEfficiencyMatrix","",\r
 221                                hEta->GetNbinsX(),\r
 222                                hEta->GetXaxis()->GetXmin(),\r
 223                                hEta->GetXaxis()->GetXmax(),\r
 224                                hPt->GetNbinsX(),\r
 225                                hPt->GetXaxis()->GetXmin(),\r
 226                                hPt->GetXaxis()->GetXmax(),\r
 227                                hPhi->GetNbinsX(),\r
 228                                hPhi->GetXaxis()->GetXmin(),\r
 229                                hPhi->GetXaxis()->GetXmax());\r
 230 \r
 231   //Efficiency in pt\r
 232   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
 233   for(Int_t i=1;i<=20;i++) {\r
 234     hPt->SetBinContent(i,epsilon[i-1]);\r
 235     hPt->SetBinError(i,0.01);\r
 236   }\r
 237   for(Int_t i=21;i<=200;i++) {\r
 238     hPt->SetBinContent(i,epsilon[19]);\r
 239     hPt->SetBinError(i,0.01);\r
 240   }\r
 241 \r
 242   //Efficiency in eta\r
 243   for(Int_t i=1;i<=hEta->GetNbinsX();i++) {\r
 244     hEta->SetBinContent(i,1.0);\r
 245     hEta->SetBinError(i,0.01);\r
 246   }\r
 247   if(fEfficiencyDropNearEtaEdges) {\r
 248     hEta->SetBinContent(1,0.7); hEta->SetBinContent(2,0.8);\r
 249     hEta->SetBinContent(3,0.9);\r
 250     hEta->SetBinContent(18,0.9); hEta->SetBinContent(19,0.8);\r
 251     hEta->SetBinContent(20,0.7);\r
 252   }\r
 253 \r
 254   //Efficiency in phi\r
 255   for(Int_t i=1;i<=hPhi->GetNbinsX();i++) {\r
 256     hPhi->SetBinContent(i,1.0);\r
 257     hPhi->SetBinError(i,0.01);\r
 258   }\r
 259   for(Int_t i=1;i<=fNumberOfInefficientSectors;i++)\r
 260     hPhi->SetBinContent(hPhi->FindBin(hPhi->GetRandom()),fInefficiencyFactorInPhi);\r
 261   for(Int_t i=1;i<=fNumberOfDeadSectors;i++)\r
 262     hPhi->SetBinContent(hPhi->FindBin(hPhi->GetRandom()),0.0);\r
 263   \r
 264   //Fill the 3D efficiency map\r
 265   for(Int_t iBinX = 1; iBinX<=fEfficiencyMatrix->GetXaxis()->GetNbins();iBinX++) {\r
 266     //cout<<"==================================="<<endl;\r
 267     for(Int_t iBinY = 1; iBinY<=fEfficiencyMatrix->GetYaxis()->GetNbins();iBinY++) {\r
 268       //cout<<"==================================="<<endl;\r
 269       for(Int_t iBinZ = 1; iBinZ<=fEfficiencyMatrix->GetZaxis()->GetNbins();iBinZ++) {\r
 270         fEfficiencyMatrix->SetBinContent(iBinX,iBinY,iBinZ,hEta->GetBinContent(iBinX)*hPt->GetBinContent(iBinY)*hPhi->GetBinContent(iBinZ));\r
 271         //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
 272       }\r
 273     }\r
 274   }\r
 275 }\r
 276 \r
 277 //________________________________________________________________________\r
 278 void AliAnalysisTaskToyModel::CreateOutputObjects() {\r
 279   // Create histograms\r
 280   // Called once\r
 281 \r
 282   // global switch disabling the reference \r
 283   // (to avoid "Replacing existing TH1" if several wagons are created in train)\r
 284   Bool_t oldStatus = TH1::AddDirectoryStatus();\r
 285   TH1::AddDirectory(kFALSE);\r
 286 \r
 287   if(!fBalance) {\r
 288     fBalance = new AliBalancePsi();\r
 289     fBalance->SetDeltaEtaMax(2.0);\r
 290   }\r
 291   if(fRunShuffling) {\r
 292     if(!fShuffledBalance) {\r
 293       fShuffledBalance = new AliBalancePsi();\r
 294       fShuffledBalance->SetDeltaEtaMax(2.0);\r
 295     }\r
 296   }\r
 297   if(fRunMixing) {\r
 298     if(!fMixedBalance) {\r
 299       fMixedBalance = new AliBalancePsi();\r
 300       fMixedBalance->SetDeltaEtaMax(2.0);\r
 301     }\r
 302   }\r
 303 \r
 304   //QA list\r
 305   fList = new TList();\r
 306   fList->SetName("listQA");\r
 307   fList->SetOwner();\r
 308 \r
 309   //Balance Function list\r
 310   fListBF = new TList();\r
 311   fListBF->SetName("listBF");\r
 312   fListBF->SetOwner();\r
 313 \r
 314   if(fRunShuffling) {\r
 315     fListBFS = new TList();\r
 316     fListBFS->SetName("listBFShuffled");\r
 317     fListBFS->SetOwner();\r
 318   }\r
 319 \r
 320   if(fRunMixing) {\r
 321     fListBFM = new TList();\r
 322     fListBFM->SetName("listBFMixed");\r
 323     fListBFM->SetOwner();\r
 324   }\r
 325 \r
 326   //==============QA================//\r
 327   //Event stats.\r
 328   TString gCutName[4] = {"Total","Offline trigger",\r
 329                          "Vertex","Analyzed"};\r
 330   fHistEventStats = new TH1F("fHistEventStats",\r
 331                              "Event statistics;;N_{events}",\r
 332                              4,0.5,4.5);\r
 333   for(Int_t i = 1; i <= 4; i++)\r
 334     fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());\r
 335   fList->Add(fHistEventStats);\r
 336 \r
 337   fHistNumberOfAcceptedParticles = new TH1F("fHistNumberOfAcceptedParticles",";N_{acc.};Entries",10000,0,10000);\r
 338   fList->Add(fHistNumberOfAcceptedParticles);\r
 339   \r
 340   fHistReactionPlane = new TH1F("fHistReactionPlane","Reaction plane angle;#Psi [rad];Entries",1000,0.,2.*TMath::Pi());\r
 341   fList->Add(fHistReactionPlane);\r
 342 \r
 343   //Pseudo-rapidity\r
 344   fHistEtaTotal = new TH1F("fHistEtaTotal","Pseudo-rapidity (full phase space);#eta;Entries",1000,-15.,15.); \r
 345   fList->Add(fHistEtaTotal);\r
 346 \r
 347   fHistEta = new TH1F("fHistEta","Pseudo-rapidity (acceptance);#eta;Entries",1000,-1.5,1.5); \r
 348   fList->Add(fHistEta);\r
 349 \r
 350   fHistEtaPhiPos = new TH2F("fHistEtaPhiPos","#eta-#phi distribution (+);#eta;#varphi (rad)",80,-2.,2.,72,-TMath::Pi()/2.,3.*TMath::Pi()/2.);\r
 351   fList->Add(fHistEtaPhiPos);\r
 352   fHistEtaPhiNeg = new TH2F("fHistEtaPhiNeg","#eta-#phi distribution (-);#eta;#varphi (rad)",80,-2.,2.,72,-TMath::Pi()/2.,3.*TMath::Pi()/2.);\r
 353   fList->Add(fHistEtaPhiNeg);\r
 354 \r
 355   //Rapidity\r
 356   fHistRapidity = new TH1F("fHistRapidity","Rapidity (acceptance);y;Entries",1000,-1.5,1.5); \r
 357   fList->Add(fHistRapidity);\r
 358   fHistRapidityPions = new TH1F("fHistRapidityPions","Rapidity (acceptance - pions);y;Entries",1000,-1.5,1.5); \r
 359   fList->Add(fHistRapidityPions);\r
 360   fHistRapidityKaons = new TH1F("fHistRapidityKaons","Rapidity (acceptance - kaons);y;Entries",1000,-1.5,1.5); \r
 361   fList->Add(fHistRapidityKaons);\r
 362   fHistRapidityProtons = new TH1F("fHistRapidityProtons","Rapidity (acceptance - protons);y;Entries",1000,-1.5,1.5); \r
 363   fList->Add(fHistRapidityProtons);\r
 364 \r
 365   //Phi\r
 366   fHistPhi = new TH1F("fHistPhi","Phi (acceptance);#phi (rad);Entries",1000,0.,2*TMath::Pi());\r
 367   fList->Add(fHistPhi);\r
 368 \r
 369   fHistPhiPions = new TH1F("fHistPhiPions","Phi (acceptance - pions);#phi (rad);Entries",1000,0.,2*TMath::Pi());\r
 370   fList->Add(fHistPhiPions);\r
 371   fHistPhiKaons = new TH1F("fHistPhiKaons","Phi (acceptance - kaons);#phi (rad);Entries",1000,0.,2*TMath::Pi());\r
 372   fList->Add(fHistPhiKaons);\r
 373   fHistPhiProtons = new TH1F("fHistPhiProtons","Phi (acceptance - protons);#phi (rad);Entries",1000,0.,2*TMath::Pi());\r
 374   fList->Add(fHistPhiProtons);\r
 375 \r
 376   //Pt\r
 377   fHistPt = new TH1F("fHistPt","Pt (acceptance);p_{t} (GeV/c);Entries",1000,0.,10.);\r
 378   fList->Add(fHistPt);\r
 379 \r
 380   fHistPtPions = new TH1F("fHistPtPions","Pt (acceptance - pions);p_{t} (GeV/c);Entries",1000,0.,10.);\r
 381   fList->Add(fHistPtPions);\r
 382   fHistPtKaons = new TH1F("fHistPtKaons","Pt (acceptance - kaons);p_{t} (GeV/c);Entries",1000,0.,10.);\r
 383   fList->Add(fHistPtKaons);\r
 384   fHistPtProtons = new TH1F("fHistPtProtons","Pt (acceptance - protons);p_{t} (GeV/c);Entries",1000,0.,10.);\r
 385   fList->Add(fHistPtProtons);\r
 386 \r
 387   if(fEfficiencyMatrix) fList->Add(fEfficiencyMatrix);\r
 388 \r
 389   //==============Balance function histograms================//\r
 390   // Initialize histograms if not done yet\r
 391   if(!fBalance->GetHistNp()){\r
 392     AliWarning("Histograms not yet initialized! --> Will be done now");\r
 393     AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
 394     fBalance->InitHistograms();\r
 395   }\r
 396 \r
 397   if(fRunShuffling) {\r
 398     if(!fShuffledBalance->GetHistNp()) {\r
 399       AliWarning("Histograms (shuffling) not yet initialized! --> Will be done now");\r
 400       AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
 401       fShuffledBalance->InitHistograms();\r
 402     }\r
 403   }\r
 404 \r
 405   fListBF->Add(fBalance->GetHistNp());\r
 406   fListBF->Add(fBalance->GetHistNn());\r
 407   fListBF->Add(fBalance->GetHistNpn());\r
 408   fListBF->Add(fBalance->GetHistNnn());\r
 409   fListBF->Add(fBalance->GetHistNpp());\r
 410   fListBF->Add(fBalance->GetHistNnp());\r
 411 \r
 412   if(fRunShuffling) {\r
 413     fListBFS->Add(fShuffledBalance->GetHistNp());\r
 414     fListBFS->Add(fShuffledBalance->GetHistNn());\r
 415     fListBFS->Add(fShuffledBalance->GetHistNpn());\r
 416     fListBFS->Add(fShuffledBalance->GetHistNnn());\r
 417     fListBFS->Add(fShuffledBalance->GetHistNpp());\r
 418     fListBFS->Add(fShuffledBalance->GetHistNnp());\r
 419   }\r
 420 \r
 421   if(fRunMixing) {\r
 422     fListBFM->Add(fMixedBalance->GetHistNp());\r
 423     fListBFM->Add(fMixedBalance->GetHistNn());\r
 424     fListBFM->Add(fMixedBalance->GetHistNpn());\r
 425     fListBFM->Add(fMixedBalance->GetHistNnn());\r
 426     fListBFM->Add(fMixedBalance->GetHistNpp());\r
 427     fListBFM->Add(fMixedBalance->GetHistNnp());\r
 428   }\r
 429 \r
 430   // Event Mixing\r
 431   if(fRunMixing){\r
 432     Int_t fMixingTracks = 2000;\r
 433     Int_t trackDepth = fMixingTracks; \r
 434     Int_t poolsize   = 1000;  // Maximum number of events, ignored in the present implemented of AliEventPoolManager\r
 435     \r
 436     // centrality bins\r
 437     // 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
 438     // Double_t* centbins        = centralityBins;\r
 439     // Int_t nCentralityBins     = sizeof(centralityBins) / sizeof(Double_t) - 1;\r
 440 \r
 441     // multiplicity bins\r
 442     Double_t multiplicityBins[] = {0,10,20,30,40,50,60,70,80,100,100000}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
 443     Double_t* multbins        = multiplicityBins;\r
 444     Int_t nMultiplicityBins     = sizeof(multiplicityBins) / sizeof(Double_t) - 1;\r
 445     \r
 446     // Zvtx bins\r
 447     Double_t vertexBins[] = {-10., -7., -5., -3., -1., 1., 3., 5., 7., 10.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
 448     Double_t* vtxbins     = vertexBins;\r
 449     Int_t nVertexBins     = sizeof(vertexBins) / sizeof(Double_t) - 1;\r
 450     \r
 451     // Event plane angle (Psi) bins\r
 452     // Double_t psiBins[] = {0.,45.,135.,215.,305.,360.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
 453     // Double_t* psibins     = psiBins;\r
 454     // Int_t nPsiBins     = sizeof(psiBins) / sizeof(Double_t) - 1;\r
 455     \r
 456     // // run the event mixing also in bins of event plane (statistics!)\r
 457     // if(fRunMixingEventPlane){\r
 458     //   if(fEventClass=="Multiplicity"){\r
 459     //  fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins, nPsiBins, psibins);\r
 460     //   }\r
 461     //   else{\r
 462     //  fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins, nPsiBins, psibins);\r
 463     //   }\r
 464     // }\r
 465     // else{\r
 466     //if(fEventClass=="Multiplicity"){\r
 467     fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins);\r
 468     //}\r
 469     //else{\r
 470     //fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);\r
 471     //}\r
 472   }\r
 473 \r
 474   TH1::AddDirectory(oldStatus);\r
 475 }\r
 476 \r
 477 //________________________________________________________________________\r
 478 void AliAnalysisTaskToyModel::Run(Int_t nEvents) {\r
 479   // Main loop\r
 480   // Called for each event\r
 481   Short_t vCharge = 0;\r
 482   Float_t vY = 0.0;\r
 483   Float_t vEta = 0.0;\r
 484   Float_t vPhi = 0.0;\r
 485   Float_t vP[3] = {0.,0.,0.};\r
 486   Float_t vPt = 0.0;\r
 487   Float_t vE = 0.0;\r
 488   Bool_t isPion = kFALSE, isKaon = kFALSE, isProton = kFALSE;\r
 489 \r
 490   Double_t gDecideCharge = 0.;\r
 491 \r
 492   if(fUseAllCharges) {\r
 493     //fPtSpectraAllCharges->SetParameter(0,fParticleMass);\r
 494     fPtSpectraAllCharges->SetParameter(0,1.05);\r
 495     fPtSpectraAllCharges->SetParameter(1,fTemperatureAllCharges);\r
 496 \r
 497     fAzimuthalAngleAllCharges->SetParameter(1,fDirectedFlowAllCharges);\r
 498     fAzimuthalAngleAllCharges->SetParameter(2,fEllipticFlowAllCharges);\r
 499     fAzimuthalAngleAllCharges->SetParameter(3,fTriangularFlowAllCharges);\r
 500     fAzimuthalAngleAllCharges->SetParameter(4,fQuandrangularFlowAllCharges);\r
 501     fAzimuthalAngleAllCharges->SetParameter(5,fPentangularFlowAllCharges);\r
 502   }\r
 503   else {\r
 504     fPtSpectraPions->SetParameter(0,fPionMass);\r
 505     fPtSpectraPions->SetParameter(1,fTemperaturePions);\r
 506     fPtSpectraKaons->SetParameter(0,fKaonMass);\r
 507     fPtSpectraKaons->SetParameter(1,fTemperatureKaons);\r
 508     fPtSpectraProtons->SetParameter(0,fProtonMass);\r
 509     fPtSpectraProtons->SetParameter(1,fTemperatureProtons);\r
 510 \r
 511     fAzimuthalAnglePions->SetParameter(1,fDirectedFlowPions);\r
 512     fAzimuthalAnglePions->SetParameter(2,fEllipticFlowPions);\r
 513     fAzimuthalAnglePions->SetParameter(3,fTriangularFlowPions);\r
 514     fAzimuthalAnglePions->SetParameter(4,fQuandrangularFlowPions);\r
 515     fAzimuthalAnglePions->SetParameter(5,fPentangularFlowPions);\r
 516 \r
 517     fAzimuthalAngleKaons->SetParameter(1,fDirectedFlowKaons);\r
 518     fAzimuthalAngleKaons->SetParameter(2,fEllipticFlowKaons);\r
 519     fAzimuthalAngleKaons->SetParameter(3,fTriangularFlowKaons);\r
 520     fAzimuthalAngleKaons->SetParameter(4,fQuandrangularFlowKaons);\r
 521     fAzimuthalAngleKaons->SetParameter(5,fPentangularFlowKaons);\r
 522 \r
 523     fAzimuthalAngleProtons->SetParameter(1,fDirectedFlowProtons);\r
 524     fAzimuthalAngleProtons->SetParameter(2,fEllipticFlowProtons);\r
 525     fAzimuthalAngleProtons->SetParameter(3,fTriangularFlowProtons);\r
 526     fAzimuthalAngleProtons->SetParameter(4,fQuandrangularFlowProtons);\r
 527     fAzimuthalAngleProtons->SetParameter(5,fPentangularFlowProtons);\r
 528   }\r
 529 \r
 530 \r
 531   for(Int_t iEvent = 0; iEvent < nEvents; iEvent++) {\r
 532     // vector holding the charges/kinematics of all tracks (charge,y,eta,phi,p0,p1,p2,pt,E)\r
 533     //vector<Double_t> *chargeVectorShuffle[9];   // this will be shuffled\r
 534     //vector<Double_t> *chargeVector[9];          // original charge\r
 535     //for(Int_t i = 0; i < 9; i++){\r
 536     //chargeVectorShuffle[i] = new vector<Double_t>;\r
 537     //chargeVector[i]        = new vector<Double_t>;\r
 538     //}\r
 539 \r
 540     // TObjArray for the accepted particles \r
 541     // (has to be done here, otherwise mxing with event pool does not work, overwriting pointers!)\r
 542     TObjArray *tracksMain = new TObjArray();\r
 543     tracksMain->SetOwner(kTRUE);\r
 544     TObjArray *tracksMixing = 0x0;\r
 545     if(fRunMixing) {\r
 546       tracksMixing = new TObjArray();\r
 547       tracksMixing->SetOwner(kTRUE);\r
 548     }\r
 549 \r
 550     tracksMain->Clear();\r
 551     if(fRunMixing) tracksMixing->Clear();\r
 552     \r
 553     fHistEventStats->Fill(1);\r
 554     fHistEventStats->Fill(2);\r
 555     fHistEventStats->Fill(3);\r
 556 \r
 557     if((iEvent%1000) == 0) \r
 558       cout<<"Event: "<<iEvent<<"/"<<nEvents<<endl;\r
 559     \r
 560     //Multiplicities\r
 561     Int_t nMultiplicity = (Int_t)(gRandom->Gaus(fTotalMultiplicityMean,fTotalMultiplicitySigma));\r
 562     Int_t nNetCharge = (Int_t)(gRandom->Gaus(fNetChargeMean,fNetChargeSigma));\r
 563     \r
 564     Int_t nGeneratedPositive = (Int_t)((nMultiplicity/2) + nNetCharge);\r
 565     Int_t nGeneratedNegative = nMultiplicity - nGeneratedPositive;\r
 566     if(fUseDebug) \r
 567       Printf("Total multiplicity: %d - Generated positive: %d - Generated negative: %d",nMultiplicity,nGeneratedPositive,nGeneratedNegative);\r
 568     \r
 569     //Randomization of the reaction plane\r
 570     fReactionPlane = 2.0*TMath::Pi()*gRandom->Rndm();\r
 571     //fReactionPlane = 0.0;\r
 572     if(fUseAllCharges) \r
 573       fAzimuthalAngleAllCharges->SetParameter(0,fReactionPlane);\r
 574     else {\r
 575       fAzimuthalAnglePions->SetParameter(0,fReactionPlane);\r
 576       fAzimuthalAngleKaons->SetParameter(0,fReactionPlane);\r
 577       fAzimuthalAngleProtons->SetParameter(0,fReactionPlane);\r
 578     }\r
 579     \r
 580     Int_t gNumberOfAcceptedParticles = 0;\r
 581     Int_t gNumberOfAcceptedPositiveParticles = 0;\r
 582     Int_t gNumberOfAcceptedNegativeParticles = 0;\r
 583     Int_t nGeneratedPions = 0, nGeneratedKaons = 0, nGeneratedProtons = 0;\r
 584  \r
 585     //Generate particles\r
 586     for(Int_t iParticleCount = 0; iParticleCount < nMultiplicity; iParticleCount++) {\r
 587       isPion = kFALSE; isKaon = kFALSE; isProton = kFALSE;\r
 588       if(fUseDebug) \r
 589         Printf("Generating positive: %d(%d)",iParticleCount+1,nGeneratedPositive);\r
 590 \r
 591       //Pseudo-rapidity sampled from a Gaussian centered @ 0\r
 592       vEta = gRandom->Gaus(0.0,4.0);\r
 593 \r
 594       //Fill QA histograms (full phase space)\r
 595       fHistEtaTotal->Fill(vEta);\r
 596 \r
 597       //Decide the charge\r
 598       gDecideCharge = gRandom->Rndm();\r
 599       if(gDecideCharge <= 1.*nGeneratedPositive/nMultiplicity)       \r
 600         vCharge = 1;\r
 601       else \r
 602         vCharge = -1;\r
 603       \r
 604       //Acceptance\r
 605       if((vEta < fEtaMin) || (vEta > fEtaMax)) continue;\r
 606 \r
 607       if(!fUseAllCharges) {\r
 608         //Decide the specie\r
 609         Double_t randomNumberSpecies = gRandom->Rndm();\r
 610         if((randomNumberSpecies >= 0.0)&&(randomNumberSpecies < fPionPercentage)) {\r
 611           nGeneratedPions += 1;\r
 612           vPt = fPtSpectraPions->GetRandom();\r
 613           vPhi = fAzimuthalAnglePions->GetRandom();\r
 614           fParticleMass = fPionMass;\r
 615           isPion = kTRUE;\r
 616         }\r
 617         else if((randomNumberSpecies >= fPionPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage)) {\r
 618           nGeneratedKaons += 1;\r
 619           vPt = fPtSpectraKaons->GetRandom();\r
 620           vPhi = fAzimuthalAngleKaons->GetRandom();\r
 621           fParticleMass = fKaonMass;\r
 622           isKaon = kTRUE;\r
 623         }\r
 624         else if((randomNumberSpecies >= fPionPercentage + fKaonPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage + fProtonPercentage)) {\r
 625           nGeneratedProtons += 1;\r
 626           vPt = fPtSpectraProtons->GetRandom();\r
 627           vPhi = fAzimuthalAngleProtons->GetRandom();\r
 628           fParticleMass = fProtonMass;\r
 629           isProton = kTRUE;\r
 630         }\r
 631       }\r
 632       else {\r
 633         vPt = fPtSpectraAllCharges->GetRandom();\r
 634         vPhi = fAzimuthalAngleAllCharges->GetRandom();\r
 635       }\r
 636       \r
 637       vP[0] = vPt*TMath::Cos(vPhi);\r
 638       vP[1] = vPt*TMath::Sin(vPhi);\r
 639       vP[2] = vPt*TMath::SinH(vEta);\r
 640       vE = TMath::Sqrt(TMath::Power(fParticleMass,2) +\r
 641                         TMath::Power(vP[0],2) +\r
 642                         TMath::Power(vP[1],2) +\r
 643                         TMath::Power(vP[2],2));\r
 644       \r
 645       vY = 0.5*TMath::Log((vE + vP[2])/(vE - vP[2]));\r
 646       \r
 647       //pt coverage\r
 648       if((vPt < fPtMin) || (vPt > fPtMax)) continue;\r
 649       //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);\r
 650 \r
 651       //acceptance filter\r
 652       if(fUseAcceptanceParameterization) {\r
 653         Double_t gRandomNumberForAcceptance = gRandom->Rndm();\r
 654         if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(vPt)) \r
 655           continue;\r
 656       }\r
 657 \r
 658       //Detector effects\r
 659       if(fSimulateDetectorEffects) {\r
 660         Double_t randomNumber = gRandom->Rndm();\r
 661         if(randomNumber > fEfficiencyMatrix->GetBinContent(fEfficiencyMatrix->FindBin(vEta,vPt,vPhi)))\r
 662           continue;\r
 663       }\r
 664 \r
 665       //Fill QA histograms (acceptance);\r
 666       if(vCharge > 0) {\r
 667         gNumberOfAcceptedPositiveParticles += 1;\r
 668         if(vPhi > 3.*TMath::Pi()/2.)\r
 669           fHistEtaPhiPos->Fill(vEta,vPhi-2.*TMath::Pi());\r
 670         else\r
 671           fHistEtaPhiPos->Fill(vEta,vPhi);\r
 672       }\r
 673       else {\r
 674         gNumberOfAcceptedNegativeParticles += 1;\r
 675         if(vPhi > 3.*TMath::Pi()/2.)\r
 676           fHistEtaPhiNeg->Fill(vEta,vPhi-2.*TMath::Pi());\r
 677         else\r
 678           fHistEtaPhiNeg->Fill(vEta,vPhi);\r
 679       }\r
 680 \r
 681       fHistEta->Fill(vEta);\r
 682       fHistRapidity->Fill(vY);\r
 683       fHistPhi->Fill(vPhi);\r
 684       fHistPt->Fill(vPt);\r
 685       if(isPion) {\r
 686         fHistRapidityPions->Fill(vY);\r
 687         fHistPhiPions->Fill(vPhi);\r
 688         fHistPtPions->Fill(vPt);\r
 689       }\r
 690       else if(isKaon) {\r
 691         fHistRapidityKaons->Fill(vY);\r
 692         fHistPhiKaons->Fill(vPhi);\r
 693         fHistPtKaons->Fill(vPt);\r
 694       }\r
 695       else if(isProton) {\r
 696         fHistRapidityProtons->Fill(vY);\r
 697         fHistPhiProtons->Fill(vPhi);\r
 698         fHistPtProtons->Fill(vPt);\r
 699       }\r
 700       gNumberOfAcceptedParticles += 1;\r
 701 \r
 702       // add the track to the TObjArray\r
 703       tracksMain->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, 1.0));  \r
 704       if(fRunMixing) \r
 705         tracksMixing->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, 1.0));  \r
 706     }//generated positive particle loop\r
 707 \r
 708     //Jets\r
 709     if(fUseJets) {\r
 710       const Int_t nAssociated = 3;\r
 711 \r
 712       Double_t gPtTrig1 = 0., gPtTrig2 = 0., gPtAssoc = 0.;\r
 713       Double_t gPhiTrig1 = 0., gPhiTrig2 = 0., gPhiAssoc = 0.;\r
 714       Double_t gEtaTrig1 = 0., gEtaTrig2 = 0., gEtaAssoc = 0.;\r
 715       Short_t gChargeTrig1 = 0, gChargeTrig2 = 0, gChargeAssoc = 0;\r
 716  \r
 717       Double_t gJetCone = 0.2;\r
 718 \r
 719       //First leading particle\r
 720       gPtTrig1 = gRandom->Uniform(3.,5.);\r
 721       gEtaTrig1 = gRandom->Uniform(-0.8,0.8);\r
 722       gPhiTrig1 = gRandom->Uniform(0.,TMath::TwoPi());\r
 723 \r
 724       //Decide the charge\r
 725       gDecideCharge = gRandom->Rndm();\r
 726       if(gDecideCharge <= 0.5)\r
 727         gChargeTrig1 = 1;\r
 728       else \r
 729         gChargeTrig1 = -1;\r
 730       \r
 731       //Acceptance\r
 732       if((gEtaTrig1 < fEtaMin) || (gEtaTrig1 > fEtaMax)) continue;\r
 733       //pt coverage\r
 734       if((gPtTrig1 < fPtMin) || (gPtTrig1 > fPtMax)) continue;\r
 735       //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);\r
 736 \r
 737       //acceptance filter\r
 738       if(fUseAcceptanceParameterization) {\r
 739         Double_t gRandomNumberForAcceptance = gRandom->Rndm();\r
 740         if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(gPtTrig1)) \r
 741           continue;\r
 742       }\r
 743 \r
 744       //Detector effects\r
 745       if(fSimulateDetectorEffects) {\r
 746         Double_t randomNumber = gRandom->Rndm();\r
 747         if(randomNumber > fEfficiencyMatrix->GetBinContent(fEfficiencyMatrix->FindBin(gEtaTrig1,gPtTrig1,gPhiTrig1)))\r
 748           continue;\r
 749       }\r
 750 \r
 751       gNumberOfAcceptedParticles += 1;\r
 752 \r
 753       // add the track to the TObjArray\r
 754       tracksMain->Add(new AliBFBasicParticle(gEtaTrig1, gPhiTrig1, gPtTrig1, gChargeTrig1, 1.0));  \r
 755       if(fRunMixing) \r
 756         tracksMixing->Add(new AliBFBasicParticle(gEtaTrig1, gPhiTrig1, gPtTrig1, gChargeTrig1, 1.0));  \r
 757 \r
 758       Int_t iAssociated = 0; \r
 759       while(iAssociated < nAssociated) {\r
 760         gPtAssoc = fPtAssoc->GetRandom();\r
 761         if(gPtAssoc < gPtTrig1) {\r
 762           gEtaAssoc = gRandom->Uniform(gEtaTrig1 - gJetCone/2.,gEtaTrig1 + gJetCone/2.);\r
 763           gPhiAssoc = gRandom->Uniform(gPhiTrig1 - gJetCone/2.,gPhiTrig1 + gJetCone/2.);\r
 764           if(gPhiAssoc < 0.) gPhiAssoc += TMath::TwoPi();\r
 765           else if(gPhiAssoc > TMath::TwoPi()) gPhiAssoc -= TMath::TwoPi();\r
 766           \r
 767           iAssociated += 1;\r
 768 \r
 769           gDecideCharge = gRandom->Rndm();\r
 770           if(gDecideCharge <= 0.5)\r
 771             gChargeAssoc = 1;\r
 772           else \r
 773             gChargeAssoc = -1;\r
 774           \r
 775           //Acceptance\r
 776           if((gEtaAssoc < fEtaMin) || (gEtaAssoc > fEtaMax)) continue;\r
 777           //pt coverage\r
 778           if((gPtAssoc < fPtMin) || (gPtAssoc > fPtMax)) continue;\r
 779           //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);\r
 780 \r
 781           //acceptance filter\r
 782           if(fUseAcceptanceParameterization) {\r
 783             Double_t gRandomNumberForAcceptance = gRandom->Rndm();\r
 784             if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(gPtAssoc)) \r
 785               continue;\r
 786           }\r
 787 \r
 788           //Detector effects\r
 789           if(fSimulateDetectorEffects) {\r
 790             Double_t randomNumber = gRandom->Rndm();\r
 791             if(randomNumber > fEfficiencyMatrix->GetBinContent(fEfficiencyMatrix->FindBin(gEtaAssoc,gPtAssoc,gPhiAssoc)))\r
 792               continue;\r
 793           }\r
 794 \r
 795           gNumberOfAcceptedParticles += 1;\r
 796 \r
 797           // add the track to the TObjArray\r
 798           tracksMain->Add(new AliBFBasicParticle(gEtaAssoc, gPhiAssoc, gPtAssoc, gChargeAssoc, 1.0));  \r
 799           if(fRunMixing) \r
 800             tracksMixing->Add(new AliBFBasicParticle(gEtaAssoc, gPhiAssoc, gPtAssoc, gChargeAssoc, 1.0));  \r
 801         }//pt,assoc < pt,trig\r
 802       }//associated\r
 803       \r
 804       //back2back\r
 805       gPtTrig2 = gPtTrig1;\r
 806       gEtaTrig2 = -gEtaTrig1;\r
 807       gPhiTrig2 = TMath::Pi() + gPhiTrig1;\r
 808       if(gPhiTrig2 < 0.) gPhiTrig2 += TMath::TwoPi();\r
 809       else if(gPhiTrig2 > TMath::TwoPi()) gPhiTrig2 -= TMath::TwoPi();\r
 810 \r
 811       //Decide the charge\r
 812       gDecideCharge = gRandom->Rndm();\r
 813       if(gDecideCharge <= 0.5)\r
 814         gChargeTrig2 = 1;\r
 815       else \r
 816         gChargeTrig2 = -1;\r
 817       \r
 818       //Acceptance\r
 819       if((gEtaTrig2 < fEtaMin) || (gEtaTrig2 > fEtaMax)) continue;\r
 820       //pt coverage\r
 821       if((gPtTrig2 < fPtMin) || (gPtTrig2 > fPtMax)) continue;\r
 822       //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);\r
 823 \r
 824       //acceptance filter\r
 825       if(fUseAcceptanceParameterization) {\r
 826         Double_t gRandomNumberForAcceptance = gRandom->Rndm();\r
 827         if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(gPtTrig2)) \r
 828           continue;\r
 829       }\r
 830 \r
 831       //Detector effects\r
 832       if(fSimulateDetectorEffects) {\r
 833         Double_t randomNumber = gRandom->Rndm();\r
 834         if(randomNumber > fEfficiencyMatrix->GetBinContent(fEfficiencyMatrix->FindBin(gEtaTrig2,gPtTrig2,gPhiTrig2)))\r
 835           continue;\r
 836       }\r
 837 \r
 838       gNumberOfAcceptedParticles += 1;\r
 839 \r
 840       // add the track to the TObjArray\r
 841       tracksMain->Add(new AliBFBasicParticle(gEtaTrig2, gPhiTrig2, gPtTrig2, gChargeTrig2, 1.0));  \r
 842       if(fRunMixing) \r
 843         tracksMixing->Add(new AliBFBasicParticle(gEtaTrig2, gPhiTrig2, gPtTrig2, gChargeTrig2, 1.0));  \r
 844 \r
 845       iAssociated = 0; \r
 846       while(iAssociated < nAssociated) {\r
 847         gPtAssoc = fPtAssoc->GetRandom();\r
 848         if(gPtAssoc < gPtTrig2) {\r
 849           gEtaAssoc = gRandom->Uniform(gEtaTrig2 - gJetCone/2.,gEtaTrig2 + gJetCone/2.);\r
 850           gPhiAssoc = gRandom->Uniform(gPhiTrig2 - gJetCone/2.,gPhiTrig2 + gJetCone/2.);\r
 851           if(gPhiAssoc < 0.) gPhiAssoc += TMath::TwoPi();\r
 852           else if(gPhiAssoc > TMath::TwoPi()) gPhiAssoc -= TMath::TwoPi();\r
 853           \r
 854           iAssociated += 1;\r
 855 \r
 856           gDecideCharge = gRandom->Rndm();\r
 857           if(gDecideCharge <= 0.5)\r
 858             gChargeAssoc = 1;\r
 859           else \r
 860             gChargeAssoc = -1;\r
 861           \r
 862           //Acceptance\r
 863           if((gEtaAssoc < fEtaMin) || (gEtaAssoc > fEtaMax)) continue;\r
 864           //pt coverage\r
 865           if((gPtAssoc < fPtMin) || (gPtAssoc > fPtMax)) continue;\r
 866           //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);\r
 867 \r
 868           //acceptance filter\r
 869           if(fUseAcceptanceParameterization) {\r
 870             Double_t gRandomNumberForAcceptance = gRandom->Rndm();\r
 871             if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(gPtAssoc)) \r
 872               continue;\r
 873           }\r
 874 \r
 875           //Detector effects\r
 876           if(fSimulateDetectorEffects) {\r
 877             Double_t randomNumber = gRandom->Rndm();\r
 878             if(randomNumber > fEfficiencyMatrix->GetBinContent(fEfficiencyMatrix->FindBin(gEtaAssoc,gPtAssoc,gPhiAssoc)))\r
 879               continue;\r
 880           }\r
 881 \r
 882           gNumberOfAcceptedParticles += 1;\r
 883 \r
 884           // add the track to the TObjArray\r
 885           tracksMain->Add(new AliBFBasicParticle(gEtaAssoc, gPhiAssoc, gPtAssoc, gChargeAssoc, 1.0));  \r
 886           if(fRunMixing) \r
 887             tracksMixing->Add(new AliBFBasicParticle(gEtaAssoc, gPhiAssoc, gPtAssoc, gChargeAssoc, 1.0));  \r
 888         }//pt,assoc < pt,trig\r
 889       }//associated\r
 890     }//Jet usage\r
 891 \r
 892 \r
 893 \r
 894 \r
 895 \r
 896 \r
 897 \r
 898 \r
 899     \r
 900     //Dynamical correlations\r
 901     Int_t nGeneratedPositiveDynamicalCorrelations = 0;\r
 902     Int_t nGeneratedNegativeDynamicalCorrelations = 0;\r
 903     /*Double_t vChargePrime = 0;\r
 904     Double_t vYPrime = 0.0;\r
 905     Double_t vEtaPrime = 0.0;\r
 906     Double_t vPhiPrime = 0.0;\r
 907     Double_t vPPrime[3] = {0.,0.,0.};\r
 908     Double_t vPtPrime = 0.0;\r
 909     Double_t vEPrime = 0.0;\r
 910     //Generate "correlated" particles \r
 911     if(fUseDynamicalCorrelations) {\r
 912       Int_t gNumberOfDynamicalCorrelations = (Int_t)(0.5*gNumberOfAcceptedParticles*fDynamicalCorrelationsPercentage);\r
 913       for(Int_t iDynamicalCorrelations = 0; iDynamicalCorrelations < gNumberOfDynamicalCorrelations; iDynamicalCorrelations++) {\r
 914         isPion = kFALSE; isKaon = kFALSE; isProton = kFALSE;\r
 915         \r
 916         //Pseudo-rapidity sampled from a Gaussian centered @ 0\r
 917         vEta = gRandom->Gaus(0.0,0.1);\r
 918         vCharge = 1.0;\r
 919         nGeneratedPositiveDynamicalCorrelations += 1;\r
 920         \r
 921         vEtaPrime = -vEta;\r
 922         vChargePrime = -1.0;\r
 923         nGeneratedNegativeDynamicalCorrelations += 1;\r
 924           \r
 925         //Acceptance\r
 926         if((vEta < fEtaMin) || (vEta > fEtaMax)) continue;\r
 927         if((vEtaPrime < fEtaMin) || (vEtaPrime > fEtaMax)) continue;\r
 928       \r
 929         if(!fUseAllCharges) {\r
 930           //Decide the specie\r
 931           Double_t randomNumberSpecies = gRandom->Rndm();\r
 932           if((randomNumberSpecies >= 0.0)&&(randomNumberSpecies < fPionPercentage)) {\r
 933             nGeneratedPions += 1;\r
 934             vPt = fPtSpectraPions->GetRandom();\r
 935             vPhi = fAzimuthalAnglePions->GetRandom();\r
 936             fParticleMass = fPionMass;\r
 937             isPion = kTRUE;\r
 938           }\r
 939           else if((randomNumberSpecies >= fPionPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage)) {\r
 940             nGeneratedKaons += 1;\r
 941             vPt = fPtSpectraKaons->GetRandom();\r
 942             vPhi = fAzimuthalAngleKaons->GetRandom();\r
 943             fParticleMass = fKaonMass;\r
 944             isKaon = kTRUE;\r
 945           }\r
 946           else if((randomNumberSpecies >= fPionPercentage + fKaonPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage + fProtonPercentage)) {\r
 947             nGeneratedProtons += 1;\r
 948             vPt = fPtSpectraProtons->GetRandom();\r
 949             vPtPrime = vPt;\r
 950             vPhi = fAzimuthalAngleProtons->GetRandom();\r
 951             fParticleMass = fProtonMass;\r
 952             isProton = kTRUE;\r
 953           }\r
 954         }\r
 955         else {\r
 956           vPt = fPtSpectraAllCharges->GetRandom();\r
 957           vPhi = fAzimuthalAngleAllCharges->GetRandom();\r
 958         }\r
 959         vPtPrime = vPt;\r
 960         vPhiPrime = vPhi;\r
 961 \r
 962         vP[0] = vPt*TMath::Cos(vPhi);\r
 963         vP[1] = vPt*TMath::Sin(vPhi);\r
 964         vP[2] = vPt*TMath::SinH(vEta);\r
 965         vE = TMath::Sqrt(TMath::Power(fParticleMass,2) +\r
 966                           TMath::Power(vP[0],2) +\r
 967                           TMath::Power(vP[1],2) +\r
 968                           TMath::Power(vP[2],2));\r
 969         \r
 970         vY = 0.5*TMath::Log((vE + vP[2])/(vE - vP[2]));\r
 971 \r
 972         vPPrime[0] = vPtPrime*TMath::Cos(vPhiPrime);\r
 973         vPPrime[1] = vPtPrime*TMath::Sin(vPhiPrime);\r
 974         vPPrime[2] = vPtPrime*TMath::SinH(vEtaPrime);\r
 975         vEPrime = TMath::Sqrt(TMath::Power(fParticleMass,2) +\r
 976                           TMath::Power(vPPrime[0],2) +\r
 977                           TMath::Power(vPPrime[1],2) +\r
 978                           TMath::Power(vPPrime[2],2));\r
 979         \r
 980         vYPrime = 0.5*TMath::Log((vEPrime + vPPrime[2])/(vEPrime - vPPrime[2]));\r
 981       \r
 982         //pt coverage\r
 983         if((vPt < fPtMin) || (vPt > fPtMax)) continue;\r
 984         if((vPtPrime < fPtMin) || (vPtPrime > fPtMax)) continue;\r
 985 \r
 986         //acceptance filter\r
 987         if(fUseAcceptanceParameterization) {\r
 988           Double_t gRandomNumberForAcceptance = gRandom->Rndm();\r
 989           if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(vPt)) \r
 990             continue;\r
 991           \r
 992           Double_t gRandomNumberForAcceptancePrime = gRandom->Rndm();\r
 993           if(gRandomNumberForAcceptancePrime > fAcceptanceParameterization->Eval(vPtPrime)) \r
 994             continue;\r
 995         }\r
 996       \r
 997         // fill charge vector (positive)\r
 998         chargeVector[0]->push_back(vCharge);\r
 999         chargeVector[1]->push_back(vY);\r
1000         chargeVector[2]->push_back(vEta);\r
1001         chargeVector[3]->push_back(vPhi);\r
1002         chargeVector[4]->push_back(vP[0]);\r
1003         chargeVector[5]->push_back(vP[1]);\r
1004         chargeVector[6]->push_back(vP[2]);\r
1005         chargeVector[7]->push_back(vPt);\r
1006         chargeVector[8]->push_back(vE);\r
1007         \r
1008         if(fRunShuffling) {\r
1009           chargeVectorShuffle[0]->push_back(vCharge);\r
1010           chargeVectorShuffle[1]->push_back(vY);\r
1011           chargeVectorShuffle[2]->push_back(vEta);\r
1012           chargeVectorShuffle[3]->push_back(vPhi);\r
1013           chargeVectorShuffle[4]->push_back(vP[0]);\r
1014           chargeVectorShuffle[5]->push_back(vP[1]);\r
1015           chargeVectorShuffle[6]->push_back(vP[2]);\r
1016           chargeVectorShuffle[7]->push_back(vPt);\r
1017           chargeVectorShuffle[8]->push_back(vE);\r
1018         }\r
1019 \r
1020         // fill charge vector (negative)\r
1021         chargeVector[0]->push_back(vChargePrime);\r
1022         chargeVector[1]->push_back(vYPrime);\r
1023         chargeVector[2]->push_back(vEtaPrime);\r
1024         chargeVector[3]->push_back(vPhiPrime);\r
1025         chargeVector[4]->push_back(vPPrime[0]);\r
1026         chargeVector[5]->push_back(vPPrime[1]);\r
1027         chargeVector[6]->push_back(vPPrime[2]);\r
1028         chargeVector[7]->push_back(vPtPrime);\r
1029         chargeVector[8]->push_back(vEPrime);\r
1030         \r
1031         if(fRunShuffling) {\r
1032           chargeVectorShuffle[0]->push_back(vChargePrime);\r
1033           chargeVectorShuffle[1]->push_back(vYPrime);\r
1034           chargeVectorShuffle[2]->push_back(vEtaPrime);\r
1035           chargeVectorShuffle[3]->push_back(vPhiPrime);\r
1036           chargeVectorShuffle[4]->push_back(vPPrime[0]);\r
1037           chargeVectorShuffle[5]->push_back(vPPrime[1]);\r
1038           chargeVectorShuffle[6]->push_back(vPPrime[2]);\r
1039           chargeVectorShuffle[7]->push_back(vPtPrime);\r
1040           chargeVectorShuffle[8]->push_back(vEPrime);\r
1041         }\r
1042 \r
1043         gNumberOfAcceptedParticles += 2;\r
1044         }//loop over the dynamical correlations\r
1045         }*/// usage of dynamical correlations\r
1046 \r
1047     if(fUseDebug) {\r
1048       Printf("=======================================================");\r
1049       Printf("Total: %d - Total positive: %d - Total negative: %d",nMultiplicity,nGeneratedPositive,nGeneratedNegative);\r
1050       Printf("Accepted positive: %d - Accepted negative: %d",gNumberOfAcceptedPositiveParticles,gNumberOfAcceptedNegativeParticles);\r
1051       Printf("Correlations: %d - Correlations positive: %d - Correlations negative: %d",nGeneratedPositiveDynamicalCorrelations+nGeneratedNegativeDynamicalCorrelations,nGeneratedPositiveDynamicalCorrelations,nGeneratedNegativeDynamicalCorrelations);\r
1052       Printf("Number of accepted particles: %d",gNumberOfAcceptedParticles);\r
1053       if(!fUseAllCharges)\r
1054         Printf("Pions: %lf - Kaons: %lf - Protons: %lf",1.*nGeneratedPions/nMultiplicity,1.*nGeneratedKaons/nMultiplicity,1.*nGeneratedProtons/nMultiplicity);\r
1055       //Printf("Calculating the balance function for %d particles",chargeVector[0]->size());\r
1056     }\r
1057 \r
1058     fHistEventStats->Fill(4);\r
1059     fHistNumberOfAcceptedParticles->Fill(gNumberOfAcceptedParticles);\r
1060     fHistReactionPlane->Fill(fReactionPlane);\r
1061 \r
1062     // Event mixing \r
1063     if (fRunMixing)\r
1064       {\r
1065         // 1. First get an event pool corresponding in mult (cent) and\r
1066         //    zvertex to the current event. Once initialized, the pool\r
1067         //    should contain nMix (reduced) events. This routine does not\r
1068         //    pre-scan the chain. The first several events of every chain\r
1069         //    will be skipped until the needed pools are filled to the\r
1070         //    specified depth. If the pool categories are not too rare, this\r
1071         //    should not be a problem. If they are rare, you could lose`\r
1072         //    statistics.\r
1073         \r
1074         // 2. Collect the whole pool's content of tracks into one TObjArray\r
1075         //    (bgTracks), which is effectively a single background super-event.\r
1076         \r
1077         // 3. The reduced and bgTracks arrays must both be passed into\r
1078         //    FillCorrelations(). Also nMix should be passed in, so a weight\r
1079         //    of 1./nMix can be applied.\r
1080         \r
1081         AliEventPool* pool = fPoolMgr->GetEventPool(1., 0.,fReactionPlane);\r
1082       \r
1083         if (!pool){\r
1084           AliFatal(Form("No pool found for centrality = %f, zVtx = %f, psi = %f", 1., 0.,fReactionPlane));\r
1085         }\r
1086         else{\r
1087           \r
1088           //pool->SetDebug(1);\r
1089           \r
1090           if (pool->IsReady() || pool->NTracksInPool() > 50000 / 10 || pool->GetCurrentNEvents() >= 5){ \r
1091             \r
1092             \r
1093             Int_t nMix = pool->GetCurrentNEvents();\r
1094             //cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << ", tracks in this event = "<<gNumberOfAcceptedParticles <<", event Plane = "<<fReactionPlane<<endl;\r
1095             \r
1096             //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(2);\r
1097             //((TH2F*) fListOfHistos->FindObject("mixedDist"))->Fill(centrality, pool->NTracksInPool());\r
1098             //if (pool->IsReady())\r
1099             //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(3);\r
1100             \r
1101             // Fill mixed-event histos here  \r
1102             for (Int_t jMix=0; jMix<nMix; jMix++) \r
1103               {\r
1104                 TObjArray* tracksMixed = pool->GetEvent(jMix);\r
1105                 fMixedBalance->CalculateBalance(fReactionPlane,tracksMain,tracksMixed,1,1.,0.);\r
1106               }\r
1107           }\r
1108           \r
1109           // Update the Event pool\r
1110           pool->UpdatePool(tracksMain);\r
1111           //pool->PrintInfo();\r
1112           \r
1113         }//pool NULL check  \r
1114       }//run mixing\r
1115     \r
1116     //Calculate the balance function\r
1117     fBalance->CalculateBalance(fReactionPlane,tracksMain,NULL,1,1.,0.);\r
1118     //if(fRunMixing)\r
1119     //  fMixedBalance->CalculateBalance(fReactionPlane,tracksMixing,NULL,1,1.,0.);           \r
1120 \r
1121   }//event loop\r
1122 }      \r
1123 \r
1124 //________________________________________________________________________\r
1125 void  AliAnalysisTaskToyModel::FinishOutput() {\r
1126   //Printf("END BF");\r
1127 \r
1128   TFile *gOutput = TFile::Open("outputToyModel.root","recreate");\r
1129   TDirectoryFile *dir = new TDirectoryFile("PWGCFEbyE.outputBalanceFunctionPsiAnalysis","PWGCFEbyE.outputBalanceFunctionPsiAnalysis");\r
1130   \r
1131   fList->SetName("listQA");\r
1132   fList->SetOwner(kTRUE);\r
1133   dir->Add(fList);\r
1134   //fList->Write("listQA",TObject::kSingleKey);\r
1135 \r
1136   if (!fBalance) {\r
1137     AliError("ERROR: fBalance not available");\r
1138     return;\r
1139   }  \r
1140   fListBF->SetName("listBF");\r
1141   fListBF->SetOwner(kTRUE);\r
1142   dir->Add(fListBF);\r
1143   //fListBF->Write("listBF",TObject::kSingleKey);\r
1144 \r
1145   if(fRunShuffling) {\r
1146     if (!fShuffledBalance) {\r
1147       AliError("ERROR: fShuffledBalance not available");\r
1148       return;\r
1149     }\r
1150     fListBFS->SetName("listBFShuffled");\r
1151     fListBFS->SetOwner(kTRUE);\r
1152     dir->Add(fListBFS);\r
1153     //fListBFS->Write("listBFShuffled",TObject::kSingleKey);\r
1154   }\r
1155 \r
1156   if(fRunMixing) {\r
1157     if (!fMixedBalance) {\r
1158       AliError("ERROR: fMixedBalance not available");\r
1159       return;\r
1160     }\r
1161     fListBFM->SetName("listBFMixed");\r
1162     fListBFM->SetOwner(kTRUE);\r
1163     dir->Add(fListBFM);\r
1164     //fListBFM->Write("listBFMixed",TObject::kSingleKey);\r
1165   }\r
1166 \r
1167   dir->Write(dir->GetName(),TObject::kSingleKey);\r
1168   gOutput->Close();\r
1169 }\r
1170 \r