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