]> 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 0f2b671b06535259f841ec8ee6765455492c498a..5255a44f73ce3da5fc0c7fa065d505bb87a8c976 100755 (executable)
@@ -20,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
@@ -126,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
@@ -138,10 +149,19 @@ AliAnalysisTaskBFPsi::AliAnalysisTaskBFPsi(const char *name)
   fExcludeResonancesInMC(kFALSE),\r
   fUseMCPdgCode(kFALSE),\r
   fPDGCodeToBeAnalyzed(-1),\r
-  fEventClass("EventPlane") {\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
@@ -256,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
@@ -447,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
@@ -457,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
@@ -484,7 +566,7 @@ 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
@@ -493,17 +575,17 @@ void AliAnalysisTaskBFPsi::UserExec(Option_t *) {
 \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
@@ -512,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
@@ -534,10 +616,10 @@ 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
@@ -584,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
@@ -595,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
@@ -633,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
@@ -651,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
@@ -672,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
@@ -700,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
@@ -733,38 +818,41 @@ Double_t AliAnalysisTaskBFPsi::GetRefMultiOrCentrality(AliVEvent *event){
     // Checks the Event cuts\r
     // Fills Event statistics histograms\r
   \r
-  Float_t fCentrality = -1.;\r
+  Float_t gCentrality = -1.;\r
   Double_t fMultiplicity = -100.;\r
   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
   if(fEventClass == "Centrality"){\r
-    Bool_t isSelectedMain = kTRUE;\r
     \r
 \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
       }//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
-      fCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
+      gCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
     }//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
   }//End if "Centrality"\r
   if(fEventClass=="Multiplicity"&&gAnalysisLevel == "ESD"){\r
-    fMultiplicity = fESDtrackCuts->GetReferenceMultiplicity(dynamic_cast<AliESDEvent*>(event), AliESDtrackCuts::kTrackletsITSTPC,0.5);\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
@@ -776,7 +864,7 @@ Double_t AliAnalysisTaskBFPsi::GetRefMultiOrCentrality(AliVEvent *event){
   if(fEventClass=="Multiplicity"){\r
     lReturnVal = fMultiplicity;\r
   }else if(fEventClass=="Centrality"){\r
-    lReturnVal = fCentrality;\r
+    lReturnVal = gCentrality;\r
   }\r
   return lReturnVal;\r
 }\r
@@ -793,12 +881,13 @@ 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
@@ -817,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
@@ -827,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
@@ -847,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
@@ -882,16 +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
+      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
@@ -1062,16 +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
-      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
@@ -1149,15 +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
-      //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
@@ -1177,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
@@ -1194,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
@@ -1212,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