From 089b43581f81d07481955e138fb15dcd60dd0f20 Mon Sep 17 00:00:00 2001
From: Gaetano Barone <gaetano.barone@cern.ch>
Date: Tue, 27 Jul 2021 17:53:12 +0200
Subject: [PATCH 01/14] Adding systematic changes for single hist method

---
 .../MuonMomentumCorrections/CMakeLists.txt    |   1 +
 .../MuonCalibrationAndSmearingTool.h          |   5 +-
 .../Root/MuonCalibrationAndSmearingTool.cxx   | 180 ++++++++++++++----
 3 files changed, 152 insertions(+), 34 deletions(-)

diff --git a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/CMakeLists.txt b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/CMakeLists.txt
index 993bc7d99f23..9d55fcedfcf0 100644
--- a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/CMakeLists.txt
+++ b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/CMakeLists.txt
@@ -66,3 +66,4 @@ endif()
 
 # Install files from the package:
 atlas_install_joboptions( share/*.py )
+atlas_install_data( data/* )
diff --git a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/MuonMomentumCorrections/MuonCalibrationAndSmearingTool.h b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/MuonMomentumCorrections/MuonCalibrationAndSmearingTool.h
index dae845658046..fe866b406f4c 100644
--- a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/MuonMomentumCorrections/MuonCalibrationAndSmearingTool.h
+++ b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/MuonMomentumCorrections/MuonCalibrationAndSmearingTool.h
@@ -42,7 +42,7 @@ namespace MCAST {
   namespace DetectorType { enum { MS = 1, ID = 2, CB = 3 }; }
   namespace SystVariation { enum { Default = 0, Down = -1, Up = 1 }; }
   namespace SagittaCorType { enum { CB=0, ID=1, ME=2, WEIGHTS=3, AUTO=4}; }
-  namespace SagittaSysType { enum { NOMINAL=0, RHO=1, BIAS=2}; }
+  namespace SagittaSysType { enum { NOMINAL=0, RHO=1, BIAS=2, DATASTAT=3}; }
   namespace MST_Categories { enum { Undefined = -1, Zero = 0, One = 1, Two = 2, Three = 3, Four = 4, Total = 5 }; }
   namespace SagittaInputHistType { enum {NOMINAL=0,SINGLE=1 };   } 
 }
@@ -130,7 +130,7 @@ class MuonCalibrationAndSmearingTool : public virtual IMuonCalibrationAndSmearin
     virtual CorrectionCode applyStatCombination( xAOD::Muon& mu, InfoHelper& muonInfo ) const;
     virtual CorrectionCode applySagittaBiasCorrectionAuto(const int DetType, xAOD::Muon& mu, bool isMC, const unsigned int SystCase, InfoHelper& muonInfo) const;
     virtual CorrectionCode CorrectForCharge(double p2, double& pt, int q, bool isMC, double p2Kin=0) const;
-    virtual CorrectionCode applySagittaBiasCorrection(const unsigned int SgCorrType, xAOD::Muon& mu, unsigned int iter, bool stop, bool isMC, InfoHelper& muonInfo) const;
+  virtual CorrectionCode applySagittaBiasCorrection(const unsigned int SgCorrType, xAOD::Muon& mu, unsigned int iter, bool stop, bool isMC, InfoHelper& muonInfo, const unsigned int SystCase=0) const;
 
 
   protected:
@@ -181,6 +181,7 @@ class MuonCalibrationAndSmearingTool : public virtual IMuonCalibrationAndSmearin
       double ScaleMS_egLoss;
       double SagittaRho;
       double SagittaBias;
+      double SagittaDataStat;
     };
 
     bool  m_useExternalSeed;
diff --git a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx
index 8ec9206cb4cc..5479cdf02489 100644
--- a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx
+++ b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx
@@ -359,7 +359,7 @@ namespace CP {
 	      ATH_MSG_INFO("Setting up the tool for sinlge histogram input for sagitta bias corrections");
 	      m_SagittaIterations.clear();
         m_SagittaIterations={1,1,1};
-        m_SagittaCorrPhaseSpace=false; 
+        m_SagittaCorrPhaseSpace=true; 
       }
 
       std::vector<std::string> trackNames; trackNames.push_back("CB"); trackNames.push_back("ID");  trackNames.push_back("ME");
@@ -375,7 +375,13 @@ namespace CP {
 	    }}
 	  else if (m_saggitaMapsInputType==MCAST::SagittaInputHistType::SINGLE){
 	    ATH_MSG_VERBOSE("Track "<<i<<" file "<< PathResolverFindCalibFile(Form("MuonMomentumCorrections/%s/%s/",m_SagittaRelease.c_str(),m_year.c_str())+ trackNames.at(i) + "_data.root"));
+
 	    MuonCalibrationAndSmearingTool::setSagittaHistogramsSingle(GetHistSingleMethod(PathResolverFindCalibFile(Form("MuonMomentumCorrections/%s/%s/",m_SagittaRelease.c_str(),m_year.c_str())+ trackNames.at(i) + "_data.root"),Form("p%s_0",trackNames.at(i).c_str())),i);
+
+	    MuonCalibrationAndSmearingTool::setSagittaHistogramsSingle(GetHistSingleMethod(PathResolverFindCalibFile(Form("MuonMomentumCorrections/%s/%s/",m_SagittaRelease.c_str(),m_year.c_str())+ trackNames.at(i) + "_data.root"),Form("p%s_1",trackNames.at(i).c_str())),i);
+
+	    MuonCalibrationAndSmearingTool::setSagittaHistogramsSingle(GetHistSingleMethod(PathResolverFindCalibFile(Form("MuonMomentumCorrections/%s/%s/",m_SagittaRelease.c_str(),m_year.c_str())+ trackNames.at(i) + "_data.root"),Form("p%s_statError",trackNames.at(i).c_str())));
+	    
 	  }
 	  else {
 	    ATH_MSG_ERROR("Unkown sagitta bias map input configuration type ");
@@ -385,10 +391,19 @@ namespace CP {
       }
 
       if(m_SagittaCorrPhaseSpace){
+	if( m_saggitaMapsInputType==MCAST::SagittaInputHistType::NOMINAL) {
         // Load the mc sagitta bias maps 
-        m_sagittaPhaseSpaceCB = std::make_unique<TProfile2D>(*GetHist( PathResolverFindCalibFile(Form("MuonMomentumCorrections/%s/outqDeltamPlots_iter%d/",m_SagittaRelease.c_str(),0) + trackNames.at(0) + "_mc_NoCorr.root"),"inclusive",m_GlobalZScales.at(0)));
-        m_sagittaPhaseSpaceID = std::make_unique<TProfile2D>(*GetHist( PathResolverFindCalibFile(Form("MuonMomentumCorrections/%s/outqDeltamPlots_iter%d/",m_SagittaRelease.c_str(),0) + trackNames.at(1) + "_mc_NoCorr.root"),"inclusive",m_GlobalZScales.at(1)));
-        m_sagittaPhaseSpaceME = std::make_unique<TProfile2D>(*GetHist( PathResolverFindCalibFile(Form("MuonMomentumCorrections/%s/outqDeltamPlots_iter%d/",m_SagittaRelease.c_str(),0) + trackNames.at(2) + "_mc_NoCorr.root"),"inclusive",m_GlobalZScales.at(2)));
+	  m_sagittaPhaseSpaceCB = std::make_unique<TProfile2D>(*GetHist( PathResolverFindCalibFile(Form("MuonMomentumCorrections/%s/outqDeltamPlots_iter%d/",m_SagittaRelease.c_str(),0) + trackNames.at(0) + "_mc_NoCorr.root"),"inclusive",m_GlobalZScales.at(0)));
+	  m_sagittaPhaseSpaceID = std::make_unique<TProfile2D>(*GetHist( PathResolverFindCalibFile(Form("MuonMomentumCorrections/%s/outqDeltamPlots_iter%d/",m_SagittaRelease.c_str(),0) + trackNames.at(1) + "_mc_NoCorr.root"),"inclusive",m_GlobalZScales.at(1)));
+	  m_sagittaPhaseSpaceME = std::make_unique<TProfile2D>(*GetHist( PathResolverFindCalibFile(Form("MuonMomentumCorrections/%s/outqDeltamPlots_iter%d/",m_SagittaRelease.c_str(),0) + trackNames.at(2) + "_mc_NoCorr.root"),"inclusive",m_GlobalZScales.at(2)));
+	}
+	else if (m_saggitaMapsInputType==MCAST::SagittaInputHistType::SINGLE){
+	  
+	  m_sagittaPhaseSpaceCB = std::make_unique<TProfile2D>(*(GetHistSingleMethod(PathResolverFindCalibFile(Form("MuonMomentumCorrections/%s/%s/",m_SagittaRelease.c_str(),m_year.c_str())+ trackNames.at(0) + "_mc.root"),Form("p%s_0",trackNames.at(0).c_str()))));
+	  m_sagittaPhaseSpaceID =  std::make_unique<TProfile2D>(*(GetHistSingleMethod(PathResolverFindCalibFile(Form("MuonMomentumCorrections/%s/%s/",m_SagittaRelease.c_str(),m_year.c_str())+ trackNames.at(1) + "_mc.root"),Form("p%s_0",trackNames.at(1).c_str()))));
+	  m_sagittaPhaseSpaceME =  std::make_unique<TProfile2D>(*(GetHistSingleMethod(PathResolverFindCalibFile(Form("MuonMomentumCorrections/%s/%s/",m_SagittaRelease.c_str(),m_year.c_str())+ trackNames.at(2) + "_mc.root"),Form("p%s_0",trackNames.at(2).c_str()))));
+	}
+	
       }
       else{
         m_sagittaPhaseSpaceCB=nullptr;
@@ -503,10 +518,10 @@ namespace CP {
     double p2=corrM->GetBinContent(binEta,binPhi);
     
     if(m_currentParameters->SagittaBias != MCAST::SystVariation::Default){
-      if(m_currentParameters->SagittaBias == MCAST::SystVariation::Up){
+      if(m_currentParameters->SagittaBias == MCAST::SystVariation::Up || m_currentParameters->SagittaDataStat == MCAST::SystVariation::Up){
         p2=0.5*p2;
       }
-      else if (m_currentParameters->SagittaBias == MCAST::SystVariation::Down){
+      else if (m_currentParameters->SagittaBias == MCAST::SystVariation::Down || m_currentParameters->SagittaDataStat == MCAST::SystVariation::Down){
         p2=-0.5*p2;
       }
     }
@@ -530,10 +545,34 @@ namespace CP {
     return CorrectionCode::Ok;
   }
 
-  CorrectionCode MuonCalibrationAndSmearingTool::applySagittaBiasCorrection(const unsigned int SgCorrType, xAOD::Muon& mu, unsigned int iter, bool stop, bool isMC, InfoHelper& muonInfo) const {
+  CorrectionCode MuonCalibrationAndSmearingTool::applySagittaBiasCorrection(const unsigned int SgCorrType, xAOD::Muon& mu, unsigned int iter, bool stop, bool isMC, InfoHelper& muonInfo,const unsigned int SystCase) const {
 
     if(m_doSagittaMCDistortion && isMC && iter == 0) stop=true;
+    if(isMC && SystCase == MCAST::SagittaSysType::BIAS && m_saggitaMapsInputType==MCAST::SagittaInputHistType::SINGLE && iter >1 ) stop=true;
+
+    //
+    // to be checked if this needs to be enforced only for m_saggitaMapsInputType==MCAST::SagittaInputHistType::SINGLE or also for m_saggitaMapsInputType==MCAST::SagittaInputHistType::NOMINAL
+
+    // for SagittaInputHistType::NOMINAL the iterations can not ever exceed the number of set iterations
+    if (m_saggitaMapsInputType==MCAST::SagittaInputHistType::NOMINAL ) {
+      if (   SgCorrType==MCAST::SagittaCorType::CB  &&   iter > m_SagittaIterations[0] ) 
+	return CorrectionCode::Ok;
+      if (   SgCorrType==MCAST::SagittaCorType::ID  &&   iter > m_SagittaIterations[1] ) 
+	return CorrectionCode::Ok;
+      if (   SgCorrType==MCAST::SagittaCorType::ME  &&   iter > m_SagittaIterations[2] ) 
+	return CorrectionCode::Ok;
+    }
+    // for SagittaInputHistType::SINGLE they can exceed them only for the BIAS and DATASTAT systematics 
+    if (m_saggitaMapsInputType==MCAST::SagittaInputHistType::SINGLE && (SystCase==MCAST::SagittaSysType::NOMINAL|| SystCase==MCAST::SagittaSysType::RHO)) {
+      if (   SgCorrType==MCAST::SagittaCorType::CB  &&   iter >= m_SagittaIterations[0] ) 
+	return CorrectionCode::Ok;
+      if (   SgCorrType==MCAST::SagittaCorType::ID  &&   iter >= m_SagittaIterations[1] ) 
+	return CorrectionCode::Ok;
+      if (   SgCorrType==MCAST::SagittaCorType::ME  &&   iter >= m_SagittaIterations[2] ) 
+	return CorrectionCode::Ok;
+    }
 
+    
     ATH_MSG_DEBUG("Sagitta correction method "<< SgCorrType <<" iter "<<iter);
 
     if(iter == 0){
@@ -610,7 +649,7 @@ namespace CP {
       muonInfo.sagitta_calibrated_ptcb=muonInfo.ptcb;
       iter++;
       if(corr != CorrectionCode::Ok) return corr;
-      if(!stop)  return applySagittaBiasCorrection(MCAST::SagittaCorType::CB, mu, iter, stop, isMC, muonInfo);
+      if(!stop)  return applySagittaBiasCorrection(MCAST::SagittaCorType::CB, mu, iter, stop, isMC, muonInfo,SystCase);
     }
 
     else if(SgCorrType == MCAST::SagittaCorType::ID){
@@ -620,7 +659,7 @@ namespace CP {
       muonInfo.sagitta_calibrated_ptid=muonInfo.ptid;
       iter++;
       if(corr != CorrectionCode::Ok) return corr;
-      if(!stop)return applySagittaBiasCorrection(MCAST::SagittaCorType::ID, mu, iter, stop, isMC, muonInfo);
+      if(!stop)return applySagittaBiasCorrection(MCAST::SagittaCorType::ID, mu, iter, stop, isMC, muonInfo,SystCase);
     }
 
     else if(SgCorrType == MCAST::SagittaCorType::ME){
@@ -630,7 +669,7 @@ namespace CP {
       muonInfo.sagitta_calibrated_ptms=muonInfo.ptms;
       iter++;
       if(corr != CorrectionCode::Ok) return corr;
-      if(!stop) return applySagittaBiasCorrection(MCAST::SagittaCorType::ME, mu, iter, stop, isMC, muonInfo);
+      if(!stop) return applySagittaBiasCorrection(MCAST::SagittaCorType::ME, mu, iter, stop, isMC, muonInfo,SystCase);
     }
 
     else if( SgCorrType == MCAST::SagittaCorType::WEIGHTS){
@@ -691,7 +730,8 @@ namespace CP {
       double tmpDeltaID=0;
       double tmpDeltaMS=0;
 
-      CorrectionCode idOK=applySagittaBiasCorrection(MCAST::SagittaCorType::ID, mu,0,false, isMC, muonInfo);
+      // GB has a question here: shoulnd't the number of iterations here be itersID and itersME instead of 0 ?
+      CorrectionCode idOK=applySagittaBiasCorrection(MCAST::SagittaCorType::ID, mu,0,false, isMC, muonInfo,SystCase);
       TLorentzVector lvIDCorr; lvIDCorr.SetPtEtaPhiM(muonInfo.ptid,muonInfo.eta,muonInfo.phi,mu.m()/1e3);
       if(idOK == CorrectionCode::Ok && tmpPtID!=0 ) tmpDeltaID = ( -tmpPtID +lvIDCorr.Pt() )/ tmpPtID  ;
       else tmpDeltaID=0;
@@ -702,7 +742,7 @@ namespace CP {
       double stored_ptid = muonInfo.ptid;
 
 
-      CorrectionCode meOK=applySagittaBiasCorrection(MCAST::SagittaCorType::ME, mu,0, false, isMC, muonInfo);
+      CorrectionCode meOK=applySagittaBiasCorrection(MCAST::SagittaCorType::ME, mu,0, false, isMC, muonInfo,SystCase);
       TLorentzVector lvMECorr;  lvMECorr.SetPtEtaPhiM(muonInfo.ptms,muonInfo.eta,muonInfo.phi,mu.m()/1e3);
       if(meOK == CorrectionCode::Ok && tmpPtMS!=0 ) tmpDeltaMS = ( -tmpPtMS + lvMECorr.Pt()/1e3 ) /tmpPtMS  ;
       else tmpDeltaMS=0;
@@ -805,30 +845,56 @@ namespace CP {
     }
 
     unsigned int itersCB=0;
-    if(SystCase == MCAST::SagittaSysType::BIAS && m_SagittaIterations.at(0) > 1)
-      itersCB= m_SagittaIterations.at(0) - 1;
-
     unsigned int itersID=0;
-
-    if(SystCase == MCAST::SagittaSysType::BIAS && m_SagittaIterations.at(1) > 1)
-      itersID= m_SagittaIterations.at(1) - 1;
-
     unsigned int itersME=0;
-    if(SystCase == MCAST::SagittaSysType::BIAS && m_SagittaIterations.at(2) > 1)
-      itersME= m_SagittaIterations.at(2) - 1;
+    if( m_saggitaMapsInputType==MCAST::SagittaInputHistType::NOMINAL) {
+      if(SystCase == MCAST::SagittaSysType::BIAS && m_SagittaIterations.at(0) > 1)
+	itersCB= m_SagittaIterations.at(0) - 1;
+      
+      if(SystCase == MCAST::SagittaSysType::BIAS && m_SagittaIterations.at(1) > 1)
+	itersID= m_SagittaIterations.at(1) - 1;
+      
+      if(SystCase == MCAST::SagittaSysType::BIAS && m_SagittaIterations.at(2) > 1)
+	itersME= m_SagittaIterations.at(2) - 1;
+      
+      // In this case this systematic should return the nominal
+      if( SystCase ==  MCAST::SagittaSysType::DATASTAT ){
+	itersCB=0;
+	itersID=0;
+	itersME=0;
+      }
+    }
 
-    // In case one distrots the MC iterations are set to 1. Systamtics willl be calculated based on the full effect.
-    if(m_doSagittaMCDistortion && isMC){
-      itersCB =0; itersID =0;  itersME=0;
+    else if (m_saggitaMapsInputType==MCAST::SagittaInputHistType::SINGLE){
+      if(SystCase == MCAST::SagittaSysType::BIAS && m_SagittaIterations.at(0) > 0)
+	itersCB= m_SagittaIterations.at(0) + 1;
+      
+      if(SystCase == MCAST::SagittaSysType::BIAS && m_SagittaIterations.at(1) > 0)
+	itersID= m_SagittaIterations.at(1) + 1;
+      
+      if(SystCase == MCAST::SagittaSysType::BIAS && m_SagittaIterations.at(2) > 0)
+	itersME= m_SagittaIterations.at(2) + 1;
+      
+      if ( SystCase ==  MCAST::SagittaSysType::DATASTAT && m_SagittaIterations.at(0) > 0)
+	itersCB= m_SagittaIterations.at(0) + 2;
+      if ( SystCase ==  MCAST::SagittaSysType::DATASTAT && m_SagittaIterations.at(1) > 0)
+	itersID= m_SagittaIterations.at(1) + 2;
+      if ( SystCase ==  MCAST::SagittaSysType::DATASTAT && m_SagittaIterations.at(2) > 0)
+	itersME= m_SagittaIterations.at(2) + 2;
     }
+    
+      // In case one distrots the MC iterations are set to 1. Systamtics willl be calculated based on the full effect.
+      if(m_doSagittaMCDistortion && isMC){
+	itersCB =0; itersID =0;  itersME=0;
+      }
 
 
     if(DetType == MCAST::DetectorType::ID){  // Correct the ID momentum
-      return applySagittaBiasCorrection(MCAST::SagittaCorType::ID, mu, itersID, false, isMC, muonInfo);
+      return applySagittaBiasCorrection(MCAST::SagittaCorType::ID, mu, itersID, false, isMC, muonInfo,SystCase);
     }
 
     else if(DetType == MCAST::DetectorType::MS){  // Correct the ME momentum
-      return applySagittaBiasCorrection(MCAST::SagittaCorType::ME, mu, itersME, false, isMC, muonInfo);
+      return applySagittaBiasCorrection(MCAST::SagittaCorType::ME, mu, itersME, false, isMC, muonInfo,SystCase);
     }
     else if(DetType == MCAST::DetectorType::CB){  // Correct the CB momentum;
       if( muonInfo.ptcb == 0) {
@@ -882,12 +948,12 @@ namespace CP {
         ATH_MSG_VERBOSE("Applying sagitta correction for Standalone");
         rho=0;
         if(muonInfo.ptid == 0  && muonInfo.ptms != 0 )  {
-          if(applySagittaBiasCorrection(MCAST::SagittaCorType::ME, mu, itersME, false, isMC, muonInfo)!=CorrectionCode::Ok){
+          if(applySagittaBiasCorrection(MCAST::SagittaCorType::ME, mu, itersME, false, isMC, muonInfo,SystCase)!=CorrectionCode::Ok){
             return CP::CorrectionCode::Error;                                                
           }                                                                                  
         }                                                                                    
         else if(muonInfo.ptid != 0  && muonInfo.ptms == 0 )  {                               
-          if(applySagittaBiasCorrection(MCAST::SagittaCorType::ID, mu, itersID, false, isMC, muonInfo)!=CorrectionCode::Ok){
+          if(applySagittaBiasCorrection(MCAST::SagittaCorType::ID, mu, itersID, false, isMC, muonInfo,SystCase)!=CorrectionCode::Ok){
             return CP::CorrectionCode::Error;
           }
         }
@@ -906,7 +972,7 @@ namespace CP {
 
       ATH_MSG_VERBOSE("Applying CB sagitta correction");
 
-      if(applySagittaBiasCorrection(MCAST::SagittaCorType::CB, mu, itersCB, false, isMC, muonInfo)!=CP::CorrectionCode::Ok){
+      if(applySagittaBiasCorrection(MCAST::SagittaCorType::CB, mu, itersCB, false, isMC, muonInfo,SystCase)!=CP::CorrectionCode::Ok){
         return CP::CorrectionCode::Error;
       }
 
@@ -917,7 +983,7 @@ namespace CP {
 
       ATH_MSG_VERBOSE("Applying Weighted sagitta correction");
 
-      if(applySagittaBiasCorrection(MCAST::SagittaCorType::WEIGHTS, mu, itersCB, false, isMC, muonInfo)!=CP::CorrectionCode::Ok){
+      if(applySagittaBiasCorrection(MCAST::SagittaCorType::WEIGHTS, mu, itersCB, false, isMC, muonInfo,SystCase)!=CP::CorrectionCode::Ok){
         return CP::CorrectionCode::Error;
       }
       else {
@@ -1472,7 +1538,13 @@ namespace CP {
     // Sagitta correction resid bias
     result.insert( SystematicVariation( "MUON_SAGITTA_RESBIAS", 1 ) );
     result.insert( SystematicVariation( "MUON_SAGITTA_RESBIAS", -1 ) );
-
+    
+    if( m_saggitaMapsInputType==MCAST::SagittaInputHistType::SINGLE ){
+      // Sagitta correction resid bias
+      result.insert( SystematicVariation( "MUON_SAGITTA_DATASTAT", 1 ) );
+      result.insert( SystematicVariation( "MUON_SAGITTA_DATASTAT", -1 ) );
+    }
+    
     return result;
 
   }
@@ -1506,7 +1578,7 @@ namespace CP {
     param.ScaleMS_egLoss  = MCAST::SystVariation::Default;
     param.SagittaRho      = MCAST::SystVariation::Default;
     param.SagittaBias     = MCAST::SystVariation::Default;
-
+    param.SagittaDataStat = MCAST::SystVariation::Default;
     // ID systematics
     SystematicVariation syst = systConfig.getSystematicByBaseName( "MUON_ID" );
 
@@ -1518,6 +1590,7 @@ namespace CP {
       param.ScaleMS_egLoss= MCAST::SystVariation::Default;
       param.SagittaRho  = MCAST::SystVariation::Default;
       param.SagittaBias = MCAST::SystVariation::Default;
+      param.SagittaDataStat = MCAST::SystVariation::Default;
     }
     else if( syst == SystematicVariation( "MUON_ID", -1 ) ) {
       param.SmearTypeMS = MCAST::SystVariation::Default;
@@ -1527,6 +1600,7 @@ namespace CP {
       param.ScaleMS_egLoss= MCAST::SystVariation::Default;
       param.SagittaRho  = MCAST::SystVariation::Default;
       param.SagittaBias = MCAST::SystVariation::Default;
+      param.SagittaDataStat = MCAST::SystVariation::Default;
     }
     else if( !syst.empty() ) return SystematicCode::Unsupported;
 
@@ -1541,6 +1615,7 @@ namespace CP {
       param.ScaleMS_egLoss= MCAST::SystVariation::Default;
       param.SagittaRho  = MCAST::SystVariation::Default;
       param.SagittaBias = MCAST::SystVariation::Default;
+      param.SagittaDataStat = MCAST::SystVariation::Default;
     }
     else if( syst == SystematicVariation( "MUON_MS", -1 ) ) {
       param.SmearTypeMS = MCAST::SystVariation::Down;
@@ -1550,6 +1625,7 @@ namespace CP {
       param.ScaleMS_egLoss= MCAST::SystVariation::Default;
       param.SagittaRho  = MCAST::SystVariation::Default;
       param.SagittaBias = MCAST::SystVariation::Default;
+      param.SagittaDataStat = MCAST::SystVariation::Default;
     }
     else if( !syst.empty() ) return SystematicCode::Unsupported;
 
@@ -1564,6 +1640,7 @@ namespace CP {
       param.ScaleMS_egLoss= MCAST::SystVariation::Down;
       param.SagittaRho  = MCAST::SystVariation::Default;
       param.SagittaBias = MCAST::SystVariation::Default;
+      param.SagittaDataStat = MCAST::SystVariation::Default;
     }
     else if( syst == SystematicVariation( "MUON_SCALE", -1 ) ) {
       param.SmearTypeMS = MCAST::SystVariation::Default;
@@ -1573,6 +1650,7 @@ namespace CP {
       param.ScaleMS_egLoss= MCAST::SystVariation::Up;
       param.SagittaRho  = MCAST::SystVariation::Default;
       param.SagittaBias = MCAST::SystVariation::Default;
+      param.SagittaDataStat = MCAST::SystVariation::Default;
     }
     else if( !syst.empty() ) return SystematicCode::Unsupported;
 
@@ -1587,6 +1665,7 @@ namespace CP {
       param.ScaleMS_egLoss= MCAST::SystVariation::Default;
       param.SagittaRho  = MCAST::SystVariation::Default;
       param.SagittaBias = MCAST::SystVariation::Default;
+      param.SagittaDataStat = MCAST::SystVariation::Default;
     }
     else if( syst == SystematicVariation( "MUON_SCALE_ID", -1 ) ) {
       param.SmearTypeMS = MCAST::SystVariation::Default;
@@ -1596,6 +1675,7 @@ namespace CP {
       param.ScaleMS_egLoss= MCAST::SystVariation::Default;
       param.SagittaRho  = MCAST::SystVariation::Default;
       param.SagittaBias = MCAST::SystVariation::Default;
+      param.SagittaDataStat = MCAST::SystVariation::Default;
     }
     else if( !syst.empty() ) return SystematicCode::Unsupported;
 
@@ -1610,6 +1690,7 @@ namespace CP {
       param.ScaleMS_egLoss= MCAST::SystVariation::Default;
       param.SagittaRho  = MCAST::SystVariation::Default;
       param.SagittaBias = MCAST::SystVariation::Default;
+      param.SagittaDataStat = MCAST::SystVariation::Default;
     }
     else if( syst == SystematicVariation( "MUON_SCALE_MS", -1 ) ) {
       param.SmearTypeMS = MCAST::SystVariation::Default;
@@ -1619,6 +1700,7 @@ namespace CP {
       param.ScaleMS_egLoss= MCAST::SystVariation::Default;
       param.SagittaRho  = MCAST::SystVariation::Default;
       param.SagittaBias = MCAST::SystVariation::Default;
+      param.SagittaDataStat = MCAST::SystVariation::Default;
     }
     else if( !syst.empty() ) return SystematicCode::Unsupported;
 
@@ -1633,6 +1715,7 @@ namespace CP {
       param.ScaleMS_egLoss= MCAST::SystVariation::Down;
       param.SagittaRho  = MCAST::SystVariation::Default;
       param.SagittaBias = MCAST::SystVariation::Default;
+      param.SagittaDataStat = MCAST::SystVariation::Default;
     }
     else if( syst == SystematicVariation( "MUON_SCALE_MS_ELOSS", -1 ) ) {
       param.SmearTypeMS = MCAST::SystVariation::Default;
@@ -1642,6 +1725,7 @@ namespace CP {
       param.ScaleMS_egLoss= MCAST::SystVariation::Up;
       param.SagittaRho  = MCAST::SystVariation::Default;
       param.SagittaBias = MCAST::SystVariation::Default;
+      param.SagittaDataStat = MCAST::SystVariation::Default;
     }
     else if( !syst.empty() ) return SystematicCode::Unsupported;
 
@@ -1657,6 +1741,7 @@ namespace CP {
       param.ScaleMS_egLoss= MCAST::SystVariation::Default;
       param.SagittaRho  = MCAST::SystVariation::Down;
       param.SagittaBias = MCAST::SystVariation::Default;
+      param.SagittaDataStat = MCAST::SystVariation::Default;
     }
     else if( syst == SystematicVariation( "MUON_SAGITTA_RHO", -1 ) ) {
       param.SmearTypeMS = MCAST::SystVariation::Default;
@@ -1666,6 +1751,7 @@ namespace CP {
       param.ScaleMS_egLoss= MCAST::SystVariation::Default;
       param.SagittaRho  = MCAST::SystVariation::Up;
       param.SagittaBias = MCAST::SystVariation::Default;
+      param.SagittaDataStat = MCAST::SystVariation::Default;
     }
     else if( !syst.empty() ) return SystematicCode::Unsupported;
 
@@ -1681,6 +1767,7 @@ namespace CP {
       param.ScaleMS_egLoss= MCAST::SystVariation::Default;
       param.SagittaRho  = MCAST::SystVariation::Default;
       param.SagittaBias = MCAST::SystVariation::Down;
+      param.SagittaDataStat = MCAST::SystVariation::Default;
     }
     else if( syst == SystematicVariation( "MUON_SAGITTA_RESBIAS", -1 ) ) {
       param.SmearTypeMS = MCAST::SystVariation::Default;
@@ -1690,10 +1777,39 @@ namespace CP {
       param.ScaleMS_egLoss= MCAST::SystVariation::Default;
       param.SagittaRho  = MCAST::SystVariation::Default;
       param.SagittaBias = MCAST::SystVariation::Up;
+      param.SagittaDataStat = MCAST::SystVariation::Default;
+    }
+    else if( !syst.empty() ) return SystematicCode::Unsupported;
+
+
+    // Sagitta Residual Bias systematics
+    syst = systConfig.getSystematicByBaseName( "MUON_SAGITTA_DATASTAT" );
+
+    if( syst == SystematicVariation( "MUON_SAGITTA_DATASTAT", 1 ) ) {
+      param.SmearTypeMS = MCAST::SystVariation::Default;
+      param.SmearTypeID = MCAST::SystVariation::Default;
+      param.ScaleID       = MCAST::SystVariation::Default;
+      param.ScaleMS_scale = MCAST::SystVariation::Default;
+      param.ScaleMS_egLoss= MCAST::SystVariation::Default;
+      param.SagittaRho  = MCAST::SystVariation::Default;
+      param.SagittaBias = MCAST::SystVariation::Default;
+      param.SagittaDataStat = MCAST::SystVariation::Up;
+    }
+    else if( syst == SystematicVariation( "MUON_SAGITTA_DATASTAT", -1 ) ) {
+      param.SmearTypeMS = MCAST::SystVariation::Default;
+      param.SmearTypeID = MCAST::SystVariation::Default;
+      param.ScaleID       = MCAST::SystVariation::Default;
+      param.ScaleMS_scale = MCAST::SystVariation::Default;
+      param.ScaleMS_egLoss= MCAST::SystVariation::Default;
+      param.SagittaRho  = MCAST::SystVariation::Default;
+      param.SagittaBias = MCAST::SystVariation::Default;
+      param.SagittaDataStat = MCAST::SystVariation::Down;
     }
     else if( !syst.empty() ) return SystematicCode::Unsupported;
 
 
+    
+
     //
     ATH_MSG_DEBUG( "Systematic variation's parameters, SmearTypeID: "   << param.SmearTypeID );
     ATH_MSG_DEBUG( "Systematic variation's parameters, SmearTypeMS: "   << param.SmearTypeMS );
-- 
GitLab


From 37c45c06494f5f23aa921657f72528fda611f640 Mon Sep 17 00:00:00 2001
From: Syed Haider Abidi <syed.haider.abidi@cern.ch>
Date: Thu, 29 Jul 2021 18:49:24 +0200
Subject: [PATCH 02/14] clean up of verbose statement

---
 .../Root/MuonCalibrationAndSmearingTool.cxx   | 70 +++----------------
 1 file changed, 8 insertions(+), 62 deletions(-)

diff --git a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx
index 5479cdf02489..a5d5803f21eb 100644
--- a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx
+++ b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx
@@ -1035,13 +1035,6 @@ namespace CP {
     InfoHelper muonInfo;
 
     // Set pt ID:
-    ATH_MSG_VERBOSE( "Retrieving ElementLink to ID TrackParticle..." );
-    ATH_MSG_VERBOSE( "Setting Pt  [ID]: if no track available, set to 0..." );
-    ATH_MSG_VERBOSE( "mu.isAvailable< ElementLink< xAOD::TrackParticleContainer > >( \"inDetTrackParticleLink\" ) = " << mu.isAvailable< ElementLink< xAOD::TrackParticleContainer > >( "inDetTrackParticleLink" ) );
-    ATH_MSG_VERBOSE( "( mu.inDetTrackParticleLink() == nullptr ) = " << ( mu.inDetTrackParticleLink() == nullptr ) );
-    ATH_MSG_VERBOSE( "mu.inDetTrackParticleLink() = " << mu.inDetTrackParticleLink() );
-    ATH_MSG_VERBOSE( "( mu.inDetTrackParticleLink() ).isValid() = " << ( mu.inDetTrackParticleLink() ).isValid() );
-
     if( ( mu.inDetTrackParticleLink() ).isValid() ) {
             const ElementLink< xAOD::TrackParticleContainer >& id_track = mu.inDetTrackParticleLink();
       muonInfo.ptid = ( !id_track ) ? 0. : ( *id_track )->pt() * MeVtoGeV;
@@ -1051,12 +1044,6 @@ namespace CP {
     }
 
     // Set pt MS:
-    ATH_MSG_VERBOSE( "Retrieving ElementLink to MS TrackParticle..." );
-    ATH_MSG_VERBOSE( "Setting Pt  [MS]: if no track available, set to 0..." );
-    ATH_MSG_VERBOSE( "mu.isAvailable< ElementLink< xAOD::TrackParticleContainer > >( \"extrapolatedMuonSpectrometerTrackParticleLink\" ) = " << mu.isAvailable< ElementLink< xAOD::TrackParticleContainer > >( "extrapolatedMuonSpectrometerTrackParticleLink" ) );
-    ATH_MSG_VERBOSE( "( mu.extrapolatedMuonSpectrometerTrackParticleLink() == nullptr ) = " << ( mu.extrapolatedMuonSpectrometerTrackParticleLink() == nullptr ) );
-    ATH_MSG_VERBOSE( "mu.extrapolatedMuonSpectrometerTrackParticleLink() = " << mu.extrapolatedMuonSpectrometerTrackParticleLink() );
-    ATH_MSG_VERBOSE( "( mu.extrapolatedMuonSpectrometerTrackParticleLink() ).isValid() = " << ( mu.extrapolatedMuonSpectrometerTrackParticleLink() ).isValid() );
     if( ( mu.extrapolatedMuonSpectrometerTrackParticleLink() ).isValid() ) {
           const ElementLink< xAOD::TrackParticleContainer >& ms_track = mu.extrapolatedMuonSpectrometerTrackParticleLink();
       muonInfo.ptms = ( !ms_track ) ? 0. : ( *ms_track )->pt() * MeVtoGeV;
@@ -1067,13 +1054,6 @@ namespace CP {
     }
 
     // Set pt CB:
-    ATH_MSG_VERBOSE( "Retrieving ElementLink to CB TrackParticle..." );
-    ATH_MSG_VERBOSE( "Setting Pt  [CB]: if no track available, set to 0..." );
-    ATH_MSG_VERBOSE( "mu.isAvailable< ElementLink< xAOD::TrackParticleContainer > >( \"primaryTrackParticleLink\" ) = " << mu.isAvailable< ElementLink< xAOD::TrackParticleContainer > >( "primaryTrackParticleLink" ) );
-    ATH_MSG_VERBOSE( "( mu.primaryTrackParticleLink() == nullptr ) = " << ( mu.primaryTrackParticleLink() == nullptr ) );
-    ATH_MSG_VERBOSE( "mu.primaryTrackParticleLink() = " << mu.primaryTrackParticleLink() );
-    ATH_MSG_VERBOSE( "( mu.primaryTrackParticleLink() ).isValid() = " << ( mu.primaryTrackParticleLink() ).isValid() );
-
     if( ( mu.primaryTrackParticleLink() ).isValid() ) {
       const ElementLink< xAOD::TrackParticleContainer >& cb_track = mu.primaryTrackParticleLink();
       muonInfo.ptcb = ( !cb_track ) ? 0. : ( *cb_track )->pt() * MeVtoGeV;
@@ -1096,9 +1076,7 @@ namespace CP {
       muonInfo.charge = mu.charge();
     }
    
-    ATH_MSG_VERBOSE( "Setting Eta [CB]..." );
     muonInfo.eta = mu.eta();
-    ATH_MSG_VERBOSE( "Setting Phi [CB]..." );
     muonInfo.phi = mu.phi();
 
     // Getting detector region
@@ -1143,7 +1121,6 @@ namespace CP {
 
     // Retrieve the event information:
     const xAOD::EventInfo* evtInfo = 0;
-    ATH_MSG_VERBOSE( "Retrieving EventInfo from the EventStore..." );
     if( evtStore()->retrieve( evtInfo, "EventInfo" ).isFailure() ) {
       ATH_MSG_ERROR( "No EventInfo object could be retrieved" );
       ATH_MSG_ERROR( "Random number generation not configured correctly, impossible to determine if dealing with data or MC" );
@@ -1175,7 +1152,7 @@ namespace CP {
         mu.setP4( muonInfo.ptid * GeVtoMeV, muonInfo.eta, muonInfo.phi );
       }
       else {
-	mu.auxdata< float >( "MuonSpectrometerPt" ) = muonInfo.ptms * GeVtoMeV;
+	       mu.auxdata< float >( "MuonSpectrometerPt" ) = muonInfo.ptms * GeVtoMeV;
       }
 
       // SAF specifics
@@ -2023,7 +2000,7 @@ namespace CP {
     double tmpDelta = 1.;
     double outPtID = muonInfo.ptid, outPtMS = muonInfo.ptms, outPtCB = muonInfo.ptcb;
     if( DetType == MCAST::DetectorType::ID ) {
-      ATH_MSG_VERBOSE( "Checking s1_ID = " << scaleID );
+      ATH_MSG_VERBOSE( "Checking s1_ID = " << 1 - scaleID );
       ATH_MSG_VERBOSE( "Double-Checking outPtID = " << outPtID << " at first..." );
       if( outPtID == 0. ) return outPtID;
       //Load the ID scale and smearing that you will use
@@ -2032,7 +2009,7 @@ namespace CP {
       tmpDelta = ( int( inSmearID ) == DEFAULT_INIT_VAL ) ? 1. + muonInfo.smearDeltaID : 1. + inSmearID;
       ATH_MSG_VERBOSE( "Double-Checking inSmearID = " << inSmearID );
       ATH_MSG_VERBOSE( "Double-Checking muonInfo.smearDeltaID = " << muonInfo.smearDeltaID );
-      ATH_MSG_VERBOSE( "Double-Checking tmpDelta = " << tmpDelta );
+      ATH_MSG_VERBOSE( "Double-Checking denominator of smeaering = " << tmpDelta );
       if( m_Trel < MCAST::Release::Rel17_2_Sum13 ) {
         if( m_Tsmear == MCAST::SmearingType::Pt )  outPtID = outPtID * tmpDelta;
         if( m_Tsmear == MCAST::SmearingType::QoverPt ) outPtID = ( tmpDelta == 0 ) ? MCAST_MAX_PT : outPtID / tmpDelta;
@@ -2055,7 +2032,7 @@ namespace CP {
       tmpDelta = ( int( inSmearMS ) == DEFAULT_INIT_VAL ) ? 1. + muonInfo.smearDeltaMS : 1. + inSmearMS;
       ATH_MSG_VERBOSE( "Double-Checking inSmearMS = " << inSmearMS );
       ATH_MSG_VERBOSE( "Double-Checking muonInfo.smearDeltaMS = " << muonInfo.smearDeltaMS );
-      ATH_MSG_VERBOSE( "Double-Checking tmpDelta = " << tmpDelta );
+      ATH_MSG_VERBOSE( "Double-Checking denominator of smeaering = " << tmpDelta );
       //In these releases the smearing was applied to the pt before the scale
       if( m_Trel < MCAST::Release::Rel17_2_Sum13 ) {
         if( m_Tsmear == MCAST::SmearingType::Pt )  outPtMS = outPtMS * tmpDelta;
@@ -2453,6 +2430,7 @@ namespace CP {
         return 0.;
       }
       else {
+        ATH_MSG_VERBOSE( "MS: Using p0MS = " << m_p0_MS[muonInfo.detRegion] << " and p1MS = " << m_p1_MS[muonInfo.detRegion] << " and p2MS = " << m_p2_MS[muonInfo.detRegion] );
         smear = m_p0_MS[muonInfo.detRegion]*muonInfo.g0/muonInfo.ptms + m_p1_MS[muonInfo.detRegion]*muonInfo.g1 + m_p2_MS[muonInfo.detRegion]*muonInfo.g2*muonInfo.ptms;
         return smear;
       }
@@ -2590,14 +2568,6 @@ namespace CP {
     double loc_ptcb = 0;
 
     // Set pt ID:
-
-    ATH_MSG_VERBOSE( "Retrieving ElementLink to ID TrackParticle..." );
-    ATH_MSG_VERBOSE( "Setting Pt  [ID]: if no track available, set to 0..." );
-    ATH_MSG_VERBOSE( "mu.isAvailable< ElementLink< xAOD::TrackParticleContainer > >( \"inDetTrackParticleLink\" ) = " << mu.isAvailable< ElementLink< xAOD::TrackParticleContainer > >( "inDetTrackParticleLink" ) );
-    ATH_MSG_VERBOSE( "( mu.inDetTrackParticleLink() == nullptr ) = " << ( mu.inDetTrackParticleLink() == nullptr ) );
-    ATH_MSG_VERBOSE( "mu.inDetTrackParticleLink() = " << mu.inDetTrackParticleLink() );
-    ATH_MSG_VERBOSE( "( mu.inDetTrackParticleLink() ).isValid() = " << ( mu.inDetTrackParticleLink() ).isValid() );
-
     if( ( mu.inDetTrackParticleLink() ).isValid() ) {
         const ElementLink< xAOD::TrackParticleContainer >& id_track = mu.inDetTrackParticleLink();
         loc_ptid = ( !id_track ) ? 0. : ( *id_track )->pt() * MeVtoGeV;
@@ -2608,13 +2578,6 @@ namespace CP {
     }
 
     // Set pt MS:
-
-    ATH_MSG_VERBOSE( "Retrieving ElementLink to MS TrackParticle..." );
-    ATH_MSG_VERBOSE( "Setting Pt  [MS]: if no track available, set to 0..." );
-    ATH_MSG_VERBOSE( "mu.isAvailable< ElementLink< xAOD::TrackParticleContainer > >( \"extrapolatedMuonSpectrometerTrackParticleLink\" ) = " << mu.isAvailable< ElementLink< xAOD::TrackParticleContainer > >( "extrapolatedMuonSpectrometerTrackParticleLink" ) );
-    ATH_MSG_VERBOSE( "( mu.extrapolatedMuonSpectrometerTrackParticleLink() == nullptr ) = " << ( mu.extrapolatedMuonSpectrometerTrackParticleLink() == nullptr ) );
-    ATH_MSG_VERBOSE( "mu.extrapolatedMuonSpectrometerTrackParticleLink() = " << mu.extrapolatedMuonSpectrometerTrackParticleLink() );
-    ATH_MSG_VERBOSE( "( mu.extrapolatedMuonSpectrometerTrackParticleLink() ).isValid() = " << ( mu.extrapolatedMuonSpectrometerTrackParticleLink() ).isValid() );
     if( ( mu.extrapolatedMuonSpectrometerTrackParticleLink() ).isValid() ) {
         const ElementLink< xAOD::TrackParticleContainer >& ms_track = mu.extrapolatedMuonSpectrometerTrackParticleLink();
         loc_ptms = ( !ms_track ) ? 0. : ( *ms_track )->pt() * MeVtoGeV;
@@ -2624,15 +2587,7 @@ namespace CP {
       loc_ptms = 0.;
     }
 
-    // Set pt CB:
-
-    ATH_MSG_VERBOSE( "Retrieving ElementLink to CB TrackParticle..." );
-    ATH_MSG_VERBOSE( "Setting Pt  [CB]: if no track available, set to 0..." );
-    ATH_MSG_VERBOSE( "mu.isAvailable< ElementLink< xAOD::TrackParticleContainer > >( \"primaryTrackParticleLink\" ) = " << mu.isAvailable< ElementLink< xAOD::TrackParticleContainer > >( "primaryTrackParticleLink" ) );
-    ATH_MSG_VERBOSE( "( mu.primaryTrackParticleLink() == nullptr ) = " << ( mu.primaryTrackParticleLink() == nullptr ) );
-    ATH_MSG_VERBOSE( "mu.primaryTrackParticleLink() = " << mu.primaryTrackParticleLink() );
-    ATH_MSG_VERBOSE( "( mu.primaryTrackParticleLink() ).isValid() = " << ( mu.primaryTrackParticleLink() ).isValid() );
-
+    // Set pt CB;
     if( ( mu.primaryTrackParticleLink() ).isValid() ) {
         const ElementLink< xAOD::TrackParticleContainer >& cb_track = mu.primaryTrackParticleLink();
         loc_ptcb = ( !cb_track ) ? 0. : ( *cb_track )->pt() * MeVtoGeV;
@@ -2659,14 +2614,6 @@ namespace CP {
       loc_detRegion = GetRegion( loc_eta, loc_phi );
     }
 
-    ATH_MSG_VERBOSE( "Getting Expected Resolution: " );
-    ATH_MSG_VERBOSE( " detRegion = "<<loc_detRegion);
-    ATH_MSG_VERBOSE( " ptid      = "<<loc_ptid);
-    ATH_MSG_VERBOSE( " ptms      = "<<loc_ptms);
-    ATH_MSG_VERBOSE( " ptcb      = "<<loc_ptcb);
-    ATH_MSG_VERBOSE( " m_nb_regions = "<<m_nb_regions);
-    ATH_MSG_VERBOSE( " mc = "<<mc);
-
     // Expected resolution in data (or unsmeared MC if second argument is true)
     bool useTan2 = true;
     bool useTan  = false;
@@ -2887,7 +2834,6 @@ namespace CP {
 
   int MuonCalibrationAndSmearingTool::GetRegion( const double eta, const double phi ) const {
 
-    ATH_MSG_VERBOSE( "eta,phi = " << eta << "," << phi );
     int ret_k =-1;
     for( int k=0; k < m_nb_regions; ++k ) {
       if( eta>=m_eta_min[k] && eta<m_eta_max[k] ) {
@@ -2909,8 +2855,8 @@ namespace CP {
       ATH_MSG_DEBUG( "Region corresponding to Eta=" << eta << ", Phi=" << phi << " NOT FOUND!" );
       return -1;
     }
-    if( m_doMacroRegions ) {
-      ATH_MSG_VERBOSE( "INDEX = " << m_MacroRegionIdxMap.find( ret_k )->second );
+    if( m_doMacroRegions ) 
+    {
       return m_MacroRegionIdxMap.find( ret_k )->second;
     }
     return ret_k;
-- 
GitLab


From b4369f83fbda7e2893d664dd46c1c2979a770c97 Mon Sep 17 00:00:00 2001
From: Syed Haider Abidi <syed.haider.abidi@cern.ch>
Date: Mon, 2 Aug 2021 11:15:03 +0200
Subject: [PATCH 03/14] Clean up of systematics implementation

---
 .../Root/MuonCalibrationAndSmearingTool.cxx   | 542 ++++++++++--------
 1 file changed, 299 insertions(+), 243 deletions(-)

diff --git a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx
index a5d5803f21eb..1f8ea57a9b48 100644
--- a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx
+++ b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx
@@ -173,15 +173,15 @@ namespace CP {
     m_IterWeight(tool.m_IterWeight),
     m_SagittaIterations(tool.m_SagittaIterations),
     m_GlobalZScales(tool.m_GlobalZScales){
-    declareProperty( "Year", m_year = "Data16" );
-    declareProperty( "Algo", m_algo = "muons" );
-    declareProperty( "SmearingType", m_type = "q_pT" );
-    declareProperty( "Release", m_release = "Recs2020_03_03" );
-    declareProperty( "ToroidOff", m_toroidOff = false );
+    declareProperty("Year", m_year = "Data16" );
+    declareProperty("Algo", m_algo = "muons" );
+    declareProperty("SmearingType", m_type = "q_pT" );
+    declareProperty("Release", m_release = "Recs2020_03_03" );
+    declareProperty("ToroidOff", m_toroidOff = false );
     declareProperty("AddExtraDecorations", m_extra_decorations = false );
     declareProperty("doExtraSmearing", m_extra_highpt_smearing = false );
     declareProperty("do2StationsHighPt", m_2stations_highpt_smearing = false );
-    declareProperty( "FilesPath", m_FilesPath = "" );
+    declareProperty("FilesPath", m_FilesPath = "" );
     declareProperty("sagittaMapsInputType", m_saggitaMapsInputType=tool.m_saggitaMapsInputType);
     declareProperty("sagittaMapUnitConversion",m_sagittaMapUnitConversion=tool.m_sagittaMapUnitConversion);
   }
@@ -365,50 +365,53 @@ namespace CP {
       std::vector<std::string> trackNames; trackNames.push_back("CB"); trackNames.push_back("ID");  trackNames.push_back("ME");
 
 
-      for( unsigned int i=0; i<m_SagittaIterations.size(); i++){
+      for( unsigned int i=0; i<m_SagittaIterations.size(); i++)
+      {
         ATH_MSG_VERBOSE("Case "<<i<<" track Name "<<trackNames.at(i)<<" and iterations "<<m_SagittaIterations.at(i));
         
-	  if( m_saggitaMapsInputType==MCAST::SagittaInputHistType::NOMINAL) {
-	    for( unsigned int j=0; j< m_SagittaIterations.at(i) ;  j++){
-	      ATH_MSG_VERBOSE("Track "<<i<<" file "<< PathResolverFindCalibFile(Form("MuonMomentumCorrections/%s/outqDeltamPlots_iter%d/",m_SagittaRelease.c_str(),j) + trackNames.at(i) + "_data.root"));
-	      MuonCalibrationAndSmearingTool::setSagittaHistogramsSingle(GetHist( PathResolverFindCalibFile(Form("MuonMomentumCorrections/%s/outqDeltamPlots_iter%d/",m_SagittaRelease.c_str(),j) + trackNames.at(i) + "_data.root"),"inclusive",m_GlobalZScales.at(i)),i);
-	    }}
-	  else if (m_saggitaMapsInputType==MCAST::SagittaInputHistType::SINGLE){
-	    ATH_MSG_VERBOSE("Track "<<i<<" file "<< PathResolverFindCalibFile(Form("MuonMomentumCorrections/%s/%s/",m_SagittaRelease.c_str(),m_year.c_str())+ trackNames.at(i) + "_data.root"));
-
-	    MuonCalibrationAndSmearingTool::setSagittaHistogramsSingle(GetHistSingleMethod(PathResolverFindCalibFile(Form("MuonMomentumCorrections/%s/%s/",m_SagittaRelease.c_str(),m_year.c_str())+ trackNames.at(i) + "_data.root"),Form("p%s_0",trackNames.at(i).c_str())),i);
-
-	    MuonCalibrationAndSmearingTool::setSagittaHistogramsSingle(GetHistSingleMethod(PathResolverFindCalibFile(Form("MuonMomentumCorrections/%s/%s/",m_SagittaRelease.c_str(),m_year.c_str())+ trackNames.at(i) + "_data.root"),Form("p%s_1",trackNames.at(i).c_str())),i);
-
-	    MuonCalibrationAndSmearingTool::setSagittaHistogramsSingle(GetHistSingleMethod(PathResolverFindCalibFile(Form("MuonMomentumCorrections/%s/%s/",m_SagittaRelease.c_str(),m_year.c_str())+ trackNames.at(i) + "_data.root"),Form("p%s_statError",trackNames.at(i).c_str())));
-	    
-	  }
-	  else {
-	    ATH_MSG_ERROR("Unkown sagitta bias map input configuration type ");
-	    return StatusCode::FAILURE;
-	  }
+        if( m_saggitaMapsInputType==MCAST::SagittaInputHistType::NOMINAL) 
+        {
+          for( unsigned int j=0; j< m_SagittaIterations.at(i) ;  j++)
+          {
+            ATH_MSG_VERBOSE("Track "<<i<<" file "<< PathResolverFindCalibFile(Form("MuonMomentumCorrections/%s/outqDeltamPlots_iter%d/",m_SagittaRelease.c_str(),j) + trackNames.at(i) + "_data.root"));
+            MuonCalibrationAndSmearingTool::setSagittaHistogramsSingle(GetHist( PathResolverFindCalibFile(Form("MuonMomentumCorrections/%s/outqDeltamPlots_iter%d/",m_SagittaRelease.c_str(),j) + trackNames.at(i) + "_data.root"),"inclusive",m_GlobalZScales.at(i)),i);
+          }
+        }
+        else if (m_saggitaMapsInputType==MCAST::SagittaInputHistType::SINGLE)
+        {
+          ATH_MSG_VERBOSE("Track "<<i<<" file "<< PathResolverFindCalibFile(Form("MuonMomentumCorrections/%s/%s/",m_SagittaRelease.c_str(),m_year.c_str())+ trackNames.at(i) + "_data.root"));
+
+          MuonCalibrationAndSmearingTool::setSagittaHistogramsSingle(GetHistSingleMethod(PathResolverFindCalibFile(Form("MuonMomentumCorrections/%s/%s/",m_SagittaRelease.c_str(),m_year.c_str())+ trackNames.at(i) + "_data.root"),Form("p%s_0",trackNames.at(i).c_str())),        i);
+          MuonCalibrationAndSmearingTool::setSagittaHistogramsSingle(GetHistSingleMethod(PathResolverFindCalibFile(Form("MuonMomentumCorrections/%s/%s/",m_SagittaRelease.c_str(),m_year.c_str())+ trackNames.at(i) + "_data.root"),Form("p%s_1",trackNames.at(i).c_str())),        i);
+          MuonCalibrationAndSmearingTool::setSagittaHistogramsSingle(GetHistSingleMethod(PathResolverFindCalibFile(Form("MuonMomentumCorrections/%s/%s/",m_SagittaRelease.c_str(),m_year.c_str())+ trackNames.at(i) + "_data.root"),Form("p%s_statError",trackNames.at(i).c_str())),i);  
+        }
+        else {
+          ATH_MSG_ERROR("Unkown sagitta bias map input configuration type ");
+          return StatusCode::FAILURE;
+        }
 	
       }
 
       if(m_SagittaCorrPhaseSpace){
-	if( m_saggitaMapsInputType==MCAST::SagittaInputHistType::NOMINAL) {
-        // Load the mc sagitta bias maps 
-	  m_sagittaPhaseSpaceCB = std::make_unique<TProfile2D>(*GetHist( PathResolverFindCalibFile(Form("MuonMomentumCorrections/%s/outqDeltamPlots_iter%d/",m_SagittaRelease.c_str(),0) + trackNames.at(0) + "_mc_NoCorr.root"),"inclusive",m_GlobalZScales.at(0)));
-	  m_sagittaPhaseSpaceID = std::make_unique<TProfile2D>(*GetHist( PathResolverFindCalibFile(Form("MuonMomentumCorrections/%s/outqDeltamPlots_iter%d/",m_SagittaRelease.c_str(),0) + trackNames.at(1) + "_mc_NoCorr.root"),"inclusive",m_GlobalZScales.at(1)));
-	  m_sagittaPhaseSpaceME = std::make_unique<TProfile2D>(*GetHist( PathResolverFindCalibFile(Form("MuonMomentumCorrections/%s/outqDeltamPlots_iter%d/",m_SagittaRelease.c_str(),0) + trackNames.at(2) + "_mc_NoCorr.root"),"inclusive",m_GlobalZScales.at(2)));
-	}
-	else if (m_saggitaMapsInputType==MCAST::SagittaInputHistType::SINGLE){
-	  
-	  m_sagittaPhaseSpaceCB = std::make_unique<TProfile2D>(*(GetHistSingleMethod(PathResolverFindCalibFile(Form("MuonMomentumCorrections/%s/%s/",m_SagittaRelease.c_str(),m_year.c_str())+ trackNames.at(0) + "_mc.root"),Form("p%s_0",trackNames.at(0).c_str()))));
-	  m_sagittaPhaseSpaceID =  std::make_unique<TProfile2D>(*(GetHistSingleMethod(PathResolverFindCalibFile(Form("MuonMomentumCorrections/%s/%s/",m_SagittaRelease.c_str(),m_year.c_str())+ trackNames.at(1) + "_mc.root"),Form("p%s_0",trackNames.at(1).c_str()))));
-	  m_sagittaPhaseSpaceME =  std::make_unique<TProfile2D>(*(GetHistSingleMethod(PathResolverFindCalibFile(Form("MuonMomentumCorrections/%s/%s/",m_SagittaRelease.c_str(),m_year.c_str())+ trackNames.at(2) + "_mc.root"),Form("p%s_0",trackNames.at(2).c_str()))));
-	}
-	
+        if( m_saggitaMapsInputType==MCAST::SagittaInputHistType::NOMINAL) 
+        {
+          // Load the mc sagitta bias maps 
+          m_sagittaPhaseSpaceCB = std::make_unique<TProfile2D>(*GetHist( PathResolverFindCalibFile(Form("MuonMomentumCorrections/%s/outqDeltamPlots_iter%d/",m_SagittaRelease.c_str(),0) + trackNames.at(0) + "_mc_NoCorr.root"),"inclusive",m_GlobalZScales.at(0)));
+          m_sagittaPhaseSpaceID = std::make_unique<TProfile2D>(*GetHist( PathResolverFindCalibFile(Form("MuonMomentumCorrections/%s/outqDeltamPlots_iter%d/",m_SagittaRelease.c_str(),0) + trackNames.at(1) + "_mc_NoCorr.root"),"inclusive",m_GlobalZScales.at(1)));
+          m_sagittaPhaseSpaceME = std::make_unique<TProfile2D>(*GetHist( PathResolverFindCalibFile(Form("MuonMomentumCorrections/%s/outqDeltamPlots_iter%d/",m_SagittaRelease.c_str(),0) + trackNames.at(2) + "_mc_NoCorr.root"),"inclusive",m_GlobalZScales.at(2)));
+        }
+        else if (m_saggitaMapsInputType==MCAST::SagittaInputHistType::SINGLE)
+        {
+          m_sagittaPhaseSpaceCB = std::make_unique<TProfile2D>(*(GetHistSingleMethod(PathResolverFindCalibFile(Form("MuonMomentumCorrections/%s/%s/",m_SagittaRelease.c_str(),m_year.c_str())+ trackNames.at(0) + "_mc.root"),Form("p%s_0",trackNames.at(0).c_str()))));
+          m_sagittaPhaseSpaceID = std::make_unique<TProfile2D>(*(GetHistSingleMethod(PathResolverFindCalibFile(Form("MuonMomentumCorrections/%s/%s/",m_SagittaRelease.c_str(),m_year.c_str())+ trackNames.at(1) + "_mc.root"),Form("p%s_0",trackNames.at(1).c_str()))));
+          m_sagittaPhaseSpaceME = std::make_unique<TProfile2D>(*(GetHistSingleMethod(PathResolverFindCalibFile(Form("MuonMomentumCorrections/%s/%s/",m_SagittaRelease.c_str(),m_year.c_str())+ trackNames.at(2) + "_mc.root"),Form("p%s_0",trackNames.at(2).c_str()))));
+        }
       }
-      else{
-        m_sagittaPhaseSpaceCB=nullptr;
-        m_sagittaPhaseSpaceID=nullptr;
-        m_sagittaPhaseSpaceME=nullptr;
+      else
+      {
+        m_sagittaPhaseSpaceCB = nullptr;
+        m_sagittaPhaseSpaceID = nullptr;
+        m_sagittaPhaseSpaceME = nullptr;
       }
 
       // Set configuration in case only systematic uncertainty is used. 
@@ -517,65 +520,77 @@ namespace CP {
     }
     double p2=corrM->GetBinContent(binEta,binPhi);
     
-    if(m_currentParameters->SagittaBias != MCAST::SystVariation::Default){
-      if(m_currentParameters->SagittaBias == MCAST::SystVariation::Up || m_currentParameters->SagittaDataStat == MCAST::SystVariation::Up){
-        p2=0.5*p2;
-      }
-      else if (m_currentParameters->SagittaBias == MCAST::SystVariation::Down || m_currentParameters->SagittaDataStat == MCAST::SystVariation::Down){
-        p2=-0.5*p2;
-      }
+    if(m_currentParameters->SagittaBias == MCAST::SystVariation::Up || m_currentParameters->SagittaDataStat == MCAST::SystVariation::Up)
+    {
+      p2 = 0.5*p2;
+    }
+    else if (m_currentParameters->SagittaBias == MCAST::SystVariation::Down || m_currentParameters->SagittaDataStat == MCAST::SystVariation::Down)
+    {
+      p2 = -0.5*p2;
     }
+    
     return p2;
   }
 
   CorrectionCode MuonCalibrationAndSmearingTool::CorrectForCharge(double p2, double& pt, int q, bool isMC,double p2Kin) const {
     
-    if( q==0 ) {
+    if( q==0 ) 
+    {
       ATH_MSG_DEBUG("Muon charge is 0");
       return CorrectionCode::OutOfValidityRange;
     }
+
     double corrPt=pt;
-    if(isMC)
-      corrPt = corrPt/(1 - q*p2*m_sagittaMapUnitConversion*corrPt);
-    else{
-      p2=p2-p2Kin;
+    if(isMC) corrPt = corrPt/(1 - q*p2*m_sagittaMapUnitConversion*corrPt);
+    else
+    {
+      // Remove the kinematic term from MC
+      p2 = p2 - p2Kin;
       corrPt = corrPt/(1 + q*p2*m_sagittaMapUnitConversion*corrPt);
     }
-    pt=corrPt;
+    ATH_MSG_DEBUG("CorrectForCharge - in pT: "<<pt<<" corrPt: "<<corrPt);
+    pt = corrPt;
     return CorrectionCode::Ok;
   }
 
   CorrectionCode MuonCalibrationAndSmearingTool::applySagittaBiasCorrection(const unsigned int SgCorrType, xAOD::Muon& mu, unsigned int iter, bool stop, bool isMC, InfoHelper& muonInfo,const unsigned int SystCase) const {
 
-    if(m_doSagittaMCDistortion && isMC && iter == 0) stop=true;
-    if(isMC && SystCase == MCAST::SagittaSysType::BIAS && m_saggitaMapsInputType==MCAST::SagittaInputHistType::SINGLE && iter >1 ) stop=true;
+
+    ATH_MSG_DEBUG("applySagittaBiasCorrection:: asking for iter "<<iter<<" max CB iter: "<<m_SagittaIterations[0]<<" max ID iter: "<<m_SagittaIterations[1]<<" max ME iter: "<<m_SagittaIterations[2]);
+    if(m_doSagittaMCDistortion && isMC && iter == 0) stop = true;
+
+    // Special case for Bias
+    if(isMC && (SystCase == MCAST::SagittaSysType::BIAS)     && (m_saggitaMapsInputType==MCAST::SagittaInputHistType::SINGLE) && (iter >= 1) ) stop = true;
+    if(isMC && (SystCase == MCAST::SagittaSysType::DATASTAT) && (m_saggitaMapsInputType==MCAST::SagittaInputHistType::SINGLE) && (iter >= 1) ) stop = true;
 
     //
     // to be checked if this needs to be enforced only for m_saggitaMapsInputType==MCAST::SagittaInputHistType::SINGLE or also for m_saggitaMapsInputType==MCAST::SagittaInputHistType::NOMINAL
 
     // for SagittaInputHistType::NOMINAL the iterations can not ever exceed the number of set iterations
-    if (m_saggitaMapsInputType==MCAST::SagittaInputHistType::NOMINAL ) {
-      if (   SgCorrType==MCAST::SagittaCorType::CB  &&   iter > m_SagittaIterations[0] ) 
-	return CorrectionCode::Ok;
-      if (   SgCorrType==MCAST::SagittaCorType::ID  &&   iter > m_SagittaIterations[1] ) 
-	return CorrectionCode::Ok;
-      if (   SgCorrType==MCAST::SagittaCorType::ME  &&   iter > m_SagittaIterations[2] ) 
-	return CorrectionCode::Ok;
+    if (m_saggitaMapsInputType==MCAST::SagittaInputHistType::NOMINAL ) 
+    {
+      if (   SgCorrType==MCAST::SagittaCorType::CB  &&   iter > m_SagittaIterations[0] )  return CorrectionCode::Ok;
+      if (   SgCorrType==MCAST::SagittaCorType::ID  &&   iter > m_SagittaIterations[1] )  return CorrectionCode::Ok;
+      if (   SgCorrType==MCAST::SagittaCorType::ME  &&   iter > m_SagittaIterations[2] )  return CorrectionCode::Ok;
     }
     // for SagittaInputHistType::SINGLE they can exceed them only for the BIAS and DATASTAT systematics 
-    if (m_saggitaMapsInputType==MCAST::SagittaInputHistType::SINGLE && (SystCase==MCAST::SagittaSysType::NOMINAL|| SystCase==MCAST::SagittaSysType::RHO)) {
-      if (   SgCorrType==MCAST::SagittaCorType::CB  &&   iter >= m_SagittaIterations[0] ) 
-	return CorrectionCode::Ok;
-      if (   SgCorrType==MCAST::SagittaCorType::ID  &&   iter >= m_SagittaIterations[1] ) 
-	return CorrectionCode::Ok;
-      if (   SgCorrType==MCAST::SagittaCorType::ME  &&   iter >= m_SagittaIterations[2] ) 
-	return CorrectionCode::Ok;
+    if (m_saggitaMapsInputType==MCAST::SagittaInputHistType::SINGLE && (SystCase==MCAST::SagittaSysType::NOMINAL || SystCase==MCAST::SagittaSysType::RHO)) 
+    {
+      if (   SgCorrType==MCAST::SagittaCorType::CB  &&   iter >= m_SagittaIterations[0] )  return CorrectionCode::Ok;
+      if (   SgCorrType==MCAST::SagittaCorType::ID  &&   iter >= m_SagittaIterations[1] )  return CorrectionCode::Ok;
+      if (   SgCorrType==MCAST::SagittaCorType::ME  &&   iter >= m_SagittaIterations[2] )  return CorrectionCode::Ok;
     }
 
     
-    ATH_MSG_DEBUG("Sagitta correction method "<< SgCorrType <<" iter "<<iter);
-
-    if(iter == 0){
+         if ( SgCorrType==MCAST::SagittaCorType::CB)      ATH_MSG_DEBUG("Sagitta correction method CB iter "<<iter);
+    else if ( SgCorrType==MCAST::SagittaCorType::ID)      ATH_MSG_DEBUG("Sagitta correction method ID iter "<<iter);
+    else if ( SgCorrType==MCAST::SagittaCorType::ME)      ATH_MSG_DEBUG("Sagitta correction method ME iter "<<iter);
+    else if ( SgCorrType==MCAST::SagittaCorType::WEIGHTS) ATH_MSG_DEBUG("Sagitta correction method WEIGHTS iter "<<iter);
+    else if ( SgCorrType==MCAST::SagittaCorType::AUTO)    ATH_MSG_DEBUG("Sagitta correction method AUTO iter "<<iter);
+    else ATH_MSG_DEBUG("Sagitta correction method "<< SgCorrType <<" iter "<<iter);
+
+    if(iter == 0)
+    {
       if( ( mu.primaryTrackParticleLink() ).isValid() ) {
         const ElementLink< xAOD::TrackParticleContainer >& cb_track = mu.primaryTrackParticleLink();
         muonInfo.ptcb = ( !cb_track ) ? 0. : ( *cb_track )->pt() * MeVtoGeV;
@@ -602,11 +617,11 @@ namespace CP {
         const ElementLink< xAOD::TrackParticleContainer >& ms_track = mu.extrapolatedMuonSpectrometerTrackParticleLink();
         muonInfo.ptms = ( !ms_track ) ? 0. : ( *ms_track )->pt() * MeVtoGeV;
       }
-      else
-        muonInfo.ptms = 0.;
+      else muonInfo.ptms = 0.;
     }
 
-    if(muonInfo.ptid==0 && muonInfo.ptms ==0){
+    if(muonInfo.ptid==0 && muonInfo.ptms ==0)
+    {
       ATH_MSG_DEBUG("ME and ID momenta == 0");
       return CorrectionCode::OutOfValidityRange;
     }
@@ -628,12 +643,12 @@ namespace CP {
     double p2PhaseSpaceID=0.0;
     double p2PhaseSpaceME=0.0;
     
-    if(m_SagittaCorrPhaseSpace && m_sagittaPhaseSpaceCB!=nullptr)
-      p2PhaseSpaceCB=m_SagittaCorrPhaseSpace ? sagitta(m_sagittaPhaseSpaceCB.get(), lvCB):0.0;
-    if(m_SagittaCorrPhaseSpace && m_sagittaPhaseSpaceID!=nullptr)
-      p2PhaseSpaceID=m_SagittaCorrPhaseSpace ? sagitta(m_sagittaPhaseSpaceID.get(), lvID):0.0;
-    if(m_SagittaCorrPhaseSpace && m_sagittaPhaseSpaceME!=nullptr)
-      p2PhaseSpaceME=m_SagittaCorrPhaseSpace ? sagitta(m_sagittaPhaseSpaceME.get(), lvME):0.0;
+    // Kinematic terms from MC
+    if(m_SagittaCorrPhaseSpace && m_sagittaPhaseSpaceCB!=nullptr) p2PhaseSpaceCB = m_SagittaCorrPhaseSpace ? sagitta(m_sagittaPhaseSpaceCB.get(), lvCB):0.0;
+    if(m_SagittaCorrPhaseSpace && m_sagittaPhaseSpaceID!=nullptr) p2PhaseSpaceID = m_SagittaCorrPhaseSpace ? sagitta(m_sagittaPhaseSpaceID.get(), lvID):0.0;
+    if(m_SagittaCorrPhaseSpace && m_sagittaPhaseSpaceME!=nullptr) p2PhaseSpaceME = m_SagittaCorrPhaseSpace ? sagitta(m_sagittaPhaseSpaceME.get(), lvME):0.0;
+
+    ATH_MSG_VERBOSE("iter "<<iter<< " CB kin term: "<<p2PhaseSpaceCB<<" ID kin term: "<<p2PhaseSpaceID<<" ME kin term: "<<p2PhaseSpaceME<<" eta "<<muonInfo.eta<<" phi "<<muonInfo.phi<<" charge "<< muonInfo.charge);
 
 
 
@@ -642,96 +657,77 @@ namespace CP {
       return CorrectionCode::OutOfValidityRange;
     }
 
-    if(SgCorrType==MCAST::SagittaCorType::CB) {
+    if(SgCorrType==MCAST::SagittaCorType::CB) 
+    {
       if(muonInfo.ptcb == 0 ) return CorrectionCode::Ok;
+      ATH_MSG_VERBOSE("CB saggita at "<<iter<< " size - "<<m_sagittasCB.size());
       if( iter >=  m_sagittasCB.size())  return CorrectionCode::Ok;
-      CorrectionCode corr = CorrectForCharge( sagitta(m_sagittasCB.at(iter).get(), lvCB)*m_IterWeight, muonInfo.ptcb, q, isMC,p2PhaseSpaceCB);
+      ATH_MSG_VERBOSE("CB saggita at "<<iter<< ": "<<sagitta(m_sagittasCB.at(iter).get(), lvCB)*m_IterWeight);
+
+      CorrectionCode corr = CorrectForCharge( sagitta(m_sagittasCB.at(iter).get(), lvCB)*m_IterWeight, muonInfo.ptcb, q, isMC, p2PhaseSpaceCB);
       muonInfo.sagitta_calibrated_ptcb=muonInfo.ptcb;
       iter++;
       if(corr != CorrectionCode::Ok) return corr;
       if(!stop)  return applySagittaBiasCorrection(MCAST::SagittaCorType::CB, mu, iter, stop, isMC, muonInfo,SystCase);
     }
 
-    else if(SgCorrType == MCAST::SagittaCorType::ID){
+    else if(SgCorrType == MCAST::SagittaCorType::ID)
+    {
       if(muonInfo.ptid == 0 ) return CorrectionCode::Ok;
+      ATH_MSG_VERBOSE("ID saggita at "<<iter<< " size - "<<m_sagittasID.size());
       if( iter >= m_sagittasID.size())  return CorrectionCode::Ok;
       CorrectionCode corr = CorrectForCharge(  sagitta(m_sagittasID.at(iter).get(), lvID)*m_IterWeight, muonInfo.ptid, q, isMC,p2PhaseSpaceID);
+      ATH_MSG_VERBOSE("ID saggita at "<<iter<< ": "<<sagitta(m_sagittasID.at(iter).get(), lvID)*m_IterWeight);
       muonInfo.sagitta_calibrated_ptid=muonInfo.ptid;
       iter++;
       if(corr != CorrectionCode::Ok) return corr;
       if(!stop)return applySagittaBiasCorrection(MCAST::SagittaCorType::ID, mu, iter, stop, isMC, muonInfo,SystCase);
     }
 
-    else if(SgCorrType == MCAST::SagittaCorType::ME){
+    else if(SgCorrType == MCAST::SagittaCorType::ME)
+    {
       if(muonInfo.ptms == 0 ) return CorrectionCode::Ok;
+      ATH_MSG_VERBOSE("ME saggita at "<<iter<< " size - "<<m_sagittasME.size());
       if( iter >=  m_sagittasME.size() )  return CorrectionCode::Ok;
       CorrectionCode corr = CorrectForCharge(sagitta(m_sagittasME.at(iter).get(), lvME)*m_IterWeight, muonInfo.ptms, q, isMC,p2PhaseSpaceME);
+      ATH_MSG_VERBOSE("ME saggita at "<<iter<< ": "<<sagitta(m_sagittasME.at(iter).get(), lvME)*m_IterWeight);
+
       muonInfo.sagitta_calibrated_ptms=muonInfo.ptms;
       iter++;
       if(corr != CorrectionCode::Ok) return corr;
       if(!stop) return applySagittaBiasCorrection(MCAST::SagittaCorType::ME, mu, iter, stop, isMC, muonInfo,SystCase);
     }
 
-    else if( SgCorrType == MCAST::SagittaCorType::WEIGHTS){
+    else if( SgCorrType == MCAST::SagittaCorType::WEIGHTS)
+    {
       const xAOD::TrackParticle* CB_track  = mu.trackParticle( xAOD::Muon::CombinedTrackParticle );
       const xAOD::TrackParticle* ID_track  = mu.trackParticle( xAOD::Muon::InnerDetectorTrackParticle );
-      const xAOD::TrackParticle* ME_track  =  mu.trackParticle( xAOD::Muon::ExtrapolatedMuonSpectrometerTrackParticle );
-
-      double CBqOverPE = 1e10;
-      if( ID_track != nullptr && ME_track != nullptr )
-        CBqOverPE=std::pow(std::sqrt( ID_track->definingParametersCovMatrix()( 4, 4 ) + ME_track->definingParametersCovMatrix()( 4, 4 ))/1e3,2);
-
-      double IDqOverPE =  1e10;
-      if(ID_track!=nullptr)
-        IDqOverPE=std::pow( ID_track->definingParametersCovMatrix()( 4, 4 )/1e3,2);
-
-      double MEqOverPE = 1e10;
-      if(ME_track!=nullptr)
-        MEqOverPE = std::pow( ME_track->definingParametersCovMatrix()( 4, 4 )/1e3,2);
-
+      const xAOD::TrackParticle* ME_track  = mu.trackParticle( xAOD::Muon::ExtrapolatedMuonSpectrometerTrackParticle );
 
       CalcCBWeights( mu, muonInfo );
 
-
-      double ptTilde = (1/IDqOverPE*muonInfo.ptid + 1/MEqOverPE*muonInfo.ptms)/(1/IDqOverPE+1/MEqOverPE);
-      double pTilde = ((1/IDqOverPE)*(1/muonInfo.ptid) + (1/MEqOverPE)*(1/muonInfo.ptms))/(1/IDqOverPE+1/MEqOverPE);
-      double deltaPTilde  =(1/muonInfo.ptcb-pTilde)/(1/muonInfo.ptcb);
-
-      double deltaID=std::abs((1/muonInfo.ptid-1/muonInfo.ptcb)/(1/muonInfo.ptcb));
-      double deltaMS=std::abs((1/muonInfo.ptms-1/muonInfo.ptcb)/(1/muonInfo.ptcb));
-      if(deltaID == 0 ) deltaID=1e-6;
-      if(deltaMS == 0 ) deltaMS=1e-6;
-
-      bool dump=true;
-
-      if(iter==0){
-
-        ATH_MSG_DEBUG("eta "<<muonInfo.eta<<" phi  "<<muonInfo.phi);
-        ATH_MSG_DEBUG(" ptCB: "<<muonInfo.ptcb<<" --> "<<ptTilde<<" diff "<< (muonInfo.ptcb-ptTilde)*100/muonInfo.ptcb<<" %");
-        ATH_MSG_DEBUG(" 1/pT: "<<1/muonInfo.ptcb<<" --> "<<pTilde<<" diff "<< deltaPTilde *100 <<" %"<< "1/pttilde "<< 1/ptTilde);
-        ATH_MSG_DEBUG(" deltaID "<<(muonInfo.ptid-muonInfo.ptcb)*100/muonInfo.ptcb<<" % delta MS "<<(muonInfo.ptms-muonInfo.ptcb)*100/muonInfo.ptcb<<" % ");
-        ATH_MSG_DEBUG(" sigma(q/p) CB "<<CBqOverPE*100<<" ID "<<IDqOverPE*100<<" ME "<<MEqOverPE*100);
-      }
-
-
-      if( iter ==  m_sagittasCB.size() ) {
-        if(dump)
-          ATH_MSG_DEBUG("--> " << muonInfo.ptcb);
+      if( iter == m_sagittasCB.size() ) 
+      {
         return CorrectionCode::Ok;
       }
 
-
-      //float sagittaID=iter >= m_sagittasID->size() ? 0 : sagitta(m_sagittasID->at(iter),lvID)-p2PhaseSpaceID;
-      //float sagittaME=iter >= m_sagittasME->size() ? 0 : sagitta(m_sagittasME->at(iter),lvME)-p2PhaseSpaceME;
-
       double tmpPtID = lvID.Pt();   //muonInfo.ptid;
       double tmpPtMS = lvME.Pt();   //muonInfo.ptms;
 
       double tmpDeltaID=0;
       double tmpDeltaMS=0;
 
-      // GB has a question here: shoulnd't the number of iterations here be itersID and itersME instead of 0 ?
-      CorrectionCode idOK=applySagittaBiasCorrection(MCAST::SagittaCorType::ID, mu,0,false, isMC, muonInfo,SystCase);
+      int itersID = 0;
+      int itersME = 0;
+
+      if ((SystCase != MCAST::SagittaSysType::NOMINAL && SystCase != MCAST::SagittaSysType::RHO))
+      {
+        itersID = iter;
+        itersME = iter;
+      }
+
+
+      CorrectionCode idOK=applySagittaBiasCorrection(MCAST::SagittaCorType::ID, mu, itersID, false, isMC, muonInfo, SystCase);
       TLorentzVector lvIDCorr; lvIDCorr.SetPtEtaPhiM(muonInfo.ptid,muonInfo.eta,muonInfo.phi,mu.m()/1e3);
       if(idOK == CorrectionCode::Ok && tmpPtID!=0 ) tmpDeltaID = ( -tmpPtID +lvIDCorr.Pt() )/ tmpPtID  ;
       else tmpDeltaID=0;
@@ -742,7 +738,7 @@ namespace CP {
       double stored_ptid = muonInfo.ptid;
 
 
-      CorrectionCode meOK=applySagittaBiasCorrection(MCAST::SagittaCorType::ME, mu,0, false, isMC, muonInfo,SystCase);
+      CorrectionCode meOK=applySagittaBiasCorrection(MCAST::SagittaCorType::ME, mu, itersME, false, isMC, muonInfo, SystCase);
       TLorentzVector lvMECorr;  lvMECorr.SetPtEtaPhiM(muonInfo.ptms,muonInfo.eta,muonInfo.phi,mu.m()/1e3);
       if(meOK == CorrectionCode::Ok && tmpPtMS!=0 ) tmpDeltaMS = ( -tmpPtMS + lvMECorr.Pt()/1e3 ) /tmpPtMS  ;
       else tmpDeltaMS=0;
@@ -753,14 +749,6 @@ namespace CP {
       double stored_ptms = muonInfo.ptms;
 
 
-      //double CBsagitta3 = (1/(IDqOverPE*deltaID) * sagittaID
-      //                   + 1/(MEqOverPE*deltaMS) * sagittaME)/(1/(IDqOverPE*deltaID)+1/(MEqOverPE*deltaMS));
-
-     
-
-      double simpleCombPt  = (1/(IDqOverPE*deltaID) * lvIDCorr.Pt() + 
-			      1/(MEqOverPE*deltaMS) * lvMECorr.Pt())/(1/(IDqOverPE*deltaID)+1/(MEqOverPE*deltaMS)); 
-
       // Calculate the stat combination before sagitta bias: 
      
       double chi2Nom=-999;
@@ -824,7 +812,8 @@ namespace CP {
       double statCombPtNom = std::sin(parsCBNom[3])/std::abs(parsCBNom[4]) * MeVtoGeV;
       double statCombPt    = std::sin(parsCB[3])/std::abs(parsCB[4]) * MeVtoGeV;
       muonInfo.ptcb =  muonInfo.ptcb * (1  +  (statCombPt-statCombPtNom)/statCombPtNom ) ;
-      ATH_MSG_VERBOSE(" Poor man's combination "<<simpleCombPt<<" Stat comb "<<statCombPt<<" Stat comb nom "<<" statCombPtNom "<<statCombPtNom ); 
+      // ATH_MSG_VERBOSE(" Poor man's combination "<<simpleCombPt<<" Stat comb "<<statCombPt<<" Stat comb nom "<<" statCombPtNom "<<statCombPtNom ); 
+      ATH_MSG_VERBOSE(" Stat comb "<<statCombPt<<" Stat comb nom "<<" statCombPtNom "<<statCombPtNom ); 
 
       muonInfo.ptid = stored_ptid;
       muonInfo.ptms = stored_ptms;
@@ -847,69 +836,78 @@ namespace CP {
     unsigned int itersCB=0;
     unsigned int itersID=0;
     unsigned int itersME=0;
-    if( m_saggitaMapsInputType==MCAST::SagittaInputHistType::NOMINAL) {
-      if(SystCase == MCAST::SagittaSysType::BIAS && m_SagittaIterations.at(0) > 1)
-	itersCB= m_SagittaIterations.at(0) - 1;
-      
-      if(SystCase == MCAST::SagittaSysType::BIAS && m_SagittaIterations.at(1) > 1)
-	itersID= m_SagittaIterations.at(1) - 1;
-      
-      if(SystCase == MCAST::SagittaSysType::BIAS && m_SagittaIterations.at(2) > 1)
-	itersME= m_SagittaIterations.at(2) - 1;
+    if( m_saggitaMapsInputType==MCAST::SagittaInputHistType::NOMINAL) 
+    {
+      if(SystCase == MCAST::SagittaSysType::BIAS && m_SagittaIterations.at(0) > 1) itersCB = m_SagittaIterations.at(0) - 1;
+      if(SystCase == MCAST::SagittaSysType::BIAS && m_SagittaIterations.at(1) > 1) itersID = m_SagittaIterations.at(1) - 1;
+      if(SystCase == MCAST::SagittaSysType::BIAS && m_SagittaIterations.at(2) > 1) itersME = m_SagittaIterations.at(2) - 1;
       
       // In this case this systematic should return the nominal
-      if( SystCase ==  MCAST::SagittaSysType::DATASTAT ){
-	itersCB=0;
-	itersID=0;
-	itersME=0;
+      if( SystCase ==  MCAST::SagittaSysType::DATASTAT )
+      {
+        itersCB=0;
+        itersID=0;
+        itersME=0;
       }
     }
 
-    else if (m_saggitaMapsInputType==MCAST::SagittaInputHistType::SINGLE){
-      if(SystCase == MCAST::SagittaSysType::BIAS && m_SagittaIterations.at(0) > 0)
-	itersCB= m_SagittaIterations.at(0) + 1;
-      
-      if(SystCase == MCAST::SagittaSysType::BIAS && m_SagittaIterations.at(1) > 0)
-	itersID= m_SagittaIterations.at(1) + 1;
-      
-      if(SystCase == MCAST::SagittaSysType::BIAS && m_SagittaIterations.at(2) > 0)
-	itersME= m_SagittaIterations.at(2) + 1;
-      
-      if ( SystCase ==  MCAST::SagittaSysType::DATASTAT && m_SagittaIterations.at(0) > 0)
-	itersCB= m_SagittaIterations.at(0) + 2;
-      if ( SystCase ==  MCAST::SagittaSysType::DATASTAT && m_SagittaIterations.at(1) > 0)
-	itersID= m_SagittaIterations.at(1) + 2;
-      if ( SystCase ==  MCAST::SagittaSysType::DATASTAT && m_SagittaIterations.at(2) > 0)
-	itersME= m_SagittaIterations.at(2) + 2;
+    else if (m_saggitaMapsInputType==MCAST::SagittaInputHistType::SINGLE)
+    {
+      if(SystCase == MCAST::SagittaSysType::BIAS && m_SagittaIterations.at(0) > 0) itersCB = m_SagittaIterations.at(0) + 0;
+      if(SystCase == MCAST::SagittaSysType::BIAS && m_SagittaIterations.at(1) > 0) itersID = m_SagittaIterations.at(1) + 0;
+      if(SystCase == MCAST::SagittaSysType::BIAS && m_SagittaIterations.at(2) > 0) itersME = m_SagittaIterations.at(2) + 0;
+      if(SystCase == MCAST::SagittaSysType::DATASTAT && m_SagittaIterations.at(0) > 0) itersCB = m_SagittaIterations.at(0) + 1;
+      if(SystCase == MCAST::SagittaSysType::DATASTAT && m_SagittaIterations.at(1) > 0) itersID = m_SagittaIterations.at(1) + 1;
+      if(SystCase == MCAST::SagittaSysType::DATASTAT && m_SagittaIterations.at(2) > 0) itersME = m_SagittaIterations.at(2) + 1;
     }
     
-      // In case one distrots the MC iterations are set to 1. Systamtics willl be calculated based on the full effect.
-      if(m_doSagittaMCDistortion && isMC){
-	itersCB =0; itersID =0;  itersME=0;
-      }
+    // In case one distrots the MC iterations are set to 1. Systamtics willl be calculated based on the full effect.
+    if(m_doSagittaMCDistortion && isMC)
+    {
+      itersCB =0; itersID =0;  itersME=0;
+    }
+    ATH_MSG_VERBOSE("applySagittaBiasCorrectionAuto: itersCB: "<<itersCB<<" itersID: "<<itersID<<" itersME: "<<itersME);
 
 
-    if(DetType == MCAST::DetectorType::ID){  // Correct the ID momentum
+
+    if(DetType == MCAST::DetectorType::ID)
+    { 
+      ATH_MSG_VERBOSE("applySagittaBiasCorrectionAuto: Doing ID correction"); 
+      // Correct the ID momentum
       return applySagittaBiasCorrection(MCAST::SagittaCorType::ID, mu, itersID, false, isMC, muonInfo,SystCase);
     }
 
-    else if(DetType == MCAST::DetectorType::MS){  // Correct the ME momentum
+    else if(DetType == MCAST::DetectorType::MS)
+    { 
+      ATH_MSG_VERBOSE("applySagittaBiasCorrectionAuto: Doing MS correction"); 
+      // Correct the ME momentum
       return applySagittaBiasCorrection(MCAST::SagittaCorType::ME, mu, itersME, false, isMC, muonInfo,SystCase);
     }
-    else if(DetType == MCAST::DetectorType::CB){  // Correct the CB momentum;
-      if( muonInfo.ptcb == 0) {
+    else if(DetType == MCAST::DetectorType::CB)
+    {  
+
+      ATH_MSG_VERBOSE("applySagittaBiasCorrectionAuto: Doing CB correction"); 
+
+      // Correct the CB momentum, in case CB momentume is ZERO
+      // Does this even ever come here?
+      if( muonInfo.ptcb == 0) 
+      {
         ATH_MSG_VERBOSE("Combined pt = 0 correcting separtly ID and ME");
-        if(muonInfo.ptid !=0 && muonInfo.ptms !=0){
+        if(muonInfo.ptid !=0 && muonInfo.ptms !=0)
+        {
           if( applySagittaBiasCorrectionAuto(MCAST::DetectorType::ID, mu, isMC, SystCase, muonInfo) != CorrectionCode::Ok &&
               applySagittaBiasCorrectionAuto(MCAST::DetectorType::MS, mu, isMC, SystCase, muonInfo) != CorrectionCode::Ok ) return CorrectionCode::Error;
         }
-        else if(muonInfo.ptid !=0 ){
+        else if(muonInfo.ptid !=0 )
+        {
           if (applySagittaBiasCorrectionAuto(MCAST::DetectorType::ID, mu, isMC, SystCase, muonInfo) != CorrectionCode::Ok) return CorrectionCode::Error;
         }
-        else if(muonInfo.ptms !=0 ){
+        else if(muonInfo.ptms !=0 )
+        {
           if (applySagittaBiasCorrectionAuto(MCAST::DetectorType::MS, mu, isMC, SystCase, muonInfo) != CorrectionCode::Ok) return CorrectionCode::Error;
         }
-        else {
+        else 
+        {
           return CP::CorrectionCode::Ok;
         }
       }
@@ -918,41 +916,52 @@ namespace CP {
       double rho = m_useFixedRho ? m_fixedRho : 0.0; 
       double nom_central(central), nom_sigmas(sigmas), nom_rho(rho);
 
-      bool isSystematic = (SystCase == MCAST::SagittaSysType::RHO);
+      bool isRhoSystematic = (SystCase == MCAST::SagittaSysType::RHO);
 
-      if(isSystematic) {
+      if(isRhoSystematic) 
+      {
         double sigmaID = ExpectedResolution( MCAST::DetectorType::ID, mu, true ) * muonInfo.ptcb;
         double sigmaMS = ExpectedResolution( MCAST::DetectorType::MS, mu, true ) * muonInfo.ptcb;
         double denominator = muonInfo.ptcb * std::sqrt( sigmaID*sigmaID + sigmaMS*sigmaMS );
-        double res= denominator ? std::sqrt( 2. ) * sigmaID * sigmaMS / denominator : 0.;
+        double res = denominator ? std::sqrt( 2. ) * sigmaID * sigmaMS / denominator : 0.;
 
-        if(m_currentParameters->SagittaRho==MCAST::SystVariation::Up){
-          central=nom_central + std::abs(0.5 * res  * nom_central);
+        if(m_currentParameters->SagittaRho==MCAST::SystVariation::Up)
+        {
+          central = nom_central + std::abs(0.5 * res  * nom_central);
         }
-        else if(m_currentParameters->SagittaRho==MCAST::SystVariation::Down){
-          central=nom_central - std::abs(0.5 * res  * nom_central);
+        else if(m_currentParameters->SagittaRho==MCAST::SystVariation::Down)
+        {
+          central = nom_central - std::abs(0.5 * res  * nom_central);
         }
       }
       
-      if(!m_useFixedRho){
-        sigmas = (std::abs(muonInfo.ptcb - central)/width);
-        nom_sigmas = (std::abs(muonInfo.ptcb - nom_central)/width);
-        rho = 1./sigmas;
-        nom_rho = 1./nom_sigmas;
-        if(sigmas < 1.) rho = 1.;
+      if(!m_useFixedRho)
+      {
+        sigmas      = (std::abs(muonInfo.ptcb - central)/width);
+        nom_sigmas  = (std::abs(muonInfo.ptcb - nom_central)/width);
+        rho         = 1./sigmas;
+        nom_rho     = 1./nom_sigmas;
+
+        if(sigmas < 1.)     rho = 1.;
         if(nom_sigmas < 1.) nom_rho = 1.;
       }
+
+      ATH_MSG_VERBOSE("rho: "<<rho<<" nom rho: "<<nom_rho);
+
       
       // For standalone muons and Silicon associated fowrad do not use the combined
-      if( muonInfo.ptid ==0 || muonInfo.ptms ==0){
+      if( muonInfo.ptid ==0 || muonInfo.ptms ==0)
+      {
         ATH_MSG_VERBOSE("Applying sagitta correction for Standalone");
         rho=0;
-        if(muonInfo.ptid == 0  && muonInfo.ptms != 0 )  {
+        if(muonInfo.ptid == 0  && muonInfo.ptms != 0 )  
+        {
           if(applySagittaBiasCorrection(MCAST::SagittaCorType::ME, mu, itersME, false, isMC, muonInfo,SystCase)!=CorrectionCode::Ok){
             return CP::CorrectionCode::Error;                                                
           }                                                                                  
         }                                                                                    
-        else if(muonInfo.ptid != 0  && muonInfo.ptms == 0 )  {                               
+        else if(muonInfo.ptid != 0  && muonInfo.ptms == 0 )  
+        {                               
           if(applySagittaBiasCorrection(MCAST::SagittaCorType::ID, mu, itersID, false, isMC, muonInfo,SystCase)!=CorrectionCode::Ok){
             return CP::CorrectionCode::Error;
           }
@@ -962,36 +971,50 @@ namespace CP {
           return    CP::CorrectionCode::OutOfValidityRange;
         }
 
-        ATH_MSG_VERBOSE("Final pt "<<muonInfo.ptcb);
         return CP::CorrectionCode::Ok;
       }
 
-      double origPt=muonInfo.ptcb;
-      double ptCB=muonInfo.ptcb;
-      double ptWeight=muonInfo.ptcb;
+      double origPtid = muonInfo.ptid; // Back up cuz the code will modify the MuonInfo.ptcb var
+      double origPtms = muonInfo.ptms; // Back up cuz the code will modify the MuonInfo.ptcb var
+      double origPt   = muonInfo.ptcb; // Back up cuz the code will modify the MuonInfo.ptcb var
+      double ptCB     = muonInfo.ptcb;
+      double ptWeight = muonInfo.ptcb;
 
       ATH_MSG_VERBOSE("Applying CB sagitta correction");
 
-      if(applySagittaBiasCorrection(MCAST::SagittaCorType::CB, mu, itersCB, false, isMC, muonInfo,SystCase)!=CP::CorrectionCode::Ok){
-        return CP::CorrectionCode::Error;
+      if(applySagittaBiasCorrection(MCAST::SagittaCorType::CB, mu, itersCB, false, isMC, muonInfo,SystCase) == CP::CorrectionCode::Ok)
+      {
+        ptCB          = muonInfo.ptcb;
+        muonInfo.ptcb = origPt;
+        
+        // only over write if we are doing rho sys, for dataset and bias we want this to change
+        if(isRhoSystematic)
+        {
+          muonInfo.ptid = origPtid;
+          muonInfo.ptms = origPtms;
+        }
       }
+      else  return CP::CorrectionCode::Error;
 
-      else {
-        ptCB = muonInfo.ptcb;
-        muonInfo.ptcb=origPt;
-      }
 
       ATH_MSG_VERBOSE("Applying Weighted sagitta correction");
 
-      if(applySagittaBiasCorrection(MCAST::SagittaCorType::WEIGHTS, mu, itersCB, false, isMC, muonInfo,SystCase)!=CP::CorrectionCode::Ok){
-        return CP::CorrectionCode::Error;
-      }
-      else {
-        ptWeight =  muonInfo.ptcb;
-        muonInfo.ptcb=origPt;
+      if(applySagittaBiasCorrection(MCAST::SagittaCorType::WEIGHTS, mu, itersCB, false, isMC, muonInfo,SystCase) == CP::CorrectionCode::Ok)
+      {
+        ptWeight      = muonInfo.ptcb;
+        muonInfo.ptcb = origPt;
+        // only over write if we are doing rho sys, for dataset and bias we want this to change
+        if(isRhoSystematic)
+        {
+          muonInfo.ptid = origPtid;
+          muonInfo.ptms = origPtms; 
+        }
       }
+      else  return CP::CorrectionCode::Error;
 
-      if(m_useFixedRho){
+
+      if(m_useFixedRho)
+      {
         ATH_MSG_VERBOSE("Using fixed rho value "<<m_fixedRho);
         rho=m_fixedRho;
       }
@@ -1002,16 +1025,18 @@ namespace CP {
 
       // Rescaling momentum to make sure it is consistent with nominal value (needed for systematics only)
       // Systematics are evaluated around the vanilla momentum value, hence the need to rescale
-      if(isSystematic) {
+      if(isRhoSystematic) 
+      {
         // Nominal case, to determine scaling factor
         double nom_ptcb = nom_rho*ptCB + (1-nom_rho)*ptWeight;
-        double ptcb_ratio = nom_ptcb / origPt;
+        double ptcb_ratio = origPt / nom_ptcb;
         // This specific (systematic) case
         double ptcb = rho*ptCB + (1-rho)*ptWeight;
         // Rescaling
         muonInfo.ptcb = ptcb * ptcb_ratio; 
       }
-      else {
+      else 
+      {
         muonInfo.ptcb = rho*ptCB + (1-rho)*ptWeight;
       }
 
@@ -1136,10 +1161,13 @@ namespace CP {
       }
 
       // Sagitta Correction  specifics
-      if(m_doSagittaCorrection){
+      if(m_doSagittaCorrection)
+      {
         CorrectionCode sgCode = applySagittaBiasCorrectionAuto(MCAST::DetectorType::CB, mu, false, MCAST::SagittaSysType::NOMINAL, muonInfo);
         if(sgCode!=CorrectionCode::Ok) return sgCode;
-        ATH_MSG_VERBOSE("Rho weighting details: stored_pt = " << muonInfo.ptcb );
+        ATH_MSG_VERBOSE("Rho weighting details: stored_pt CB = " << muonInfo.ptcb );
+        ATH_MSG_VERBOSE("Rho weighting details: stored_pt ID = " << muonInfo.ptid );
+        ATH_MSG_VERBOSE("Rho weighting details: stored_pt ME = " << muonInfo.ptms );
         mu.setP4( muonInfo.ptcb * GeVtoMeV, muonInfo.eta, muonInfo.phi );
         ATH_MSG_VERBOSE("Rho weighting details: really_stored_pt = " << mu.pt() );
       }
@@ -1183,8 +1211,8 @@ namespace CP {
     }
 
     CalcCBWeights( mu, muonInfo );
-    ATH_MSG_VERBOSE( "Checking Weights - weightID: " << muonInfo.weightID << " - std::abs( weightID - 1 ): " << std::abs( muonInfo.weightID - 1 ) );
-    ATH_MSG_VERBOSE( "Checking Weights - weightMS: " << muonInfo.weightMS << " - std::abs( weightMS - 1 ): " << std::abs( muonInfo.weightMS - 1 ) );
+    ATH_MSG_VERBOSE( "Checking Weights - weightID: " << muonInfo.weightID );
+    ATH_MSG_VERBOSE( "Checking Weights - weightMS: " << muonInfo.weightMS );
     muonInfo.smearDeltaCB = muonInfo.smearDeltaID * muonInfo.weightID + muonInfo.smearDeltaMS * muonInfo.weightMS;
 
     // Calibrate the pt of the muon:
@@ -1192,37 +1220,62 @@ namespace CP {
     double res_msPt = GeVtoMeV * CalculatePt( MCAST::DetectorType::MS, muonInfo.smearDeltaID, muonInfo.smearDeltaMS, m_currentParameters->ScaleID, m_currentParameters->ScaleMS_scale, m_currentParameters->ScaleMS_egLoss, muonInfo );
     double res_cbPt = GeVtoMeV * CalculatePt( MCAST::DetectorType::CB, muonInfo.smearDeltaID, muonInfo.smearDeltaMS, m_currentParameters->ScaleID, m_currentParameters->ScaleMS_scale, m_currentParameters->ScaleMS_egLoss, muonInfo );
 
+
+    ATH_MSG_VERBOSE( "Calibrated pT ID after correction: " << res_idPt );
+    ATH_MSG_VERBOSE( "Calibrated pT MS after correction: " << res_msPt );
+    ATH_MSG_VERBOSE( "Calibrated pT CB after correction: " << res_cbPt );
+
     if( (mu.muonType()!=xAOD::Muon::SiliconAssociatedForwardMuon) && 
         (m_doSagittaCorrection || m_doSagittaMCDistortion) &&
-        (m_currentParameters->SagittaRho != MCAST::SystVariation::Default || m_currentParameters->SagittaBias != MCAST::SystVariation::Default) ){
-      ATH_MSG_VERBOSE( "Systematic uncertainties for sagitta bias "<< muonInfo.ptcb << res_idPt);
-
+        (m_currentParameters->SagittaRho != MCAST::SystVariation::Default || 
+          m_currentParameters->SagittaBias != MCAST::SystVariation::Default || 
+          m_currentParameters->SagittaDataStat != MCAST::SystVariation::Default) )
+    {
+      ATH_MSG_VERBOSE( "Systematic uncertainties for sagitta bias ");
+      // Overwrtie calibrated with uncalibrated ones 
       muonInfo.ptid = res_idPt * MeVtoGeV;
       muonInfo.ptms = res_msPt * MeVtoGeV;
       muonInfo.ptcb = res_cbPt * MeVtoGeV;
 
-
       CorrectionCode sgCode=CorrectionCode::Ok;
-      if(m_currentParameters->SagittaRho != MCAST::SystVariation::Default){
-        if(m_currentParameters->SagittaRho == MCAST::SystVariation::Up || m_currentParameters->SagittaRho == MCAST::SystVariation::Down){
-          sgCode=applySagittaBiasCorrectionAuto(MCAST::DetectorType::CB, mu, true, MCAST::SagittaSysType::RHO, muonInfo);
+
+
+
+      if(m_currentParameters->SagittaRho != MCAST::SystVariation::Default)
+      {
+        if(m_currentParameters->SagittaRho == MCAST::SystVariation::Up || m_currentParameters->SagittaRho == MCAST::SystVariation::Down)
+        {
+          sgCode = applySagittaBiasCorrectionAuto(MCAST::DetectorType::CB, mu, true, MCAST::SagittaSysType::RHO, muonInfo);
         }
       }
 
-      else if(m_currentParameters->SagittaBias != MCAST::SystVariation::Default){
-        if(m_currentParameters->SagittaBias == MCAST::SystVariation::Up || m_currentParameters->SagittaBias == MCAST::SystVariation::Down){
-          float scale_factor = muonInfo.ptcb / (mu.pt() * MeVtoGeV);
-          sgCode=applySagittaBiasCorrectionAuto(MCAST::DetectorType::CB, mu, true, MCAST::SagittaSysType::BIAS, muonInfo);
-          muonInfo.ptcb = muonInfo.ptcb * scale_factor;
+      else if(m_currentParameters->SagittaBias != MCAST::SystVariation::Default)
+      {
+        if(m_currentParameters->SagittaBias == MCAST::SystVariation::Up || m_currentParameters->SagittaBias == MCAST::SystVariation::Down)
+        {
+          sgCode = applySagittaBiasCorrectionAuto(MCAST::DetectorType::CB, mu, true, MCAST::SagittaSysType::BIAS, muonInfo);
         }
       }
+      else if(m_currentParameters->SagittaDataStat != MCAST::SystVariation::Default)
+      {
+        if(m_currentParameters->SagittaDataStat == MCAST::SystVariation::Up || m_currentParameters->SagittaDataStat == MCAST::SystVariation::Down)
+        {
+          sgCode=applySagittaBiasCorrectionAuto(MCAST::DetectorType::CB, mu, true, MCAST::SagittaSysType::DATASTAT, muonInfo);
+        }
+      }
+
 
       if(sgCode!=CorrectionCode::Ok)
         return sgCode;
 
-      res_idPt=muonInfo.ptid * GeVtoMeV;
-      res_msPt=muonInfo.ptms * GeVtoMeV;
-      res_cbPt=muonInfo.ptcb * GeVtoMeV;
+      res_idPt = muonInfo.ptid * GeVtoMeV;
+      res_msPt = muonInfo.ptms * GeVtoMeV;
+      res_cbPt = muonInfo.ptcb * GeVtoMeV;
+
+      ATH_MSG_VERBOSE( "After ID saggita correction: " << muonInfo.ptid );
+      ATH_MSG_VERBOSE( "After MS saggita correction: " << muonInfo.ptms );
+      ATH_MSG_VERBOSE( "After CB saggita correction: " << muonInfo.ptcb );
+
     }
 
     // Sagitta Distrotion  specifics
@@ -1242,11 +1295,13 @@ namespace CP {
       mu.auxdata< float >( "MuonSpectrometerPt" ) = 0.;
       mu.setP4( res_idPt, muonInfo.eta, muonInfo.phi );
     }
-    else {
+    else 
+    {
       mu.auxdata< float >( "MuonSpectrometerPt" ) = res_msPt;
       double cbPT=0;
 
-      if( ( mu.primaryTrackParticleLink() ).isValid() ) {
+      if( ( mu.primaryTrackParticleLink() ).isValid() ) 
+      {
         const ElementLink< xAOD::TrackParticleContainer >& cb_track = mu.primaryTrackParticleLink();
         cbPT=(*cb_track)->pt() * MeVtoGeV;
       }
@@ -1516,7 +1571,8 @@ namespace CP {
     result.insert( SystematicVariation( "MUON_SAGITTA_RESBIAS", 1 ) );
     result.insert( SystematicVariation( "MUON_SAGITTA_RESBIAS", -1 ) );
     
-    if( m_saggitaMapsInputType==MCAST::SagittaInputHistType::SINGLE ){
+    if( m_saggitaMapsInputType==MCAST::SagittaInputHistType::SINGLE )
+    {
       // Sagitta correction resid bias
       result.insert( SystematicVariation( "MUON_SAGITTA_DATASTAT", 1 ) );
       result.insert( SystematicVariation( "MUON_SAGITTA_DATASTAT", -1 ) );
-- 
GitLab


From 89fbee3858ff857bc3a24ff00da31e23d8f47abd Mon Sep 17 00:00:00 2001
From: Syed Haider Abidi <syed.haider.abidi@cern.ch>
Date: Mon, 2 Aug 2021 13:02:53 +0200
Subject: [PATCH 04/14] MCAST update

---
 .../util/MCAST_Tester.cxx                     | 49 ++++++++++++++-----
 1 file changed, 36 insertions(+), 13 deletions(-)

diff --git a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/util/MCAST_Tester.cxx b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/util/MCAST_Tester.cxx
index 908114aed6c7..29faab5e2158 100644
--- a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/util/MCAST_Tester.cxx
+++ b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/util/MCAST_Tester.cxx
@@ -125,11 +125,11 @@ int main( int argc, char* argv[] ) {
   //::: create the tool handle
   asg::AnaToolHandle<CP::IMuonCalibrationAndSmearingTool> corrTool; //!
   corrTool.setTypeAndName("CP::MuonCalibrationAndSmearingTool/MuonCorrectionTool");
-    //::: set the properties
+  //::: set the properties
   corrTool.setProperty("Year",                  "Data17" );
   //   corrTool.setProperty("Algo",                  "muons" );
   //   corrTool.setProperty("SmearingType",          "q_pT" );
-  corrTool.setProperty("Release",               "Recs2020_03_03" );
+  corrTool.setProperty("Release",               "Recs2021_04_18_Prelim" );
   //corrTool.setProperty("Release",               "Recs2018_05_20" );     
   //   corrTool.setProperty("ToroidOff",             false );
   //   corrTool.setProperty("FilesPath",             "" );
@@ -137,7 +137,7 @@ int main( int argc, char* argv[] ) {
   //   corrTool.setProperty("MinCombPt",             300.0);
   corrTool.setProperty("SagittaCorr",           true);
   //corrTool.setProperty("SagittaRelease",        "sagittaBiasDataAll_30_07_18");
-  corrTool.setProperty("SagittaRelease",        "sagittaBiasDataAll_SingleHistMethod_B_07_04_2021");
+  corrTool.setProperty("SagittaRelease",        "sagittaBiasDataAll_SingleHistMethod_B_17_06_2021_v2");
   corrTool.setProperty("doSagittaMCDistortion", false);
   corrTool.setProperty("SagittaCorrPhaseSpace", false);
   corrTool.setProperty("SagittaIterWeight", 1);
@@ -152,6 +152,23 @@ int main( int argc, char* argv[] ) {
   corrTool.setProperty("sagittaMapUnitConversion",1); 
   //   corrTool.setProperty("useExternalSeed" ,      false);
   //   corrTool.setProperty("externalSeed" ,         0);
+
+
+  // corrTool.setProperty("Year",                  "Data17" );
+  // corrTool.setProperty("Release",               "Recs2021_04_18_Prelim" );
+  // corrTool.setProperty("StatComb",              false);
+  // corrTool.setProperty("SagittaCorr",           true);
+  // corrTool.setProperty("SagittaRelease",        "sagittaBiasDataAll_03_02_19_Data17");
+  // corrTool.setProperty("doSagittaMCDistortion", false);
+  // corrTool.setProperty("SagittaCorrPhaseSpace", true);
+  // corrTool.setProperty("SagittaIterWeight", 1);
+  // corrTool.setProperty("fixedRho",              0.0);
+  // corrTool.setProperty("useFixedRho",           false);
+  // corrTool.setProperty("noEigenDecor" ,         false);
+
+
+  corrTool.setProperty( "OutputLevel", MSG::VERBOSE);
+
   //::: retrieve the tool
   corrTool.retrieve();
 
@@ -260,14 +277,15 @@ int main( int argc, char* argv[] ) {
 
     //Info( APP_NAME, "Number of muons: %i", static_cast< int >( muons->size() ) );
 
-    // create a shallow copy of the muons container
-    std::pair< xAOD::MuonContainer*, xAOD::ShallowAuxContainer* > muons_shallowCopy = xAOD::shallowCopyContainer( *muons );
-
-    xAOD::MuonContainer* muonsCorr = muons_shallowCopy.first;
-
     //::: Loop over systematics
-    for( sysListItr = sysList.begin(); sysListItr != sysList.end(); ++sysListItr ) {
+    for( sysListItr = sysList.begin(); sysListItr != sysList.end(); ++sysListItr ) 
+    {
+      // create a shallow copy of the muons container
+      std::pair< xAOD::MuonContainer*, xAOD::ShallowAuxContainer* > muons_shallowCopy = xAOD::shallowCopyContainer( *muons );
 
+      xAOD::MuonContainer* muonsCorr = muons_shallowCopy.first;
+
+      Info( APP_NAME, "-----------------------------------------------------------");
       Info( APP_NAME, "Looking at %s systematic", ( sysListItr->name() ).c_str() );
 
       //::: Check if systematic is applicable to the correction tool
@@ -277,7 +295,8 @@ int main( int argc, char* argv[] ) {
 
       //for( int i=-0; i<1e3 ; i++) {
       //::: Loop over muon container
-      for( auto muon: *muonsCorr ) {
+      for( auto muon: *muonsCorr ) 
+      {
 
         //::: Select "good" muons:
         // if( ! selTool.passedIDCuts( **mu_itr ) ) {
@@ -325,7 +344,7 @@ int main( int argc, char* argv[] ) {
         }
 
 
-        Info( APP_NAME, "--> CB %g, ID %g, ME %g, author: %d, type: %d",ptCB,ptID,ptME,muon->author(),muon->muonType());
+        Info( APP_NAME, "--> CB %g, ID %g, ME %g, author: %d, type: %d",ptCB/1e3,ptID/1e3,ptME/1e3,muon->author(),muon->muonType());
 
         // either use the correctedCopy call or correct the muon object itself
         if( useCorrectedCopy ) {
@@ -338,7 +357,7 @@ int main( int argc, char* argv[] ) {
           CorrPtCB = mu->pt();
           CorrPtID = mu->auxdata< float >( "InnerDetectorPt" );
           CorrPtMS = mu->auxdata< float >( "MuonSpectrometerPt" );
-          Info( APP_NAME, "Calibrated muon: eta = %g, phi = %g, pt(CB) = %g, pt(ID) = %g, pt(MS) = %g", mu->eta(), mu->phi(), mu->pt(), mu->auxdata< float >( "InnerDetectorPt" ), mu->auxdata< float >( "MuonSpectrometerPt" ) );
+          Info( APP_NAME, "Calibrated muon: eta = %g, phi = %g, pt(CB) = %g, pt(ID) = %g, pt(MS) = %g", mu->eta(), mu->phi(), mu->pt()/1e3, mu->auxdata< float >( "InnerDetectorPt" )/1e3, mu->auxdata< float >( "MuonSpectrometerPt" )/1e3 );
           sysTreeMap[ *sysListItr ]->Fill();
           //::: Delete the calibrated muon:
           delete mu;
@@ -354,11 +373,15 @@ int main( int argc, char* argv[] ) {
           ExpResoCB = corrTool->expectedResolution( "CB", *muon, true );
           ExpResoID = corrTool->expectedResolution( "ID", *muon, true );
           ExpResoMS = corrTool->expectedResolution( "MS", *muon, true );
-          Info( APP_NAME, "Calibrated muon: eta = %g, phi = %g, pt(CB) = %g, pt(ID) = %g, pt(MS) = %g", muon->eta(), muon->phi(), muon->pt()/1e3,CorrPtID,CorrPtMS);
+          Info( APP_NAME, "Calibrated muon: eta = %g, phi = %g, pt(CB) = %g, pt(ID) = %g, pt(MS) = %g", muon->eta(), muon->phi(), muon->pt()/1e3,CorrPtID/1e3,CorrPtMS/1e3);
           Info( APP_NAME, " expReso : ExpResoCB = %g , ExpResoID = %g , ExpResoMS = %g", ExpResoCB, ExpResoID, ExpResoMS);
           sysTreeMap[ *sysListItr ]->Fill();
         }
       }
+      Info( APP_NAME, "-----------------------------------------------------------");
+
+      delete muons_shallowCopy.first;
+      delete muons_shallowCopy.second;
       //}
     }
 
-- 
GitLab


From b8e2a9260ec92c664ba1c6fe12fc241db5a1c492 Mon Sep 17 00:00:00 2001
From: Syed Haider Abidi <syed.haider.abidi@cern.ch>
Date: Mon, 2 Aug 2021 15:12:51 +0200
Subject: [PATCH 05/14] add nom pT to tester

---
 .../util/MCAST_Tester.cxx                     | 63 +++++++++++++++----
 1 file changed, 52 insertions(+), 11 deletions(-)

diff --git a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/util/MCAST_Tester.cxx b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/util/MCAST_Tester.cxx
index 29faab5e2158..abb096d98d3f 100644
--- a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/util/MCAST_Tester.cxx
+++ b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/util/MCAST_Tester.cxx
@@ -6,6 +6,7 @@
 #include <memory>
 #include <cstdlib>
 #include <string>
+#include <map>
 #include "boost/unordered_map.hpp"
 
 // ROOT include(s):
@@ -167,7 +168,10 @@ int main( int argc, char* argv[] ) {
   // corrTool.setProperty("noEigenDecor" ,         false);
 
 
-  corrTool.setProperty( "OutputLevel", MSG::VERBOSE);
+  bool isDebug = false;
+  if(nEvents >= 0 || Ievent >= 0) isDebug = true;
+
+  if(isDebug) corrTool.setProperty( "OutputLevel", MSG::VERBOSE);
 
   //::: retrieve the tool
   corrTool.retrieve();
@@ -220,9 +224,11 @@ int main( int argc, char* argv[] ) {
 
   // branches to be stored
   Float_t InitPtCB( 0. ), InitPtID( 0. ), InitPtMS( 0. );
+  Float_t NomPtCB( 0. ), NomPtID( 0. ), NomPtMS( 0. );
   Float_t CorrPtCB( 0. ), CorrPtID( 0. ), CorrPtMS( 0. );
   Float_t Eta( 0. ), Phi( 0. ), Charge( 0. );
   Float_t ExpResoCB( 0. ), ExpResoID( 0. ), ExpResoMS( 0. );
+  long long unsigned int eventNum (0);
 
   // output file
   TFile* outputFile = TFile::Open( "output.root", "recreate" );
@@ -242,12 +248,17 @@ int main( int argc, char* argv[] ) {
     sysTree->Branch( "CorrPtCB", &CorrPtCB, "CorrPtCB/F" );
     sysTree->Branch( "CorrPtID", &CorrPtID, "CorrPtID/F" );
     sysTree->Branch( "CorrPtMS", &CorrPtMS, "CorrPtMS/F" );
+    sysTree->Branch(  "NomPtCB",  &NomPtCB,  "NomPtCB/F" );
+    sysTree->Branch(  "NomPtID",  &NomPtID,  "NomPtID/F" );
+    sysTree->Branch(  "NomPtMS",  &NomPtMS,  "NomPtMS/F" );
+
     sysTree->Branch( "Eta", &Eta, "Eta/F" );
     sysTree->Branch( "Phi", &Phi, "Phi/F" );
     sysTree->Branch( "Charge", &Charge, "Charge/F" );
     sysTree->Branch( "ExpResoCB", &ExpResoCB, "ExpResoCB/F" );
     sysTree->Branch( "ExpResoID", &ExpResoID, "ExpResoID/F" );
     sysTree->Branch( "ExpResoMS", &ExpResoMS, "ExpResoMS/F" );
+    sysTree->Branch( "eventNum", &eventNum );
 
     sysTreeMap[ *sysListItr ] = sysTree;
   }
@@ -270,13 +281,15 @@ int main( int argc, char* argv[] ) {
       continue;
     }
 
+    eventNum = evtInfo->eventNumber();
+
     //::: Get the Muons from the event:
     const xAOD::MuonContainer* muons = 0;
     if( event.retrieve( muons, "Muons" ) != StatusCode::SUCCESS)
       continue ;
 
     //Info( APP_NAME, "Number of muons: %i", static_cast< int >( muons->size() ) );
-
+    std::map<int, float> cbpt, idpt, mspt;
     //::: Loop over systematics
     for( sysListItr = sysList.begin(); sysListItr != sysList.end(); ++sysListItr ) 
     {
@@ -285,14 +298,19 @@ int main( int argc, char* argv[] ) {
 
       xAOD::MuonContainer* muonsCorr = muons_shallowCopy.first;
 
-      Info( APP_NAME, "-----------------------------------------------------------");
-      Info( APP_NAME, "Looking at %s systematic", ( sysListItr->name() ).c_str() );
+      if(isDebug)
+      {
+        Info( APP_NAME, "-----------------------------------------------------------");
+        Info( APP_NAME, "Looking at %s systematic", ( sysListItr->name() ).c_str() );
+      }
 
       //::: Check if systematic is applicable to the correction tool
       if( corrTool->applySystematicVariation( *sysListItr ) != CP::SystematicCode::Ok ) {
         Error( APP_NAME, "Cannot configure muon calibration tool for systematic" );
       }
 
+      int counter = 0;
+
       //for( int i=-0; i<1e3 ; i++) {
       //::: Loop over muon container
       for( auto muon: *muonsCorr ) 
@@ -322,7 +340,7 @@ int main( int argc, char* argv[] ) {
         Charge = muon->charge();
 
         //::: Print some info about the selected muon:
-        Info( APP_NAME, "Selected muon: eta = %g, phi = %g, pt = %g", muon->eta(), muon->phi(), muon->pt()/1e3 );
+        if(isDebug) Info( APP_NAME, "Selected muon: eta = %g, phi = %g, pt = %g", muon->eta(), muon->phi(), muon->pt()/1e3 );
 
         float ptCB = 0 ;
         if(muon->primaryTrackParticleLink().isValid()){
@@ -330,7 +348,7 @@ int main( int argc, char* argv[] ) {
           ptCB = (!cb_track) ? 0:(*cb_track)->pt();
         }
         else {
-          Info( APP_NAME, "Missing primary track particle link for --> CB %g, author: %d, type: %d",ptCB,muon->author(),muon->muonType());
+          if(isDebug) Info( APP_NAME, "Missing primary track particle link for --> CB %g, author: %d, type: %d",ptCB,muon->author(),muon->muonType());
         }
         float ptID = 0;
         if(muon->inDetTrackParticleLink().isValid()){
@@ -344,7 +362,7 @@ int main( int argc, char* argv[] ) {
         }
 
 
-        Info( APP_NAME, "--> CB %g, ID %g, ME %g, author: %d, type: %d",ptCB/1e3,ptID/1e3,ptME/1e3,muon->author(),muon->muonType());
+        if(isDebug) Info( APP_NAME, "--> CB %g, ID %g, ME %g, author: %d, type: %d",ptCB/1e3,ptID/1e3,ptME/1e3,muon->author(),muon->muonType());
 
         // either use the correctedCopy call or correct the muon object itself
         if( useCorrectedCopy ) {
@@ -357,7 +375,18 @@ int main( int argc, char* argv[] ) {
           CorrPtCB = mu->pt();
           CorrPtID = mu->auxdata< float >( "InnerDetectorPt" );
           CorrPtMS = mu->auxdata< float >( "MuonSpectrometerPt" );
-          Info( APP_NAME, "Calibrated muon: eta = %g, phi = %g, pt(CB) = %g, pt(ID) = %g, pt(MS) = %g", mu->eta(), mu->phi(), mu->pt()/1e3, mu->auxdata< float >( "InnerDetectorPt" )/1e3, mu->auxdata< float >( "MuonSpectrometerPt" )/1e3 );
+
+          if(sysListItr->name().size()==0)
+          {
+            cbpt[counter] = mu->pt();
+            idpt[counter] = mu->auxdata< float >( "InnerDetectorPt" );
+            mspt[counter] = mu->auxdata< float >( "MuonSpectrometerPt" );           
+          }
+          NomPtCB = cbpt[counter];
+          NomPtID = idpt[counter];
+          NomPtMS = mspt[counter];
+
+          if(isDebug) Info( APP_NAME, "Calibrated muon: eta = %g, phi = %g, pt(CB) = %g, pt(ID) = %g, pt(MS) = %g", mu->eta(), mu->phi(), mu->pt()/1e3, mu->auxdata< float >( "InnerDetectorPt" )/1e3, mu->auxdata< float >( "MuonSpectrometerPt" )/1e3 );
           sysTreeMap[ *sysListItr ]->Fill();
           //::: Delete the calibrated muon:
           delete mu;
@@ -373,12 +402,24 @@ int main( int argc, char* argv[] ) {
           ExpResoCB = corrTool->expectedResolution( "CB", *muon, true );
           ExpResoID = corrTool->expectedResolution( "ID", *muon, true );
           ExpResoMS = corrTool->expectedResolution( "MS", *muon, true );
-          Info( APP_NAME, "Calibrated muon: eta = %g, phi = %g, pt(CB) = %g, pt(ID) = %g, pt(MS) = %g", muon->eta(), muon->phi(), muon->pt()/1e3,CorrPtID/1e3,CorrPtMS/1e3);
-          Info( APP_NAME, " expReso : ExpResoCB = %g , ExpResoID = %g , ExpResoMS = %g", ExpResoCB, ExpResoID, ExpResoMS);
+
+          if(sysListItr->name().size()==0)
+          {
+            cbpt[counter] = muon->pt();
+            idpt[counter] = muon->auxdata< float >( "InnerDetectorPt" );
+            mspt[counter] = muon->auxdata< float >( "MuonSpectrometerPt" );           
+          }
+          NomPtCB = cbpt[counter];
+          NomPtID = idpt[counter];
+          NomPtMS = mspt[counter];
+
+          if(isDebug) Info( APP_NAME, "Calibrated muon: eta = %g, phi = %g, pt(CB) = %g, pt(ID) = %g, pt(MS) = %g", muon->eta(), muon->phi(), muon->pt()/1e3,CorrPtID/1e3,CorrPtMS/1e3);
+          if(isDebug) Info( APP_NAME, " expReso : ExpResoCB = %g , ExpResoID = %g , ExpResoMS = %g", ExpResoCB, ExpResoID, ExpResoMS);
           sysTreeMap[ *sysListItr ]->Fill();
         }
+        counter++;
       }
-      Info( APP_NAME, "-----------------------------------------------------------");
+      if(isDebug) Info( APP_NAME, "-----------------------------------------------------------");
 
       delete muons_shallowCopy.first;
       delete muons_shallowCopy.second;
-- 
GitLab


From ce58d2eefe86aafe7971b6eb9623a19e9a78b688 Mon Sep 17 00:00:00 2001
From: Syed Haider Abidi <syed.haider.abidi@cern.ch>
Date: Tue, 3 Aug 2021 18:07:20 +0200
Subject: [PATCH 06/14] systematic fix for saggita bias and addition of corner
 cases where ID/MS tracks are missing

---
 .../MuonCalibrationAndSmearingTool.h          |   3 +
 .../Root/MuonCalibrationAndSmearingTool.cxx   | 119 +++++++++++++++---
 .../util/MCAST_Tester.cxx                     |  21 ++--
 3 files changed, 117 insertions(+), 26 deletions(-)

diff --git a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/MuonMomentumCorrections/MuonCalibrationAndSmearingTool.h b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/MuonMomentumCorrections/MuonCalibrationAndSmearingTool.h
index fe866b406f4c..e2248a8752c3 100644
--- a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/MuonMomentumCorrections/MuonCalibrationAndSmearingTool.h
+++ b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/MuonMomentumCorrections/MuonCalibrationAndSmearingTool.h
@@ -100,6 +100,9 @@ class MuonCalibrationAndSmearingTool : public virtual IMuonCalibrationAndSmearin
       double smearDeltaCB = 0;
       double smearDeltaCBOnly = 0;
       int    sel_category = -1;
+      double uncorrected_ptcb = 0;
+      double uncorrected_ptid = 0;
+      double uncorrected_ptms = 0;
     };
 
   public:
diff --git a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx
index 1f8ea57a9b48..942171c445b2 100644
--- a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx
+++ b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx
@@ -620,9 +620,20 @@ namespace CP {
       else muonInfo.ptms = 0.;
     }
 
-    if(muonInfo.ptid==0 && muonInfo.ptms ==0)
+    if(muonInfo.ptcb == 0 && SgCorrType==MCAST::SagittaCorType::CB)
     {
-      ATH_MSG_DEBUG("ME and ID momenta == 0");
+      ATH_MSG_DEBUG("CB momenta == 0");
+      return CorrectionCode::OutOfValidityRange;
+    }
+    if(muonInfo.ptid == 0 && SgCorrType==MCAST::SagittaCorType::ID)
+    {
+      ATH_MSG_DEBUG("ID momenta == 0");
+      return CorrectionCode::OutOfValidityRange;
+    }
+
+    if(muonInfo.ptms == 0 && SgCorrType==MCAST::SagittaCorType::ME)
+    {
+      ATH_MSG_DEBUG("ME momenta == 0");
       return CorrectionCode::OutOfValidityRange;
     }
 
@@ -740,7 +751,7 @@ namespace CP {
 
       CorrectionCode meOK=applySagittaBiasCorrection(MCAST::SagittaCorType::ME, mu, itersME, false, isMC, muonInfo, SystCase);
       TLorentzVector lvMECorr;  lvMECorr.SetPtEtaPhiM(muonInfo.ptms,muonInfo.eta,muonInfo.phi,mu.m()/1e3);
-      if(meOK == CorrectionCode::Ok && tmpPtMS!=0 ) tmpDeltaMS = ( -tmpPtMS + lvMECorr.Pt()/1e3 ) /tmpPtMS  ;
+      if(meOK == CorrectionCode::Ok && tmpPtMS!=0 ) tmpDeltaMS = ( -tmpPtMS + lvMECorr.Pt() ) /tmpPtMS  ;
       else tmpDeltaMS=0;
       ATH_MSG_VERBOSE( "Shift MS "<<tmpDeltaMS );
       // Now modify the ME covariance matrix 
@@ -948,9 +959,29 @@ namespace CP {
 
       ATH_MSG_VERBOSE("rho: "<<rho<<" nom rho: "<<nom_rho);
 
+      double origPtid = muonInfo.ptid; // Back up cuz the code will modify the MuonInfo.ptcb var
+      double origPtms = muonInfo.ptms; // Back up cuz the code will modify the MuonInfo.ptcb var
+      double origPt   = muonInfo.ptcb; // Back up cuz the code will modify the MuonInfo.ptcb var
+      double ptCB     = muonInfo.ptcb;
+      double ptWeight = muonInfo.ptcb;
+      double ptWeight_unCorr = muonInfo.ptcb;
       
+      // Special case for CB track with missing ID and MS links
+      if(muonInfo.ptcb != 0 && muonInfo.ptid == 0 && muonInfo.ptms == 0)
+      {
+        if(applySagittaBiasCorrection(MCAST::SagittaCorType::CB, mu, itersCB, false, isMC, muonInfo,SystCase) != CP::CorrectionCode::Ok) return CP::CorrectionCode::Error;  
+
+        // Scale for rho systematic
+        if(isRhoSystematic)
+        {
+          muonInfo.ptcb = origPt * muonInfo.ptcb/muonInfo.uncorrected_ptcb;
+        }
+
+        return CP::CorrectionCode::Ok;
+      }
+
       // For standalone muons and Silicon associated fowrad do not use the combined
-      if( muonInfo.ptid ==0 || muonInfo.ptms ==0)
+      if( (muonInfo.ptid ==0 || muonInfo.ptms == 0))
       {
         ATH_MSG_VERBOSE("Applying sagitta correction for Standalone");
         rho=0;
@@ -958,13 +989,25 @@ namespace CP {
         {
           if(applySagittaBiasCorrection(MCAST::SagittaCorType::ME, mu, itersME, false, isMC, muonInfo,SystCase)!=CorrectionCode::Ok){
             return CP::CorrectionCode::Error;                                                
-          }                                                                                  
+          }    
+          // Scale for rho systematic
+          if(isRhoSystematic)
+          {
+            muonInfo.ptms = origPtms * muonInfo.ptms/muonInfo.uncorrected_ptms;
+          }
+
         }                                                                                    
         else if(muonInfo.ptid != 0  && muonInfo.ptms == 0 )  
         {                               
           if(applySagittaBiasCorrection(MCAST::SagittaCorType::ID, mu, itersID, false, isMC, muonInfo,SystCase)!=CorrectionCode::Ok){
             return CP::CorrectionCode::Error;
           }
+          // Scale for rho systematic
+          if(isRhoSystematic)
+          {
+            muonInfo.ptid = origPtid * muonInfo.ptid/muonInfo.uncorrected_ptid;
+          }
+
         }
         else {
           ATH_MSG_DEBUG("All momenta are 0");
@@ -974,11 +1017,23 @@ namespace CP {
         return CP::CorrectionCode::Ok;
       }
 
-      double origPtid = muonInfo.ptid; // Back up cuz the code will modify the MuonInfo.ptcb var
-      double origPtms = muonInfo.ptms; // Back up cuz the code will modify the MuonInfo.ptcb var
-      double origPt   = muonInfo.ptcb; // Back up cuz the code will modify the MuonInfo.ptcb var
-      double ptCB     = muonInfo.ptcb;
-      double ptWeight = muonInfo.ptcb;
+
+      ATH_MSG_VERBOSE("Doing the WEIGHT combination to have the value for later");
+      if(applySagittaBiasCorrection(MCAST::SagittaCorType::WEIGHTS, mu, 1000, false, isMC, muonInfo,SystCase) == CP::CorrectionCode::Ok)
+      {
+        ptWeight_unCorr          = muonInfo.ptcb;
+        muonInfo.ptcb = origPt;
+        
+        // only over write if we are doing rho sys, for dataset and bias we want this to change
+        if(isRhoSystematic)
+        {
+          muonInfo.ptid = origPtid;
+          muonInfo.ptms = origPtms;
+        }
+      }
+      else  return CP::CorrectionCode::Error;
+
+
 
       ATH_MSG_VERBOSE("Applying CB sagitta correction");
 
@@ -1034,6 +1089,20 @@ namespace CP {
         double ptcb = rho*ptCB + (1-rho)*ptWeight;
         // Rescaling
         muonInfo.ptcb = ptcb * ptcb_ratio; 
+      }
+      else if((SystCase == MCAST::SagittaSysType::DATASTAT) || (SystCase == MCAST::SagittaSysType::BIAS))
+      {
+        double nom_ptcb = nom_rho*ptCB + (1-nom_rho)*ptWeight;
+        double org_ptcb = nom_rho*origPt + (1-nom_rho)*ptWeight_unCorr;
+
+        muonInfo.ptcb = origPt * nom_ptcb/org_ptcb; 
+
+        ATH_MSG_VERBOSE("corrected ptCB: "<<ptCB<<" corrected ptWeight: "<<ptWeight);
+        ATH_MSG_VERBOSE("uncorrected ptCB: "<<origPt<<" uncorrected ptWeight: "<<ptWeight_unCorr);
+        ATH_MSG_VERBOSE("combined corrected pt: "<<nom_ptcb<<" combined uncorrected pt: "<<org_ptcb);
+        ATH_MSG_VERBOSE("origPt: "<<origPt<<" corrected pt: "<<muonInfo.ptcb);
+
+
       }
       else 
       {
@@ -1063,29 +1132,35 @@ namespace CP {
     if( ( mu.inDetTrackParticleLink() ).isValid() ) {
             const ElementLink< xAOD::TrackParticleContainer >& id_track = mu.inDetTrackParticleLink();
       muonInfo.ptid = ( !id_track ) ? 0. : ( *id_track )->pt() * MeVtoGeV;
+      muonInfo.uncorrected_ptid = muonInfo.ptid;
     } else {
       ATH_MSG_VERBOSE("The ID track link is not valid - setting pT(ID)=0");
       muonInfo.ptid = 0.;
+      muonInfo.uncorrected_ptid = 0.;
     }
 
     // Set pt MS:
     if( ( mu.extrapolatedMuonSpectrometerTrackParticleLink() ).isValid() ) {
           const ElementLink< xAOD::TrackParticleContainer >& ms_track = mu.extrapolatedMuonSpectrometerTrackParticleLink();
       muonInfo.ptms = ( !ms_track ) ? 0. : ( *ms_track )->pt() * MeVtoGeV;
+      muonInfo.uncorrected_ptms = muonInfo.ptms;
     }
     else{
       ATH_MSG_VERBOSE("No link to extrapolatedMuonSpectrometerTrackParticleLink setting pT(MS)=0");
       muonInfo.ptms = 0.;
+      muonInfo.uncorrected_ptms = 0.;
     }
 
     // Set pt CB:
     if( ( mu.primaryTrackParticleLink() ).isValid() ) {
       const ElementLink< xAOD::TrackParticleContainer >& cb_track = mu.primaryTrackParticleLink();
       muonInfo.ptcb = ( !cb_track ) ? 0. : ( *cb_track )->pt() * MeVtoGeV;
+      muonInfo.uncorrected_ptcb = muonInfo.ptcb;
     }
     else{
       ATH_MSG_VERBOSE("The ID+MS combined track link is not valid - setting pT(ID)=0");
       muonInfo.ptcb = 0.;
+      muonInfo.uncorrected_ptcb = 0;
     }
 
     // Set the charge
@@ -1127,11 +1202,13 @@ namespace CP {
       if( cbCode==CorrectionCode::Ok){
         if(m_doNotUseAMGMATRIXDECOR){
           muonInfo.ptcb = std::sin(muonInfo.cbParsA[3])/std::abs(muonInfo.cbParsA[4]) * MeVtoGeV;
+          muonInfo.uncorrected_ptcb = muonInfo.ptcb;
         }
         else {
           if(!mu.isAvailable < AmgVector(5) >( "StatCombCBPars" )) return CorrectionCode::Error;
           AmgVector(5) parsCB = mu.auxdata < AmgVector(5) >( "StatCombCBPars" );
           muonInfo.ptcb = std::sin(parsCB[3])/std::abs(parsCB[4]) * MeVtoGeV;
+          muonInfo.uncorrected_ptcb = muonInfo.ptcb;
         }
       }
     }
@@ -1225,6 +1302,8 @@ namespace CP {
     ATH_MSG_VERBOSE( "Calibrated pT MS after correction: " << res_msPt );
     ATH_MSG_VERBOSE( "Calibrated pT CB after correction: " << res_cbPt );
 
+    CorrectionCode sgCode = CorrectionCode::Ok;
+
     if( (mu.muonType()!=xAOD::Muon::SiliconAssociatedForwardMuon) && 
         (m_doSagittaCorrection || m_doSagittaMCDistortion) &&
         (m_currentParameters->SagittaRho != MCAST::SystVariation::Default || 
@@ -1237,10 +1316,6 @@ namespace CP {
       muonInfo.ptms = res_msPt * MeVtoGeV;
       muonInfo.ptcb = res_cbPt * MeVtoGeV;
 
-      CorrectionCode sgCode=CorrectionCode::Ok;
-
-
-
       if(m_currentParameters->SagittaRho != MCAST::SystVariation::Default)
       {
         if(m_currentParameters->SagittaRho == MCAST::SystVariation::Up || m_currentParameters->SagittaRho == MCAST::SystVariation::Down)
@@ -1265,7 +1340,9 @@ namespace CP {
       }
 
 
-      if(sgCode!=CorrectionCode::Ok)
+      // This needs to be Error, so that if there is out of validity send by the saggitta bias, atleast the smearing corrections are applied
+      // The out of validity is returned afterwards
+      if(sgCode == CorrectionCode::Error)
         return sgCode;
 
       res_idPt = muonInfo.ptid * GeVtoMeV;
@@ -1334,6 +1411,8 @@ namespace CP {
     ATH_MSG_DEBUG( "Checking Output Muon Info - Pt_MS: " << acc_me_pt(mu));
     ATH_MSG_DEBUG( "Checking Output Muon Info - Pt_CB: " << mu.pt());
 
+    // If saggita was out of validity, return it here
+    if(sgCode == CorrectionCode::OutOfValidityRange) return sgCode;
     // Return gracefully:
     return CorrectionCode::Ok;
   }
@@ -1347,12 +1426,18 @@ namespace CP {
       muonInfo.ptid = inTrk.pt() * MeVtoGeV;
       muonInfo.ptms = 0.;
       muonInfo.ptcb = muonInfo.ptid;
+      muonInfo.uncorrected_ptcb = muonInfo.ptcb;
+      muonInfo.uncorrected_ptid = muonInfo.ptid;
+      muonInfo.uncorrected_ptms = muonInfo.ptms;
       muonInfo.weightID = 1.;
       muonInfo.weightMS = 0.;
     } else if ( DetType == MCAST::DetectorType::MS){ //MS-trk only correction
       muonInfo.ptid = 0.;
       muonInfo.ptms = inTrk.pt() * MeVtoGeV;
       muonInfo.ptcb = muonInfo.ptms;
+      muonInfo.uncorrected_ptcb = muonInfo.ptcb;
+      muonInfo.uncorrected_ptid = muonInfo.ptms;
+      muonInfo.uncorrected_ptid = muonInfo.ptms;
       muonInfo.weightID = 0.;
       muonInfo.weightMS = 1.;
     } else {
@@ -2118,12 +2203,16 @@ namespace CP {
         if( m_Tsmear == MCAST::SmearingType::Pt )  outPtCB = outPtCB * tmpDelta;
         if( m_Tsmear == MCAST::SmearingType::QoverPt ) outPtCB = ( tmpDelta == 0 ) ? MCAST_MAX_PT : outPtCB / tmpDelta;
       }
+      ATH_MSG_VERBOSE( "Org pT " << outPtCB );
+      ATH_MSG_VERBOSE( "Checking scale corr for CB = " << scaleCB );
+      ATH_MSG_VERBOSE( "Checking smearing corr for CB = " << tmpDelta );
 
       outPtCB = ScaleApply( std::abs(outPtCB), scaleCB, 0., muonInfo );
       if( m_Trel >= MCAST::Release::Rel17_2_Sum13 ) {
         if( m_Tsmear == MCAST::SmearingType::Pt )  outPtCB = outPtCB * tmpDelta;
         if( m_Tsmear == MCAST::SmearingType::QoverPt ) outPtCB = ( tmpDelta == 0 ) ? MCAST_MAX_PT : outPtCB / tmpDelta;
       }
+      ATH_MSG_VERBOSE( "Corrected pT " << outPtCB );
 
       return outPtCB;
     }
diff --git a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/util/MCAST_Tester.cxx b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/util/MCAST_Tester.cxx
index abb096d98d3f..8a1170c137ca 100644
--- a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/util/MCAST_Tester.cxx
+++ b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/util/MCAST_Tester.cxx
@@ -151,6 +151,7 @@ int main( int argc, char* argv[] ) {
   corrTool.setProperty("noEigenDecor" ,         false);
   corrTool.setProperty("sagittaMapsInputType" ,         1);
   corrTool.setProperty("sagittaMapUnitConversion",1); 
+  corrTool.setProperty("systematicCorrelationScheme", "Decorr_Scale"); 
   //   corrTool.setProperty("useExternalSeed" ,      false);
   //   corrTool.setProperty("externalSeed" ,         0);
 
@@ -187,9 +188,9 @@ int main( int argc, char* argv[] ) {
   selTool.setTypeAndName("CP::MuonSelectionTool/MuonSelectionTool");
 
   //::: set the properties
-//   selTool.setProperty( "MaxEta",              2.5 );
-//   selTool.setProperty( "MuQuality",           1 );   //corresponds to 0=Tight, 1=Medium, 2=Loose, 3=VeryLoose, 4=HighPt, 5=LowPtEfficiency
-//   selTool.setProperty( "ToroidOff",           false );
+  selTool.setProperty( "MaxEta",              2.5 );
+  selTool.setProperty( "MuQuality",           (int)xAOD::Muon::Loose );   //corresponds to 0=Tight, 1=Medium, 2=Loose, 3=VeryLoose, 4=HighPt, 5=LowPtEfficiency
+  // selTool.setProperty( "ToroidOff",           false );
 //   selTool.setProperty( "TurnOffMomCorr",      false );
 //   selTool.setProperty( "CalibrationRelease",  "PreRec2016_2016-04-13" );
 //   selTool.setProperty( "ExpertDevelopMode",   false );
@@ -315,12 +316,11 @@ int main( int argc, char* argv[] ) {
       //::: Loop over muon container
       for( auto muon: *muonsCorr ) 
       {
-
         //::: Select "good" muons:
-        // if( ! selTool.passedIDCuts( **mu_itr ) ) {
-        //   Info( APP_NAME, "This muon doesn't pass the ID hits quality cuts");
-        //   continue;
-        // }
+        if( ! selTool->accept( *muon ) ) {
+          if(isDebug) Info( APP_NAME, "This muon doesn't pass the ID hits quality cuts");
+          continue;
+        }
 
         //::: Should be using correctedCopy here, testing behaviour of applyCorrection though
         InitPtCB = muon->pt();
@@ -418,17 +418,16 @@ int main( int argc, char* argv[] ) {
           sysTreeMap[ *sysListItr ]->Fill();
         }
         counter++;
+        // break;
       }
       if(isDebug) Info( APP_NAME, "-----------------------------------------------------------");
 
       delete muons_shallowCopy.first;
       delete muons_shallowCopy.second;
-      //}
     }
 
     //::: Close with a message:
-    //if(entry %  1000 ==0 )
-      Info( APP_NAME, "===>>>  done processing event #%i, run #%i %i events processed so far  <<<===", static_cast< int >( evtInfo->eventNumber() ), static_cast< int >( evtInfo->runNumber() ), static_cast< int >( entry + 1 ) );
+    if(entry %  1000 ==0 ) Info( APP_NAME, "===>>>  done processing event #%i, run #%i %i events processed so far  <<<===", static_cast< int >( evtInfo->eventNumber() ), static_cast< int >( evtInfo->runNumber() ), static_cast< int >( entry + 1 ) );
   }
 
   for( sysListItr = sysList.begin(); sysListItr != sysList.end(); ++sysListItr ) {
-- 
GitLab


From 531b3360bba793e45d8c57ab699e560cc410e42b Mon Sep 17 00:00:00 2001
From: Syed Haider Abidi <syed.haider.abidi@cern.ch>
Date: Thu, 5 Aug 2021 19:06:03 +0200
Subject: [PATCH 07/14] clean up and adding back rho variation to the resbias
 sys

---
 .../MuonCalibrationAndSmearingTool.h          |   7 +-
 .../Root/MuonCalibrationAndSmearingTool.cxx   | 446 ++++++++++--------
 .../util/MCAST_Tester.cxx                     |   2 +-
 3 files changed, 253 insertions(+), 202 deletions(-)

diff --git a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/MuonMomentumCorrections/MuonCalibrationAndSmearingTool.h b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/MuonMomentumCorrections/MuonCalibrationAndSmearingTool.h
index e2248a8752c3..7357b022a3b4 100644
--- a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/MuonMomentumCorrections/MuonCalibrationAndSmearingTool.h
+++ b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/MuonMomentumCorrections/MuonCalibrationAndSmearingTool.h
@@ -124,8 +124,8 @@ class MuonCalibrationAndSmearingTool : public virtual IMuonCalibrationAndSmearin
     // Expert method to apply the MC correction on a modifyable trackParticle for ID- or MS-only corrections
     virtual CorrectionCode applyCorrectionTrkOnly( xAOD::TrackParticle& inTrk, const int DetType ) const;
 
-    virtual CorrectionCode applyStatCombination( const ElementLink< xAOD::TrackParticleContainer >& inDetTrackParticle,
-                                                 const ElementLink< xAOD::TrackParticleContainer >& extrTrackParticle ,
+    virtual CorrectionCode applyStatCombination( AmgVector(5) parsID, AmgSymMatrix(5) covID,
+                                                 AmgVector(5) parsMS, AmgSymMatrix(5) covMS,
                                                  int charge,
                                                  AmgVector(5)& parsCB,
                                                  AmgSymMatrix(5)& covCB,
@@ -187,6 +187,9 @@ class MuonCalibrationAndSmearingTool : public virtual IMuonCalibrationAndSmearin
       double SagittaDataStat;
     };
 
+    bool m_expertMode;
+    bool m_expertMode_isData;
+
     bool  m_useExternalSeed;
     int   m_externalSeed;
 
diff --git a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx
index 942171c445b2..46649d24eaf9 100644
--- a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx
+++ b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx
@@ -72,6 +72,9 @@ namespace CP {
     declareProperty("sagittaMapsInputType", m_saggitaMapsInputType=MCAST::SagittaInputHistType::NOMINAL);
     declareProperty("sagittaMapUnitConversion",m_sagittaMapUnitConversion=1e-3);
     declareProperty("systematicCorrelationScheme", m_sysScheme = "Corr_Scale");
+    declareProperty("expertMode", m_expertMode = false);
+    declareProperty("expertSetting_isData", m_expertMode_isData = false);
+
     m_SagittaIterations.push_back(0);
     m_SagittaIterations.push_back(0);
     m_SagittaIterations.push_back(0);
@@ -591,33 +594,53 @@ namespace CP {
 
     if(iter == 0)
     {
-      if( ( mu.primaryTrackParticleLink() ).isValid() ) {
-        const ElementLink< xAOD::TrackParticleContainer >& cb_track = mu.primaryTrackParticleLink();
-        muonInfo.ptcb = ( !cb_track ) ? 0. : ( *cb_track )->pt() * MeVtoGeV;
+      if(m_expertMode)
+      {
+        static const SG::AuxElement::Accessor<float> id_pt("expert_ptid");
+        static const SG::AuxElement::Accessor<float> ms_pt("expert_ptms");
+        static const SG::AuxElement::Accessor<float> cb_pt("expert_ptcb");
+
+        muonInfo.charge = mu.charge();
+        muonInfo.ptid = id_pt(mu) * MeVtoGeV;
+        muonInfo.uncorrected_ptid = muonInfo.ptid;
+
+        muonInfo.ptms = ms_pt(mu) * MeVtoGeV;
+        muonInfo.uncorrected_ptms = muonInfo.ptms;
+
+        muonInfo.ptcb = cb_pt(mu) * MeVtoGeV;
+        muonInfo.uncorrected_ptcb = muonInfo.ptcb;
+
       }
-      else muonInfo.ptcb = 0.;
-      if( m_useStatComb && muonInfo.ptcb > m_StatCombPtThreshold && isBadMuon(mu, muonInfo)) {
-        if(m_doNotUseAMGMATRIXDECOR){
-          muonInfo.ptcb = std::sin(muonInfo.cbParsA[3])/std::abs(muonInfo.cbParsA[4]) * MeVtoGeV;
-        }
-        else {
-          if(!mu.isAvailable < AmgVector(5) >( "StatCombCBPars" )) return CorrectionCode::Error;
-          AmgVector(5) parsCB = mu.auxdata < AmgVector(5) >( "StatCombCBPars" );
-          muonInfo.ptcb = std::sin(parsCB[3])/std::abs(parsCB[4]) * MeVtoGeV;
+      else
+      {
+        if( ( mu.primaryTrackParticleLink() ).isValid() ) {
+          const ElementLink< xAOD::TrackParticleContainer >& cb_track = mu.primaryTrackParticleLink();
+          muonInfo.ptcb = ( !cb_track ) ? 0. : ( *cb_track )->pt() * MeVtoGeV;
+        }
+        else muonInfo.ptcb = 0.;
+        if( m_useStatComb && muonInfo.ptcb > m_StatCombPtThreshold && isBadMuon(mu, muonInfo)) {
+          if(m_doNotUseAMGMATRIXDECOR){
+            muonInfo.ptcb = std::sin(muonInfo.cbParsA[3])/std::abs(muonInfo.cbParsA[4]) * MeVtoGeV;
+          }
+          else {
+            if(!mu.isAvailable < AmgVector(5) >( "StatCombCBPars" )) return CorrectionCode::Error;
+            AmgVector(5) parsCB = mu.auxdata < AmgVector(5) >( "StatCombCBPars" );
+            muonInfo.ptcb = std::sin(parsCB[3])/std::abs(parsCB[4]) * MeVtoGeV;
+          }
         }
-      }
 
-      if( ( mu.inDetTrackParticleLink() ).isValid() ) {
-        const ElementLink< xAOD::TrackParticleContainer >& id_track = mu.inDetTrackParticleLink();
-        muonInfo.ptid = ( !id_track ) ? 0. : ( *id_track )->pt() * MeVtoGeV;
-      }
-      else muonInfo.ptid = 0.;
+        if( ( mu.inDetTrackParticleLink() ).isValid() ) {
+          const ElementLink< xAOD::TrackParticleContainer >& id_track = mu.inDetTrackParticleLink();
+          muonInfo.ptid = ( !id_track ) ? 0. : ( *id_track )->pt() * MeVtoGeV;
+        }
+        else muonInfo.ptid = 0.;
 
-      if( ( mu.extrapolatedMuonSpectrometerTrackParticleLink() ).isValid() ) {
-        const ElementLink< xAOD::TrackParticleContainer >& ms_track = mu.extrapolatedMuonSpectrometerTrackParticleLink();
-        muonInfo.ptms = ( !ms_track ) ? 0. : ( *ms_track )->pt() * MeVtoGeV;
+        if( ( mu.extrapolatedMuonSpectrometerTrackParticleLink() ).isValid() ) {
+          const ElementLink< xAOD::TrackParticleContainer >& ms_track = mu.extrapolatedMuonSpectrometerTrackParticleLink();
+          muonInfo.ptms = ( !ms_track ) ? 0. : ( *ms_track )->pt() * MeVtoGeV;
+        }
+        else muonInfo.ptms = 0.;
       }
-      else muonInfo.ptms = 0.;
     }
 
     if(muonInfo.ptcb == 0 && SgCorrType==MCAST::SagittaCorType::CB)
@@ -738,13 +761,54 @@ namespace CP {
       }
 
 
+      // Calculate the stat combination before sagitta bias: 
+     
+      double chi2Nom=-999;
+      AmgVector(5) parsCBNom, parsCB;
+      AmgSymMatrix(5) covCBNom, covCB;
+      int charge = mu.charge();
+
+      AmgVector(5) parsID;
+      AmgVector(5) parsMS;
+      AmgSymMatrix(5) covID;
+      AmgSymMatrix(5) covMS;
+
+      if(m_expertMode)
+      {
+        parsCBNom = mu.auxdata<AmgVector(5)>("CBParam");
+        parsID = mu.auxdata<AmgVector(5)>("IDParam");
+        parsMS = mu.auxdata<AmgVector(5)>("MSParam");
+        covCBNom = mu.auxdata<AmgSymMatrix(5)>("CBCov");
+        covID = mu.auxdata<AmgSymMatrix(5)>("IDCov");
+        covMS = mu.auxdata<AmgSymMatrix(5)>("MSCov");
+      }
+      else
+      {
+        parsCBNom = CB_track->definingParameters();
+        parsID = ID_track->definingParameters();
+        parsMS = ME_track->definingParameters();
+        covCBNom = CB_track->definingParametersCovMatrix();
+        covID = ID_track->definingParametersCovMatrix();
+        covMS = ME_track->definingParametersCovMatrix();
+      }
+
+
+      CorrectionCode NominalCorrCode=applyStatCombination(parsID, covID,
+                                                          parsMS, covMS,
+                                                          charge,
+                                                          parsCBNom,
+                                                          covCBNom,
+                                                          chi2Nom);
+      if(NominalCorrCode!=CorrectionCode::Ok) return NominalCorrCode;
+
+
+
       CorrectionCode idOK=applySagittaBiasCorrection(MCAST::SagittaCorType::ID, mu, itersID, false, isMC, muonInfo, SystCase);
       TLorentzVector lvIDCorr; lvIDCorr.SetPtEtaPhiM(muonInfo.ptid,muonInfo.eta,muonInfo.phi,mu.m()/1e3);
       if(idOK == CorrectionCode::Ok && tmpPtID!=0 ) tmpDeltaID = ( -tmpPtID +lvIDCorr.Pt() )/ tmpPtID  ;
       else tmpDeltaID=0;
       ATH_MSG_VERBOSE( "Shift ID "<<tmpDeltaID );
       // Now modify the ID covariance matrix
-      AmgVector(5) parsID = ID_track->definingParameters();
       parsID[4]=1.0 / (lvIDCorr.P()*1e3); 
       double stored_ptid = muonInfo.ptid;
 
@@ -755,75 +819,22 @@ namespace CP {
       else tmpDeltaMS=0;
       ATH_MSG_VERBOSE( "Shift MS "<<tmpDeltaMS );
       // Now modify the ME covariance matrix 
-      AmgVector(5) parsMS = ME_track->definingParameters();
       parsMS[4]=1.0 / (lvMECorr.P()*1e3);
       double stored_ptms = muonInfo.ptms;
 
 
-      // Calculate the stat combination before sagitta bias: 
-     
-      double chi2Nom=-999;
-      AmgVector(5) parsCBNom=CB_track->definingParameters();
-      AmgSymMatrix(5) covCBNom=CB_track->definingParametersCovMatrix();
-      int charge = mu.charge();
-      const ElementLink< xAOD::TrackParticleContainer >& ms_track = mu.extrapolatedMuonSpectrometerTrackParticleLink();
-      const ElementLink< xAOD::TrackParticleContainer >& id_track=mu.inDetTrackParticleLink();
-      CorrectionCode NominalCorrCode=applyStatCombination(id_track,
-                                                          ms_track,
-                                                          charge,
-                                                          parsCBNom,
-                                                          covCBNom,
-                                                          chi2Nom);
-      if(NominalCorrCode!=CorrectionCode::Ok) return NominalCorrCode;
+      CorrectionCode SysCorrCode=applyStatCombination(parsID, covID,
+                                                      parsMS, covMS,
+                                                      charge,
+                                                      parsCB,
+                                                      covCB,
+                                                      chi2Nom);
+      if(NominalCorrCode!=CorrectionCode::Ok) return SysCorrCode;
           
-      
-      // Perform the statistical combination 
-      AmgSymMatrix(5) covID = ID_track->definingParametersCovMatrix();
-      AmgSymMatrix(5) covMS = ME_track->definingParametersCovMatrix();
-      
-      const AmgSymMatrix(5)  weightID = covID.inverse();
-      if  ( weightID.determinant() == 0 ){
-        ATH_MSG_WARNING( " ID weight matrix computation failed     " ) ;
-        return CorrectionCode::Error;
-      }
-      
-      const AmgSymMatrix(5)  weightMS = covMS.inverse();
-      if  ( weightMS.determinant() == 0 ){
-        ATH_MSG_WARNING( "weightMS computation failed      " ) ;
-        return CorrectionCode::Error;
-      }
-      
-      AmgSymMatrix(5) weightCB = weightID + weightMS ;
-      AmgSymMatrix(5) covCB = weightCB.inverse();
-      if (covCB.determinant() == 0){
-        ATH_MSG_WARNING( " Inversion of weightCB failed " ) ;
-        return CorrectionCode::Error;
-      }
-
-      AmgSymMatrix(5) covSum = covID + covMS ;
-      AmgSymMatrix(5) invCovSum = covSum.inverse();
-      if (invCovSum.determinant() == 0){
-        ATH_MSG_WARNING( " Inversion of covSum failed " ) ;
-        return CorrectionCode::Error;
-      }
-      double  diffPhi = parsMS[2] - parsID[2] ;
-      if(diffPhi>M_PI)       parsMS[2] -= 2.*M_PI;
-      else if(diffPhi<-M_PI) parsMS[2] += 2.*M_PI;
-
-      //:: Chisquare calculation. Not used at this stage. To be decided if included as decoration in next round of reccomendations. 
-      //AmgVector(5) diffPars = parsID - parsMS;
-      //double chi2 = diffPars.transpose() * invCovSum * diffPars;
-      //chi2 = chi2/5. ;
-
-      AmgVector(5) parsCB = covCB * ( weightID * parsID + weightMS * parsMS ) ;
-      parsCB[4] *= muonInfo.charge;
-
-      if(parsCB[2]>M_PI)       parsCB[2] -= 2.*M_PI;
-      else if(parsCB[2]<-M_PI) parsCB[2] += 2.*M_PI;
-      double statCombPtNom = std::sin(parsCBNom[3])/std::abs(parsCBNom[4]) * MeVtoGeV;
+       double statCombPtNom = std::sin(parsCBNom[3])/std::abs(parsCBNom[4]) * MeVtoGeV;
       double statCombPt    = std::sin(parsCB[3])/std::abs(parsCB[4]) * MeVtoGeV;
       muonInfo.ptcb =  muonInfo.ptcb * (1  +  (statCombPt-statCombPtNom)/statCombPtNom ) ;
-      // ATH_MSG_VERBOSE(" Poor man's combination "<<simpleCombPt<<" Stat comb "<<statCombPt<<" Stat comb nom "<<" statCombPtNom "<<statCombPtNom ); 
+
       ATH_MSG_VERBOSE(" Stat comb "<<statCombPt<<" Stat comb nom "<<" statCombPtNom "<<statCombPtNom ); 
 
       muonInfo.ptid = stored_ptid;
@@ -1090,7 +1101,7 @@ namespace CP {
         // Rescaling
         muonInfo.ptcb = ptcb * ptcb_ratio; 
       }
-      else if((SystCase == MCAST::SagittaSysType::DATASTAT) || (SystCase == MCAST::SagittaSysType::BIAS))
+      else if((SystCase == MCAST::SagittaSysType::DATASTAT))
       {
         double nom_ptcb = nom_rho*ptCB + (1-nom_rho)*ptWeight;
         double org_ptcb = nom_rho*origPt + (1-nom_rho)*ptWeight_unCorr;
@@ -1123,57 +1134,75 @@ namespace CP {
 
     ATH_MSG_VERBOSE( "Muon Type = " << mu.muonType() << " ( 0: Combined, 1: StandAlone, 2: SegmentTagged, 3: CaloTagged )" );
     ATH_MSG_VERBOSE( "Muon Author = " << mu.author() );
-    ATH_MSG_VERBOSE( "Passes MST Selection = " << m_MuonSelectionTool->accept(mu));
 
     // Make InfoHelper for passing muon info between functions
     InfoHelper muonInfo;
+    if(m_expertMode)
+    {
+      static const SG::AuxElement::Accessor<float> id_pt("expert_ptid");
+      static const SG::AuxElement::Accessor<float> ms_pt("expert_ptms");
+      static const SG::AuxElement::Accessor<float> cb_pt("expert_ptcb");
 
-    // Set pt ID:
-    if( ( mu.inDetTrackParticleLink() ).isValid() ) {
-            const ElementLink< xAOD::TrackParticleContainer >& id_track = mu.inDetTrackParticleLink();
-      muonInfo.ptid = ( !id_track ) ? 0. : ( *id_track )->pt() * MeVtoGeV;
+      muonInfo.charge = mu.charge();
+      muonInfo.ptid = id_pt(mu) * MeVtoGeV;
       muonInfo.uncorrected_ptid = muonInfo.ptid;
-    } else {
-      ATH_MSG_VERBOSE("The ID track link is not valid - setting pT(ID)=0");
-      muonInfo.ptid = 0.;
-      muonInfo.uncorrected_ptid = 0.;
-    }
 
-    // Set pt MS:
-    if( ( mu.extrapolatedMuonSpectrometerTrackParticleLink() ).isValid() ) {
-          const ElementLink< xAOD::TrackParticleContainer >& ms_track = mu.extrapolatedMuonSpectrometerTrackParticleLink();
-      muonInfo.ptms = ( !ms_track ) ? 0. : ( *ms_track )->pt() * MeVtoGeV;
+      muonInfo.ptms = ms_pt(mu) * MeVtoGeV;
       muonInfo.uncorrected_ptms = muonInfo.ptms;
-    }
-    else{
-      ATH_MSG_VERBOSE("No link to extrapolatedMuonSpectrometerTrackParticleLink setting pT(MS)=0");
-      muonInfo.ptms = 0.;
-      muonInfo.uncorrected_ptms = 0.;
-    }
 
-    // Set pt CB:
-    if( ( mu.primaryTrackParticleLink() ).isValid() ) {
-      const ElementLink< xAOD::TrackParticleContainer >& cb_track = mu.primaryTrackParticleLink();
-      muonInfo.ptcb = ( !cb_track ) ? 0. : ( *cb_track )->pt() * MeVtoGeV;
+      muonInfo.ptcb = cb_pt(mu) * MeVtoGeV;
       muonInfo.uncorrected_ptcb = muonInfo.ptcb;
+
     }
-    else{
-      ATH_MSG_VERBOSE("The ID+MS combined track link is not valid - setting pT(ID)=0");
-      muonInfo.ptcb = 0.;
-      muonInfo.uncorrected_ptcb = 0;
-    }
+    else
+    {
+      // Set pt ID:
+      if( ( mu.inDetTrackParticleLink() ).isValid() ) {
+              const ElementLink< xAOD::TrackParticleContainer >& id_track = mu.inDetTrackParticleLink();
+        muonInfo.ptid = ( !id_track ) ? 0. : ( *id_track )->pt() * MeVtoGeV;
+        muonInfo.uncorrected_ptid = muonInfo.ptid;
+      } else {
+        ATH_MSG_VERBOSE("The ID track link is not valid - setting pT(ID)=0");
+        muonInfo.ptid = 0.;
+        muonInfo.uncorrected_ptid = 0.;
+      }
 
-    // Set the charge
-    if( mu.muonType() == xAOD::Muon::SiliconAssociatedForwardMuon ) {
-      if( mu.isAvailable<ElementLink<xAOD::TrackParticleContainer > > ("combinedTrackParticleLink")){
-        if(mu.combinedTrackParticleLink().isValid()){
-          const ElementLink< xAOD::TrackParticleContainer >& cb_track = mu.combinedTrackParticleLink();
-          muonInfo.charge = ( *cb_track )->charge();
+      // Set pt MS:
+      if( ( mu.extrapolatedMuonSpectrometerTrackParticleLink() ).isValid() ) {
+            const ElementLink< xAOD::TrackParticleContainer >& ms_track = mu.extrapolatedMuonSpectrometerTrackParticleLink();
+        muonInfo.ptms = ( !ms_track ) ? 0. : ( *ms_track )->pt() * MeVtoGeV;
+        muonInfo.uncorrected_ptms = muonInfo.ptms;
+      }
+      else{
+        ATH_MSG_VERBOSE("No link to extrapolatedMuonSpectrometerTrackParticleLink setting pT(MS)=0");
+        muonInfo.ptms = 0.;
+        muonInfo.uncorrected_ptms = 0.;
+      }
+
+      // Set pt CB:
+      if( ( mu.primaryTrackParticleLink() ).isValid() ) {
+        const ElementLink< xAOD::TrackParticleContainer >& cb_track = mu.primaryTrackParticleLink();
+        muonInfo.ptcb = ( !cb_track ) ? 0. : ( *cb_track )->pt() * MeVtoGeV;
+        muonInfo.uncorrected_ptcb = muonInfo.ptcb;
+      }
+      else{
+        ATH_MSG_VERBOSE("The ID+MS combined track link is not valid - setting pT(ID)=0");
+        muonInfo.ptcb = 0.;
+        muonInfo.uncorrected_ptcb = 0;
+      }
+
+      // Set the charge
+      if( mu.muonType() == xAOD::Muon::SiliconAssociatedForwardMuon ) {
+        if( mu.isAvailable<ElementLink<xAOD::TrackParticleContainer > > ("combinedTrackParticleLink")){
+          if(mu.combinedTrackParticleLink().isValid()){
+            const ElementLink< xAOD::TrackParticleContainer >& cb_track = mu.combinedTrackParticleLink();
+            muonInfo.charge = ( *cb_track )->charge();
+          }
         }
       }
-    }
-    else{
-      muonInfo.charge = mu.charge();
+      else{
+        muonInfo.charge = mu.charge();
+      }
     }
    
     muonInfo.eta = mu.eta();
@@ -1221,16 +1250,22 @@ namespace CP {
       }
     }
 
-    // Retrieve the event information:
-    const xAOD::EventInfo* evtInfo = 0;
-    if( evtStore()->retrieve( evtInfo, "EventInfo" ).isFailure() ) {
-      ATH_MSG_ERROR( "No EventInfo object could be retrieved" );
-      ATH_MSG_ERROR( "Random number generation not configured correctly, impossible to determine if dealing with data or MC" );
-      return CorrectionCode::Error;
+    bool isData = false;
+    if(m_expertMode) isData = m_expertMode_isData;
+    else
+    {
+      // Retrieve the event information:
+      const xAOD::EventInfo* evtInfo = 0;
+      if( evtStore()->retrieve( evtInfo, "EventInfo" ).isFailure() ) {
+        ATH_MSG_ERROR( "No EventInfo object could be retrieved" );
+        ATH_MSG_ERROR( "Random number generation not configured correctly, impossible to determine if dealing with data or MC" );
+        return CorrectionCode::Error;
+      }
+      isData = !evtInfo->eventType( xAOD::EventInfo::IS_SIMULATION );
     }
-    ATH_MSG_VERBOSE( "Checking Simulation flag: " << evtInfo->eventType( xAOD::EventInfo::IS_SIMULATION ) );
-
-    if( !evtInfo->eventType( xAOD::EventInfo::IS_SIMULATION ) ) {
+    ATH_MSG_VERBOSE( "Checking Simulation flag: " << !isData);
+    if( isData ) {
+      ATH_MSG_VERBOSE( "Doing data corrections");
 
       // Statistical combiantion specifics
       if(m_useStatComb){
@@ -1373,16 +1408,9 @@ namespace CP {
       mu.setP4( res_idPt, muonInfo.eta, muonInfo.phi );
     }
     else 
-    {
+    { 
       mu.auxdata< float >( "MuonSpectrometerPt" ) = res_msPt;
-      double cbPT=0;
-
-      if( ( mu.primaryTrackParticleLink() ).isValid() ) 
-      {
-        const ElementLink< xAOD::TrackParticleContainer >& cb_track = mu.primaryTrackParticleLink();
-        cbPT=(*cb_track)->pt() * MeVtoGeV;
-      }
-      ATH_MSG_VERBOSE("CB pt "<<cbPT<<" stat comb "<<muonInfo.ptcb<<" corrected pt "<<res_cbPt * MeVtoGeV);
+      ATH_MSG_VERBOSE("CB stat comb "<<muonInfo.ptcb<<" corrected pt "<<res_cbPt * MeVtoGeV);
       mu.setP4( res_cbPt, muonInfo.eta, muonInfo.phi );
     }
 
@@ -1501,11 +1529,14 @@ namespace CP {
     if( !m_useExternalSeed ) {
       // Get Event Number, Retrieve the event information:
       const xAOD::EventInfo* evtInfo = 0;
-      ATH_MSG_VERBOSE( "Retrieving EventInfo from the EventStore..." );
-      if( evtStore()->retrieve( evtInfo, "EventInfo" ).isFailure() ) {
-        ATH_MSG_ERROR( "No EventInfo object could be retrieved" );
-        ATH_MSG_ERROR( "Random number generation not configured correctly, impossible to determine if dealing with data or MC" );
-        return CorrectionCode::Error;
+      if(!m_expertMode) // expert mode is running for validation where this information is not avaliable
+      {
+        ATH_MSG_VERBOSE( "Retrieving EventInfo from the EventStore..." );
+        if( evtStore()->retrieve( evtInfo, "EventInfo" ).isFailure() ) {
+          ATH_MSG_ERROR( "No EventInfo object could be retrieved" );
+          ATH_MSG_ERROR( "Random number generation not configured correctly, impossible to determine if dealing with data or MC" );
+          return CorrectionCode::Error;
+        }
       }
 
       const unsigned long long eventNumber = evtInfo ? evtInfo->eventNumber() : 0;
@@ -2707,39 +2738,52 @@ namespace CP {
   double MuonCalibrationAndSmearingTool::ExpectedResolution( const int DetType,const xAOD::Muon& mu, const bool mc ) const {
 
     // get the pt measurements from the xAOD::Muon object
-
     double loc_ptid = 0;
     double loc_ptms = 0;
     double loc_ptcb = 0;
 
-    // Set pt ID:
-    if( ( mu.inDetTrackParticleLink() ).isValid() ) {
-        const ElementLink< xAOD::TrackParticleContainer >& id_track = mu.inDetTrackParticleLink();
-        loc_ptid = ( !id_track ) ? 0. : ( *id_track )->pt() * MeVtoGeV;
-    }
-    else{
-      ATH_MSG_VERBOSE("The ID track link is not valid - setting pT(ID)=0");
-      loc_ptid = 0.;
-    }
+    if(m_expertMode)
+    {
+      static const SG::AuxElement::Accessor<float> id_pt("expert_ptid");
+      static const SG::AuxElement::Accessor<float> ms_pt("expert_ptms");
+      static const SG::AuxElement::Accessor<float> cb_pt("expert_ptcb");
 
-    // Set pt MS:
-    if( ( mu.extrapolatedMuonSpectrometerTrackParticleLink() ).isValid() ) {
-        const ElementLink< xAOD::TrackParticleContainer >& ms_track = mu.extrapolatedMuonSpectrometerTrackParticleLink();
-        loc_ptms = ( !ms_track ) ? 0. : ( *ms_track )->pt() * MeVtoGeV;
-    }
-    else{
-      ATH_MSG_VERBOSE("No link to extrapolatedMuonSpectrometerTrackParticleLink setting pT(MS)=0");
-      loc_ptms = 0.;
-    }
+      loc_ptid = id_pt(mu) * MeVtoGeV;
+      loc_ptms = ms_pt(mu) * MeVtoGeV;
+      loc_ptcb = cb_pt(mu) * MeVtoGeV;
 
-    // Set pt CB;
-    if( ( mu.primaryTrackParticleLink() ).isValid() ) {
-        const ElementLink< xAOD::TrackParticleContainer >& cb_track = mu.primaryTrackParticleLink();
-        loc_ptcb = ( !cb_track ) ? 0. : ( *cb_track )->pt() * MeVtoGeV;
     }
-    else{
-      ATH_MSG_VERBOSE("The ID+MS combined track link is not valid - setting pT(ID)=0");
-      loc_ptcb = 0.;
+    else
+    {
+      // Set pt ID:
+      if( ( mu.inDetTrackParticleLink() ).isValid() ) {
+          const ElementLink< xAOD::TrackParticleContainer >& id_track = mu.inDetTrackParticleLink();
+          loc_ptid = ( !id_track ) ? 0. : ( *id_track )->pt() * MeVtoGeV;
+      }
+      else{
+        ATH_MSG_VERBOSE("The ID track link is not valid - setting pT(ID)=0");
+        loc_ptid = 0.;
+      }
+
+      // Set pt MS:
+      if( ( mu.extrapolatedMuonSpectrometerTrackParticleLink() ).isValid() ) {
+          const ElementLink< xAOD::TrackParticleContainer >& ms_track = mu.extrapolatedMuonSpectrometerTrackParticleLink();
+          loc_ptms = ( !ms_track ) ? 0. : ( *ms_track )->pt() * MeVtoGeV;
+      }
+      else{
+        ATH_MSG_VERBOSE("No link to extrapolatedMuonSpectrometerTrackParticleLink setting pT(MS)=0");
+        loc_ptms = 0.;
+      }
+
+      // Set pt CB;
+      if( ( mu.primaryTrackParticleLink() ).isValid() ) {
+          const ElementLink< xAOD::TrackParticleContainer >& cb_track = mu.primaryTrackParticleLink();
+          loc_ptcb = ( !cb_track ) ? 0. : ( *cb_track )->pt() * MeVtoGeV;
+      }
+      else{
+        ATH_MSG_VERBOSE("The ID+MS combined track link is not valid - setting pT(ID)=0");
+        loc_ptcb = 0.;
+      }
     }
 
     // obtain detRegion given the muon (eta,phi)
@@ -3386,12 +3430,37 @@ namespace CP {
       }}
 
     const xAOD::TrackParticle*       cbtrack = mu.trackParticle( xAOD::Muon::CombinedTrackParticle );
-    AmgVector(5) parsCB=cbtrack->definingParameters();
-    AmgSymMatrix(5) covCB=cbtrack->definingParametersCovMatrix();
+    AmgVector(5) parsCB ;
+    AmgSymMatrix(5) covCB;
     double chi2=-999;
 
+    AmgVector(5) parsID;
+    AmgVector(5) parsMS;
+    AmgSymMatrix(5) covID;
+    AmgSymMatrix(5) covMS;
 
-    if(applyStatCombination(id_track,ms_track,charge,parsCB,covCB,chi2)!=CorrectionCode::Ok)
+    if(m_expertMode)
+    {
+      parsCB = mu.auxdata<AmgVector(5)>("CBParam");
+      parsID = mu.auxdata<AmgVector(5)>("IDParam");
+      parsMS = mu.auxdata<AmgVector(5)>("MSParam");
+      covCB = mu.auxdata<AmgSymMatrix(5)>("CBCov");
+      covID = mu.auxdata<AmgSymMatrix(5)>("IDCov");
+      covMS = mu.auxdata<AmgSymMatrix(5)>("MSCov");
+    }
+    else
+    {
+      parsCB = cbtrack->definingParameters();
+      parsID = (*id_track)->definingParameters();
+      parsMS = (*ms_track)->definingParameters();
+      covCB = cbtrack->definingParametersCovMatrix();
+      covID = (*id_track)->definingParametersCovMatrix();
+      covMS = (*ms_track)->definingParametersCovMatrix();
+    }
+
+
+
+    if(applyStatCombination(parsID, covID, parsMS, covMS,charge,parsCB,covCB,chi2)!=CorrectionCode::Ok)
       return CorrectionCode::Error;
 
     muonInfo.cbParsA.clear();
@@ -3417,21 +3486,17 @@ namespace CP {
     return CorrectionCode::Ok;
   }
 
-  CorrectionCode MuonCalibrationAndSmearingTool::applyStatCombination( const ElementLink< xAOD::TrackParticleContainer >& inDetTrackParticle,
-                                                                       const ElementLink< xAOD::TrackParticleContainer >& extrTrackParticle,
+  CorrectionCode MuonCalibrationAndSmearingTool::applyStatCombination( AmgVector(5) parsID, AmgSymMatrix(5) covID,
+                                                                       AmgVector(5) parsMS, AmgSymMatrix(5) covMS,
                                                                        int charge,
                                                                        AmgVector(5)& parsCB,
                                                                        AmgSymMatrix(5)& covCB,
                                                                        double& chi2) const {
 
     chi2   = 1e20;
-    AmgVector(5) parsID = (*inDetTrackParticle)->definingParameters();;
     parsID[4] = std::abs(parsID[4]);
-
-    AmgVector(5) parsMS = (*extrTrackParticle)->definingParameters();;
     parsMS[4] = std::abs(parsMS[4]);
 
-    AmgSymMatrix(5) covID = (*inDetTrackParticle)->definingParametersCovMatrix();
 
     const AmgSymMatrix(5)  weightID = covID.inverse();
     if  ( weightID.determinant() == 0 ){
@@ -3439,23 +3504,6 @@ namespace CP {
       return CorrectionCode::Error;
     }
 
-    AmgSymMatrix(5) covMS = (*extrTrackParticle)->definingParametersCovMatrix();
-
-    /*
-    // Error inflation to account for ID-MS misaligment
-    double dSigma[5] = { 1.0, 2.0/sin(parsMS[3]), 0.001, 0.001, 0.0};
-    double factor[5];
-    for ( int i = 0; i<5 ; i++ ) {
-      factor[i] = std::sqrt((covMS(i,i)+dSigma[i]*dSigma[i])/covMS(i,i));
-    }
-    for ( int i = 0 ; i<5 ; i++ ) {
-      for ( int j = 0 ; j<5 ; j++ ) {
-        if ( i==j ) {
-          covMS(i,j) = covMS(i,j)*factor[i]*factor[j];
-        }
-      }
-      }*/
-
     const AmgSymMatrix(5)  weightMS = covMS.inverse();
     if  ( weightMS.determinant() == 0 ){
       ATH_MSG_WARNING( "weightMS computation failed      " ) ;
diff --git a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/util/MCAST_Tester.cxx b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/util/MCAST_Tester.cxx
index 8a1170c137ca..9a0fa5755aca 100644
--- a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/util/MCAST_Tester.cxx
+++ b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/util/MCAST_Tester.cxx
@@ -151,7 +151,7 @@ int main( int argc, char* argv[] ) {
   corrTool.setProperty("noEigenDecor" ,         false);
   corrTool.setProperty("sagittaMapsInputType" ,         1);
   corrTool.setProperty("sagittaMapUnitConversion",1); 
-  corrTool.setProperty("systematicCorrelationScheme", "Decorr_Scale"); 
+  // corrTool.setProperty("systematicCorrelationScheme", "Decorr_Scale"); 
   //   corrTool.setProperty("useExternalSeed" ,      false);
   //   corrTool.setProperty("externalSeed" ,         0);
 
-- 
GitLab


From 7a173e8807bcd0da600ab3c1f3fc291e6d65b5c9 Mon Sep 17 00:00:00 2001
From: Syed Haider Abidi <syed.haider.abidi@cern.ch>
Date: Wed, 15 Sep 2021 17:20:01 +0200
Subject: [PATCH 08/14] final config for saggita setup and sys

---
 .../MuonCalibrationAndSmearingTool.h          |  1 +
 .../Root/MuonCalibrationAndSmearingTool.cxx   | 34 +++++++++++++------
 .../util/MCAST_Tester.cxx                     | 11 +++---
 3 files changed, 28 insertions(+), 18 deletions(-)

diff --git a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/MuonMomentumCorrections/MuonCalibrationAndSmearingTool.h b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/MuonMomentumCorrections/MuonCalibrationAndSmearingTool.h
index 7357b022a3b4..44d6c310e021 100644
--- a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/MuonMomentumCorrections/MuonCalibrationAndSmearingTool.h
+++ b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/MuonMomentumCorrections/MuonCalibrationAndSmearingTool.h
@@ -199,6 +199,7 @@ class MuonCalibrationAndSmearingTool : public virtual IMuonCalibrationAndSmearin
     bool m_2stations_highpt_smearing;
     bool m_extra_decorations;
     bool m_toroidOff;
+    double m_extraRebiasSys;
     int m_Tsmear;
     int m_Tdata;
     int m_Trel;
diff --git a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx
index 46649d24eaf9..911311028263 100644
--- a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx
+++ b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx
@@ -56,7 +56,7 @@ namespace CP {
     declareProperty("MinCombPt", m_StatCombPtThreshold=300.0);
     declareProperty("HighPtSystThr", m_HighPtSystThreshold=300.0);
     declareProperty("SagittaCorr", m_doSagittaCorrection = false);
-    declareProperty("SagittaRelease", m_SagittaRelease = "sagittaBiasDataAll_03_02_19");
+    declareProperty("SagittaRelease", m_SagittaRelease = "sagittaBiasDataAll_15_09_2021");
     declareProperty("doSagittaMCDistortion",m_doSagittaMCDistortion=true);
     declareProperty("SagittaCorrPhaseSpace",m_SagittaCorrPhaseSpace=false);
     declareProperty("SagittaIterWeight", m_IterWeight=0.5);
@@ -64,16 +64,17 @@ namespace CP {
     declareProperty("sgItersID",m_sgItersID=4);
     declareProperty("sgItersME",m_sgItersME=4);
     declareProperty("sgIetrsManual",m_sgIetrsManual=false);
-    declareProperty("fixedRho",m_fixedRho=1.0);
-    declareProperty("useFixedRho",m_useFixedRho=false);
+    declareProperty("fixedRho", m_fixedRho=1.0);
+    declareProperty("useFixedRho",m_useFixedRho=true);
     declareProperty("noEigenDecor" ,m_doNotUseAMGMATRIXDECOR=false);
     declareProperty("useExternalSeed" ,m_useExternalSeed=false);
     declareProperty("externalSeed" ,m_externalSeed=0);
-    declareProperty("sagittaMapsInputType", m_saggitaMapsInputType=MCAST::SagittaInputHistType::NOMINAL);
-    declareProperty("sagittaMapUnitConversion",m_sagittaMapUnitConversion=1e-3);
+    declareProperty("sagittaMapsInputType", m_saggitaMapsInputType=MCAST::SagittaInputHistType::SINGLE);
+    declareProperty("sagittaMapUnitConversion",m_sagittaMapUnitConversion=1);
     declareProperty("systematicCorrelationScheme", m_sysScheme = "Corr_Scale");
     declareProperty("expertMode", m_expertMode = false);
     declareProperty("expertSetting_isData", m_expertMode_isData = false);
+    declareProperty("expertSetting_extraRebiasSys", m_extraRebiasSys = 0.00003);
 
     m_SagittaIterations.push_back(0);
     m_SagittaIterations.push_back(0);
@@ -525,12 +526,22 @@ namespace CP {
     
     if(m_currentParameters->SagittaBias == MCAST::SystVariation::Up || m_currentParameters->SagittaDataStat == MCAST::SystVariation::Up)
     {
-      p2 = 0.5*p2;
+      p2 = p2;
     }
     else if (m_currentParameters->SagittaBias == MCAST::SystVariation::Down || m_currentParameters->SagittaDataStat == MCAST::SystVariation::Down)
     {
-      p2 = -0.5*p2;
+      p2 = -p2;
     }
+
+    if(m_currentParameters->SagittaBias == MCAST::SystVariation::Up)
+    {
+      p2 += m_extraRebiasSys;
+    }
+    else if (m_currentParameters->SagittaBias == MCAST::SystVariation::Down)
+    {
+      p2 -= m_extraRebiasSys;
+    }
+
     
     return p2;
   }
@@ -1678,10 +1689,11 @@ namespace CP {
       result.insert( SystematicVariation( "MUON_SCALE_MS_ELOSS", -1 ) );
     }
     // Sagitta correction rho
-    //if(!m_useFixedRho){
-    result.insert( SystematicVariation( "MUON_SAGITTA_RHO", 1 ) );
-    result.insert( SystematicVariation( "MUON_SAGITTA_RHO", -1 ) );
-    //}
+    if(!m_useFixedRho)
+    {
+      result.insert( SystematicVariation( "MUON_SAGITTA_RHO", 1 ) );
+      result.insert( SystematicVariation( "MUON_SAGITTA_RHO", -1 ) );
+    }
 
     // Sagitta correction resid bias
     result.insert( SystematicVariation( "MUON_SAGITTA_RESBIAS", 1 ) );
diff --git a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/util/MCAST_Tester.cxx b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/util/MCAST_Tester.cxx
index 9a0fa5755aca..20fa7d7384e3 100644
--- a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/util/MCAST_Tester.cxx
+++ b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/util/MCAST_Tester.cxx
@@ -131,14 +131,11 @@ int main( int argc, char* argv[] ) {
   //   corrTool.setProperty("Algo",                  "muons" );
   //   corrTool.setProperty("SmearingType",          "q_pT" );
   corrTool.setProperty("Release",               "Recs2021_04_18_Prelim" );
-  //corrTool.setProperty("Release",               "Recs2018_05_20" );     
-  //   corrTool.setProperty("ToroidOff",             false );
-  //   corrTool.setProperty("FilesPath",             "" );
+
   corrTool.setProperty("StatComb",              false);
   //   corrTool.setProperty("MinCombPt",             300.0);
   corrTool.setProperty("SagittaCorr",           true);
-  //corrTool.setProperty("SagittaRelease",        "sagittaBiasDataAll_30_07_18");
-  corrTool.setProperty("SagittaRelease",        "sagittaBiasDataAll_SingleHistMethod_B_17_06_2021_v2");
+  corrTool.setProperty("SagittaRelease",        "sagittaBiasDataAll_15_09_2021");
   corrTool.setProperty("doSagittaMCDistortion", false);
   corrTool.setProperty("SagittaCorrPhaseSpace", false);
   corrTool.setProperty("SagittaIterWeight", 1);
@@ -146,8 +143,8 @@ int main( int argc, char* argv[] ) {
   //   corrTool.setProperty("sgItersID",             11);
   //   corrTool.setProperty("sgItersME",             11);
   //   corrTool.setProperty("sgIetrsMamual",         false);
-  corrTool.setProperty("fixedRho",              0.0);
-  corrTool.setProperty("useFixedRho",           false);
+  corrTool.setProperty("fixedRho",              1.0);
+  corrTool.setProperty("useFixedRho",           true);
   corrTool.setProperty("noEigenDecor" ,         false);
   corrTool.setProperty("sagittaMapsInputType" ,         1);
   corrTool.setProperty("sagittaMapUnitConversion",1); 
-- 
GitLab


From e4a30b922948db2ad55d5f34a64e3c879dd3fe28 Mon Sep 17 00:00:00 2001
From: Haider Abidi <syed.haider.abidi@cern.ch>
Date: Wed, 6 Oct 2021 18:02:02 +0200
Subject: [PATCH 09/14] typo fix

---
 .../Root/MuonCalibrationAndSmearingTool.cxx                   | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx
index 911311028263..f75d92b5b146 100644
--- a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx
+++ b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx
@@ -840,9 +840,9 @@ namespace CP {
                                                       parsCB,
                                                       covCB,
                                                       chi2Nom);
-      if(NominalCorrCode!=CorrectionCode::Ok) return SysCorrCode;
+      if(SysCorrCode!=CorrectionCode::Ok) return SysCorrCode;
           
-       double statCombPtNom = std::sin(parsCBNom[3])/std::abs(parsCBNom[4]) * MeVtoGeV;
+      double statCombPtNom = std::sin(parsCBNom[3])/std::abs(parsCBNom[4]) * MeVtoGeV;
       double statCombPt    = std::sin(parsCB[3])/std::abs(parsCB[4]) * MeVtoGeV;
       muonInfo.ptcb =  muonInfo.ptcb * (1  +  (statCombPt-statCombPtNom)/statCombPtNom ) ;
 
-- 
GitLab


From 9f6bb0b90d293f7a9b17700b0a58b259193e138f Mon Sep 17 00:00:00 2001
From: Syed Haider Abidi <syed.haider.abidi@cern.ch>
Date: Mon, 8 Nov 2021 11:14:20 +0100
Subject: [PATCH 10/14] support for new calibration

---
 .../MuonCalibrationAndSmearingTool.h          |  16 +-
 .../Root/MuonCalibrationAndSmearingTool.cxx   | 278 ++++++++++++++----
 .../util/MCAST_Tester.cxx                     |   3 +-
 3 files changed, 240 insertions(+), 57 deletions(-)

diff --git a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/MuonMomentumCorrections/MuonCalibrationAndSmearingTool.h b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/MuonMomentumCorrections/MuonCalibrationAndSmearingTool.h
index 44d6c310e021..1afc9c72864f 100644
--- a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/MuonMomentumCorrections/MuonCalibrationAndSmearingTool.h
+++ b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/MuonMomentumCorrections/MuonCalibrationAndSmearingTool.h
@@ -99,6 +99,7 @@ class MuonCalibrationAndSmearingTool : public virtual IMuonCalibrationAndSmearin
       double smearDeltaID = 0;
       double smearDeltaCB = 0;
       double smearDeltaCBOnly = 0;
+      double smearDeltaCBDirect = 0;
       int    sel_category = -1;
       double uncorrected_ptcb = 0;
       double uncorrected_ptid = 0;
@@ -145,8 +146,8 @@ class MuonCalibrationAndSmearingTool : public virtual IMuonCalibrationAndSmearin
     float        GetRegionInnerEta( const int r_i ) const; //Return Eta closer to the origin
     std::string  GetRegionName( const int r_i ) const;
     std::string  GetRegionName( const double eta, const double phi ) const;
-    double GetSmearing( int DetType, InfoHelper& muonInfo ) const;
-    double GetSystVariation( int DetType, double var, InfoHelper& muonInfo ) const;
+    double GetSmearing( int DetType, InfoHelper& muonInfo, bool doDirectCB ) const;
+    double GetSystVariation( int DetType, double var, InfoHelper& muonInfo, bool doDirectCB ) const;
     StatusCode SetInfoHelperCorConsts(InfoHelper& inMuonInfo) const;
     void CalcCBWeights( xAOD::Muon&, InfoHelper& muonInfo ) const;
     double CalculatePt( const int DetType, const double inSmearID, const double inSmearMS, const double scaleVarID, const double scaleMS_scale, const double scaleMS_egLoss, InfoHelper& muonInfo ) const;
@@ -199,12 +200,13 @@ class MuonCalibrationAndSmearingTool : public virtual IMuonCalibrationAndSmearin
     bool m_2stations_highpt_smearing;
     bool m_extra_decorations;
     bool m_toroidOff;
-    double m_extraRebiasSys;
     int m_Tsmear;
     int m_Tdata;
     int m_Trel;
     int m_Talgo;
+    double m_extraRebiasSys;
     double m_useNsigmaForICombine;
+    bool m_doDirectCBCalib;
     std::vector<double> m_scale_ID, m_enLoss_MS, m_scale_MS, m_scale_CB;
 
     //sys variations (stat error added in quadrature), one if it's simmetrized, 2 if Up != Dw.
@@ -222,6 +224,14 @@ class MuonCalibrationAndSmearingTool : public virtual IMuonCalibrationAndSmearin
     std::vector<double> m_SUp_p1_ID, m_SUp_p2_ID, m_SUp_p2_ID_TAN, m_SUp_p0_MS, m_SUp_p1_MS, m_SUp_p2_MS;
     std::vector<double> m_SDw_p1_ID, m_SDw_p2_ID, m_SDw_p2_ID_TAN, m_SDw_p0_MS, m_SDw_p1_MS, m_SDw_p2_MS;
     std::vector<double> m_MC_p1_ID, m_MC_p2_ID, m_MC_p2_ID_TAN, m_MC_p0_MS, m_MC_p1_MS, m_MC_p2_MS;
+
+
+    std::vector<double> m_S_0_CB, m_SUp_0_CB, m_SDw_0_CB;
+    std::vector<double> m_S_1_CB, m_SUp_1_CB, m_SDw_1_CB;
+    std::vector<double> m_R_0_CB, m_RUp_0_CB, m_RDw_0_CB;
+    std::vector<double> m_R_1_CB, m_RUp_1_CB, m_RDw_1_CB;
+    std::vector<double> m_R_2_CB, m_RUp_2_CB, m_RDw_2_CB;
+
     // Special "p2" systematics and corrections for non-three-station muons
     // Maps have two keys: detector region and category
     std::map<std::pair<int, int>, std::pair<double, double> > m_extra_p1_p2_MS_AlignedOnly, m_extra_p1_p2_MS_AlignedAndCorrected, m_extra_p1_p2_MS_Misaligned;
diff --git a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx
index 911311028263..cb4be124157f 100644
--- a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx
+++ b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx
@@ -39,14 +39,14 @@ namespace CP {
     m_doSagittaCorrection(false),
     m_doSagittaMCDistortion(false),
     m_IterWeight(0.5),
-    m_SagittaRelease("sagittaBiasDataAll_03_02_19"),
+    m_SagittaRelease("sagittaBiasDataAll_15_09_2021"),
     
     m_MuonSelectionTool("") {
     
     declareProperty("Year", m_year = "Data16" );
     declareProperty("Algo", m_algo = "muons" );
     declareProperty("SmearingType", m_type = "q_pT" );
-    declareProperty("Release", m_release = "Recs2020_03_03" );
+    declareProperty("Release", m_release = "Recs2021_11_04" );
     declareProperty("ToroidOff", m_toroidOff = false );
     declareProperty("AddExtraDecorations", m_extra_decorations = false );
     declareProperty("doExtraSmearing", m_extra_highpt_smearing = false );
@@ -55,16 +55,16 @@ namespace CP {
     declareProperty("StatComb", m_useStatComb = false);
     declareProperty("MinCombPt", m_StatCombPtThreshold=300.0);
     declareProperty("HighPtSystThr", m_HighPtSystThreshold=300.0);
-    declareProperty("SagittaCorr", m_doSagittaCorrection = false);
+    declareProperty("SagittaCorr", m_doSagittaCorrection = true);
     declareProperty("SagittaRelease", m_SagittaRelease = "sagittaBiasDataAll_15_09_2021");
-    declareProperty("doSagittaMCDistortion",m_doSagittaMCDistortion=true);
+    declareProperty("doSagittaMCDistortion",m_doSagittaMCDistortion=false);
     declareProperty("SagittaCorrPhaseSpace",m_SagittaCorrPhaseSpace=false);
-    declareProperty("SagittaIterWeight", m_IterWeight=0.5);
+    declareProperty("SagittaIterWeight", m_IterWeight=1.0);
     declareProperty("sgItersCB",m_sgItersCB=4);
     declareProperty("sgItersID",m_sgItersID=4);
     declareProperty("sgItersME",m_sgItersME=4);
     declareProperty("sgIetrsManual",m_sgIetrsManual=false);
-    declareProperty("fixedRho", m_fixedRho=1.0);
+    declareProperty("fixedRho",m_fixedRho=1.0);
     declareProperty("useFixedRho",m_useFixedRho=true);
     declareProperty("noEigenDecor" ,m_doNotUseAMGMATRIXDECOR=false);
     declareProperty("useExternalSeed" ,m_useExternalSeed=false);
@@ -74,14 +74,15 @@ namespace CP {
     declareProperty("systematicCorrelationScheme", m_sysScheme = "Corr_Scale");
     declareProperty("expertMode", m_expertMode = false);
     declareProperty("expertSetting_isData", m_expertMode_isData = false);
-    declareProperty("expertSetting_extraRebiasSys", m_extraRebiasSys = 0.00003);
+    declareProperty("expert_extraRebiasSys", m_extraRebiasSys = 0.00003);
+    declareProperty("expert_doDirectCBCalib", m_doDirectCBCalib = false);
 
     m_SagittaIterations.push_back(0);
     m_SagittaIterations.push_back(0);
     m_SagittaIterations.push_back(0);
-    m_sagittaPhaseSpaceCB=nullptr;
-    m_sagittaPhaseSpaceID=nullptr;
-    m_sagittaPhaseSpaceME=nullptr;
+    m_sagittaPhaseSpaceCB = nullptr;
+    m_sagittaPhaseSpaceID = nullptr;
+    m_sagittaPhaseSpaceME = nullptr;
   }
 
   MuonCalibrationAndSmearingTool::MuonCalibrationAndSmearingTool( const MuonCalibrationAndSmearingTool& tool ) :
@@ -230,6 +231,7 @@ namespace CP {
     ATH_CHECK(m_MuonSelectionTool.setProperty("MaxEta", 2.7));
     ATH_CHECK(m_MuonSelectionTool.setProperty("MuQuality", 1));
     ATH_CHECK(m_MuonSelectionTool.setProperty("TurnOffMomCorr", true));
+    ATH_CHECK(m_MuonSelectionTool.setProperty("OutputLevel", MSG::FATAL));
     ATH_CHECK(m_MuonSelectionTool.retrieve());
 
     std::string regionsPath;
@@ -526,11 +528,11 @@ namespace CP {
     
     if(m_currentParameters->SagittaBias == MCAST::SystVariation::Up || m_currentParameters->SagittaDataStat == MCAST::SystVariation::Up)
     {
-      p2 = p2;
+      p2 = 0.5*p2;
     }
     else if (m_currentParameters->SagittaBias == MCAST::SystVariation::Down || m_currentParameters->SagittaDataStat == MCAST::SystVariation::Down)
     {
-      p2 = -p2;
+      p2 = -0.5*p2;
     }
 
     if(m_currentParameters->SagittaBias == MCAST::SystVariation::Up)
@@ -745,10 +747,6 @@ namespace CP {
 
     else if( SgCorrType == MCAST::SagittaCorType::WEIGHTS)
     {
-      const xAOD::TrackParticle* CB_track  = mu.trackParticle( xAOD::Muon::CombinedTrackParticle );
-      const xAOD::TrackParticle* ID_track  = mu.trackParticle( xAOD::Muon::InnerDetectorTrackParticle );
-      const xAOD::TrackParticle* ME_track  = mu.trackParticle( xAOD::Muon::ExtrapolatedMuonSpectrometerTrackParticle );
-
       CalcCBWeights( mu, muonInfo );
 
       if( iter == m_sagittasCB.size() ) 
@@ -795,6 +793,9 @@ namespace CP {
       }
       else
       {
+        const xAOD::TrackParticle* CB_track  = mu.trackParticle( xAOD::Muon::CombinedTrackParticle );
+        const xAOD::TrackParticle* ID_track  = mu.trackParticle( xAOD::Muon::InnerDetectorTrackParticle );
+        const xAOD::TrackParticle* ME_track  = mu.trackParticle( xAOD::Muon::ExtrapolatedMuonSpectrometerTrackParticle );
         parsCBNom = CB_track->definingParameters();
         parsID = ID_track->definingParameters();
         parsMS = ME_track->definingParameters();
@@ -1600,29 +1601,34 @@ namespace CP {
     // Getting smearing values
     // MS
     if( m_currentParameters->SmearTypeMS == MCAST::SystVariation::Default ) {
-      inMuonInfo.smearDeltaMS = GetSmearing( MCAST::DetectorType::MS, inMuonInfo );
-      inMuonInfo.smearDeltaCBOnly = GetSmearing( MCAST::DetectorType::CB, inMuonInfo );
+      inMuonInfo.smearDeltaMS = GetSmearing( MCAST::DetectorType::MS, inMuonInfo, false );
+      inMuonInfo.smearDeltaCBOnly = GetSmearing( MCAST::DetectorType::CB, inMuonInfo, false );
+      if(m_doDirectCBCalib) inMuonInfo.smearDeltaCBDirect = GetSmearing( MCAST::DetectorType::CB, inMuonInfo, true );
     }
     else if( m_currentParameters->SmearTypeMS == MCAST::SystVariation::Up ) {
-      inMuonInfo.smearDeltaMS = GetSystVariation( MCAST::DetectorType::MS, 1., inMuonInfo );
-      inMuonInfo.smearDeltaCBOnly = GetSystVariation( MCAST::DetectorType::CB, 1., inMuonInfo );
+      inMuonInfo.smearDeltaMS = GetSystVariation( MCAST::DetectorType::MS, 1., inMuonInfo, false );
+      inMuonInfo.smearDeltaCBOnly = GetSystVariation( MCAST::DetectorType::CB, 1., inMuonInfo, false );
+      if(m_doDirectCBCalib) inMuonInfo.smearDeltaCBDirect = GetSystVariation( MCAST::DetectorType::CB, 1, inMuonInfo, true );
+
+
     }
     else if( m_currentParameters->SmearTypeMS == MCAST::SystVariation::Down ) {
-      inMuonInfo.smearDeltaMS = GetSystVariation( MCAST::DetectorType::MS, -1., inMuonInfo );
-      inMuonInfo.smearDeltaCBOnly = GetSystVariation( MCAST::DetectorType::CB, -1., inMuonInfo );
+      inMuonInfo.smearDeltaMS = GetSystVariation( MCAST::DetectorType::MS, -1., inMuonInfo, false );
+      inMuonInfo.smearDeltaCBOnly = GetSystVariation( MCAST::DetectorType::CB, -1., inMuonInfo, false );
+      if(m_doDirectCBCalib) inMuonInfo.smearDeltaCBDirect = GetSystVariation( MCAST::DetectorType::CB, -1, inMuonInfo, true );
     }
     else {
       ATH_MSG_ERROR( "Invalid value for m_currentParameters->SmearTypeMS" );
     }
     // ID
     if( m_currentParameters->SmearTypeID == MCAST::SystVariation::Default ) {
-      inMuonInfo.smearDeltaID = GetSmearing( MCAST::DetectorType::ID, inMuonInfo );
+      inMuonInfo.smearDeltaID = GetSmearing( MCAST::DetectorType::ID, inMuonInfo, false );
     }
     else if( m_currentParameters->SmearTypeID == MCAST::SystVariation::Up ) {
-      inMuonInfo.smearDeltaID = GetSystVariation( MCAST::DetectorType::ID, 1., inMuonInfo );
+      inMuonInfo.smearDeltaID = GetSystVariation( MCAST::DetectorType::ID, 1., inMuonInfo, false );
     }
     else if( m_currentParameters->SmearTypeID == MCAST::SystVariation::Down ) {
-      inMuonInfo.smearDeltaID = GetSystVariation( MCAST::DetectorType::ID, -1., inMuonInfo );
+      inMuonInfo.smearDeltaID = GetSystVariation( MCAST::DetectorType::ID, -1., inMuonInfo, false );
     }
     else {
       ATH_MSG_ERROR( "Invalid value for m_currentParameters->SmearTypeID" );
@@ -1689,11 +1695,10 @@ namespace CP {
       result.insert( SystematicVariation( "MUON_SCALE_MS_ELOSS", -1 ) );
     }
     // Sagitta correction rho
-    if(!m_useFixedRho)
-    {
-      result.insert( SystematicVariation( "MUON_SAGITTA_RHO", 1 ) );
-      result.insert( SystematicVariation( "MUON_SAGITTA_RHO", -1 ) );
-    }
+    //if(!m_useFixedRho){
+    result.insert( SystematicVariation( "MUON_SAGITTA_RHO", 1 ) );
+    result.insert( SystematicVariation( "MUON_SAGITTA_RHO", -1 ) );
+    //}
 
     // Sagitta correction resid bias
     result.insert( SystematicVariation( "MUON_SAGITTA_RESBIAS", 1 ) );
@@ -2111,6 +2116,12 @@ namespace CP {
     else if (rel == "Recs2021_04_18_Prelim") {
         m_Trel = MCAST::Release::Recs2021_07_01;
     }
+    else if (rel == "Recs2021_07_01") {
+        m_Trel = MCAST::Release::Recs2021_07_01;
+    }
+    else if (rel == "Recs2021_11_04") {
+        m_Trel = MCAST::Release::Recs2021_07_01;
+    }
     else 
     {
         m_Trel = MCAST::Release::Recs2021_07_01;
@@ -2135,7 +2146,7 @@ namespace CP {
 
   double MuonCalibrationAndSmearingTool::CalculatePt( const int DetType, const double inSmearID, const double inSmearMS, const double scaleVarID, const double scaleMS_scale, const double scaleMS_egLoss, InfoHelper& muonInfo ) const {
 
-    double scaleID = 0., enLossCorrMS = 0., scaleMS = 0., scaleCB = 0.;//initialize all to 0
+    double scaleID = 0., enLossCorrMS = 0., enLossCorrCB = 0., scaleMS = 0., scaleCB = 0.;//initialize all to 0
     // These are alternative scale corrections (KPKM,KC,K,C) they are != 0. if Tscale != SCALE_DEFAULT.
     //double s1_ID = 0., s2_ID = 0., s1_MS = 0., s2_MS = 0., s1_CB = 0., s2_CB = 0.;//Description of these is in ApplyScale
 
@@ -2176,9 +2187,25 @@ namespace CP {
     if( m_scale_CB[muonInfo.detRegion] != -1 ) {
       scaleCB = m_scale_CB[muonInfo.detRegion] + scaleVarID * m_scaleSyst_CB[muonInfo.detRegion];
     }
-    else {
+    else 
+    {
       if( muonInfo.ptms ) scaleCB = ( scaleID * muonInfo.weightID + ( scaleMS + enLossCorrMS / muonInfo.ptms ) * muonInfo.weightMS );
-      else         scaleCB = scaleID; // was scaleID * muonInfo.weightID
+      else                scaleCB = scaleID; // was scaleID * muonInfo.weightID
+
+      if(m_doDirectCBCalib)
+      {
+        // TODO:: Fix the NPs
+        scaleCB       = scaleMS_scale > 0.  ? m_SUp_1_CB[muonInfo.detRegion] : m_SDw_1_CB[muonInfo.detRegion];
+        enLossCorrCB  = scaleMS_egLoss > 0. ? m_SUp_0_CB[muonInfo.detRegion] : m_SDw_0_CB[muonInfo.detRegion];
+
+        // ATH_MSG_VERBOSE( "CB Nom Scale:  = " << m_S_1_CB[muonInfo.detRegion] << " enLossCorrCB: " << m_S_0_CB[muonInfo.detRegion] );
+        // ATH_MSG_VERBOSE( "CB sys Scale:  = " << scaleMS_scale << " "<< scaleCB << " enLossCorrCB: " << scaleMS_egLoss <<" "<< enLossCorrCB );
+
+        scaleCB      =  m_S_1_CB[muonInfo.detRegion] + scaleMS_scale * scaleCB;
+        enLossCorrCB =  m_S_0_CB[muonInfo.detRegion] + scaleMS_egLoss * enLossCorrCB;
+
+        scaleCB += 1;
+      }
     }
 
     double tmpDelta = 1.;
@@ -2232,14 +2259,23 @@ namespace CP {
       ATH_MSG_VERBOSE( "Double-Checking outPtMS = " << outPtMS << " at fourth..." );
       return outPtMS;
     }
-    if( DetType == MCAST::DetectorType::CB ) {
-      if( int( inSmearID ) != DEFAULT_INIT_VAL && int( inSmearMS ) != DEFAULT_INIT_VAL ) {
+    if( DetType == MCAST::DetectorType::CB ) 
+    {
+      if(m_doDirectCBCalib)
+      {
+        tmpDelta = 1 + muonInfo.smearDeltaCBDirect;
+      }
+      else if( int( inSmearID ) != DEFAULT_INIT_VAL && int( inSmearMS ) != DEFAULT_INIT_VAL ) 
+      {
         tmpDelta = 1. + inSmearID * muonInfo.weightID + inSmearMS * muonInfo.weightMS;
       }
-      else if( int( inSmearID ) == DEFAULT_INIT_VAL && int( inSmearMS ) == DEFAULT_INIT_VAL ) {
+      else if( int( inSmearID ) == DEFAULT_INIT_VAL && int( inSmearMS ) == DEFAULT_INIT_VAL ) 
+      {
         tmpDelta = 1. + muonInfo.smearDeltaCB;
       }
-      else {
+      else 
+      {
+
       }
       //In these releases the smearing was applied to the pt before the scale
       if( m_Trel < MCAST::Release::Rel17_2_Sum13 ) {
@@ -2247,10 +2283,11 @@ namespace CP {
         if( m_Tsmear == MCAST::SmearingType::QoverPt ) outPtCB = ( tmpDelta == 0 ) ? MCAST_MAX_PT : outPtCB / tmpDelta;
       }
       ATH_MSG_VERBOSE( "Org pT " << outPtCB );
+      ATH_MSG_VERBOSE( "Checking energy loss for cb = " << enLossCorrCB );
       ATH_MSG_VERBOSE( "Checking scale corr for CB = " << scaleCB );
       ATH_MSG_VERBOSE( "Checking smearing corr for CB = " << tmpDelta );
 
-      outPtCB = ScaleApply( std::abs(outPtCB), scaleCB, 0., muonInfo );
+      outPtCB = ScaleApply( std::abs(outPtCB), scaleCB, enLossCorrCB, muonInfo );
       if( m_Trel >= MCAST::Release::Rel17_2_Sum13 ) {
         if( m_Tsmear == MCAST::SmearingType::Pt )  outPtCB = outPtCB * tmpDelta;
         if( m_Tsmear == MCAST::SmearingType::QoverPt ) outPtCB = ( tmpDelta == 0 ) ? MCAST_MAX_PT : outPtCB / tmpDelta;
@@ -2274,7 +2311,10 @@ namespace CP {
     std::string scale_val;
     // Check if FilesPath defined: if so override other configurations (advanced user setting, for debugging within MCP)
     if ( m_FilesPath == "" ) {
-      if (m_Trel >= MCAST::Release::Recs2020_03_03) {
+      if (m_Trel >= MCAST::Release::Recs2021_07_01) {
+        scale_val = PathResolverFindCalibFile( "MuonMomentumCorrections/" + m_release + "/Scale_" + m_algo + "_" + m_year + ".dat" );
+      }
+      else if (m_Trel >= MCAST::Release::Recs2020_03_03) {
         scale_val = PathResolverFindCalibFile( "MuonMomentumCorrections/" + m_release + "/Scale_" + m_algo + "_" + m_year + "_" + m_release + ".dat" );
       }
       else if( m_Trel >= MCAST::Release::PreRec ) {
@@ -2288,7 +2328,10 @@ namespace CP {
       }
     }
     else {
-      if (m_Trel >= MCAST::Release::Recs2020_03_03) {
+      if (m_Trel >= MCAST::Release::Recs2021_07_01) {
+        scale_val = m_FilesPath + m_release + "/Scale_" + m_algo + "_" + m_year  + ".dat";
+      }
+      else if (m_Trel >= MCAST::Release::Recs2020_03_03) {
         scale_val = m_FilesPath + m_release + "/Scale_" + m_algo + "_" + m_year + "_" + m_release + ".dat";
       }
       else {
@@ -2358,7 +2401,10 @@ namespace CP {
     std::string data_val;
     // Check if FilesPath defined: if so override other configurations (advanced user setting, for debugging within MCP)
     if ( m_FilesPath == "" ) {
-      if (m_Trel >= MCAST::Release::Recs2020_03_03) {
+      if (m_Trel >= MCAST::Release::Recs2021_07_01) {
+        data_val = PathResolverFindCalibFile( "MuonMomentumCorrections/" + m_release + "/Smearing_" + m_algo + "_" + m_year + ".dat" );
+      }
+      else if (m_Trel >= MCAST::Release::Recs2020_03_03) {
         data_val = PathResolverFindCalibFile( "MuonMomentumCorrections/" + m_release + "/Smearing_" + m_algo + "_" + m_year + "_" + m_release + ".dat" );
       }
       else if( m_Trel >= MCAST::Release::PreRec ) {
@@ -2369,7 +2415,10 @@ namespace CP {
       }
     }
     else {
-      if (m_Trel >= MCAST::Release::Recs2020_03_03) {
+      if (m_Trel >= MCAST::Release::Recs2021_07_01) {
+        data_val = m_FilesPath + m_release + "/Smearing_" + m_algo + "_" + m_year + ".dat";
+      }
+      else if (m_Trel >= MCAST::Release::Recs2020_03_03) {
         data_val = m_FilesPath + m_release + "/Smearing_" + m_algo + "_" + m_year + "_" + m_release + ".dat";
       }
       else {
@@ -2471,7 +2520,10 @@ namespace CP {
     std::string mc_val;
     // Check if FilesPath defined: if so override other configurations (advanced user setting, for debugging within MCP)
     if ( m_FilesPath == "" ) {
-      if (m_Trel >= MCAST::Release::Recs2020_03_03) {
+      if (m_Trel >= MCAST::Release::Recs2021_07_01) {
+        mc_val = PathResolverFindCalibFile( "MuonMomentumCorrections/" + m_release + "/MC_values_" + m_algo + "_" + m_year + ".dat" );
+      }
+      else if (m_Trel >= MCAST::Release::Recs2020_03_03) {
         mc_val = PathResolverFindCalibFile( "MuonMomentumCorrections/" + m_release + "/MC_values_" + m_algo + "_" + m_year + "_" + m_release + ".dat" );
       }
       else if ( m_Trel >= MCAST::Release::PreRec ) {
@@ -2485,7 +2537,10 @@ namespace CP {
       }
     }
     else {
-      if (m_Trel >= MCAST::Release::Recs2020_03_03) {
+      if (m_Trel >= MCAST::Release::Recs2021_07_01) {
+        mc_val = m_FilesPath + m_release + "/MC_values_" + m_algo + "_" + m_year + ".dat"; 
+      }
+      else if (m_Trel >= MCAST::Release::Recs2020_03_03) {
         mc_val = m_FilesPath + m_release + "/MC_values_" + m_algo + "_" + m_year + "_" + m_release + ".dat"; 
       }
       else {
@@ -2536,9 +2591,15 @@ namespace CP {
       for(auto& kv: files_and_maps) {
         std::string extra_file;
         if(m_FilesPath == "")
-          extra_file = PathResolverFindCalibFile("MuonMomentumCorrections/" + m_release + "/" + kv.first + "_" + m_algo + "_" + m_year + "_" + m_release + ".dat");
+        { 
+          if (m_Trel >= MCAST::Release::Recs2021_07_01) extra_file = PathResolverFindCalibFile("MuonMomentumCorrections/" + m_release + "/" + kv.first + "_" + m_algo + "_" + m_year + ".dat");
+                                                   else extra_file = PathResolverFindCalibFile("MuonMomentumCorrections/" + m_release + "/" + kv.first + "_" + m_algo + "_" + m_year + "_" + m_release + ".dat");
+        }
         else
-          extra_file = m_FilesPath + m_release + "/" + kv.first + "_" + m_algo + "_" + m_year + "_" + m_release + ".dat";
+        {
+          if (m_Trel >= MCAST::Release::Recs2021_07_01) extra_file = m_FilesPath + m_release + "/" + kv.first + "_" + m_algo + "_" + m_year + ".dat";
+                                                   else extra_file = m_FilesPath + m_release + "/" + kv.first + "_" + m_algo + "_" + m_year + "_" + m_release + ".dat";
+        }
         ATH_MSG_DEBUG( "Checking Files - Extra smearing for high pt muons: " << extra_file );
         std::map<std::pair<int, int>, std::pair<double, double> >* this_map = kv.second;
 
@@ -2583,11 +2644,83 @@ namespace CP {
       }
     }
 
+
+    if(m_doDirectCBCalib)
+    {
+      std::string CBFile = "";
+      if ( m_FilesPath == "" ) 
+      {
+          if (m_Trel >= MCAST::Release::Recs2021_07_01) CBFile = PathResolverFindCalibFile( "MuonMomentumCorrections/" + m_release + "/CB_" + m_algo + "_" + m_year + ".dat" );
+                                                   else CBFile = PathResolverFindCalibFile( "MuonMomentumCorrections/" + m_release + "/CB_" + m_algo + "_" + m_year + "_" + m_release + ".dat" );
+      }
+      else 
+      {
+          if (m_Trel >= MCAST::Release::Recs2021_07_01) CBFile = m_FilesPath + m_release + "/CB_" + m_algo + "_" + m_year + ".dat";
+                                                   else CBFile = m_FilesPath + m_release + "/CB_" + m_algo + "_" + m_year + "_" + m_release + ".dat";
+      }
+      ATH_MSG_INFO( "Checking Files - Data: " << CBFile );
+
+      InValues.open( CBFile.c_str() );
+      i = 0;
+      if( !InValues.good() ) {
+        ATH_MSG_ERROR( "File " << CBFile << " not found!!" );
+        return StatusCode::FAILURE;
+      }
+      else 
+      {
+        while( InValues.good() && i < m_nb_regions ) 
+        {
+          tmpval = 0;
+          if( i == 0 ) {
+            getline( InValues,tmpname );
+          }
+
+          InValues>>tmpval;
+          m_S_0_CB.push_back( tmpval );
+          InValues>>tmpval;
+          m_SUp_0_CB.push_back( tmpval );
+          InValues>>tmpval;
+          m_SDw_0_CB.push_back( tmpval );
+
+          InValues>>tmpval;
+          m_S_1_CB.push_back( tmpval );
+          InValues>>tmpval;
+          m_SUp_1_CB.push_back( tmpval );
+          InValues>>tmpval;
+          m_SDw_1_CB.push_back( tmpval );
+
+
+          InValues>>tmpval;
+          m_R_0_CB.push_back( tmpval );
+          InValues>>tmpval;
+          m_RUp_0_CB.push_back( tmpval );
+          InValues>>tmpval;
+          m_RDw_0_CB.push_back( tmpval );
+
+          InValues>>tmpval;
+          m_R_1_CB.push_back( tmpval );
+          InValues>>tmpval;
+          m_RUp_1_CB.push_back( tmpval );
+          InValues>>tmpval;
+          m_RDw_1_CB.push_back( tmpval );
+
+          InValues>>tmpval;
+          m_R_2_CB.push_back( tmpval );
+          InValues>>tmpval;
+          m_RUp_2_CB.push_back( tmpval );
+          InValues>>tmpval;
+          m_RDw_2_CB.push_back( tmpval );
+        }
+      }
+
+
+    }
+
     return StatusCode::SUCCESS;
 
   }
 
-  double MuonCalibrationAndSmearingTool::GetSmearing( int DetType, InfoHelper& muonInfo ) const {
+  double MuonCalibrationAndSmearingTool::GetSmearing( int DetType, InfoHelper& muonInfo, bool doDirectCB ) const {
 
     bool useTan2 = true;
     bool useTan  = false;
@@ -2599,7 +2732,9 @@ namespace CP {
     }    
     if ( muonInfo.detRegion < 0 || muonInfo.detRegion >= m_nb_regions ) return 0; //++++++ HOW TO IMPROVE THIS CHECK?!
     double smear = 0.;
-    if ( DetType == MCAST::DetectorType::CB ) {
+    if ( DetType == MCAST::DetectorType::CB && !doDirectCB) 
+    {
+      // This is for High Pt
       ATH_MSG_VERBOSE("[Direct CB Smearing] Map Key: " << muonInfo.detRegion << ", " << ConvertToMacroCategory(muonInfo.sel_category));
       std::pair<int, int> key = std::make_pair(muonInfo.detRegion, ConvertToMacroCategory(muonInfo.sel_category)); 
       if((m_extra_p1_p2_MS_Misaligned.find(key) != m_extra_p1_p2_MS_Misaligned.end()) and (m_extra_p1_p2_MS_AlignedAndCorrected.find(key) != m_extra_p1_p2_MS_AlignedAndCorrected.end())) {
@@ -2613,7 +2748,20 @@ namespace CP {
       }
       return smear; 
     }
-    if ( DetType == MCAST::DetectorType::MS ) {
+    else if ( DetType == MCAST::DetectorType::CB && doDirectCB ) 
+    {
+      if( muonInfo.ptcb == 0 ) 
+      {
+        return 0.;
+      }
+      else {
+        ATH_MSG_VERBOSE( "CB: Using p0CB = " << m_R_0_CB[muonInfo.detRegion] << " and p1CB = " << m_R_1_CB[muonInfo.detRegion] << " and p2CB = " << m_R_2_CB[muonInfo.detRegion] );
+        smear = m_R_0_CB[muonInfo.detRegion]*muonInfo.g0/muonInfo.ptcb + m_R_1_CB[muonInfo.detRegion]*muonInfo.g1 + m_R_2_CB[muonInfo.detRegion]*muonInfo.g2*muonInfo.ptcb;
+        return smear;
+      }
+    }    
+    else if ( DetType == MCAST::DetectorType::MS ) 
+    {
       if( muonInfo.ptms == 0 ) {
         return 0.;
       }
@@ -3300,7 +3448,7 @@ namespace CP {
 
   }
 
-  double MuonCalibrationAndSmearingTool::GetSystVariation( int DetType, double var, InfoHelper& muonInfo ) const {
+  double MuonCalibrationAndSmearingTool::GetSystVariation( int DetType, double var, InfoHelper& muonInfo, bool doDirectCB ) const {
 
     if( var != 1. && var != -1 ) {
       ATH_MSG_WARNING( "Systematic variation is not +/- 1 sigma, returning 0." );
@@ -3310,8 +3458,8 @@ namespace CP {
 
     double p0_MS_var = 0., p1_MS_var = 0., p2_MS_var = 0., p1_ID_var = 0., p2_ID_var = 0., p2_ID_TAN_var = 0.;
     double newSmear = 0.;
-    if ( DetType == MCAST::DetectorType::CB ) {
-        double smear = GetSmearing(DetType, muonInfo);
+    if ( DetType == MCAST::DetectorType::CB && !doDirectCB ) {
+        double smear = GetSmearing(DetType, muonInfo, doDirectCB);
         std::pair<int, int> key = std::make_pair(muonInfo.detRegion, ConvertToMacroCategory(muonInfo.sel_category)); 
         ATH_MSG_VERBOSE("[Direct CB Smearing] Map Key: " << muonInfo.detRegion << ", " << muonInfo.sel_category);
         if((m_extra_p1_p2_MS_Misaligned.find(key) != m_extra_p1_p2_MS_Misaligned.end()) and (m_extra_p1_p2_MS_AlignedOnly.find(key) != m_extra_p1_p2_MS_AlignedOnly.end())) {
@@ -3358,7 +3506,31 @@ namespace CP {
         newSmear = ( p0_MS_var * muonInfo.g0 / muonInfo.ptms + p1_MS_var * muonInfo.g1 + p2_MS_var * muonInfo.g2 * muonInfo.ptms );
         return newSmear;
       }
-    } else if( DetType == MCAST::DetectorType::ID ) {
+    } 
+    if( DetType == MCAST::DetectorType::CB && doDirectCB ) 
+    {
+      if( muonInfo.ptcb == 0 ) 
+      {
+        return 0;
+      } 
+      else 
+      {
+        p0_MS_var = var > 0. ? m_RUp_0_CB[ muonInfo.detRegion ] : m_RDw_0_CB[ muonInfo.detRegion ];
+        p1_MS_var = var > 0. ? m_RUp_1_CB[ muonInfo.detRegion ] : m_RDw_1_CB[ muonInfo.detRegion ];
+        p2_MS_var = var > 0. ? m_RUp_2_CB[ muonInfo.detRegion ] : m_RDw_2_CB[ muonInfo.detRegion ];
+        
+        p0_MS_var = m_R_0_CB[ muonInfo.detRegion ] + var * p0_MS_var;
+        p1_MS_var = m_R_1_CB[ muonInfo.detRegion ] + var * p1_MS_var;
+        p2_MS_var = m_R_2_CB[ muonInfo.detRegion ] + var * p2_MS_var;
+
+        if( p0_MS_var < 0. ) p0_MS_var = 0.; //Truncation of unphysical variations
+        if( p1_MS_var < 0. ) p1_MS_var = 0.;
+        if( p2_MS_var < 0. ) p2_MS_var = 0.;
+        newSmear = ( p0_MS_var * muonInfo.g0 / muonInfo.ptcb + p1_MS_var * muonInfo.g1 + p2_MS_var * muonInfo.g2 * muonInfo.ptcb );
+        return newSmear;
+      }
+    } 
+    if( DetType == MCAST::DetectorType::ID ) {
       if( m_Trel < MCAST::Release::Rel17_2_Sum13 ) {
         p1_ID_var     = std::pow( m_E_p1_ID[ muonInfo.detRegion ] * m_E_p1_ID[ muonInfo.detRegion ] + m_S_p1_ID[ muonInfo.detRegion ] * m_S_p1_ID[ muonInfo.detRegion ], 0.5 );
         p2_ID_var     = std::pow( m_E_p2_ID[ muonInfo.detRegion ] * m_E_p2_ID[ muonInfo.detRegion ] + m_S_p2_ID[ muonInfo.detRegion ] * m_S_p2_ID[ muonInfo.detRegion ], 0.5 );
diff --git a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/util/MCAST_Tester.cxx b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/util/MCAST_Tester.cxx
index 20fa7d7384e3..36effe0d7874 100644
--- a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/util/MCAST_Tester.cxx
+++ b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/util/MCAST_Tester.cxx
@@ -130,7 +130,7 @@ int main( int argc, char* argv[] ) {
   corrTool.setProperty("Year",                  "Data17" );
   //   corrTool.setProperty("Algo",                  "muons" );
   //   corrTool.setProperty("SmearingType",          "q_pT" );
-  corrTool.setProperty("Release",               "Recs2021_04_18_Prelim" );
+  corrTool.setProperty("Release",               "Recs2021_11_04" );
 
   corrTool.setProperty("StatComb",              false);
   //   corrTool.setProperty("MinCombPt",             300.0);
@@ -170,6 +170,7 @@ int main( int argc, char* argv[] ) {
   if(nEvents >= 0 || Ievent >= 0) isDebug = true;
 
   if(isDebug) corrTool.setProperty( "OutputLevel", MSG::VERBOSE);
+  Info( APP_NAME, "is debug: ", isDebug);
 
   //::: retrieve the tool
   corrTool.retrieve();
-- 
GitLab


From bf13810d2c1273d3ec5ef49e4431fa65cbdca846 Mon Sep 17 00:00:00 2001
From: Syed Haider Abidi <syed.haider.abidi@cern.ch>
Date: Mon, 8 Nov 2021 14:11:57 +0100
Subject: [PATCH 11/14] default values to recover the current recommendation

---
 .../Root/MuonCalibrationAndSmearingTool.cxx   | 20 +++---
 .../util/MCAST_Tester.cxx                     | 67 ++++++++++---------
 2 files changed, 45 insertions(+), 42 deletions(-)

diff --git a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx
index cb4be124157f..780afb0ef55d 100644
--- a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx
+++ b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx
@@ -46,7 +46,7 @@ namespace CP {
     declareProperty("Year", m_year = "Data16" );
     declareProperty("Algo", m_algo = "muons" );
     declareProperty("SmearingType", m_type = "q_pT" );
-    declareProperty("Release", m_release = "Recs2021_11_04" );
+    declareProperty("Release", m_release = "Recs2020_03_03" ); // new Recs2021_11_04
     declareProperty("ToroidOff", m_toroidOff = false );
     declareProperty("AddExtraDecorations", m_extra_decorations = false );
     declareProperty("doExtraSmearing", m_extra_highpt_smearing = false );
@@ -55,26 +55,26 @@ namespace CP {
     declareProperty("StatComb", m_useStatComb = false);
     declareProperty("MinCombPt", m_StatCombPtThreshold=300.0);
     declareProperty("HighPtSystThr", m_HighPtSystThreshold=300.0);
-    declareProperty("SagittaCorr", m_doSagittaCorrection = true);
-    declareProperty("SagittaRelease", m_SagittaRelease = "sagittaBiasDataAll_15_09_2021");
-    declareProperty("doSagittaMCDistortion",m_doSagittaMCDistortion=false);
-    declareProperty("SagittaCorrPhaseSpace",m_SagittaCorrPhaseSpace=false);
-    declareProperty("SagittaIterWeight", m_IterWeight=1.0);
+    declareProperty("SagittaCorr", m_doSagittaCorrection = false); // new true
+    declareProperty("SagittaRelease", m_SagittaRelease = "sagittaBiasDataAll_03_02_19"); // new sagittaBiasDataAll_15_09_2021
+    declareProperty("doSagittaMCDistortion",m_doSagittaMCDistortion=true); // new false
+    declareProperty("SagittaCorrPhaseSpace",m_SagittaCorrPhaseSpace=false); 
+    declareProperty("SagittaIterWeight", m_IterWeight=0.5); // new 1.0
     declareProperty("sgItersCB",m_sgItersCB=4);
     declareProperty("sgItersID",m_sgItersID=4);
     declareProperty("sgItersME",m_sgItersME=4);
     declareProperty("sgIetrsManual",m_sgIetrsManual=false);
     declareProperty("fixedRho",m_fixedRho=1.0);
-    declareProperty("useFixedRho",m_useFixedRho=true);
+    declareProperty("useFixedRho",m_useFixedRho=false); // new true
     declareProperty("noEigenDecor" ,m_doNotUseAMGMATRIXDECOR=false);
     declareProperty("useExternalSeed" ,m_useExternalSeed=false);
     declareProperty("externalSeed" ,m_externalSeed=0);
-    declareProperty("sagittaMapsInputType", m_saggitaMapsInputType=MCAST::SagittaInputHistType::SINGLE);
-    declareProperty("sagittaMapUnitConversion",m_sagittaMapUnitConversion=1);
+    declareProperty("sagittaMapsInputType", m_saggitaMapsInputType=MCAST::SagittaInputHistType::NOMINAL); // new MCAST::SagittaInputHistType::SINGLE
+    declareProperty("sagittaMapUnitConversion",m_sagittaMapUnitConversion=1e-3); // new calib 1 
     declareProperty("systematicCorrelationScheme", m_sysScheme = "Corr_Scale");
     declareProperty("expertMode", m_expertMode = false);
     declareProperty("expertSetting_isData", m_expertMode_isData = false);
-    declareProperty("expert_extraRebiasSys", m_extraRebiasSys = 0.00003);
+    declareProperty("expert_extraRebiasSys", m_extraRebiasSys = 0.0); // new calib 0.00003
     declareProperty("expert_doDirectCBCalib", m_doDirectCBCalib = false);
 
     m_SagittaIterations.push_back(0);
diff --git a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/util/MCAST_Tester.cxx b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/util/MCAST_Tester.cxx
index 36effe0d7874..b82a8d90f7e0 100644
--- a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/util/MCAST_Tester.cxx
+++ b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/util/MCAST_Tester.cxx
@@ -126,44 +126,47 @@ int main( int argc, char* argv[] ) {
   //::: create the tool handle
   asg::AnaToolHandle<CP::IMuonCalibrationAndSmearingTool> corrTool; //!
   corrTool.setTypeAndName("CP::MuonCalibrationAndSmearingTool/MuonCorrectionTool");
-  //::: set the properties
-  corrTool.setProperty("Year",                  "Data17" );
-  //   corrTool.setProperty("Algo",                  "muons" );
-  //   corrTool.setProperty("SmearingType",          "q_pT" );
-  corrTool.setProperty("Release",               "Recs2021_11_04" );
-
-  corrTool.setProperty("StatComb",              false);
-  //   corrTool.setProperty("MinCombPt",             300.0);
-  corrTool.setProperty("SagittaCorr",           true);
-  corrTool.setProperty("SagittaRelease",        "sagittaBiasDataAll_15_09_2021");
-  corrTool.setProperty("doSagittaMCDistortion", false);
-  corrTool.setProperty("SagittaCorrPhaseSpace", false);
-  corrTool.setProperty("SagittaIterWeight", 1);
-  //   corrTool.setProperty("sgItersCB",             11);
-  //   corrTool.setProperty("sgItersID",             11);
-  //   corrTool.setProperty("sgItersME",             11);
-  //   corrTool.setProperty("sgIetrsMamual",         false);
-  corrTool.setProperty("fixedRho",              1.0);
-  corrTool.setProperty("useFixedRho",           true);
-  corrTool.setProperty("noEigenDecor" ,         false);
-  corrTool.setProperty("sagittaMapsInputType" ,         1);
-  corrTool.setProperty("sagittaMapUnitConversion",1); 
-  // corrTool.setProperty("systematicCorrelationScheme", "Decorr_Scale"); 
-  //   corrTool.setProperty("useExternalSeed" ,      false);
-  //   corrTool.setProperty("externalSeed" ,         0);
-
-
+  // New tool settings
+  // //::: set the properties
   // corrTool.setProperty("Year",                  "Data17" );
-  // corrTool.setProperty("Release",               "Recs2021_04_18_Prelim" );
+  // //   corrTool.setProperty("Algo",                  "muons" );
+  // //   corrTool.setProperty("SmearingType",          "q_pT" );
+  // corrTool.setProperty("Release",               "Recs2021_11_04" );
+
   // corrTool.setProperty("StatComb",              false);
+  // //   corrTool.setProperty("MinCombPt",             300.0);
   // corrTool.setProperty("SagittaCorr",           true);
-  // corrTool.setProperty("SagittaRelease",        "sagittaBiasDataAll_03_02_19_Data17");
+  // corrTool.setProperty("SagittaRelease",        "sagittaBiasDataAll_15_09_2021");
   // corrTool.setProperty("doSagittaMCDistortion", false);
-  // corrTool.setProperty("SagittaCorrPhaseSpace", true);
+  // corrTool.setProperty("SagittaCorrPhaseSpace", false);
   // corrTool.setProperty("SagittaIterWeight", 1);
-  // corrTool.setProperty("fixedRho",              0.0);
-  // corrTool.setProperty("useFixedRho",           false);
+  // //   corrTool.setProperty("sgItersCB",             11);
+  // //   corrTool.setProperty("sgItersID",             11);
+  // //   corrTool.setProperty("sgItersME",             11);
+  // //   corrTool.setProperty("sgIetrsMamual",         false);
+  // corrTool.setProperty("fixedRho",              1.0);
+  // corrTool.setProperty("useFixedRho",           true);
   // corrTool.setProperty("noEigenDecor" ,         false);
+  // corrTool.setProperty("sagittaMapsInputType" ,         1);
+  // corrTool.setProperty("sagittaMapUnitConversion",1); 
+  // // corrTool.setProperty("systematicCorrelationScheme", "Decorr_Scale"); 
+  // //   corrTool.setProperty("useExternalSeed" ,      false);
+  // //   corrTool.setProperty("externalSeed" ,         0);
+
+
+  corrTool.setProperty("Year",                  "Data17" );
+  corrTool.setProperty("Release", "Recs2020_03_03");
+  corrTool.setProperty("SagittaRelease", "sagittaBiasDataAll_03_02_19_Data17");
+  corrTool.setProperty("StatComb", false);
+  corrTool.setProperty("SagittaCorr", true);
+  corrTool.setProperty("doSagittaMCDistortion", false);
+  corrTool.setProperty("SagittaCorrPhaseSpace", true);
+  // corrTool.setProperty("SagittaIterWeight", 0.5);
+  // corrTool.setProperty("fixedRho", 1);
+  // corrTool.setProperty("useFixedRho", false);
+  // corrTool.setProperty("sagittaMapsInputType", 0);
+  // corrTool.setProperty("sagittaMapUnitConversion", 1e-3);
+  // corrTool.setProperty("systematicCorrelationScheme", "Corr_Scale");
 
 
   bool isDebug = false;
-- 
GitLab


From b8f8df60c62b8d34f3318c8af498ec4106c5182a Mon Sep 17 00:00:00 2001
From: Haider Abidi <syed.haider.abidi@cern.ch>
Date: Tue, 9 Nov 2021 11:34:46 +0100
Subject: [PATCH 12/14] Update CMakeLists.txt

---
 .../MuonID/MuonIDAnalysis/MuonMomentumCorrections/CMakeLists.txt | 1 -
 1 file changed, 1 deletion(-)

diff --git a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/CMakeLists.txt b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/CMakeLists.txt
index 9d55fcedfcf0..993bc7d99f23 100644
--- a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/CMakeLists.txt
+++ b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/CMakeLists.txt
@@ -66,4 +66,3 @@ endif()
 
 # Install files from the package:
 atlas_install_joboptions( share/*.py )
-atlas_install_data( data/* )
-- 
GitLab


From a174cc6aa961073c86febf9e0b97a486bd36313b Mon Sep 17 00:00:00 2001
From: Haider Abidi <syed.haider.abidi@cern.ch>
Date: Tue, 9 Nov 2021 11:38:14 +0100
Subject: [PATCH 13/14] L2 comments

---
 .../Root/MuonCalibrationAndSmearingTool.cxx                | 7 ++-----
 1 file changed, 2 insertions(+), 5 deletions(-)

diff --git a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx
index 3c04db082fd6..d3c68c519a54 100644
--- a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx
+++ b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/Root/MuonCalibrationAndSmearingTool.cxx
@@ -2198,9 +2198,6 @@ namespace CP {
         scaleCB       = scaleMS_scale > 0.  ? m_SUp_1_CB[muonInfo.detRegion] : m_SDw_1_CB[muonInfo.detRegion];
         enLossCorrCB  = scaleMS_egLoss > 0. ? m_SUp_0_CB[muonInfo.detRegion] : m_SDw_0_CB[muonInfo.detRegion];
 
-        // ATH_MSG_VERBOSE( "CB Nom Scale:  = " << m_S_1_CB[muonInfo.detRegion] << " enLossCorrCB: " << m_S_0_CB[muonInfo.detRegion] );
-        // ATH_MSG_VERBOSE( "CB sys Scale:  = " << scaleMS_scale << " "<< scaleCB << " enLossCorrCB: " << scaleMS_egLoss <<" "<< enLossCorrCB );
-
         scaleCB      =  m_S_1_CB[muonInfo.detRegion] + scaleMS_scale * scaleCB;
         enLossCorrCB =  m_S_0_CB[muonInfo.detRegion] + scaleMS_egLoss * enLossCorrCB;
 
@@ -2220,7 +2217,7 @@ namespace CP {
       tmpDelta = ( int( inSmearID ) == DEFAULT_INIT_VAL ) ? 1. + muonInfo.smearDeltaID : 1. + inSmearID;
       ATH_MSG_VERBOSE( "Double-Checking inSmearID = " << inSmearID );
       ATH_MSG_VERBOSE( "Double-Checking muonInfo.smearDeltaID = " << muonInfo.smearDeltaID );
-      ATH_MSG_VERBOSE( "Double-Checking denominator of smeaering = " << tmpDelta );
+      ATH_MSG_VERBOSE( "Double-Checking denominator of smearing = " << tmpDelta );
       if( m_Trel < MCAST::Release::Rel17_2_Sum13 ) {
         if( m_Tsmear == MCAST::SmearingType::Pt )  outPtID = outPtID * tmpDelta;
         if( m_Tsmear == MCAST::SmearingType::QoverPt ) outPtID = ( tmpDelta == 0 ) ? MCAST_MAX_PT : outPtID / tmpDelta;
@@ -2243,7 +2240,7 @@ namespace CP {
       tmpDelta = ( int( inSmearMS ) == DEFAULT_INIT_VAL ) ? 1. + muonInfo.smearDeltaMS : 1. + inSmearMS;
       ATH_MSG_VERBOSE( "Double-Checking inSmearMS = " << inSmearMS );
       ATH_MSG_VERBOSE( "Double-Checking muonInfo.smearDeltaMS = " << muonInfo.smearDeltaMS );
-      ATH_MSG_VERBOSE( "Double-Checking denominator of smeaering = " << tmpDelta );
+      ATH_MSG_VERBOSE( "Double-Checking denominator of smearing = " << tmpDelta );
       //In these releases the smearing was applied to the pt before the scale
       if( m_Trel < MCAST::Release::Rel17_2_Sum13 ) {
         if( m_Tsmear == MCAST::SmearingType::Pt )  outPtMS = outPtMS * tmpDelta;
-- 
GitLab


From 2006c71bf5e0b59fe88b1b9db95cd98320a836b1 Mon Sep 17 00:00:00 2001
From: Haider Abidi <syed.haider.abidi@cern.ch>
Date: Tue, 9 Nov 2021 13:48:35 +0100
Subject: [PATCH 14/14] Update MCAST_Tester.cxx

---
 .../MuonIDAnalysis/MuonMomentumCorrections/util/MCAST_Tester.cxx | 1 -
 1 file changed, 1 deletion(-)

diff --git a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/util/MCAST_Tester.cxx b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/util/MCAST_Tester.cxx
index b82a8d90f7e0..0c236dc56c1b 100644
--- a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/util/MCAST_Tester.cxx
+++ b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/util/MCAST_Tester.cxx
@@ -173,7 +173,6 @@ int main( int argc, char* argv[] ) {
   if(nEvents >= 0 || Ievent >= 0) isDebug = true;
 
   if(isDebug) corrTool.setProperty( "OutputLevel", MSG::VERBOSE);
-  Info( APP_NAME, "is debug: ", isDebug);
 
   //::: retrieve the tool
   corrTool.retrieve();
-- 
GitLab