X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;ds=sidebyside;f=PWGCF%2FEBYE%2FBalanceFunctions%2FAliAnalysisTaskBFPsi.cxx;h=5255a44f73ce3da5fc0c7fa065d505bb87a8c976;hb=951cf00e4e35f6433ef1450aef81e12165487c57;hp=3d36a1d7a0b0d38cec42db6ea5d0310603f36af6;hpb=3ce45c12467b09f7d2e1b42cd8b9eaf844d82a51;p=u%2Fmrichter%2FAliRoot.git

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