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