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