diff --git a/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/MuonMomentumCorrections/MuonCalibrationAndSmearingTool.h b/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/MuonMomentumCorrections/MuonCalibrationAndSmearingTool.h index dae845658046fda333fbafb779c63727d2e42ee9..1afc9c72864f4562a9d7702a0d8607e29e49f4cb 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 }; } } @@ -99,7 +99,11 @@ 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; + double uncorrected_ptms = 0; }; public: @@ -121,8 +125,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, @@ -130,7 +134,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: @@ -142,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; @@ -181,8 +185,12 @@ class MuonCalibrationAndSmearingTool : public virtual IMuonCalibrationAndSmearin double ScaleMS_egLoss; double SagittaRho; double SagittaBias; + double SagittaDataStat; }; + bool m_expertMode; + bool m_expertMode_isData; + bool m_useExternalSeed; int m_externalSeed; @@ -196,7 +204,9 @@ class MuonCalibrationAndSmearingTool : public virtual IMuonCalibrationAndSmearin 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. @@ -214,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 8ec9206cb4cc7f8183bb084c3e6d153448218a1c..d3c68c519a54213a325ca3391f2b729fdc8b44f2 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 = "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,29 +55,34 @@ 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("SagittaRelease", m_SagittaRelease = "sagittaBiasDataAll_03_02_19"); - declareProperty("doSagittaMCDistortion",m_doSagittaMCDistortion=true); - declareProperty("SagittaCorrPhaseSpace",m_SagittaCorrPhaseSpace=false); - declareProperty("SagittaIterWeight", m_IterWeight=0.5); + 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=false); + 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::NOMINAL); - declareProperty("sagittaMapUnitConversion",m_sagittaMapUnitConversion=1e-3); + 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.0); // new calib 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 ) : @@ -173,15 +178,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); } @@ -226,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; @@ -359,41 +365,59 @@ 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"); - 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); - } - 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){ - // 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))); + 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. @@ -502,73 +526,150 @@ namespace CP { } double p2=corrM->GetBinContent(binEta,binPhi); - if(m_currentParameters->SagittaBias != MCAST::SystVariation::Default){ - if(m_currentParameters->SagittaBias == MCAST::SystVariation::Up){ - p2=0.5*p2; - } - else if (m_currentParameters->SagittaBias == 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; + } + + if(m_currentParameters->SagittaBias == MCAST::SystVariation::Up) + { + p2 += m_extraRebiasSys; + } + else if (m_currentParameters->SagittaBias == MCAST::SystVariation::Down) + { + p2 -= m_extraRebiasSys; } + + 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 { + 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; - ATH_MSG_DEBUG("Sagitta correction method "<< SgCorrType <<" iter "<<iter); + 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; + } + // 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 ( 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(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; - if(iter == 0){ - 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.primaryTrackParticleLink() ).isValid() ) { + const ElementLink< xAOD::TrackParticleContainer >& cb_track = mu.primaryTrackParticleLink(); + muonInfo.ptcb = ( !cb_track ) ? 0. : ( *cb_track )->pt() * 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 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.ptid==0 && muonInfo.ptms ==0){ - ATH_MSG_DEBUG("ME and ID momenta == 0"); + if(muonInfo.ptcb == 0 && SgCorrType==MCAST::SagittaCorType::CB) + { + 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; } @@ -589,12 +690,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); @@ -603,188 +704,150 @@ 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); + 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); + 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); + if(!stop) return applySagittaBiasCorrection(MCAST::SagittaCorType::ME, mu, iter, stop, isMC, muonInfo,SystCase); } - 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); + else if( SgCorrType == MCAST::SagittaCorType::WEIGHTS) + { + CalcCBWeights( mu, muonInfo ); - double MEqOverPE = 1e10; - if(ME_track!=nullptr) - MEqOverPE = std::pow( ME_track->definingParametersCovMatrix()( 4, 4 )/1e3,2); + if( iter == m_sagittasCB.size() ) + { + return CorrectionCode::Ok; + } + double tmpPtID = lvID.Pt(); //muonInfo.ptid; + double tmpPtMS = lvME.Pt(); //muonInfo.ptms; - CalcCBWeights( mu, muonInfo ); + double tmpDeltaID=0; + double tmpDeltaMS=0; + int itersID = 0; + int itersME = 0; - 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); + if ((SystCase != MCAST::SagittaSysType::NOMINAL && SystCase != MCAST::SagittaSysType::RHO)) + { + itersID = iter; + itersME = iter; + } - 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; + // Calculate the stat combination before sagitta bias: + + double chi2Nom=-999; + AmgVector(5) parsCBNom, parsCB; + AmgSymMatrix(5) covCBNom, covCB; + int charge = mu.charge(); - if(iter==0){ + AmgVector(5) parsID; + AmgVector(5) parsMS; + AmgSymMatrix(5) covID; + AmgSymMatrix(5) covMS; - 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(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"); } - - - if( iter == m_sagittasCB.size() ) { - if(dump) - ATH_MSG_DEBUG("--> " << muonInfo.ptcb); - return CorrectionCode::Ok; + 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(); + covCBNom = CB_track->definingParametersCovMatrix(); + covID = ID_track->definingParametersCovMatrix(); + covMS = ME_track->definingParametersCovMatrix(); } - //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; + CorrectionCode NominalCorrCode=applyStatCombination(parsID, covID, + parsMS, covMS, + charge, + parsCBNom, + covCBNom, + chi2Nom); + if(NominalCorrCode!=CorrectionCode::Ok) return NominalCorrCode; - double tmpPtID = lvID.Pt(); //muonInfo.ptid; - double tmpPtMS = lvME.Pt(); //muonInfo.ptms; - double tmpDeltaID=0; - double tmpDeltaMS=0; - CorrectionCode idOK=applySagittaBiasCorrection(MCAST::SagittaCorType::ID, mu,0,false, isMC, muonInfo); + 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; - CorrectionCode meOK=applySagittaBiasCorrection(MCAST::SagittaCorType::ME, mu,0, false, isMC, muonInfo); + 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 - AmgVector(5) parsMS = ME_track->definingParameters(); parsMS[4]=1.0 / (lvMECorr.P()*1e3); 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; - 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(SysCorrCode!=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 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; muonInfo.ptms = stored_ptms; @@ -805,45 +868,80 @@ 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; + } + } + 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){ + 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 - return applySagittaBiasCorrection(MCAST::SagittaCorType::ID, mu, itersID, false, isMC, muonInfo); + + 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 - return applySagittaBiasCorrection(MCAST::SagittaCorType::ME, mu, itersME, false, isMC, muonInfo); + 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; } } @@ -852,80 +950,149 @@ 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); + + 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; - if(muonInfo.ptid == 0 && muonInfo.ptms != 0 ) { - if(applySagittaBiasCorrection(MCAST::SagittaCorType::ME, mu, itersME, false, isMC, muonInfo)!=CorrectionCode::Ok){ + if(muonInfo.ptid == 0 && muonInfo.ptms != 0 ) + { + 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)!=CorrectionCode::Ok){ + 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"); 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; + + 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"); - if(applySagittaBiasCorrection(MCAST::SagittaCorType::CB, mu, itersCB, false, isMC, muonInfo)!=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)!=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; } @@ -936,16 +1103,32 @@ 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 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; + + 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 + { muonInfo.ptcb = rho*ptCB + (1-rho)*ptWeight; } @@ -963,76 +1146,78 @@ 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: - 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() ); + muonInfo.charge = mu.charge(); + muonInfo.ptid = id_pt(mu) * MeVtoGeV; + muonInfo.uncorrected_ptid = muonInfo.ptid; - if( ( mu.inDetTrackParticleLink() ).isValid() ) { - const ElementLink< xAOD::TrackParticleContainer >& id_track = mu.inDetTrackParticleLink(); - muonInfo.ptid = ( !id_track ) ? 0. : ( *id_track )->pt() * MeVtoGeV; - } else { - ATH_MSG_VERBOSE("The ID track link is not valid - setting pT(ID)=0"); - muonInfo.ptid = 0.; - } + muonInfo.ptms = ms_pt(mu) * MeVtoGeV; + muonInfo.uncorrected_ptms = muonInfo.ptms; + + muonInfo.ptcb = cb_pt(mu) * MeVtoGeV; + muonInfo.uncorrected_ptcb = muonInfo.ptcb; - // 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; - } - else{ - ATH_MSG_VERBOSE("No link to extrapolatedMuonSpectrometerTrackParticleLink setting pT(MS)=0"); - muonInfo.ptms = 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 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 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.; + } - if( ( mu.primaryTrackParticleLink() ).isValid() ) { - const ElementLink< xAOD::TrackParticleContainer >& cb_track = mu.primaryTrackParticleLink(); - muonInfo.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"); - muonInfo.ptcb = 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(); + // 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(); + } } - ATH_MSG_VERBOSE( "Setting Eta [CB]..." ); muonInfo.eta = mu.eta(); - ATH_MSG_VERBOSE( "Setting Phi [CB]..." ); muonInfo.phi = mu.phi(); // Getting detector region @@ -1058,11 +1243,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; } } } @@ -1075,17 +1262,22 @@ 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" ); - 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){ @@ -1093,10 +1285,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() ); } @@ -1109,7 +1304,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 @@ -1140,8 +1335,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: @@ -1149,37 +1344,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 ); + + CorrectionCode sgCode = CorrectionCode::Ok; + 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) + // 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; - 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 @@ -1199,15 +1419,10 @@ 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() ) { - 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 ); } @@ -1236,6 +1451,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; } @@ -1249,12 +1466,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 { @@ -1318,11 +1541,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; @@ -1375,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" ); @@ -1472,7 +1703,14 @@ 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 +1744,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 +1756,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 +1766,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 +1781,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 +1791,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 +1806,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 +1816,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 +1831,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 +1841,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 +1856,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 +1866,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 +1881,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 +1891,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 +1907,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 +1917,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 +1933,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 +1943,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 ); @@ -1834,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; @@ -1858,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 @@ -1899,15 +2187,28 @@ 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]; + + 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.; 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 @@ -1916,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 tmpDelta = " << 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; @@ -1939,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 tmpDelta = " << 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; @@ -1955,26 +2256,40 @@ 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 ) { 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 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; } + ATH_MSG_VERBOSE( "Corrected pT " << outPtCB ); return outPtCB; } @@ -1993,7 +2308,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 ) { @@ -2007,7 +2325,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 { @@ -2077,7 +2398,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 ) { @@ -2088,7 +2412,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 { @@ -2190,7 +2517,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 ) { @@ -2204,7 +2534,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 { @@ -2255,9 +2588,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; @@ -2302,11 +2641,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; @@ -2318,7 +2729,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())) { @@ -2332,11 +2745,25 @@ 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.; } 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; } @@ -2468,62 +2895,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: - - 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; - } - 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: + loc_ptid = id_pt(mu) * MeVtoGeV; + loc_ptms = ms_pt(mu) * MeVtoGeV; + loc_ptcb = cb_pt(mu) * MeVtoGeV; - 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; - } - else{ - ATH_MSG_VERBOSE("No link to extrapolatedMuonSpectrometerTrackParticleLink setting pT(MS)=0"); - loc_ptms = 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 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 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.; + } - 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.; + // 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) @@ -2543,14 +2960,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; @@ -2771,7 +3180,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] ) { @@ -2793,8 +3201,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; @@ -3037,7 +3445,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." ); @@ -3047,8 +3455,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())) { @@ -3095,7 +3503,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 ); @@ -3179,12 +3611,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(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(id_track,ms_track,charge,parsCB,covCB,chi2)!=CorrectionCode::Ok) + + if(applyStatCombination(parsID, covID, parsMS, covMS,charge,parsCB,covCB,chi2)!=CorrectionCode::Ok) return CorrectionCode::Error; muonInfo.cbParsA.clear(); @@ -3210,21 +3667,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 ){ @@ -3232,23 +3685,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 908114aed6c7b2c62f1851fc12d2c282eb488ea9..0c236dc56c1b1dcd0f2f815c836ff674aec2d030 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): @@ -125,33 +126,54 @@ int main( int argc, char* argv[] ) { //::: create the tool handle asg::AnaToolHandle<CP::IMuonCalibrationAndSmearingTool> corrTool; //! corrTool.setTypeAndName("CP::MuonCalibrationAndSmearingTool/MuonCorrectionTool"); - //::: set the properties + // New tool settings + // //::: 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); + + corrTool.setProperty("Year", "Data17" ); - // corrTool.setProperty("Algo", "muons" ); - // corrTool.setProperty("SmearingType", "q_pT" ); - corrTool.setProperty("Release", "Recs2020_03_03" ); - //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_07_04_2021"); + 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", false); - corrTool.setProperty("SagittaIterWeight", 1); - // corrTool.setProperty("sgItersCB", 11); - // corrTool.setProperty("sgItersID", 11); - // corrTool.setProperty("sgItersME", 11); - // corrTool.setProperty("sgIetrsMamual", false); - corrTool.setProperty("fixedRho", 0.0); - corrTool.setProperty("useFixedRho", false); - corrTool.setProperty("noEigenDecor" , false); - corrTool.setProperty("sagittaMapsInputType" , 1); - corrTool.setProperty("sagittaMapUnitConversion",1); - // corrTool.setProperty("useExternalSeed" , false); - // corrTool.setProperty("externalSeed" , 0); + 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; + if(nEvents >= 0 || Ievent >= 0) isDebug = true; + + if(isDebug) corrTool.setProperty( "OutputLevel", MSG::VERBOSE); + //::: retrieve the tool corrTool.retrieve(); @@ -166,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 ); @@ -203,9 +225,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" ); @@ -225,12 +249,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; } @@ -253,37 +282,45 @@ 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() ) ); - - // 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; - + std::map<int, float> cbpt, idpt, mspt; //::: 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 ); - Info( APP_NAME, "Looking at %s systematic", ( sysListItr->name() ).c_str() ); + xAOD::MuonContainer* muonsCorr = muons_shallowCopy.first; + + 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 ) { - + 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(); @@ -303,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()){ @@ -311,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()){ @@ -325,7 +362,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()); + 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 ) { @@ -338,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(), mu->auxdata< float >( "InnerDetectorPt" ), mu->auxdata< float >( "MuonSpectrometerPt" ) ); + + 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; @@ -354,17 +402,32 @@ 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, " 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++; + // 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 ) {