]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskBFPsi.cxx
moving from AliGenHijingEventHeader to AliCollisionGeometry also for Event plane...
[u/mrichter/AliRoot.git] / PWGCF / EBYE / BalanceFunctions / AliAnalysisTaskBFPsi.cxx
index 3d36a1d7a0b0d38cec42db6ea5d0310603f36af6..5255a44f73ce3da5fc0c7fa065d505bb87a8c976 100755 (executable)
@@ -5,7 +5,6 @@
 #include "TGraphErrors.h"\r
 #include "TH1F.h"\r
 #include "TH2F.h"\r
-#include "TH3F.h"\r
 #include "TH2D.h"                  \r
 #include "TH3D.h"\r
 #include "TArrayF.h"\r
@@ -21,8 +20,8 @@
 #include "AliAODEvent.h"\r
 #include "AliAODTrack.h"\r
 #include "AliAODInputHandler.h"\r
+#include "AliCollisionGeometry.h"\r
 #include "AliGenEventHeader.h"\r
-#include "AliGenHijingEventHeader.h"\r
 #include "AliMCEventHandler.h"\r
 #include "AliMCEvent.h"\r
 #include "AliStack.h"\r
@@ -80,8 +79,6 @@ AliAnalysisTaskBFPsi::AliAnalysisTaskBFPsi(const char *name)
   fHistEta(0),\r
   fHistRapidity(0),\r
   fHistPhi(0),\r
-  fHistEtaPhiPos(0),\r
-  fHistEtaPhiNeg(0),\r
   fHistPhiBefore(0),\r
   fHistPhiAfter(0),\r
   fHistPhiPos(0),\r
@@ -129,8 +126,19 @@ AliAnalysisTaskBFPsi::AliAnalysisTaskBFPsi(const char *name)
   nAODtrackCutBit(128),\r
   fPtMin(0.3),\r
   fPtMax(1.5),\r
+  fPtMinForCorrections(0.3),//=================================correction\r
+  fPtMaxForCorrections(1.5),//=================================correction\r
+  fPtBinForCorrections(36), //=================================correction\r
   fEtaMin(-0.8),\r
   fEtaMax(-0.8),\r
+  fEtaMinForCorrections(-0.8),//=================================correction\r
+  fEtaMaxForCorrections(0.8),//=================================correction\r
+  fEtaBinForCorrections(16), //=================================correction\r
+  fPhiMin(0.),\r
+  fPhiMax(360.),\r
+  fPhiMinForCorrections(0.),//=================================correction\r
+  fPhiMaxForCorrections(360.),//=================================correction\r
+  fPhiBinForCorrections(100), //=================================correction\r
   fDCAxyCut(-1),\r
   fDCAzCut(-1),\r
   fTPCchi2Cut(-1),\r
@@ -140,10 +148,20 @@ AliAnalysisTaskBFPsi::AliAnalysisTaskBFPsi(const char *name)
   fUseFlowAfterBurner(kFALSE),\r
   fExcludeResonancesInMC(kFALSE),\r
   fUseMCPdgCode(kFALSE),\r
-  fPDGCodeToBeAnalyzed(-1) {\r
+  fPDGCodeToBeAnalyzed(-1),\r
+  fEventClass("EventPlane") \r
+{\r
   // Constructor\r
   // Define input and output slots here\r
   // Input slot #0 works with a TChain\r
+\r
+  //======================================================correction\r
+  for (Int_t i=0; i<=kCENTRALITY; i++){\r
+    fHistMatrixCorrectionPlus[i] = NULL; \r
+    fHistMatrixCorrectionMinus[i] = NULL; \r
+  }\r
+  //=====================================================correction\r
+\r
   DefineInput(0, TChain::Class());\r
   // Output slot #0 writes into a TH1 container\r
   DefineOutput(1, TList::Class());\r
@@ -174,8 +192,6 @@ AliAnalysisTaskBFPsi::~AliAnalysisTaskBFPsi() {
   // delete fHistPt;\r
   // delete fHistEta;\r
   // delete fHistPhi;\r
-  // delete fHistEtaPhiPos;\r
-  // delete fHistEtaPhiNeg;\r
   // delete fHistV0M;\r
 }\r
 \r
@@ -192,6 +208,7 @@ void AliAnalysisTaskBFPsi::UserCreateOutputObjects() {
   if(!fBalance) {\r
     fBalance = new AliBalancePsi();\r
     fBalance->SetAnalysisLevel("ESD");\r
+    fBalance->SetEventClass(fEventClass);\r
     //fBalance->SetNumberOfBins(-1,16);\r
     //fBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);\r
   }\r
@@ -199,6 +216,7 @@ void AliAnalysisTaskBFPsi::UserCreateOutputObjects() {
     if(!fShuffledBalance) {\r
       fShuffledBalance = new AliBalancePsi();\r
       fShuffledBalance->SetAnalysisLevel("ESD");\r
+      fShuffledBalance->SetEventClass(fEventClass);\r
       //fShuffledBalance->SetNumberOfBins(-1,16);\r
       //fShuffledBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);\r
     }\r
@@ -207,6 +225,7 @@ void AliAnalysisTaskBFPsi::UserCreateOutputObjects() {
     if(!fMixedBalance) {\r
       fMixedBalance = new AliBalancePsi();\r
       fMixedBalance->SetAnalysisLevel("ESD");\r
+      fMixedBalance->SetEventClass(fEventClass);\r
     }\r
   }\r
 \r
@@ -257,10 +276,10 @@ void AliAnalysisTaskBFPsi::UserCreateOutputObjects() {
     fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());\r
   fList->Add(fHistCentStats);\r
 \r
-  fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",130,0,130);\r
+  fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",1025,0,1025);\r
   fList->Add(fHistTriggerStats);\r
 \r
-  fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TrackFilterBit;N_{events}",130,0,130);\r
+  fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TrackFilterBit;N_{events}",16,0,16);\r
   fList->Add(fHistTrackStats);\r
 \r
   fHistNumberOfAcceptedTracks = new TH2F("fHistNumberOfAcceptedTracks",";N_{acc.};Centrality percentile;Entries",4001,-0.5,4000.5,220,-5,105);\r
@@ -293,10 +312,6 @@ void AliAnalysisTaskBFPsi::UserCreateOutputObjects() {
   fList->Add(fHistRapidity);\r
   fHistPhi  = new TH2F("fHistPhi","#phi distribution;#phi (rad);Centrality percentile",200,0.0,2.*TMath::Pi(),220,-5,105);\r
   fList->Add(fHistPhi);\r
-  fHistEtaPhiPos  = new TH3F("fHistEtaPhiPos","#eta-#phi distribution (+);#eta;#phi (rad);Centrality percentile",80,-2,2,72,0.0,2.*TMath::Pi(),220,-5,105);\r
-  fList->Add(fHistEtaPhiPos);\r
-  fHistEtaPhiNeg  = new TH3F("fHistEtaPhiNeg","#eta-#phi distribution (-);#eta;#phi (rad);Centrality percentile",80,-2,2,72,0.0,2.*TMath::Pi(),220,-5,105);\r
-  fList->Add(fHistEtaPhiNeg);\r
   fHistPhiBefore  = new TH2F("fHistPhiBefore","#phi distribution;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);\r
   fList->Add(fHistPhiBefore);\r
   fHistPhiAfter  = new TH2F("fHistPhiAfter","#phi distribution;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);\r
@@ -360,7 +375,7 @@ void AliAnalysisTaskBFPsi::UserCreateOutputObjects() {
     fListBFM->Add(fMixedBalance->GetHistNnn());\r
     fListBFM->Add(fMixedBalance->GetHistNpp());\r
     fListBFM->Add(fMixedBalance->GetHistNnp());\r
-  }  \r
+  }\r
   //}\r
 \r
 \r
@@ -452,7 +467,69 @@ void AliAnalysisTaskBFPsi::UserCreateOutputObjects() {
   if(fUsePID) PostData(5, fHistListPIDQA);       //PID\r
 \r
   TH1::AddDirectory(oldStatus);\r
+}\r
+\r
 \r
+//________________________________________________________________________\r
+void AliAnalysisTaskBFPsi::SetInputCorrection(const char* filename, \r
+                                             const char* gCollSystem) {\r
+  //Open files that will be used for correction\r
+  TString gCollidingSystem = gCollSystem;\r
+\r
+  //Open the input file\r
+  TFile *f = TFile::Open(filename);\r
+  if(!f->IsOpen()) {\r
+    Printf("File not found!!!");\r
+    return;\r
+  }\r
+  \r
+  //Get the TDirectoryFile and TList\r
+  TDirectoryFile *dir = dynamic_cast<TDirectoryFile *>(f->Get("PWGCFEbyE.outputBalanceFunctionEfficiencyAnalysis"));\r
+  if(!dir) {\r
+    Printf("TDirectoryFile not found!!!");\r
+    return;\r
+  }\r
+  \r
+  TString listEffName = "";\r
+  for (Int_t iCent = 0; iCent < kCENTRALITY; iCent++) {\r
+    listEffName = "listEfficiencyBF_";  \r
+    listEffName += (TString)((Int_t)(centralityArrayForPbPb[iCent]));\r
+    listEffName += "_";\r
+    listEffName += (TString)((Int_t)(centralityArrayForPbPb[iCent+1]));\r
+    TList *list = (TList*)dir->Get(listEffName.Data());\r
+    if(!list) {\r
+      Printf("TList Efficiency not found!!!");\r
+      return;\r
+    }\r
+    \r
+    TString histoName = "fHistMatrixCorrectionPlus";\r
+    fHistMatrixCorrectionPlus[iCent]= dynamic_cast<TH3D *>(list->FindObject(histoName.Data()));\r
+    if(!fHistMatrixCorrectionPlus[iCent]) {\r
+      Printf("fHist not found!!!");\r
+      return;\r
+    }\r
+    \r
+    histoName = "fHistMatrixCorrectionMinus";\r
+    fHistMatrixCorrectionMinus[iCent] = dynamic_cast<TH3D *>(list->FindObject(histoName.Data()));  \r
+    if(!fHistMatrixCorrectionMinus[iCent]) {\r
+      Printf("fHist not found!!!");\r
+      return; \r
+    }\r
+  }//loop over centralities: ONLY the PbPb case is covered\r
+\r
+  if(fHistMatrixCorrectionPlus[0]){\r
+    fEtaMinForCorrections = fHistMatrixCorrectionPlus[0]->GetXaxis()->GetXmin();\r
+    fEtaMaxForCorrections = fHistMatrixCorrectionPlus[0]->GetXaxis()->GetXmax();\r
+    fEtaBinForCorrections = fHistMatrixCorrectionPlus[0]->GetNbinsX();\r
+    \r
+    fPtMinForCorrections = fHistMatrixCorrectionPlus[0]->GetYaxis()->GetXmin();\r
+    fPtMaxForCorrections = fHistMatrixCorrectionPlus[0]->GetYaxis()->GetXmax();\r
+    fPtBinForCorrections = fHistMatrixCorrectionPlus[0]->GetNbinsY();\r
+    \r
+    fPhiMinForCorrections = fHistMatrixCorrectionPlus[0]->GetZaxis()->GetXmin();\r
+    fPhiMaxForCorrections = fHistMatrixCorrectionPlus[0]->GetZaxis()->GetXmax();\r
+    fPhiBinForCorrections = fHistMatrixCorrectionPlus[0]->GetNbinsZ();\r
+  }\r
 }\r
 \r
 //________________________________________________________________________\r
@@ -462,10 +539,10 @@ void AliAnalysisTaskBFPsi::UserExec(Option_t *) {
 \r
   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
   Int_t gNumberOfAcceptedTracks = 0;\r
-  Double_t fCentrality          = -1.;\r
+  Double_t gCentrality          = -1.;\r
   Double_t gReactionPlane       = -1.; \r
   Float_t bSign = 0.;\r
-\r
+  \r
   // get the event (for generator level: MCEvent())\r
   AliVEvent* eventMain = NULL;\r
   if(gAnalysisLevel == "MC") {\r
@@ -481,7 +558,6 @@ void AliAnalysisTaskBFPsi::UserExec(Option_t *) {
     AliError("eventMain not available");\r
     return;\r
   }\r
-\r
   \r
   // PID Response task active?\r
   if(fUsePID) {\r
@@ -490,23 +566,26 @@ void AliAnalysisTaskBFPsi::UserExec(Option_t *) {
   }\r
   \r
   // check event cuts and fill event histograms\r
-  if((fCentrality = IsEventAccepted(eventMain)) < 0){\r
+  if((gCentrality = IsEventAccepted(eventMain)) < 0){\r
     return;\r
   }\r
+  \r
+  //Compute Multiplicity or Centrality variable\r
+  Double_t lMultiplicityVar = GetRefMultiOrCentrality( eventMain );\r
 \r
   // get the reaction plane\r
   gReactionPlane = GetEventPlane(eventMain);\r
-  fHistEventPlane->Fill(gReactionPlane,fCentrality);\r
+  fHistEventPlane->Fill(gReactionPlane,gCentrality);\r
   if(gReactionPlane < 0){\r
     return;\r
   }\r
   \r
   // get the accepted tracks in main event\r
-  TObjArray *tracksMain = GetAcceptedTracks(eventMain,fCentrality,gReactionPlane);\r
+  TObjArray *tracksMain = GetAcceptedTracks(eventMain,gCentrality,gReactionPlane);\r
   gNumberOfAcceptedTracks = tracksMain->GetEntriesFast();\r
 \r
   //multiplicity cut (used in pp)\r
-  fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedTracks,fCentrality);\r
+  fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedTracks,gCentrality);\r
   if(fUseMultiplicity) {\r
     if((gNumberOfAcceptedTracks < fNumberOfAcceptedTracksMin)||(gNumberOfAcceptedTracks > fNumberOfAcceptedTracksMax))\r
       return;\r
@@ -515,7 +594,7 @@ void AliAnalysisTaskBFPsi::UserExec(Option_t *) {
   // store charges of all accepted tracks, shuffle and reassign (two extra loops!)\r
   TObjArray* tracksShuffled = NULL;\r
   if(fRunShuffling){\r
-    tracksShuffled = GetShuffledTracks(tracksMain);\r
+    tracksShuffled = GetShuffledTracks(tracksMain,gCentrality);\r
   }\r
   \r
   // Event mixing \r
@@ -537,15 +616,15 @@ void AliAnalysisTaskBFPsi::UserExec(Option_t *) {
       //    FillCorrelations(). Also nMix should be passed in, so a weight\r
       //    of 1./nMix can be applied.\r
       \r
-      AliEventPool* pool = fPoolMgr->GetEventPool(fCentrality, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane);\r
+      AliEventPool* pool = fPoolMgr->GetEventPool(gCentrality, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane);\r
       \r
       if (!pool){\r
-       AliFatal(Form("No pool found for centrality = %f, zVtx = %f, psi = %f", fCentrality, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane));\r
+       AliFatal(Form("No pool found for centrality = %f, zVtx = %f, psi = %f", gCentrality, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane));\r
       }\r
       else{\r
        \r
        //pool->SetDebug(1);\r
-       \r
+\r
        if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 || pool->GetCurrentNEvents() >= 5){ \r
          \r
          \r
@@ -561,7 +640,7 @@ void AliAnalysisTaskBFPsi::UserExec(Option_t *) {
          for (Int_t jMix=0; jMix<nMix; jMix++) \r
            {\r
              TObjArray* tracksMixed = pool->GetEvent(jMix);\r
-             fMixedBalance->CalculateBalance(gReactionPlane,tracksMain,tracksMixed,bSign); \r
+             fMixedBalance->CalculateBalance(gReactionPlane,tracksMain,tracksMixed,bSign,lMultiplicityVar);\r
            }\r
        }\r
        \r
@@ -573,11 +652,11 @@ void AliAnalysisTaskBFPsi::UserExec(Option_t *) {
     }//run mixing\r
   \r
   // calculate balance function\r
-  fBalance->CalculateBalance(gReactionPlane,tracksMain,NULL,bSign);\r
+  fBalance->CalculateBalance(gReactionPlane,tracksMain,NULL,bSign,lMultiplicityVar);\r
   \r
   // calculate shuffled balance function\r
   if(fRunShuffling && tracksShuffled != NULL) {\r
-    fShuffledBalance->CalculateBalance(gReactionPlane,tracksShuffled,NULL,bSign);\r
+    fShuffledBalance->CalculateBalance(gReactionPlane,tracksShuffled,NULL,bSign,lMultiplicityVar);\r
   }\r
 }      \r
 \r
@@ -587,10 +666,10 @@ Double_t AliAnalysisTaskBFPsi::IsEventAccepted(AliVEvent *event){
   // Fills Event statistics histograms\r
   \r
   Bool_t isSelectedMain = kTRUE;\r
-  Float_t fCentrality = -1.;\r
+  Float_t gCentrality = -1.;\r
   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
 \r
-  fHistEventStats->Fill(1,fCentrality); //all events\r
+  fHistEventStats->Fill(1,gCentrality); //all events\r
 \r
   // Event trigger bits\r
   fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());\r
@@ -598,14 +677,14 @@ Double_t AliAnalysisTaskBFPsi::IsEventAccepted(AliVEvent *event){
     isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
   \r
   if(isSelectedMain) {\r
-    fHistEventStats->Fill(2,fCentrality); //triggered events\r
+    fHistEventStats->Fill(2,gCentrality); //triggered events\r
     \r
     //Centrality stuff \r
     if(fUseCentrality) {\r
       if(gAnalysisLevel == "AOD") { //centrality in AOD header\r
        AliAODHeader *header = (AliAODHeader*) event->GetHeader();\r
        if(header){\r
-         fCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
+         gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
          \r
          // QA for centrality estimators\r
          fHistCentStats->Fill(0.,header->GetCentralityP()->GetCentralityPercentile("V0M"));\r
@@ -636,7 +715,7 @@ Double_t AliAnalysisTaskBFPsi::IsEventAccepted(AliVEvent *event){
 \r
       else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){ // centrality class for ESDs or MC-ESDs\r
        AliCentrality *centrality = event->GetCentrality();\r
-       fCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
+       gCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
 \r
        // QA for centrality estimators\r
        fHistCentStats->Fill(0.,centrality->GetCentralityPercentile("V0M"));\r
@@ -654,14 +733,16 @@ Double_t AliAnalysisTaskBFPsi::IsEventAccepted(AliVEvent *event){
       }//ESD\r
       else if(gAnalysisLevel == "MC"){\r
        Double_t gImpactParameter = 0.;\r
-       AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(dynamic_cast<AliMCEvent*>(event)->GenEventHeader());\r
-       if(headerH){\r
-         gImpactParameter = headerH->ImpactParameter();\r
-         fCentrality      = gImpactParameter;\r
-       }//MC header\r
+       if(dynamic_cast<AliMCEvent*>(event)){\r
+         AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(dynamic_cast<AliMCEvent*>(event)->GenEventHeader());\r
+         if(headerH){\r
+           gImpactParameter = headerH->ImpactParameter();\r
+           gCentrality      = gImpactParameter;\r
+         }//MC header\r
+       }//MC event cast\r
       }//MC\r
       else{\r
-       fCentrality = -1.;\r
+       gCentrality = -1.;\r
       }\r
     }\r
     \r
@@ -675,18 +756,19 @@ Double_t AliAnalysisTaskBFPsi::IsEventAccepted(AliVEvent *event){
        //gVertexArray.At(0),\r
        //gVertexArray.At(1),\r
        //gVertexArray.At(2));\r
-       fHistEventStats->Fill(3,fCentrality); //events with a proper vertex\r
+       fHistEventStats->Fill(3,gCentrality); //events with a proper vertex\r
        if(TMath::Abs(gVertexArray.At(0)) < fVxMax) {\r
          if(TMath::Abs(gVertexArray.At(1)) < fVyMax) {\r
            if(TMath::Abs(gVertexArray.At(2)) < fVzMax) {\r
-             fHistEventStats->Fill(4,fCentrality); //analayzed events\r
+             fHistEventStats->Fill(4,gCentrality); //analayzed events\r
              fHistVx->Fill(gVertexArray.At(0));\r
              fHistVy->Fill(gVertexArray.At(1));\r
-             fHistVz->Fill(gVertexArray.At(2),fCentrality);\r
+             fHistVz->Fill(gVertexArray.At(2),gCentrality);\r
              \r
              // take only events inside centrality class\r
-             if((fImpactParameterMin < fCentrality) && (fImpactParameterMax > fCentrality)){\r
-               return fCentrality;         \r
+             if((fImpactParameterMin < gCentrality) && (fImpactParameterMax > gCentrality)){\r
+               fHistEventStats->Fill(5,gCentrality); //events with correct centrality\r
+               return gCentrality;         \r
              }//centrality class\r
            }//Vz cut\r
          }//Vy cut\r
@@ -703,19 +785,19 @@ Double_t AliAnalysisTaskBFPsi::IsEventAccepted(AliVEvent *event){
        vertex->GetCovarianceMatrix(fCov);\r
        if(vertex->GetNContributors() > 0) {\r
          if(fCov[5] != 0) {\r
-           fHistEventStats->Fill(3,fCentrality); //events with a proper vertex\r
+           fHistEventStats->Fill(3,gCentrality); //events with a proper vertex\r
            if(TMath::Abs(vertex->GetX()) < fVxMax) {\r
              if(TMath::Abs(vertex->GetY()) < fVyMax) {\r
                if(TMath::Abs(vertex->GetZ()) < fVzMax) {\r
-                 fHistEventStats->Fill(4,fCentrality); //analyzed events\r
+                 fHistEventStats->Fill(4,gCentrality); //analyzed events\r
                  fHistVx->Fill(vertex->GetX());\r
                  fHistVy->Fill(vertex->GetY());\r
-                 fHistVz->Fill(vertex->GetZ(),fCentrality);\r
+                 fHistVz->Fill(vertex->GetZ(),gCentrality);\r
                  \r
                  // take only events inside centrality class\r
-                 if((fCentrality > fCentralityPercentileMin) && (fCentrality < fCentralityPercentileMax)){\r
-                   fHistEventStats->Fill(5,fCentrality); //events with correct centrality\r
-                   return fCentrality;         \r
+                 if((gCentrality > fCentralityPercentileMin) && (gCentrality < fCentralityPercentileMax)){\r
+                   fHistEventStats->Fill(5,gCentrality); //events with correct centrality\r
+                   return gCentrality;         \r
                  }//centrality class\r
                }//Vz cut\r
              }//Vy cut\r
@@ -730,6 +812,63 @@ Double_t AliAnalysisTaskBFPsi::IsEventAccepted(AliVEvent *event){
   return -1;\r
 }\r
 \r
+\r
+//________________________________________________________________________\r
+Double_t AliAnalysisTaskBFPsi::GetRefMultiOrCentrality(AliVEvent *event){\r
+    // Checks the Event cuts\r
+    // Fills Event statistics histograms\r
+  \r
+  Float_t gCentrality = -1.;\r
+  Double_t fMultiplicity = -100.;\r
+  TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
+  if(fEventClass == "Centrality"){\r
+    \r
+\r
+    if(gAnalysisLevel == "AOD") { //centrality in AOD header\r
+      AliAODHeader *header = (AliAODHeader*) event->GetHeader();\r
+      if(header){\r
+        gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
+      }//AOD header\r
+    }//AOD\r
+    \r
+    else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){ // centrality class for ESDs or MC-ESDs\r
+      AliCentrality *centrality = event->GetCentrality();\r
+      gCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
+    }//ESD\r
+    else if(gAnalysisLevel == "MC"){\r
+      Double_t gImpactParameter = 0.;\r
+      if(dynamic_cast<AliMCEvent*>(event)){\r
+       AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(dynamic_cast<AliMCEvent*>(event)->GenEventHeader());      \r
+       if(headerH){\r
+         gImpactParameter = headerH->ImpactParameter();\r
+         gCentrality      = gImpactParameter;\r
+       }//MC header\r
+      }//MC event cast\r
+    }//MC\r
+    else{\r
+      gCentrality = -1.;\r
+    }\r
+  }//End if "Centrality"\r
+  if(fEventClass=="Multiplicity"&&gAnalysisLevel == "ESD"){\r
+    if(dynamic_cast<AliESDEvent*>(event)){\r
+      fMultiplicity = fESDtrackCuts->GetReferenceMultiplicity(dynamic_cast<AliESDEvent*>(event), AliESDtrackCuts::kTrackletsITSTPC,0.5);\r
+    }//AliESDevent cast\r
+  }\r
+  if(fEventClass=="Multiplicity"&&gAnalysisLevel != "ESD"){\r
+    AliAODHeader *header = (AliAODHeader*) event->GetHeader();\r
+    if(header){\r
+      fMultiplicity = header->GetRefMultiplicity();\r
+    }//AOD header\r
+  }\r
+  Double_t lReturnVal = -100;\r
+  if(fEventClass=="Multiplicity"){\r
+    lReturnVal = fMultiplicity;\r
+  }else if(fEventClass=="Centrality"){\r
+    lReturnVal = gCentrality;\r
+  }\r
+  return lReturnVal;\r
+}\r
+\r
 //________________________________________________________________________\r
 Double_t AliAnalysisTaskBFPsi::GetEventPlane(AliVEvent *event){\r
   // Get the event plane\r
@@ -742,16 +881,17 @@ Double_t AliAnalysisTaskBFPsi::GetEventPlane(AliVEvent *event){
 \r
   //MC: from reaction plane\r
   if(gAnalysisLevel == "MC"){\r
-   \r
-    AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(dynamic_cast<AliMCEvent*>(event)->GenEventHeader());\r
-    if (headerH) {\r
-      gReactionPlane = headerH->ReactionPlaneAngle();\r
-      //gReactionPlane *= TMath::RadToDeg();\r
-    }\r
+    if(dynamic_cast<AliMCEvent*>(event)){\r
+      AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(dynamic_cast<AliMCEvent*>(event)->GenEventHeader());    \r
+      if (headerH) {\r
+       gReactionPlane = headerH->ReactionPlaneAngle();\r
+       //gReactionPlane *= TMath::RadToDeg();\r
+      }//MC header\r
+    }//MC event cast\r
   }//MC\r
   \r
   // AOD,ESD,ESDMC: from VZERO Event Plane\r
-  else{  \r
+  else{\r
    \r
     AliEventplane *ep = event->GetEventplane();\r
     if(ep){ \r
@@ -766,7 +906,62 @@ Double_t AliAnalysisTaskBFPsi::GetEventPlane(AliVEvent *event){
 }\r
 \r
 //________________________________________________________________________\r
-TObjArray* AliAnalysisTaskBFPsi::GetAcceptedTracks(AliVEvent *event, Double_t fCentrality, Double_t gReactionPlane){\r
+Double_t AliAnalysisTaskBFPsi::GetTrackbyTrackCorrectionMatrix( Double_t vEta, \r
+                                                               Double_t vPhi, \r
+                                                               Double_t vPt, \r
+                                                               Short_t vCharge, \r
+                                                               Double_t gCentrality) {\r
+  // -- Get efficiency correction of particle dependent on (eta, phi, pt, charge, centrality) \r
+\r
+  Double_t correction = 1.;\r
+  Int_t binEta = 0, binPt = 0, binPhi = 0;\r
+\r
+  //Printf("EtaMAx: %lf - EtaMin: %lf - EtaBin: %lf", fEtaMaxForCorrections,fEtaMinForCorrections,fEtaBinForCorrections);\r
+  if(fEtaBinForCorrections != 0) {\r
+    Double_t widthEta = (fEtaMaxForCorrections - fEtaMinForCorrections)/fEtaBinForCorrections;\r
+    if(fEtaMaxForCorrections != fEtaMinForCorrections) \r
+      binEta = (Int_t)(vEta/widthEta)+1;\r
+  }\r
+\r
+  if(fPtBinForCorrections != 0) {\r
+    Double_t widthPt = (fPtMaxForCorrections - fPtMinForCorrections)/fPtBinForCorrections;\r
+    if(fPtMaxForCorrections != fPtMinForCorrections) \r
+      binPt = (Int_t)(vPt/widthPt) + 1;\r
+  }\r
\r
+  if(fPhiBinForCorrections != 0) {\r
+    Double_t widthPhi = (fPhiMaxForCorrections - fPhiMinForCorrections)/fPhiBinForCorrections;\r
+    if(fPhiMaxForCorrections != fPhiMinForCorrections) \r
+      binPhi = (Int_t)(vPhi/widthPhi)+ 1;\r
+  }\r
+\r
+  Int_t gCentralityInt = 1;\r
+  for (Int_t i=0; i<=kCENTRALITY; i++){\r
+    if((centralityArrayForPbPb[i] <= gCentrality)&&(gCentrality <= centralityArrayForPbPb[i+1]))\r
+      gCentralityInt = i;\r
+  }\r
+  \r
+  if(fHistMatrixCorrectionPlus[gCentralityInt]){\r
+    if (vCharge > 0) {\r
+      correction = fHistMatrixCorrectionPlus[gCentralityInt]->GetBinContent(fHistMatrixCorrectionPlus[gCentralityInt]->GetBin(binEta, binPt, binPhi));\r
+    }\r
+    if (vCharge < 0) {\r
+      correction = fHistMatrixCorrectionMinus[gCentralityInt]->GetBinContent(fHistMatrixCorrectionMinus[gCentralityInt]->GetBin(binEta, binPt, binPhi));\r
+    }\r
+  }\r
+  else {\r
+    correction = 1.;\r
+  }\r
+  if (correction == 0.) { \r
+    AliError(Form("Should not happen : bin content = 0. >> eta: %.2f | phi : %.2f | pt : %.2f | cent %d",vEta, vPhi, vPt, gCentralityInt)); \r
+    correction = 1.; \r
+  } \r
+    \r
+  return correction;\r
+}\r
+\r
+//________________________________________________________________________\r
+TObjArray* AliAnalysisTaskBFPsi::GetAcceptedTracks(AliVEvent *event, Double_t gCentrality, Double_t gReactionPlane){\r
   // Returns TObjArray with tracks after all track cuts (only for AOD!)\r
   // Fills QA histograms\r
 \r
@@ -776,7 +971,7 @@ TObjArray* AliAnalysisTaskBFPsi::GetAcceptedTracks(AliVEvent *event, Double_t fC
   TObjArray* tracksAccepted = new TObjArray;\r
   tracksAccepted->SetOwner(kTRUE);\r
 \r
-  Double_t vCharge;\r
+  Short_t vCharge;\r
   Double_t vEta;\r
   Double_t vY;\r
   Double_t vPhi;\r
@@ -796,7 +991,10 @@ TObjArray* AliAnalysisTaskBFPsi::GetAcceptedTracks(AliVEvent *event, Double_t fC
       \r
       // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C\r
       // take only TPC only tracks \r
-      fHistTrackStats->Fill(aodTrack->GetFilterMap());\r
+      //fHistTrackStats->Fill(aodTrack->GetFilterMap());\r
+      for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){\r
+       fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));\r
+      }\r
       if(!aodTrack->TestFilterBit(nAODtrackCutBit)) continue;\r
       \r
       vCharge = aodTrack->Charge();\r
@@ -831,18 +1029,19 @@ TObjArray* AliAnalysisTaskBFPsi::GetAcceptedTracks(AliVEvent *event, Double_t fC
       // fill QA histograms\r
       fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());\r
       fHistDCA->Fill(dcaZ,dcaXY);\r
-      fHistChi2->Fill(aodTrack->Chi2perNDF(),fCentrality);\r
-      fHistPt->Fill(vPt,fCentrality);\r
-      fHistEta->Fill(vEta,fCentrality);\r
-      fHistRapidity->Fill(vY,fCentrality);\r
-      if(vCharge > 0) fHistPhiPos->Fill(vPhi,fCentrality);\r
-      else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,fCentrality);\r
-      fHistPhi->Fill(vPhi,fCentrality);\r
-      if(vCharge > 0)      fHistEtaPhiPos->Fill(vEta,vPhi,fCentrality);\r
-      else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,fCentrality);\r
+      fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality);\r
+      fHistPt->Fill(vPt,gCentrality);\r
+      fHistEta->Fill(vEta,gCentrality);\r
+      fHistRapidity->Fill(vY,gCentrality);\r
+      if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
+      else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
+      fHistPhi->Fill(vPhi,gCentrality);\r
+      \r
+      //=======================================correction\r
+      Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  \r
       \r
       // add the track to the TObjArray\r
-      tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, 1.*vCharge));\r
+      tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));  \r
     }//track loop\r
   }// AOD analysis\r
 \r
@@ -1013,17 +1212,19 @@ TObjArray* AliAnalysisTaskBFPsi::GetAcceptedTracks(AliVEvent *event, Double_t fC
       vPt     = trackTPC->Pt();\r
       fHistClus->Fill(trackTPC->GetITSclusters(0),nClustersTPC);\r
       fHistDCA->Fill(b[1],b[0]);\r
-      fHistChi2->Fill(chi2PerClusterTPC,fCentrality);\r
-      fHistPt->Fill(vPt,fCentrality);\r
-      fHistEta->Fill(vEta,fCentrality);\r
-      fHistPhi->Fill(vPhi,fCentrality);\r
-      if(vCharge > 0)      fHistEtaPhiPos->Fill(vEta,vPhi,fCentrality);\r
-      else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,fCentrality);      fHistRapidity->Fill(vY,fCentrality);\r
-      if(vCharge > 0) fHistPhiPos->Fill(vPhi,fCentrality);\r
-      else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,fCentrality);\r
+      fHistChi2->Fill(chi2PerClusterTPC,gCentrality);\r
+      fHistPt->Fill(vPt,gCentrality);\r
+      fHistEta->Fill(vEta,gCentrality);\r
+      fHistPhi->Fill(vPhi,gCentrality);\r
+      fHistRapidity->Fill(vY,gCentrality);\r
+      if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
+      else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
       \r
+      //=======================================correction\r
+      Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  \r
+\r
       // add the track to the TObjArray\r
-      tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge));      \r
+      tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));   \r
 \r
       delete trackTPC;\r
     }//track loop\r
@@ -1101,16 +1302,15 @@ TObjArray* AliAnalysisTaskBFPsi::GetAcceptedTracks(AliVEvent *event, Double_t fC
       vPhi    = track->Phi();\r
       //Printf("phi (before): %lf",vPhi);\r
       \r
-      fHistPt->Fill(vPt,fCentrality);\r
-      fHistEta->Fill(vEta,fCentrality);\r
-      fHistPhi->Fill(vPhi,fCentrality);\r
-      if(vCharge > 0)      fHistEtaPhiPos->Fill(vEta,vPhi,fCentrality);\r
-      else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,fCentrality);      //fHistPhi->Fill(vPhi*TMath::RadToDeg(),fCentrality);\r
-      fHistRapidity->Fill(vY,fCentrality);\r
-      //if(vCharge > 0) fHistPhiPos->Fill(vPhi*TMath::RadToDeg(),fCentrality);\r
-      //else if(vCharge < 0) fHistPhiNeg->Fill(vPhi*TMath::RadToDeg(),fCentrality);\r
-      if(vCharge > 0) fHistPhiPos->Fill(vPhi,fCentrality);\r
-      else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,fCentrality);\r
+      fHistPt->Fill(vPt,gCentrality);\r
+      fHistEta->Fill(vEta,gCentrality);\r
+      fHistPhi->Fill(vPhi,gCentrality);\r
+      //fHistPhi->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r
+      fHistRapidity->Fill(vY,gCentrality);\r
+      //if(vCharge > 0) fHistPhiPos->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r
+      //else if(vCharge < 0) fHistPhiNeg->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r
+      if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
+      else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
       \r
       //Flow after burner\r
       if(fUseFlowAfterBurner) {\r
@@ -1130,16 +1330,20 @@ TObjArray* AliAnalysisTaskBFPsi::GetAcceptedTracks(AliVEvent *event, Double_t fC
        //Printf("phi (after): %lf\n",vPhi);\r
        Double_t vDeltaphiBefore = phi0 - gReactionPlane*TMath::DegToRad();\r
        if(vDeltaphiBefore < 0) vDeltaphiBefore += 2*TMath::Pi();\r
-       fHistPhiBefore->Fill(vDeltaphiBefore,fCentrality);\r
+       fHistPhiBefore->Fill(vDeltaphiBefore,gCentrality);\r
        \r
        Double_t vDeltaphiAfter = vPhi - gReactionPlane*TMath::DegToRad();\r
        if(vDeltaphiAfter < 0) vDeltaphiAfter += 2*TMath::Pi();\r
-       fHistPhiAfter->Fill(vDeltaphiAfter,fCentrality);\r
+       fHistPhiAfter->Fill(vDeltaphiAfter,gCentrality);\r
+      \r
       }\r
       \r
       //vPhi *= TMath::RadToDeg();\r
-                \r
-      tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge));\r
+\r
+      //=======================================correction\r
+      Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  \r
+      \r
+      tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction)); \r
       \r
     } //track loop\r
   }//MC\r
@@ -1147,7 +1351,7 @@ TObjArray* AliAnalysisTaskBFPsi::GetAcceptedTracks(AliVEvent *event, Double_t fC
   return tracksAccepted;  \r
 }\r
 //________________________________________________________________________\r
-TObjArray* AliAnalysisTaskBFPsi::GetShuffledTracks(TObjArray *tracks){\r
+TObjArray* AliAnalysisTaskBFPsi::GetShuffledTracks(TObjArray *tracks, Double_t gCentrality){\r
   // Clones TObjArray and returns it with tracks after shuffling the charges\r
 \r
   TObjArray* tracksShuffled = new TObjArray;\r
@@ -1165,7 +1369,9 @@ TObjArray* AliAnalysisTaskBFPsi::GetShuffledTracks(TObjArray *tracks){
   \r
   for(Int_t i = 0; i < tracks->GetEntriesFast(); i++){\r
     AliVParticle* track = (AliVParticle*) tracks->At(i);\r
-    tracksShuffled->Add(new AliBFBasicParticle(track->Eta(), track->Phi(), track->Pt(),chargeVector->at(i)));\r
+    //==============================correction\r
+    Double_t correction = GetTrackbyTrackCorrectionMatrix(track->Eta(), track->Phi(),track->Pt(), chargeVector->at(i), gCentrality);\r
+    tracksShuffled->Add(new AliBFBasicParticle(track->Eta(), track->Phi(), track->Pt(),chargeVector->at(i), correction));\r
   }\r
 \r
   delete chargeVector;\r