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