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 "AliAODMCParticle.h" \r
  25 #include "AliCollisionGeometry.h"\r
  26 #include "AliGenEventHeader.h"\r
  27 #include "AliMCEventHandler.h"\r
  28 #include "AliMCEvent.h"\r
  29 #include "AliStack.h"\r
  30 #include "AliESDtrackCuts.h"\r
  31 #include "AliEventplane.h"\r
  32 #include "AliTHn.h"    \r
  33 #include "AliLog.h"\r
  34 #include "AliAnalysisUtils.h"\r
  35 \r
  36 #include "AliEventPoolManager.h"           \r
  37 \r
  38 #include "AliPID.h"                \r
  39 #include "AliPIDResponse.h"        \r
  40 #include "AliPIDCombined.h"        \r
  41 \r
  42 #include "AliAnalysisTaskBFPsi.h"\r
  43 #include "AliBalancePsi.h"\r
  44 #include "AliAnalysisTaskTriggeredBF.h"\r
  45 \r
  46 \r
  47 // Analysis task for the BF vs Psi code\r
  48 // Authors: Panos.Christakoglou@nikhef.nl\r
  49 \r
  50 using std::cout;\r
  51 using std::endl;\r
  52 \r
  53 ClassImp(AliAnalysisTaskBFPsi)\r
  54 \r
  55 //________________________________________________________________________\r
  56 AliAnalysisTaskBFPsi::AliAnalysisTaskBFPsi(const char *name) \r
  57 : AliAnalysisTaskSE(name),\r
  58   fDebugLevel(kFALSE),\r
  59   fArrayMC(0), //+++++++++++++\r
  60   fBalance(0),\r
  61   fRunShuffling(kFALSE),\r
  62   fShuffledBalance(0),\r
  63   fRunMixing(kFALSE),\r
  64   fRunMixingEventPlane(kFALSE),\r
  65   fMixingTracks(50000),\r
  66   fMixedBalance(0),\r
  67   fPoolMgr(0),\r
  68   fList(0),\r
  69   fListBF(0),\r
  70   fListBFS(0),\r
  71   fListBFM(0),\r
  72   fHistListPIDQA(0),\r
  73   fHistEventStats(0),\r
  74   fHistCentStats(0),\r
  75   fHistCentStatsUsed(0),\r
  76   fHistTriggerStats(0),\r
  77   fHistTrackStats(0),\r
  78   fHistVx(0),\r
  79   fHistVy(0),\r
  80   fHistVz(0),\r
  81   fHistTPCvsVZEROMultiplicity(0),\r
  82   fHistVZEROSignal(0),\r
  83   fHistEventPlane(0),\r
  84   fHistClus(0),\r
  85   fHistDCA(0),\r
  86   fHistChi2(0),\r
  87   fHistPt(0),\r
  88   fHistEta(0),\r
  89   fHistRapidity(0),\r
  90   fHistPhi(0),\r
  91   fHistEtaPhiPos(0),             \r
  92   fHistEtaPhiNeg(0), \r
  93   fHistPhiBefore(0),\r
  94   fHistPhiAfter(0),\r
  95   fHistPhiPos(0),\r
  96   fHistPhiNeg(0),\r
  97   fHistV0M(0),\r
  98   fHistRefTracks(0),\r
  99   fHistdEdxVsPTPCbeforePID(NULL),\r
 100   fHistBetavsPTOFbeforePID(NULL), \r
 101   fHistProbTPCvsPtbeforePID(NULL), \r
 102   fHistProbTOFvsPtbeforePID(NULL), \r
 103   fHistProbTPCTOFvsPtbeforePID(NULL),\r
 104   fHistNSigmaTPCvsPtbeforePID(NULL), \r
 105   fHistNSigmaTOFvsPtbeforePID(NULL), \r
 106   fHistdEdxVsPTPCafterPID(NULL),\r
 107   fHistBetavsPTOFafterPID(NULL), \r
 108   fHistProbTPCvsPtafterPID(NULL), \r
 109   fHistProbTOFvsPtafterPID(NULL), \r
 110   fHistProbTPCTOFvsPtafterPID(NULL),\r
 111   fHistNSigmaTPCvsPtafterPID(NULL), \r
 112   fHistNSigmaTOFvsPtafterPID(NULL),  \r
 113   fCentralityArrayBinsForCorrections(kCENTRALITY),\r
 114   fPIDResponse(0x0),\r
 115   fPIDCombined(0x0),\r
 116   fParticleOfInterest(kPion),\r
 117   fPidDetectorConfig(kTPCTOF),\r
 118   fUsePID(kFALSE),\r
 119   fUsePIDnSigma(kTRUE),\r
 120   fUsePIDPropabilities(kFALSE), \r
 121   fPIDNSigma(3.),\r
 122   fMinAcceptedPIDProbability(0.8),\r
 123   fElectronRejection(kFALSE),\r
 124   fElectronOnlyRejection(kFALSE),\r
 125   fElectronRejectionNSigma(-1.),\r
 126   fElectronRejectionMinPt(0.),\r
 127   fElectronRejectionMaxPt(1000.),\r
 128   fESDtrackCuts(0),\r
 129   fCentralityEstimator("V0M"),\r
 130   fUseCentrality(kFALSE),\r
 131   fCentralityPercentileMin(0.), \r
 132   fCentralityPercentileMax(5.),\r
 133   fImpactParameterMin(0.),\r
 134   fImpactParameterMax(20.),\r
 135   fMultiplicityEstimator("V0M"),\r
 136   fUseMultiplicity(kFALSE),\r
 137   fNumberOfAcceptedTracksMin(0),\r
 138   fNumberOfAcceptedTracksMax(10000),\r
 139   fHistNumberOfAcceptedTracks(0),\r
 140   fHistMultiplicity(0),\r
 141   fUseOfflineTrigger(kFALSE),\r
 142   fCheckFirstEventInChunk(kFALSE),\r
 143   fCheckPileUp(kFALSE),\r
 144   fUseMCforKinematics(kFALSE),\r
 145   fVxMax(0.3),\r
 146   fVyMax(0.3),\r
 147   fVzMax(10.),\r
 148   fnAODtrackCutBit(128),\r
 149   fPtMin(0.3),\r
 150   fPtMax(1.5),\r
 151   fPtMinForCorrections(0.3),//=================================correction\r
 152   fPtMaxForCorrections(1.5),//=================================correction\r
 153   fPtBinForCorrections(36), //=================================correction\r
 154   fEtaMin(-0.8),\r
 155   fEtaMax(-0.8),\r
 156   fEtaMinForCorrections(-0.8),//=================================correction\r
 157   fEtaMaxForCorrections(0.8),//=================================correction\r
 158   fEtaBinForCorrections(16), //=================================correction\r
 159   fPhiMin(0.),\r
 160   fPhiMax(360.),\r
 161   fPhiMinForCorrections(0.),//=================================correction\r
 162   fPhiMaxForCorrections(360.),//=================================correction\r
 163   fPhiBinForCorrections(100), //=================================correction\r
 164   fDCAxyCut(-1),\r
 165   fDCAzCut(-1),\r
 166   fTPCchi2Cut(-1),\r
 167   fNClustersTPCCut(-1),\r
 168   fAcceptanceParameterization(0),\r
 169   fDifferentialV2(0),\r
 170   fUseFlowAfterBurner(kFALSE),\r
 171   fExcludeResonancesInMC(kFALSE),\r
 172   fExcludeElectronsInMC(kFALSE),\r
 173   fUseMCPdgCode(kFALSE),\r
 174   fPDGCodeToBeAnalyzed(-1),\r
 175   fEventClass("EventPlane") \r
 176 {\r
 177   // Constructor\r
 178   // Define input and output slots here\r
 179   // Input slot #0 works with a TChain\r
 180 \r
 181   //======================================================correction\r
 182   for (Int_t i=0; i<kCENTRALITY; i++){\r
 183     fHistCorrectionPlus[i] = NULL; \r
 184     fHistCorrectionMinus[i] = NULL; \r
 185     fCentralityArrayForCorrections[i] = -1.;\r
 186   }\r
 187   //=====================================================correction\r
 188 \r
 189   DefineInput(0, TChain::Class());\r
 190   // Output slot #0 writes into a TH1 container\r
 191   DefineOutput(1, TList::Class());\r
 192   DefineOutput(2, TList::Class());\r
 193   DefineOutput(3, TList::Class());\r
 194   DefineOutput(4, TList::Class());\r
 195   DefineOutput(5, TList::Class());\r
 196 }\r
 197 \r
 198 //________________________________________________________________________\r
 199 AliAnalysisTaskBFPsi::~AliAnalysisTaskBFPsi() {\r
 200 \r
 201   // delete fBalance; \r
 202   // delete fShuffledBalance; \r
 203   // delete fList;\r
 204   // delete fListBF; \r
 205   // delete fListBFS;\r
 206 \r
 207   // delete fHistEventStats; \r
 208   // delete fHistTrackStats; \r
 209   // delete fHistVx; \r
 210   // delete fHistVy; \r
 211   // delete fHistVz; \r
 212 \r
 213   // delete fHistClus;\r
 214   // delete fHistDCA;\r
 215   // delete fHistChi2;\r
 216   // delete fHistPt;\r
 217   // delete fHistEta;\r
 218   // delete fHistPhi;\r
 219   // delete fHistEtaPhiPos;                      \r
 220   // delete fHistEtaPhiNeg;\r
 221   // delete fHistV0M;\r
 222 }\r
 223 \r
 224 //________________________________________________________________________\r
 225 void AliAnalysisTaskBFPsi::UserCreateOutputObjects() {\r
 226   // Create histograms\r
 227   // Called once\r
 228 \r
 229   // global switch disabling the reference \r
 230   // (to avoid "Replacing existing TH1" if several wagons are created in train)\r
 231   Bool_t oldStatus = TH1::AddDirectoryStatus();\r
 232   TH1::AddDirectory(kFALSE);\r
 233 \r
 234   if(!fBalance) {\r
 235     fBalance = new AliBalancePsi();\r
 236     fBalance->SetAnalysisLevel("ESD");\r
 237     fBalance->SetEventClass(fEventClass);\r
 238     //fBalance->SetNumberOfBins(-1,16);\r
 239     //fBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);\r
 240   }\r
 241   if(fRunShuffling) {\r
 242     if(!fShuffledBalance) {\r
 243       fShuffledBalance = new AliBalancePsi();\r
 244       fShuffledBalance->SetAnalysisLevel("ESD");\r
 245       fShuffledBalance->SetEventClass(fEventClass);\r
 246       //fShuffledBalance->SetNumberOfBins(-1,16);\r
 247       //fShuffledBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);\r
 248     }\r
 249   }\r
 250   if(fRunMixing) {\r
 251     if(!fMixedBalance) {\r
 252       fMixedBalance = new AliBalancePsi();\r
 253       fMixedBalance->SetAnalysisLevel("ESD");\r
 254       fMixedBalance->SetEventClass(fEventClass);\r
 255     }\r
 256   }\r
 257 \r
 258   //QA list\r
 259   fList = new TList();\r
 260   fList->SetName("listQA");\r
 261   fList->SetOwner();\r
 262 \r
 263   //Balance Function list\r
 264   fListBF = new TList();\r
 265   fListBF->SetName("listBF");\r
 266   fListBF->SetOwner();\r
 267 \r
 268   if(fRunShuffling) {\r
 269     fListBFS = new TList();\r
 270     fListBFS->SetName("listBFShuffled");\r
 271     fListBFS->SetOwner();\r
 272   }\r
 273 \r
 274   if(fRunMixing) {\r
 275     fListBFM = new TList();\r
 276     fListBFM->SetName("listTriggeredBFMixed");\r
 277     fListBFM->SetOwner();\r
 278   }\r
 279 \r
 280   //PID QA list\r
 281   if(fUsePID || fElectronRejection) {\r
 282     fHistListPIDQA = new TList();\r
 283     fHistListPIDQA->SetName("listQAPID");\r
 284     fHistListPIDQA->SetOwner();\r
 285   }\r
 286 \r
 287   //Event stats.\r
 288   TString gCutName[7] = {"Total","Offline trigger",\r
 289                          "Vertex","Analyzed","sel. Centrality","Not1stEvInChunk","No Pile-Up"};\r
 290   fHistEventStats = new TH2F("fHistEventStats",\r
 291                              "Event statistics;;Centrality percentile;N_{events}",\r
 292                              7,0.5,7.5,220,-5,105);\r
 293   for(Int_t i = 1; i <= 7; i++)\r
 294     fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());\r
 295   fList->Add(fHistEventStats);\r
 296 \r
 297   TString gCentName[13] = {"V0M","V0A","V0C","FMD","TRK","TKL","CL0","CL1","ZNA","ZPA","V0MvsFMD","TKLvsV0M","ZEMvsZDC"};\r
 298   fHistCentStats = new TH2F("fHistCentStats",\r
 299                              "Centrality statistics;;Cent percentile",\r
 300                             13,-0.5,12.5,220,-5,105);\r
 301   for(Int_t i = 1; i <= 13; i++){\r
 302     fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());\r
 303     //fHistCentStatsUsed->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());  //++++++++++++++++++++++\r
 304   }\r
 305   fList->Add(fHistCentStats);\r
 306 \r
 307   fHistCentStatsUsed = new TH2F("fHistCentStatsUsed","Centrality statistics;;Cent percentile", 1,-0.5,0.5,220,-5,105); //++++++++++++++++++++++\r
 308   fHistCentStatsUsed->GetXaxis()->SetBinLabel(1,fCentralityEstimator.Data());  //++++++++++++++++++++++\r
 309   fList->Add(fHistCentStatsUsed); //++++++++++++++++++++++\r
 310 \r
 311   fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",1025,0,1025);\r
 312   fList->Add(fHistTriggerStats);\r
 313 \r
 314   fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TrackFilterBit;N_{events}",16,0,16);\r
 315   fList->Add(fHistTrackStats);\r
 316 \r
 317   fHistNumberOfAcceptedTracks = new TH2F("fHistNumberOfAcceptedTracks",";N_{acc.};Centrality percentile;Entries",4001,-0.5,4000.5,220,-5,105);\r
 318   fList->Add(fHistNumberOfAcceptedTracks);\r
 319 \r
 320   fHistMultiplicity = new TH1F("fHistMultiplicity",";N_{ch.};Entries",30001,-0.5,30000.5);\r
 321   fList->Add(fHistMultiplicity);\r
 322 \r
 323   // Vertex distributions\r
 324   fHistVx = new TH1F("fHistVx","Primary vertex distribution - x coordinate;V_{x} (cm);Entries",100,-0.5,0.5);\r
 325   fList->Add(fHistVx);\r
 326   fHistVy = new TH1F("fHistVy","Primary vertex distribution - y coordinate;V_{y} (cm);Entries",100,-0.5,0.5);\r
 327   fList->Add(fHistVy);\r
 328   fHistVz = new TH2F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Centrality percentile;Entries",100,-20.,20.,220,-5,105);\r
 329   fList->Add(fHistVz);\r
 330 \r
 331   //TPC vs VZERO multiplicity\r
 332   fHistTPCvsVZEROMultiplicity = new TH2F("fHistTPCvsVZEROMultiplicity","VZERO vs TPC multiplicity",3001,-0.5,30000.5,4001,-0.5,4000.5);\r
 333   if(fMultiplicityEstimator == "V0A") \r
 334     fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO multiplicity (a.u.)");\r
 335   else if(fMultiplicityEstimator == "V0C") \r
 336     fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO multiplicity (a.u.)");\r
 337   else \r
 338     fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO multiplicity (a.u.)");\r
 339   fList->Add(fHistTPCvsVZEROMultiplicity);\r
 340 \r
 341   fHistVZEROSignal = new TH2F("fHistVZEROSignal","VZERO signal vs VZERO channel;VZERO channel; Signal (a.u.)",64,0.5,64.5,3001,-0.5,30000.5);\r
 342   fList->Add(fHistVZEROSignal);\r
 343 \r
 344   //Event plane\r
 345   fHistEventPlane = new TH2F("fHistEventPlane",";#Psi_{2} [deg.];Centrality percentile;Counts",100,0,360.,220,-5,105);\r
 346   fList->Add(fHistEventPlane);\r
 347 \r
 348   // QA histograms\r
 349   fHistClus = new TH2F("fHistClus","# Cluster (TPC vs. ITS)",10,0,10,200,0,200);\r
 350   fList->Add(fHistClus);\r
 351   fHistChi2 = new TH2F("fHistChi2","Chi2/NDF distribution;#chi^{2}/ndf;Centrality percentile",200,0,10,220,-5,105);\r
 352   fList->Add(fHistChi2);\r
 353   fHistDCA  = new TH2F("fHistDCA","DCA (xy vs. z)",400,-5,5,400,-5,5); \r
 354   fList->Add(fHistDCA);\r
 355   fHistPt   = new TH2F("fHistPt","p_{T} distribution;p_{T} (GeV/c);Centrality percentile",200,0,10,220,-5,105);\r
 356   fList->Add(fHistPt);\r
 357   fHistEta  = new TH2F("fHistEta","#eta distribution;#eta;Centrality percentile",200,-2,2,220,-5,105);\r
 358   fList->Add(fHistEta);\r
 359   fHistRapidity  = new TH2F("fHistRapidity","y distribution;y;Centrality percentile",200,-2,2,220,-5,105);\r
 360   fList->Add(fHistRapidity);\r
 361   fHistPhi  = new TH2F("fHistPhi","#phi distribution;#phi (rad);Centrality percentile",200,0.0,2.*TMath::Pi(),220,-5,105);\r
 362   fList->Add(fHistPhi);\r
 363   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
 364   fList->Add(fHistEtaPhiPos);                    \r
 365   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
 366   fList->Add(fHistEtaPhiNeg);\r
 367   fHistPhiBefore  = new TH2F("fHistPhiBefore","#phi distribution;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);\r
 368   fList->Add(fHistPhiBefore);\r
 369   fHistPhiAfter  = new TH2F("fHistPhiAfter","#phi distribution;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);\r
 370   fList->Add(fHistPhiAfter);\r
 371   fHistPhiPos  = new TH2F("fHistPhiPos","#phi distribution for positive particles;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);\r
 372   fList->Add(fHistPhiPos);\r
 373   fHistPhiNeg  = new TH2F("fHistPhiNeg","#phi distribution for negative particles;#phi;Centrality percentile",200,0.,2.*TMath::Pi(),220,-5,105);\r
 374   fList->Add(fHistPhiNeg);\r
 375   fHistV0M  = new TH2F("fHistV0M","V0 Multiplicity C vs. A",500, 0, 20000, 500, 0, 20000);\r
 376   fList->Add(fHistV0M);\r
 377   TString gRefTrackName[6] = {"tracks","tracksPos","tracksNeg","tracksTPConly","clusITS0","clusITS1"};\r
 378   fHistRefTracks  = new TH2F("fHistRefTracks","Nr of Ref tracks/event vs. ref track estimator;;Nr of tracks",6, 0, 6, 400, 0, 20000);\r
 379   for(Int_t i = 1; i <= 6; i++)\r
 380     fHistRefTracks->GetXaxis()->SetBinLabel(i,gRefTrackName[i-1].Data());\r
 381   fList->Add(fHistRefTracks);\r
 382 \r
 383   // QA histograms for different cuts\r
 384   fList->Add(fBalance->GetQAHistHBTbefore());\r
 385   fList->Add(fBalance->GetQAHistHBTafter());\r
 386   fList->Add(fBalance->GetQAHistConversionbefore());\r
 387   fList->Add(fBalance->GetQAHistConversionafter());\r
 388   fList->Add(fBalance->GetQAHistPsiMinusPhi());\r
 389   fList->Add(fBalance->GetQAHistResonancesBefore());\r
 390   fList->Add(fBalance->GetQAHistResonancesRho());\r
 391   fList->Add(fBalance->GetQAHistResonancesK0());\r
 392   fList->Add(fBalance->GetQAHistResonancesLambda());\r
 393   fList->Add(fBalance->GetQAHistQbefore());\r
 394   fList->Add(fBalance->GetQAHistQafter());\r
 395 \r
 396   // Balance function histograms\r
 397   // Initialize histograms if not done yet\r
 398   if(!fBalance->GetHistNp()){\r
 399     AliWarning("Histograms not yet initialized! --> Will be done now");\r
 400     AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
 401     fBalance->InitHistograms();\r
 402   }\r
 403 \r
 404   if(fRunShuffling) {\r
 405     if(!fShuffledBalance->GetHistNp()) {\r
 406       AliWarning("Histograms (shuffling) not yet initialized! --> Will be done now");\r
 407       AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
 408       fShuffledBalance->InitHistograms();\r
 409     }\r
 410   }\r
 411 \r
 412   //for(Int_t a = 0; a < ANALYSIS_TYPES; a++){\r
 413   fListBF->Add(fBalance->GetHistNp());\r
 414   fListBF->Add(fBalance->GetHistNn());\r
 415   fListBF->Add(fBalance->GetHistNpn());\r
 416   fListBF->Add(fBalance->GetHistNnn());\r
 417   fListBF->Add(fBalance->GetHistNpp());\r
 418   fListBF->Add(fBalance->GetHistNnp());\r
 419 \r
 420   if(fRunShuffling) {\r
 421     fListBFS->Add(fShuffledBalance->GetHistNp());\r
 422     fListBFS->Add(fShuffledBalance->GetHistNn());\r
 423     fListBFS->Add(fShuffledBalance->GetHistNpn());\r
 424     fListBFS->Add(fShuffledBalance->GetHistNnn());\r
 425     fListBFS->Add(fShuffledBalance->GetHistNpp());\r
 426     fListBFS->Add(fShuffledBalance->GetHistNnp());\r
 427   }  \r
 428 \r
 429   if(fRunMixing) {\r
 430     fListBFM->Add(fMixedBalance->GetHistNp());\r
 431     fListBFM->Add(fMixedBalance->GetHistNn());\r
 432     fListBFM->Add(fMixedBalance->GetHistNpn());\r
 433     fListBFM->Add(fMixedBalance->GetHistNnn());\r
 434     fListBFM->Add(fMixedBalance->GetHistNpp());\r
 435     fListBFM->Add(fMixedBalance->GetHistNnp());\r
 436   }\r
 437   //}\r
 438 \r
 439 \r
 440   // Event Mixing\r
 441   if(fRunMixing){\r
 442     Int_t trackDepth = fMixingTracks; \r
 443     Int_t poolsize   = 1000;  // Maximum number of events, ignored in the present implemented of AliEventPoolManager\r
 444     \r
 445     // centrality bins\r
 446     //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
 447     Double_t centralityBins[] = {0.,1.,2.,3.,4.,5.,7.,10.,15.,20.,25.,30.,35.,40.,45.,50.,55.,60.,65.,70.,75.,80.};\r
 448     Double_t* centbins        = centralityBins;\r
 449     Int_t nCentralityBins     = sizeof(centralityBins) / sizeof(Double_t) - 1;\r
 450 \r
 451     // multiplicity bins\r
 452     Double_t multiplicityBins[] = {0,10,20,30,40,50,60,70,80,100,100000}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
 453     Double_t* multbins        = multiplicityBins;\r
 454     Int_t nMultiplicityBins     = sizeof(multiplicityBins) / sizeof(Double_t) - 1;\r
 455     \r
 456     // Zvtx bins\r
 457     Double_t vertexBins[] = {-10., -7., -5., -3., -1., 1., 3., 5., 7., 10.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
 458     Double_t* vtxbins     = vertexBins;\r
 459     Int_t nVertexBins     = sizeof(vertexBins) / sizeof(Double_t) - 1;\r
 460     \r
 461     // Event plane angle (Psi) bins\r
 462     Double_t psiBins[] = {0.,45.,135.,215.,305.,360.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
 463     Double_t* psibins     = psiBins;\r
 464     Int_t nPsiBins     = sizeof(psiBins) / sizeof(Double_t) - 1;\r
 465     \r
 466     // run the event mixing also in bins of event plane (statistics!)\r
 467     if(fRunMixingEventPlane){\r
 468       if(fEventClass=="Multiplicity"){\r
 469         fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins, nPsiBins, psibins);\r
 470       }\r
 471       else{\r
 472         fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins, nPsiBins, psibins);\r
 473       }\r
 474     }\r
 475     else{\r
 476       if(fEventClass=="Multiplicity"){\r
 477         fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins);\r
 478       }\r
 479       else{\r
 480         fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);\r
 481       }\r
 482     }\r
 483   }\r
 484 \r
 485   if(fESDtrackCuts) fList->Add(fESDtrackCuts);\r
 486 \r
 487   //====================PID========================//\r
 488   if(fUsePID) {\r
 489     fPIDCombined = new AliPIDCombined();\r
 490     fPIDCombined->SetDefaultTPCPriors();\r
 491 \r
 492     fHistdEdxVsPTPCbeforePID = new TH2D ("dEdxVsPTPCbefore","dEdxVsPTPCbefore", 1000, -10.0, 10.0, 1000, 0, 1000); \r
 493     fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePID); //addition \r
 494     \r
 495     fHistBetavsPTOFbeforePID = new TH2D ("BetavsPTOFbefore","BetavsPTOFbefore", 1000, -10.0, 10., 1000, 0, 1.2); \r
 496     fHistListPIDQA->Add(fHistBetavsPTOFbeforePID); //addition\r
 497     \r
 498     fHistProbTPCvsPtbeforePID = new TH2D ("ProbTPCvsPtbefore","ProbTPCvsPtbefore", 1000, -10.0,10.0, 1000, 0, 2.0); \r
 499     fHistListPIDQA->Add(fHistProbTPCvsPtbeforePID); //addition \r
 500     \r
 501     fHistProbTOFvsPtbeforePID = new TH2D ("ProbTOFvsPtbefore","ProbTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0); \r
 502     fHistListPIDQA->Add(fHistProbTOFvsPtbeforePID); //addition \r
 503 \r
 504     fHistProbTPCTOFvsPtbeforePID =new TH2D ("ProbTPCTOFvsPtbefore","ProbTPCTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0); \r
 505     fHistListPIDQA->Add(fHistProbTPCTOFvsPtbeforePID); //addition \r
 506     \r
 507     fHistNSigmaTPCvsPtbeforePID = new TH2D ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore", 1000, -10, 10, 1000, 0, 500); \r
 508     fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePID); //addition \r
 509     \r
 510     fHistNSigmaTOFvsPtbeforePID = new TH2D ("NSigmaTOFvsPtbefore","NSigmaTOFvsPtbefore", 1000, -10, 10, 1000, 0, 500); \r
 511     fHistListPIDQA->Add(fHistNSigmaTOFvsPtbeforePID); //addition \r
 512     \r
 513     fHistdEdxVsPTPCafterPID = new TH2D ("dEdxVsPTPCafter","dEdxVsPTPCafter", 1000, -10, 10, 1000, 0, 1000); \r
 514     fHistListPIDQA->Add(fHistdEdxVsPTPCafterPID); //addition \r
 515     \r
 516     fHistBetavsPTOFafterPID = new TH2D ("BetavsPTOFafter","BetavsPTOFafter", 1000, -10, 10, 1000, 0, 1.2); \r
 517     fHistListPIDQA->Add(fHistBetavsPTOFafterPID); //addition \r
 518     \r
 519     fHistProbTPCvsPtafterPID = new TH2D ("ProbTPCvsPtafter","ProbTPCvsPtafter", 1000, -10, 10, 1000, 0, 2); \r
 520     fHistListPIDQA->Add(fHistProbTPCvsPtafterPID); //addition \r
 521   \r
 522     fHistProbTOFvsPtafterPID = new TH2D ("ProbTOFvsPtafter","ProbTOFvsPtafter", 1000,  -10, 10, 1000, 0, 2); \r
 523     fHistListPIDQA->Add(fHistProbTOFvsPtafterPID); //addition  \r
 524     \r
 525     fHistProbTPCTOFvsPtafterPID =new TH2D ("ProbTPCTOFvsPtafter","ProbTPCTOFvsPtafter", 1000, -50, 50, 1000, 0, 2.0); \r
 526     fHistListPIDQA->Add(fHistProbTPCTOFvsPtafterPID); //addition \r
 527 \r
 528     fHistNSigmaTPCvsPtafterPID = new TH2D ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter", 1000, -10, 10, 1000, 0, 500); \r
 529     fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPID); //addition  \r
 530     \r
 531     fHistNSigmaTOFvsPtafterPID = new TH2D ("NSigmaTOFvsPtafter","NSigmaTOFvsPtafter", 1000, -10, 10, 1000, 0, 500); \r
 532     fHistListPIDQA->Add(fHistNSigmaTOFvsPtafterPID); //addition \r
 533   }\r
 534 \r
 535   // for electron rejection only TPC nsigma histograms\r
 536   if(!fUsePID && fElectronRejection) {\r
 537  \r
 538     fHistdEdxVsPTPCbeforePID = new TH2D ("dEdxVsPTPCbefore","dEdxVsPTPCbefore", 1000, -10.0, 10.0, 1000, 0, 1000); \r
 539     fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePID); //addition \r
 540     \r
 541     fHistNSigmaTPCvsPtbeforePID = new TH2D ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore", 1000, -10, 10, 1000, 0, 500); \r
 542     fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePID); //addition \r
 543     \r
 544     fHistdEdxVsPTPCafterPID = new TH2D ("dEdxVsPTPCafter","dEdxVsPTPCafter", 1000, -10, 10, 1000, 0, 1000); \r
 545     fHistListPIDQA->Add(fHistdEdxVsPTPCafterPID); //addition \r
 546 \r
 547     fHistNSigmaTPCvsPtafterPID = new TH2D ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter", 1000, -10, 10, 1000, 0, 500); \r
 548     fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPID); //addition  \r
 549   }\r
 550   //====================PID========================//\r
 551 \r
 552   // Post output data.\r
 553   PostData(1, fList);\r
 554   PostData(2, fListBF);\r
 555   if(fRunShuffling) PostData(3, fListBFS);\r
 556   if(fRunMixing) PostData(4, fListBFM);\r
 557   if(fUsePID || fElectronRejection) PostData(5, fHistListPIDQA);       //PID\r
 558 \r
 559   TH1::AddDirectory(oldStatus);\r
 560 }\r
 561 \r
 562 \r
 563 //________________________________________________________________________\r
 564 void AliAnalysisTaskBFPsi::SetInputCorrection(TString filename, \r
 565                                               Int_t nCentralityBins, \r
 566                                               Double_t *centralityArrayForCorrections) {\r
 567   //Open files that will be used for correction\r
 568   fCentralityArrayBinsForCorrections = nCentralityBins;\r
 569   for (Int_t i=0; i<nCentralityBins; i++)\r
 570     fCentralityArrayForCorrections[i] = centralityArrayForCorrections[i];\r
 571 \r
 572   // No file specified -> run without corrections\r
 573   if(!filename.Contains(".root")) {\r
 574     AliInfo(Form("No correction file specified (= %s) --> run without corrections",filename.Data()));\r
 575     return;\r
 576   }\r
 577 \r
 578   //Open the input file\r
 579   TFile *f = TFile::Open(filename);\r
 580   if(!f->IsOpen()) {\r
 581     AliInfo(Form("File %s not found --> run without corrections",filename.Data()));\r
 582     return;\r
 583   }\r
 584     \r
 585   //TString listEffName = "";\r
 586   for (Int_t iCent = 0; iCent < fCentralityArrayBinsForCorrections-1; iCent++) {    \r
 587     //Printf("iCent %d:",iCent);    \r
 588     TString histoName = "fHistCorrectionPlus";\r
 589     histoName += Form("%d-%d",(Int_t)(fCentralityArrayForCorrections[iCent]),(Int_t)(fCentralityArrayForCorrections[iCent+1]));\r
 590     fHistCorrectionPlus[iCent]= dynamic_cast<TH3F *>(f->Get(histoName.Data()));\r
 591     if(!fHistCorrectionPlus[iCent]) {\r
 592       AliError(Form("fHist %s not found!!!",histoName.Data()));\r
 593       return;\r
 594     }\r
 595     \r
 596     histoName = "fHistCorrectionMinus";\r
 597     histoName += Form("%d-%d",(Int_t)(fCentralityArrayForCorrections[iCent]),(Int_t)(fCentralityArrayForCorrections[iCent+1]));\r
 598     fHistCorrectionMinus[iCent] = dynamic_cast<TH3F *>(f->Get(histoName.Data())); \r
 599     if(!fHistCorrectionMinus[iCent]) {\r
 600       AliError(Form("fHist %s not found!!!",histoName.Data()));\r
 601       return; \r
 602     }\r
 603   }//loop over centralities: ONLY the PbPb case is covered\r
 604 \r
 605   if(fHistCorrectionPlus[0]){\r
 606     fEtaMinForCorrections = fHistCorrectionPlus[0]->GetXaxis()->GetXmin();\r
 607     fEtaMaxForCorrections = fHistCorrectionPlus[0]->GetXaxis()->GetXmax();\r
 608     fEtaBinForCorrections = fHistCorrectionPlus[0]->GetNbinsX();\r
 609     \r
 610     fPtMinForCorrections = fHistCorrectionPlus[0]->GetYaxis()->GetXmin();\r
 611     fPtMaxForCorrections = fHistCorrectionPlus[0]->GetYaxis()->GetXmax();\r
 612     fPtBinForCorrections = fHistCorrectionPlus[0]->GetNbinsY();\r
 613     \r
 614     fPhiMinForCorrections = fHistCorrectionPlus[0]->GetZaxis()->GetXmin();\r
 615     fPhiMaxForCorrections = fHistCorrectionPlus[0]->GetZaxis()->GetXmax();\r
 616     fPhiBinForCorrections = fHistCorrectionPlus[0]->GetNbinsZ();\r
 617   }\r
 618 }\r
 619 \r
 620 //________________________________________________________________________\r
 621 void AliAnalysisTaskBFPsi::UserExec(Option_t *) {\r
 622   // Main loop\r
 623   // Called for each event\r
 624 \r
 625   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
 626   Int_t gNumberOfAcceptedTracks = 0;\r
 627   Double_t gCentrality          = -1.;\r
 628   Double_t gReactionPlane       = -1.; \r
 629   Float_t bSign = 0.;\r
 630   \r
 631   // get the event (for generator level: MCEvent())\r
 632   AliVEvent* eventMain = NULL;\r
 633   if(gAnalysisLevel == "MC") {\r
 634     eventMain = dynamic_cast<AliVEvent*>(MCEvent()); \r
 635   }\r
 636   else{\r
 637     eventMain = dynamic_cast<AliVEvent*>(InputEvent()); \r
 638     \r
 639     // for HBT like cuts need magnetic field sign\r
 640     bSign = (eventMain->GetMagneticField() > 0) ? 1 : -1;\r
 641   }\r
 642   if(!eventMain) {\r
 643     AliError("eventMain not available");\r
 644     return;\r
 645   }\r
 646   \r
 647   // PID Response task active?\r
 648   if(fUsePID || fElectronRejection) {\r
 649     fPIDResponse = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();\r
 650     if (!fPIDResponse) AliFatal("This Task needs the PID response attached to the inputHandler");\r
 651   }\r
 652   \r
 653   // check event cuts and fill event histograms\r
 654   if((gCentrality = IsEventAccepted(eventMain)) < 0){\r
 655     return;\r
 656   }\r
 657   \r
 658   //Compute Multiplicity or Centrality variable\r
 659   Double_t lMultiplicityVar = GetRefMultiOrCentrality( eventMain );\r
 660 \r
 661   // get the reaction plane\r
 662   if(fEventClass != "Multiplicity") {\r
 663     gReactionPlane = GetEventPlane(eventMain);\r
 664     fHistEventPlane->Fill(gReactionPlane,lMultiplicityVar);\r
 665     if(gReactionPlane < 0){\r
 666       return;\r
 667     }\r
 668   }\r
 669   \r
 670   // get the accepted tracks in main event\r
 671   TObjArray *tracksMain = GetAcceptedTracks(eventMain,lMultiplicityVar,gReactionPlane);\r
 672   gNumberOfAcceptedTracks = tracksMain->GetEntriesFast();\r
 673 \r
 674   //multiplicity cut (used in pp)\r
 675   fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedTracks,lMultiplicityVar);\r
 676 \r
 677   // store charges of all accepted tracks,shuffle and reassign(two extra loops)\r
 678   TObjArray* tracksShuffled = NULL;\r
 679   if(fRunShuffling){\r
 680     tracksShuffled = GetShuffledTracks(tracksMain,lMultiplicityVar);\r
 681   }\r
 682   \r
 683   // Event mixing \r
 684   if (fRunMixing)\r
 685     {\r
 686       // 1. First get an event pool corresponding in mult (cent) and\r
 687       //    zvertex to the current event. Once initialized, the pool\r
 688       //    should contain nMix (reduced) events. This routine does not\r
 689       //    pre-scan the chain. The first several events of every chain\r
 690       //    will be skipped until the needed pools are filled to the\r
 691       //    specified depth. If the pool categories are not too rare, this\r
 692       //    should not be a problem. If they are rare, you could lose`\r
 693       //    statistics.\r
 694       \r
 695       // 2. Collect the whole pool's content of tracks into one TObjArray\r
 696       //    (bgTracks), which is effectively a single background super-event.\r
 697       \r
 698       // 3. The reduced and bgTracks arrays must both be passed into\r
 699       //    FillCorrelations(). Also nMix should be passed in, so a weight\r
 700       //    of 1./nMix can be applied.\r
 701       \r
 702       AliEventPool* pool = fPoolMgr->GetEventPool(lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane);\r
 703       \r
 704       if (!pool){\r
 705         AliFatal(Form("No pool found for centrality = %f, zVtx = %f, psi = %f", lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane));\r
 706       }\r
 707       else{\r
 708         \r
 709         //pool->SetDebug(1);\r
 710 \r
 711         if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 || pool->GetCurrentNEvents() >= 5){ \r
 712           \r
 713           \r
 714           Int_t nMix = pool->GetCurrentNEvents();\r
 715           //cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << endl;\r
 716           \r
 717           //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(2);\r
 718           //((TH2F*) fListOfHistos->FindObject("mixedDist"))->Fill(centrality, pool->NTracksInPool());\r
 719           //if (pool->IsReady())\r
 720           //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(3);\r
 721           \r
 722           // Fill mixed-event histos here  \r
 723           for (Int_t jMix=0; jMix<nMix; jMix++) \r
 724             {\r
 725               TObjArray* tracksMixed = pool->GetEvent(jMix);\r
 726               fMixedBalance->CalculateBalance(gReactionPlane,tracksMain,tracksMixed,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());\r
 727             }\r
 728         }\r
 729         \r
 730         // Update the Event pool\r
 731         pool->UpdatePool(tracksMain);\r
 732         //pool->PrintInfo();\r
 733         \r
 734       }//pool NULL check  \r
 735     }//run mixing\r
 736   \r
 737   // calculate balance function\r
 738   fBalance->CalculateBalance(gReactionPlane,tracksMain,NULL,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());\r
 739   \r
 740   // calculate shuffled balance function\r
 741   if(fRunShuffling && tracksShuffled != NULL) {\r
 742     fShuffledBalance->CalculateBalance(gReactionPlane,tracksShuffled,NULL,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());\r
 743   }\r
 744 }      \r
 745 \r
 746 //________________________________________________________________________\r
 747 Double_t AliAnalysisTaskBFPsi::IsEventAccepted(AliVEvent *event){\r
 748   // Checks the Event cuts\r
 749   // Fills Event statistics histograms\r
 750   \r
 751   Bool_t isSelectedMain = kTRUE;\r
 752   Float_t gCentrality = -1.;\r
 753   Float_t gRefMultiplicity = -1.;\r
 754   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
 755 \r
 756   AliMCEvent *mcevent = dynamic_cast<AliMCEvent*>(event);\r
 757 \r
 758   fHistEventStats->Fill(1,gCentrality); //all events\r
 759 \r
 760   // check first event in chunk (is not needed for new reconstructions)\r
 761   if(fCheckFirstEventInChunk){\r
 762     AliAnalysisUtils ut;\r
 763     if(ut.IsFirstEventInChunk(event)) \r
 764       return -1.;\r
 765     fHistEventStats->Fill(6,gCentrality); \r
 766   }\r
 767 \r
 768   // check for pile-up event\r
 769   if(fCheckPileUp){\r
 770     AliAnalysisUtils ut;\r
 771     ut.SetUseMVPlpSelection(kTRUE);\r
 772     ut.SetUseOutOfBunchPileUp(kTRUE);\r
 773     if(ut.IsPileUpEvent(event))\r
 774       return -1.;\r
 775     fHistEventStats->Fill(7,gCentrality); \r
 776   }\r
 777 \r
 778   // Event trigger bits\r
 779   fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());\r
 780   if(fUseOfflineTrigger)\r
 781     isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
 782   \r
 783   if(isSelectedMain) {\r
 784     fHistEventStats->Fill(2,gCentrality); //triggered events\r
 785 \r
 786     //Centrality stuff \r
 787     if(fUseCentrality) {\r
 788       if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD" || gAnalysisLevel == "MCAODrec") { //centrality in AOD header   //+++++++++++\r
 789         AliAODHeader *header = (AliAODHeader*) event->GetHeader();\r
 790         if(header){\r
 791           gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
 792 \r
 793           // QA for centrality estimators\r
 794           fHistCentStats->Fill(0.,header->GetCentralityP()->GetCentralityPercentile("V0M"));\r
 795           fHistCentStats->Fill(1.,header->GetCentralityP()->GetCentralityPercentile("V0A"));\r
 796           fHistCentStats->Fill(2.,header->GetCentralityP()->GetCentralityPercentile("V0C"));\r
 797           fHistCentStats->Fill(3.,header->GetCentralityP()->GetCentralityPercentile("FMD"));\r
 798           fHistCentStats->Fill(4.,header->GetCentralityP()->GetCentralityPercentile("TRK"));\r
 799           fHistCentStats->Fill(5.,header->GetCentralityP()->GetCentralityPercentile("TKL")); \r
 800           fHistCentStats->Fill(6.,header->GetCentralityP()->GetCentralityPercentile("CL0"));\r
 801           fHistCentStats->Fill(7.,header->GetCentralityP()->GetCentralityPercentile("CL1"));\r
 802           fHistCentStats->Fill(8.,header->GetCentralityP()->GetCentralityPercentile("ZNA"));\r
 803           fHistCentStats->Fill(9.,header->GetCentralityP()->GetCentralityPercentile("ZPA"));\r
 804           fHistCentStats->Fill(10.,header->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));\r
 805           fHistCentStats->Fill(11.,header->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));\r
 806           fHistCentStats->Fill(12.,header->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));\r
 807 \r
 808           // Centrality estimator USED   ++++++++++++++++++++++++++++++\r
 809           fHistCentStatsUsed->Fill(0.,header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data()));\r
 810           \r
 811           // centrality QA (V0M)\r
 812           fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());\r
 813           \r
 814           // centrality QA (reference tracks)\r
 815           fHistRefTracks->Fill(0.,header->GetRefMultiplicity());\r
 816           fHistRefTracks->Fill(1.,header->GetRefMultiplicityPos());\r
 817           fHistRefTracks->Fill(2.,header->GetRefMultiplicityNeg());\r
 818           fHistRefTracks->Fill(3.,header->GetTPConlyRefMultiplicity());\r
 819           fHistRefTracks->Fill(4.,header->GetNumberOfITSClusters(0));\r
 820           fHistRefTracks->Fill(5.,header->GetNumberOfITSClusters(1));\r
 821           fHistRefTracks->Fill(6.,header->GetNumberOfITSClusters(2));\r
 822           fHistRefTracks->Fill(7.,header->GetNumberOfITSClusters(3));\r
 823           fHistRefTracks->Fill(8.,header->GetNumberOfITSClusters(4));\r
 824         }//AOD header\r
 825       }//AOD\r
 826 \r
 827       else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){ // centrality class for ESDs or MC-ESDs\r
 828         AliCentrality *centrality = event->GetCentrality();\r
 829         gCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
 830 \r
 831         // QA for centrality estimators\r
 832         fHistCentStats->Fill(0.,centrality->GetCentralityPercentile("V0M"));\r
 833         fHistCentStats->Fill(1.,centrality->GetCentralityPercentile("V0A"));\r
 834         fHistCentStats->Fill(2.,centrality->GetCentralityPercentile("V0C"));\r
 835         fHistCentStats->Fill(3.,centrality->GetCentralityPercentile("FMD"));\r
 836         fHistCentStats->Fill(4.,centrality->GetCentralityPercentile("TRK"));\r
 837         fHistCentStats->Fill(5.,centrality->GetCentralityPercentile("TKL"));\r
 838         fHistCentStats->Fill(6.,centrality->GetCentralityPercentile("CL0"));\r
 839         fHistCentStats->Fill(7.,centrality->GetCentralityPercentile("CL1"));\r
 840         fHistCentStats->Fill(8.,centrality->GetCentralityPercentile("ZNA"));\r
 841         fHistCentStats->Fill(9.,centrality->GetCentralityPercentile("ZPA"));\r
 842         fHistCentStats->Fill(10.,centrality->GetCentralityPercentile("V0MvsFMD"));\r
 843         fHistCentStats->Fill(11.,centrality->GetCentralityPercentile("TKLvsV0M"));\r
 844         fHistCentStats->Fill(12.,centrality->GetCentralityPercentile("ZEMvsZDC"));\r
 845 \r
 846         // Centrality estimator USED   ++++++++++++++++++++++++++++++\r
 847         fHistCentStatsUsed->Fill(0.,centrality->GetCentralityPercentile(fCentralityEstimator.Data()));\r
 848 \r
 849         // centrality QA (V0M)\r
 850         fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());\r
 851       }//ESD\r
 852       else if(gAnalysisLevel == "MC"){\r
 853         Double_t gImpactParameter = 0.;\r
 854         if(mcevent) {\r
 855           AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(dynamic_cast<AliMCEvent*>(mcevent)->GenEventHeader());\r
 856           if(headerH){\r
 857             gImpactParameter = headerH->ImpactParameter();\r
 858             gCentrality      = gImpactParameter;\r
 859           }//MC header\r
 860         }//MC event cast\r
 861       }//MC\r
 862       else{\r
 863         gCentrality = -1.;\r
 864       }\r
 865     }//centrality\r
 866 \r
 867     //Multiplicity stuff \r
 868     if(fUseMultiplicity)\r
 869       gRefMultiplicity = GetRefMultiOrCentrality(event);\r
 870 \r
 871     // Event Vertex MC\r
 872     if(gAnalysisLevel == "MC") {\r
 873       if(!event) {\r
 874         AliError("mcEvent not available");\r
 875         return 0x0;\r
 876       }\r
 877       \r
 878       if(mcevent){\r
 879         AliGenEventHeader *header = dynamic_cast<AliGenEventHeader*>(mcevent->GenEventHeader());\r
 880         if(header){  \r
 881           TArrayF gVertexArray;\r
 882           header->PrimaryVertex(gVertexArray);\r
 883           //Printf("Vertex: %lf (x) - %lf (y) - %lf (z)",\r
 884           //gVertexArray.At(0),\r
 885           //gVertexArray.At(1),\r
 886           //gVertexArray.At(2));\r
 887           if(fUseMultiplicity) \r
 888             fHistEventStats->Fill(3,gRefMultiplicity); //events with a proper vertex\r
 889           else \r
 890             fHistEventStats->Fill(3,gCentrality); //events with a proper vertex\r
 891           if(TMath::Abs(gVertexArray.At(0)) < fVxMax) {\r
 892             if(TMath::Abs(gVertexArray.At(1)) < fVyMax) {\r
 893               if(TMath::Abs(gVertexArray.At(2)) < fVzMax) {\r
 894                 if(fUseMultiplicity) \r
 895                   fHistEventStats->Fill(4,gRefMultiplicity);//analyzed events\r
 896                 else \r
 897                   fHistEventStats->Fill(4,gCentrality); //analyzed events\r
 898                 fHistVx->Fill(gVertexArray.At(0));\r
 899                 fHistVy->Fill(gVertexArray.At(1));\r
 900                 fHistVz->Fill(gVertexArray.At(2),gCentrality);\r
 901                 \r
 902                 // take only events inside centrality class\r
 903                 if(fUseCentrality) {\r
 904                   if((fImpactParameterMin < gCentrality) && (fImpactParameterMax > gCentrality)){\r
 905                     fHistEventStats->Fill(5,gCentrality); //events with correct centrality\r
 906                     return gCentrality;     \r
 907                   }//centrality class\r
 908                 }\r
 909                 // take events only within the same multiplicity class\r
 910                 else if(fUseMultiplicity) {\r
 911                   if((gRefMultiplicity > fNumberOfAcceptedTracksMin)||(gRefMultiplicity < fNumberOfAcceptedTracksMax)) {\r
 912                     fHistEventStats->Fill(5,gRefMultiplicity); //events with correct multiplicity\r
 913                     return gRefMultiplicity;\r
 914                   }\r
 915                 }//multiplicity range\r
 916               }//Vz cut\r
 917             }//Vy cut\r
 918           }//Vx cut\r
 919         }//header    \r
 920       }//MC event object\r
 921     }//MC\r
 922     \r
 923     // Event Vertex AOD, ESD, ESDMC\r
 924     else{\r
 925       const AliVVertex *vertex = event->GetPrimaryVertex();\r
 926       \r
 927       if(vertex) {\r
 928         Double32_t fCov[6];\r
 929         vertex->GetCovarianceMatrix(fCov);\r
 930         if(vertex->GetNContributors() > 0) {\r
 931           if(fCov[5] != 0) {\r
 932             if(fUseMultiplicity) \r
 933               fHistEventStats->Fill(3,gRefMultiplicity); //proper vertex\r
 934             else if(fUseCentrality)\r
 935               fHistEventStats->Fill(3,gCentrality); //proper vertex\r
 936             if(TMath::Abs(vertex->GetX()) < fVxMax) {\r
 937               if(TMath::Abs(vertex->GetY()) < fVyMax) {\r
 938                 if(TMath::Abs(vertex->GetZ()) < fVzMax) {\r
 939                   if(fUseMultiplicity) {\r
 940                     //if(fDebugLevel)\r
 941                     //Printf("Filling 4th bin of fHistEventStats for multiplicity %lf",gRefMultiplicity);\r
 942                     fHistEventStats->Fill(4,gRefMultiplicity);//analyzed events\r
 943                   }\r
 944                   else if(fUseCentrality) {\r
 945                     //cout<<"Filling 4 for centrality..."<<endl;\r
 946                     fHistEventStats->Fill(4,gCentrality); //analyzed events\r
 947                   }\r
 948                   fHistVx->Fill(vertex->GetX());\r
 949                   fHistVy->Fill(vertex->GetY());\r
 950                   fHistVz->Fill(vertex->GetZ(),gCentrality);\r
 951                   \r
 952                   // take only events inside centrality class\r
 953                   if(fUseCentrality) {\r
 954                     //cout<<"Centrality..."<<endl;\r
 955                     if((gCentrality > fCentralityPercentileMin) && (gCentrality < fCentralityPercentileMax)){\r
 956                       fHistEventStats->Fill(5,gCentrality); //events with correct centrality\r
 957                       return gCentrality;               \r
 958                     }//centrality class\r
 959                   }\r
 960                   // take events only within the same multiplicity class\r
 961                   else if(fUseMultiplicity) {\r
 962                     //if(fDebugLevel) \r
 963                     //Printf("N(min): %.0f, N(max): %.0f - N(ref): %.0f",fNumberOfAcceptedTracksMin,\r
 964                     //fNumberOfAcceptedTracksMax,gRefMultiplicity);\r
 965 \r
 966                     if((gRefMultiplicity > fNumberOfAcceptedTracksMin)||(gRefMultiplicity < fNumberOfAcceptedTracksMax)) {\r
 967                       fHistEventStats->Fill(5,gRefMultiplicity); //events with correct multiplicity\r
 968                       return gRefMultiplicity;\r
 969                     }\r
 970                   }//multiplicity range\r
 971                 }//Vz cut\r
 972               }//Vy cut\r
 973             }//Vx cut\r
 974           }//proper vertex resolution\r
 975         }//proper number of contributors\r
 976       }//vertex object valid\r
 977     }//triggered event \r
 978   }//AOD,ESD,ESDMC\r
 979   \r
 980   // in all other cases return -1 (event not accepted)\r
 981   return -1;\r
 982 }\r
 983 \r
 984 \r
 985 //________________________________________________________________________\r
 986 Double_t AliAnalysisTaskBFPsi::GetRefMultiOrCentrality(AliVEvent *event){\r
 987     // Checks the Event cuts\r
 988     // Fills Event statistics histograms\r
 989   \r
 990   Float_t gCentrality = -1.;\r
 991   Double_t gMultiplicity = 0.;\r
 992   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
 993 \r
 994   if(fEventClass == "Centrality"){\r
 995     if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD" || gAnalysisLevel == "MCAODrec" ) { //centrality in AOD header  //++++++++++++++\r
 996       AliAODHeader *header = (AliAODHeader*) event->GetHeader();\r
 997       if(header){\r
 998         gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
 999       }//AOD header\r
1000     }//AOD\r
1001     \r
1002     else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){ // centrality class for ESDs or MC-ESDs\r
1003       AliCentrality *centrality = event->GetCentrality();\r
1004       gCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
1005     }//ESD\r
1006     else if(gAnalysisLevel == "MC"){\r
1007       Double_t gImpactParameter = 0.;\r
1008       AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);\r
1009       if(gMCEvent){\r
1010         AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader());      \r
1011         if(headerH){\r
1012           gImpactParameter = headerH->ImpactParameter();\r
1013           gCentrality      = gImpactParameter;\r
1014         }//MC header\r
1015       }//MC event cast\r
1016     }//MC\r
1017     else{\r
1018       gCentrality = -1.;\r
1019     }\r
1020   }//End if "Centrality"\r
1021   if(fEventClass=="Multiplicity"&&gAnalysisLevel == "ESD"){\r
1022     AliESDEvent* gESDEvent = dynamic_cast<AliESDEvent*>(event);\r
1023     if(gESDEvent){\r
1024       gMultiplicity = fESDtrackCuts->GetReferenceMultiplicity(gESDEvent, AliESDtrackCuts::kTrackletsITSTPC,0.5);\r
1025       fHistMultiplicity->Fill(gMultiplicity);\r
1026     }//AliESDevent cast\r
1027   }\r
1028   else if(fEventClass=="Multiplicity"&&gAnalysisLevel == "AOD"){\r
1029     AliAODHeader *header = (AliAODHeader*) event->GetHeader();\r
1030     if ((fMultiplicityEstimator == "V0M")||\r
1031         (fMultiplicityEstimator == "V0A")||\r
1032         (fMultiplicityEstimator == "V0C") ||\r
1033         (fMultiplicityEstimator == "TPC")) {\r
1034       gMultiplicity = GetReferenceMultiplicityFromAOD(event);\r
1035       if(fDebugLevel) Printf("Reference multiplicity (calculated): %.0f",gMultiplicity);\r
1036     }\r
1037     else {\r
1038       if(header)\r
1039         gMultiplicity = header->GetRefMultiplicity();\r
1040       if(fDebugLevel) Printf("Reference multiplicity (AOD header): %.0f",gMultiplicity);\r
1041     }\r
1042     fHistMultiplicity->Fill(gMultiplicity);\r
1043   }\r
1044   else if((fEventClass=="Multiplicity")&&(gAnalysisLevel == "MC")) {\r
1045     AliMCEvent* gMCEvent = dynamic_cast<AliMCEvent*>(event);\r
1046     //Calculating the multiplicity as the number of charged primaries\r
1047     //within \pm 0.8 in eta and pT > 0.1 GeV/c\r
1048     for(Int_t iParticle = 0; iParticle < gMCEvent->GetNumberOfPrimaries(); iParticle++) {\r
1049       AliMCParticle* track = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iParticle));\r
1050       if (!track) {\r
1051         AliError(Form("Could not receive particle %d", iParticle));\r
1052         continue;\r
1053       }\r
1054       \r
1055       //exclude non stable particles\r
1056       if(!(gMCEvent->IsPhysicalPrimary(iParticle))) continue;\r
1057       if(track->Pt() < 0.1)  continue;\r
1058 \r
1059       //++++++++++++++++\r
1060       if (fMultiplicityEstimator == "V0M") {\r
1061         if((track->Eta() > 5.1 || track->Eta() < 2.8)&&(track->Eta() < -3.7 || track->Eta() > -1.7)) \r
1062         continue;}\r
1063       else if (fMultiplicityEstimator == "V0A") {\r
1064         if(track->Eta() > 5.1 || track->Eta() < 2.8)  continue;}\r
1065       else if (fMultiplicityEstimator == "V0C") {\r
1066         if(track->Eta() > -1.7 || track->Eta() < -3.7)  continue;}\r
1067       else if (fMultiplicityEstimator == "TPC") {\r
1068         if(track->Eta() > 0.9 || track->Eta() < -0.9)  continue;}\r
1069       else{\r
1070         if(track->Eta() < fEtaMin || track->Eta() > fEtaMax)  continue;\r
1071       }\r
1072       //++++++++++++++++\r
1073 \r
1074       if(track->Charge() == 0) continue;\r
1075 \r
1076       gMultiplicity += 1;\r
1077     }//loop over primaries\r
1078     fHistMultiplicity->Fill(gMultiplicity);\r
1079   }//MC mode & multiplicity class\r
1080   \r
1081   Double_t lReturnVal = -100;\r
1082   if(fEventClass=="Multiplicity"){\r
1083     lReturnVal = gMultiplicity;\r
1084   }else if(fEventClass=="Centrality"){\r
1085     lReturnVal = gCentrality;\r
1086   }\r
1087   return lReturnVal;\r
1088 }\r
1089 \r
1090 //________________________________________________________________________\r
1091 Double_t AliAnalysisTaskBFPsi::GetReferenceMultiplicityFromAOD(AliVEvent *event){\r
1092   //Function that returns the reference multiplicity from AODs (data or reco MC)\r
1093   //Different ref. mult. implemented: V0M, V0A, V0C, TPC\r
1094   Double_t gRefMultiplicity = 0., gRefMultiplicityTPC = 0., gRefMultiplicityVZERO = 0.;\r
1095 \r
1096   // Loop over tracks in event\r
1097   for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r
1098     AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));\r
1099     if (!aodTrack) {\r
1100       AliError(Form("Could not receive track %d", iTracks));\r
1101       continue;\r
1102     }\r
1103     \r
1104     // AOD track cuts\r
1105     if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;\r
1106             \r
1107     if(aodTrack->Charge() == 0) continue;\r
1108     // Kinematics cuts from ESD track cuts\r
1109     if( aodTrack->Pt() < fPtMin || aodTrack->Pt() > fPtMax)      continue;\r
1110     if( aodTrack->Eta() < fEtaMin || aodTrack->Eta() > fEtaMax)  continue;\r
1111     \r
1112     //=================PID (so far only for electron rejection)==========================//\r
1113     if(fElectronRejection) {\r
1114       // get the electron nsigma\r
1115       Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));\r
1116       \r
1117       // check only for given momentum range\r
1118       if( aodTrack->Pt() > fElectronRejectionMinPt && aodTrack->Pt() < fElectronRejectionMaxPt ){\r
1119         //look only at electron nsigma\r
1120         if(!fElectronOnlyRejection) {\r
1121           //Make the decision based on the n-sigma of electrons\r
1122           if(nSigma < fElectronRejectionNSigma) continue;\r
1123         }\r
1124         else {\r
1125           Double_t nSigmaPions   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));\r
1126           Double_t nSigmaKaons   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));\r
1127           Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));\r
1128           \r
1129           //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)\r
1130           if(nSigma < fElectronRejectionNSigma\r
1131              && nSigmaPions   > fElectronRejectionNSigma\r
1132              && nSigmaKaons   > fElectronRejectionNSigma\r
1133              && nSigmaProtons > fElectronRejectionNSigma ) continue;\r
1134         }\r
1135       }\r
1136     }//electron rejection\r
1137     \r
1138     gRefMultiplicityTPC += 1.0;\r
1139   }// track loop\r
1140   \r
1141   //VZERO segmentation in two detectors (0-31: VZERO-C, 32-63: VZERO-A)\r
1142   for(Int_t i = 0; i < 64; i++) {\r
1143     fHistVZEROSignal->Fill(i,event->GetVZEROEqMultiplicity(i));\r
1144     \r
1145     if(fMultiplicityEstimator == "V0C") {\r
1146       if(i > 31) continue;}\r
1147     else if(fMultiplicityEstimator == "V0A") {\r
1148       if(i < 32) continue;}\r
1149     \r
1150     gRefMultiplicityVZERO += event->GetVZEROEqMultiplicity(i);\r
1151   }//loop over PMTs\r
1152   \r
1153   if(fDebugLevel) \r
1154     Printf("VZERO multiplicity: %.0f - TPC multiplicity: %.0f",gRefMultiplicityVZERO,gRefMultiplicityTPC);\r
1155 \r
1156   fHistTPCvsVZEROMultiplicity->Fill(gRefMultiplicityVZERO,gRefMultiplicityTPC);\r
1157 \r
1158   if(fMultiplicityEstimator == "TPC") \r
1159     gRefMultiplicity = gRefMultiplicityTPC;\r
1160   else if((fMultiplicityEstimator == "V0M")||\r
1161           (fMultiplicityEstimator == "V0A")||\r
1162           (fMultiplicityEstimator == "V0C")) \r
1163     gRefMultiplicity = gRefMultiplicityVZERO;\r
1164   \r
1165   return gRefMultiplicity;\r
1166 }\r
1167 \r
1168 //________________________________________________________________________\r
1169 Double_t AliAnalysisTaskBFPsi::GetEventPlane(AliVEvent *event){\r
1170   // Get the event plane\r
1171 \r
1172   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
1173 \r
1174   Float_t gVZEROEventPlane    = -10.;\r
1175   Float_t gReactionPlane      = -10.;\r
1176   Double_t qxTot = 0.0, qyTot = 0.0;\r
1177 \r
1178   //MC: from reaction plane\r
1179   if(gAnalysisLevel == "MC"){\r
1180     if(!event) {\r
1181       AliError("mcEvent not available");\r
1182       return 0x0;\r
1183     }\r
1184 \r
1185     AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);\r
1186     if(gMCEvent){\r
1187       AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader());    \r
1188       if (headerH) {\r
1189         gReactionPlane = headerH->ReactionPlaneAngle();\r
1190         //gReactionPlane *= TMath::RadToDeg();\r
1191       }//MC header\r
1192     }//MC event cast\r
1193   }//MC\r
1194   \r
1195   // AOD,ESD,ESDMC: from VZERO Event Plane\r
1196   else{\r
1197    \r
1198     AliEventplane *ep = event->GetEventplane();\r
1199     if(ep){ \r
1200       gVZEROEventPlane = ep->CalculateVZEROEventPlane(event,10,2,qxTot,qyTot);\r
1201       if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();\r
1202       //gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();\r
1203       gReactionPlane = gVZEROEventPlane;\r
1204     }\r
1205   }//AOD,ESD,ESDMC\r
1206 \r
1207   return gReactionPlane;\r
1208 }\r
1209 \r
1210 //________________________________________________________________________\r
1211 Double_t AliAnalysisTaskBFPsi::GetTrackbyTrackCorrectionMatrix( Double_t vEta, \r
1212                                                                 Double_t vPhi, \r
1213                                                                 Double_t vPt, \r
1214                                                                 Short_t vCharge, \r
1215                                                                 Double_t gCentrality) {\r
1216   // -- Get efficiency correction of particle dependent on (eta, phi, pt, charge, centrality) \r
1217 \r
1218   Double_t correction = 1.;\r
1219   Int_t binEta = 0, binPt = 0, binPhi = 0;\r
1220 \r
1221   //Printf("EtaMAx: %lf - EtaMin: %lf - EtaBin: %lf", fEtaMaxForCorrections,fEtaMinForCorrections,fEtaBinForCorrections);\r
1222   if(fEtaBinForCorrections != 0) {\r
1223     Double_t widthEta = (fEtaMaxForCorrections - fEtaMinForCorrections)/fEtaBinForCorrections;\r
1224     if(fEtaMaxForCorrections != fEtaMinForCorrections) \r
1225       binEta = (Int_t)((vEta-fEtaMinForCorrections)/widthEta)+1;\r
1226   }\r
1227 \r
1228   if(fPtBinForCorrections != 0) {\r
1229     Double_t widthPt = (fPtMaxForCorrections - fPtMinForCorrections)/fPtBinForCorrections;\r
1230     if(fPtMaxForCorrections != fPtMinForCorrections) \r
1231       binPt = (Int_t)((vPt-fPtMinForCorrections)/widthPt) + 1;\r
1232   }\r
1233  \r
1234   if(fPhiBinForCorrections != 0) {\r
1235     Double_t widthPhi = (fPhiMaxForCorrections - fPhiMinForCorrections)/fPhiBinForCorrections;\r
1236     if(fPhiMaxForCorrections != fPhiMinForCorrections) \r
1237       binPhi = (Int_t)((vPhi-fPhiMinForCorrections)/widthPhi)+ 1;\r
1238   }\r
1239 \r
1240   Int_t gCentralityInt = -1;\r
1241   for (Int_t i=0; i<fCentralityArrayBinsForCorrections-1; i++){\r
1242     if((fCentralityArrayForCorrections[i] <= gCentrality)&&(gCentrality <= fCentralityArrayForCorrections[i+1])){\r
1243       gCentralityInt = i;\r
1244       break;\r
1245     }\r
1246   }  \r
1247 \r
1248   // centrality not in array --> no correction\r
1249   if(gCentralityInt < 0){\r
1250     correction = 1.;\r
1251   }\r
1252   else{\r
1253     \r
1254     //Printf("//=============CENTRALITY=============// %d:",gCentralityInt);\r
1255     \r
1256     if(fHistCorrectionPlus[gCentralityInt]){\r
1257       if (vCharge > 0) {\r
1258         correction = fHistCorrectionPlus[gCentralityInt]->GetBinContent(fHistCorrectionPlus[gCentralityInt]->GetBin(binEta, binPt, binPhi));\r
1259         //Printf("CORRECTIONplus: %.2f | Centrality %d",correction,gCentralityInt);  \r
1260       }\r
1261       if (vCharge < 0) {\r
1262         correction = fHistCorrectionMinus[gCentralityInt]->GetBinContent(fHistCorrectionMinus[gCentralityInt]->GetBin(binEta, binPt, binPhi));\r
1263         //Printf("CORRECTIONminus: %.2f | Centrality %d",correction,gCentralityInt);\r
1264       }\r
1265     }\r
1266     else {\r
1267       correction = 1.;\r
1268     }\r
1269   }//centrality in array\r
1270   \r
1271   if (correction == 0.) { \r
1272     AliError(Form("Should not happen : bin content = 0. >> eta: %.2f | phi : %.2f | pt : %.2f | cent %d",vEta, vPhi, vPt, gCentralityInt)); \r
1273     correction = 1.; \r
1274   } \r
1275   \r
1276   return correction;\r
1277 }\r
1278 \r
1279 //________________________________________________________________________\r
1280 TObjArray* AliAnalysisTaskBFPsi::GetAcceptedTracks(AliVEvent *event, Double_t gCentrality, Double_t gReactionPlane){\r
1281   // Returns TObjArray with tracks after all track cuts (only for AOD!)\r
1282   // Fills QA histograms\r
1283 \r
1284   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
1285 \r
1286   //output TObjArray holding all good tracks\r
1287   TObjArray* tracksAccepted = new TObjArray;\r
1288   tracksAccepted->SetOwner(kTRUE);\r
1289 \r
1290   Short_t vCharge;\r
1291   Double_t vEta;\r
1292   Double_t vY;\r
1293   Double_t vPhi;\r
1294   Double_t vPt;\r
1295 \r
1296 \r
1297   if(gAnalysisLevel == "AOD") { // handling of TPC only tracks different in AOD and ESD\r
1298     // Loop over tracks in event\r
1299     for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r
1300       AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));\r
1301       if (!aodTrack) {\r
1302         AliError(Form("Could not receive track %d", iTracks));\r
1303         continue;\r
1304       }\r
1305       \r
1306       // AOD track cuts\r
1307       \r
1308       // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C\r
1309       // take only TPC only tracks \r
1310       //fHistTrackStats->Fill(aodTrack->GetFilterMap());\r
1311       for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){\r
1312         fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));\r
1313       }\r
1314       if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;\r
1315       \r
1316       vCharge = aodTrack->Charge();\r
1317       vEta    = aodTrack->Eta();\r
1318       vY      = aodTrack->Y();\r
1319       vPhi    = aodTrack->Phi();// * TMath::RadToDeg();\r
1320       vPt     = aodTrack->Pt();\r
1321       \r
1322       //===========================PID (so far only for electron rejection)===============================//                \r
1323       if(fElectronRejection) {\r
1324 \r
1325         // get the electron nsigma\r
1326         Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));\r
1327         \r
1328         //Fill QA before the PID\r
1329         fHistdEdxVsPTPCbeforePID -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());\r
1330         fHistNSigmaTPCvsPtbeforePID -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); \r
1331         //end of QA-before pid\r
1332         \r
1333         // check only for given momentum range\r
1334         if( vPt > fElectronRejectionMinPt && vPt < fElectronRejectionMaxPt ){\r
1335                   \r
1336           //look only at electron nsigma\r
1337           if(!fElectronOnlyRejection){\r
1338             \r
1339             //Make the decision based on the n-sigma of electrons\r
1340             if(nSigma < fElectronRejectionNSigma) continue;\r
1341           }\r
1342           else{\r
1343             \r
1344             Double_t nSigmaPions   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));\r
1345             Double_t nSigmaKaons   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));\r
1346             Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));\r
1347             \r
1348             //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)\r
1349             if(nSigma < fElectronRejectionNSigma\r
1350                && nSigmaPions   > fElectronRejectionNSigma\r
1351                && nSigmaKaons   > fElectronRejectionNSigma\r
1352                && nSigmaProtons > fElectronRejectionNSigma ) continue;\r
1353           }\r
1354         }\r
1355   \r
1356         //Fill QA after the PID\r
1357         fHistdEdxVsPTPCafterPID -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());\r
1358         fHistNSigmaTPCvsPtafterPID -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); \r
1359         \r
1360       }\r
1361       //===========================end of PID (so far only for electron rejection)===============================//\r
1362 \r
1363       Float_t dcaXY = aodTrack->DCA();      // this is the DCA from global track (not exactly what is cut on)\r
1364       Float_t dcaZ  = aodTrack->ZAtDCA();   // this is the DCA from global track (not exactly what is cut on)\r
1365       \r
1366       \r
1367       // Kinematics cuts from ESD track cuts\r
1368       if( vPt < fPtMin || vPt > fPtMax)      continue;\r
1369       if( vEta < fEtaMin || vEta > fEtaMax)  continue;\r
1370       \r
1371       // Extra DCA cuts (for systematic studies [!= -1])\r
1372       if( fDCAxyCut != -1 && fDCAzCut != -1){\r
1373         if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){\r
1374           continue;  // 2D cut\r
1375         }\r
1376       }\r
1377       \r
1378       // Extra TPC cuts (for systematic studies [!= -1])\r
1379       if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){\r
1380         continue;\r
1381       }\r
1382       if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){\r
1383         continue;\r
1384       }\r
1385       \r
1386       // fill QA histograms\r
1387       fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());\r
1388       fHistDCA->Fill(dcaZ,dcaXY);\r
1389       fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality);\r
1390       fHistPt->Fill(vPt,gCentrality);\r
1391       fHistEta->Fill(vEta,gCentrality);\r
1392       fHistRapidity->Fill(vY,gCentrality);\r
1393       if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
1394       else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
1395       fHistPhi->Fill(vPhi,gCentrality);\r
1396       if(vCharge > 0)      fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);                  \r
1397       else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r
1398       \r
1399       //=======================================correction\r
1400       Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  \r
1401       //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);\r
1402       \r
1403       // add the track to the TObjArray\r
1404       tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));  \r
1405     }//track loop\r
1406   }// AOD analysis\r
1407 \r
1408   //==============================================================================================================\r
1409   else if(gAnalysisLevel == "MCAOD") {\r
1410     \r
1411     AliMCEvent* mcEvent = MCEvent();\r
1412     if (!mcEvent) {\r
1413       AliError("ERROR: Could not retrieve MC event");\r
1414     }\r
1415     else{\r
1416       \r
1417       for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {\r
1418         AliAODMCParticle *aodTrack = (AliAODMCParticle*) mcEvent->GetTrack(iTracks); \r
1419         if (!aodTrack) {\r
1420           AliError(Form("ERROR: Could not receive track %d (mc loop)", iTracks));\r
1421           continue;\r
1422         }\r
1423         \r
1424         if(!aodTrack->IsPhysicalPrimary()) continue;   \r
1425         \r
1426         vCharge = aodTrack->Charge();\r
1427         vEta    = aodTrack->Eta();\r
1428         vY      = aodTrack->Y();\r
1429         vPhi    = aodTrack->Phi();// * TMath::RadToDeg();\r
1430         vPt     = aodTrack->Pt();\r
1431         \r
1432         // Kinematics cuts from ESD track cuts\r
1433         if( vPt < fPtMin || vPt > fPtMax)      continue;\r
1434         if( vEta < fEtaMin || vEta > fEtaMax)  continue;\r
1435         \r
1436         // Remove neutral tracks\r
1437         if( vCharge == 0 ) continue;\r
1438         \r
1439         //Exclude resonances\r
1440         if(fExcludeResonancesInMC) {\r
1441           \r
1442           Bool_t kExcludeParticle = kFALSE;\r
1443           Int_t gMotherIndex = aodTrack->GetMother();\r
1444           if(gMotherIndex != -1) {\r
1445             AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));\r
1446             if(motherTrack) {\r
1447               Int_t pdgCodeOfMother = motherTrack->GetPdgCode();\r
1448               //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {\r
1449               //if(pdgCodeOfMother == 113) {\r
1450               if(pdgCodeOfMother == 113  // rho0\r
1451                  || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+\r
1452                  // || pdgCodeOfMother == 221  // eta\r
1453                  // || pdgCodeOfMother == 331  // eta'\r
1454                  // || pdgCodeOfMother == 223  // omega\r
1455                  // || pdgCodeOfMother == 333  // phi\r
1456                  || pdgCodeOfMother == 311  || pdgCodeOfMother == -311 // K0\r
1457                  // || pdgCodeOfMother == 313  || pdgCodeOfMother == -313 // K0*\r
1458                  // || pdgCodeOfMother == 323  || pdgCodeOfMother == -323 // K+*\r
1459                  || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda\r
1460                  || pdgCodeOfMother == 111  // pi0 Dalitz\r
1461                  || pdgCodeOfMother == 22  // photon\r
1462                  ) {\r
1463                 kExcludeParticle = kTRUE;\r
1464               }\r
1465             }\r
1466           }\r
1467           \r
1468           //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi\r
1469           if(kExcludeParticle) continue;\r
1470         }\r
1471 \r
1472         //Exclude electrons with PDG\r
1473         if(fExcludeElectronsInMC) {\r
1474           \r
1475           if(TMath::Abs(aodTrack->GetPdgCode()) == 11) continue;\r
1476           \r
1477         }\r
1478         \r
1479         // fill QA histograms\r
1480         fHistPt->Fill(vPt,gCentrality);\r
1481         fHistEta->Fill(vEta,gCentrality);\r
1482         fHistRapidity->Fill(vY,gCentrality);\r
1483         if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
1484         else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
1485         fHistPhi->Fill(vPhi,gCentrality);\r
1486         if(vCharge > 0)      fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);                \r
1487         else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r
1488         \r
1489         //=======================================correction\r
1490         Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  \r
1491         //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);   \r
1492         \r
1493         // add the track to the TObjArray\r
1494         tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));  \r
1495       }//aodTracks\r
1496     }//MC event\r
1497   }//MCAOD\r
1498   //==============================================================================================================\r
1499 \r
1500   //==============================================================================================================\r
1501   else if(gAnalysisLevel == "MCAODrec") {\r
1502     \r
1503     /* fAOD = dynamic_cast<AliAODEvent*>(InputEvent());\r
1504     if (!fAOD) {\r
1505       printf("ERROR: fAOD not available\n");\r
1506       return;\r
1507       }*/\r
1508     \r
1509     fArrayMC = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName())); \r
1510     if (!fArrayMC) { \r
1511        AliError("No array of MC particles found !!!");\r
1512     }\r
1513 \r
1514     AliMCEvent* mcEvent = MCEvent();\r
1515     if (!mcEvent) {\r
1516        AliError("ERROR: Could not retrieve MC event");\r
1517        return tracksAccepted;\r
1518     }\r
1519      \r
1520     for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r
1521       AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));\r
1522       if (!aodTrack) {\r
1523         AliError(Form("Could not receive track %d", iTracks));\r
1524         continue;\r
1525       }\r
1526       \r
1527       for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){\r
1528         fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));\r
1529       }\r
1530       if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;\r
1531       \r
1532       vCharge = aodTrack->Charge();\r
1533       vEta    = aodTrack->Eta();\r
1534       vY      = aodTrack->Y();\r
1535       vPhi    = aodTrack->Phi();// * TMath::RadToDeg();\r
1536       vPt     = aodTrack->Pt();\r
1537       \r
1538       //===========================use MC information for Kinematics===============================//               \r
1539       if(fUseMCforKinematics){\r
1540 \r
1541         Int_t label = TMath::Abs(aodTrack->GetLabel());\r
1542         AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);\r
1543 \r
1544         if(AODmcTrack){\r
1545           vCharge = AODmcTrack->Charge();\r
1546           vEta    = AODmcTrack->Eta();\r
1547           vY      = AODmcTrack->Y();\r
1548           vPhi    = AODmcTrack->Phi();// * TMath::RadToDeg();\r
1549           vPt     = AODmcTrack->Pt();\r
1550         }\r
1551         else{\r
1552           AliDebug(1, "no MC particle for this track"); \r
1553           continue;\r
1554         }\r
1555       }\r
1556       //===========================end of use MC information for Kinematics========================//               \r
1557 \r
1558 \r
1559       //===========================PID (so far only for electron rejection)===============================//                \r
1560       if(fElectronRejection) {\r
1561 \r
1562         // get the electron nsigma\r
1563         Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));\r
1564 \r
1565         //Fill QA before the PID\r
1566         fHistdEdxVsPTPCbeforePID -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());\r
1567         fHistNSigmaTPCvsPtbeforePID -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); \r
1568         //end of QA-before pid\r
1569         \r
1570         // check only for given momentum range\r
1571         if( vPt > fElectronRejectionMinPt && vPt < fElectronRejectionMaxPt ){\r
1572                   \r
1573           //look only at electron nsigma\r
1574           if(!fElectronOnlyRejection){\r
1575             \r
1576             //Make the decision based on the n-sigma of electrons\r
1577             if(nSigma < fElectronRejectionNSigma) continue;\r
1578           }\r
1579           else{\r
1580             \r
1581             Double_t nSigmaPions   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));\r
1582             Double_t nSigmaKaons   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));\r
1583             Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));\r
1584             \r
1585             //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)\r
1586             if(nSigma < fElectronRejectionNSigma\r
1587                && nSigmaPions   > fElectronRejectionNSigma\r
1588                && nSigmaKaons   > fElectronRejectionNSigma\r
1589                && nSigmaProtons > fElectronRejectionNSigma ) continue;\r
1590           }\r
1591         }\r
1592   \r
1593         //Fill QA after the PID\r
1594         fHistdEdxVsPTPCafterPID -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());\r
1595         fHistNSigmaTPCvsPtafterPID -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); \r
1596         \r
1597       }\r
1598       //===========================end of PID (so far only for electron rejection)===============================//\r
1599       \r
1600       Float_t dcaXY = aodTrack->DCA();      // this is the DCA from global track (not exactly what is cut on)\r
1601       Float_t dcaZ  = aodTrack->ZAtDCA();   // this is the DCA from global track (not exactly what is cut on)\r
1602             \r
1603       // Kinematics cuts from ESD track cuts\r
1604       if( vPt < fPtMin || vPt > fPtMax)      continue;\r
1605       if( vEta < fEtaMin || vEta > fEtaMax)  continue;\r
1606       \r
1607       // Extra DCA cuts (for systematic studies [!= -1])\r
1608       if( fDCAxyCut != -1 && fDCAzCut != -1){\r
1609         if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){\r
1610           continue;  // 2D cut\r
1611         }\r
1612       }\r
1613       \r
1614       // Extra TPC cuts (for systematic studies [!= -1])\r
1615       if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){\r
1616         continue;\r
1617       }\r
1618       if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){\r
1619         continue;\r
1620       }\r
1621 \r
1622       //Exclude resonances\r
1623       if(fExcludeResonancesInMC) {\r
1624         \r
1625         Bool_t kExcludeParticle = kFALSE;\r
1626 \r
1627         Int_t label = TMath::Abs(aodTrack->GetLabel());\r
1628         AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);\r
1629       \r
1630         if (AODmcTrack){ \r
1631           //if (AODmcTrack->IsPhysicalPrimary()){\r
1632           \r
1633           Int_t gMotherIndex = AODmcTrack->GetMother();\r
1634           if(gMotherIndex != -1) {\r
1635             AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));\r
1636             if(motherTrack) {\r
1637               Int_t pdgCodeOfMother = motherTrack->GetPdgCode();\r
1638               if(pdgCodeOfMother == 113  // rho0\r
1639                  || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+\r
1640                  // || pdgCodeOfMother == 221  // eta\r
1641                  // || pdgCodeOfMother == 331  // eta'\r
1642                  // || pdgCodeOfMother == 223  // omega\r
1643                  // || pdgCodeOfMother == 333  // phi\r
1644                  || pdgCodeOfMother == 311  || pdgCodeOfMother == -311 // K0\r
1645                  // || pdgCodeOfMother == 313  || pdgCodeOfMother == -313 // K0*\r
1646                  // || pdgCodeOfMother == 323  || pdgCodeOfMother == -323 // K+*\r
1647                  || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda\r
1648                  || pdgCodeOfMother == 111  // pi0 Dalitz\r
1649                  || pdgCodeOfMother == 22   // photon\r
1650                  ) {\r
1651                 kExcludeParticle = kTRUE;\r
1652               }\r
1653             }\r
1654           }\r
1655         }       \r
1656         //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi\r
1657         if(kExcludeParticle) continue;\r
1658       }\r
1659 \r
1660       //Exclude electrons with PDG\r
1661       if(fExcludeElectronsInMC) {\r
1662         \r
1663         Int_t label = TMath::Abs(aodTrack->GetLabel());\r
1664         AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);\r
1665         \r
1666         if (AODmcTrack){ \r
1667           if(TMath::Abs(AODmcTrack->GetPdgCode()) == 11) continue;\r
1668         }\r
1669       }\r
1670       \r
1671       // fill QA histograms\r
1672       fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());\r
1673       fHistDCA->Fill(dcaZ,dcaXY);\r
1674       fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality);\r
1675       fHistPt->Fill(vPt,gCentrality);\r
1676       fHistEta->Fill(vEta,gCentrality);\r
1677       fHistRapidity->Fill(vY,gCentrality);\r
1678       if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
1679       else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
1680       fHistPhi->Fill(vPhi,gCentrality);\r
1681       if(vCharge > 0)      fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);                  \r
1682       else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r
1683       \r
1684       //=======================================correction\r
1685       Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  \r
1686       //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);\r
1687       \r
1688       // add the track to the TObjArray\r
1689       tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));  \r
1690     }//track loop\r
1691   }//MCAODrec\r
1692   //==============================================================================================================\r
1693 \r
1694   else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD") { // handling of TPC only tracks different in AOD and ESD\r
1695 \r
1696     AliESDtrack *trackTPC   = NULL;\r
1697     AliMCParticle *mcTrack   = NULL;\r
1698 \r
1699     AliMCEvent*  mcEvent     = NULL;\r
1700     \r
1701     // for MC ESDs use also MC information\r
1702     if(gAnalysisLevel == "MCESD"){\r
1703       mcEvent = MCEvent(); \r
1704       if (!mcEvent) {\r
1705         AliError("mcEvent not available");\r
1706         return tracksAccepted;\r
1707       }\r
1708     }\r
1709     \r
1710     // Loop over tracks in event\r
1711     for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r
1712       AliESDtrack* track = dynamic_cast<AliESDtrack *>(event->GetTrack(iTracks));\r
1713       if (!track) {\r
1714         AliError(Form("Could not receive track %d", iTracks));\r
1715         continue;\r
1716       } \r
1717 \r
1718       // for MC ESDs use also MC information --> MC track not used further???\r
1719       if(gAnalysisLevel == "MCESD"){\r
1720         Int_t label = TMath::Abs(track->GetLabel());\r
1721         if(label > mcEvent->GetNumberOfTracks()) continue;\r
1722         if(label > mcEvent->GetNumberOfPrimaries()) continue;\r
1723         \r
1724         mcTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label));\r
1725         if(!mcTrack) continue;\r
1726       }\r
1727 \r
1728       // take only TPC only tracks\r
1729       trackTPC   = new AliESDtrack();\r
1730       if(!track->FillTPCOnlyTrack(*trackTPC)) continue;\r
1731       \r
1732       //ESD track cuts\r
1733       if(fESDtrackCuts) \r
1734         if(!fESDtrackCuts->AcceptTrack(trackTPC)) continue;\r
1735       \r
1736       // fill QA histograms\r
1737       Float_t b[2];\r
1738       Float_t bCov[3];\r
1739       trackTPC->GetImpactParameters(b,bCov);\r
1740       if (bCov[0]<=0 || bCov[2]<=0) {\r
1741         AliDebug(1, "Estimated b resolution lower or equal zero!");\r
1742         bCov[0]=0; bCov[2]=0;\r
1743       }\r
1744       \r
1745       Int_t nClustersTPC = -1;\r
1746       nClustersTPC = trackTPC->GetTPCNclsIter1();   // TPC standalone\r
1747       //nClustersTPC = track->GetTPCclusters(0);   // global track\r
1748       Float_t chi2PerClusterTPC = -1;\r
1749       if (nClustersTPC!=0) {\r
1750         chi2PerClusterTPC = trackTPC->GetTPCchi2Iter1()/Float_t(nClustersTPC);      // TPC standalone\r
1751         //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);     // global track\r
1752       }\r
1753       \r
1754       //===========================PID===============================//             \r
1755       if(fUsePID) {\r
1756         Double_t prob[AliPID::kSPECIES]={0.};\r
1757         Double_t probTPC[AliPID::kSPECIES]={0.};\r
1758         Double_t probTOF[AliPID::kSPECIES]={0.};\r
1759         Double_t probTPCTOF[AliPID::kSPECIES]={0.};\r
1760         \r
1761         Double_t nSigma = 0.;\r
1762         UInt_t detUsedTPC = 0;\r
1763         UInt_t detUsedTOF = 0;\r
1764         UInt_t detUsedTPCTOF = 0;\r
1765         \r
1766         //Decide what detector configuration we want to use\r
1767         switch(fPidDetectorConfig) {\r
1768         case kTPCpid:\r
1769           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);\r
1770           nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest));\r
1771           detUsedTPC = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPC);\r
1772           for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
1773             prob[iSpecies] = probTPC[iSpecies];\r
1774           break;\r
1775         case kTOFpid:\r
1776           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);\r
1777           nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest));\r
1778           detUsedTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTOF);\r
1779           for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
1780             prob[iSpecies] = probTOF[iSpecies];\r
1781           break;\r
1782         case kTPCTOF:\r
1783           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);\r
1784           detUsedTPCTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPCTOF);\r
1785           for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
1786             prob[iSpecies] = probTPCTOF[iSpecies];\r
1787           break;\r
1788         default:\r
1789           break;\r
1790         }//end switch: define detector mask\r
1791         \r
1792         //Filling the PID QA\r
1793         Double_t tofTime = -999., length = 999., tof = -999.;\r
1794         Double_t c = TMath::C()*1.E-9;// m/ns\r
1795         Double_t beta = -999.;\r
1796         Double_t  nSigmaTOFForParticleOfInterest = -999.;\r
1797         if ( (track->IsOn(AliESDtrack::kTOFin)) &&\r
1798              (track->IsOn(AliESDtrack::kTIME))  ) { \r
1799           tofTime = track->GetTOFsignal();//in ps\r
1800           length = track->GetIntegratedLength();\r
1801           tof = tofTime*1E-3; // ns     \r
1802           \r
1803           if (tof <= 0) {\r
1804             //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");\r
1805             continue;\r
1806           }\r
1807           if (length <= 0){\r
1808             //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");\r
1809             continue;\r
1810           }\r
1811           \r
1812           length = length*0.01; // in meters\r
1813           tof = tof*c;\r
1814           beta = length/tof;\r
1815           \r
1816           nSigmaTOFForParticleOfInterest = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest);\r
1817           fHistBetavsPTOFbeforePID ->Fill(track->P()*track->Charge(),beta);\r
1818           fHistProbTOFvsPtbeforePID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);\r
1819           fHistNSigmaTOFvsPtbeforePID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);\r
1820         }//TOF signal \r
1821         \r
1822         \r
1823         Double_t  nSigmaTPCForParticleOfInterest = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest);\r
1824         fHistdEdxVsPTPCbeforePID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());\r
1825         fHistProbTPCvsPtbeforePID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); \r
1826         fHistNSigmaTPCvsPtbeforePID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); \r
1827         fHistProbTPCTOFvsPtbeforePID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);\r
1828         //end of QA-before pid\r
1829         \r
1830         if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {\r
1831           //Make the decision based on the n-sigma\r
1832           if(fUsePIDnSigma) {\r
1833             if(nSigma > fPIDNSigma) continue;}\r
1834           \r
1835           //Make the decision based on the bayesian\r
1836           else if(fUsePIDPropabilities) {\r
1837             if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;\r
1838             if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue;      \r
1839           }\r
1840           \r
1841           //Fill QA after the PID\r
1842           fHistBetavsPTOFafterPID ->Fill(track->P()*track->Charge(),beta);\r
1843           fHistProbTOFvsPtafterPID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);\r
1844           fHistNSigmaTOFvsPtafterPID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);\r
1845           \r
1846           fHistdEdxVsPTPCafterPID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());\r
1847           fHistProbTPCvsPtafterPID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); \r
1848           fHistProbTPCTOFvsPtafterPID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);\r
1849           fHistNSigmaTPCvsPtafterPID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); \r
1850         }\r
1851       }\r
1852       //===========================PID===============================//\r
1853       vCharge = trackTPC->Charge();\r
1854       vY      = trackTPC->Y();\r
1855       vEta    = trackTPC->Eta();\r
1856       vPhi    = trackTPC->Phi();// * TMath::RadToDeg();\r
1857       vPt     = trackTPC->Pt();\r
1858       fHistClus->Fill(trackTPC->GetITSclusters(0),nClustersTPC);\r
1859       fHistDCA->Fill(b[1],b[0]);\r
1860       fHistChi2->Fill(chi2PerClusterTPC,gCentrality);\r
1861       fHistPt->Fill(vPt,gCentrality);\r
1862       fHistEta->Fill(vEta,gCentrality);\r
1863       fHistPhi->Fill(vPhi,gCentrality);\r
1864       if(vCharge > 0)      fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);\r
1865       else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r
1866       fHistRapidity->Fill(vY,gCentrality);\r
1867       if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
1868       else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
1869       \r
1870       //=======================================correction\r
1871       Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  \r
1872       //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);\r
1873 \r
1874       // add the track to the TObjArray\r
1875       tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));   \r
1876 \r
1877       delete trackTPC;\r
1878     }//track loop\r
1879   }// ESD analysis\r
1880 \r
1881   else if(gAnalysisLevel == "MC"){\r
1882     if(!event) {\r
1883       AliError("mcEvent not available");\r
1884       return 0x0;\r
1885     }\r
1886 \r
1887     AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);\r
1888     if(gMCEvent) {\r
1889       // Loop over tracks in event\r
1890       for (Int_t iTracks = 0; iTracks < gMCEvent->GetNumberOfPrimaries(); iTracks++) {\r
1891         AliMCParticle* track = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iTracks));\r
1892         if (!track) {\r
1893           AliError(Form("Could not receive particle %d", iTracks));\r
1894           continue;\r
1895         }\r
1896         \r
1897         //exclude non stable particles\r
1898         if(!(gMCEvent->IsPhysicalPrimary(iTracks))) continue;\r
1899         \r
1900         vCharge = track->Charge();\r
1901         vEta    = track->Eta();\r
1902         vPt     = track->Pt();\r
1903         vY      = track->Y();\r
1904         \r
1905         if( vPt < fPtMin || vPt > fPtMax)      \r
1906           continue;\r
1907         if (!fUsePID) {\r
1908           if( vEta < fEtaMin || vEta > fEtaMax)  continue;\r
1909         }\r
1910         else if (fUsePID){\r
1911           if( vY < fEtaMin || vY > fEtaMax)  continue;\r
1912         }\r
1913 \r
1914         // Remove neutral tracks\r
1915         if( vCharge == 0 ) continue;\r
1916         \r
1917         //analyze one set of particles\r
1918         if(fUseMCPdgCode) {\r
1919           TParticle *particle = track->Particle();\r
1920           if(!particle) continue;\r
1921           \r
1922           Int_t gPdgCode = particle->GetPdgCode();\r
1923           if(TMath::Abs(fPDGCodeToBeAnalyzed) != TMath::Abs(gPdgCode)) \r
1924             continue;\r
1925         }\r
1926         \r
1927         //Use the acceptance parameterization\r
1928         if(fAcceptanceParameterization) {\r
1929           Double_t gRandomNumber = gRandom->Rndm();\r
1930           if(gRandomNumber > fAcceptanceParameterization->Eval(track->Pt())) \r
1931             continue;\r
1932         }\r
1933         \r
1934         //Exclude resonances\r
1935         if(fExcludeResonancesInMC) {\r
1936           TParticle *particle = track->Particle();\r
1937           if(!particle) continue;\r
1938           \r
1939           Bool_t kExcludeParticle = kFALSE;\r
1940           Int_t gMotherIndex = particle->GetFirstMother();\r
1941           if(gMotherIndex != -1) {\r
1942             AliMCParticle* motherTrack = dynamic_cast<AliMCParticle *>(event->GetTrack(gMotherIndex));\r
1943             if(motherTrack) {\r
1944               TParticle *motherParticle = motherTrack->Particle();\r
1945               if(motherParticle) {\r
1946                 Int_t pdgCodeOfMother = motherParticle->GetPdgCode();\r
1947                 //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {\r
1948                 if(pdgCodeOfMother == 113  // rho0\r
1949                    || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+\r
1950                    // || pdgCodeOfMother == 221  // eta\r
1951                    // || pdgCodeOfMother == 331  // eta'\r
1952                    // || pdgCodeOfMother == 223  // omega\r
1953                    // || pdgCodeOfMother == 333  // phi\r
1954                    || pdgCodeOfMother == 311  || pdgCodeOfMother == -311 // K0\r
1955                    // || pdgCodeOfMother == 313  || pdgCodeOfMother == -313 // K0*\r
1956                    // || pdgCodeOfMother == 323  || pdgCodeOfMother == -323 // K+*\r
1957                    || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda\r
1958                    || pdgCodeOfMother == 111  // pi0 Dalitz\r
1959                    ) {\r
1960                   kExcludeParticle = kTRUE;\r
1961                 }\r
1962               }\r
1963             }\r
1964           }\r
1965           \r
1966           //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi\r
1967           if(kExcludeParticle) continue;\r
1968         }\r
1969 \r
1970         //Exclude electrons with PDG\r
1971         if(fExcludeElectronsInMC) {\r
1972           \r
1973           TParticle *particle = track->Particle();\r
1974           \r
1975           if (particle){ \r
1976             if(TMath::Abs(particle->GetPdgCode()) == 11) continue;\r
1977           }\r
1978         }\r
1979       \r
1980         vPhi    = track->Phi();\r
1981         //Printf("phi (before): %lf",vPhi);\r
1982         \r
1983         fHistPt->Fill(vPt,gCentrality);\r
1984         fHistEta->Fill(vEta,gCentrality);\r
1985         fHistPhi->Fill(vPhi,gCentrality);\r
1986         if(vCharge > 0)      fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);\r
1987         else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r
1988         //fHistPhi->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r
1989         fHistRapidity->Fill(vY,gCentrality);\r
1990         //if(vCharge > 0) fHistPhiPos->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r
1991         //else if(vCharge < 0) fHistPhiNeg->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r
1992         if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
1993         else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
1994         \r
1995         //Flow after burner\r
1996         if(fUseFlowAfterBurner) {\r
1997           Double_t precisionPhi = 0.001;\r
1998           Int_t maxNumberOfIterations = 100;\r
1999           \r
2000           Double_t phi0 = vPhi;\r
2001           Double_t gV2 = fDifferentialV2->Eval(vPt);\r
2002           \r
2003           for (Int_t j = 0; j < maxNumberOfIterations; j++) {\r
2004             Double_t phiprev = vPhi;\r
2005             Double_t fl = vPhi - phi0 + gV2*TMath::Sin(2.*(vPhi - gReactionPlane*TMath::DegToRad()));\r
2006             Double_t fp = 1.0 + 2.0*gV2*TMath::Cos(2.*(vPhi - gReactionPlane*TMath::DegToRad())); \r
2007             vPhi -= fl/fp;\r
2008             if (TMath::AreEqualAbs(phiprev,vPhi,precisionPhi)) break;\r
2009           }\r
2010           //Printf("phi (after): %lf\n",vPhi);\r
2011           Double_t vDeltaphiBefore = phi0 - gReactionPlane*TMath::DegToRad();\r
2012           if(vDeltaphiBefore < 0) vDeltaphiBefore += 2*TMath::Pi();\r
2013           fHistPhiBefore->Fill(vDeltaphiBefore,gCentrality);\r
2014           \r
2015           Double_t vDeltaphiAfter = vPhi - gReactionPlane*TMath::DegToRad();\r
2016           if(vDeltaphiAfter < 0) vDeltaphiAfter += 2*TMath::Pi();\r
2017           fHistPhiAfter->Fill(vDeltaphiAfter,gCentrality);\r
2018           \r
2019         }\r
2020         \r
2021         //vPhi *= TMath::RadToDeg();\r
2022         \r
2023       //=======================================correction\r
2024       Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  \r
2025       //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);\r
2026 \r
2027         tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction)); \r
2028       } //track loop\r
2029     }//MC event object\r
2030   }//MC\r
2031   \r
2032   return tracksAccepted;  \r
2033 }\r
2034 \r
2035 //________________________________________________________________________\r
2036 TObjArray* AliAnalysisTaskBFPsi::GetShuffledTracks(TObjArray *tracks, Double_t gCentrality){\r
2037   // Clones TObjArray and returns it with tracks after shuffling the charges\r
2038 \r
2039   TObjArray* tracksShuffled = new TObjArray;\r
2040   tracksShuffled->SetOwner(kTRUE);\r
2041 \r
2042   vector<Short_t> *chargeVector = new vector<Short_t>;   //original charge of accepted tracks \r
2043 \r
2044   for (Int_t i=0; i<tracks->GetEntriesFast(); i++)\r
2045   {\r
2046     AliVParticle* track = (AliVParticle*) tracks->At(i);\r
2047     chargeVector->push_back(track->Charge());\r
2048   }  \r
2049  \r
2050   random_shuffle(chargeVector->begin(), chargeVector->end());\r
2051   \r
2052   for(Int_t i = 0; i < tracks->GetEntriesFast(); i++){\r
2053     AliVParticle* track = (AliVParticle*) tracks->At(i);\r
2054     //==============================correction\r
2055     Double_t correction = GetTrackbyTrackCorrectionMatrix(track->Eta(), track->Phi(),track->Pt(), chargeVector->at(i), gCentrality);\r
2056     //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);\r
2057     tracksShuffled->Add(new AliBFBasicParticle(track->Eta(), track->Phi(), track->Pt(),chargeVector->at(i), correction));\r
2058   }\r
2059 \r
2060   delete chargeVector;\r
2061    \r
2062   return tracksShuffled;\r
2063 }\r
2064 \r
2065 \r
2066 //________________________________________________________________________\r
2067 void  AliAnalysisTaskBFPsi::FinishTaskOutput(){\r
2068   //Printf("END BF");\r
2069 \r
2070   if (!fBalance) {\r
2071     AliError("fBalance not available");\r
2072     return;\r
2073   }  \r
2074   if(fRunShuffling) {\r
2075     if (!fShuffledBalance) {\r
2076       AliError("fShuffledBalance not available");\r
2077       return;\r
2078     }\r
2079   }\r
2080 \r
2081 }\r
2082 \r
2083 //________________________________________________________________________\r
2084 void AliAnalysisTaskBFPsi::Terminate(Option_t *) {\r
2085   // Draw result to the screen\r
2086   // Called once at the end of the query\r
2087 \r
2088   // not implemented ...\r
2089 \r
2090 }\r
2091 \r