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