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.,7.,10.,20.,30.,40.,50.,60.,70.,80.,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     // multiplicity bins\r
 408     Double_t multiplicityBins[] = {0,10,20,30,40,50,60,70,80,100,100000}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
 409     Double_t* multbins        = multiplicityBins;\r
 410     Int_t nMultiplicityBins     = sizeof(multiplicityBins) / sizeof(Double_t) - 1;\r
 411     \r
 412     // Zvtx bins\r
 413     Double_t vertexBins[] = {-10., -7., -5., -3., -1., 1., 3., 5., 7., 10.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
 414     Double_t* vtxbins     = vertexBins;\r
 415     Int_t nVertexBins     = sizeof(vertexBins) / sizeof(Double_t) - 1;\r
 416     \r
 417     // Event plane angle (Psi) bins\r
 418     Double_t psiBins[] = {0.,45.,135.,215.,305.,360.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
 419     Double_t* psibins     = psiBins;\r
 420     Int_t nPsiBins     = sizeof(psiBins) / sizeof(Double_t) - 1;\r
 421     \r
 422     // run the event mixing also in bins of event plane (statistics!)\r
 423     if(fRunMixingEventPlane){\r
 424       if(fEventClass=="Multiplicity"){\r
 425         fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins, nPsiBins, psibins);\r
 426       }\r
 427       else{\r
 428         fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins, nPsiBins, psibins);\r
 429       }\r
 430     }\r
 431     else{\r
 432       if(fEventClass=="Multiplicity"){\r
 433         fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins);\r
 434       }\r
 435       else{\r
 436         fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);\r
 437       }\r
 438     }\r
 439   }\r
 440 \r
 441   if(fESDtrackCuts) fList->Add(fESDtrackCuts);\r
 442 \r
 443   //====================PID========================//\r
 444   if(fUsePID) {\r
 445     fPIDCombined = new AliPIDCombined();\r
 446     fPIDCombined->SetDefaultTPCPriors();\r
 447 \r
 448     fHistdEdxVsPTPCbeforePID = new TH2D ("dEdxVsPTPCbefore","dEdxVsPTPCbefore", 1000, -10.0, 10.0, 1000, 0, 1000); \r
 449     fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePID); //addition \r
 450     \r
 451     fHistBetavsPTOFbeforePID = new TH2D ("BetavsPTOFbefore","BetavsPTOFbefore", 1000, -10.0, 10., 1000, 0, 1.2); \r
 452     fHistListPIDQA->Add(fHistBetavsPTOFbeforePID); //addition\r
 453     \r
 454     fHistProbTPCvsPtbeforePID = new TH2D ("ProbTPCvsPtbefore","ProbTPCvsPtbefore", 1000, -10.0,10.0, 1000, 0, 2.0); \r
 455     fHistListPIDQA->Add(fHistProbTPCvsPtbeforePID); //addition \r
 456     \r
 457     fHistProbTOFvsPtbeforePID = new TH2D ("ProbTOFvsPtbefore","ProbTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0); \r
 458     fHistListPIDQA->Add(fHistProbTOFvsPtbeforePID); //addition \r
 459 \r
 460     fHistProbTPCTOFvsPtbeforePID =new TH2D ("ProbTPCTOFvsPtbefore","ProbTPCTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0); \r
 461     fHistListPIDQA->Add(fHistProbTPCTOFvsPtbeforePID); //addition \r
 462     \r
 463     fHistNSigmaTPCvsPtbeforePID = new TH2D ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore", 1000, -10, 10, 1000, 0, 500); \r
 464     fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePID); //addition \r
 465     \r
 466     fHistNSigmaTOFvsPtbeforePID = new TH2D ("NSigmaTOFvsPtbefore","NSigmaTOFvsPtbefore", 1000, -10, 10, 1000, 0, 500); \r
 467     fHistListPIDQA->Add(fHistNSigmaTOFvsPtbeforePID); //addition \r
 468     \r
 469     fHistdEdxVsPTPCafterPID = new TH2D ("dEdxVsPTPCafter","dEdxVsPTPCafter", 1000, -10, 10, 1000, 0, 1000); \r
 470     fHistListPIDQA->Add(fHistdEdxVsPTPCafterPID); //addition \r
 471     \r
 472     fHistBetavsPTOFafterPID = new TH2D ("BetavsPTOFafter","BetavsPTOFafter", 1000, -10, 10, 1000, 0, 1.2); \r
 473     fHistListPIDQA->Add(fHistBetavsPTOFafterPID); //addition \r
 474     \r
 475     fHistProbTPCvsPtafterPID = new TH2D ("ProbTPCvsPtafter","ProbTPCvsPtafter", 1000, -10, 10, 1000, 0, 2); \r
 476     fHistListPIDQA->Add(fHistProbTPCvsPtafterPID); //addition \r
 477   \r
 478     fHistProbTOFvsPtafterPID = new TH2D ("ProbTOFvsPtafter","ProbTOFvsPtafter", 1000,  -10, 10, 1000, 0, 2); \r
 479     fHistListPIDQA->Add(fHistProbTOFvsPtafterPID); //addition  \r
 480     \r
 481     fHistProbTPCTOFvsPtafterPID =new TH2D ("ProbTPCTOFvsPtafter","ProbTPCTOFvsPtafter", 1000, -50, 50, 1000, 0, 2.0); \r
 482     fHistListPIDQA->Add(fHistProbTPCTOFvsPtafterPID); //addition \r
 483 \r
 484     fHistNSigmaTPCvsPtafterPID = new TH2D ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter", 1000, -10, 10, 1000, 0, 500); \r
 485     fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPID); //addition  \r
 486     \r
 487     fHistNSigmaTOFvsPtafterPID = new TH2D ("NSigmaTOFvsPtafter","NSigmaTOFvsPtafter", 1000, -10, 10, 1000, 0, 500); \r
 488     fHistListPIDQA->Add(fHistNSigmaTOFvsPtafterPID); //addition \r
 489   }\r
 490   //====================PID========================//\r
 491 \r
 492   // Post output data.\r
 493   PostData(1, fList);\r
 494   PostData(2, fListBF);\r
 495   if(fRunShuffling) PostData(3, fListBFS);\r
 496   if(fRunMixing) PostData(4, fListBFM);\r
 497   if(fUsePID) PostData(5, fHistListPIDQA);       //PID\r
 498 \r
 499   TH1::AddDirectory(oldStatus);\r
 500 }\r
 501 \r
 502 \r
 503 //________________________________________________________________________\r
 504 void AliAnalysisTaskBFPsi::SetInputCorrection(const char* filename, \r
 505                                               const char* gCollSystem) {\r
 506   //Open files that will be used for correction\r
 507   TString gCollidingSystem = gCollSystem;\r
 508 \r
 509   //Open the input file\r
 510   TFile *f = TFile::Open(filename);\r
 511   if(!f->IsOpen()) {\r
 512     Printf("File not found!!!");\r
 513     return;\r
 514   }\r
 515   \r
 516   //Get the TDirectoryFile and TList\r
 517   TDirectoryFile *dir = dynamic_cast<TDirectoryFile *>(f->Get("PWGCFEbyE.outputBalanceFunctionEfficiencyAnalysis"));\r
 518   if(!dir) {\r
 519     Printf("TDirectoryFile not found!!!");\r
 520     return;\r
 521   }\r
 522   \r
 523   TString listEffName = "";\r
 524   for (Int_t iCent = 0; iCent < kCENTRALITY; iCent++) {\r
 525     listEffName = "listEfficiencyBF_";  \r
 526     listEffName += (TString)((Int_t)(centralityArrayForPbPb[iCent]));\r
 527     listEffName += "_";\r
 528     listEffName += (TString)((Int_t)(centralityArrayForPbPb[iCent+1]));\r
 529     TList *list = (TList*)dir->Get(listEffName.Data());\r
 530     if(!list) {\r
 531       Printf("TList Efficiency not found!!!");\r
 532       return;\r
 533     }\r
 534     \r
 535     TString histoName = "fHistMatrixCorrectionPlus";\r
 536     fHistMatrixCorrectionPlus[iCent]= dynamic_cast<TH3D *>(list->FindObject(histoName.Data()));\r
 537     if(!fHistMatrixCorrectionPlus[iCent]) {\r
 538       Printf("fHist not found!!!");\r
 539       return;\r
 540     }\r
 541     \r
 542     histoName = "fHistMatrixCorrectionMinus";\r
 543     fHistMatrixCorrectionMinus[iCent] = dynamic_cast<TH3D *>(list->FindObject(histoName.Data()));  \r
 544     if(!fHistMatrixCorrectionMinus[iCent]) {\r
 545       Printf("fHist not found!!!");\r
 546       return; \r
 547     }\r
 548   }//loop over centralities: ONLY the PbPb case is covered\r
 549 \r
 550   if(fHistMatrixCorrectionPlus[0]){\r
 551     fEtaMinForCorrections = fHistMatrixCorrectionPlus[0]->GetXaxis()->GetXmin();\r
 552     fEtaMaxForCorrections = fHistMatrixCorrectionPlus[0]->GetXaxis()->GetXmax();\r
 553     fEtaBinForCorrections = fHistMatrixCorrectionPlus[0]->GetNbinsX();\r
 554     \r
 555     fPtMinForCorrections = fHistMatrixCorrectionPlus[0]->GetYaxis()->GetXmin();\r
 556     fPtMaxForCorrections = fHistMatrixCorrectionPlus[0]->GetYaxis()->GetXmax();\r
 557     fPtBinForCorrections = fHistMatrixCorrectionPlus[0]->GetNbinsY();\r
 558     \r
 559     fPhiMinForCorrections = fHistMatrixCorrectionPlus[0]->GetZaxis()->GetXmin();\r
 560     fPhiMaxForCorrections = fHistMatrixCorrectionPlus[0]->GetZaxis()->GetXmax();\r
 561     fPhiBinForCorrections = fHistMatrixCorrectionPlus[0]->GetNbinsZ();\r
 562   }\r
 563 }\r
 564 \r
 565 //________________________________________________________________________\r
 566 void AliAnalysisTaskBFPsi::UserExec(Option_t *) {\r
 567   // Main loop\r
 568   // Called for each event\r
 569 \r
 570   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
 571   Int_t gNumberOfAcceptedTracks = 0;\r
 572   Double_t gCentrality          = -1.;\r
 573   Double_t gReactionPlane       = -1.; \r
 574   Float_t bSign = 0.;\r
 575   \r
 576   // get the event (for generator level: MCEvent())\r
 577   AliVEvent* eventMain = NULL;\r
 578   if(gAnalysisLevel == "MC") {\r
 579     eventMain = dynamic_cast<AliVEvent*>(MCEvent()); \r
 580   }\r
 581   else{\r
 582     eventMain = dynamic_cast<AliVEvent*>(InputEvent()); \r
 583     \r
 584     // for HBT like cuts need magnetic field sign\r
 585     bSign = (eventMain->GetMagneticField() > 0) ? 1 : -1;\r
 586   }\r
 587   if(!eventMain) {\r
 588     AliError("eventMain not available");\r
 589     return;\r
 590   }\r
 591   \r
 592   // PID Response task active?\r
 593   if(fUsePID) {\r
 594     fPIDResponse = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();\r
 595     if (!fPIDResponse) AliFatal("This Task needs the PID response attached to the inputHandler");\r
 596   }\r
 597   \r
 598   // check event cuts and fill event histograms\r
 599   if((gCentrality = IsEventAccepted(eventMain)) < 0){\r
 600     return;\r
 601   }\r
 602   \r
 603   //Compute Multiplicity or Centrality variable\r
 604   Double_t lMultiplicityVar = GetRefMultiOrCentrality( eventMain );\r
 605 \r
 606   // get the reaction plane\r
 607   if(fEventClass != "Multiplicity") {\r
 608     gReactionPlane = GetEventPlane(eventMain);\r
 609     fHistEventPlane->Fill(gReactionPlane,lMultiplicityVar);\r
 610     if(gReactionPlane < 0){\r
 611       return;\r
 612     }\r
 613   }\r
 614   \r
 615   // get the accepted tracks in main event\r
 616   TObjArray *tracksMain = GetAcceptedTracks(eventMain,lMultiplicityVar,gReactionPlane);\r
 617   gNumberOfAcceptedTracks = tracksMain->GetEntriesFast();\r
 618 \r
 619   //multiplicity cut (used in pp)\r
 620   fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedTracks,lMultiplicityVar);\r
 621 \r
 622   // store charges of all accepted tracks,shuffle and reassign(two extra loops)\r
 623   TObjArray* tracksShuffled = NULL;\r
 624   if(fRunShuffling){\r
 625     tracksShuffled = GetShuffledTracks(tracksMain,lMultiplicityVar);\r
 626   }\r
 627   \r
 628   // Event mixing \r
 629   if (fRunMixing)\r
 630     {\r
 631       // 1. First get an event pool corresponding in mult (cent) and\r
 632       //    zvertex to the current event. Once initialized, the pool\r
 633       //    should contain nMix (reduced) events. This routine does not\r
 634       //    pre-scan the chain. The first several events of every chain\r
 635       //    will be skipped until the needed pools are filled to the\r
 636       //    specified depth. If the pool categories are not too rare, this\r
 637       //    should not be a problem. If they are rare, you could lose`\r
 638       //    statistics.\r
 639       \r
 640       // 2. Collect the whole pool's content of tracks into one TObjArray\r
 641       //    (bgTracks), which is effectively a single background super-event.\r
 642       \r
 643       // 3. The reduced and bgTracks arrays must both be passed into\r
 644       //    FillCorrelations(). Also nMix should be passed in, so a weight\r
 645       //    of 1./nMix can be applied.\r
 646       \r
 647       AliEventPool* pool = fPoolMgr->GetEventPool(lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane);\r
 648       \r
 649       if (!pool){\r
 650         AliFatal(Form("No pool found for centrality = %f, zVtx = %f, psi = %f", lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane));\r
 651       }\r
 652       else{\r
 653         \r
 654         //pool->SetDebug(1);\r
 655 \r
 656         if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 || pool->GetCurrentNEvents() >= 5){ \r
 657           \r
 658           \r
 659           Int_t nMix = pool->GetCurrentNEvents();\r
 660           //cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << endl;\r
 661           \r
 662           //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(2);\r
 663           //((TH2F*) fListOfHistos->FindObject("mixedDist"))->Fill(centrality, pool->NTracksInPool());\r
 664           //if (pool->IsReady())\r
 665           //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(3);\r
 666           \r
 667           // Fill mixed-event histos here  \r
 668           for (Int_t jMix=0; jMix<nMix; jMix++) \r
 669             {\r
 670               TObjArray* tracksMixed = pool->GetEvent(jMix);\r
 671               fMixedBalance->CalculateBalance(gReactionPlane,tracksMain,tracksMixed,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());\r
 672             }\r
 673         }\r
 674         \r
 675         // Update the Event pool\r
 676         pool->UpdatePool(tracksMain);\r
 677         //pool->PrintInfo();\r
 678         \r
 679       }//pool NULL check  \r
 680     }//run mixing\r
 681   \r
 682   // calculate balance function\r
 683   fBalance->CalculateBalance(gReactionPlane,tracksMain,NULL,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());\r
 684   \r
 685   // calculate shuffled balance function\r
 686   if(fRunShuffling && tracksShuffled != NULL) {\r
 687     fShuffledBalance->CalculateBalance(gReactionPlane,tracksShuffled,NULL,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());\r
 688   }\r
 689 }      \r
 690 \r
 691 //________________________________________________________________________\r
 692 Double_t AliAnalysisTaskBFPsi::IsEventAccepted(AliVEvent *event){\r
 693   // Checks the Event cuts\r
 694   // Fills Event statistics histograms\r
 695   \r
 696   Bool_t isSelectedMain = kTRUE;\r
 697   Float_t gCentrality = -1.;\r
 698   Float_t gRefMultiplicity = -1.;\r
 699   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
 700 \r
 701   fHistEventStats->Fill(1,gCentrality); //all events\r
 702 \r
 703   // Event trigger bits\r
 704   fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());\r
 705   if(fUseOfflineTrigger)\r
 706     isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
 707   \r
 708   if(isSelectedMain) {\r
 709     fHistEventStats->Fill(2,gCentrality); //triggered events\r
 710     \r
 711     //Centrality stuff \r
 712     if(fUseCentrality) {\r
 713       if(gAnalysisLevel == "AOD") { //centrality in AOD header\r
 714         AliAODHeader *header = (AliAODHeader*) event->GetHeader();\r
 715         if(header){\r
 716           gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
 717           \r
 718           // QA for centrality estimators\r
 719           fHistCentStats->Fill(0.,header->GetCentralityP()->GetCentralityPercentile("V0M"));\r
 720           fHistCentStats->Fill(1.,header->GetCentralityP()->GetCentralityPercentile("FMD"));\r
 721           fHistCentStats->Fill(2.,header->GetCentralityP()->GetCentralityPercentile("TRK"));\r
 722           fHistCentStats->Fill(3.,header->GetCentralityP()->GetCentralityPercentile("TKL"));\r
 723           fHistCentStats->Fill(4.,header->GetCentralityP()->GetCentralityPercentile("CL0"));\r
 724           fHistCentStats->Fill(5.,header->GetCentralityP()->GetCentralityPercentile("CL1"));\r
 725           fHistCentStats->Fill(6.,header->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));\r
 726           fHistCentStats->Fill(7.,header->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));\r
 727           fHistCentStats->Fill(8.,header->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));\r
 728           \r
 729           // centrality QA (V0M)\r
 730           fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());\r
 731           \r
 732           // centrality QA (reference tracks)\r
 733           fHistRefTracks->Fill(0.,header->GetRefMultiplicity());\r
 734           fHistRefTracks->Fill(1.,header->GetRefMultiplicityPos());\r
 735           fHistRefTracks->Fill(2.,header->GetRefMultiplicityNeg());\r
 736           fHistRefTracks->Fill(3.,header->GetTPConlyRefMultiplicity());\r
 737           fHistRefTracks->Fill(4.,header->GetNumberOfITSClusters(0));\r
 738           fHistRefTracks->Fill(5.,header->GetNumberOfITSClusters(1));\r
 739           fHistRefTracks->Fill(6.,header->GetNumberOfITSClusters(2));\r
 740           fHistRefTracks->Fill(7.,header->GetNumberOfITSClusters(3));\r
 741           fHistRefTracks->Fill(8.,header->GetNumberOfITSClusters(4));\r
 742         }//AOD header\r
 743       }//AOD\r
 744 \r
 745       else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){ // centrality class for ESDs or MC-ESDs\r
 746         AliCentrality *centrality = event->GetCentrality();\r
 747         gCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
 748 \r
 749         // QA for centrality estimators\r
 750         fHistCentStats->Fill(0.,centrality->GetCentralityPercentile("V0M"));\r
 751         fHistCentStats->Fill(1.,centrality->GetCentralityPercentile("FMD"));\r
 752         fHistCentStats->Fill(2.,centrality->GetCentralityPercentile("TRK"));\r
 753         fHistCentStats->Fill(3.,centrality->GetCentralityPercentile("TKL"));\r
 754         fHistCentStats->Fill(4.,centrality->GetCentralityPercentile("CL0"));\r
 755         fHistCentStats->Fill(5.,centrality->GetCentralityPercentile("CL1"));\r
 756         fHistCentStats->Fill(6.,centrality->GetCentralityPercentile("V0MvsFMD"));\r
 757         fHistCentStats->Fill(7.,centrality->GetCentralityPercentile("TKLvsV0M"));\r
 758         fHistCentStats->Fill(8.,centrality->GetCentralityPercentile("ZEMvsZDC"));\r
 759 \r
 760         // centrality QA (V0M)\r
 761         fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());\r
 762       }//ESD\r
 763       else if(gAnalysisLevel == "MC"){\r
 764         Double_t gImpactParameter = 0.;\r
 765         if(dynamic_cast<AliMCEvent*>(event)){\r
 766           AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(dynamic_cast<AliMCEvent*>(event)->GenEventHeader());\r
 767           if(headerH){\r
 768             gImpactParameter = headerH->ImpactParameter();\r
 769             gCentrality      = gImpactParameter;\r
 770           }//MC header\r
 771         }//MC event cast\r
 772       }//MC\r
 773       else{\r
 774         gCentrality = -1.;\r
 775       }\r
 776     }\r
 777 \r
 778     //Multiplicity stuff \r
 779     if(fUseMultiplicity) \r
 780       gRefMultiplicity = GetRefMultiOrCentrality(event);\r
 781     \r
 782     // Event Vertex MC\r
 783     if(gAnalysisLevel == "MC"){\r
 784       if(!event) {\r
 785         AliError("mcEvent not available");\r
 786         return 0x0;\r
 787       }\r
 788       \r
 789       if(dynamic_cast<AliMCEvent*>(event)){\r
 790         AliGenEventHeader *header = dynamic_cast<AliGenEventHeader*>(dynamic_cast<AliMCEvent*>(event)->GenEventHeader());\r
 791         if(header){  \r
 792           TArrayF gVertexArray;\r
 793           header->PrimaryVertex(gVertexArray);\r
 794           //Printf("Vertex: %lf (x) - %lf (y) - %lf (z)",\r
 795           //gVertexArray.At(0),\r
 796           //gVertexArray.At(1),\r
 797           //gVertexArray.At(2));\r
 798           if(fUseMultiplicity) \r
 799             fHistEventStats->Fill(3,gRefMultiplicity); //events with a proper vertex\r
 800           else \r
 801             fHistEventStats->Fill(3,gCentrality); //events with a proper vertex\r
 802           if(TMath::Abs(gVertexArray.At(0)) < fVxMax) {\r
 803             if(TMath::Abs(gVertexArray.At(1)) < fVyMax) {\r
 804               if(TMath::Abs(gVertexArray.At(2)) < fVzMax) {\r
 805                 if(fUseMultiplicity) \r
 806                   fHistEventStats->Fill(4,gRefMultiplicity);//analyzed events\r
 807                 else \r
 808                   fHistEventStats->Fill(4,gCentrality); //analyzed events\r
 809                 fHistVx->Fill(gVertexArray.At(0));\r
 810                 fHistVy->Fill(gVertexArray.At(1));\r
 811                 fHistVz->Fill(gVertexArray.At(2),gCentrality);\r
 812                 \r
 813                 // take only events inside centrality class\r
 814                 if(fUseCentrality) {\r
 815                   if((fImpactParameterMin < gCentrality) && (fImpactParameterMax > gCentrality)){\r
 816                     fHistEventStats->Fill(5,gCentrality); //events with correct centrality\r
 817                     return gCentrality;     \r
 818                   }//centrality class\r
 819                 }\r
 820                 // take events only within the same multiplicity class\r
 821                 else if(fUseMultiplicity) {\r
 822                   if((gRefMultiplicity > fNumberOfAcceptedTracksMin)||(gRefMultiplicity < fNumberOfAcceptedTracksMax)) {\r
 823                     fHistEventStats->Fill(5,gRefMultiplicity); //events with correct multiplicity\r
 824                     return gRefMultiplicity;\r
 825                   }\r
 826                 }//multiplicity range\r
 827               }//Vz cut\r
 828             }//Vy cut\r
 829           }//Vx cut\r
 830         }//header    \r
 831       }//MC event object\r
 832     }//MC\r
 833     \r
 834     // Event Vertex AOD, ESD, ESDMC\r
 835     else{\r
 836       const AliVVertex *vertex = event->GetPrimaryVertex();\r
 837       \r
 838       if(vertex) {\r
 839         Double32_t fCov[6];\r
 840         vertex->GetCovarianceMatrix(fCov);\r
 841         if(vertex->GetNContributors() > 0) {\r
 842           if(fCov[5] != 0) {\r
 843             if(fUseMultiplicity) \r
 844               fHistEventStats->Fill(3,gRefMultiplicity); //proper vertex\r
 845             else \r
 846               fHistEventStats->Fill(3,gCentrality); //proper vertex\r
 847             if(TMath::Abs(vertex->GetX()) < fVxMax) {\r
 848               if(TMath::Abs(vertex->GetY()) < fVyMax) {\r
 849                 if(TMath::Abs(vertex->GetZ()) < fVzMax) {\r
 850                   if(fUseMultiplicity) \r
 851                     fHistEventStats->Fill(4,gRefMultiplicity);//analyzed events\r
 852                   else \r
 853                     fHistEventStats->Fill(4,gCentrality); //analyzed events\r
 854                   fHistVx->Fill(vertex->GetX());\r
 855                   fHistVy->Fill(vertex->GetY());\r
 856                   fHistVz->Fill(vertex->GetZ(),gCentrality);\r
 857                   \r
 858                   // take only events inside centrality class\r
 859                   if(fUseCentrality) {\r
 860                     if((gCentrality > fCentralityPercentileMin) && (gCentrality < fCentralityPercentileMax)){\r
 861                       fHistEventStats->Fill(5,gCentrality); //events with correct centrality\r
 862                       return gCentrality;               \r
 863                     }//centrality class\r
 864                   }\r
 865                   // take events only within the same multiplicity class\r
 866                   else if(fUseMultiplicity) {\r
 867                     if((gRefMultiplicity > fNumberOfAcceptedTracksMin)||(gRefMultiplicity < fNumberOfAcceptedTracksMax)) {\r
 868                       fHistEventStats->Fill(5,gRefMultiplicity); //events with correct multiplicity\r
 869                       return gRefMultiplicity;\r
 870                     }\r
 871                   }//multiplicity range\r
 872                 }//Vz cut\r
 873               }//Vy cut\r
 874             }//Vx cut\r
 875           }//proper vertex resolution\r
 876         }//proper number of contributors\r
 877       }//vertex object valid\r
 878     }//triggered event \r
 879   }//AOD,ESD,ESDMC\r
 880   \r
 881   // in all other cases return -1 (event not accepted)\r
 882   return -1;\r
 883 }\r
 884 \r
 885 \r
 886 //________________________________________________________________________\r
 887 Double_t AliAnalysisTaskBFPsi::GetRefMultiOrCentrality(AliVEvent *event){\r
 888     // Checks the Event cuts\r
 889     // Fills Event statistics histograms\r
 890   \r
 891   Float_t gCentrality = -1.;\r
 892   Double_t fMultiplicity = -100.;\r
 893   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
 894 \r
 895   if(fEventClass == "Centrality"){\r
 896     if(gAnalysisLevel == "AOD") { //centrality in AOD header\r
 897       AliAODHeader *header = (AliAODHeader*) event->GetHeader();\r
 898       if(header){\r
 899         gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
 900       }//AOD header\r
 901     }//AOD\r
 902     \r
 903     else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){ // centrality class for ESDs or MC-ESDs\r
 904       AliCentrality *centrality = event->GetCentrality();\r
 905       gCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
 906     }//ESD\r
 907     else if(gAnalysisLevel == "MC"){\r
 908       Double_t gImpactParameter = 0.;\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           gImpactParameter = headerH->ImpactParameter();\r
 913           gCentrality      = gImpactParameter;\r
 914         }//MC header\r
 915       }//MC event cast\r
 916     }//MC\r
 917     else{\r
 918       gCentrality = -1.;\r
 919     }\r
 920   }//End if "Centrality"\r
 921   if(fEventClass=="Multiplicity"&&gAnalysisLevel == "ESD"){\r
 922     AliESDEvent* gESDEvent = dynamic_cast<AliESDEvent*>(event);\r
 923     if(gESDEvent){\r
 924       fMultiplicity = fESDtrackCuts->GetReferenceMultiplicity(gESDEvent, AliESDtrackCuts::kTrackletsITSTPC,0.5);\r
 925     }//AliESDevent cast\r
 926   }\r
 927   else if(fEventClass=="Multiplicity"&&gAnalysisLevel == "AOD"){\r
 928     AliAODHeader *header = (AliAODHeader*) event->GetHeader();\r
 929     if(header){\r
 930       fMultiplicity = header->GetRefMultiplicity();\r
 931     }//AOD header\r
 932   }\r
 933   else if(fEventClass=="Multiplicity"&&gAnalysisLevel == "MC") {\r
 934     AliMCEvent* gMCEvent = dynamic_cast<AliMCEvent*>(event);\r
 935     fMultiplicity = gMCEvent->GetNumberOfPrimaries();\r
 936   }\r
 937 \r
 938   Double_t lReturnVal = -100;\r
 939   if(fEventClass=="Multiplicity"){\r
 940     lReturnVal = fMultiplicity;\r
 941   }else if(fEventClass=="Centrality"){\r
 942     lReturnVal = gCentrality;\r
 943   }\r
 944   return lReturnVal;\r
 945 }\r
 946 \r
 947 //________________________________________________________________________\r
 948 Double_t AliAnalysisTaskBFPsi::GetEventPlane(AliVEvent *event){\r
 949   // Get the event plane\r
 950 \r
 951   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
 952 \r
 953   Float_t gVZEROEventPlane    = -10.;\r
 954   Float_t gReactionPlane      = -10.;\r
 955   Double_t qxTot = 0.0, qyTot = 0.0;\r
 956 \r
 957   //MC: from reaction plane\r
 958   if(gAnalysisLevel == "MC"){\r
 959     if(!event) {\r
 960       AliError("mcEvent not available");\r
 961       return 0x0;\r
 962     }\r
 963 \r
 964     AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);\r
 965     if(gMCEvent){\r
 966       AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader());    \r
 967       if (headerH) {\r
 968         gReactionPlane = headerH->ReactionPlaneAngle();\r
 969         //gReactionPlane *= TMath::RadToDeg();\r
 970       }//MC header\r
 971     }//MC event cast\r
 972   }//MC\r
 973   \r
 974   // AOD,ESD,ESDMC: from VZERO Event Plane\r
 975   else{\r
 976    \r
 977     AliEventplane *ep = event->GetEventplane();\r
 978     if(ep){ \r
 979       gVZEROEventPlane = ep->CalculateVZEROEventPlane(event,10,2,qxTot,qyTot);\r
 980       if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();\r
 981       //gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();\r
 982       gReactionPlane = gVZEROEventPlane;\r
 983     }\r
 984   }//AOD,ESD,ESDMC\r
 985 \r
 986   return gReactionPlane;\r
 987 }\r
 988 \r
 989 //________________________________________________________________________\r
 990 Double_t AliAnalysisTaskBFPsi::GetTrackbyTrackCorrectionMatrix( Double_t vEta, \r
 991                                                                 Double_t vPhi, \r
 992                                                                 Double_t vPt, \r
 993                                                                 Short_t vCharge, \r
 994                                                                 Double_t gCentrality) {\r
 995   // -- Get efficiency correction of particle dependent on (eta, phi, pt, charge, centrality) \r
 996 \r
 997   Double_t correction = 1.;\r
 998   Int_t binEta = 0, binPt = 0, binPhi = 0;\r
 999 \r
1000   //Printf("EtaMAx: %lf - EtaMin: %lf - EtaBin: %lf", fEtaMaxForCorrections,fEtaMinForCorrections,fEtaBinForCorrections);\r
1001   if(fEtaBinForCorrections != 0) {\r
1002     Double_t widthEta = (fEtaMaxForCorrections - fEtaMinForCorrections)/fEtaBinForCorrections;\r
1003     if(fEtaMaxForCorrections != fEtaMinForCorrections) \r
1004       binEta = (Int_t)(vEta/widthEta)+1;\r
1005   }\r
1006 \r
1007   if(fPtBinForCorrections != 0) {\r
1008     Double_t widthPt = (fPtMaxForCorrections - fPtMinForCorrections)/fPtBinForCorrections;\r
1009     if(fPtMaxForCorrections != fPtMinForCorrections) \r
1010       binPt = (Int_t)(vPt/widthPt) + 1;\r
1011   }\r
1012  \r
1013   if(fPhiBinForCorrections != 0) {\r
1014     Double_t widthPhi = (fPhiMaxForCorrections - fPhiMinForCorrections)/fPhiBinForCorrections;\r
1015     if(fPhiMaxForCorrections != fPhiMinForCorrections) \r
1016       binPhi = (Int_t)(vPhi/widthPhi)+ 1;\r
1017   }\r
1018 \r
1019   Int_t gCentralityInt = 1;\r
1020   for (Int_t i=0; i<kCENTRALITY; i++){\r
1021     if((centralityArrayForPbPb[i] <= gCentrality)&&(gCentrality <= centralityArrayForPbPb[i+1]))\r
1022       gCentralityInt = i;\r
1023   }\r
1024   \r
1025   if(fHistMatrixCorrectionPlus[gCentralityInt]){\r
1026     if (vCharge > 0) {\r
1027       correction = fHistMatrixCorrectionPlus[gCentralityInt]->GetBinContent(fHistMatrixCorrectionPlus[gCentralityInt]->GetBin(binEta, binPt, binPhi));\r
1028     }\r
1029     if (vCharge < 0) {\r
1030       correction = fHistMatrixCorrectionMinus[gCentralityInt]->GetBinContent(fHistMatrixCorrectionMinus[gCentralityInt]->GetBin(binEta, binPt, binPhi));\r
1031     }\r
1032   }\r
1033   else {\r
1034     correction = 1.;\r
1035   }\r
1036   if (correction == 0.) { \r
1037     AliError(Form("Should not happen : bin content = 0. >> eta: %.2f | phi : %.2f | pt : %.2f | cent %d",vEta, vPhi, vPt, gCentralityInt)); \r
1038     correction = 1.; \r
1039   } \r
1040     \r
1041   return correction;\r
1042 }\r
1043 \r
1044 //________________________________________________________________________\r
1045 TObjArray* AliAnalysisTaskBFPsi::GetAcceptedTracks(AliVEvent *event, Double_t gCentrality, Double_t gReactionPlane){\r
1046   // Returns TObjArray with tracks after all track cuts (only for AOD!)\r
1047   // Fills QA histograms\r
1048 \r
1049   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
1050 \r
1051   //output TObjArray holding all good tracks\r
1052   TObjArray* tracksAccepted = new TObjArray;\r
1053   tracksAccepted->SetOwner(kTRUE);\r
1054 \r
1055   Short_t vCharge;\r
1056   Double_t vEta;\r
1057   Double_t vY;\r
1058   Double_t vPhi;\r
1059   Double_t vPt;\r
1060 \r
1061 \r
1062   if(gAnalysisLevel == "AOD") { // handling of TPC only tracks different in AOD and ESD\r
1063     // Loop over tracks in event\r
1064     for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r
1065       AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));\r
1066       if (!aodTrack) {\r
1067         AliError(Form("Could not receive track %d", iTracks));\r
1068         continue;\r
1069       }\r
1070       \r
1071       // AOD track cuts\r
1072       \r
1073       // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C\r
1074       // take only TPC only tracks \r
1075       //fHistTrackStats->Fill(aodTrack->GetFilterMap());\r
1076       for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){\r
1077         fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));\r
1078       }\r
1079       if(!aodTrack->TestFilterBit(nAODtrackCutBit)) continue;\r
1080       \r
1081       vCharge = aodTrack->Charge();\r
1082       vEta    = aodTrack->Eta();\r
1083       vY      = aodTrack->Y();\r
1084       vPhi    = aodTrack->Phi();// * TMath::RadToDeg();\r
1085       vPt     = aodTrack->Pt();\r
1086       \r
1087       Float_t dcaXY = aodTrack->DCA();      // this is the DCA from global track (not exactly what is cut on)\r
1088       Float_t dcaZ  = aodTrack->ZAtDCA();   // this is the DCA from global track (not exactly what is cut on)\r
1089       \r
1090       \r
1091       // Kinematics cuts from ESD track cuts\r
1092       if( vPt < fPtMin || vPt > fPtMax)      continue;\r
1093       if( vEta < fEtaMin || vEta > fEtaMax)  continue;\r
1094       \r
1095       // Extra DCA cuts (for systematic studies [!= -1])\r
1096       if( fDCAxyCut != -1 && fDCAzCut != -1){\r
1097         if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){\r
1098           continue;  // 2D cut\r
1099         }\r
1100       }\r
1101       \r
1102       // Extra TPC cuts (for systematic studies [!= -1])\r
1103       if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){\r
1104         continue;\r
1105       }\r
1106       if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){\r
1107         continue;\r
1108       }\r
1109       \r
1110       // fill QA histograms\r
1111       fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());\r
1112       fHistDCA->Fill(dcaZ,dcaXY);\r
1113       fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality);\r
1114       fHistPt->Fill(vPt,gCentrality);\r
1115       fHistEta->Fill(vEta,gCentrality);\r
1116       fHistRapidity->Fill(vY,gCentrality);\r
1117       if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
1118       else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
1119       fHistPhi->Fill(vPhi,gCentrality);\r
1120       if(vCharge > 0)      fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);                  \r
1121       else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r
1122       \r
1123       //=======================================correction\r
1124       Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  \r
1125       \r
1126       // add the track to the TObjArray\r
1127       tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));  \r
1128     }//track loop\r
1129   }// AOD analysis\r
1130 \r
1131 \r
1132   else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD") { // handling of TPC only tracks different in AOD and ESD\r
1133 \r
1134     AliESDtrack *trackTPC   = NULL;\r
1135     AliMCParticle *mcTrack   = NULL;\r
1136 \r
1137     AliMCEvent*  mcEvent     = NULL;\r
1138     \r
1139     // for MC ESDs use also MC information\r
1140     if(gAnalysisLevel == "MCESD"){\r
1141       mcEvent = MCEvent(); \r
1142       if (!mcEvent) {\r
1143         AliError("mcEvent not available");\r
1144         return tracksAccepted;\r
1145       }\r
1146     }\r
1147     \r
1148     // Loop over tracks in event\r
1149     for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r
1150       AliESDtrack* track = dynamic_cast<AliESDtrack *>(event->GetTrack(iTracks));\r
1151       if (!track) {\r
1152         AliError(Form("Could not receive track %d", iTracks));\r
1153         continue;\r
1154       } \r
1155 \r
1156       // for MC ESDs use also MC information --> MC track not used further???\r
1157       if(gAnalysisLevel == "MCESD"){\r
1158         Int_t label = TMath::Abs(track->GetLabel());\r
1159         if(label > mcEvent->GetNumberOfTracks()) continue;\r
1160         if(label > mcEvent->GetNumberOfPrimaries()) continue;\r
1161         \r
1162         mcTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label));\r
1163         if(!mcTrack) continue;\r
1164       }\r
1165 \r
1166       // take only TPC only tracks\r
1167       trackTPC   = new AliESDtrack();\r
1168       if(!track->FillTPCOnlyTrack(*trackTPC)) continue;\r
1169       \r
1170       //ESD track cuts\r
1171       if(fESDtrackCuts) \r
1172         if(!fESDtrackCuts->AcceptTrack(trackTPC)) continue;\r
1173       \r
1174       // fill QA histograms\r
1175       Float_t b[2];\r
1176       Float_t bCov[3];\r
1177       trackTPC->GetImpactParameters(b,bCov);\r
1178       if (bCov[0]<=0 || bCov[2]<=0) {\r
1179         AliDebug(1, "Estimated b resolution lower or equal zero!");\r
1180         bCov[0]=0; bCov[2]=0;\r
1181       }\r
1182       \r
1183       Int_t nClustersTPC = -1;\r
1184       nClustersTPC = trackTPC->GetTPCNclsIter1();   // TPC standalone\r
1185       //nClustersTPC = track->GetTPCclusters(0);   // global track\r
1186       Float_t chi2PerClusterTPC = -1;\r
1187       if (nClustersTPC!=0) {\r
1188         chi2PerClusterTPC = trackTPC->GetTPCchi2Iter1()/Float_t(nClustersTPC);      // TPC standalone\r
1189         //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);     // global track\r
1190       }\r
1191       \r
1192       //===========================PID===============================//             \r
1193       if(fUsePID) {\r
1194         Double_t prob[AliPID::kSPECIES]={0.};\r
1195         Double_t probTPC[AliPID::kSPECIES]={0.};\r
1196         Double_t probTOF[AliPID::kSPECIES]={0.};\r
1197         Double_t probTPCTOF[AliPID::kSPECIES]={0.};\r
1198         \r
1199         Double_t nSigma = 0.;\r
1200         UInt_t detUsedTPC = 0;\r
1201         UInt_t detUsedTOF = 0;\r
1202         UInt_t detUsedTPCTOF = 0;\r
1203         \r
1204         //Decide what detector configuration we want to use\r
1205         switch(fPidDetectorConfig) {\r
1206         case kTPCpid:\r
1207           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);\r
1208           nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest));\r
1209           detUsedTPC = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPC);\r
1210           for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
1211             prob[iSpecies] = probTPC[iSpecies];\r
1212           break;\r
1213         case kTOFpid:\r
1214           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);\r
1215           nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest));\r
1216           detUsedTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTOF);\r
1217           for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
1218             prob[iSpecies] = probTOF[iSpecies];\r
1219           break;\r
1220         case kTPCTOF:\r
1221           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);\r
1222           detUsedTPCTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPCTOF);\r
1223           for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
1224             prob[iSpecies] = probTPCTOF[iSpecies];\r
1225           break;\r
1226         default:\r
1227           break;\r
1228         }//end switch: define detector mask\r
1229         \r
1230         //Filling the PID QA\r
1231         Double_t tofTime = -999., length = 999., tof = -999.;\r
1232         Double_t c = TMath::C()*1.E-9;// m/ns\r
1233         Double_t beta = -999.;\r
1234         Double_t  nSigmaTOFForParticleOfInterest = -999.;\r
1235         if ( (track->IsOn(AliESDtrack::kTOFin)) &&\r
1236              (track->IsOn(AliESDtrack::kTIME))  ) { \r
1237           tofTime = track->GetTOFsignal();//in ps\r
1238           length = track->GetIntegratedLength();\r
1239           tof = tofTime*1E-3; // ns     \r
1240           \r
1241           if (tof <= 0) {\r
1242             //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");\r
1243             continue;\r
1244           }\r
1245           if (length <= 0){\r
1246             //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");\r
1247             continue;\r
1248           }\r
1249           \r
1250           length = length*0.01; // in meters\r
1251           tof = tof*c;\r
1252           beta = length/tof;\r
1253           \r
1254           nSigmaTOFForParticleOfInterest = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest);\r
1255           fHistBetavsPTOFbeforePID ->Fill(track->P()*track->Charge(),beta);\r
1256           fHistProbTOFvsPtbeforePID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);\r
1257           fHistNSigmaTOFvsPtbeforePID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);\r
1258         }//TOF signal \r
1259         \r
1260         \r
1261         Double_t  nSigmaTPCForParticleOfInterest = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest);\r
1262         fHistdEdxVsPTPCbeforePID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());\r
1263         fHistProbTPCvsPtbeforePID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); \r
1264         fHistNSigmaTPCvsPtbeforePID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); \r
1265         fHistProbTPCTOFvsPtbeforePID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);\r
1266         //end of QA-before pid\r
1267         \r
1268         if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {\r
1269           //Make the decision based on the n-sigma\r
1270           if(fUsePIDnSigma) {\r
1271             if(nSigma > fPIDNSigma) continue;}\r
1272           \r
1273           //Make the decision based on the bayesian\r
1274           else if(fUsePIDPropabilities) {\r
1275             if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;\r
1276             if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue;      \r
1277           }\r
1278           \r
1279           //Fill QA after the PID\r
1280           fHistBetavsPTOFafterPID ->Fill(track->P()*track->Charge(),beta);\r
1281           fHistProbTOFvsPtafterPID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);\r
1282           fHistNSigmaTOFvsPtafterPID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);\r
1283           \r
1284           fHistdEdxVsPTPCafterPID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());\r
1285           fHistProbTPCvsPtafterPID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); \r
1286           fHistProbTPCTOFvsPtafterPID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);\r
1287           fHistNSigmaTPCvsPtafterPID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); \r
1288         }\r
1289       }\r
1290       //===========================PID===============================//\r
1291       vCharge = trackTPC->Charge();\r
1292       vY      = trackTPC->Y();\r
1293       vEta    = trackTPC->Eta();\r
1294       vPhi    = trackTPC->Phi();// * TMath::RadToDeg();\r
1295       vPt     = trackTPC->Pt();\r
1296       fHistClus->Fill(trackTPC->GetITSclusters(0),nClustersTPC);\r
1297       fHistDCA->Fill(b[1],b[0]);\r
1298       fHistChi2->Fill(chi2PerClusterTPC,gCentrality);\r
1299       fHistPt->Fill(vPt,gCentrality);\r
1300       fHistEta->Fill(vEta,gCentrality);\r
1301       fHistPhi->Fill(vPhi,gCentrality);\r
1302       if(vCharge > 0)      fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);\r
1303       else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r
1304       fHistRapidity->Fill(vY,gCentrality);\r
1305       if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
1306       else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
1307       \r
1308       //=======================================correction\r
1309       Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  \r
1310 \r
1311       // add the track to the TObjArray\r
1312       tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));   \r
1313 \r
1314       delete trackTPC;\r
1315     }//track loop\r
1316   }// ESD analysis\r
1317 \r
1318   else if(gAnalysisLevel == "MC"){\r
1319     if(!event) {\r
1320       AliError("mcEvent not available");\r
1321       return 0x0;\r
1322     }\r
1323 \r
1324     AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);\r
1325     if(gMCEvent) {\r
1326       // Loop over tracks in event\r
1327       for (Int_t iTracks = 0; iTracks < gMCEvent->GetNumberOfPrimaries(); iTracks++) {\r
1328         AliMCParticle* track = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iTracks));\r
1329         if (!track) {\r
1330           AliError(Form("Could not receive particle %d", iTracks));\r
1331           continue;\r
1332         }\r
1333         \r
1334         //exclude non stable particles\r
1335         if(!(gMCEvent->IsPhysicalPrimary(iTracks))) continue;\r
1336         \r
1337         vEta    = track->Eta();\r
1338         vPt     = track->Pt();\r
1339         vY      = track->Y();\r
1340         \r
1341         if( vPt < fPtMin || vPt > fPtMax)      \r
1342           continue;\r
1343         if (!fUsePID) {\r
1344           if( vEta < fEtaMin || vEta > fEtaMax)  continue;\r
1345         }\r
1346         else if (fUsePID){\r
1347           if( vY < fEtaMin || vY > fEtaMax)  continue;\r
1348         }\r
1349         \r
1350         //analyze one set of particles\r
1351         if(fUseMCPdgCode) {\r
1352           TParticle *particle = track->Particle();\r
1353           if(!particle) continue;\r
1354           \r
1355           Int_t gPdgCode = particle->GetPdgCode();\r
1356           if(TMath::Abs(fPDGCodeToBeAnalyzed) != TMath::Abs(gPdgCode)) \r
1357             continue;\r
1358         }\r
1359         \r
1360         //Use the acceptance parameterization\r
1361         if(fAcceptanceParameterization) {\r
1362           Double_t gRandomNumber = gRandom->Rndm();\r
1363           if(gRandomNumber > fAcceptanceParameterization->Eval(track->Pt())) \r
1364             continue;\r
1365         }\r
1366         \r
1367         //Exclude resonances\r
1368         if(fExcludeResonancesInMC) {\r
1369           TParticle *particle = track->Particle();\r
1370           if(!particle) continue;\r
1371           \r
1372           Bool_t kExcludeParticle = kFALSE;\r
1373           Int_t gMotherIndex = particle->GetFirstMother();\r
1374           if(gMotherIndex != -1) {\r
1375             AliMCParticle* motherTrack = dynamic_cast<AliMCParticle *>(event->GetTrack(gMotherIndex));\r
1376             if(motherTrack) {\r
1377               TParticle *motherParticle = motherTrack->Particle();\r
1378               if(motherParticle) {\r
1379                 Int_t pdgCodeOfMother = motherParticle->GetPdgCode();\r
1380                 //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {\r
1381                 if(pdgCodeOfMother == 113) {\r
1382                   kExcludeParticle = kTRUE;\r
1383                 }\r
1384               }\r
1385             }\r
1386           }\r
1387           \r
1388           //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi\r
1389           if(kExcludeParticle) continue;\r
1390         }\r
1391         \r
1392         vCharge = track->Charge();\r
1393         vPhi    = track->Phi();\r
1394         //Printf("phi (before): %lf",vPhi);\r
1395         \r
1396         fHistPt->Fill(vPt,gCentrality);\r
1397         fHistEta->Fill(vEta,gCentrality);\r
1398         fHistPhi->Fill(vPhi,gCentrality);\r
1399         if(vCharge > 0)      fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);\r
1400         else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r
1401         //fHistPhi->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r
1402         fHistRapidity->Fill(vY,gCentrality);\r
1403         //if(vCharge > 0) fHistPhiPos->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r
1404         //else if(vCharge < 0) fHistPhiNeg->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r
1405         if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
1406         else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
1407         \r
1408         //Flow after burner\r
1409         if(fUseFlowAfterBurner) {\r
1410           Double_t precisionPhi = 0.001;\r
1411           Int_t maxNumberOfIterations = 100;\r
1412           \r
1413           Double_t phi0 = vPhi;\r
1414           Double_t gV2 = fDifferentialV2->Eval(vPt);\r
1415           \r
1416           for (Int_t j = 0; j < maxNumberOfIterations; j++) {\r
1417             Double_t phiprev = vPhi;\r
1418             Double_t fl = vPhi - phi0 + gV2*TMath::Sin(2.*(vPhi - gReactionPlane*TMath::DegToRad()));\r
1419             Double_t fp = 1.0 + 2.0*gV2*TMath::Cos(2.*(vPhi - gReactionPlane*TMath::DegToRad())); \r
1420             vPhi -= fl/fp;\r
1421             if (TMath::AreEqualAbs(phiprev,vPhi,precisionPhi)) break;\r
1422           }\r
1423           //Printf("phi (after): %lf\n",vPhi);\r
1424           Double_t vDeltaphiBefore = phi0 - gReactionPlane*TMath::DegToRad();\r
1425           if(vDeltaphiBefore < 0) vDeltaphiBefore += 2*TMath::Pi();\r
1426           fHistPhiBefore->Fill(vDeltaphiBefore,gCentrality);\r
1427           \r
1428           Double_t vDeltaphiAfter = vPhi - gReactionPlane*TMath::DegToRad();\r
1429           if(vDeltaphiAfter < 0) vDeltaphiAfter += 2*TMath::Pi();\r
1430           fHistPhiAfter->Fill(vDeltaphiAfter,gCentrality);\r
1431           \r
1432         }\r
1433         \r
1434         //vPhi *= TMath::RadToDeg();\r
1435         \r
1436         //=======================================correction\r
1437         Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge,gCentrality);  \r
1438 \r
1439         tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction)); \r
1440       } //track loop\r
1441     }//MC event object\r
1442   }//MC\r
1443   \r
1444   return tracksAccepted;  \r
1445 }\r
1446 \r
1447 //________________________________________________________________________\r
1448 TObjArray* AliAnalysisTaskBFPsi::GetShuffledTracks(TObjArray *tracks, Double_t gCentrality){\r
1449   // Clones TObjArray and returns it with tracks after shuffling the charges\r
1450 \r
1451   TObjArray* tracksShuffled = new TObjArray;\r
1452   tracksShuffled->SetOwner(kTRUE);\r
1453 \r
1454   vector<Short_t> *chargeVector = new vector<Short_t>;   //original charge of accepted tracks \r
1455 \r
1456   for (Int_t i=0; i<tracks->GetEntriesFast(); i++)\r
1457   {\r
1458     AliVParticle* track = (AliVParticle*) tracks->At(i);\r
1459     chargeVector->push_back(track->Charge());\r
1460   }  \r
1461  \r
1462   random_shuffle(chargeVector->begin(), chargeVector->end());\r
1463   \r
1464   for(Int_t i = 0; i < tracks->GetEntriesFast(); i++){\r
1465     AliVParticle* track = (AliVParticle*) tracks->At(i);\r
1466     //==============================correction\r
1467     Double_t correction = GetTrackbyTrackCorrectionMatrix(track->Eta(), track->Phi(),track->Pt(), chargeVector->at(i), gCentrality);\r
1468     tracksShuffled->Add(new AliBFBasicParticle(track->Eta(), track->Phi(), track->Pt(),chargeVector->at(i), correction));\r
1469   }\r
1470 \r
1471   delete chargeVector;\r
1472    \r
1473   return tracksShuffled;\r
1474 }\r
1475 \r
1476 \r
1477 //________________________________________________________________________\r
1478 void  AliAnalysisTaskBFPsi::FinishTaskOutput(){\r
1479   //Printf("END BF");\r
1480 \r
1481   if (!fBalance) {\r
1482     AliError("fBalance not available");\r
1483     return;\r
1484   }  \r
1485   if(fRunShuffling) {\r
1486     if (!fShuffledBalance) {\r
1487       AliError("fShuffledBalance not available");\r
1488       return;\r
1489     }\r
1490   }\r
1491 \r
1492 }\r
1493 \r
1494 //________________________________________________________________________\r
1495 void AliAnalysisTaskBFPsi::Terminate(Option_t *) {\r
1496   // Draw result to the screen\r
1497   // Called once at the end of the query\r
1498 \r
1499   // not implemented ...\r
1500 \r
1501 }\r