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