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("V0A"),\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-A multiplicity (a.u.)");\r
 336   else if(fMultiplicityEstimator == "V0C") \r
 337     fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO-C 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 lMultiplicityVar     = -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((lMultiplicityVar = IsEventAccepted(eventMain)) < 0){\r
 676     return;\r
 677   }\r
 678   \r
 679   // get the reaction plane\r
 680   if(fEventClass != "Multiplicity") {\r
 681     gReactionPlane = GetEventPlane(eventMain);\r
 682     fHistEventPlane->Fill(gReactionPlane,lMultiplicityVar);\r
 683     if(gReactionPlane < 0){\r
 684       return;\r
 685     }\r
 686   }\r
 687   \r
 688   // get the accepted tracks in main event\r
 689   TObjArray *tracksMain = GetAcceptedTracks(eventMain,lMultiplicityVar,gReactionPlane);\r
 690   gNumberOfAcceptedTracks = tracksMain->GetEntriesFast();\r
 691 \r
 692   //multiplicity cut (used in pp)\r
 693   fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedTracks,lMultiplicityVar);\r
 694 \r
 695   // store charges of all accepted tracks,shuffle and reassign(two extra loops)\r
 696   TObjArray* tracksShuffled = NULL;\r
 697   if(fRunShuffling){\r
 698     tracksShuffled = GetShuffledTracks(tracksMain,lMultiplicityVar);\r
 699   }\r
 700   \r
 701   // Event mixing \r
 702   if (fRunMixing)\r
 703     {\r
 704       // 1. First get an event pool corresponding in mult (cent) and\r
 705       //    zvertex to the current event. Once initialized, the pool\r
 706       //    should contain nMix (reduced) events. This routine does not\r
 707       //    pre-scan the chain. The first several events of every chain\r
 708       //    will be skipped until the needed pools are filled to the\r
 709       //    specified depth. If the pool categories are not too rare, this\r
 710       //    should not be a problem. If they are rare, you could lose`\r
 711       //    statistics.\r
 712       \r
 713       // 2. Collect the whole pool's content of tracks into one TObjArray\r
 714       //    (bgTracks), which is effectively a single background super-event.\r
 715       \r
 716       // 3. The reduced and bgTracks arrays must both be passed into\r
 717       //    FillCorrelations(). Also nMix should be passed in, so a weight\r
 718       //    of 1./nMix can be applied.\r
 719       \r
 720       AliEventPool* pool = fPoolMgr->GetEventPool(lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane);\r
 721       \r
 722       if (!pool){\r
 723         AliFatal(Form("No pool found for centrality = %f, zVtx = %f, psi = %f", lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane));\r
 724       }\r
 725       else{\r
 726         \r
 727         //pool->SetDebug(1);\r
 728 \r
 729         if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 || pool->GetCurrentNEvents() >= 5){ \r
 730           \r
 731           \r
 732           Int_t nMix = pool->GetCurrentNEvents();\r
 733           //cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << endl;\r
 734           \r
 735           //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(2);\r
 736           //((TH2F*) fListOfHistos->FindObject("mixedDist"))->Fill(centrality, pool->NTracksInPool());\r
 737           //if (pool->IsReady())\r
 738           //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(3);\r
 739           \r
 740           // Fill mixed-event histos here  \r
 741           for (Int_t jMix=0; jMix<nMix; jMix++) \r
 742             {\r
 743               TObjArray* tracksMixed = pool->GetEvent(jMix);\r
 744               fMixedBalance->CalculateBalance(gReactionPlane,tracksMain,tracksMixed,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());\r
 745             }\r
 746         }\r
 747         \r
 748         // Update the Event pool\r
 749         pool->UpdatePool(tracksMain);\r
 750         //pool->PrintInfo();\r
 751         \r
 752       }//pool NULL check  \r
 753     }//run mixing\r
 754   \r
 755   // calculate balance function\r
 756   fBalance->CalculateBalance(gReactionPlane,tracksMain,NULL,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());\r
 757   \r
 758   // calculate shuffled balance function\r
 759   if(fRunShuffling && tracksShuffled != NULL) {\r
 760     fShuffledBalance->CalculateBalance(gReactionPlane,tracksShuffled,NULL,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());\r
 761   }\r
 762 }      \r
 763 \r
 764 //________________________________________________________________________\r
 765 Double_t AliAnalysisTaskBFPsi::IsEventAccepted(AliVEvent *event){\r
 766   // Checks the Event cuts\r
 767   // Fills Event statistics histograms\r
 768   \r
 769   Bool_t isSelectedMain = kTRUE;\r
 770   Float_t gRefMultiplicity = -1.;\r
 771   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
 772 \r
 773   AliMCEvent *mcevent = dynamic_cast<AliMCEvent*>(event);\r
 774 \r
 775   fHistEventStats->Fill(1,gRefMultiplicity); //all events\r
 776 \r
 777   // check first event in chunk (is not needed for new reconstructions)\r
 778   if(fCheckFirstEventInChunk){\r
 779     AliAnalysisUtils ut;\r
 780     if(ut.IsFirstEventInChunk(event)) \r
 781       return -1.;\r
 782     fHistEventStats->Fill(6,gRefMultiplicity); \r
 783   }\r
 784 \r
 785   // check for pile-up event\r
 786   if(fCheckPileUp){\r
 787     AliAnalysisUtils ut;\r
 788     ut.SetUseMVPlpSelection(kTRUE);\r
 789     ut.SetUseOutOfBunchPileUp(kTRUE);\r
 790     if(ut.IsPileUpEvent(event))\r
 791       return -1.;\r
 792     fHistEventStats->Fill(7,gRefMultiplicity); \r
 793   }\r
 794 \r
 795   // Event trigger bits\r
 796   fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());\r
 797   if(fUseOfflineTrigger)\r
 798     isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
 799   \r
 800   if(isSelectedMain) {\r
 801     fHistEventStats->Fill(2,gRefMultiplicity); //triggered events\r
 802  \r
 803     // Event Vertex MC\r
 804     if(gAnalysisLevel == "MC") {\r
 805       if(!event) {\r
 806         AliError("mcEvent not available");\r
 807         return 0x0;\r
 808       }\r
 809       \r
 810       if(mcevent){\r
 811         AliGenEventHeader *header = dynamic_cast<AliGenEventHeader*>(mcevent->GenEventHeader());\r
 812         if(header){  \r
 813           TArrayF gVertexArray;\r
 814           header->PrimaryVertex(gVertexArray);\r
 815           //Printf("Vertex: %lf (x) - %lf (y) - %lf (z)",\r
 816           //gVertexArray.At(0),\r
 817           //gVertexArray.At(1),\r
 818           //gVertexArray.At(2));\r
 819           fHistEventStats->Fill(3,gRefMultiplicity); //events with a proper vertex\r
 820           if(TMath::Abs(gVertexArray.At(0)) < fVxMax) {\r
 821             if(TMath::Abs(gVertexArray.At(1)) < fVyMax) {\r
 822               if(TMath::Abs(gVertexArray.At(2)) < fVzMax) {\r
 823                 fHistEventStats->Fill(4,gRefMultiplicity);//analyzed events\r
 824 \r
 825                 // get the reference multiplicty or centrality\r
 826                 gRefMultiplicity = GetRefMultiOrCentrality(event);\r
 827 \r
 828                 fHistVx->Fill(gVertexArray.At(0));\r
 829                 fHistVy->Fill(gVertexArray.At(1));\r
 830                 fHistVz->Fill(gVertexArray.At(2),gRefMultiplicity);\r
 831                 \r
 832                 // take only events inside centrality class\r
 833                 if(fUseCentrality) {\r
 834                   if((fImpactParameterMin < gRefMultiplicity) && (fImpactParameterMax > gRefMultiplicity)){\r
 835                     fHistEventStats->Fill(5,gRefMultiplicity); //events with correct centrality\r
 836                     return gRefMultiplicity;        \r
 837                   }//centrality class\r
 838                 }\r
 839                 // take events only within the same multiplicity class\r
 840                 else if(fUseMultiplicity) {\r
 841                   if((gRefMultiplicity > fNumberOfAcceptedTracksMin) && (gRefMultiplicity < fNumberOfAcceptedTracksMax)) {\r
 842                     fHistEventStats->Fill(5,gRefMultiplicity); //events with correct multiplicity\r
 843                     return gRefMultiplicity;\r
 844                   }\r
 845                 }//multiplicity range\r
 846               }//Vz cut\r
 847             }//Vy cut\r
 848           }//Vx cut\r
 849         }//header    \r
 850       }//MC event object\r
 851     }//MC\r
 852     \r
 853     // Event Vertex AOD, ESD, ESDMC\r
 854     else{\r
 855       const AliVVertex *vertex = event->GetPrimaryVertex();\r
 856       \r
 857       if(vertex) {\r
 858         Double32_t fCov[6];\r
 859         vertex->GetCovarianceMatrix(fCov);\r
 860         if(vertex->GetNContributors() > 0) {\r
 861           if(fCov[5] != 0) {\r
 862             fHistEventStats->Fill(3,gRefMultiplicity); //proper vertex\r
 863             if(TMath::Abs(vertex->GetX()) < fVxMax) {\r
 864               if(TMath::Abs(vertex->GetY()) < fVyMax) {\r
 865                 if(TMath::Abs(vertex->GetZ()) < fVzMax) {\r
 866                   fHistEventStats->Fill(4,gRefMultiplicity);//analyzed events\r
 867 \r
 868                   // get the reference multiplicty or centrality\r
 869                   gRefMultiplicity = GetRefMultiOrCentrality(event);\r
 870                   \r
 871                   fHistVx->Fill(vertex->GetX());\r
 872                   fHistVy->Fill(vertex->GetY());\r
 873                   fHistVz->Fill(vertex->GetZ(),gRefMultiplicity);\r
 874                   \r
 875                   // take only events inside centrality class\r
 876                   if(fUseCentrality) {\r
 877                     if((gRefMultiplicity > fCentralityPercentileMin) && (gRefMultiplicity < fCentralityPercentileMax)){\r
 878                       fHistEventStats->Fill(5,gRefMultiplicity); //events with correct centrality\r
 879                       return gRefMultiplicity;          \r
 880                     }//centrality class\r
 881                   }\r
 882                   // take events only within the same multiplicity class\r
 883                   else if(fUseMultiplicity) {\r
 884                     //if(fDebugLevel) \r
 885                     //Printf("N(min): %.0f, N(max): %.0f - N(ref): %.0f",fNumberOfAcceptedTracksMin,\r
 886                     //fNumberOfAcceptedTracksMax,gRefMultiplicity);\r
 887 \r
 888                     if((gRefMultiplicity > fNumberOfAcceptedTracksMin) && (gRefMultiplicity < fNumberOfAcceptedTracksMax)) {\r
 889                       fHistEventStats->Fill(5,gRefMultiplicity); //events with correct multiplicity\r
 890                       return gRefMultiplicity;\r
 891                     }\r
 892                   }//multiplicity range\r
 893                 }//Vz cut\r
 894               }//Vy cut\r
 895             }//Vx cut\r
 896           }//proper vertex resolution\r
 897         }//proper number of contributors\r
 898       }//vertex object valid\r
 899     }//triggered event \r
 900   }//AOD,ESD,ESDMC\r
 901   \r
 902   // in all other cases return -1 (event not accepted)\r
 903   return -1;\r
 904 }\r
 905 \r
 906 \r
 907 //________________________________________________________________________\r
 908 Double_t AliAnalysisTaskBFPsi::GetRefMultiOrCentrality(AliVEvent *event){\r
 909     // Checks the Event cuts\r
 910     // Fills Event statistics histograms\r
 911   \r
 912   Float_t gCentrality = -1.;\r
 913   Double_t gMultiplicity = -1.;\r
 914   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
 915 \r
 916 \r
 917   // calculate centrality always (not only in centrality mode)\r
 918   if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD" || gAnalysisLevel == "MCAODrec" ) { //centrality in AOD header  //++++++++++++++\r
 919     AliAODHeader *header = (AliAODHeader*) event->GetHeader();\r
 920     if(header){\r
 921       gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
 922 \r
 923       // QA for centrality estimators\r
 924       fHistCentStats->Fill(0.,header->GetCentralityP()->GetCentralityPercentile("V0M"));\r
 925       fHistCentStats->Fill(1.,header->GetCentralityP()->GetCentralityPercentile("V0A"));\r
 926       fHistCentStats->Fill(2.,header->GetCentralityP()->GetCentralityPercentile("V0C"));\r
 927       fHistCentStats->Fill(3.,header->GetCentralityP()->GetCentralityPercentile("FMD"));\r
 928       fHistCentStats->Fill(4.,header->GetCentralityP()->GetCentralityPercentile("TRK"));\r
 929       fHistCentStats->Fill(5.,header->GetCentralityP()->GetCentralityPercentile("TKL")); \r
 930       fHistCentStats->Fill(6.,header->GetCentralityP()->GetCentralityPercentile("CL0"));\r
 931       fHistCentStats->Fill(7.,header->GetCentralityP()->GetCentralityPercentile("CL1"));\r
 932       fHistCentStats->Fill(8.,header->GetCentralityP()->GetCentralityPercentile("ZNA"));\r
 933       fHistCentStats->Fill(9.,header->GetCentralityP()->GetCentralityPercentile("ZPA"));\r
 934       fHistCentStats->Fill(10.,header->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));\r
 935       fHistCentStats->Fill(11.,header->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));\r
 936       fHistCentStats->Fill(12.,header->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));\r
 937       \r
 938       // Centrality estimator USED   ++++++++++++++++++++++++++++++\r
 939       fHistCentStatsUsed->Fill(0.,header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data()));\r
 940       \r
 941       // centrality QA (V0M)\r
 942       fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());\r
 943       \r
 944       // centrality QA (reference tracks)\r
 945       fHistRefTracks->Fill(0.,header->GetRefMultiplicity());\r
 946       fHistRefTracks->Fill(1.,header->GetRefMultiplicityPos());\r
 947       fHistRefTracks->Fill(2.,header->GetRefMultiplicityNeg());\r
 948       fHistRefTracks->Fill(3.,header->GetTPConlyRefMultiplicity());\r
 949       fHistRefTracks->Fill(4.,header->GetNumberOfITSClusters(0));\r
 950       fHistRefTracks->Fill(5.,header->GetNumberOfITSClusters(1));\r
 951       fHistRefTracks->Fill(6.,header->GetNumberOfITSClusters(2));\r
 952       fHistRefTracks->Fill(7.,header->GetNumberOfITSClusters(3));\r
 953       fHistRefTracks->Fill(8.,header->GetNumberOfITSClusters(4));\r
 954 \r
 955     }//AOD header\r
 956   }//AOD\r
 957   \r
 958   else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){ // centrality class for ESDs or MC-ESDs\r
 959     AliCentrality *centrality = event->GetCentrality();\r
 960     gCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
 961 \r
 962     // QA for centrality estimators\r
 963     fHistCentStats->Fill(0.,centrality->GetCentralityPercentile("V0M"));\r
 964     fHistCentStats->Fill(1.,centrality->GetCentralityPercentile("V0A"));\r
 965     fHistCentStats->Fill(2.,centrality->GetCentralityPercentile("V0C"));\r
 966     fHistCentStats->Fill(3.,centrality->GetCentralityPercentile("FMD"));\r
 967     fHistCentStats->Fill(4.,centrality->GetCentralityPercentile("TRK"));\r
 968     fHistCentStats->Fill(5.,centrality->GetCentralityPercentile("TKL"));\r
 969     fHistCentStats->Fill(6.,centrality->GetCentralityPercentile("CL0"));\r
 970     fHistCentStats->Fill(7.,centrality->GetCentralityPercentile("CL1"));\r
 971     fHistCentStats->Fill(8.,centrality->GetCentralityPercentile("ZNA"));\r
 972     fHistCentStats->Fill(9.,centrality->GetCentralityPercentile("ZPA"));\r
 973     fHistCentStats->Fill(10.,centrality->GetCentralityPercentile("V0MvsFMD"));\r
 974     fHistCentStats->Fill(11.,centrality->GetCentralityPercentile("TKLvsV0M"));\r
 975     fHistCentStats->Fill(12.,centrality->GetCentralityPercentile("ZEMvsZDC"));\r
 976     \r
 977     // Centrality estimator USED   ++++++++++++++++++++++++++++++\r
 978     fHistCentStatsUsed->Fill(0.,centrality->GetCentralityPercentile(fCentralityEstimator.Data()));\r
 979     \r
 980     // centrality QA (V0M)\r
 981     fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());\r
 982   }//ESD\r
 983 \r
 984   else if(gAnalysisLevel == "MC"){\r
 985     Double_t gImpactParameter = 0.;\r
 986     AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);\r
 987     if(gMCEvent){\r
 988       AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader());      \r
 989       if(headerH){\r
 990         gImpactParameter = headerH->ImpactParameter();\r
 991         gCentrality      = gImpactParameter;\r
 992       }//MC header\r
 993     }//MC event cast\r
 994   }//MC\r
 995 \r
 996   else{\r
 997     gCentrality = -1.;\r
 998   }\r
 999   \r
1000   // calculate reference multiplicity always (not only in multiplicity mode)\r
1001   if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){\r
1002     AliESDEvent* gESDEvent = dynamic_cast<AliESDEvent*>(event);\r
1003     if(gESDEvent){\r
1004       gMultiplicity = fESDtrackCuts->GetReferenceMultiplicity(gESDEvent, AliESDtrackCuts::kTrackletsITSTPC,0.5);\r
1005       fHistMultiplicity->Fill(gMultiplicity);\r
1006     }//AliESDevent cast\r
1007   }//ESD mode\r
1008 \r
1009   else if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD" || gAnalysisLevel == "MCAODrec" ){\r
1010     AliAODHeader *header = (AliAODHeader*) event->GetHeader();\r
1011     if ((fMultiplicityEstimator == "V0M")||\r
1012         (fMultiplicityEstimator == "V0A")||\r
1013         (fMultiplicityEstimator == "V0C") ||\r
1014         (fMultiplicityEstimator == "TPC")) {\r
1015       gMultiplicity = GetReferenceMultiplicityFromAOD(event);\r
1016       if(fDebugLevel) Printf("Reference multiplicity (calculated): %.0f",gMultiplicity);\r
1017     }\r
1018     else {\r
1019       if(header)\r
1020         gMultiplicity = header->GetRefMultiplicity();\r
1021       if(fDebugLevel) Printf("Reference multiplicity (AOD header): %.0f",gMultiplicity);\r
1022     }\r
1023     fHistMultiplicity->Fill(gMultiplicity);\r
1024   }//AOD mode\r
1025   else if(gAnalysisLevel == "MC") {\r
1026     AliMCEvent* gMCEvent = dynamic_cast<AliMCEvent*>(event);\r
1027     //Calculating the multiplicity as the number of charged primaries\r
1028     //within \pm 0.8 in eta and pT > 0.1 GeV/c\r
1029     for(Int_t iParticle = 0; iParticle < gMCEvent->GetNumberOfPrimaries(); iParticle++) {\r
1030       AliMCParticle* track = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iParticle));\r
1031       if (!track) {\r
1032         AliError(Form("Could not receive particle %d", iParticle));\r
1033         continue;\r
1034       }\r
1035       \r
1036       //exclude non stable particles\r
1037       if(!(gMCEvent->IsPhysicalPrimary(iParticle))) continue;\r
1038       \r
1039       //++++++++++++++++\r
1040       if (fMultiplicityEstimator == "V0M") {\r
1041         if((track->Eta() > 5.1 || track->Eta() < 2.8)&&(track->Eta() < -3.7 || track->Eta() > -1.7)) \r
1042           continue;}\r
1043       else if (fMultiplicityEstimator == "V0A") {\r
1044         if(track->Eta() > 5.1 || track->Eta() < 2.8)  continue;}\r
1045       else if (fMultiplicityEstimator == "V0C") {\r
1046         if(track->Eta() > -1.7 || track->Eta() < -3.7)  continue;}\r
1047       else if (fMultiplicityEstimator == "TPC") {\r
1048         if(track->Eta() < fEtaMin || track->Eta() > fEtaMax)  continue;\r
1049         if(track->Pt() < fPtMin || track->Pt() > fPtMax)  continue;\r
1050       }\r
1051       else{\r
1052         if(track->Pt() < fPtMin || track->Pt() > fPtMax)  continue;\r
1053         if(track->Eta() < fEtaMin || track->Eta() > fEtaMax)  continue;\r
1054       }\r
1055       //++++++++++++++++\r
1056       \r
1057       if(track->Charge() == 0) continue;\r
1058       \r
1059       gMultiplicity += 1;\r
1060     }//loop over primaries\r
1061     fHistMultiplicity->Fill(gMultiplicity);\r
1062   }//MC mode\r
1063   else{\r
1064     gMultiplicity = -1;\r
1065   }\r
1066   \r
1067 \r
1068   // decide what should be returned only here\r
1069   Double_t lReturnVal = -100;\r
1070   if(fEventClass=="Multiplicity"){\r
1071     lReturnVal = gMultiplicity;\r
1072   }else if(fEventClass=="Centrality"){\r
1073     lReturnVal = gCentrality;\r
1074   }\r
1075   return lReturnVal;\r
1076 }\r
1077 \r
1078 //________________________________________________________________________\r
1079 Double_t AliAnalysisTaskBFPsi::GetReferenceMultiplicityFromAOD(AliVEvent *event){\r
1080   //Function that returns the reference multiplicity from AODs (data or reco MC)\r
1081   //Different ref. mult. implemented: V0M, V0A, V0C, TPC\r
1082   Double_t gRefMultiplicity = 0., gRefMultiplicityTPC = 0., gRefMultiplicityVZERO = 0.;\r
1083 \r
1084   // Loop over tracks in event\r
1085   for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r
1086     AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));\r
1087     if (!aodTrack) {\r
1088       AliError(Form("Could not receive track %d", iTracks));\r
1089       continue;\r
1090     }\r
1091     \r
1092     // AOD track cuts\r
1093     if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;\r
1094             \r
1095     if(aodTrack->Charge() == 0) continue;\r
1096     // Kinematics cuts from ESD track cuts\r
1097     if( aodTrack->Pt() < fPtMin || aodTrack->Pt() > fPtMax)      continue;\r
1098     if( aodTrack->Eta() < fEtaMin || aodTrack->Eta() > fEtaMax)  continue;\r
1099     \r
1100     //=================PID (so far only for electron rejection)==========================//\r
1101     if(fElectronRejection) {\r
1102       // get the electron nsigma\r
1103       Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));\r
1104       \r
1105       // check only for given momentum range\r
1106       if( aodTrack->Pt() > fElectronRejectionMinPt && aodTrack->Pt() < fElectronRejectionMaxPt ){\r
1107         //look only at electron nsigma\r
1108         if(!fElectronOnlyRejection) {\r
1109           //Make the decision based on the n-sigma of electrons\r
1110           if(nSigma < fElectronRejectionNSigma) continue;\r
1111         }\r
1112         else {\r
1113           Double_t nSigmaPions   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));\r
1114           Double_t nSigmaKaons   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));\r
1115           Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));\r
1116           \r
1117           //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)\r
1118           if(nSigma < fElectronRejectionNSigma\r
1119              && nSigmaPions   > fElectronRejectionNSigma\r
1120              && nSigmaKaons   > fElectronRejectionNSigma\r
1121              && nSigmaProtons > fElectronRejectionNSigma ) continue;\r
1122         }\r
1123       }\r
1124     }//electron rejection\r
1125     \r
1126     gRefMultiplicityTPC += 1.0;\r
1127   }// track loop\r
1128   \r
1129   //VZERO segmentation in two detectors (0-31: VZERO-C, 32-63: VZERO-A)\r
1130   for(Int_t i = 0; i < 64; i++) {\r
1131     fHistVZEROSignal->Fill(i,event->GetVZEROEqMultiplicity(i));\r
1132     \r
1133     if(fMultiplicityEstimator == "V0C") {\r
1134       if(i > 31) continue;}\r
1135     else if(fMultiplicityEstimator == "V0A") {\r
1136       if(i < 32) continue;}\r
1137     \r
1138     gRefMultiplicityVZERO += event->GetVZEROEqMultiplicity(i);\r
1139   }//loop over PMTs\r
1140   \r
1141   if(fDebugLevel) \r
1142     Printf("VZERO multiplicity: %.0f - TPC multiplicity: %.0f",gRefMultiplicityVZERO,gRefMultiplicityTPC);\r
1143 \r
1144   fHistTPCvsVZEROMultiplicity->Fill(gRefMultiplicityVZERO,gRefMultiplicityTPC);\r
1145 \r
1146   if(fMultiplicityEstimator == "TPC") \r
1147     gRefMultiplicity = gRefMultiplicityTPC;\r
1148   else if((fMultiplicityEstimator == "V0M")||\r
1149           (fMultiplicityEstimator == "V0A")||\r
1150           (fMultiplicityEstimator == "V0C")) \r
1151     gRefMultiplicity = gRefMultiplicityVZERO;\r
1152   \r
1153   return gRefMultiplicity;\r
1154 }\r
1155 \r
1156 //________________________________________________________________________\r
1157 Double_t AliAnalysisTaskBFPsi::GetEventPlane(AliVEvent *event){\r
1158   // Get the event plane\r
1159 \r
1160   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
1161 \r
1162   Float_t gVZEROEventPlane    = -10.;\r
1163   Float_t gReactionPlane      = -10.;\r
1164   Double_t qxTot = 0.0, qyTot = 0.0;\r
1165 \r
1166   //MC: from reaction plane\r
1167   if(gAnalysisLevel == "MC"){\r
1168     if(!event) {\r
1169       AliError("mcEvent not available");\r
1170       return 0x0;\r
1171     }\r
1172 \r
1173     AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);\r
1174     if(gMCEvent){\r
1175       AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader());    \r
1176       if (headerH) {\r
1177         gReactionPlane = headerH->ReactionPlaneAngle();\r
1178         //gReactionPlane *= TMath::RadToDeg();\r
1179       }//MC header\r
1180     }//MC event cast\r
1181   }//MC\r
1182   \r
1183   // AOD,ESD,ESDMC: from VZERO Event Plane\r
1184   else{\r
1185    \r
1186     AliEventplane *ep = event->GetEventplane();\r
1187     if(ep){ \r
1188       gVZEROEventPlane = ep->CalculateVZEROEventPlane(event,10,2,qxTot,qyTot);\r
1189       if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();\r
1190       //gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();\r
1191       gReactionPlane = gVZEROEventPlane;\r
1192     }\r
1193   }//AOD,ESD,ESDMC\r
1194 \r
1195   return gReactionPlane;\r
1196 }\r
1197 \r
1198 //________________________________________________________________________\r
1199 Double_t AliAnalysisTaskBFPsi::GetTrackbyTrackCorrectionMatrix( Double_t vEta, \r
1200                                                                 Double_t vPhi, \r
1201                                                                 Double_t vPt, \r
1202                                                                 Short_t vCharge, \r
1203                                                                 Double_t gCentrality) {\r
1204   // -- Get efficiency correction of particle dependent on (eta, phi, pt, charge, centrality) \r
1205 \r
1206   Double_t correction = 1.;\r
1207   Int_t binEta = 0, binPt = 0, binPhi = 0;\r
1208 \r
1209   //Printf("EtaMAx: %lf - EtaMin: %lf - EtaBin: %lf", fEtaMaxForCorrections,fEtaMinForCorrections,fEtaBinForCorrections);\r
1210   if(fEtaBinForCorrections != 0) {\r
1211     Double_t widthEta = (fEtaMaxForCorrections - fEtaMinForCorrections)/fEtaBinForCorrections;\r
1212     if(fEtaMaxForCorrections != fEtaMinForCorrections) \r
1213       binEta = (Int_t)((vEta-fEtaMinForCorrections)/widthEta)+1;\r
1214   }\r
1215 \r
1216   if(fPtBinForCorrections != 0) {\r
1217     Double_t widthPt = (fPtMaxForCorrections - fPtMinForCorrections)/fPtBinForCorrections;\r
1218     if(fPtMaxForCorrections != fPtMinForCorrections) \r
1219       binPt = (Int_t)((vPt-fPtMinForCorrections)/widthPt) + 1;\r
1220   }\r
1221  \r
1222   if(fPhiBinForCorrections != 0) {\r
1223     Double_t widthPhi = (fPhiMaxForCorrections - fPhiMinForCorrections)/fPhiBinForCorrections;\r
1224     if(fPhiMaxForCorrections != fPhiMinForCorrections) \r
1225       binPhi = (Int_t)((vPhi-fPhiMinForCorrections)/widthPhi)+ 1;\r
1226   }\r
1227 \r
1228   Int_t gCentralityInt = -1;\r
1229   for (Int_t i=0; i<fCentralityArrayBinsForCorrections-1; i++){\r
1230     if((fCentralityArrayForCorrections[i] <= gCentrality)&&(gCentrality <= fCentralityArrayForCorrections[i+1])){\r
1231       gCentralityInt = i;\r
1232       break;\r
1233     }\r
1234   }  \r
1235 \r
1236   // centrality not in array --> no correction\r
1237   if(gCentralityInt < 0){\r
1238     correction = 1.;\r
1239   }\r
1240   else{\r
1241     \r
1242     //Printf("//=============CENTRALITY=============// %d:",gCentralityInt);\r
1243     \r
1244     if(fHistCorrectionPlus[gCentralityInt]){\r
1245       if (vCharge > 0) {\r
1246         correction = fHistCorrectionPlus[gCentralityInt]->GetBinContent(fHistCorrectionPlus[gCentralityInt]->GetBin(binEta, binPt, binPhi));\r
1247         //Printf("CORRECTIONplus: %.2f | Centrality %d",correction,gCentralityInt);  \r
1248       }\r
1249       if (vCharge < 0) {\r
1250         correction = fHistCorrectionMinus[gCentralityInt]->GetBinContent(fHistCorrectionMinus[gCentralityInt]->GetBin(binEta, binPt, binPhi));\r
1251         //Printf("CORRECTIONminus: %.2f | Centrality %d",correction,gCentralityInt);\r
1252       }\r
1253     }\r
1254     else {\r
1255       correction = 1.;\r
1256     }\r
1257   }//centrality in array\r
1258   \r
1259   if (correction == 0.) { \r
1260     AliError(Form("Should not happen : bin content = 0. >> eta: %.2f | phi : %.2f | pt : %.2f | cent %d",vEta, vPhi, vPt, gCentralityInt)); \r
1261     correction = 1.; \r
1262   } \r
1263   \r
1264   return correction;\r
1265 }\r
1266 \r
1267 //________________________________________________________________________\r
1268 TObjArray* AliAnalysisTaskBFPsi::GetAcceptedTracks(AliVEvent *event, Double_t gCentrality, Double_t gReactionPlane){\r
1269   // Returns TObjArray with tracks after all track cuts (only for AOD!)\r
1270   // Fills QA histograms\r
1271 \r
1272   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
1273 \r
1274   //output TObjArray holding all good tracks\r
1275   TObjArray* tracksAccepted = new TObjArray;\r
1276   tracksAccepted->SetOwner(kTRUE);\r
1277 \r
1278   Short_t vCharge;\r
1279   Double_t vEta;\r
1280   Double_t vY;\r
1281   Double_t vPhi;\r
1282   Double_t vPt;\r
1283 \r
1284 \r
1285   if(gAnalysisLevel == "AOD") { // handling of TPC only tracks different in AOD and ESD\r
1286     // Loop over tracks in event\r
1287     for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r
1288       AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));\r
1289       if (!aodTrack) {\r
1290         AliError(Form("Could not receive track %d", iTracks));\r
1291         continue;\r
1292       }\r
1293       \r
1294       // AOD track cuts\r
1295       \r
1296       // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C\r
1297       // take only TPC only tracks \r
1298       //fHistTrackStats->Fill(aodTrack->GetFilterMap());\r
1299       for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){\r
1300         fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));\r
1301       }\r
1302       if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;\r
1303       \r
1304       vCharge = aodTrack->Charge();\r
1305       vEta    = aodTrack->Eta();\r
1306       vY      = aodTrack->Y();\r
1307       vPhi    = aodTrack->Phi();// * TMath::RadToDeg();\r
1308       vPt     = aodTrack->Pt();\r
1309       \r
1310       //===========================PID (so far only for electron rejection)===============================//                \r
1311       if(fElectronRejection) {\r
1312 \r
1313         // get the electron nsigma\r
1314         Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));\r
1315         \r
1316         //Fill QA before the PID\r
1317         fHistdEdxVsPTPCbeforePID -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());\r
1318         fHistNSigmaTPCvsPtbeforePID -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); \r
1319         //end of QA-before pid\r
1320         \r
1321         // check only for given momentum range\r
1322         if( vPt > fElectronRejectionMinPt && vPt < fElectronRejectionMaxPt ){\r
1323                   \r
1324           //look only at electron nsigma\r
1325           if(!fElectronOnlyRejection){\r
1326             \r
1327             //Make the decision based on the n-sigma of electrons\r
1328             if(nSigma < fElectronRejectionNSigma) continue;\r
1329           }\r
1330           else{\r
1331             \r
1332             Double_t nSigmaPions   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));\r
1333             Double_t nSigmaKaons   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));\r
1334             Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));\r
1335             \r
1336             //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)\r
1337             if(nSigma < fElectronRejectionNSigma\r
1338                && nSigmaPions   > fElectronRejectionNSigma\r
1339                && nSigmaKaons   > fElectronRejectionNSigma\r
1340                && nSigmaProtons > fElectronRejectionNSigma ) continue;\r
1341           }\r
1342         }\r
1343   \r
1344         //Fill QA after the PID\r
1345         fHistdEdxVsPTPCafterPID -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());\r
1346         fHistNSigmaTPCvsPtafterPID -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); \r
1347         \r
1348       }\r
1349       //===========================end of PID (so far only for electron rejection)===============================//\r
1350 \r
1351       Float_t dcaXY = aodTrack->DCA();      // this is the DCA from global track (not exactly what is cut on)\r
1352       Float_t dcaZ  = aodTrack->ZAtDCA();   // this is the DCA from global track (not exactly what is cut on)\r
1353       \r
1354       \r
1355       // Kinematics cuts from ESD track cuts\r
1356       if( vPt < fPtMin || vPt > fPtMax)      continue;\r
1357       if( vEta < fEtaMin || vEta > fEtaMax)  continue;\r
1358       \r
1359       // Extra DCA cuts (for systematic studies [!= -1])\r
1360       if( fDCAxyCut != -1 && fDCAzCut != -1){\r
1361         if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){\r
1362           continue;  // 2D cut\r
1363         }\r
1364       }\r
1365       \r
1366       // Extra TPC cuts (for systematic studies [!= -1])\r
1367       if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){\r
1368         continue;\r
1369       }\r
1370       if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){\r
1371         continue;\r
1372       }\r
1373       \r
1374       // fill QA histograms\r
1375       fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());\r
1376       fHistDCA->Fill(dcaZ,dcaXY);\r
1377       fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality);\r
1378       fHistPt->Fill(vPt,gCentrality);\r
1379       fHistEta->Fill(vEta,gCentrality);\r
1380       fHistRapidity->Fill(vY,gCentrality);\r
1381       if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
1382       else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
1383       fHistPhi->Fill(vPhi,gCentrality);\r
1384       if(vCharge > 0)      fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);                  \r
1385       else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r
1386       \r
1387       //=======================================correction\r
1388       Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  \r
1389       //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);\r
1390       \r
1391       // add the track to the TObjArray\r
1392       tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));  \r
1393     }//track loop\r
1394   }// AOD analysis\r
1395 \r
1396   //==============================================================================================================\r
1397   else if(gAnalysisLevel == "MCAOD") {\r
1398     \r
1399     AliMCEvent* mcEvent = MCEvent();\r
1400     if (!mcEvent) {\r
1401       AliError("ERROR: Could not retrieve MC event");\r
1402     }\r
1403     else{\r
1404       \r
1405       for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {\r
1406         AliAODMCParticle *aodTrack = (AliAODMCParticle*) mcEvent->GetTrack(iTracks); \r
1407         if (!aodTrack) {\r
1408           AliError(Form("ERROR: Could not receive track %d (mc loop)", iTracks));\r
1409           continue;\r
1410         }\r
1411         \r
1412         if(!aodTrack->IsPhysicalPrimary()) continue;   \r
1413         \r
1414         vCharge = aodTrack->Charge();\r
1415         vEta    = aodTrack->Eta();\r
1416         vY      = aodTrack->Y();\r
1417         vPhi    = aodTrack->Phi();// * TMath::RadToDeg();\r
1418         vPt     = aodTrack->Pt();\r
1419         \r
1420         // Kinematics cuts from ESD track cuts\r
1421         if( vPt < fPtMin || vPt > fPtMax)      continue;\r
1422         if( vEta < fEtaMin || vEta > fEtaMax)  continue;\r
1423         \r
1424         // Remove neutral tracks\r
1425         if( vCharge == 0 ) continue;\r
1426         \r
1427         //Exclude resonances\r
1428         if(fExcludeResonancesInMC) {\r
1429           \r
1430           Bool_t kExcludeParticle = kFALSE;\r
1431           Int_t gMotherIndex = aodTrack->GetMother();\r
1432           if(gMotherIndex != -1) {\r
1433             AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));\r
1434             if(motherTrack) {\r
1435               Int_t pdgCodeOfMother = motherTrack->GetPdgCode();\r
1436               //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {\r
1437               //if(pdgCodeOfMother == 113) {\r
1438               if(pdgCodeOfMother == 113  // rho0\r
1439                  || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+\r
1440                  // || pdgCodeOfMother == 221  // eta\r
1441                  // || pdgCodeOfMother == 331  // eta'\r
1442                  // || pdgCodeOfMother == 223  // omega\r
1443                  // || pdgCodeOfMother == 333  // phi\r
1444                  || pdgCodeOfMother == 311  || pdgCodeOfMother == -311 // K0\r
1445                  // || pdgCodeOfMother == 313  || pdgCodeOfMother == -313 // K0*\r
1446                  // || pdgCodeOfMother == 323  || pdgCodeOfMother == -323 // K+*\r
1447                  || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda\r
1448                  || pdgCodeOfMother == 111  // pi0 Dalitz\r
1449                  || pdgCodeOfMother == 22  // photon\r
1450                  ) {\r
1451                 kExcludeParticle = kTRUE;\r
1452               }\r
1453             }\r
1454           }\r
1455           \r
1456           //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi\r
1457           if(kExcludeParticle) continue;\r
1458         }\r
1459 \r
1460         //Exclude electrons with PDG\r
1461         if(fExcludeElectronsInMC) {\r
1462           \r
1463           if(TMath::Abs(aodTrack->GetPdgCode()) == 11) continue;\r
1464           \r
1465         }\r
1466         \r
1467         // fill QA histograms\r
1468         fHistPt->Fill(vPt,gCentrality);\r
1469         fHistEta->Fill(vEta,gCentrality);\r
1470         fHistRapidity->Fill(vY,gCentrality);\r
1471         if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
1472         else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
1473         fHistPhi->Fill(vPhi,gCentrality);\r
1474         if(vCharge > 0)      fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);                \r
1475         else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r
1476         \r
1477         //=======================================correction\r
1478         Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  \r
1479         //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);   \r
1480         \r
1481         // add the track to the TObjArray\r
1482         tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));  \r
1483       }//aodTracks\r
1484     }//MC event\r
1485   }//MCAOD\r
1486   //==============================================================================================================\r
1487 \r
1488   //==============================================================================================================\r
1489   else if(gAnalysisLevel == "MCAODrec") {\r
1490     \r
1491     /* fAOD = dynamic_cast<AliAODEvent*>(InputEvent());\r
1492     if (!fAOD) {\r
1493       printf("ERROR: fAOD not available\n");\r
1494       return;\r
1495       }*/\r
1496     \r
1497     fArrayMC = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName())); \r
1498     if (!fArrayMC) { \r
1499        AliError("No array of MC particles found !!!");\r
1500     }\r
1501 \r
1502     AliMCEvent* mcEvent = MCEvent();\r
1503     if (!mcEvent) {\r
1504        AliError("ERROR: Could not retrieve MC event");\r
1505        return tracksAccepted;\r
1506     }\r
1507      \r
1508     for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r
1509       AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));\r
1510       if (!aodTrack) {\r
1511         AliError(Form("Could not receive track %d", iTracks));\r
1512         continue;\r
1513       }\r
1514       \r
1515       for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){\r
1516         fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));\r
1517       }\r
1518       if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;\r
1519       \r
1520       vCharge = aodTrack->Charge();\r
1521       vEta    = aodTrack->Eta();\r
1522       vY      = aodTrack->Y();\r
1523       vPhi    = aodTrack->Phi();// * TMath::RadToDeg();\r
1524       vPt     = aodTrack->Pt();\r
1525       \r
1526       //===========================use MC information for Kinematics===============================//               \r
1527       if(fUseMCforKinematics){\r
1528 \r
1529         Int_t label = TMath::Abs(aodTrack->GetLabel());\r
1530         AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);\r
1531 \r
1532         if(AODmcTrack){\r
1533           vCharge = AODmcTrack->Charge();\r
1534           vEta    = AODmcTrack->Eta();\r
1535           vY      = AODmcTrack->Y();\r
1536           vPhi    = AODmcTrack->Phi();// * TMath::RadToDeg();\r
1537           vPt     = AODmcTrack->Pt();\r
1538         }\r
1539         else{\r
1540           AliDebug(1, "no MC particle for this track"); \r
1541           continue;\r
1542         }\r
1543       }\r
1544       //===========================end of use MC information for Kinematics========================//               \r
1545 \r
1546 \r
1547       //===========================PID (so far only for electron rejection)===============================//                \r
1548       if(fElectronRejection) {\r
1549 \r
1550         // get the electron nsigma\r
1551         Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));\r
1552 \r
1553         //Fill QA before the PID\r
1554         fHistdEdxVsPTPCbeforePID -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());\r
1555         fHistNSigmaTPCvsPtbeforePID -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); \r
1556         //end of QA-before pid\r
1557         \r
1558         // check only for given momentum range\r
1559         if( vPt > fElectronRejectionMinPt && vPt < fElectronRejectionMaxPt ){\r
1560                   \r
1561           //look only at electron nsigma\r
1562           if(!fElectronOnlyRejection){\r
1563             \r
1564             //Make the decision based on the n-sigma of electrons\r
1565             if(nSigma < fElectronRejectionNSigma) continue;\r
1566           }\r
1567           else{\r
1568             \r
1569             Double_t nSigmaPions   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));\r
1570             Double_t nSigmaKaons   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));\r
1571             Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));\r
1572             \r
1573             //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)\r
1574             if(nSigma < fElectronRejectionNSigma\r
1575                && nSigmaPions   > fElectronRejectionNSigma\r
1576                && nSigmaKaons   > fElectronRejectionNSigma\r
1577                && nSigmaProtons > fElectronRejectionNSigma ) continue;\r
1578           }\r
1579         }\r
1580   \r
1581         //Fill QA after the PID\r
1582         fHistdEdxVsPTPCafterPID -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());\r
1583         fHistNSigmaTPCvsPtafterPID -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); \r
1584         \r
1585       }\r
1586       //===========================end of PID (so far only for electron rejection)===============================//\r
1587       \r
1588       Float_t dcaXY = aodTrack->DCA();      // this is the DCA from global track (not exactly what is cut on)\r
1589       Float_t dcaZ  = aodTrack->ZAtDCA();   // this is the DCA from global track (not exactly what is cut on)\r
1590             \r
1591       // Kinematics cuts from ESD track cuts\r
1592       if( vPt < fPtMin || vPt > fPtMax)      continue;\r
1593       if( vEta < fEtaMin || vEta > fEtaMax)  continue;\r
1594       \r
1595       // Extra DCA cuts (for systematic studies [!= -1])\r
1596       if( fDCAxyCut != -1 && fDCAzCut != -1){\r
1597         if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){\r
1598           continue;  // 2D cut\r
1599         }\r
1600       }\r
1601       \r
1602       // Extra TPC cuts (for systematic studies [!= -1])\r
1603       if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){\r
1604         continue;\r
1605       }\r
1606       if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){\r
1607         continue;\r
1608       }\r
1609 \r
1610       //Exclude resonances\r
1611       if(fExcludeResonancesInMC) {\r
1612         \r
1613         Bool_t kExcludeParticle = kFALSE;\r
1614 \r
1615         Int_t label = TMath::Abs(aodTrack->GetLabel());\r
1616         AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);\r
1617       \r
1618         if (AODmcTrack){ \r
1619           //if (AODmcTrack->IsPhysicalPrimary()){\r
1620           \r
1621           Int_t gMotherIndex = AODmcTrack->GetMother();\r
1622           if(gMotherIndex != -1) {\r
1623             AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));\r
1624             if(motherTrack) {\r
1625               Int_t pdgCodeOfMother = motherTrack->GetPdgCode();\r
1626               if(pdgCodeOfMother == 113  // rho0\r
1627                  || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+\r
1628                  // || pdgCodeOfMother == 221  // eta\r
1629                  // || pdgCodeOfMother == 331  // eta'\r
1630                  // || pdgCodeOfMother == 223  // omega\r
1631                  // || pdgCodeOfMother == 333  // phi\r
1632                  || pdgCodeOfMother == 311  || pdgCodeOfMother == -311 // K0\r
1633                  // || pdgCodeOfMother == 313  || pdgCodeOfMother == -313 // K0*\r
1634                  // || pdgCodeOfMother == 323  || pdgCodeOfMother == -323 // K+*\r
1635                  || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda\r
1636                  || pdgCodeOfMother == 111  // pi0 Dalitz\r
1637                  || pdgCodeOfMother == 22   // photon\r
1638                  ) {\r
1639                 kExcludeParticle = kTRUE;\r
1640               }\r
1641             }\r
1642           }\r
1643         }       \r
1644         //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi\r
1645         if(kExcludeParticle) continue;\r
1646       }\r
1647 \r
1648       //Exclude electrons with PDG\r
1649       if(fExcludeElectronsInMC) {\r
1650         \r
1651         Int_t label = TMath::Abs(aodTrack->GetLabel());\r
1652         AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);\r
1653         \r
1654         if (AODmcTrack){ \r
1655           if(TMath::Abs(AODmcTrack->GetPdgCode()) == 11) continue;\r
1656         }\r
1657       }\r
1658       \r
1659       // fill QA histograms\r
1660       fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());\r
1661       fHistDCA->Fill(dcaZ,dcaXY);\r
1662       fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality);\r
1663       fHistPt->Fill(vPt,gCentrality);\r
1664       fHistEta->Fill(vEta,gCentrality);\r
1665       fHistRapidity->Fill(vY,gCentrality);\r
1666       if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
1667       else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
1668       fHistPhi->Fill(vPhi,gCentrality);\r
1669       if(vCharge > 0)      fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);                  \r
1670       else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r
1671       \r
1672       //=======================================correction\r
1673       Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  \r
1674       //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);\r
1675       \r
1676       // add the track to the TObjArray\r
1677       tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));  \r
1678     }//track loop\r
1679   }//MCAODrec\r
1680   //==============================================================================================================\r
1681 \r
1682   else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD") { // handling of TPC only tracks different in AOD and ESD\r
1683 \r
1684     AliESDtrack *trackTPC   = NULL;\r
1685     AliMCParticle *mcTrack   = NULL;\r
1686 \r
1687     AliMCEvent*  mcEvent     = NULL;\r
1688     \r
1689     // for MC ESDs use also MC information\r
1690     if(gAnalysisLevel == "MCESD"){\r
1691       mcEvent = MCEvent(); \r
1692       if (!mcEvent) {\r
1693         AliError("mcEvent not available");\r
1694         return tracksAccepted;\r
1695       }\r
1696     }\r
1697     \r
1698     // Loop over tracks in event\r
1699     for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r
1700       AliESDtrack* track = dynamic_cast<AliESDtrack *>(event->GetTrack(iTracks));\r
1701       if (!track) {\r
1702         AliError(Form("Could not receive track %d", iTracks));\r
1703         continue;\r
1704       } \r
1705 \r
1706       // for MC ESDs use also MC information --> MC track not used further???\r
1707       if(gAnalysisLevel == "MCESD"){\r
1708         Int_t label = TMath::Abs(track->GetLabel());\r
1709         if(label > mcEvent->GetNumberOfTracks()) continue;\r
1710         if(label > mcEvent->GetNumberOfPrimaries()) continue;\r
1711         \r
1712         mcTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label));\r
1713         if(!mcTrack) continue;\r
1714       }\r
1715 \r
1716       // take only TPC only tracks\r
1717       trackTPC   = new AliESDtrack();\r
1718       if(!track->FillTPCOnlyTrack(*trackTPC)) continue;\r
1719       \r
1720       //ESD track cuts\r
1721       if(fESDtrackCuts) \r
1722         if(!fESDtrackCuts->AcceptTrack(trackTPC)) continue;\r
1723       \r
1724       // fill QA histograms\r
1725       Float_t b[2];\r
1726       Float_t bCov[3];\r
1727       trackTPC->GetImpactParameters(b,bCov);\r
1728       if (bCov[0]<=0 || bCov[2]<=0) {\r
1729         AliDebug(1, "Estimated b resolution lower or equal zero!");\r
1730         bCov[0]=0; bCov[2]=0;\r
1731       }\r
1732       \r
1733       Int_t nClustersTPC = -1;\r
1734       nClustersTPC = trackTPC->GetTPCNclsIter1();   // TPC standalone\r
1735       //nClustersTPC = track->GetTPCclusters(0);   // global track\r
1736       Float_t chi2PerClusterTPC = -1;\r
1737       if (nClustersTPC!=0) {\r
1738         chi2PerClusterTPC = trackTPC->GetTPCchi2Iter1()/Float_t(nClustersTPC);      // TPC standalone\r
1739         //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);     // global track\r
1740       }\r
1741       \r
1742       //===========================PID===============================//             \r
1743       if(fUsePID) {\r
1744         Double_t prob[AliPID::kSPECIES]={0.};\r
1745         Double_t probTPC[AliPID::kSPECIES]={0.};\r
1746         Double_t probTOF[AliPID::kSPECIES]={0.};\r
1747         Double_t probTPCTOF[AliPID::kSPECIES]={0.};\r
1748         \r
1749         Double_t nSigma = 0.;\r
1750         UInt_t detUsedTPC = 0;\r
1751         UInt_t detUsedTOF = 0;\r
1752         UInt_t detUsedTPCTOF = 0;\r
1753         \r
1754         //Decide what detector configuration we want to use\r
1755         switch(fPidDetectorConfig) {\r
1756         case kTPCpid:\r
1757           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);\r
1758           nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest));\r
1759           detUsedTPC = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPC);\r
1760           for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
1761             prob[iSpecies] = probTPC[iSpecies];\r
1762           break;\r
1763         case kTOFpid:\r
1764           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);\r
1765           nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest));\r
1766           detUsedTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTOF);\r
1767           for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
1768             prob[iSpecies] = probTOF[iSpecies];\r
1769           break;\r
1770         case kTPCTOF:\r
1771           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);\r
1772           detUsedTPCTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPCTOF);\r
1773           for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
1774             prob[iSpecies] = probTPCTOF[iSpecies];\r
1775           break;\r
1776         default:\r
1777           break;\r
1778         }//end switch: define detector mask\r
1779         \r
1780         //Filling the PID QA\r
1781         Double_t tofTime = -999., length = 999., tof = -999.;\r
1782         Double_t c = TMath::C()*1.E-9;// m/ns\r
1783         Double_t beta = -999.;\r
1784         Double_t  nSigmaTOFForParticleOfInterest = -999.;\r
1785         if ( (track->IsOn(AliESDtrack::kTOFin)) &&\r
1786              (track->IsOn(AliESDtrack::kTIME))  ) { \r
1787           tofTime = track->GetTOFsignal();//in ps\r
1788           length = track->GetIntegratedLength();\r
1789           tof = tofTime*1E-3; // ns     \r
1790           \r
1791           if (tof <= 0) {\r
1792             //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");\r
1793             continue;\r
1794           }\r
1795           if (length <= 0){\r
1796             //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");\r
1797             continue;\r
1798           }\r
1799           \r
1800           length = length*0.01; // in meters\r
1801           tof = tof*c;\r
1802           beta = length/tof;\r
1803           \r
1804           nSigmaTOFForParticleOfInterest = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest);\r
1805           fHistBetavsPTOFbeforePID ->Fill(track->P()*track->Charge(),beta);\r
1806           fHistProbTOFvsPtbeforePID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);\r
1807           fHistNSigmaTOFvsPtbeforePID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);\r
1808         }//TOF signal \r
1809         \r
1810         \r
1811         Double_t  nSigmaTPCForParticleOfInterest = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest);\r
1812         fHistdEdxVsPTPCbeforePID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());\r
1813         fHistProbTPCvsPtbeforePID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); \r
1814         fHistNSigmaTPCvsPtbeforePID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); \r
1815         fHistProbTPCTOFvsPtbeforePID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);\r
1816         //end of QA-before pid\r
1817         \r
1818         if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {\r
1819           //Make the decision based on the n-sigma\r
1820           if(fUsePIDnSigma) {\r
1821             if(nSigma > fPIDNSigma) continue;}\r
1822           \r
1823           //Make the decision based on the bayesian\r
1824           else if(fUsePIDPropabilities) {\r
1825             if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;\r
1826             if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue;      \r
1827           }\r
1828           \r
1829           //Fill QA after the PID\r
1830           fHistBetavsPTOFafterPID ->Fill(track->P()*track->Charge(),beta);\r
1831           fHistProbTOFvsPtafterPID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);\r
1832           fHistNSigmaTOFvsPtafterPID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);\r
1833           \r
1834           fHistdEdxVsPTPCafterPID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());\r
1835           fHistProbTPCvsPtafterPID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); \r
1836           fHistProbTPCTOFvsPtafterPID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);\r
1837           fHistNSigmaTPCvsPtafterPID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); \r
1838         }\r
1839       }\r
1840       //===========================PID===============================//\r
1841       vCharge = trackTPC->Charge();\r
1842       vY      = trackTPC->Y();\r
1843       vEta    = trackTPC->Eta();\r
1844       vPhi    = trackTPC->Phi();// * TMath::RadToDeg();\r
1845       vPt     = trackTPC->Pt();\r
1846       fHistClus->Fill(trackTPC->GetITSclusters(0),nClustersTPC);\r
1847       fHistDCA->Fill(b[1],b[0]);\r
1848       fHistChi2->Fill(chi2PerClusterTPC,gCentrality);\r
1849       fHistPt->Fill(vPt,gCentrality);\r
1850       fHistEta->Fill(vEta,gCentrality);\r
1851       fHistPhi->Fill(vPhi,gCentrality);\r
1852       if(vCharge > 0)      fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);\r
1853       else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r
1854       fHistRapidity->Fill(vY,gCentrality);\r
1855       if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
1856       else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
1857       \r
1858       //=======================================correction\r
1859       Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  \r
1860       //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);\r
1861 \r
1862       // add the track to the TObjArray\r
1863       tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));   \r
1864 \r
1865       delete trackTPC;\r
1866     }//track loop\r
1867   }// ESD analysis\r
1868 \r
1869   else if(gAnalysisLevel == "MC"){\r
1870     if(!event) {\r
1871       AliError("mcEvent not available");\r
1872       return 0x0;\r
1873     }\r
1874 \r
1875     AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);\r
1876     if(gMCEvent) {\r
1877       // Loop over tracks in event\r
1878       for (Int_t iTracks = 0; iTracks < gMCEvent->GetNumberOfPrimaries(); iTracks++) {\r
1879         AliMCParticle* track = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iTracks));\r
1880         if (!track) {\r
1881           AliError(Form("Could not receive particle %d", iTracks));\r
1882           continue;\r
1883         }\r
1884         \r
1885         //exclude non stable particles\r
1886         if(!(gMCEvent->IsPhysicalPrimary(iTracks))) continue;\r
1887         \r
1888         vCharge = track->Charge();\r
1889         vEta    = track->Eta();\r
1890         vPt     = track->Pt();\r
1891         vY      = track->Y();\r
1892         \r
1893         if( vPt < fPtMin || vPt > fPtMax)      \r
1894           continue;\r
1895         if (!fUsePID) {\r
1896           if( vEta < fEtaMin || vEta > fEtaMax)  continue;\r
1897         }\r
1898         else if (fUsePID){\r
1899           if( vY < fEtaMin || vY > fEtaMax)  continue;\r
1900         }\r
1901 \r
1902         // Remove neutral tracks\r
1903         if( vCharge == 0 ) continue;\r
1904         \r
1905         //analyze one set of particles\r
1906         if(fUseMCPdgCode) {\r
1907           TParticle *particle = track->Particle();\r
1908           if(!particle) continue;\r
1909           \r
1910           Int_t gPdgCode = particle->GetPdgCode();\r
1911           if(TMath::Abs(fPDGCodeToBeAnalyzed) != TMath::Abs(gPdgCode)) \r
1912             continue;\r
1913         }\r
1914         \r
1915         //Use the acceptance parameterization\r
1916         if(fAcceptanceParameterization) {\r
1917           Double_t gRandomNumber = gRandom->Rndm();\r
1918           if(gRandomNumber > fAcceptanceParameterization->Eval(track->Pt())) \r
1919             continue;\r
1920         }\r
1921         \r
1922         //Exclude resonances\r
1923         if(fExcludeResonancesInMC) {\r
1924           TParticle *particle = track->Particle();\r
1925           if(!particle) continue;\r
1926           \r
1927           Bool_t kExcludeParticle = kFALSE;\r
1928           Int_t gMotherIndex = particle->GetFirstMother();\r
1929           if(gMotherIndex != -1) {\r
1930             AliMCParticle* motherTrack = dynamic_cast<AliMCParticle *>(event->GetTrack(gMotherIndex));\r
1931             if(motherTrack) {\r
1932               TParticle *motherParticle = motherTrack->Particle();\r
1933               if(motherParticle) {\r
1934                 Int_t pdgCodeOfMother = motherParticle->GetPdgCode();\r
1935                 //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {\r
1936                 if(pdgCodeOfMother == 113  // rho0\r
1937                    || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+\r
1938                    // || pdgCodeOfMother == 221  // eta\r
1939                    // || pdgCodeOfMother == 331  // eta'\r
1940                    // || pdgCodeOfMother == 223  // omega\r
1941                    // || pdgCodeOfMother == 333  // phi\r
1942                    || pdgCodeOfMother == 311  || pdgCodeOfMother == -311 // K0\r
1943                    // || pdgCodeOfMother == 313  || pdgCodeOfMother == -313 // K0*\r
1944                    // || pdgCodeOfMother == 323  || pdgCodeOfMother == -323 // K+*\r
1945                    || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda\r
1946                    || pdgCodeOfMother == 111  // pi0 Dalitz\r
1947                    ) {\r
1948                   kExcludeParticle = kTRUE;\r
1949                 }\r
1950               }\r
1951             }\r
1952           }\r
1953           \r
1954           //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi\r
1955           if(kExcludeParticle) continue;\r
1956         }\r
1957 \r
1958         //Exclude electrons with PDG\r
1959         if(fExcludeElectronsInMC) {\r
1960           \r
1961           TParticle *particle = track->Particle();\r
1962           \r
1963           if (particle){ \r
1964             if(TMath::Abs(particle->GetPdgCode()) == 11) continue;\r
1965           }\r
1966         }\r
1967       \r
1968         vPhi    = track->Phi();\r
1969         //Printf("phi (before): %lf",vPhi);\r
1970         \r
1971         fHistPt->Fill(vPt,gCentrality);\r
1972         fHistEta->Fill(vEta,gCentrality);\r
1973         fHistPhi->Fill(vPhi,gCentrality);\r
1974         if(vCharge > 0)      fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);\r
1975         else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r
1976         //fHistPhi->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r
1977         fHistRapidity->Fill(vY,gCentrality);\r
1978         //if(vCharge > 0) fHistPhiPos->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r
1979         //else if(vCharge < 0) fHistPhiNeg->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r
1980         if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
1981         else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
1982         \r
1983         //Flow after burner\r
1984         if(fUseFlowAfterBurner) {\r
1985           Double_t precisionPhi = 0.001;\r
1986           Int_t maxNumberOfIterations = 100;\r
1987           \r
1988           Double_t phi0 = vPhi;\r
1989           Double_t gV2 = fDifferentialV2->Eval(vPt);\r
1990           \r
1991           for (Int_t j = 0; j < maxNumberOfIterations; j++) {\r
1992             Double_t phiprev = vPhi;\r
1993             Double_t fl = vPhi - phi0 + gV2*TMath::Sin(2.*(vPhi - gReactionPlane*TMath::DegToRad()));\r
1994             Double_t fp = 1.0 + 2.0*gV2*TMath::Cos(2.*(vPhi - gReactionPlane*TMath::DegToRad())); \r
1995             vPhi -= fl/fp;\r
1996             if (TMath::AreEqualAbs(phiprev,vPhi,precisionPhi)) break;\r
1997           }\r
1998           //Printf("phi (after): %lf\n",vPhi);\r
1999           Double_t vDeltaphiBefore = phi0 - gReactionPlane*TMath::DegToRad();\r
2000           if(vDeltaphiBefore < 0) vDeltaphiBefore += 2*TMath::Pi();\r
2001           fHistPhiBefore->Fill(vDeltaphiBefore,gCentrality);\r
2002           \r
2003           Double_t vDeltaphiAfter = vPhi - gReactionPlane*TMath::DegToRad();\r
2004           if(vDeltaphiAfter < 0) vDeltaphiAfter += 2*TMath::Pi();\r
2005           fHistPhiAfter->Fill(vDeltaphiAfter,gCentrality);\r
2006           \r
2007         }\r
2008         \r
2009         //vPhi *= TMath::RadToDeg();\r
2010         \r
2011       //=======================================correction\r
2012       Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  \r
2013       //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);\r
2014 \r
2015         tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction)); \r
2016       } //track loop\r
2017     }//MC event object\r
2018   }//MC\r
2019   \r
2020   return tracksAccepted;  \r
2021 }\r
2022 \r
2023 //________________________________________________________________________\r
2024 TObjArray* AliAnalysisTaskBFPsi::GetShuffledTracks(TObjArray *tracks, Double_t gCentrality){\r
2025   // Clones TObjArray and returns it with tracks after shuffling the charges\r
2026 \r
2027   TObjArray* tracksShuffled = new TObjArray;\r
2028   tracksShuffled->SetOwner(kTRUE);\r
2029 \r
2030   vector<Short_t> *chargeVector = new vector<Short_t>;   //original charge of accepted tracks \r
2031 \r
2032   for (Int_t i=0; i<tracks->GetEntriesFast(); i++)\r
2033   {\r
2034     AliVParticle* track = (AliVParticle*) tracks->At(i);\r
2035     chargeVector->push_back(track->Charge());\r
2036   }  \r
2037  \r
2038   random_shuffle(chargeVector->begin(), chargeVector->end());\r
2039   \r
2040   for(Int_t i = 0; i < tracks->GetEntriesFast(); i++){\r
2041     AliVParticle* track = (AliVParticle*) tracks->At(i);\r
2042     //==============================correction\r
2043     Double_t correction = GetTrackbyTrackCorrectionMatrix(track->Eta(), track->Phi(),track->Pt(), chargeVector->at(i), gCentrality);\r
2044     //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);\r
2045     tracksShuffled->Add(new AliBFBasicParticle(track->Eta(), track->Phi(), track->Pt(),chargeVector->at(i), correction));\r
2046   }\r
2047 \r
2048   delete chargeVector;\r
2049    \r
2050   return tracksShuffled;\r
2051 }\r
2052 \r
2053 \r
2054 //________________________________________________________________________\r
2055 void  AliAnalysisTaskBFPsi::FinishTaskOutput(){\r
2056   //Printf("END BF");\r
2057 \r
2058   if (!fBalance) {\r
2059     AliError("fBalance not available");\r
2060     return;\r
2061   }  \r
2062   if(fRunShuffling) {\r
2063     if (!fShuffledBalance) {\r
2064       AliError("fShuffledBalance not available");\r
2065       return;\r
2066     }\r
2067   }\r
2068 \r
2069 }\r
2070 \r
2071 //________________________________________________________________________\r
2072 void AliAnalysisTaskBFPsi::Terminate(Option_t *) {\r
2073   // Draw result to the screen\r
2074   // Called once at the end of the query\r
2075 \r
2076   // not implemented ...\r
2077 \r
2078 }\r
2079 \r