]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG3/muondep/AliAnalysisTaskMuonResolution.cxx
Fixed warnings from FC
[u/mrichter/AliRoot.git] / PWG3 / muondep / AliAnalysisTaskMuonResolution.cxx
index 93a4681334f0ad95d481a1bc6110c3470af3a1df..d06d09ad695b5d16b15dd8a930ecb1d25ed0ebdc 100644 (file)
@@ -13,6 +13,8 @@
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
+/* $Id$ */
+
 // ROOT includes
 #include <TFile.h>
 #include <TTree.h>
@@ -25,6 +27,9 @@
 #include <Riostream.h>
 #include <TString.h>
 #include <TGeoManager.h>
+#include <TList.h>
+#include <TObjString.h>
+#include <TRegexp.h>
 
 // STEER includes
 #include "AliESDEvent.h"
 #include "AliGeomManager.h"
 
 // ANALYSIS includes
-#include "AliAnalysisTaskSE.h"
 #include "AliAnalysisDataSlot.h"
 #include "AliAnalysisManager.h"
 #include "AliInputEventHandler.h"
 #include "AliAnalysisTaskMuonResolution.h"
+#include "AliCentrality.h"
 
 // MUON includes
 #include "AliMUONCDB.h"
 #include "AliMUONTrackParam.h"
 #include "AliMUONTrackExtrap.h"
 #include "AliMUONVCluster.h"
+#include "AliMUONVDigit.h"
 #include "AliMUONGeometryTransformer.h"
 #include "AliMUONGeometryModuleTransformer.h"
 #include "AliMUONGeometryDetElement.h"
 #include "AliMpDEIterator.h"
+#include "AliMpSegmentation.h"
+#include "AliMpVSegmentation.h"
+#include "AliMpConstants.h"
+#include "AliMpDDLStore.h"
+#include "AliMpPad.h"
+#include "AliMpDetElement.h"
+#include "AliMpCathodType.h"
 
 #ifndef SafeDelete
 #define SafeDelete(x) if (x != NULL) { delete x; x = NULL; }
@@ -67,23 +80,35 @@ AliAnalysisTaskMuonResolution::AliAnalysisTaskMuonResolution() :
   AliAnalysisTaskSE(), 
   fResiduals(NULL),
   fResidualsVsP(NULL),
+  fResidualsVsCent(NULL),
   fLocalChi2(NULL),
   fChamberRes(NULL),
   fTrackRes(NULL),
   fCanvases(NULL),
+  fDefaultStorage(""),
   fNEvents(0),
+  fShowProgressBar(kFALSE),
+  fPrintClResPerCh(kFALSE),
+  fPrintClResPerDE(kFALSE),
+  fGaus(NULL),
   fMinMomentum(0.),
   fSelectPhysics(kFALSE),
   fMatchTrig(kFALSE),
+  fApplyAccCut(kFALSE),
+  fSelectTrigger(kFALSE),
+  fTriggerMask(0),
   fExtrapMode(1),
   fCorrectForSystematics(kTRUE),
+  fRemoveMonoCathCl(kFALSE),
+  fCheckAllPads(kFALSE),
   fOCDBLoaded(kFALSE),
   fNDE(0),
   fReAlign(kFALSE),
   fOldAlignStorage(""),
   fNewAlignStorage(""),
   fOldGeoTransformer(NULL),
-  fNewGeoTransformer(NULL)
+  fNewGeoTransformer(NULL),
+  fSelectTriggerClass(NULL)
 {
   /// Default constructor
 }
@@ -93,28 +118,42 @@ AliAnalysisTaskMuonResolution::AliAnalysisTaskMuonResolution(const char *name) :
   AliAnalysisTaskSE(name), 
   fResiduals(NULL),
   fResidualsVsP(NULL),
+  fResidualsVsCent(NULL),
   fLocalChi2(NULL),
   fChamberRes(NULL),
   fTrackRes(NULL),
   fCanvases(NULL),
+  fDefaultStorage("raw://"),
   fNEvents(0),
+  fShowProgressBar(kFALSE),
+  fPrintClResPerCh(kFALSE),
+  fPrintClResPerDE(kFALSE),
+  fGaus(NULL),
   fMinMomentum(0.),
   fSelectPhysics(kFALSE),
   fMatchTrig(kFALSE),
+  fApplyAccCut(kFALSE),
+  fSelectTrigger(kFALSE),
+  fTriggerMask(0),
   fExtrapMode(1),
   fCorrectForSystematics(kTRUE),
+  fRemoveMonoCathCl(kFALSE),
+  fCheckAllPads(kFALSE),
   fOCDBLoaded(kFALSE),
   fNDE(0),
   fReAlign(kFALSE),
   fOldAlignStorage(""),
   fNewAlignStorage(""),
   fOldGeoTransformer(NULL),
-  fNewGeoTransformer(NULL)
+  fNewGeoTransformer(NULL),
+  fSelectTriggerClass(NULL)
 {
   /// Constructor
   
   for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) SetStartingResolution(i, -1., -1.);
   
+  FitResiduals();
+  
   // Output slot #1 writes into a TObjArray container
   DefineOutput(1,TObjArray::Class());
   // Output slot #2 writes into a TObjArray container
@@ -125,19 +164,26 @@ AliAnalysisTaskMuonResolution::AliAnalysisTaskMuonResolution(const char *name) :
   DefineOutput(4,TObjArray::Class());
   // Output slot #5 writes into a TObjArray container
   DefineOutput(5,TObjArray::Class());
+  // Output slot #5 writes into a TObjArray container
+  DefineOutput(6,TObjArray::Class());
 }
 
 //________________________________________________________________________
 AliAnalysisTaskMuonResolution::~AliAnalysisTaskMuonResolution()
 {
   /// Destructor
-  SafeDelete(fResiduals);
-  SafeDelete(fResidualsVsP);
+  if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
+    SafeDelete(fResiduals);
+    SafeDelete(fResidualsVsP);
+    SafeDelete(fResidualsVsCent);
+    SafeDelete(fTrackRes);
+  }
   SafeDelete(fChamberRes);
-  SafeDelete(fTrackRes);
   SafeDelete(fCanvases);
+  SafeDelete(fGaus);
   SafeDelete(fOldGeoTransformer);
   SafeDelete(fNewGeoTransformer);
+  SafeDelete(fSelectTriggerClass);
 }
 
 //___________________________________________________________________________
@@ -148,10 +194,21 @@ void AliAnalysisTaskMuonResolution::UserCreateOutputObjects()
   // do it once the OCDB has been loaded (i.e. from NotifyRun())
   if (!fOCDBLoaded) return;
   
+  // set the list of trigger classes that can be selected to fill histograms (in case the physics selection is not used)
+  fSelectTriggerClass = new TList();
+  fSelectTriggerClass->SetOwner();
+  fSelectTriggerClass->AddLast(new TObjString(" CINT1B-ABCE-NOPF-ALL ")); fSelectTriggerClass->Last()->SetUniqueID(AliVEvent::kMB);
+  fSelectTriggerClass->AddLast(new TObjString(" CMUS1B-ABCE-NOPF-MUON ")); fSelectTriggerClass->Last()->SetUniqueID(AliVEvent::kMUON);
+  fSelectTriggerClass->AddLast(new TObjString(" CINT1-B-")); fSelectTriggerClass->Last()->SetUniqueID(AliVEvent::kMB);
+  fSelectTriggerClass->AddLast(new TObjString(" CMUS1-B-")); fSelectTriggerClass->Last()->SetUniqueID(AliVEvent::kMUON);
+  fSelectTriggerClass->AddLast(new TObjString(" CSH1-B-")); fSelectTriggerClass->Last()->SetUniqueID(AliVEvent::kHighMult);
+  
   fResiduals = new TObjArray(1000);
   fResiduals->SetOwner();
   fResidualsVsP = new TObjArray(1000);
   fResidualsVsP->SetOwner();
+  fResidualsVsCent = new TObjArray(1000);
+  fResidualsVsCent->SetOwner();
   fTrackRes = new TObjArray(1000);
   fTrackRes->SetOwner();
   TH2F* h2;
@@ -164,10 +221,12 @@ void AliAnalysisTaskMuonResolution::UserCreateOutputObjects()
     if (recoParam->GetDefaultBendingReso(i) > maxSigma[1]) maxSigma[1] = recoParam->GetDefaultBendingReso(i);
   }
   const char* axes[2] = {"X", "Y"};
-  const Int_t nBins = 2000;
+  const Int_t nBins = 5000;
   const Int_t nSigma = 10;
   const Int_t pNBins = 20;
   const Double_t pEdges[2] = {0., 50.};
+  Int_t nCentBins = 12;
+  Double_t centRange[2] = {-10., 110.};
   
   for (Int_t ia = 0; ia < 2; ia++) {
     
@@ -175,25 +234,28 @@ void AliAnalysisTaskMuonResolution::UserCreateOutputObjects()
     
     // List of residual histos
     h2 = new TH2F(Form("hResidual%sPerCh_ClusterIn",axes[ia]), Form("cluster-track residual-%s distribution per chamber (cluster attached to the track);chamber ID;#Delta_{%s} (cm)",axes[ia],axes[ia]), 10, 0.5, 10.5, nBins, -maxRes, maxRes);
-    fResiduals->AddAtAndExpand(h2, kResidualPerCh_ClusterIn+ia);
+    fResiduals->AddAtAndExpand(h2, kResidualPerChClusterIn+ia);
     h2 = new TH2F(Form("hResidual%sPerCh_ClusterOut",axes[ia]), Form("cluster-track residual-%s distribution per chamber (cluster not attached to the track);chamber ID;#Delta_{%s} (cm)",axes[ia],axes[ia]), 10, 0.5, 10.5, nBins, -2.*maxRes, 2.*maxRes);
-    fResiduals->AddAtAndExpand(h2, kResidualPerCh_ClusterOut+ia);
+    fResiduals->AddAtAndExpand(h2, kResidualPerChClusterOut+ia);
     
     h2 = new TH2F(Form("hResidual%sPerHalfCh_ClusterIn",axes[ia]), Form("cluster-track residual-%s distribution per half chamber (cluster attached to the track);half chamber ID;#Delta_{%s} (cm)",axes[ia],axes[ia]), 20, 0.5, 20.5, nBins, -maxRes, maxRes);
-    fResiduals->AddAtAndExpand(h2, kResidualPerHalfCh_ClusterIn+ia);
+    for (Int_t i = 0; i < 10; i++) { h2->GetXaxis()->SetBinLabel(2*i+1, Form("%d-I",i+1)); h2->GetXaxis()->SetBinLabel(2*i+2, Form("%d-O",i+1)); }
+    fResiduals->AddAtAndExpand(h2, kResidualPerHalfChClusterIn+ia);
     h2 = new TH2F(Form("hResidual%sPerHalfCh_ClusterOut",axes[ia]), Form("cluster-track residual-%s distribution per half chamber (cluster not attached to the track);half chamber ID;#Delta_{%s} (cm)",axes[ia],axes[ia]), 20, 0.5, 20.5, nBins, -2.*maxRes, 2.*maxRes);
-    fResiduals->AddAtAndExpand(h2, kResidualPerHalfCh_ClusterOut+ia);
+    for (Int_t i = 0; i < 10; i++) { h2->GetXaxis()->SetBinLabel(2*i+1, Form("%d-I",i+1)); h2->GetXaxis()->SetBinLabel(2*i+2, Form("%d-O",i+1)); }
+    fResiduals->AddAtAndExpand(h2, kResidualPerHalfChClusterOut+ia);
     
-    h2 = new TH2F(Form("hResidual%sPerDE_ClusterIn",axes[ia]), Form("cluster-track residual-%s distribution per DE (cluster not attached to the track);DE ID;#Delta_{%s} (cm)",axes[ia],axes[ia]), fNDE, 0.5, fNDE+0.5, nBins, -maxRes, maxRes);
+    h2 = new TH2F(Form("hResidual%sPerDE_ClusterIn",axes[ia]), Form("cluster-track residual-%s distribution per DE (cluster attached to the track);DE ID;#Delta_{%s} (cm)",axes[ia],axes[ia]), fNDE, 0.5, fNDE+0.5, nBins, -maxRes, maxRes);
     for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
-    fResiduals->AddAtAndExpand(h2, kResidualPerDE_ClusterIn+ia);
+    fResiduals->AddAtAndExpand(h2, kResidualPerDEClusterIn+ia);
     h2 = new TH2F(Form("hResidual%sPerDE_ClusterOut",axes[ia]), Form("cluster-track residual-%s distribution per DE (cluster not attached to the track);DE ID;#Delta_{%s} (cm)",axes[ia],axes[ia]), fNDE, 0.5, fNDE+0.5, nBins, -2.*maxRes, 2.*maxRes);
     for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
-    fResiduals->AddAtAndExpand(h2, kResidualPerDE_ClusterOut+ia);
+    fResiduals->AddAtAndExpand(h2, kResidualPerDEClusterOut+ia);
     
     h2 = new TH2F(Form("hTrackRes%sPerCh",axes[ia]), Form("track #sigma_{%s} per Ch;chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]), 10, 0.5, 10.5, nBins, 0., maxRes);
     fResiduals->AddAtAndExpand(h2, kTrackResPerCh+ia);
     h2 = new TH2F(Form("hTrackRes%sPerHalfCh",axes[ia]), Form("track #sigma_{%s} per half Ch;half chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]), 20, 0.5, 20.5, nBins, 0., maxRes);
+    for (Int_t i = 0; i < 10; i++) { h2->GetXaxis()->SetBinLabel(2*i+1, Form("%d-I",i+1)); h2->GetXaxis()->SetBinLabel(2*i+2, Form("%d-O",i+1)); }
     fResiduals->AddAtAndExpand(h2, kTrackResPerHalfCh+ia);
     h2 = new TH2F(Form("hTrackRes%sPerDE",axes[ia]), Form("track #sigma_{%s} per DE;DE ID;#sigma_{%s} (cm)",axes[ia],axes[ia]), fNDE, 0.5, fNDE+0.5, nBins, 0., maxRes);
     for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
@@ -202,6 +264,7 @@ void AliAnalysisTaskMuonResolution::UserCreateOutputObjects()
     h2 = new TH2F(Form("hMCS%sPerCh",axes[ia]), Form("MCS %s-dispersion of extrapolated track per Ch;chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]), 10, 0.5, 10.5, nBins, 0., 0.2);
     fResiduals->AddAtAndExpand(h2, kMCSPerCh+ia);
     h2 = new TH2F(Form("hMCS%sPerHalfCh",axes[ia]), Form("MCS %s-dispersion of extrapolated track per half Ch;half chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]), 20, 0.5, 20.5, nBins, 0., 0.2);
+    for (Int_t i = 0; i < 10; i++) { h2->GetXaxis()->SetBinLabel(2*i+1, Form("%d-I",i+1)); h2->GetXaxis()->SetBinLabel(2*i+2, Form("%d-O",i+1)); }
     fResiduals->AddAtAndExpand(h2, kMCSPerHalfCh+ia);
     h2 = new TH2F(Form("hMCS%sPerDE",axes[ia]), Form("MCS %s-dispersion of extrapolated track per DE;DE ID;#sigma_{%s} (cm)",axes[ia],axes[ia]), fNDE, 0.5, fNDE+0.5, nBins, 0., 0.2);
     for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
@@ -210,14 +273,29 @@ void AliAnalysisTaskMuonResolution::UserCreateOutputObjects()
     h2 = new TH2F(Form("hClusterRes2%sPerCh",axes[ia]), Form("cluster #sigma_{%s}^{2} per Ch;chamber ID;#sigma_{%s}^{2} (cm^{2})",axes[ia],axes[ia]), 10, 0.5, 10.5, nSigma*nBins, -0.1*maxRes*maxRes, maxRes*maxRes);
     fResiduals->AddAtAndExpand(h2, kClusterRes2PerCh+ia);
     
-    // List of residual vs. p histos
+    // List of residual vs. p or vs. centrality histos
     for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
       h2 = new TH2F(Form("hResidual%sInCh%dVsP_ClusterIn",axes[ia],i+1), Form("cluster-track residual-%s distribution in chamber %d versus momentum (cluster attached to the track);p (GeV/c^{2});#Delta_{%s} (cm)",axes[ia],i+1,axes[ia]), pNBins, pEdges[0], pEdges[1], nBins, -maxRes, maxRes);
-      fResidualsVsP->AddAtAndExpand(h2, kResidualInChVsP_ClusterIn+10*ia+i);
+      fResidualsVsP->AddAtAndExpand(h2, kResidualInChVsPClusterIn+10*ia+i);
       h2 = new TH2F(Form("hResidual%sInCh%dVsP_ClusterOut",axes[ia],i+1), Form("cluster-track residual-%s distribution in chamber %d versus momentum (cluster not attached to the track);p (GeV/c^{2});#Delta_{%s} (cm)",axes[ia],i+1,axes[ia]), pNBins, pEdges[0], pEdges[1], nBins, -2.*maxRes, 2.*maxRes);
-      fResidualsVsP->AddAtAndExpand(h2, kResidualInChVsP_ClusterOut+10*ia+i);
+      fResidualsVsP->AddAtAndExpand(h2, kResidualInChVsPClusterOut+10*ia+i);
+      
+      h2 = new TH2F(Form("hResidual%sInCh%dVsCent_ClusterIn",axes[ia],i+1), Form("cluster-track residual-%s distribution in chamber %d versus centrality (cluster attached to the track);centrality (%%);#Delta_{%s} (cm)",axes[ia],i+1,axes[ia]), nCentBins, centRange[0], centRange[1], nBins, -maxRes, maxRes);
+      fResidualsVsCent->AddAtAndExpand(h2, kResidualInChVsCentClusterIn+10*ia+i);
+      h2 = new TH2F(Form("hResidual%sInCh%dVsCent_ClusterOut",axes[ia],i+1), Form("cluster-track residual-%s distribution in chamber %d versus centrality (cluster not attached to the track);centrality (%%);#Delta_{%s} (cm)",axes[ia],i+1,axes[ia]), nCentBins, centRange[0], centRange[1], nBins, -2.*maxRes, 2.*maxRes);
+      fResidualsVsCent->AddAtAndExpand(h2, kResidualInChVsCentClusterOut+10*ia+i);
     }
     
+    h2 = new TH2F(Form("hResidual%sVsP_ClusterIn",axes[ia]), Form("cluster-track residual-%s distribution integrated over chambers versus momentum (cluster attached to the track);p (GeV/c^{2});#Delta_{%s} (cm)",axes[ia],axes[ia]), pNBins, pEdges[0], pEdges[1], nBins, -maxRes, maxRes);
+    fResidualsVsP->AddAtAndExpand(h2, kResidualVsPClusterIn+ia);
+    h2 = new TH2F(Form("hResidual%sVsP_ClusterOut",axes[ia]), Form("cluster-track residual-%s distribution integrated over chambers versus momentum (cluster not attached to the track);p (GeV/c^{2});#Delta_{%s} (cm)",axes[ia],axes[ia]), pNBins, pEdges[0], pEdges[1], nBins, -2.*maxRes, 2.*maxRes);
+    fResidualsVsP->AddAtAndExpand(h2, kResidualVsPClusterOut+ia);
+    
+    h2 = new TH2F(Form("hResidual%sVsCent_ClusterIn",axes[ia]), Form("cluster-track residual-%s distribution integrated over chambers versus centrality (cluster attached to the track);centrality (%%);#Delta_{%s} (cm)",axes[ia],axes[ia]), nCentBins, centRange[0], centRange[1], nBins, -maxRes, maxRes);
+    fResidualsVsCent->AddAtAndExpand(h2, kResidualVsCentClusterIn+ia);
+    h2 = new TH2F(Form("hResidual%sVsCent_ClusterOut",axes[ia]), Form("cluster-track residual-%s distribution integrated over chambers versus centrality (cluster not attached to the track);centrality (%%);#Delta_{%s} (cm)",axes[ia],axes[ia]), nCentBins, centRange[0], centRange[1], nBins, -2.*maxRes, 2.*maxRes);
+    fResidualsVsCent->AddAtAndExpand(h2, kResidualVsCentClusterOut+ia);
+    
     // local chi2
     h2 = new TH2F(Form("hLocalChi2%sPerCh",axes[ia]), Form("local chi2-%s distribution per chamber;chamber ID;local #chi^{2}_{%s}", axes[ia], axes[ia]), 10, 0.5, 10.5, 1000, 0., 25.);
     fResiduals->AddAtAndExpand(h2, kLocalChi2PerCh+ia);
@@ -252,7 +330,8 @@ void AliAnalysisTaskMuonResolution::UserCreateOutputObjects()
   // Post data at least once per task to ensure data synchronisation (required for merging)
   PostData(1, fResiduals);
   PostData(2, fResidualsVsP);
-  PostData(5, fTrackRes);
+  PostData(3, fResidualsVsCent);
+  PostData(6, fTrackRes);
 }
 
 //________________________________________________________________________
@@ -266,7 +345,19 @@ void AliAnalysisTaskMuonResolution::UserExec(Option_t *)
   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
   if (!esd) return;
   
-  if ((++fNEvents)%100 == 0) cout<<"\rEvent processing... "<<fNEvents<<"\r"<<flush;
+  if (fShowProgressBar && (++fNEvents)%100 == 0) cout<<"\rEvent processing... "<<fNEvents<<"\r"<<flush;
+  
+  // skip events that do not pass the physics selection if required
+  UInt_t triggerWord = (fInputHandler) ? fInputHandler->IsEventSelected() : 0;
+  if (fSelectPhysics && triggerWord == 0) return;
+  
+  // skip events that do not pass the trigger selection if required
+  TString firedTriggerClasses = esd->GetFiredTriggerClasses();
+  if (!fSelectPhysics) triggerWord = BuildTriggerWord(firedTriggerClasses);
+  if (fSelectTrigger && (triggerWord & fTriggerMask) == 0) return;
+  
+  // get the centrality
+  Float_t centrality = esd->GetCentrality()->GetCentralityPercentileUnchecked("V0M");
   
   // get tracker to refit
   AliMUONVTrackReconstructor* tracker = AliMUONESDInterface::GetTracker();
@@ -281,12 +372,14 @@ void AliAnalysisTaskMuonResolution::UserExec(Option_t *)
     // skip ghost tracks
     if (!esdTrack->ContainTrackerData()) continue;
     
-    // skip tracks that do not pass the physics selection if required
-    if (fSelectPhysics && fInputHandler && !fInputHandler->IsEventSelected()) continue;
-    
     // skip tracks not matched with trigger if required
     if (fMatchTrig && !esdTrack->ContainTriggerData()) continue;
     
+    // skip tracks that do not pass the acceptance cuts if required
+    Double_t thetaAbs = TMath::ATan(esdTrack->GetRAtAbsorberEnd()/505.) * TMath::RadToDeg();
+    Double_t eta = esdTrack->Eta();
+    if (fApplyAccCut && (thetaAbs < 2. || thetaAbs > 10. || eta < -4. || eta > -2.5)) continue;
+    
     // skip low momentum tracks
     if (esdTrack->PUncorrected() < fMinMomentum) continue;
     
@@ -360,6 +453,19 @@ void AliAnalysisTaskMuonResolution::UserExec(Option_t *)
       AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->UncheckedAt(iCluster));
       AliMUONTrackParam* previousTrackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->Before(trackParam));
       AliMUONTrackParam* nextTrackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->After(trackParam));
+      if (!previousTrackParam && !nextTrackParam) {
+       AliError(Form("unable to find a cluster neither before nor after the one at position %d !?!", iCluster));
+       track.RecursiveDump();
+       break;
+      }
+      
+      // remove mono-cathod clusters if required
+      if (fRemoveMonoCathCl) {
+       Bool_t hasBending, hasNonBending;
+       if (fCheckAllPads) CheckPads(trackParam->GetClusterPtr(), hasBending, hasNonBending);
+       else CheckPadsBelow(trackParam->GetClusterPtr(), hasBending, hasNonBending);
+       if (!hasBending || !hasNonBending) continue;
+      }
       
       // save current trackParam and remove it from the track
       AliMUONTrackParam currentTrackParam(*trackParam);
@@ -399,24 +505,29 @@ void AliAnalysisTaskMuonResolution::UserExec(Option_t *)
        if (tracker->RefitTrack(track, kFALSE)) {
          
          // fill histograms of residuals with cluster still attached to the track
-         ((TH2F*)fResiduals->UncheckedAt(kResidualPerCh_ClusterIn))->Fill(chId+1, deltaX);
-         ((TH2F*)fResiduals->UncheckedAt(kResidualPerCh_ClusterIn+1))->Fill(chId+1, deltaY);
-         ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfCh_ClusterIn))->Fill(halfChId+1, deltaX);
-         ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfCh_ClusterIn+1))->Fill(halfChId+1, deltaY);
-         ((TH2F*)fResiduals->UncheckedAt(kResidualPerDE_ClusterIn))->Fill(fDEIndices[deId], deltaX);
-         ((TH2F*)fResiduals->UncheckedAt(kResidualPerDE_ClusterIn+1))->Fill(fDEIndices[deId], deltaY);
-         ((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsP_ClusterIn+chId))->Fill(pUncorr, deltaX);
-         ((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsP_ClusterIn+10+chId))->Fill(pUncorr, deltaY);
+         ((TH2F*)fResiduals->UncheckedAt(kResidualPerChClusterIn))->Fill(chId+1, deltaX);
+         ((TH2F*)fResiduals->UncheckedAt(kResidualPerChClusterIn+1))->Fill(chId+1, deltaY);
+         ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfChClusterIn))->Fill(halfChId+1, deltaX);
+         ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfChClusterIn+1))->Fill(halfChId+1, deltaY);
+         ((TH2F*)fResiduals->UncheckedAt(kResidualPerDEClusterIn))->Fill(fDEIndices[deId], deltaX);
+         ((TH2F*)fResiduals->UncheckedAt(kResidualPerDEClusterIn+1))->Fill(fDEIndices[deId], deltaY);
+         ((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsPClusterIn+chId))->Fill(pUncorr, deltaX);
+         ((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsPClusterIn+10+chId))->Fill(pUncorr, deltaY);
+         ((TH2F*)fResidualsVsP->UncheckedAt(kResidualVsPClusterIn))->Fill(pUncorr, deltaX);
+         ((TH2F*)fResidualsVsP->UncheckedAt(kResidualVsPClusterIn+1))->Fill(pUncorr, deltaY);
+         ((TH2F*)fResidualsVsCent->UncheckedAt(kResidualInChVsCentClusterIn+chId))->Fill(centrality, deltaX);
+         ((TH2F*)fResidualsVsCent->UncheckedAt(kResidualInChVsCentClusterIn+10+chId))->Fill(centrality, deltaY);
+         ((TH2F*)fResidualsVsCent->UncheckedAt(kResidualVsCentClusterIn))->Fill(centrality, deltaX);
+         ((TH2F*)fResidualsVsCent->UncheckedAt(kResidualVsCentClusterIn+1))->Fill(centrality, deltaY);
          
          // find the track parameters closest to the current cluster position
          Double_t dZWithPrevious = (previousTrackParam) ? TMath::Abs(previousTrackParam->GetClusterPtr()->GetZ() - cluster->GetZ()) : FLT_MAX;
          Int_t previousChId = (previousTrackParam) ? previousTrackParam->GetClusterPtr()->GetChamberId() : -1;
          Double_t dZWithNext = (nextTrackParam) ? TMath::Abs(nextTrackParam->GetClusterPtr()->GetZ() - cluster->GetZ()) : FLT_MAX;
-         AliMUONTrackParam* startingTrackParam = 0x0;
-         if ((fExtrapMode == 0 && dZWithPrevious < dZWithNext) ||
+         AliMUONTrackParam* startingTrackParam = (nextTrackParam) ? nextTrackParam : previousTrackParam;
+         if ((fExtrapMode == 0 && previousTrackParam && dZWithPrevious < dZWithNext) ||
              (fExtrapMode == 1 && previousTrackParam && !(chId/2 == 2 && previousChId/2 == 1) &&
               !(chId/2 == 3 && previousChId/2 == 2))) startingTrackParam = previousTrackParam;
-         else startingTrackParam = nextTrackParam;
          
          // reset current parameters
          currentTrackParam.SetParameters(startingTrackParam->GetParameters());
@@ -445,14 +556,20 @@ void AliAnalysisTaskMuonResolution::UserExec(Option_t *)
            deltaY = cluster->GetY() - currentTrackParam.GetBendingCoor();
            
            // fill histograms with cluster not attached to the track
-           ((TH2F*)fResiduals->UncheckedAt(kResidualPerCh_ClusterOut))->Fill(chId+1, deltaX);
-           ((TH2F*)fResiduals->UncheckedAt(kResidualPerCh_ClusterOut+1))->Fill(chId+1, deltaY);
-           ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfCh_ClusterOut))->Fill(halfChId+1, deltaX);
-           ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfCh_ClusterOut+1))->Fill(halfChId+1, deltaY);
-           ((TH2F*)fResiduals->UncheckedAt(kResidualPerDE_ClusterOut))->Fill(fDEIndices[deId], deltaX);
-           ((TH2F*)fResiduals->UncheckedAt(kResidualPerDE_ClusterOut+1))->Fill(fDEIndices[deId], deltaY);
-           ((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsP_ClusterOut+chId))->Fill(pUncorr, deltaX);
-           ((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsP_ClusterOut+10+chId))->Fill(pUncorr, deltaY);
+           ((TH2F*)fResiduals->UncheckedAt(kResidualPerChClusterOut))->Fill(chId+1, deltaX);
+           ((TH2F*)fResiduals->UncheckedAt(kResidualPerChClusterOut+1))->Fill(chId+1, deltaY);
+           ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfChClusterOut))->Fill(halfChId+1, deltaX);
+           ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfChClusterOut+1))->Fill(halfChId+1, deltaY);
+           ((TH2F*)fResiduals->UncheckedAt(kResidualPerDEClusterOut))->Fill(fDEIndices[deId], deltaX);
+           ((TH2F*)fResiduals->UncheckedAt(kResidualPerDEClusterOut+1))->Fill(fDEIndices[deId], deltaY);
+           ((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsPClusterOut+chId))->Fill(pUncorr, deltaX);
+           ((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsPClusterOut+10+chId))->Fill(pUncorr, deltaY);
+           ((TH2F*)fResidualsVsP->UncheckedAt(kResidualVsPClusterOut))->Fill(pUncorr, deltaX);
+           ((TH2F*)fResidualsVsP->UncheckedAt(kResidualVsPClusterOut+1))->Fill(pUncorr, deltaY);
+           ((TH2F*)fResidualsVsCent->UncheckedAt(kResidualInChVsCentClusterOut+chId))->Fill(centrality, deltaX);
+           ((TH2F*)fResidualsVsCent->UncheckedAt(kResidualInChVsCentClusterOut+10+chId))->Fill(centrality, deltaY);
+           ((TH2F*)fResidualsVsCent->UncheckedAt(kResidualVsCentClusterOut))->Fill(centrality, deltaX);
+           ((TH2F*)fResidualsVsCent->UncheckedAt(kResidualVsCentClusterOut+1))->Fill(centrality, deltaY);
            ((TH2F*)fResiduals->UncheckedAt(kTrackResPerCh))->Fill(chId+1, TMath::Sqrt(trackResX2));
            ((TH2F*)fResiduals->UncheckedAt(kTrackResPerCh+1))->Fill(chId+1, TMath::Sqrt(trackResY2));
            ((TH2F*)fResiduals->UncheckedAt(kTrackResPerHalfCh))->Fill(halfChId+1, TMath::Sqrt(trackResX2));
@@ -483,7 +600,8 @@ void AliAnalysisTaskMuonResolution::UserExec(Option_t *)
   // Post final data. It will be written to a file with option "RECREATE"
   PostData(1, fResiduals);
   PostData(2, fResidualsVsP);
-  PostData(5, fTrackRes);
+  PostData(3, fResidualsVsCent);
+  PostData(6, fTrackRes);
 }
 
 //________________________________________________________________________
@@ -494,6 +612,7 @@ void AliAnalysisTaskMuonResolution::NotifyRun()
   if (fOCDBLoaded) return;
   
   AliCDBManager* cdbm = AliCDBManager::Instance();
+  cdbm->SetDefaultStorage(fDefaultStorage.Data());
   cdbm->SetRun(fCurrentRunNumber);
   
   if (!AliMUONCDB::LoadField()) return;
@@ -503,8 +622,6 @@ void AliAnalysisTaskMuonResolution::NotifyRun()
   AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam();
   if (!recoParam) return;
   
-  fOCDBLoaded = kTRUE;
-  
   AliMUONESDInterface::ResetTracker(recoParam);
   
   for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
@@ -527,7 +644,7 @@ void AliAnalysisTaskMuonResolution::NotifyRun()
   
   if (fReAlign) {
     
-    // recover default storage name
+    // recover default storage full name (raw:// cannot be used to set specific storage)
     TString defaultStorage(cdbm->GetDefaultStorage()->GetType());
     if (defaultStorage == "alien") defaultStorage += Form("://folder=%s", cdbm->GetDefaultStorage()->GetBaseFolder().Data());
     else defaultStorage += Form("://%s", cdbm->GetDefaultStorage()->GetBaseFolder().Data());
@@ -564,19 +681,27 @@ void AliAnalysisTaskMuonResolution::NotifyRun()
     
     // load geometry for track extrapolation to vertex
     if (cdbm->GetEntryCache()->Contains("GRP/Geometry/Data")) cdbm->UnloadFromCache("GRP/Geometry/Data");
+    if (cdbm->GetEntryCache()->Contains("MUON/Align/Data")) cdbm->UnloadFromCache("MUON/Align/Data");
     if (AliGeomManager::GetGeometry()) AliGeomManager::GetGeometry()->UnlockGeometry();
     AliGeomManager::LoadGeometry();
     if (!AliGeomManager::GetGeometry()) return;
+    AliGeomManager::ApplyAlignObjsFromCDB("MUON");
+    fNewGeoTransformer = new AliMUONGeometryTransformer();
+    fNewGeoTransformer->LoadGeometryData();
     
   }
   
-  // print starting chamber resolution
-  printf("\nstarting chamber resolution:\n");
-  printf(" - non-bending:");
-  for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %5.3f":", %5.3f",fClusterResNB[i]);
-  printf("\n -     bending:");
-  for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %6.4f":", %6.4f",fClusterResB[i]);
-  printf("\n\n");
+  // print starting chamber resolution if required
+  if (fPrintClResPerCh) {
+    printf("\nstarting chamber resolution:\n");
+    printf(" - non-bending:");
+    for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %5.3f":", %5.3f",fClusterResNB[i]);
+    printf("\n -     bending:");
+    for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %6.4f":", %6.4f",fClusterResB[i]);
+    printf("\n\n");
+  }
+  
+  fOCDBLoaded = kTRUE;
   
   UserCreateOutputObjects();
   
@@ -590,7 +715,8 @@ void AliAnalysisTaskMuonResolution::Terminate(Option_t *)
   // recover output objects
   fResiduals = static_cast<TObjArray*> (GetOutputData(1));
   fResidualsVsP = static_cast<TObjArray*> (GetOutputData(2));
-  fTrackRes = static_cast<TObjArray*> (GetOutputData(5));
+  fResidualsVsCent = static_cast<TObjArray*> (GetOutputData(3));
+  fTrackRes = static_cast<TObjArray*> (GetOutputData(6));
   if (!fResiduals || !fResidualsVsP || !fTrackRes) return;
   
   // summary graphs
@@ -603,7 +729,7 @@ void AliAnalysisTaskMuonResolution::Terminate(Option_t *)
   
   const char* axes[2] = {"X", "Y"};
   Double_t newClusterRes[2][10], newClusterResErr[2][10];
-  fNDE = ((TH2F*)fResiduals->UncheckedAt(kResidualPerDE_ClusterIn))->GetXaxis()->GetNbins();
+  fNDE = ((TH2F*)fResiduals->UncheckedAt(kResidualPerDEClusterIn))->GetXaxis()->GetNbins();
   
   for (Int_t ia = 0; ia < 2; ia++) {
     
@@ -611,51 +737,51 @@ void AliAnalysisTaskMuonResolution::Terminate(Option_t *)
     g->SetName(Form("gResidual%sPerChMean_ClusterIn",axes[ia]));
     g->SetTitle(Form("cluster-track residual-%s per Ch: mean (cluster in);chamber ID;<#Delta_{%s}> (cm)",axes[ia],axes[ia]));
     g->SetMarkerStyle(kFullDotLarge);
-    fChamberRes->AddAtAndExpand(g, kResidualPerChMean_ClusterIn+ia);
+    fChamberRes->AddAtAndExpand(g, kResidualPerChMeanClusterIn+ia);
     g = new TGraphErrors(AliMUONConstants::NTrackingCh());
     g->SetName(Form("gResidual%sPerChMean_ClusterOut",axes[ia]));
     g->SetTitle(Form("cluster-track residual-%s per Ch: mean (cluster out);chamber ID;<#Delta_{%s}> (cm)",axes[ia],axes[ia]));
     g->SetMarkerStyle(kFullDotLarge);
-    fChamberRes->AddAtAndExpand(g, kResidualPerChMean_ClusterOut+ia);
+    fChamberRes->AddAtAndExpand(g, kResidualPerChMeanClusterOut+ia);
     
     g = new TGraphErrors(2*AliMUONConstants::NTrackingCh());
     g->SetName(Form("gResidual%sPerHalfChMean_ClusterIn",axes[ia]));
     g->SetTitle(Form("cluster-track residual-%s per half Ch: mean (cluster in);half chamber ID;<#Delta_{%s}> (cm)",axes[ia],axes[ia]));
     g->SetMarkerStyle(kFullDotLarge);
-    fChamberRes->AddAtAndExpand(g, kResidualPerHalfChMean_ClusterIn+ia);
+    fChamberRes->AddAtAndExpand(g, kResidualPerHalfChMeanClusterIn+ia);
     g = new TGraphErrors(2*AliMUONConstants::NTrackingCh());
     g->SetName(Form("gResidual%sPerHalfChMean_ClusterOut",axes[ia]));
     g->SetTitle(Form("cluster-track residual-%s per half Ch: mean (cluster out);half chamber ID;<#Delta_{%s}> (cm)",axes[ia],axes[ia]));
     g->SetMarkerStyle(kFullDotLarge);
-    fChamberRes->AddAtAndExpand(g, kResidualPerHalfChMean_ClusterOut+ia);
+    fChamberRes->AddAtAndExpand(g, kResidualPerHalfChMeanClusterOut+ia);
     
     g = new TGraphErrors(fNDE);
     g->SetName(Form("gResidual%sPerDEMean_ClusterIn",axes[ia]));
     g->SetTitle(Form("cluster-track residual-%s per DE: mean (cluster in);DE ID;<#Delta_{%s}> (cm)",axes[ia],axes[ia]));
     g->SetMarkerStyle(kFullDotLarge);
-    fChamberRes->AddAtAndExpand(g, kResidualPerDEMean_ClusterIn+ia);
+    fChamberRes->AddAtAndExpand(g, kResidualPerDEMeanClusterIn+ia);
     g = new TGraphErrors(fNDE);
     g->SetName(Form("gResidual%sPerDEMean_ClusterOut",axes[ia]));
     g->SetTitle(Form("cluster-track residual-%s per DE: mean (cluster out);DE ID;<#Delta_{%s}> (cm)",axes[ia],axes[ia]));
     g->SetMarkerStyle(kFullDotLarge);
-    fChamberRes->AddAtAndExpand(g, kResidualPerDEMean_ClusterOut+ia);
+    fChamberRes->AddAtAndExpand(g, kResidualPerDEMeanClusterOut+ia);
     
     g = new TGraphErrors(AliMUONConstants::NTrackingCh());
     g->SetName(Form("gResidual%sPerChSigma_ClusterIn",axes[ia]));
     g->SetTitle(Form("cluster-track residual-%s per Ch: sigma (cluster in);chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
     g->SetMarkerStyle(kFullDotLarge);
-    fChamberRes->AddAtAndExpand(g, kResidualPerChSigma_ClusterIn+ia);
+    fChamberRes->AddAtAndExpand(g, kResidualPerChSigmaClusterIn+ia);
     g = new TGraphErrors(AliMUONConstants::NTrackingCh());
     g->SetName(Form("gResidual%sPerChSigma_ClusterOut",axes[ia]));
     g->SetTitle(Form("cluster-track residual-%s per Ch: sigma (cluster out);chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
     g->SetMarkerStyle(kFullDotLarge);
-    fChamberRes->AddAtAndExpand(g, kResidualPerChSigma_ClusterOut+ia);
+    fChamberRes->AddAtAndExpand(g, kResidualPerChSigmaClusterOut+ia);
     
     g = new TGraphErrors(AliMUONConstants::NTrackingCh());
     g->SetName(Form("gResidual%sPerChDispersion_ClusterOut",axes[ia]));
     g->SetTitle(Form("cluster-track residual-%s per Ch: dispersion (cluster out);chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
     g->SetMarkerStyle(kFullDotLarge);
-    fChamberRes->AddAtAndExpand(g, kResidualPerChDispersion_ClusterOut+ia);
+    fChamberRes->AddAtAndExpand(g, kResidualPerChDispersionClusterOut+ia);
     
     g = new TGraphErrors(AliMUONConstants::NTrackingCh());
     g->SetName(Form("gCombinedResidual%sPerChSigma",axes[ia]));
@@ -677,7 +803,7 @@ void AliAnalysisTaskMuonResolution::Terminate(Option_t *)
     
     mg = new TMultiGraph(Form("mgCombinedResidual%sSigmaVsP",axes[ia]),Form("cluster %s-resolution per chamber versus momentum;p (GeV/c^{2});#sigma_{%s} (cm)",axes[ia],axes[ia]));
     for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
-      g = new TGraphErrors(((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsP_ClusterIn+10*ia+i))->GetNbinsX());
+      g = new TGraphErrors(((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsPClusterIn+10*ia+i))->GetNbinsX());
       g->SetName(Form("gRes%sVsP_ch%d",axes[ia],i+1));
       g->SetMarkerStyle(kFullDotMedium);
       g->SetMarkerColor(i+1+i/9);
@@ -685,6 +811,28 @@ void AliAnalysisTaskMuonResolution::Terminate(Option_t *)
     }
     fChamberRes->AddAtAndExpand(mg, kCombinedResidualSigmaVsP+ia);
     
+    g = new TGraphErrors(((TH2F*)fResidualsVsP->UncheckedAt(kResidualVsPClusterIn+ia))->GetNbinsX());
+    g->SetName(Form("gCombinedResidual%sSigmaVsP",axes[ia]));
+    g->SetTitle(Form("cluster %s-resolution integrated over chambers versus momentum: sigma;p (GeV/c^{2});#sigma_{%s} (cm)",axes[ia],axes[ia]));
+    g->SetMarkerStyle(kFullDotLarge);
+    fChamberRes->AddAtAndExpand(g, kCombinedResidualAllChSigmaVsP+ia);
+    
+    mg = new TMultiGraph(Form("mgCombinedResidual%sSigmaVsCent",axes[ia]),Form("cluster %s-resolution per chamber versus centrality;centrality (%%);#sigma_{%s} (cm)",axes[ia],axes[ia]));
+    for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
+      g = new TGraphErrors(((TH2F*)fResidualsVsCent->UncheckedAt(kResidualInChVsCentClusterIn+10*ia+i))->GetNbinsX());
+      g->SetName(Form("gRes%sVsCent_ch%d",axes[ia],i+1));
+      g->SetMarkerStyle(kFullDotMedium);
+      g->SetMarkerColor(i+1+i/9);
+      mg->Add(g,"p");
+    }
+    fChamberRes->AddAtAndExpand(mg, kCombinedResidualSigmaVsCent+ia);
+    
+    g = new TGraphErrors(((TH2F*)fResidualsVsCent->UncheckedAt(kResidualVsCentClusterIn+ia))->GetNbinsX());
+    g->SetName(Form("gCombinedResidual%sSigmaVsCent",axes[ia]));
+    g->SetTitle(Form("cluster %s-resolution integrated over chambers versus centrality: sigma;centrality (%%);#sigma_{%s} (cm)",axes[ia],axes[ia]));
+    g->SetMarkerStyle(kFullDotLarge);
+    fChamberRes->AddAtAndExpand(g, kCombinedResidualAllChSigmaVsCent+ia);
+    
     g = new TGraphErrors(AliMUONConstants::NTrackingCh());
     g->SetName(Form("gTrackRes%sPerCh",axes[ia]));
     g->SetTitle(Form("track <#sigma_{%s}> per Ch;chamber ID;#sigma_{%s} (cm)",axes[ia],axes[ia]));
@@ -733,30 +881,42 @@ void AliAnalysisTaskMuonResolution::Terminate(Option_t *)
     g->SetMarkerStyle(kFullDotLarge);
     fLocalChi2->AddAtAndExpand(g, kLocalChi2PerDEMean+ia);
     
+    // compute residual mean and dispersion integrated over chambers versus p
+    FillSigmaClusterVsP((TH2F*)fResidualsVsP->UncheckedAt(kResidualVsPClusterIn+ia),
+                       (TH2F*)fResidualsVsP->UncheckedAt(kResidualVsPClusterOut+ia),
+                       (TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualAllChSigmaVsP+ia));
+    
+    // compute residual mean and dispersion integrated over chambers versus centrality
+    FillSigmaClusterVsCent((TH2F*)fResidualsVsCent->UncheckedAt(kResidualVsCentClusterIn+ia),
+                          (TH2F*)fResidualsVsCent->UncheckedAt(kResidualVsCentClusterOut+ia),
+                          (TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualAllChSigmaVsCent+ia));
+    
     // compute residual mean and dispersion and averaged local chi2 per chamber and half chamber
-    Double_t meanIn, meanInErr, meanOut, meanOutErr, sigmaIn, sigmaInErr, sigmaOut, sigmaOutErr;
+    Double_t meanIn, meanInErr, meanOut, meanOutErr, sigma, sigmaIn, sigmaInErr, sigmaOut, sigmaOutErr;
     Double_t sigmaTrack, sigmaTrackErr, sigmaMCS, sigmaMCSErr, clusterRes, clusterResErr, sigmaCluster, sigmaClusterErr;
     for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
       
       // method 1
-      TH1D *tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerCh_ClusterIn+ia))->ProjectionY("tmp",i+1,i+1,"e");
-      GetMean(tmp, meanIn, meanInErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChMean_ClusterIn+ia), i, i+1);
-      GetRMS(tmp, sigmaIn, sigmaInErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChSigma_ClusterIn+ia), i, i+1);
+      TH1D *tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerChClusterIn+ia))->ProjectionY("tmp",i+1,i+1,"e");
+      GetMean(tmp, meanIn, meanInErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChMeanClusterIn+ia), i, i+1);
+      GetRMS(tmp, sigmaIn, sigmaInErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChSigmaClusterIn+ia), i, i+1);
       delete tmp;
       
-      tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerCh_ClusterOut+ia))->ProjectionY("tmp",i+1,i+1,"e");
-      GetMean(tmp, meanOut, meanOutErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChMean_ClusterOut+ia), i, i+1);
-      GetRMS(tmp, sigmaOut, sigmaOutErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChSigma_ClusterOut+ia), i, i+1);
+      tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerChClusterOut+ia))->ProjectionY("tmp",i+1,i+1,"e");
+      GetMean(tmp, meanOut, meanOutErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChMeanClusterOut+ia), i, i+1);
+      GetRMS(tmp, sigmaOut, sigmaOutErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChSigmaClusterOut+ia), i, i+1);
       delete tmp;
       
       if (fCorrectForSystematics) {
-       sigmaIn = TMath::Sqrt(sigmaIn*sigmaIn + meanIn*meanIn);
-       sigmaInErr = (sigmaIn>0) ? TMath::Sqrt(sigmaIn*sigmaIn*sigmaInErr*sigmaInErr + meanIn*meanIn*meanInErr*meanInErr) / sigmaIn : 0.;
-       sigmaOut = TMath::Sqrt(sigmaOut*sigmaOut + meanOut*meanOut);
-       sigmaOutErr = (sigmaOut>0) ? TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + meanOut*meanOut*meanOutErr*meanOutErr) / sigmaOut : 0.;
+       sigma = TMath::Sqrt(sigmaIn*sigmaIn + meanIn*meanIn);
+       sigmaInErr = (sigma>0) ? TMath::Sqrt(sigmaIn*sigmaIn*sigmaInErr*sigmaInErr + meanIn*meanIn*meanInErr*meanInErr) / sigma : 0.;
+       sigmaIn = sigma;
+       sigma = TMath::Sqrt(sigmaOut*sigmaOut + meanOut*meanOut);
+       sigmaOutErr = (sigma>0) ? TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + meanOut*meanOut*meanOutErr*meanOutErr) / sigma : 0.;
+       sigmaOut = sigma;
       }
-      ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChDispersion_ClusterOut+ia))->SetPoint(i, i+1, sigmaOut);
-      ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChDispersion_ClusterOut+ia))->SetPointError(i, 0., sigmaOutErr);
+      ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChDispersionClusterOut+ia))->SetPoint(i, i+1, sigmaOut);
+      ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChDispersionClusterOut+ia))->SetPointError(i, 0., sigmaOutErr);
       
       clusterRes = TMath::Sqrt(sigmaIn*sigmaOut);
       //      clusterResErr = (clusterRes > 0.) ? 0.5 * TMath::Sqrt(sigmaInErr*sigmaInErr*sigmaOut*sigmaOut + sigmaIn*sigmaIn*sigmaOutErr*sigmaOutErr) / clusterRes : 0.;
@@ -768,11 +928,11 @@ void AliAnalysisTaskMuonResolution::Terminate(Option_t *)
       
       // method 2
       tmp = ((TH2F*)fResiduals->UncheckedAt(kTrackResPerCh+ia))->ProjectionY("tmp",i+1,i+1,"e");
-      GetMean(tmp, sigmaTrack, sigmaTrackErr, (TGraphErrors*)fChamberRes->UncheckedAt(kTrackResPerChMean+ia), i, i+1, kFALSE);
+      GetMean(tmp, sigmaTrack, sigmaTrackErr, (TGraphErrors*)fChamberRes->UncheckedAt(kTrackResPerChMean+ia), i, i+1, kFALSE, kFALSE);
       delete tmp;
       
       tmp = ((TH2F*)fResiduals->UncheckedAt(kMCSPerCh+ia))->ProjectionY("tmp",i+1,i+1,"e");
-      GetMean(tmp, sigmaMCS, sigmaMCSErr, (TGraphErrors*)fChamberRes->UncheckedAt(kMCSPerChMean+ia), i, i+1, kFALSE);
+      GetMean(tmp, sigmaMCS, sigmaMCSErr, (TGraphErrors*)fChamberRes->UncheckedAt(kMCSPerChMean+ia), i, i+1, kFALSE, kFALSE);
       delete tmp;
       
       sigmaCluster = sigmaOut*sigmaOut - sigmaTrack*sigmaTrack;
@@ -800,30 +960,37 @@ void AliAnalysisTaskMuonResolution::Terminate(Option_t *)
       delete tmp;
       
       // method 1 versus p
-      FillSigmaClusterVsP((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsP_ClusterIn+10*ia+i),
-                         (TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsP_ClusterOut+10*ia+i),
+      FillSigmaClusterVsP((TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsPClusterIn+10*ia+i),
+                         (TH2F*)fResidualsVsP->UncheckedAt(kResidualInChVsPClusterOut+10*ia+i),
                          (TGraphErrors*)((TMultiGraph*)fChamberRes->UncheckedAt(kCombinedResidualSigmaVsP+ia))->GetListOfGraphs()->FindObject(Form("gRes%sVsP_ch%d",axes[ia],i+1)));
       
+      // method 1 versus centrality
+      FillSigmaClusterVsCent((TH2F*)fResidualsVsCent->UncheckedAt(kResidualInChVsCentClusterIn+10*ia+i),
+                            (TH2F*)fResidualsVsCent->UncheckedAt(kResidualInChVsCentClusterOut+10*ia+i),
+                            (TGraphErrors*)((TMultiGraph*)fChamberRes->UncheckedAt(kCombinedResidualSigmaVsCent+ia))->GetListOfGraphs()->FindObject(Form("gRes%sVsCent_ch%d",axes[ia],i+1)));
+      
       // compute residual mean and dispersion per half chamber
       for (Int_t j = 0; j < 2; j++) {
        Int_t k = 2*i+j;
        
        // method 1
-       tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfCh_ClusterIn+ia))->ProjectionY("tmp",k+1,k+1,"e");
-       GetMean(tmp, meanIn, meanInErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMean_ClusterIn+ia), k, k+1);
+       tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfChClusterIn+ia))->ProjectionY("tmp",k+1,k+1,"e");
+       GetMean(tmp, meanIn, meanInErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMeanClusterIn+ia), k, k+1);
        GetRMS(tmp, sigmaIn, sigmaInErr);
        delete tmp;
        
-       tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfCh_ClusterOut+ia))->ProjectionY("tmp",k+1,k+1,"e");
-       GetMean(tmp, meanOut, meanOutErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMean_ClusterOut+ia), k, k+1);
+       tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfChClusterOut+ia))->ProjectionY("tmp",k+1,k+1,"e");
+       GetMean(tmp, meanOut, meanOutErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMeanClusterOut+ia), k, k+1);
        GetRMS(tmp, sigmaOut, sigmaOutErr);
        delete tmp;
        
        if (fCorrectForSystematics) {
-         sigmaIn = TMath::Sqrt(sigmaIn*sigmaIn + meanIn*meanIn);
-         sigmaInErr = (sigmaIn>0) ? TMath::Sqrt(sigmaIn*sigmaIn*sigmaInErr*sigmaInErr + meanIn*meanIn*meanInErr*meanInErr) / sigmaIn : 0.;
-         sigmaOut = TMath::Sqrt(sigmaOut*sigmaOut + meanOut*meanOut);
-         sigmaOutErr = (sigmaOut>0) ? TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + meanOut*meanOut*meanOutErr*meanOutErr) / sigmaOut : 0.;
+         sigma = TMath::Sqrt(sigmaIn*sigmaIn + meanIn*meanIn);
+         sigmaInErr = (sigma>0) ? TMath::Sqrt(sigmaIn*sigmaIn*sigmaInErr*sigmaInErr + meanIn*meanIn*meanInErr*meanInErr) / sigma : 0.;
+         sigmaIn = sigma;
+         sigma = TMath::Sqrt(sigmaOut*sigmaOut + meanOut*meanOut);
+         sigmaOutErr = (sigma>0) ? TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + meanOut*meanOut*meanOutErr*meanOutErr) / sigma : 0.;
+         sigmaOut = sigma;
        }
        
        clusterRes = TMath::Sqrt(sigmaIn*sigmaOut);
@@ -834,11 +1001,11 @@ void AliAnalysisTaskMuonResolution::Terminate(Option_t *)
        
        // method 2
        tmp = ((TH2F*)fResiduals->UncheckedAt(kTrackResPerHalfCh+ia))->ProjectionY("tmp",k+1,k+1,"e");
-       GetMean(tmp, sigmaTrack, sigmaTrackErr, 0x0, 0, 0, kFALSE);
+       GetMean(tmp, sigmaTrack, sigmaTrackErr, 0x0, 0, 0, kFALSE, kFALSE);
        delete tmp;
        
        tmp = ((TH2F*)fResiduals->UncheckedAt(kMCSPerHalfCh+ia))->ProjectionY("tmp",k+1,k+1,"e");
-       GetMean(tmp, sigmaMCS, sigmaMCSErr, 0x0, 0, 0, kFALSE);
+       GetMean(tmp, sigmaMCS, sigmaMCSErr, 0x0, 0, 0, kFALSE, kFALSE);
        delete tmp;
        
        sigmaCluster = sigmaOut*sigmaOut - sigmaTrack*sigmaTrack;
@@ -866,21 +1033,23 @@ void AliAnalysisTaskMuonResolution::Terminate(Option_t *)
     for (Int_t i = 0; i < fNDE; i++) {
       
       // method 1
-      TH1D *tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerDE_ClusterIn+ia))->ProjectionY("tmp",i+1,i+1,"e");
-      GetMean(tmp, meanIn, meanInErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMean_ClusterIn+ia), i, i+1);
+      TH1D *tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerDEClusterIn+ia))->ProjectionY("tmp",i+1,i+1,"e");
+      GetMean(tmp, meanIn, meanInErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMeanClusterIn+ia), i, i+1);
       GetRMS(tmp, sigmaIn, sigmaInErr);
       delete tmp;
       
-      tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerDE_ClusterOut+ia))->ProjectionY("tmp",i+1,i+1,"e");
-      GetMean(tmp, meanOut, meanOutErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMean_ClusterOut+ia), i, i+1);
+      tmp = ((TH2F*)fResiduals->UncheckedAt(kResidualPerDEClusterOut+ia))->ProjectionY("tmp",i+1,i+1,"e");
+      GetMean(tmp, meanOut, meanOutErr, (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMeanClusterOut+ia), i, i+1);
       GetRMS(tmp, sigmaOut, sigmaOutErr);
       delete tmp;
       
       if (fCorrectForSystematics) {
-       sigmaIn = TMath::Sqrt(sigmaIn*sigmaIn + meanIn*meanIn);
-       sigmaInErr = (sigmaIn>0) ? TMath::Sqrt(sigmaIn*sigmaIn*sigmaInErr*sigmaInErr + meanIn*meanIn*meanInErr*meanInErr) / sigmaIn : 0.;
-       sigmaOut = TMath::Sqrt(sigmaOut*sigmaOut + meanOut*meanOut);
-       sigmaOutErr = (sigmaOut>0) ? TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + meanOut*meanOut*meanOutErr*meanOutErr) / sigmaOut : 0.;
+       sigma = TMath::Sqrt(sigmaIn*sigmaIn + meanIn*meanIn);
+       sigmaInErr = (sigma>0) ? TMath::Sqrt(sigmaIn*sigmaIn*sigmaInErr*sigmaInErr + meanIn*meanIn*meanInErr*meanInErr) / sigma : 0.;
+       sigmaIn = sigma;
+       sigma = TMath::Sqrt(sigmaOut*sigmaOut + meanOut*meanOut);
+       sigmaOutErr = (sigma>0) ? TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + meanOut*meanOut*meanOutErr*meanOutErr) / sigma : 0.;
+       sigmaOut = sigma;
       }
       
       clusterRes = TMath::Sqrt(sigmaIn*sigmaOut);
@@ -891,11 +1060,11 @@ void AliAnalysisTaskMuonResolution::Terminate(Option_t *)
       
       // method 2
       tmp = ((TH2F*)fResiduals->UncheckedAt(kTrackResPerDE+ia))->ProjectionY("tmp",i+1,i+1,"e");
-      GetMean(tmp, sigmaTrack, sigmaTrackErr, 0x0, 0, 0, kFALSE);
+      GetMean(tmp, sigmaTrack, sigmaTrackErr, 0x0, 0, 0, kFALSE, kFALSE);
       delete tmp;
       
       tmp = ((TH2F*)fResiduals->UncheckedAt(kMCSPerDE+ia))->ProjectionY("tmp",i+1,i+1,"e");
-      GetMean(tmp, sigmaMCS, sigmaMCSErr, 0x0, 0, 0, kFALSE);
+      GetMean(tmp, sigmaMCS, sigmaMCSErr, 0x0, 0, 0, kFALSE, kFALSE);
       delete tmp;
       
       sigmaCluster = sigmaOut*sigmaOut - sigmaTrack*sigmaTrack;
@@ -917,17 +1086,31 @@ void AliAnalysisTaskMuonResolution::Terminate(Option_t *)
       
     }
     
-    // set graph labels
-    TAxis* xAxis = ((TH2F*)fResiduals->UncheckedAt(kResidualPerDE_ClusterOut+ia))->GetXaxis();
-    ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMean_ClusterIn+ia))->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
-    ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMean_ClusterOut+ia))->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
+    // set half-chamber graph labels
+    TAxis* xAxis = ((TH2F*)fResiduals->UncheckedAt(kResidualPerHalfChClusterIn+ia))->GetXaxis();
+    ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMeanClusterIn+ia))->GetXaxis()->Set(20, 0.5, 20.5);
+    ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMeanClusterOut+ia))->GetXaxis()->Set(20, 0.5, 20.5);
+    ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerHalfChSigma+ia))->GetXaxis()->Set(20, 0.5, 20.5);
+    ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerHalfCh+ia))->GetXaxis()->Set(20, 0.5, 20.5);
+    for (Int_t i = 1; i <= 20; i++) {
+      const char* label = xAxis->GetBinLabel(i);
+      ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMeanClusterIn+ia))->GetXaxis()->SetBinLabel(i, label);
+      ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMeanClusterOut+ia))->GetXaxis()->SetBinLabel(i, label);
+      ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerHalfChSigma+ia))->GetXaxis()->SetBinLabel(i, label);
+      ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerHalfCh+ia))->GetXaxis()->SetBinLabel(i, label);
+    }
+    
+    // set DE graph labels
+    xAxis = ((TH2F*)fResiduals->UncheckedAt(kResidualPerDEClusterOut+ia))->GetXaxis();
+    ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMeanClusterIn+ia))->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
+    ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMeanClusterOut+ia))->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
     ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerDESigma+ia))->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
     ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerDE+ia))->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
     ((TGraphErrors*)fLocalChi2->UncheckedAt(kLocalChi2PerDEMean+ia))->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
     for (Int_t i = 1; i <= fNDE; i++) {
       const char* label = xAxis->GetBinLabel(i);
-      ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMean_ClusterIn+ia))->GetXaxis()->SetBinLabel(i, label);
-      ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMean_ClusterOut+ia))->GetXaxis()->SetBinLabel(i, label);
+      ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMeanClusterIn+ia))->GetXaxis()->SetBinLabel(i, label);
+      ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMeanClusterOut+ia))->GetXaxis()->SetBinLabel(i, label);
       ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerDESigma+ia))->GetXaxis()->SetBinLabel(i, label);
       ((TGraphErrors*)fChamberRes->UncheckedAt(kClusterResPerDE+ia))->GetXaxis()->SetBinLabel(i, label);
       ((TGraphErrors*)fLocalChi2->UncheckedAt(kLocalChi2PerDEMean+ia))->GetXaxis()->SetBinLabel(i, label);
@@ -983,12 +1166,12 @@ void AliAnalysisTaskMuonResolution::Terminate(Option_t *)
   cResPerCh->Divide(4,2);
   for (Int_t ia = 0; ia < 2; ia++) {
     cResPerCh->cd(1+4*ia);
-    g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChMean_ClusterOut+ia);
+    g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChMeanClusterOut+ia);
     g->Draw("ap");
     g->SetMarkerColor(2);
     g->SetLineColor(2);
     if (ia == 0) lResPerChMean->AddEntry(g,"cluster out","PL");
-    g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChMean_ClusterIn+ia);
+    g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChMeanClusterIn+ia);
     g->Draw("p");
     g->SetMarkerColor(4);
     g->SetLineColor(4);
@@ -996,13 +1179,13 @@ void AliAnalysisTaskMuonResolution::Terminate(Option_t *)
     if (ia == 0) lResPerChMean->Draw();
     else lResPerChMean->DrawClone();
     cResPerCh->cd(2+4*ia);
-    g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChSigma_ClusterOut+ia);
+    g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChSigmaClusterOut+ia);
     g->Draw("ap");
     g->SetMinimum(0.);
     g->SetMarkerColor(2);
     g->SetLineColor(2);
     if (ia == 0) lResPerChSigma1->AddEntry(g,"cluster out","PL");
-    g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChSigma_ClusterIn+ia);
+    g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChSigmaClusterIn+ia);
     g->Draw("p");
     g->SetMarkerColor(4);
     g->SetLineColor(4);
@@ -1020,7 +1203,7 @@ void AliAnalysisTaskMuonResolution::Terminate(Option_t *)
     if (ia == 0) lResPerChSigma1->Draw();
     else lResPerChSigma1->DrawClone();
     cResPerCh->cd(3+4*ia);
-    g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChDispersion_ClusterOut+ia);
+    g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChDispersionClusterOut+ia);
     g->Draw("ap");
     g->SetMinimum(0.);
     g->SetMarkerColor(2);
@@ -1056,11 +1239,11 @@ void AliAnalysisTaskMuonResolution::Terminate(Option_t *)
   cResPerHalfCh->Divide(2,2);
   for (Int_t ia = 0; ia < 2; ia++) {
     cResPerHalfCh->cd(1+2*ia);
-    g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMean_ClusterOut+ia);
+    g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMeanClusterOut+ia);
     g->Draw("ap");
     g->SetMarkerColor(2);
     g->SetLineColor(2);
-    g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMean_ClusterIn+ia);
+    g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerHalfChMeanClusterIn+ia);
     g->Draw("p");
     g->SetMarkerColor(4);
     g->SetLineColor(4);
@@ -1081,11 +1264,11 @@ void AliAnalysisTaskMuonResolution::Terminate(Option_t *)
   cResPerDE->Divide(1,4);
   for (Int_t ia = 0; ia < 2; ia++) {
     cResPerDE->cd(1+ia);
-    g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMean_ClusterOut+ia);
+    g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMeanClusterOut+ia);
     g->Draw("ap");
     g->SetMarkerColor(2);
     g->SetLineColor(2);
-    g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMean_ClusterIn+ia);
+    g = (TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerDEMeanClusterIn+ia);
     g->Draw("p");
     g->SetMarkerColor(4);
     g->SetLineColor(4);
@@ -1111,17 +1294,44 @@ void AliAnalysisTaskMuonResolution::Terminate(Option_t *)
   }
   fCanvases->AddAtAndExpand(cResPerChVsP, kResPerChVsP);
   
+  TCanvas* cResPerChVsCent = new TCanvas("cResPerChVsCent","cResPerChVsCent");
+  cResPerChVsCent->Divide(1,2);
+  for (Int_t ia = 0; ia < 2; ia++) {
+    cResPerChVsCent->cd(1+ia);
+    mg = (TMultiGraph*)fChamberRes->UncheckedAt(kCombinedResidualSigmaVsCent+ia);
+    mg->Draw("ap");
+  }
+  fCanvases->AddAtAndExpand(cResPerChVsCent, kResPerChVsCent);
+  
   // print results
-  printf("\nchamber resolution:\n");
-  printf(" - non-bending:");
-  for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %5.3f":", %5.3f",newClusterRes[0][i]);
-  printf("\n -     bending:");
-  for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %6.4f":", %6.4f",newClusterRes[1][i]);
-  printf("\n\n");
+  if (fPrintClResPerCh) {
+    printf("\nchamber resolution:\n");
+    printf(" - non-bending:");
+    for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %5.3f":", %5.3f",newClusterRes[0][i]);
+    printf("\n -     bending:");
+    for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %6.4f":", %6.4f",newClusterRes[1][i]);
+    printf("\n\n");
+  }
+  
+  if (fPrintClResPerDE) {
+    Double_t iDE, clRes;
+    printf("\nDE resolution:\n");
+    printf(" - non-bending:");
+    for (Int_t i = 0; i < fNDE; i++) {
+      ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerDESigma))->GetPoint(i, iDE, clRes);
+      printf((i==0)?" %5.3f":", %5.3f", clRes);
+    }
+    printf("\n -     bending:");
+    for (Int_t i = 0; i < fNDE; i++) {
+      ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerDESigma+1))->GetPoint(i, iDE, clRes);
+      printf((i==0)?" %6.4f":", %6.4f", clRes);
+    }
+    printf("\n\n");
+  }
   
   // Post final data.
-  PostData(3, fLocalChi2);
-  PostData(4, fChamberRes);
+  PostData(4, fLocalChi2);
+  PostData(5, fChamberRes);
 }
 
 //________________________________________________________________________
@@ -1202,39 +1412,115 @@ void AliAnalysisTaskMuonResolution::ZoomRight(TH1* h, Double_t fractionCut)
 }
 
 //________________________________________________________________________
-void AliAnalysisTaskMuonResolution::GetMean(TH1* h, Double_t& mean, Double_t& meanErr, TGraphErrors* g, Int_t i, Double_t x, Bool_t zoom)
+void AliAnalysisTaskMuonResolution::GetMean(TH1* h, Double_t& mean, Double_t& meanErr, TGraphErrors* g, Int_t i, Double_t x, Bool_t zoom, Bool_t enableFit)
 {
-  /// Fill graph with the mean value of the histogram and the corresponding error (zooming if required)
-  Int_t firstBin = h->GetXaxis()->GetFirst();
-  Int_t lastBin = h->GetXaxis()->GetLast();
-  if (zoom) Zoom(h);
-  mean = (h->GetEntries() > fgkMinEntries) ? h->GetMean() : 0.;
-  meanErr = (h->GetEntries() > fgkMinEntries) ? h->GetMeanError() : 0.;
+  /// Fill graph with the mean value and the corresponding error (zooming if required)
+  
+  if (h->GetEntries() < fgkMinEntries) { // not enough entries
+    
+    mean = 0.;
+    meanErr = 0.;
+    
+  } else if (enableFit && fGaus) { // take the mean of a gaussian fit
+    
+    fGaus->SetParameters(h->GetEntries(), 0., 0.1);
+    
+    if (h->GetUniqueID() != 10) {
+      
+      // first fit
+      h->Fit("fGaus", "WWNQ");
+      
+      // rebin histo
+      Int_t rebin = TMath::Max(static_cast<Int_t>(0.2*fGaus->GetParameter(2)/h->GetBinWidth(1)),1);
+      while (h->GetNbinsX()%rebin!=0) rebin--;
+      h->Rebin(rebin);
+      
+      // use the unique ID to remember that this histogram has already been rebinned
+      h->SetUniqueID(10);
+      
+    }
+    
+    // second fit
+    h->Fit("fGaus","NQ");
+    
+    mean = fGaus->GetParameter(1);
+    meanErr = fGaus->GetParError(1);
+    
+  } else { // take the mean of the distribution
+    
+    if (zoom) Zoom(h);
+    
+    mean = h->GetMean();
+    meanErr = h->GetMeanError();
+    
+    if (zoom) h->GetXaxis()->SetRange(0,0);
+    
+  }
+  
+  // fill graph if required
   if (g) {
     g->SetPoint(i, x, mean);
     g->SetPointError(i, 0., meanErr);
   }
-  if (zoom) h->GetXaxis()->SetRange(firstBin,lastBin);
+  
 }
 
 //________________________________________________________________________
 void AliAnalysisTaskMuonResolution::GetRMS(TH1* h, Double_t& rms, Double_t& rmsErr, TGraphErrors* g, Int_t i, Double_t x, Bool_t zoom)
 {
-  /// Return the RMS of the histogram and the corresponding error (zooming if required) and fill graph if !=0x0
-  Int_t firstBin = h->GetXaxis()->GetFirst();
-  Int_t lastBin = h->GetXaxis()->GetLast();
-  if (zoom) Zoom(h);
-  rms = (h->GetEntries() > fgkMinEntries) ? h->GetRMS() : 0.;
-  rmsErr = (h->GetEntries() > fgkMinEntries) ? h->GetRMSError() : 0.;
+  /// Return the dispersion value and the corresponding error (zooming if required) and fill graph if !=0x0
+  
+  if (h->GetEntries() < fgkMinEntries) { // not enough entries
+    
+    rms = 0.;
+    rmsErr = 0.;
+    
+  } else if (fGaus) { // take the sigma of a gaussian fit
+    
+    fGaus->SetParameters(h->GetEntries(), 0., 0.1);
+    
+    if (h->GetUniqueID() != 10) {
+      
+      // first fit
+      h->Fit("fGaus", "WWNQ");
+      
+      // rebin histo
+      Int_t rebin = TMath::Max(static_cast<Int_t>(0.2*fGaus->GetParameter(2)/h->GetBinWidth(1)),1);
+      while (h->GetNbinsX()%rebin!=0) rebin--;
+      h->Rebin(rebin);
+      
+      // use the unique ID to remember that this histogram has already been rebinned
+      h->SetUniqueID(10);
+      
+    }
+    
+    // second fit
+    h->Fit("fGaus","NQ");
+    
+    rms = fGaus->GetParameter(2);
+    rmsErr = fGaus->GetParError(2);
+    
+  } else { // take the RMS of the distribution
+    
+    if (zoom) Zoom(h);
+    
+    rms = h->GetRMS();
+    rmsErr = h->GetRMSError();
+    
+    if (zoom) h->GetXaxis()->SetRange(0,0);
+    
+  }
+  
+  // fill graph if required
   if (g) {
     g->SetPoint(i, x, rms);
     g->SetPointError(i, 0., rmsErr);
   }
-  if (zoom) h->GetXaxis()->SetRange(firstBin,lastBin);
+  
 }
 
 //________________________________________________________________________
-void AliAnalysisTaskMuonResolution::FillSigmaClusterVsP(TH2* hIn, TH2* hOut, TGraphErrors* g, Bool_t zoom)
+void AliAnalysisTaskMuonResolution::FillSigmaClusterVsP(const TH2* hIn, const TH2* hOut, TGraphErrors* g, Bool_t zoom)
 {
   /// Fill graph with cluster resolution from combined residuals with cluster in/out (zooming if required)
   Double_t sigmaIn, sigmaInErr, sigmaOut, sigmaOutErr, clusterRes, clusterResErr;
@@ -1248,12 +1534,35 @@ void AliAnalysisTaskMuonResolution::FillSigmaClusterVsP(TH2* hIn, TH2* hOut, TGr
     Double_t p = 0.5 * (hIn->GetBinLowEdge(j) + hIn->GetBinLowEdge(j+1));
     Double_t pErr = p - hIn->GetBinLowEdge(j);
     clusterRes = TMath::Sqrt(sigmaIn*sigmaOut);
-    clusterResErr = (clusterRes > 0.) ? 0.5 * TMath::Sqrt(sigmaInErr*sigmaInErr*sigmaOut*sigmaOut + sigmaIn*sigmaIn*sigmaOutErr*sigmaOutErr) / clusterRes : 0.;
+    //clusterResErr = (clusterRes > 0.) ? 0.5 * TMath::Sqrt(sigmaInErr*sigmaInErr*sigmaOut*sigmaOut + sigmaIn*sigmaIn*sigmaOutErr*sigmaOutErr) / clusterRes : 0.;
+    clusterResErr = TMath::Sqrt(sigmaInErr*sigmaOutErr);
     g->SetPoint(j, p, clusterRes);
     g->SetPointError(j, pErr, clusterResErr);
   }
 }
 
+//________________________________________________________________________
+void AliAnalysisTaskMuonResolution::FillSigmaClusterVsCent(const TH2* hIn, const TH2* hOut, TGraphErrors* g, Bool_t zoom)
+{
+  /// Fill graph with cluster resolution from combined residuals with cluster in/out (zooming if required)
+  Double_t sigmaIn, sigmaInErr, sigmaOut, sigmaOutErr, clusterRes, clusterResErr;
+  for (Int_t j = 1; j <= hIn->GetNbinsX(); j++) {
+    TH1D* tmp = hIn->ProjectionY("tmp",j,j,"e");
+    GetRMS(tmp, sigmaIn, sigmaInErr, 0x0, 0, 0., zoom);
+    delete tmp;
+    tmp = hOut->ProjectionY("tmp",j,j,"e");
+    GetRMS(tmp, sigmaOut, sigmaOutErr, 0x0, 0, 0., zoom);
+    delete tmp;
+    Double_t cent = 0.5 * (hIn->GetBinLowEdge(j) + hIn->GetBinLowEdge(j+1));
+    Double_t centErr = cent - hIn->GetBinLowEdge(j);
+    clusterRes = TMath::Sqrt(sigmaIn*sigmaOut);
+    //clusterResErr = (clusterRes > 0.) ? 0.5 * TMath::Sqrt(sigmaInErr*sigmaInErr*sigmaOut*sigmaOut + sigmaIn*sigmaIn*sigmaOutErr*sigmaOutErr) / clusterRes : 0.;
+    clusterResErr = TMath::Sqrt(sigmaInErr*sigmaOutErr);
+    g->SetPoint(j, cent, clusterRes);
+    g->SetPointError(j, centErr, clusterResErr);
+  }
+}
+
 //__________________________________________________________________________
 void AliAnalysisTaskMuonResolution::Cov2CovP(const AliMUONTrackParam &param, TMatrixD &covP)
 {
@@ -1286,3 +1595,103 @@ void AliAnalysisTaskMuonResolution::Cov2CovP(const AliMUONTrackParam &param, TMa
   covP.Mult(jacob,tmp);
 }
 
+//__________________________________________________________________________
+UInt_t AliAnalysisTaskMuonResolution::BuildTriggerWord(const TString& firedTriggerClasses)
+{
+  /// build the trigger word from the fired trigger classes and the list of selectable trigger
+  
+  UInt_t word = 0;
+  
+  TObjString* trigClasseName = 0x0;
+  TIter nextTrigger(fSelectTriggerClass);
+  while ((trigClasseName = static_cast<TObjString*>(nextTrigger()))) {
+    
+    TRegexp GenericTriggerClasseName(trigClasseName->String());
+    if (firedTriggerClasses.Contains(GenericTriggerClasseName)) word |= trigClasseName->GetUniqueID();
+    
+  }
+  
+  return word;
+}
+
+//__________________________________________________________________________
+void AliAnalysisTaskMuonResolution::CheckPads(AliMUONVCluster *cl, Bool_t &hasBending, Bool_t &hasNonBending) const
+{
+  /// Check that this cluster contains pads on both cathods
+  
+  // reset
+  hasBending = kFALSE;
+  hasNonBending = kFALSE;
+  
+  // loop over digits contained in the cluster
+  for (Int_t iDigit = 0; iDigit < cl->GetNDigits(); iDigit++) {
+    
+    Int_t manuId = AliMUONVDigit::ManuId(cl->GetDigitId(iDigit));
+    
+    // check the location of the manu the digit belongs to
+    if (manuId > 0) {
+      if (manuId & AliMpConstants::ManuMask(AliMp::kNonBendingPlane)) hasNonBending = kTRUE;
+      else hasBending = kTRUE;
+    }
+    
+  }
+  
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskMuonResolution::CheckPadsBelow(AliMUONVCluster *cl, Bool_t &hasBending, Bool_t &hasNonBending) const
+{
+  /// Check that this cluster contains pads on both cathods just under its position
+  
+  // reset
+  hasBending = kFALSE;
+  hasNonBending = kFALSE;
+  
+  // get the cathod corresponding to the bending/non-bending plane
+  Int_t deId = cl->GetDetElemId();
+  AliMpDetElement* de = AliMpDDLStore::Instance()->GetDetElement(deId, kFALSE);
+  if (!de) return;
+  AliMp::CathodType cath1 = de->GetCathodType(AliMp::kBendingPlane); 
+  AliMp::CathodType cath2 = de->GetCathodType(AliMp::kNonBendingPlane); 
+  
+  // get the corresponding segmentation
+  const AliMpVSegmentation* seg1 = AliMpSegmentation::Instance()->GetMpSegmentation(deId, cath1);
+  const AliMpVSegmentation* seg2 = AliMpSegmentation::Instance()->GetMpSegmentation(deId, cath2);
+  if (!seg1 || !seg2) return;
+  
+  // get local coordinate of the cluster
+  Double_t lX,lY,lZ;
+  Double_t gX = cl->GetX();
+  Double_t gY = cl->GetY();
+  Double_t gZ = cl->GetZ();
+  fNewGeoTransformer->Global2Local(deId,gX,gY,gZ,lX,lY,lZ);
+  
+  // find pads below the cluster
+  AliMpPad pad1 = seg1->PadByPosition(lX, lY, kFALSE);
+  AliMpPad pad2 = seg2->PadByPosition(lX, lY, kFALSE);
+  
+  // build their ID if pads are valid
+  UInt_t padId1 = (pad1.IsValid()) ? AliMUONVDigit::BuildUniqueID(deId, pad1.GetManuId(), pad1.GetManuChannel(), cath1) : 0;
+  UInt_t padId2 = (pad2.IsValid()) ? AliMUONVDigit::BuildUniqueID(deId, pad2.GetManuId(), pad2.GetManuChannel(), cath2) : 0;
+  
+  // check if the cluster contains these pads 
+  for (Int_t iDigit = 0; iDigit < cl->GetNDigits(); iDigit++) {
+    
+    UInt_t digitId = cl->GetDigitId(iDigit);
+    
+    if (digitId == padId1) {
+      
+      hasBending = kTRUE;
+      if (hasNonBending) break;
+      
+    } else if (digitId == padId2) {
+      
+      hasNonBending = kTRUE;
+      if (hasBending) break;
+      
+    }
+    
+  }
+  
+}
+