PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskBFPsi.cxx

   1 #include "TChain.h"\r
   2 #include "TList.h"\r
   3 #include "TCanvas.h"\r
   4 #include "TLorentzVector.h"\r
   5 #include "TGraphErrors.h"\r
   6 #include "TH1F.h"\r
   7 #include "TH2F.h"\r
   8 #include "TH2D.h"                  \r
   9 #include "TH3D.h"\r
  10 #include "TArrayF.h"\r
  11 #include "TF1.h"\r
  12 #include "TRandom.h"\r
  13 \r
  14 #include "AliAnalysisTaskSE.h"\r
  15 #include "AliAnalysisManager.h"\r
  16 \r
  17 #include "AliESDVertex.h"\r
  18 #include "AliESDEvent.h"\r
  19 #include "AliESDInputHandler.h"\r
  20 #include "AliAODEvent.h"\r
  21 #include "AliAODTrack.h"\r
  22 #include "AliAODInputHandler.h"\r
  23 #include "AliCollisionGeometry.h"\r
  24 #include "AliGenEventHeader.h"\r
  25 #include "AliMCEventHandler.h"\r
  26 #include "AliMCEvent.h"\r
  27 #include "AliStack.h"\r
  28 #include "AliESDtrackCuts.h"\r
  29 #include "AliEventplane.h"\r
  30 #include "AliTHn.h"    \r
  31 \r
  32 #include "AliEventPoolManager.h"           \r
  33 \r
  34 #include "AliPID.h"                \r
  35 #include "AliPIDResponse.h"        \r
  36 #include "AliPIDCombined.h"        \r
  37 \r
  38 #include "AliAnalysisTaskBFPsi.h"\r
  39 #include "AliBalancePsi.h"\r
  40 #include "AliAnalysisTaskTriggeredBF.h"\r
  41 \r
  42 \r
  43 // Analysis task for the BF vs Psi code\r
  44 // Authors: Panos.Christakoglou@nikhef.nl\r
  45 \r
  46 using std::cout;\r
  47 using std::endl;\r
  48 \r
  49 ClassImp(AliAnalysisTaskBFPsi)\r
  50 \r
  51 //________________________________________________________________________\r
  52 AliAnalysisTaskBFPsi::AliAnalysisTaskBFPsi(const char *name) \r
  53 : AliAnalysisTaskSE(name), \r
  54   fBalance(0),\r
  55   fRunShuffling(kFALSE),\r
  56   fShuffledBalance(0),\r
  57   fRunMixing(kFALSE),\r
  58   fRunMixingEventPlane(kFALSE),\r
  59   fMixingTracks(50000),\r
  60   fMixedBalance(0),\r
  61   fPoolMgr(0),\r
  62   fList(0),\r
  63   fListBF(0),\r
  64   fListBFS(0),\r
  65   fListBFM(0),\r
  66   fHistListPIDQA(0),\r
  67   fHistEventStats(0),\r
  68   fHistCentStats(0),\r
  69   fHistTriggerStats(0),\r
  70   fHistTrackStats(0),\r
  71   fHistVx(0),\r
  72   fHistVy(0),\r
  73   fHistVz(0),\r
  74   fHistEventPlane(0),\r
  75   fHistClus(0),\r
  76   fHistDCA(0),\r
  77   fHistChi2(0),\r
  78   fHistPt(0),\r
  79   fHistEta(0),\r
  80   fHistRapidity(0),\r
  81   fHistPhi(0),\r
  82   fHistPhiBefore(0),\r
  83   fHistPhiAfter(0),\r
  84   fHistPhiPos(0),\r
  85   fHistPhiNeg(0),\r
  86   fHistV0M(0),\r
  87   fHistRefTracks(0),\r
  88   fHistdEdxVsPTPCbeforePID(NULL),\r
  89   fHistBetavsPTOFbeforePID(NULL), \r
  90   fHistProbTPCvsPtbeforePID(NULL), \r
  91   fHistProbTOFvsPtbeforePID(NULL), \r
  92   fHistProbTPCTOFvsPtbeforePID(NULL),\r
  93   fHistNSigmaTPCvsPtbeforePID(NULL), \r
  94   fHistNSigmaTOFvsPtbeforePID(NULL), \r
  95   fHistdEdxVsPTPCafterPID(NULL),\r
  96   fHistBetavsPTOFafterPID(NULL), \r
  97   fHistProbTPCvsPtafterPID(NULL), \r
  98   fHistProbTOFvsPtafterPID(NULL), \r
  99   fHistProbTPCTOFvsPtafterPID(NULL),\r
 100   fHistNSigmaTPCvsPtafterPID(NULL), \r
 101   fHistNSigmaTOFvsPtafterPID(NULL),  \r
 102   fPIDResponse(0x0),\r
 103   fPIDCombined(0x0),\r
 104   fParticleOfInterest(kPion),\r
 105   fPidDetectorConfig(kTPCTOF),\r
 106   fUsePID(kFALSE),\r
 107   fUsePIDnSigma(kTRUE),\r
 108   fUsePIDPropabilities(kFALSE), \r
 109   fPIDNSigma(3.),\r
 110   fMinAcceptedPIDProbability(0.8),\r
 111   fESDtrackCuts(0),\r
 112   fCentralityEstimator("V0M"),\r
 113   fUseCentrality(kFALSE),\r
 114   fCentralityPercentileMin(0.), \r
 115   fCentralityPercentileMax(5.),\r
 116   fImpactParameterMin(0.),\r
 117   fImpactParameterMax(20.),\r
 118   fUseMultiplicity(kFALSE),\r
 119   fNumberOfAcceptedTracksMin(0),\r
 120   fNumberOfAcceptedTracksMax(10000),\r
 121   fHistNumberOfAcceptedTracks(0),\r
 122   fUseOfflineTrigger(kFALSE),\r
 123   fVxMax(0.3),\r
 124   fVyMax(0.3),\r
 125   fVzMax(10.),\r
 126   nAODtrackCutBit(128),\r
 127   fPtMin(0.3),\r
 128   fPtMax(1.5),\r
 129   fPtMinForCorrections(0.3),//=================================correction\r
 130   fPtMaxForCorrections(1.5),//=================================correction\r
 131   fPtBinForCorrections(36), //=================================correction\r
 132   fEtaMin(-0.8),\r
 133   fEtaMax(-0.8),\r
 134   fEtaMinForCorrections(-0.8),//=================================correction\r
 135   fEtaMaxForCorrections(0.8),//=================================correction\r
 136   fEtaBinForCorrections(16), //=================================correction\r
 137   fPhiMin(0.),\r
 138   fPhiMax(360.),\r
 139   fPhiMinForCorrections(0.),//=================================correction\r
 140   fPhiMaxForCorrections(360.),//=================================correction\r
 141   fPhiBinForCorrections(100), //=================================correction\r
 142   fDCAxyCut(-1),\r
 143   fDCAzCut(-1),\r
 144   fTPCchi2Cut(-1),\r
 145   fNClustersTPCCut(-1),\r
 146   fAcceptanceParameterization(0),\r
 147   fDifferentialV2(0),\r
 148   fUseFlowAfterBurner(kFALSE),\r
 149   fExcludeResonancesInMC(kFALSE),\r
 150   fUseMCPdgCode(kFALSE),\r
 151   fPDGCodeToBeAnalyzed(-1),\r
 152   fEventClass("EventPlane") \r
 153 {\r
 154   // Constructor\r
 155   // Define input and output slots here\r
 156   // Input slot #0 works with a TChain\r
 157 \r
 158   //======================================================correction\r
 159   for (Int_t i=0; i<kCENTRALITY; i++){\r
 160     fHistMatrixCorrectionPlus[i] = NULL; \r
 161     fHistMatrixCorrectionMinus[i] = NULL; \r
 162   }\r
 163   //=====================================================correction\r
 164 \r
 165   DefineInput(0, TChain::Class());\r
 166   // Output slot #0 writes into a TH1 container\r
 167   DefineOutput(1, TList::Class());\r
 168   DefineOutput(2, TList::Class());\r
 169   DefineOutput(3, TList::Class());\r
 170   DefineOutput(4, TList::Class());\r
 171   DefineOutput(5, TList::Class());\r
 172 }\r
 173 \r
 174 //________________________________________________________________________\r
 175 AliAnalysisTaskBFPsi::~AliAnalysisTaskBFPsi() {\r
 176 \r
 177   // delete fBalance; \r
 178   // delete fShuffledBalance; \r
 179   // delete fList;\r
 180   // delete fListBF; \r
 181   // delete fListBFS;\r
 182 \r
 183   // delete fHistEventStats; \r
 184   // delete fHistTrackStats; \r
 185   // delete fHistVx; \r
 186   // delete fHistVy; \r
 187   // delete fHistVz; \r
 188 \r
 189   // delete fHistClus;\r
 190   // delete fHistDCA;\r
 191   // delete fHistChi2;\r
 192   // delete fHistPt;\r
 193   // delete fHistEta;\r
 194   // delete fHistPhi;\r
 195   // delete fHistV0M;\r
 196 }\r
 197 \r
 198 //________________________________________________________________________\r
 199 void AliAnalysisTaskBFPsi::UserCreateOutputObjects() {\r
 200   // Create histograms\r
 201   // Called once\r
 202 \r
 203   // global switch disabling the reference \r
 204   // (to avoid "Replacing existing TH1" if several wagons are created in train)\r
 205   Bool_t oldStatus = TH1::AddDirectoryStatus();\r
 206   TH1::AddDirectory(kFALSE);\r
 207 \r
 208   if(!fBalance) {\r
 209     fBalance = new AliBalancePsi();\r
 210     fBalance->SetAnalysisLevel("ESD");\r
 211     fBalance->SetEventClass(fEventClass);\r
 212     //fBalance->SetNumberOfBins(-1,16);\r
 213     //fBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);\r
 214   }\r
 215   if(fRunShuffling) {\r
 216     if(!fShuffledBalance) {\r
 217       fShuffledBalance = new AliBalancePsi();\r
 218       fShuffledBalance->SetAnalysisLevel("ESD");\r
 219       fShuffledBalance->SetEventClass(fEventClass);\r
 220       //fShuffledBalance->SetNumberOfBins(-1,16);\r
 221       //fShuffledBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);\r
 222     }\r
 223   }\r
 224   if(fRunMixing) {\r
 225     if(!fMixedBalance) {\r
 226       fMixedBalance = new AliBalancePsi();\r
 227       fMixedBalance->SetAnalysisLevel("ESD");\r
 228       fMixedBalance->SetEventClass(fEventClass);\r
 229     }\r
 230   }\r
 231 \r
 232   //QA list\r
 233   fList = new TList();\r
 234   fList->SetName("listQA");\r
 235   fList->SetOwner();\r
 236 \r
 237   //Balance Function list\r
 238   fListBF = new TList();\r
 239   fListBF->SetName("listBF");\r
 240   fListBF->SetOwner();\r
 241 \r
 242   if(fRunShuffling) {\r
 243     fListBFS = new TList();\r
 244     fListBFS->SetName("listBFShuffled");\r
 245     fListBFS->SetOwner();\r
 246   }\r
 247 \r
 248   if(fRunMixing) {\r
 249     fListBFM = new TList();\r
 250     fListBFM->SetName("listTriggeredBFMixed");\r
 251     fListBFM->SetOwner();\r
 252   }\r
 253 \r
 254   //PID QA list\r
 255   if(fUsePID) {\r
 256     fHistListPIDQA = new TList();\r
 257     fHistListPIDQA->SetName("listQAPID");\r
 258     fHistListPIDQA->SetOwner();\r
 259   }\r
 260 \r
 261   //Event stats.\r
 262   TString gCutName[5] = {"Total","Offline trigger",\r
 263                          "Vertex","Analyzed","sel. Centrality"};\r
 264   fHistEventStats = new TH2F("fHistEventStats",\r
 265                              "Event statistics;;Centrality percentile;N_{events}",\r
 266                              5,0.5,5.5,220,-5,105);\r
 267   for(Int_t i = 1; i <= 5; i++)\r
 268     fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());\r
 269   fList->Add(fHistEventStats);\r
 270 \r
 271   TString gCentName[9] = {"V0M","FMD","TRK","TKL","CL0","CL1","V0MvsFMD","TKLvsV0M","ZEMvsZDC"};\r
 272   fHistCentStats = new TH2F("fHistCentStats",\r
 273                              "Centrality statistics;;Cent percentile",\r
 274                             9,-0.5,8.5,220,-5,105);\r
 275   for(Int_t i = 1; i <= 9; i++)\r
 276     fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());\r
 277   fList->Add(fHistCentStats);\r
 278 \r
 279   fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",1025,0,1025);\r
 280   fList->Add(fHistTriggerStats);\r
 281 \r
 282   fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TrackFilterBit;N_{events}",16,0,16);\r
 283   fList->Add(fHistTrackStats);\r
 284 \r
 285   fHistNumberOfAcceptedTracks = new TH2F("fHistNumberOfAcceptedTracks",";N_{acc.};Centrality percentile;Entries",4001,-0.5,4000.5,220,-5,105);\r
 286   fList->Add(fHistNumberOfAcceptedTracks);\r
 287 \r
 288   // Vertex distributions\r
 289   fHistVx = new TH1F("fHistVx","Primary vertex distribution - x coordinate;V_{x} (cm);Entries",100,-0.5,0.5);\r
 290   fList->Add(fHistVx);\r
 291   fHistVy = new TH1F("fHistVy","Primary vertex distribution - y coordinate;V_{y} (cm);Entries",100,-0.5,0.5);\r
 292   fList->Add(fHistVy);\r
 293   fHistVz = new TH2F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Centrality percentile;Entries",100,-20.,20.,220,-5,105);\r
 294   fList->Add(fHistVz);\r
 295 \r
 296   //Event plane\r
 297   fHistEventPlane = new TH2F("fHistEventPlane",";#Psi_{2} [deg.];Centrality percentile;Counts",100,0,360.,220,-5,105);\r
 298   fList->Add(fHistEventPlane);\r
 299 \r
 300   // QA histograms\r
 301   fHistClus = new TH2F("fHistClus","# Cluster (TPC vs. ITS)",10,0,10,200,0,200);\r
 302   fList->Add(fHistClus);\r
 303   fHistChi2 = new TH2F("fHistChi2","Chi2/NDF distribution;#chi^{2}/ndf;Centrality percentile",200,0,10,220,-5,105);\r
 304   fList->Add(fHistChi2);\r
 305   fHistDCA  = new TH2F("fHistDCA","DCA (xy vs. z)",400,-5,5,400,-5,5); \r
 306   fList->Add(fHistDCA);\r
 307   fHistPt   = new TH2F("fHistPt","p_{T} distribution;p_{T} (GeV/c);Centrality percentile",200,0,10,220,-5,105);\r
 308   fList->Add(fHistPt);\r
 309   fHistEta  = new TH2F("fHistEta","#eta distribution;#eta;Centrality percentile",200,-2,2,220,-5,105);\r
 310   fList->Add(fHistEta);\r
 311   fHistRapidity  = new TH2F("fHistRapidity","y distribution;y;Centrality percentile",200,-2,2,220,-5,105);\r
 312   fList->Add(fHistRapidity);\r
 313   fHistPhi  = new TH2F("fHistPhi","#phi distribution;#phi (rad);Centrality percentile",200,0.0,2.*TMath::Pi(),220,-5,105);\r
 314   fList->Add(fHistPhi);\r
 315   fHistPhiBefore  = new TH2F("fHistPhiBefore","#phi distribution;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);\r
 316   fList->Add(fHistPhiBefore);\r
 317   fHistPhiAfter  = new TH2F("fHistPhiAfter","#phi distribution;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);\r
 318   fList->Add(fHistPhiAfter);\r
 319   fHistPhiPos  = new TH2F("fHistPhiPos","#phi distribution for positive particles;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);\r
 320   fList->Add(fHistPhiPos);\r
 321   fHistPhiNeg  = new TH2F("fHistPhiNeg","#phi distribution for negative particles;#phi;Centrality percentile",200,0.,2.*TMath::Pi(),220,-5,105);\r
 322   fList->Add(fHistPhiNeg);\r
 323   fHistV0M  = new TH2F("fHistV0M","V0 Multiplicity C vs. A",500, 0, 20000, 500, 0, 20000);\r
 324   fList->Add(fHistV0M);\r
 325   TString gRefTrackName[6] = {"tracks","tracksPos","tracksNeg","tracksTPConly","clusITS0","clusITS1"};\r
 326   fHistRefTracks  = new TH2F("fHistRefTracks","Nr of Ref tracks/event vs. ref track estimator;;Nr of tracks",6, 0, 6, 400, 0, 20000);\r
 327   for(Int_t i = 1; i <= 6; i++)\r
 328     fHistRefTracks->GetXaxis()->SetBinLabel(i,gRefTrackName[i-1].Data());\r
 329   fList->Add(fHistRefTracks);\r
 330 \r
 331   // QA histograms for HBTinspired and Conversion cuts\r
 332   fList->Add(fBalance->GetQAHistHBTbefore());\r
 333   fList->Add(fBalance->GetQAHistHBTafter());\r
 334   fList->Add(fBalance->GetQAHistConversionbefore());\r
 335   fList->Add(fBalance->GetQAHistConversionafter());\r
 336   fList->Add(fBalance->GetQAHistPsiMinusPhi());\r
 337 \r
 338   // Balance function histograms\r
 339   // Initialize histograms if not done yet\r
 340   if(!fBalance->GetHistNp()){\r
 341     AliWarning("Histograms not yet initialized! --> Will be done now");\r
 342     AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
 343     fBalance->InitHistograms();\r
 344   }\r
 345 \r
 346   if(fRunShuffling) {\r
 347     if(!fShuffledBalance->GetHistNp()) {\r
 348       AliWarning("Histograms (shuffling) not yet initialized! --> Will be done now");\r
 349       AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
 350       fShuffledBalance->InitHistograms();\r
 351     }\r
 352   }\r
 353 \r
 354   //for(Int_t a = 0; a < ANALYSIS_TYPES; a++){\r
 355   fListBF->Add(fBalance->GetHistNp());\r
 356   fListBF->Add(fBalance->GetHistNn());\r
 357   fListBF->Add(fBalance->GetHistNpn());\r
 358   fListBF->Add(fBalance->GetHistNnn());\r
 359   fListBF->Add(fBalance->GetHistNpp());\r
 360   fListBF->Add(fBalance->GetHistNnp());\r
 361 \r
 362   if(fRunShuffling) {\r
 363     fListBFS->Add(fShuffledBalance->GetHistNp());\r
 364     fListBFS->Add(fShuffledBalance->GetHistNn());\r
 365     fListBFS->Add(fShuffledBalance->GetHistNpn());\r
 366     fListBFS->Add(fShuffledBalance->GetHistNnn());\r
 367     fListBFS->Add(fShuffledBalance->GetHistNpp());\r
 368     fListBFS->Add(fShuffledBalance->GetHistNnp());\r
 369   }  \r
 370 \r
 371   if(fRunMixing) {\r
 372     fListBFM->Add(fMixedBalance->GetHistNp());\r
 373     fListBFM->Add(fMixedBalance->GetHistNn());\r
 374     fListBFM->Add(fMixedBalance->GetHistNpn());\r
 375     fListBFM->Add(fMixedBalance->GetHistNnn());\r
 376     fListBFM->Add(fMixedBalance->GetHistNpp());\r
 377     fListBFM->Add(fMixedBalance->GetHistNnp());\r
 378   }\r
 379   //}\r
 380 \r
 381 \r
 382   // Event Mixing\r
 383   if(fRunMixing){\r
 384     Int_t trackDepth = fMixingTracks; \r
 385     Int_t poolsize   = 1000;  // Maximum number of events, ignored in the present implemented of AliEventPoolManager\r
 386     \r
 387     // centrality bins\r
 388     Double_t centralityBins[] = {0.,1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,15.,20.,25.,30.,35.,40.,45.,50.,55.,60.,65.,70.,75.,80.,90.,100.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
 389     Double_t* centbins        = centralityBins;\r
 390     Int_t nCentralityBins     = sizeof(centralityBins) / sizeof(Double_t) - 1;\r
 391     \r
 392     // Zvtx bins\r
 393     Double_t vertexBins[] = {-10., -7., -5., -3., -1., 1., 3., 5., 7., 10.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
 394     Double_t* vtxbins     = vertexBins;\r
 395     Int_t nVertexBins     = sizeof(vertexBins) / sizeof(Double_t) - 1;\r
 396     \r
 397     // Event plane angle (Psi) bins\r
 398     Double_t psiBins[] = {0.,45.,135.,215.,305.,360.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
 399     Double_t* psibins     = psiBins;\r
 400     Int_t nPsiBins     = sizeof(psiBins) / sizeof(Double_t) - 1;\r
 401     \r
 402     // run the event mixing also in bins of event plane (statistics!)\r
 403     if(fRunMixingEventPlane){\r
 404       fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins, nPsiBins, psibins);\r
 405     }\r
 406     else{\r
 407       fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);\r
 408     }\r
 409   }\r
 410 \r
 411   if(fESDtrackCuts) fList->Add(fESDtrackCuts);\r
 412 \r
 413   //====================PID========================//\r
 414   if(fUsePID) {\r
 415     fPIDCombined = new AliPIDCombined();\r
 416     fPIDCombined->SetDefaultTPCPriors();\r
 417 \r
 418     fHistdEdxVsPTPCbeforePID = new TH2D ("dEdxVsPTPCbefore","dEdxVsPTPCbefore", 1000, -10.0, 10.0, 1000, 0, 1000); \r
 419     fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePID); //addition \r
 420     \r
 421     fHistBetavsPTOFbeforePID = new TH2D ("BetavsPTOFbefore","BetavsPTOFbefore", 1000, -10.0, 10., 1000, 0, 1.2); \r
 422     fHistListPIDQA->Add(fHistBetavsPTOFbeforePID); //addition\r
 423     \r
 424     fHistProbTPCvsPtbeforePID = new TH2D ("ProbTPCvsPtbefore","ProbTPCvsPtbefore", 1000, -10.0,10.0, 1000, 0, 2.0); \r
 425     fHistListPIDQA->Add(fHistProbTPCvsPtbeforePID); //addition \r
 426     \r
 427     fHistProbTOFvsPtbeforePID = new TH2D ("ProbTOFvsPtbefore","ProbTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0); \r
 428     fHistListPIDQA->Add(fHistProbTOFvsPtbeforePID); //addition \r
 429 \r
 430     fHistProbTPCTOFvsPtbeforePID =new TH2D ("ProbTPCTOFvsPtbefore","ProbTPCTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0); \r
 431     fHistListPIDQA->Add(fHistProbTPCTOFvsPtbeforePID); //addition \r
 432     \r
 433     fHistNSigmaTPCvsPtbeforePID = new TH2D ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore", 1000, -10, 10, 1000, 0, 500); \r
 434     fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePID); //addition \r
 435     \r
 436     fHistNSigmaTOFvsPtbeforePID = new TH2D ("NSigmaTOFvsPtbefore","NSigmaTOFvsPtbefore", 1000, -10, 10, 1000, 0, 500); \r
 437     fHistListPIDQA->Add(fHistNSigmaTOFvsPtbeforePID); //addition \r
 438     \r
 439     fHistdEdxVsPTPCafterPID = new TH2D ("dEdxVsPTPCafter","dEdxVsPTPCafter", 1000, -10, 10, 1000, 0, 1000); \r
 440     fHistListPIDQA->Add(fHistdEdxVsPTPCafterPID); //addition \r
 441     \r
 442     fHistBetavsPTOFafterPID = new TH2D ("BetavsPTOFafter","BetavsPTOFafter", 1000, -10, 10, 1000, 0, 1.2); \r
 443     fHistListPIDQA->Add(fHistBetavsPTOFafterPID); //addition \r
 444     \r
 445     fHistProbTPCvsPtafterPID = new TH2D ("ProbTPCvsPtafter","ProbTPCvsPtafter", 1000, -10, 10, 1000, 0, 2); \r
 446     fHistListPIDQA->Add(fHistProbTPCvsPtafterPID); //addition \r
 447   \r
 448     fHistProbTOFvsPtafterPID = new TH2D ("ProbTOFvsPtafter","ProbTOFvsPtafter", 1000,  -10, 10, 1000, 0, 2); \r
 449     fHistListPIDQA->Add(fHistProbTOFvsPtafterPID); //addition  \r
 450     \r
 451     fHistProbTPCTOFvsPtafterPID =new TH2D ("ProbTPCTOFvsPtafter","ProbTPCTOFvsPtafter", 1000, -50, 50, 1000, 0, 2.0); \r
 452     fHistListPIDQA->Add(fHistProbTPCTOFvsPtafterPID); //addition \r
 453 \r
 454     fHistNSigmaTPCvsPtafterPID = new TH2D ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter", 1000, -10, 10, 1000, 0, 500); \r
 455     fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPID); //addition  \r
 456     \r
 457     fHistNSigmaTOFvsPtafterPID = new TH2D ("NSigmaTOFvsPtafter","NSigmaTOFvsPtafter", 1000, -10, 10, 1000, 0, 500); \r
 458     fHistListPIDQA->Add(fHistNSigmaTOFvsPtafterPID); //addition \r
 459   }\r
 460   //====================PID========================//\r
 461 \r
 462   // Post output data.\r
 463   PostData(1, fList);\r
 464   PostData(2, fListBF);\r
 465   if(fRunShuffling) PostData(3, fListBFS);\r
 466   if(fRunMixing) PostData(4, fListBFM);\r
 467   if(fUsePID) PostData(5, fHistListPIDQA);       //PID\r
 468 \r
 469   TH1::AddDirectory(oldStatus);\r
 470 }\r
 471 \r
 472 \r
 473 //________________________________________________________________________\r
 474 void AliAnalysisTaskBFPsi::SetInputCorrection(const char* filename, \r
 475                                               const char* gCollSystem) {\r
 476   //Open files that will be used for correction\r
 477   TString gCollidingSystem = gCollSystem;\r
 478 \r
 479   //Open the input file\r
 480   TFile *f = TFile::Open(filename);\r
 481   if(!f->IsOpen()) {\r
 482     Printf("File not found!!!");\r
 483     return;\r
 484   }\r
 485   \r
 486   //Get the TDirectoryFile and TList\r
 487   TDirectoryFile *dir = dynamic_cast<TDirectoryFile *>(f->Get("PWGCFEbyE.outputBalanceFunctionEfficiencyAnalysis"));\r
 488   if(!dir) {\r
 489     Printf("TDirectoryFile not found!!!");\r
 490     return;\r
 491   }\r
 492   \r
 493   TString listEffName = "";\r
 494   for (Int_t iCent = 0; iCent < kCENTRALITY; iCent++) {\r
 495     listEffName = "listEfficiencyBF_";  \r
 496     listEffName += (TString)((Int_t)(centralityArrayForPbPb[iCent]));\r
 497     listEffName += "_";\r
 498     listEffName += (TString)((Int_t)(centralityArrayForPbPb[iCent+1]));\r
 499     TList *list = (TList*)dir->Get(listEffName.Data());\r
 500     if(!list) {\r
 501       Printf("TList Efficiency not found!!!");\r
 502       return;\r
 503     }\r
 504     \r
 505     TString histoName = "fHistMatrixCorrectionPlus";\r
 506     fHistMatrixCorrectionPlus[iCent]= dynamic_cast<TH3D *>(list->FindObject(histoName.Data()));\r
 507     if(!fHistMatrixCorrectionPlus[iCent]) {\r
 508       Printf("fHist not found!!!");\r
 509       return;\r
 510     }\r
 511     \r
 512     histoName = "fHistMatrixCorrectionMinus";\r
 513     fHistMatrixCorrectionMinus[iCent] = dynamic_cast<TH3D *>(list->FindObject(histoName.Data()));  \r
 514     if(!fHistMatrixCorrectionMinus[iCent]) {\r
 515       Printf("fHist not found!!!");\r
 516       return; \r
 517     }\r
 518   }//loop over centralities: ONLY the PbPb case is covered\r
 519 \r
 520   if(fHistMatrixCorrectionPlus[0]){\r
 521     fEtaMinForCorrections = fHistMatrixCorrectionPlus[0]->GetXaxis()->GetXmin();\r
 522     fEtaMaxForCorrections = fHistMatrixCorrectionPlus[0]->GetXaxis()->GetXmax();\r
 523     fEtaBinForCorrections = fHistMatrixCorrectionPlus[0]->GetNbinsX();\r
 524     \r
 525     fPtMinForCorrections = fHistMatrixCorrectionPlus[0]->GetYaxis()->GetXmin();\r
 526     fPtMaxForCorrections = fHistMatrixCorrectionPlus[0]->GetYaxis()->GetXmax();\r
 527     fPtBinForCorrections = fHistMatrixCorrectionPlus[0]->GetNbinsY();\r
 528     \r
 529     fPhiMinForCorrections = fHistMatrixCorrectionPlus[0]->GetZaxis()->GetXmin();\r
 530     fPhiMaxForCorrections = fHistMatrixCorrectionPlus[0]->GetZaxis()->GetXmax();\r
 531     fPhiBinForCorrections = fHistMatrixCorrectionPlus[0]->GetNbinsZ();\r
 532   }\r
 533 }\r
 534 \r
 535 //________________________________________________________________________\r
 536 void AliAnalysisTaskBFPsi::UserExec(Option_t *) {\r
 537   // Main loop\r
 538   // Called for each event\r
 539 \r
 540   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
 541   Int_t gNumberOfAcceptedTracks = 0;\r
 542   Double_t gCentrality          = -1.;\r
 543   Double_t gReactionPlane       = -1.; \r
 544   Float_t bSign = 0.;\r
 545   \r
 546   // get the event (for generator level: MCEvent())\r
 547   AliVEvent* eventMain = NULL;\r
 548   if(gAnalysisLevel == "MC") {\r
 549     eventMain = dynamic_cast<AliVEvent*>(MCEvent()); \r
 550   }\r
 551   else{\r
 552     eventMain = dynamic_cast<AliVEvent*>(InputEvent()); \r
 553     \r
 554     // for HBT like cuts need magnetic field sign\r
 555     bSign = (eventMain->GetMagneticField() > 0) ? 1 : -1;\r
 556   }\r
 557   if(!eventMain) {\r
 558     AliError("eventMain not available");\r
 559     return;\r
 560   }\r
 561   \r
 562   // PID Response task active?\r
 563   if(fUsePID) {\r
 564     fPIDResponse = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();\r
 565     if (!fPIDResponse) AliFatal("This Task needs the PID response attached to the inputHandler");\r
 566   }\r
 567   \r
 568   // check event cuts and fill event histograms\r
 569   if((gCentrality = IsEventAccepted(eventMain)) < 0){\r
 570     return;\r
 571   }\r
 572   \r
 573   //Compute Multiplicity or Centrality variable\r
 574   Double_t lMultiplicityVar = GetRefMultiOrCentrality( eventMain );\r
 575 \r
 576   // get the reaction plane\r
 577   gReactionPlane = GetEventPlane(eventMain);\r
 578   fHistEventPlane->Fill(gReactionPlane,gCentrality);\r
 579   if(gReactionPlane < 0){\r
 580     return;\r
 581   }\r
 582   \r
 583   // get the accepted tracks in main event\r
 584   TObjArray *tracksMain = GetAcceptedTracks(eventMain,gCentrality,gReactionPlane);\r
 585   gNumberOfAcceptedTracks = tracksMain->GetEntriesFast();\r
 586 \r
 587   //multiplicity cut (used in pp)\r
 588   fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedTracks,gCentrality);\r
 589   if(fUseMultiplicity) {\r
 590     if((gNumberOfAcceptedTracks < fNumberOfAcceptedTracksMin)||(gNumberOfAcceptedTracks > fNumberOfAcceptedTracksMax))\r
 591       return;\r
 592   }\r
 593 \r
 594   // store charges of all accepted tracks, shuffle and reassign (two extra loops!)\r
 595   TObjArray* tracksShuffled = NULL;\r
 596   if(fRunShuffling){\r
 597     tracksShuffled = GetShuffledTracks(tracksMain,gCentrality);\r
 598   }\r
 599   \r
 600   // Event mixing \r
 601   if (fRunMixing)\r
 602     {\r
 603       // 1. First get an event pool corresponding in mult (cent) and\r
 604       //    zvertex to the current event. Once initialized, the pool\r
 605       //    should contain nMix (reduced) events. This routine does not\r
 606       //    pre-scan the chain. The first several events of every chain\r
 607       //    will be skipped until the needed pools are filled to the\r
 608       //    specified depth. If the pool categories are not too rare, this\r
 609       //    should not be a problem. If they are rare, you could lose`\r
 610       //    statistics.\r
 611       \r
 612       // 2. Collect the whole pool's content of tracks into one TObjArray\r
 613       //    (bgTracks), which is effectively a single background super-event.\r
 614       \r
 615       // 3. The reduced and bgTracks arrays must both be passed into\r
 616       //    FillCorrelations(). Also nMix should be passed in, so a weight\r
 617       //    of 1./nMix can be applied.\r
 618       \r
 619       AliEventPool* pool = fPoolMgr->GetEventPool(gCentrality, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane);\r
 620       \r
 621       if (!pool){\r
 622         AliFatal(Form("No pool found for centrality = %f, zVtx = %f, psi = %f", gCentrality, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane));\r
 623       }\r
 624       else{\r
 625         \r
 626         //pool->SetDebug(1);\r
 627 \r
 628         if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 || pool->GetCurrentNEvents() >= 5){ \r
 629           \r
 630           \r
 631           Int_t nMix = pool->GetCurrentNEvents();\r
 632           //cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << endl;\r
 633           \r
 634           //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(2);\r
 635           //((TH2F*) fListOfHistos->FindObject("mixedDist"))->Fill(centrality, pool->NTracksInPool());\r
 636           //if (pool->IsReady())\r
 637           //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(3);\r
 638           \r
 639           // Fill mixed-event histos here  \r
 640           for (Int_t jMix=0; jMix<nMix; jMix++) \r
 641             {\r
 642               TObjArray* tracksMixed = pool->GetEvent(jMix);\r
 643               fMixedBalance->CalculateBalance(gReactionPlane,tracksMain,tracksMixed,bSign,lMultiplicityVar);\r
 644             }\r
 645         }\r
 646         \r
 647         // Update the Event pool\r
 648         pool->UpdatePool(tracksMain);\r
 649         //pool->PrintInfo();\r
 650         \r
 651       }//pool NULL check  \r
 652     }//run mixing\r
 653   \r
 654   // calculate balance function\r
 655   fBalance->CalculateBalance(gReactionPlane,tracksMain,NULL,bSign,lMultiplicityVar);\r
 656   \r
 657   // calculate shuffled balance function\r
 658   if(fRunShuffling && tracksShuffled != NULL) {\r
 659     fShuffledBalance->CalculateBalance(gReactionPlane,tracksShuffled,NULL,bSign,lMultiplicityVar);\r
 660   }\r
 661 }      \r
 662 \r
 663 //________________________________________________________________________\r
 664 Double_t AliAnalysisTaskBFPsi::IsEventAccepted(AliVEvent *event){\r
 665   // Checks the Event cuts\r
 666   // Fills Event statistics histograms\r
 667   \r
 668   Bool_t isSelectedMain = kTRUE;\r
 669   Float_t gCentrality = -1.;\r
 670   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
 671 \r
 672   fHistEventStats->Fill(1,gCentrality); //all events\r
 673 \r
 674   // Event trigger bits\r
 675   fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());\r
 676   if(fUseOfflineTrigger)\r
 677     isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
 678   \r
 679   if(isSelectedMain) {\r
 680     fHistEventStats->Fill(2,gCentrality); //triggered events\r
 681     \r
 682     //Centrality stuff \r
 683     if(fUseCentrality) {\r
 684       if(gAnalysisLevel == "AOD") { //centrality in AOD header\r
 685         AliAODHeader *header = (AliAODHeader*) event->GetHeader();\r
 686         if(header){\r
 687           gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
 688           \r
 689           // QA for centrality estimators\r
 690           fHistCentStats->Fill(0.,header->GetCentralityP()->GetCentralityPercentile("V0M"));\r
 691           fHistCentStats->Fill(1.,header->GetCentralityP()->GetCentralityPercentile("FMD"));\r
 692           fHistCentStats->Fill(2.,header->GetCentralityP()->GetCentralityPercentile("TRK"));\r
 693           fHistCentStats->Fill(3.,header->GetCentralityP()->GetCentralityPercentile("TKL"));\r
 694           fHistCentStats->Fill(4.,header->GetCentralityP()->GetCentralityPercentile("CL0"));\r
 695           fHistCentStats->Fill(5.,header->GetCentralityP()->GetCentralityPercentile("CL1"));\r
 696           fHistCentStats->Fill(6.,header->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));\r
 697           fHistCentStats->Fill(7.,header->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));\r
 698           fHistCentStats->Fill(8.,header->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));\r
 699           \r
 700           // centrality QA (V0M)\r
 701           fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());\r
 702           \r
 703           // centrality QA (reference tracks)\r
 704           fHistRefTracks->Fill(0.,header->GetRefMultiplicity());\r
 705           fHistRefTracks->Fill(1.,header->GetRefMultiplicityPos());\r
 706           fHistRefTracks->Fill(2.,header->GetRefMultiplicityNeg());\r
 707           fHistRefTracks->Fill(3.,header->GetTPConlyRefMultiplicity());\r
 708           fHistRefTracks->Fill(4.,header->GetNumberOfITSClusters(0));\r
 709           fHistRefTracks->Fill(5.,header->GetNumberOfITSClusters(1));\r
 710           fHistRefTracks->Fill(6.,header->GetNumberOfITSClusters(2));\r
 711           fHistRefTracks->Fill(7.,header->GetNumberOfITSClusters(3));\r
 712           fHistRefTracks->Fill(8.,header->GetNumberOfITSClusters(4));\r
 713         }//AOD header\r
 714       }//AOD\r
 715 \r
 716       else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){ // centrality class for ESDs or MC-ESDs\r
 717         AliCentrality *centrality = event->GetCentrality();\r
 718         gCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
 719 \r
 720         // QA for centrality estimators\r
 721         fHistCentStats->Fill(0.,centrality->GetCentralityPercentile("V0M"));\r
 722         fHistCentStats->Fill(1.,centrality->GetCentralityPercentile("FMD"));\r
 723         fHistCentStats->Fill(2.,centrality->GetCentralityPercentile("TRK"));\r
 724         fHistCentStats->Fill(3.,centrality->GetCentralityPercentile("TKL"));\r
 725         fHistCentStats->Fill(4.,centrality->GetCentralityPercentile("CL0"));\r
 726         fHistCentStats->Fill(5.,centrality->GetCentralityPercentile("CL1"));\r
 727         fHistCentStats->Fill(6.,centrality->GetCentralityPercentile("V0MvsFMD"));\r
 728         fHistCentStats->Fill(7.,centrality->GetCentralityPercentile("TKLvsV0M"));\r
 729         fHistCentStats->Fill(8.,centrality->GetCentralityPercentile("ZEMvsZDC"));\r
 730 \r
 731         // centrality QA (V0M)\r
 732         fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());\r
 733       }//ESD\r
 734       else if(gAnalysisLevel == "MC"){\r
 735         Double_t gImpactParameter = 0.;\r
 736         if(dynamic_cast<AliMCEvent*>(event)){\r
 737           AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(dynamic_cast<AliMCEvent*>(event)->GenEventHeader());\r
 738           if(headerH){\r
 739             gImpactParameter = headerH->ImpactParameter();\r
 740             gCentrality      = gImpactParameter;\r
 741           }//MC header\r
 742         }//MC event cast\r
 743       }//MC\r
 744       else{\r
 745         gCentrality = -1.;\r
 746       }\r
 747     }\r
 748     \r
 749     // Event Vertex MC\r
 750     if(gAnalysisLevel == "MC"){\r
 751       if(!event) {\r
 752         AliError("mcEvent not available");\r
 753         return 0x0;\r
 754       }\r
 755       \r
 756       AliGenEventHeader *header = dynamic_cast<AliMCEvent*>(event)->GenEventHeader();\r
 757       if(header){  \r
 758         TArrayF gVertexArray;\r
 759         header->PrimaryVertex(gVertexArray);\r
 760         //Printf("Vertex: %lf (x) - %lf (y) - %lf (z)",\r
 761         //gVertexArray.At(0),\r
 762         //gVertexArray.At(1),\r
 763         //gVertexArray.At(2));\r
 764         fHistEventStats->Fill(3,gCentrality); //events with a proper vertex\r
 765         if(TMath::Abs(gVertexArray.At(0)) < fVxMax) {\r
 766           if(TMath::Abs(gVertexArray.At(1)) < fVyMax) {\r
 767             if(TMath::Abs(gVertexArray.At(2)) < fVzMax) {\r
 768               fHistEventStats->Fill(4,gCentrality); //analayzed events\r
 769               fHistVx->Fill(gVertexArray.At(0));\r
 770               fHistVy->Fill(gVertexArray.At(1));\r
 771               fHistVz->Fill(gVertexArray.At(2),gCentrality);\r
 772               \r
 773               // take only events inside centrality class\r
 774               if((fImpactParameterMin < gCentrality) && (fImpactParameterMax > gCentrality)){\r
 775                 fHistEventStats->Fill(5,gCentrality); //events with correct centrality\r
 776                 return gCentrality;         \r
 777               }//centrality class\r
 778             }//Vz cut\r
 779           }//Vy cut\r
 780         }//Vx cut\r
 781       }//header    \r
 782     }//MC\r
 783     \r
 784     // Event Vertex AOD, ESD, ESDMC\r
 785     else{\r
 786       const AliVVertex *vertex = event->GetPrimaryVertex();\r
 787       \r
 788       if(vertex) {\r
 789         Double32_t fCov[6];\r
 790         vertex->GetCovarianceMatrix(fCov);\r
 791         if(vertex->GetNContributors() > 0) {\r
 792           if(fCov[5] != 0) {\r
 793             fHistEventStats->Fill(3,gCentrality); //events with a proper vertex\r
 794             if(TMath::Abs(vertex->GetX()) < fVxMax) {\r
 795               if(TMath::Abs(vertex->GetY()) < fVyMax) {\r
 796                 if(TMath::Abs(vertex->GetZ()) < fVzMax) {\r
 797                   fHistEventStats->Fill(4,gCentrality); //analyzed events\r
 798                   fHistVx->Fill(vertex->GetX());\r
 799                   fHistVy->Fill(vertex->GetY());\r
 800                   fHistVz->Fill(vertex->GetZ(),gCentrality);\r
 801                   \r
 802                   // take only events inside centrality class\r
 803                   if((gCentrality > fCentralityPercentileMin) && (gCentrality < fCentralityPercentileMax)){\r
 804                     fHistEventStats->Fill(5,gCentrality); //events with correct centrality\r
 805                     return gCentrality;         \r
 806                   }//centrality class\r
 807                 }//Vz cut\r
 808               }//Vy cut\r
 809             }//Vx cut\r
 810           }//proper vertex resolution\r
 811         }//proper number of contributors\r
 812       }//vertex object valid\r
 813     }//triggered event \r
 814   }//AOD,ESD,ESDMC\r
 815   \r
 816   // in all other cases return -1 (event not accepted)\r
 817   return -1;\r
 818 }\r
 819 \r
 820 \r
 821 //________________________________________________________________________\r
 822 Double_t AliAnalysisTaskBFPsi::GetRefMultiOrCentrality(AliVEvent *event){\r
 823     // Checks the Event cuts\r
 824     // Fills Event statistics histograms\r
 825   \r
 826   Float_t gCentrality = -1.;\r
 827   Double_t fMultiplicity = -100.;\r
 828   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
 829   if(fEventClass == "Centrality"){\r
 830     \r
 831 \r
 832     if(gAnalysisLevel == "AOD") { //centrality in AOD header\r
 833       AliAODHeader *header = (AliAODHeader*) event->GetHeader();\r
 834       if(header){\r
 835         gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
 836       }//AOD header\r
 837     }//AOD\r
 838     \r
 839     else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){ // centrality class for ESDs or MC-ESDs\r
 840       AliCentrality *centrality = event->GetCentrality();\r
 841       gCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
 842     }//ESD\r
 843     else if(gAnalysisLevel == "MC"){\r
 844       Double_t gImpactParameter = 0.;\r
 845       if(dynamic_cast<AliMCEvent*>(event)){\r
 846         AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(dynamic_cast<AliMCEvent*>(event)->GenEventHeader());      \r
 847         if(headerH){\r
 848           gImpactParameter = headerH->ImpactParameter();\r
 849           gCentrality      = gImpactParameter;\r
 850         }//MC header\r
 851       }//MC event cast\r
 852     }//MC\r
 853     else{\r
 854       gCentrality = -1.;\r
 855     }\r
 856   }//End if "Centrality"\r
 857   if(fEventClass=="Multiplicity"&&gAnalysisLevel == "ESD"){\r
 858     if(dynamic_cast<AliESDEvent*>(event)){\r
 859       fMultiplicity = fESDtrackCuts->GetReferenceMultiplicity(dynamic_cast<AliESDEvent*>(event), AliESDtrackCuts::kTrackletsITSTPC,0.5);\r
 860     }//AliESDevent cast\r
 861   }\r
 862   if(fEventClass=="Multiplicity"&&gAnalysisLevel != "ESD"){\r
 863     AliAODHeader *header = (AliAODHeader*) event->GetHeader();\r
 864     if(header){\r
 865       fMultiplicity = header->GetRefMultiplicity();\r
 866     }//AOD header\r
 867   }\r
 868   Double_t lReturnVal = -100;\r
 869   if(fEventClass=="Multiplicity"){\r
 870     lReturnVal = fMultiplicity;\r
 871   }else if(fEventClass=="Centrality"){\r
 872     lReturnVal = gCentrality;\r
 873   }\r
 874   return lReturnVal;\r
 875 }\r
 876 \r
 877 //________________________________________________________________________\r
 878 Double_t AliAnalysisTaskBFPsi::GetEventPlane(AliVEvent *event){\r
 879   // Get the event plane\r
 880 \r
 881   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
 882 \r
 883   Float_t gVZEROEventPlane    = -10.;\r
 884   Float_t gReactionPlane      = -10.;\r
 885   Double_t qxTot = 0.0, qyTot = 0.0;\r
 886 \r
 887   //MC: from reaction plane\r
 888   if(gAnalysisLevel == "MC"){\r
 889     if(!event) {\r
 890       AliError("mcEvent not available");\r
 891       return 0x0;\r
 892     }\r
 893 \r
 894     if(dynamic_cast<AliMCEvent*>(event)){\r
 895       AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(dynamic_cast<AliMCEvent*>(event)->GenEventHeader());    \r
 896       if (headerH) {\r
 897         gReactionPlane = headerH->ReactionPlaneAngle();\r
 898         //gReactionPlane *= TMath::RadToDeg();\r
 899       }//MC header\r
 900     }//MC event cast\r
 901   }//MC\r
 902   \r
 903   // AOD,ESD,ESDMC: from VZERO Event Plane\r
 904   else{\r
 905    \r
 906     AliEventplane *ep = event->GetEventplane();\r
 907     if(ep){ \r
 908       gVZEROEventPlane = ep->CalculateVZEROEventPlane(event,10,2,qxTot,qyTot);\r
 909       if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();\r
 910       //gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();\r
 911       gReactionPlane = gVZEROEventPlane;\r
 912     }\r
 913   }//AOD,ESD,ESDMC\r
 914 \r
 915   return gReactionPlane;\r
 916 }\r
 917 \r
 918 //________________________________________________________________________\r
 919 Double_t AliAnalysisTaskBFPsi::GetTrackbyTrackCorrectionMatrix( Double_t vEta, \r
 920                                                                 Double_t vPhi, \r
 921                                                                 Double_t vPt, \r
 922                                                                 Short_t vCharge, \r
 923                                                                 Double_t gCentrality) {\r
 924   // -- Get efficiency correction of particle dependent on (eta, phi, pt, charge, centrality) \r
 925 \r
 926   Double_t correction = 1.;\r
 927   Int_t binEta = 0, binPt = 0, binPhi = 0;\r
 928 \r
 929   //Printf("EtaMAx: %lf - EtaMin: %lf - EtaBin: %lf", fEtaMaxForCorrections,fEtaMinForCorrections,fEtaBinForCorrections);\r
 930   if(fEtaBinForCorrections != 0) {\r
 931     Double_t widthEta = (fEtaMaxForCorrections - fEtaMinForCorrections)/fEtaBinForCorrections;\r
 932     if(fEtaMaxForCorrections != fEtaMinForCorrections) \r
 933       binEta = (Int_t)(vEta/widthEta)+1;\r
 934   }\r
 935 \r
 936   if(fPtBinForCorrections != 0) {\r
 937     Double_t widthPt = (fPtMaxForCorrections - fPtMinForCorrections)/fPtBinForCorrections;\r
 938     if(fPtMaxForCorrections != fPtMinForCorrections) \r
 939       binPt = (Int_t)(vPt/widthPt) + 1;\r
 940   }\r
 941  \r
 942   if(fPhiBinForCorrections != 0) {\r
 943     Double_t widthPhi = (fPhiMaxForCorrections - fPhiMinForCorrections)/fPhiBinForCorrections;\r
 944     if(fPhiMaxForCorrections != fPhiMinForCorrections) \r
 945       binPhi = (Int_t)(vPhi/widthPhi)+ 1;\r
 946   }\r
 947 \r
 948   Int_t gCentralityInt = 1;\r
 949   for (Int_t i=0; i<kCENTRALITY; i++){\r
 950     if((centralityArrayForPbPb[i] <= gCentrality)&&(gCentrality <= centralityArrayForPbPb[i+1]))\r
 951       gCentralityInt = i;\r
 952   }\r
 953   \r
 954   if(fHistMatrixCorrectionPlus[gCentralityInt]){\r
 955     if (vCharge > 0) {\r
 956       correction = fHistMatrixCorrectionPlus[gCentralityInt]->GetBinContent(fHistMatrixCorrectionPlus[gCentralityInt]->GetBin(binEta, binPt, binPhi));\r
 957     }\r
 958     if (vCharge < 0) {\r
 959       correction = fHistMatrixCorrectionMinus[gCentralityInt]->GetBinContent(fHistMatrixCorrectionMinus[gCentralityInt]->GetBin(binEta, binPt, binPhi));\r
 960     }\r
 961   }\r
 962   else {\r
 963     correction = 1.;\r
 964   }\r
 965   if (correction == 0.) { \r
 966     AliError(Form("Should not happen : bin content = 0. >> eta: %.2f | phi : %.2f | pt : %.2f | cent %d",vEta, vPhi, vPt, gCentralityInt)); \r
 967     correction = 1.; \r
 968   } \r
 969     \r
 970   return correction;\r
 971 }\r
 972 \r
 973 //________________________________________________________________________\r
 974 TObjArray* AliAnalysisTaskBFPsi::GetAcceptedTracks(AliVEvent *event, Double_t gCentrality, Double_t gReactionPlane){\r
 975   // Returns TObjArray with tracks after all track cuts (only for AOD!)\r
 976   // Fills QA histograms\r
 977 \r
 978   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
 979 \r
 980   //output TObjArray holding all good tracks\r
 981   TObjArray* tracksAccepted = new TObjArray;\r
 982   tracksAccepted->SetOwner(kTRUE);\r
 983 \r
 984   Short_t vCharge;\r
 985   Double_t vEta;\r
 986   Double_t vY;\r
 987   Double_t vPhi;\r
 988   Double_t vPt;\r
 989 \r
 990 \r
 991   if(gAnalysisLevel == "AOD") { // handling of TPC only tracks different in AOD and ESD\r
 992     // Loop over tracks in event\r
 993     for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r
 994       AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));\r
 995       if (!aodTrack) {\r
 996         AliError(Form("Could not receive track %d", iTracks));\r
 997         continue;\r
 998       }\r
 999       \r
1000       // AOD track cuts\r
1001       \r
1002       // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C\r
1003       // take only TPC only tracks \r
1004       //fHistTrackStats->Fill(aodTrack->GetFilterMap());\r
1005       for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){\r
1006         fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));\r
1007       }\r
1008       if(!aodTrack->TestFilterBit(nAODtrackCutBit)) continue;\r
1009       \r
1010       vCharge = aodTrack->Charge();\r
1011       vEta    = aodTrack->Eta();\r
1012       vY      = aodTrack->Y();\r
1013       vPhi    = aodTrack->Phi();// * TMath::RadToDeg();\r
1014       vPt     = aodTrack->Pt();\r
1015       \r
1016       Float_t dcaXY = aodTrack->DCA();      // this is the DCA from global track (not exactly what is cut on)\r
1017       Float_t dcaZ  = aodTrack->ZAtDCA();   // this is the DCA from global track (not exactly what is cut on)\r
1018       \r
1019       \r
1020       // Kinematics cuts from ESD track cuts\r
1021       if( vPt < fPtMin || vPt > fPtMax)      continue;\r
1022       if( vEta < fEtaMin || vEta > fEtaMax)  continue;\r
1023       \r
1024       // Extra DCA cuts (for systematic studies [!= -1])\r
1025       if( fDCAxyCut != -1 && fDCAzCut != -1){\r
1026         if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){\r
1027           continue;  // 2D cut\r
1028         }\r
1029       }\r
1030       \r
1031       // Extra TPC cuts (for systematic studies [!= -1])\r
1032       if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){\r
1033         continue;\r
1034       }\r
1035       if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){\r
1036         continue;\r
1037       }\r
1038       \r
1039       // fill QA histograms\r
1040       fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());\r
1041       fHistDCA->Fill(dcaZ,dcaXY);\r
1042       fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality);\r
1043       fHistPt->Fill(vPt,gCentrality);\r
1044       fHistEta->Fill(vEta,gCentrality);\r
1045       fHistRapidity->Fill(vY,gCentrality);\r
1046       if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
1047       else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
1048       fHistPhi->Fill(vPhi,gCentrality);\r
1049       \r
1050       //=======================================correction\r
1051       Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  \r
1052       \r
1053       // add the track to the TObjArray\r
1054       tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));  \r
1055     }//track loop\r
1056   }// AOD analysis\r
1057 \r
1058 \r
1059   else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD") { // handling of TPC only tracks different in AOD and ESD\r
1060 \r
1061     AliESDtrack *trackTPC   = NULL;\r
1062     AliMCParticle *mcTrack   = NULL;\r
1063 \r
1064     AliMCEvent*  mcEvent     = NULL;\r
1065     \r
1066     // for MC ESDs use also MC information\r
1067     if(gAnalysisLevel == "MCESD"){\r
1068       mcEvent = MCEvent(); \r
1069       if (!mcEvent) {\r
1070         AliError("mcEvent not available");\r
1071         return tracksAccepted;\r
1072       }\r
1073     }\r
1074     \r
1075     // Loop over tracks in event\r
1076     for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r
1077       AliESDtrack* track = dynamic_cast<AliESDtrack *>(event->GetTrack(iTracks));\r
1078       if (!track) {\r
1079         AliError(Form("Could not receive track %d", iTracks));\r
1080         continue;\r
1081       } \r
1082 \r
1083       // for MC ESDs use also MC information --> MC track not used further???\r
1084       if(gAnalysisLevel == "MCESD"){\r
1085         Int_t label = TMath::Abs(track->GetLabel());\r
1086         if(label > mcEvent->GetNumberOfTracks()) continue;\r
1087         if(label > mcEvent->GetNumberOfPrimaries()) continue;\r
1088         \r
1089         mcTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label));\r
1090         if(!mcTrack) continue;\r
1091       }\r
1092 \r
1093       // take only TPC only tracks\r
1094       trackTPC   = new AliESDtrack();\r
1095       if(!track->FillTPCOnlyTrack(*trackTPC)) continue;\r
1096       \r
1097       //ESD track cuts\r
1098       if(fESDtrackCuts) \r
1099         if(!fESDtrackCuts->AcceptTrack(trackTPC)) continue;\r
1100       \r
1101       // fill QA histograms\r
1102       Float_t b[2];\r
1103       Float_t bCov[3];\r
1104       trackTPC->GetImpactParameters(b,bCov);\r
1105       if (bCov[0]<=0 || bCov[2]<=0) {\r
1106         AliDebug(1, "Estimated b resolution lower or equal zero!");\r
1107         bCov[0]=0; bCov[2]=0;\r
1108       }\r
1109       \r
1110       Int_t nClustersTPC = -1;\r
1111       nClustersTPC = trackTPC->GetTPCNclsIter1();   // TPC standalone\r
1112       //nClustersTPC = track->GetTPCclusters(0);   // global track\r
1113       Float_t chi2PerClusterTPC = -1;\r
1114       if (nClustersTPC!=0) {\r
1115         chi2PerClusterTPC = trackTPC->GetTPCchi2Iter1()/Float_t(nClustersTPC);      // TPC standalone\r
1116         //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);     // global track\r
1117       }\r
1118       \r
1119       //===========================PID===============================//             \r
1120       if(fUsePID) {\r
1121         Double_t prob[AliPID::kSPECIES]={0.};\r
1122         Double_t probTPC[AliPID::kSPECIES]={0.};\r
1123         Double_t probTOF[AliPID::kSPECIES]={0.};\r
1124         Double_t probTPCTOF[AliPID::kSPECIES]={0.};\r
1125         \r
1126         Double_t nSigma = 0.;\r
1127         UInt_t detUsedTPC = 0;\r
1128         UInt_t detUsedTOF = 0;\r
1129         UInt_t detUsedTPCTOF = 0;\r
1130         \r
1131         //Decide what detector configuration we want to use\r
1132         switch(fPidDetectorConfig) {\r
1133         case kTPCpid:\r
1134           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);\r
1135           nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest));\r
1136           detUsedTPC = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPC);\r
1137           for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
1138             prob[iSpecies] = probTPC[iSpecies];\r
1139           break;\r
1140         case kTOFpid:\r
1141           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);\r
1142           nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest));\r
1143           detUsedTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTOF);\r
1144           for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
1145             prob[iSpecies] = probTOF[iSpecies];\r
1146           break;\r
1147         case kTPCTOF:\r
1148           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);\r
1149           detUsedTPCTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPCTOF);\r
1150           for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
1151             prob[iSpecies] = probTPCTOF[iSpecies];\r
1152           break;\r
1153         default:\r
1154           break;\r
1155         }//end switch: define detector mask\r
1156         \r
1157         //Filling the PID QA\r
1158         Double_t tofTime = -999., length = 999., tof = -999.;\r
1159         Double_t c = TMath::C()*1.E-9;// m/ns\r
1160         Double_t beta = -999.;\r
1161         Double_t  nSigmaTOFForParticleOfInterest = -999.;\r
1162         if ( (track->IsOn(AliESDtrack::kTOFin)) &&\r
1163              (track->IsOn(AliESDtrack::kTIME))  ) { \r
1164           tofTime = track->GetTOFsignal();//in ps\r
1165           length = track->GetIntegratedLength();\r
1166           tof = tofTime*1E-3; // ns     \r
1167           \r
1168           if (tof <= 0) {\r
1169             //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");\r
1170             continue;\r
1171           }\r
1172           if (length <= 0){\r
1173             //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");\r
1174             continue;\r
1175           }\r
1176           \r
1177           length = length*0.01; // in meters\r
1178           tof = tof*c;\r
1179           beta = length/tof;\r
1180           \r
1181           nSigmaTOFForParticleOfInterest = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest);\r
1182           fHistBetavsPTOFbeforePID ->Fill(track->P()*track->Charge(),beta);\r
1183           fHistProbTOFvsPtbeforePID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);\r
1184           fHistNSigmaTOFvsPtbeforePID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);\r
1185         }//TOF signal \r
1186         \r
1187         \r
1188         Double_t  nSigmaTPCForParticleOfInterest = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest);\r
1189         fHistdEdxVsPTPCbeforePID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());\r
1190         fHistProbTPCvsPtbeforePID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); \r
1191         fHistNSigmaTPCvsPtbeforePID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); \r
1192         fHistProbTPCTOFvsPtbeforePID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);\r
1193         //end of QA-before pid\r
1194         \r
1195         if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {\r
1196           //Make the decision based on the n-sigma\r
1197           if(fUsePIDnSigma) {\r
1198             if(nSigma > fPIDNSigma) continue;}\r
1199           \r
1200           //Make the decision based on the bayesian\r
1201           else if(fUsePIDPropabilities) {\r
1202             if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;\r
1203             if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue;      \r
1204           }\r
1205           \r
1206           //Fill QA after the PID\r
1207           fHistBetavsPTOFafterPID ->Fill(track->P()*track->Charge(),beta);\r
1208           fHistProbTOFvsPtafterPID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);\r
1209           fHistNSigmaTOFvsPtafterPID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);\r
1210           \r
1211           fHistdEdxVsPTPCafterPID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());\r
1212           fHistProbTPCvsPtafterPID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); \r
1213           fHistProbTPCTOFvsPtafterPID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);\r
1214           fHistNSigmaTPCvsPtafterPID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); \r
1215         }\r
1216       }\r
1217       //===========================PID===============================//\r
1218       vCharge = trackTPC->Charge();\r
1219       vY      = trackTPC->Y();\r
1220       vEta    = trackTPC->Eta();\r
1221       vPhi    = trackTPC->Phi();// * TMath::RadToDeg();\r
1222       vPt     = trackTPC->Pt();\r
1223       fHistClus->Fill(trackTPC->GetITSclusters(0),nClustersTPC);\r
1224       fHistDCA->Fill(b[1],b[0]);\r
1225       fHistChi2->Fill(chi2PerClusterTPC,gCentrality);\r
1226       fHistPt->Fill(vPt,gCentrality);\r
1227       fHistEta->Fill(vEta,gCentrality);\r
1228       fHistPhi->Fill(vPhi,gCentrality);\r
1229       fHistRapidity->Fill(vY,gCentrality);\r
1230       if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
1231       else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
1232       \r
1233       //=======================================correction\r
1234       Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  \r
1235 \r
1236       // add the track to the TObjArray\r
1237       tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));   \r
1238 \r
1239       delete trackTPC;\r
1240     }//track loop\r
1241   }// ESD analysis\r
1242 \r
1243   else if(gAnalysisLevel == "MC"){\r
1244     if(!event) {\r
1245       AliError("mcEvent not available");\r
1246       return 0x0;\r
1247     }\r
1248 \r
1249     // Loop over tracks in event\r
1250     for (Int_t iTracks = 0; iTracks < dynamic_cast<AliMCEvent*>(event)->GetNumberOfPrimaries(); iTracks++) {\r
1251       AliMCParticle* track = dynamic_cast<AliMCParticle *>(event->GetTrack(iTracks));\r
1252       if (!track) {\r
1253         AliError(Form("Could not receive particle %d", iTracks));\r
1254         continue;\r
1255       }\r
1256             \r
1257       //exclude non stable particles\r
1258       if(!(dynamic_cast<AliMCEvent*>(event)->IsPhysicalPrimary(iTracks))) continue;\r
1259 \r
1260       vEta    = track->Eta();\r
1261       vPt     = track->Pt();\r
1262       vY      = track->Y();\r
1263       \r
1264       if( vPt < fPtMin || vPt > fPtMax)      \r
1265         continue;\r
1266       if (!fUsePID) {\r
1267         if( vEta < fEtaMin || vEta > fEtaMax)  continue;\r
1268       }\r
1269       else if (fUsePID){\r
1270         if( vY < fEtaMin || vY > fEtaMax)  continue;\r
1271       }\r
1272       \r
1273       //analyze one set of particles\r
1274       if(fUseMCPdgCode) {\r
1275         TParticle *particle = track->Particle();\r
1276         if(!particle) continue;\r
1277         \r
1278         Int_t gPdgCode = particle->GetPdgCode();\r
1279         if(TMath::Abs(fPDGCodeToBeAnalyzed) != TMath::Abs(gPdgCode)) \r
1280           continue;\r
1281       }\r
1282       \r
1283       //Use the acceptance parameterization\r
1284       if(fAcceptanceParameterization) {\r
1285         Double_t gRandomNumber = gRandom->Rndm();\r
1286         if(gRandomNumber > fAcceptanceParameterization->Eval(track->Pt())) \r
1287           continue;\r
1288       }\r
1289       \r
1290       //Exclude resonances\r
1291       if(fExcludeResonancesInMC) {\r
1292         TParticle *particle = track->Particle();\r
1293         if(!particle) continue;\r
1294         \r
1295         Bool_t kExcludeParticle = kFALSE;\r
1296         Int_t gMotherIndex = particle->GetFirstMother();\r
1297         if(gMotherIndex != -1) {\r
1298           AliMCParticle* motherTrack = dynamic_cast<AliMCParticle *>(event->GetTrack(gMotherIndex));\r
1299           if(motherTrack) {\r
1300             TParticle *motherParticle = motherTrack->Particle();\r
1301             if(motherParticle) {\r
1302               Int_t pdgCodeOfMother = motherParticle->GetPdgCode();\r
1303               //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {\r
1304               if(pdgCodeOfMother == 113) {\r
1305                 kExcludeParticle = kTRUE;\r
1306               }\r
1307             }\r
1308           }\r
1309         }\r
1310         \r
1311         //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi\r
1312         if(kExcludeParticle) continue;\r
1313       }\r
1314       \r
1315       vCharge = track->Charge();\r
1316       vPhi    = track->Phi();\r
1317       //Printf("phi (before): %lf",vPhi);\r
1318       \r
1319       fHistPt->Fill(vPt,gCentrality);\r
1320       fHistEta->Fill(vEta,gCentrality);\r
1321       fHistPhi->Fill(vPhi,gCentrality);\r
1322       //fHistPhi->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r
1323       fHistRapidity->Fill(vY,gCentrality);\r
1324       //if(vCharge > 0) fHistPhiPos->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r
1325       //else if(vCharge < 0) fHistPhiNeg->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r
1326       if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
1327       else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
1328       \r
1329       //Flow after burner\r
1330       if(fUseFlowAfterBurner) {\r
1331         Double_t precisionPhi = 0.001;\r
1332         Int_t maxNumberOfIterations = 100;\r
1333         \r
1334         Double_t phi0 = vPhi;\r
1335         Double_t gV2 = fDifferentialV2->Eval(vPt);\r
1336         \r
1337         for (Int_t j = 0; j < maxNumberOfIterations; j++) {\r
1338           Double_t phiprev = vPhi;\r
1339           Double_t fl = vPhi - phi0 + gV2*TMath::Sin(2.*(vPhi - gReactionPlane*TMath::DegToRad()));\r
1340           Double_t fp = 1.0 + 2.0*gV2*TMath::Cos(2.*(vPhi - gReactionPlane*TMath::DegToRad())); \r
1341           vPhi -= fl/fp;\r
1342           if (TMath::AreEqualAbs(phiprev,vPhi,precisionPhi)) break;\r
1343         }\r
1344         //Printf("phi (after): %lf\n",vPhi);\r
1345         Double_t vDeltaphiBefore = phi0 - gReactionPlane*TMath::DegToRad();\r
1346         if(vDeltaphiBefore < 0) vDeltaphiBefore += 2*TMath::Pi();\r
1347         fHistPhiBefore->Fill(vDeltaphiBefore,gCentrality);\r
1348         \r
1349         Double_t vDeltaphiAfter = vPhi - gReactionPlane*TMath::DegToRad();\r
1350         if(vDeltaphiAfter < 0) vDeltaphiAfter += 2*TMath::Pi();\r
1351         fHistPhiAfter->Fill(vDeltaphiAfter,gCentrality);\r
1352       \r
1353       }\r
1354       \r
1355       //vPhi *= TMath::RadToDeg();\r
1356 \r
1357       //=======================================correction\r
1358       Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  \r
1359       \r
1360       tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction)); \r
1361       \r
1362     } //track loop\r
1363   }//MC\r
1364   \r
1365   return tracksAccepted;  \r
1366 }\r
1367 //________________________________________________________________________\r
1368 TObjArray* AliAnalysisTaskBFPsi::GetShuffledTracks(TObjArray *tracks, Double_t gCentrality){\r
1369   // Clones TObjArray and returns it with tracks after shuffling the charges\r
1370 \r
1371   TObjArray* tracksShuffled = new TObjArray;\r
1372   tracksShuffled->SetOwner(kTRUE);\r
1373 \r
1374   vector<Short_t> *chargeVector = new vector<Short_t>;   //original charge of accepted tracks \r
1375 \r
1376   for (Int_t i=0; i<tracks->GetEntriesFast(); i++)\r
1377   {\r
1378     AliVParticle* track = (AliVParticle*) tracks->At(i);\r
1379     chargeVector->push_back(track->Charge());\r
1380   }  \r
1381  \r
1382   random_shuffle(chargeVector->begin(), chargeVector->end());\r
1383   \r
1384   for(Int_t i = 0; i < tracks->GetEntriesFast(); i++){\r
1385     AliVParticle* track = (AliVParticle*) tracks->At(i);\r
1386     //==============================correction\r
1387     Double_t correction = GetTrackbyTrackCorrectionMatrix(track->Eta(), track->Phi(),track->Pt(), chargeVector->at(i), gCentrality);\r
1388     tracksShuffled->Add(new AliBFBasicParticle(track->Eta(), track->Phi(), track->Pt(),chargeVector->at(i), correction));\r
1389   }\r
1390 \r
1391   delete chargeVector;\r
1392    \r
1393   return tracksShuffled;\r
1394 }\r
1395 \r
1396 \r
1397 //________________________________________________________________________\r
1398 void  AliAnalysisTaskBFPsi::FinishTaskOutput(){\r
1399   //Printf("END BF");\r
1400 \r
1401   if (!fBalance) {\r
1402     AliError("fBalance not available");\r
1403     return;\r
1404   }  \r
1405   if(fRunShuffling) {\r
1406     if (!fShuffledBalance) {\r
1407       AliError("fShuffledBalance not available");\r
1408       return;\r
1409     }\r
1410   }\r
1411 \r
1412 }\r
1413 \r
1414 //________________________________________________________________________\r
1415 void AliAnalysisTaskBFPsi::Terminate(Option_t *) {\r
1416   // Draw result to the screen\r
1417   // Called once at the end of the query\r
1418 \r
1419   // not implemented ...\r
1420 \r
1421 }\r