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