PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskBFPsi.cxx

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