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