PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskBFPsi.cxx

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