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