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       AliGenEventHeader *header = dynamic_cast<AliMCEvent*>(event)->GenEventHeader();\r
 752       if(header){  \r
 753         TArrayF gVertexArray;\r
 754         header->PrimaryVertex(gVertexArray);\r
 755         //Printf("Vertex: %lf (x) - %lf (y) - %lf (z)",\r
 756         //gVertexArray.At(0),\r
 757         //gVertexArray.At(1),\r
 758         //gVertexArray.At(2));\r
 759         fHistEventStats->Fill(3,gCentrality); //events with a proper vertex\r
 760         if(TMath::Abs(gVertexArray.At(0)) < fVxMax) {\r
 761           if(TMath::Abs(gVertexArray.At(1)) < fVyMax) {\r
 762             if(TMath::Abs(gVertexArray.At(2)) < fVzMax) {\r
 763               fHistEventStats->Fill(4,gCentrality); //analayzed events\r
 764               fHistVx->Fill(gVertexArray.At(0));\r
 765               fHistVy->Fill(gVertexArray.At(1));\r
 766               fHistVz->Fill(gVertexArray.At(2),gCentrality);\r
 767               \r
 768               // take only events inside centrality class\r
 769               if((fImpactParameterMin < gCentrality) && (fImpactParameterMax > gCentrality)){\r
 770                 fHistEventStats->Fill(5,gCentrality); //events with correct centrality\r
 771                 return gCentrality;         \r
 772               }//centrality class\r
 773             }//Vz cut\r
 774           }//Vy cut\r
 775         }//Vx cut\r
 776       }//header    \r
 777     }//MC\r
 778     \r
 779     // Event Vertex AOD, ESD, ESDMC\r
 780     else{\r
 781       const AliVVertex *vertex = event->GetPrimaryVertex();\r
 782       \r
 783       if(vertex) {\r
 784         Double32_t fCov[6];\r
 785         vertex->GetCovarianceMatrix(fCov);\r
 786         if(vertex->GetNContributors() > 0) {\r
 787           if(fCov[5] != 0) {\r
 788             fHistEventStats->Fill(3,gCentrality); //events with a proper vertex\r
 789             if(TMath::Abs(vertex->GetX()) < fVxMax) {\r
 790               if(TMath::Abs(vertex->GetY()) < fVyMax) {\r
 791                 if(TMath::Abs(vertex->GetZ()) < fVzMax) {\r
 792                   fHistEventStats->Fill(4,gCentrality); //analyzed events\r
 793                   fHistVx->Fill(vertex->GetX());\r
 794                   fHistVy->Fill(vertex->GetY());\r
 795                   fHistVz->Fill(vertex->GetZ(),gCentrality);\r
 796                   \r
 797                   // take only events inside centrality class\r
 798                   if((gCentrality > fCentralityPercentileMin) && (gCentrality < fCentralityPercentileMax)){\r
 799                     fHistEventStats->Fill(5,gCentrality); //events with correct centrality\r
 800                     return gCentrality;         \r
 801                   }//centrality class\r
 802                 }//Vz cut\r
 803               }//Vy cut\r
 804             }//Vx cut\r
 805           }//proper vertex resolution\r
 806         }//proper number of contributors\r
 807       }//vertex object valid\r
 808     }//triggered event \r
 809   }//AOD,ESD,ESDMC\r
 810   \r
 811   // in all other cases return -1 (event not accepted)\r
 812   return -1;\r
 813 }\r
 814 \r
 815 \r
 816 //________________________________________________________________________\r
 817 Double_t AliAnalysisTaskBFPsi::GetRefMultiOrCentrality(AliVEvent *event){\r
 818     // Checks the Event cuts\r
 819     // Fills Event statistics histograms\r
 820   \r
 821   Float_t gCentrality = -1.;\r
 822   Double_t fMultiplicity = -100.;\r
 823   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
 824   if(fEventClass == "Centrality"){\r
 825     \r
 826 \r
 827     if(gAnalysisLevel == "AOD") { //centrality in AOD header\r
 828       AliAODHeader *header = (AliAODHeader*) event->GetHeader();\r
 829       if(header){\r
 830         gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
 831       }//AOD header\r
 832     }//AOD\r
 833     \r
 834     else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){ // centrality class for ESDs or MC-ESDs\r
 835       AliCentrality *centrality = event->GetCentrality();\r
 836       gCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
 837     }//ESD\r
 838     else if(gAnalysisLevel == "MC"){\r
 839       Double_t gImpactParameter = 0.;\r
 840       if(dynamic_cast<AliMCEvent*>(event)){\r
 841         AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(dynamic_cast<AliMCEvent*>(event)->GenEventHeader());      \r
 842         if(headerH){\r
 843           gImpactParameter = headerH->ImpactParameter();\r
 844           gCentrality      = gImpactParameter;\r
 845         }//MC header\r
 846       }//MC event cast\r
 847     }//MC\r
 848     else{\r
 849       gCentrality = -1.;\r
 850     }\r
 851   }//End if "Centrality"\r
 852   if(fEventClass=="Multiplicity"&&gAnalysisLevel == "ESD"){\r
 853     if(dynamic_cast<AliESDEvent*>(event)){\r
 854       fMultiplicity = fESDtrackCuts->GetReferenceMultiplicity(dynamic_cast<AliESDEvent*>(event), AliESDtrackCuts::kTrackletsITSTPC,0.5);\r
 855     }//AliESDevent cast\r
 856   }\r
 857   if(fEventClass=="Multiplicity"&&gAnalysisLevel != "ESD"){\r
 858     AliAODHeader *header = (AliAODHeader*) event->GetHeader();\r
 859     if(header){\r
 860       fMultiplicity = header->GetRefMultiplicity();\r
 861     }//AOD header\r
 862   }\r
 863   Double_t lReturnVal = -100;\r
 864   if(fEventClass=="Multiplicity"){\r
 865     lReturnVal = fMultiplicity;\r
 866   }else if(fEventClass=="Centrality"){\r
 867     lReturnVal = gCentrality;\r
 868   }\r
 869   return lReturnVal;\r
 870 }\r
 871 \r
 872 //________________________________________________________________________\r
 873 Double_t AliAnalysisTaskBFPsi::GetEventPlane(AliVEvent *event){\r
 874   // Get the event plane\r
 875 \r
 876   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
 877 \r
 878   Float_t gVZEROEventPlane    = -10.;\r
 879   Float_t gReactionPlane      = -10.;\r
 880   Double_t qxTot = 0.0, qyTot = 0.0;\r
 881 \r
 882   //MC: from reaction plane\r
 883   if(gAnalysisLevel == "MC"){\r
 884     if(dynamic_cast<AliMCEvent*>(event)){\r
 885       AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(dynamic_cast<AliMCEvent*>(event)->GenEventHeader());    \r
 886       if (headerH) {\r
 887         gReactionPlane = headerH->ReactionPlaneAngle();\r
 888         //gReactionPlane *= TMath::RadToDeg();\r
 889       }//MC header\r
 890     }//MC event cast\r
 891   }//MC\r
 892   \r
 893   // AOD,ESD,ESDMC: from VZERO Event Plane\r
 894   else{\r
 895    \r
 896     AliEventplane *ep = event->GetEventplane();\r
 897     if(ep){ \r
 898       gVZEROEventPlane = ep->CalculateVZEROEventPlane(event,10,2,qxTot,qyTot);\r
 899       if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();\r
 900       //gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();\r
 901       gReactionPlane = gVZEROEventPlane;\r
 902     }\r
 903   }//AOD,ESD,ESDMC\r
 904 \r
 905   return gReactionPlane;\r
 906 }\r
 907 \r
 908 //________________________________________________________________________\r
 909 Double_t AliAnalysisTaskBFPsi::GetTrackbyTrackCorrectionMatrix( Double_t vEta, \r
 910                                                                 Double_t vPhi, \r
 911                                                                 Double_t vPt, \r
 912                                                                 Short_t vCharge, \r
 913                                                                 Double_t gCentrality) {\r
 914   // -- Get efficiency correction of particle dependent on (eta, phi, pt, charge, centrality) \r
 915 \r
 916   Double_t correction = 1.;\r
 917   Int_t binEta = 0, binPt = 0, binPhi = 0;\r
 918 \r
 919   //Printf("EtaMAx: %lf - EtaMin: %lf - EtaBin: %lf", fEtaMaxForCorrections,fEtaMinForCorrections,fEtaBinForCorrections);\r
 920   if(fEtaBinForCorrections != 0) {\r
 921     Double_t widthEta = (fEtaMaxForCorrections - fEtaMinForCorrections)/fEtaBinForCorrections;\r
 922     if(fEtaMaxForCorrections != fEtaMinForCorrections) \r
 923       binEta = (Int_t)(vEta/widthEta)+1;\r
 924   }\r
 925 \r
 926   if(fPtBinForCorrections != 0) {\r
 927     Double_t widthPt = (fPtMaxForCorrections - fPtMinForCorrections)/fPtBinForCorrections;\r
 928     if(fPtMaxForCorrections != fPtMinForCorrections) \r
 929       binPt = (Int_t)(vPt/widthPt) + 1;\r
 930   }\r
 931  \r
 932   if(fPhiBinForCorrections != 0) {\r
 933     Double_t widthPhi = (fPhiMaxForCorrections - fPhiMinForCorrections)/fPhiBinForCorrections;\r
 934     if(fPhiMaxForCorrections != fPhiMinForCorrections) \r
 935       binPhi = (Int_t)(vPhi/widthPhi)+ 1;\r
 936   }\r
 937 \r
 938   Int_t gCentralityInt = 1;\r
 939   for (Int_t i=0; i<=kCENTRALITY; i++){\r
 940     if((centralityArrayForPbPb[i] <= gCentrality)&&(gCentrality <= centralityArrayForPbPb[i+1]))\r
 941       gCentralityInt = i;\r
 942   }\r
 943   \r
 944   if(fHistMatrixCorrectionPlus[gCentralityInt]){\r
 945     if (vCharge > 0) {\r
 946       correction = fHistMatrixCorrectionPlus[gCentralityInt]->GetBinContent(fHistMatrixCorrectionPlus[gCentralityInt]->GetBin(binEta, binPt, binPhi));\r
 947     }\r
 948     if (vCharge < 0) {\r
 949       correction = fHistMatrixCorrectionMinus[gCentralityInt]->GetBinContent(fHistMatrixCorrectionMinus[gCentralityInt]->GetBin(binEta, binPt, binPhi));\r
 950     }\r
 951   }\r
 952   else {\r
 953     correction = 1.;\r
 954   }\r
 955   if (correction == 0.) { \r
 956     AliError(Form("Should not happen : bin content = 0. >> eta: %.2f | phi : %.2f | pt : %.2f | cent %d",vEta, vPhi, vPt, gCentralityInt)); \r
 957     correction = 1.; \r
 958   } \r
 959     \r
 960   return correction;\r
 961 }\r
 962 \r
 963 //________________________________________________________________________\r
 964 TObjArray* AliAnalysisTaskBFPsi::GetAcceptedTracks(AliVEvent *event, Double_t gCentrality, Double_t gReactionPlane){\r
 965   // Returns TObjArray with tracks after all track cuts (only for AOD!)\r
 966   // Fills QA histograms\r
 967 \r
 968   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
 969 \r
 970   //output TObjArray holding all good tracks\r
 971   TObjArray* tracksAccepted = new TObjArray;\r
 972   tracksAccepted->SetOwner(kTRUE);\r
 973 \r
 974   Short_t vCharge;\r
 975   Double_t vEta;\r
 976   Double_t vY;\r
 977   Double_t vPhi;\r
 978   Double_t vPt;\r
 979 \r
 980 \r
 981   if(gAnalysisLevel == "AOD") { // handling of TPC only tracks different in AOD and ESD\r
 982     // Loop over tracks in event\r
 983     for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r
 984       AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));\r
 985       if (!aodTrack) {\r
 986         AliError(Form("Could not receive track %d", iTracks));\r
 987         continue;\r
 988       }\r
 989       \r
 990       // AOD track cuts\r
 991       \r
 992       // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C\r
 993       // take only TPC only tracks \r
 994       //fHistTrackStats->Fill(aodTrack->GetFilterMap());\r
 995       for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){\r
 996         fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));\r
 997       }\r
 998       if(!aodTrack->TestFilterBit(nAODtrackCutBit)) continue;\r
 999       \r
1000       vCharge = aodTrack->Charge();\r
1001       vEta    = aodTrack->Eta();\r
1002       vY      = aodTrack->Y();\r
1003       vPhi    = aodTrack->Phi();// * TMath::RadToDeg();\r
1004       vPt     = aodTrack->Pt();\r
1005       \r
1006       Float_t dcaXY = aodTrack->DCA();      // this is the DCA from global track (not exactly what is cut on)\r
1007       Float_t dcaZ  = aodTrack->ZAtDCA();   // this is the DCA from global track (not exactly what is cut on)\r
1008       \r
1009       \r
1010       // Kinematics cuts from ESD track cuts\r
1011       if( vPt < fPtMin || vPt > fPtMax)      continue;\r
1012       if( vEta < fEtaMin || vEta > fEtaMax)  continue;\r
1013       \r
1014       // Extra DCA cuts (for systematic studies [!= -1])\r
1015       if( fDCAxyCut != -1 && fDCAzCut != -1){\r
1016         if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){\r
1017           continue;  // 2D cut\r
1018         }\r
1019       }\r
1020       \r
1021       // Extra TPC cuts (for systematic studies [!= -1])\r
1022       if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){\r
1023         continue;\r
1024       }\r
1025       if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){\r
1026         continue;\r
1027       }\r
1028       \r
1029       // fill QA histograms\r
1030       fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());\r
1031       fHistDCA->Fill(dcaZ,dcaXY);\r
1032       fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality);\r
1033       fHistPt->Fill(vPt,gCentrality);\r
1034       fHistEta->Fill(vEta,gCentrality);\r
1035       fHistRapidity->Fill(vY,gCentrality);\r
1036       if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
1037       else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
1038       fHistPhi->Fill(vPhi,gCentrality);\r
1039       \r
1040       //=======================================correction\r
1041       Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  \r
1042       \r
1043       // add the track to the TObjArray\r
1044       tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));  \r
1045     }//track loop\r
1046   }// AOD analysis\r
1047 \r
1048 \r
1049   else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD") { // handling of TPC only tracks different in AOD and ESD\r
1050 \r
1051     AliESDtrack *trackTPC   = NULL;\r
1052     AliMCParticle *mcTrack   = NULL;\r
1053 \r
1054     AliMCEvent*  mcEvent     = NULL;\r
1055     \r
1056     // for MC ESDs use also MC information\r
1057     if(gAnalysisLevel == "MCESD"){\r
1058       mcEvent = MCEvent(); \r
1059       if (!mcEvent) {\r
1060         AliError("mcEvent not available");\r
1061         return tracksAccepted;\r
1062       }\r
1063     }\r
1064     \r
1065     // Loop over tracks in event\r
1066     for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r
1067       AliESDtrack* track = dynamic_cast<AliESDtrack *>(event->GetTrack(iTracks));\r
1068       if (!track) {\r
1069         AliError(Form("Could not receive track %d", iTracks));\r
1070         continue;\r
1071       } \r
1072 \r
1073       // for MC ESDs use also MC information --> MC track not used further???\r
1074       if(gAnalysisLevel == "MCESD"){\r
1075         Int_t label = TMath::Abs(track->GetLabel());\r
1076         if(label > mcEvent->GetNumberOfTracks()) continue;\r
1077         if(label > mcEvent->GetNumberOfPrimaries()) continue;\r
1078         \r
1079         mcTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label));\r
1080         if(!mcTrack) continue;\r
1081       }\r
1082 \r
1083       // take only TPC only tracks\r
1084       trackTPC   = new AliESDtrack();\r
1085       if(!track->FillTPCOnlyTrack(*trackTPC)) continue;\r
1086       \r
1087       //ESD track cuts\r
1088       if(fESDtrackCuts) \r
1089         if(!fESDtrackCuts->AcceptTrack(trackTPC)) continue;\r
1090       \r
1091       // fill QA histograms\r
1092       Float_t b[2];\r
1093       Float_t bCov[3];\r
1094       trackTPC->GetImpactParameters(b,bCov);\r
1095       if (bCov[0]<=0 || bCov[2]<=0) {\r
1096         AliDebug(1, "Estimated b resolution lower or equal zero!");\r
1097         bCov[0]=0; bCov[2]=0;\r
1098       }\r
1099       \r
1100       Int_t nClustersTPC = -1;\r
1101       nClustersTPC = trackTPC->GetTPCNclsIter1();   // TPC standalone\r
1102       //nClustersTPC = track->GetTPCclusters(0);   // global track\r
1103       Float_t chi2PerClusterTPC = -1;\r
1104       if (nClustersTPC!=0) {\r
1105         chi2PerClusterTPC = trackTPC->GetTPCchi2Iter1()/Float_t(nClustersTPC);      // TPC standalone\r
1106         //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);     // global track\r
1107       }\r
1108       \r
1109       //===========================PID===============================//             \r
1110       if(fUsePID) {\r
1111         Double_t prob[AliPID::kSPECIES]={0.};\r
1112         Double_t probTPC[AliPID::kSPECIES]={0.};\r
1113         Double_t probTOF[AliPID::kSPECIES]={0.};\r
1114         Double_t probTPCTOF[AliPID::kSPECIES]={0.};\r
1115         \r
1116         Double_t nSigma = 0.;\r
1117         UInt_t detUsedTPC = 0;\r
1118         UInt_t detUsedTOF = 0;\r
1119         UInt_t detUsedTPCTOF = 0;\r
1120         \r
1121         //Decide what detector configuration we want to use\r
1122         switch(fPidDetectorConfig) {\r
1123         case kTPCpid:\r
1124           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);\r
1125           nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest));\r
1126           detUsedTPC = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPC);\r
1127           for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
1128             prob[iSpecies] = probTPC[iSpecies];\r
1129           break;\r
1130         case kTOFpid:\r
1131           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);\r
1132           nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest));\r
1133           detUsedTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTOF);\r
1134           for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
1135             prob[iSpecies] = probTOF[iSpecies];\r
1136           break;\r
1137         case kTPCTOF:\r
1138           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);\r
1139           detUsedTPCTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPCTOF);\r
1140           for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
1141             prob[iSpecies] = probTPCTOF[iSpecies];\r
1142           break;\r
1143         default:\r
1144           break;\r
1145         }//end switch: define detector mask\r
1146         \r
1147         //Filling the PID QA\r
1148         Double_t tofTime = -999., length = 999., tof = -999.;\r
1149         Double_t c = TMath::C()*1.E-9;// m/ns\r
1150         Double_t beta = -999.;\r
1151         Double_t  nSigmaTOFForParticleOfInterest = -999.;\r
1152         if ( (track->IsOn(AliESDtrack::kTOFin)) &&\r
1153              (track->IsOn(AliESDtrack::kTIME))  ) { \r
1154           tofTime = track->GetTOFsignal();//in ps\r
1155           length = track->GetIntegratedLength();\r
1156           tof = tofTime*1E-3; // ns     \r
1157           \r
1158           if (tof <= 0) {\r
1159             //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");\r
1160             continue;\r
1161           }\r
1162           if (length <= 0){\r
1163             //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");\r
1164             continue;\r
1165           }\r
1166           \r
1167           length = length*0.01; // in meters\r
1168           tof = tof*c;\r
1169           beta = length/tof;\r
1170           \r
1171           nSigmaTOFForParticleOfInterest = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest);\r
1172           fHistBetavsPTOFbeforePID ->Fill(track->P()*track->Charge(),beta);\r
1173           fHistProbTOFvsPtbeforePID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);\r
1174           fHistNSigmaTOFvsPtbeforePID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);\r
1175         }//TOF signal \r
1176         \r
1177         \r
1178         Double_t  nSigmaTPCForParticleOfInterest = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest);\r
1179         fHistdEdxVsPTPCbeforePID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());\r
1180         fHistProbTPCvsPtbeforePID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); \r
1181         fHistNSigmaTPCvsPtbeforePID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); \r
1182         fHistProbTPCTOFvsPtbeforePID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);\r
1183         //end of QA-before pid\r
1184         \r
1185         if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {\r
1186           //Make the decision based on the n-sigma\r
1187           if(fUsePIDnSigma) {\r
1188             if(nSigma > fPIDNSigma) continue;}\r
1189           \r
1190           //Make the decision based on the bayesian\r
1191           else if(fUsePIDPropabilities) {\r
1192             if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;\r
1193             if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue;      \r
1194           }\r
1195           \r
1196           //Fill QA after the PID\r
1197           fHistBetavsPTOFafterPID ->Fill(track->P()*track->Charge(),beta);\r
1198           fHistProbTOFvsPtafterPID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);\r
1199           fHistNSigmaTOFvsPtafterPID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);\r
1200           \r
1201           fHistdEdxVsPTPCafterPID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());\r
1202           fHistProbTPCvsPtafterPID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); \r
1203           fHistProbTPCTOFvsPtafterPID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);\r
1204           fHistNSigmaTPCvsPtafterPID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); \r
1205         }\r
1206       }\r
1207       //===========================PID===============================//\r
1208       vCharge = trackTPC->Charge();\r
1209       vY      = trackTPC->Y();\r
1210       vEta    = trackTPC->Eta();\r
1211       vPhi    = trackTPC->Phi();// * TMath::RadToDeg();\r
1212       vPt     = trackTPC->Pt();\r
1213       fHistClus->Fill(trackTPC->GetITSclusters(0),nClustersTPC);\r
1214       fHistDCA->Fill(b[1],b[0]);\r
1215       fHistChi2->Fill(chi2PerClusterTPC,gCentrality);\r
1216       fHistPt->Fill(vPt,gCentrality);\r
1217       fHistEta->Fill(vEta,gCentrality);\r
1218       fHistPhi->Fill(vPhi,gCentrality);\r
1219       fHistRapidity->Fill(vY,gCentrality);\r
1220       if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
1221       else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
1222       \r
1223       //=======================================correction\r
1224       Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  \r
1225 \r
1226       // add the track to the TObjArray\r
1227       tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));   \r
1228 \r
1229       delete trackTPC;\r
1230     }//track loop\r
1231   }// ESD analysis\r
1232 \r
1233   else if(gAnalysisLevel == "MC"){\r
1234     \r
1235     // Loop over tracks in event\r
1236     for (Int_t iTracks = 0; iTracks < dynamic_cast<AliMCEvent*>(event)->GetNumberOfPrimaries(); iTracks++) {\r
1237       AliMCParticle* track = dynamic_cast<AliMCParticle *>(event->GetTrack(iTracks));\r
1238       if (!track) {\r
1239         AliError(Form("Could not receive particle %d", iTracks));\r
1240         continue;\r
1241       }\r
1242             \r
1243       //exclude non stable particles\r
1244       if(!(dynamic_cast<AliMCEvent*>(event)->IsPhysicalPrimary(iTracks))) continue;\r
1245 \r
1246       vEta    = track->Eta();\r
1247       vPt     = track->Pt();\r
1248       vY      = track->Y();\r
1249       \r
1250       if( vPt < fPtMin || vPt > fPtMax)      \r
1251         continue;\r
1252       if (!fUsePID) {\r
1253         if( vEta < fEtaMin || vEta > fEtaMax)  continue;\r
1254       }\r
1255       else if (fUsePID){\r
1256         if( vY < fEtaMin || vY > fEtaMax)  continue;\r
1257       }\r
1258       \r
1259       //analyze one set of particles\r
1260       if(fUseMCPdgCode) {\r
1261         TParticle *particle = track->Particle();\r
1262         if(!particle) continue;\r
1263         \r
1264         Int_t gPdgCode = particle->GetPdgCode();\r
1265         if(TMath::Abs(fPDGCodeToBeAnalyzed) != TMath::Abs(gPdgCode)) \r
1266           continue;\r
1267       }\r
1268       \r
1269       //Use the acceptance parameterization\r
1270       if(fAcceptanceParameterization) {\r
1271         Double_t gRandomNumber = gRandom->Rndm();\r
1272         if(gRandomNumber > fAcceptanceParameterization->Eval(track->Pt())) \r
1273           continue;\r
1274       }\r
1275       \r
1276       //Exclude resonances\r
1277       if(fExcludeResonancesInMC) {\r
1278         TParticle *particle = track->Particle();\r
1279         if(!particle) continue;\r
1280         \r
1281         Bool_t kExcludeParticle = kFALSE;\r
1282         Int_t gMotherIndex = particle->GetFirstMother();\r
1283         if(gMotherIndex != -1) {\r
1284           AliMCParticle* motherTrack = dynamic_cast<AliMCParticle *>(event->GetTrack(gMotherIndex));\r
1285           if(motherTrack) {\r
1286             TParticle *motherParticle = motherTrack->Particle();\r
1287             if(motherParticle) {\r
1288               Int_t pdgCodeOfMother = motherParticle->GetPdgCode();\r
1289               //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {\r
1290               if(pdgCodeOfMother == 113) {\r
1291                 kExcludeParticle = kTRUE;\r
1292               }\r
1293             }\r
1294           }\r
1295         }\r
1296         \r
1297         //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi\r
1298         if(kExcludeParticle) continue;\r
1299       }\r
1300       \r
1301       vCharge = track->Charge();\r
1302       vPhi    = track->Phi();\r
1303       //Printf("phi (before): %lf",vPhi);\r
1304       \r
1305       fHistPt->Fill(vPt,gCentrality);\r
1306       fHistEta->Fill(vEta,gCentrality);\r
1307       fHistPhi->Fill(vPhi,gCentrality);\r
1308       //fHistPhi->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r
1309       fHistRapidity->Fill(vY,gCentrality);\r
1310       //if(vCharge > 0) fHistPhiPos->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r
1311       //else if(vCharge < 0) fHistPhiNeg->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r
1312       if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
1313       else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
1314       \r
1315       //Flow after burner\r
1316       if(fUseFlowAfterBurner) {\r
1317         Double_t precisionPhi = 0.001;\r
1318         Int_t maxNumberOfIterations = 100;\r
1319         \r
1320         Double_t phi0 = vPhi;\r
1321         Double_t gV2 = fDifferentialV2->Eval(vPt);\r
1322         \r
1323         for (Int_t j = 0; j < maxNumberOfIterations; j++) {\r
1324           Double_t phiprev = vPhi;\r
1325           Double_t fl = vPhi - phi0 + gV2*TMath::Sin(2.*(vPhi - gReactionPlane*TMath::DegToRad()));\r
1326           Double_t fp = 1.0 + 2.0*gV2*TMath::Cos(2.*(vPhi - gReactionPlane*TMath::DegToRad())); \r
1327           vPhi -= fl/fp;\r
1328           if (TMath::AreEqualAbs(phiprev,vPhi,precisionPhi)) break;\r
1329         }\r
1330         //Printf("phi (after): %lf\n",vPhi);\r
1331         Double_t vDeltaphiBefore = phi0 - gReactionPlane*TMath::DegToRad();\r
1332         if(vDeltaphiBefore < 0) vDeltaphiBefore += 2*TMath::Pi();\r
1333         fHistPhiBefore->Fill(vDeltaphiBefore,gCentrality);\r
1334         \r
1335         Double_t vDeltaphiAfter = vPhi - gReactionPlane*TMath::DegToRad();\r
1336         if(vDeltaphiAfter < 0) vDeltaphiAfter += 2*TMath::Pi();\r
1337         fHistPhiAfter->Fill(vDeltaphiAfter,gCentrality);\r
1338       \r
1339       }\r
1340       \r
1341       //vPhi *= TMath::RadToDeg();\r
1342 \r
1343       //=======================================correction\r
1344       Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  \r
1345       \r
1346       tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction)); \r
1347       \r
1348     } //track loop\r
1349   }//MC\r
1350   \r
1351   return tracksAccepted;  \r
1352 }\r
1353 //________________________________________________________________________\r
1354 TObjArray* AliAnalysisTaskBFPsi::GetShuffledTracks(TObjArray *tracks, Double_t gCentrality){\r
1355   // Clones TObjArray and returns it with tracks after shuffling the charges\r
1356 \r
1357   TObjArray* tracksShuffled = new TObjArray;\r
1358   tracksShuffled->SetOwner(kTRUE);\r
1359 \r
1360   vector<Short_t> *chargeVector = new vector<Short_t>;   //original charge of accepted tracks \r
1361 \r
1362   for (Int_t i=0; i<tracks->GetEntriesFast(); i++)\r
1363   {\r
1364     AliVParticle* track = (AliVParticle*) tracks->At(i);\r
1365     chargeVector->push_back(track->Charge());\r
1366   }  \r
1367  \r
1368   random_shuffle(chargeVector->begin(), chargeVector->end());\r
1369   \r
1370   for(Int_t i = 0; i < tracks->GetEntriesFast(); i++){\r
1371     AliVParticle* track = (AliVParticle*) tracks->At(i);\r
1372     //==============================correction\r
1373     Double_t correction = GetTrackbyTrackCorrectionMatrix(track->Eta(), track->Phi(),track->Pt(), chargeVector->at(i), gCentrality);\r
1374     tracksShuffled->Add(new AliBFBasicParticle(track->Eta(), track->Phi(), track->Pt(),chargeVector->at(i), correction));\r
1375   }\r
1376 \r
1377   delete chargeVector;\r
1378    \r
1379   return tracksShuffled;\r
1380 }\r
1381 \r
1382 \r
1383 //________________________________________________________________________\r
1384 void  AliAnalysisTaskBFPsi::FinishTaskOutput(){\r
1385   //Printf("END BF");\r
1386 \r
1387   if (!fBalance) {\r
1388     AliError("fBalance not available");\r
1389     return;\r
1390   }  \r
1391   if(fRunShuffling) {\r
1392     if (!fShuffledBalance) {\r
1393       AliError("fShuffledBalance not available");\r
1394       return;\r
1395     }\r
1396   }\r
1397 \r
1398 }\r
1399 \r
1400 //________________________________________________________________________\r
1401 void AliAnalysisTaskBFPsi::Terminate(Option_t *) {\r
1402   // Draw result to the screen\r
1403   // Called once at the end of the query\r
1404 \r
1405   // not implemented ...\r
1406 \r
1407 }\r