Skip to content

GitLab

Commit c7243abd authored by Nils Erik Krumnack's avatar Nils Erik Krumnack
Browse files

Merge branch 'MMC_sagittaBiasSystematics' into '21.2' 

Technical developments for new MCP calibration recommendations

See merge request !45435
parents 0d496eff 3e4c57b5
Expand all Hide whitespace changes
Inline Side-by-side
PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/MuonMomentumCorrections/MuonCalibrationAndSmearingTool.h
View file @ c7243abd
......@@ -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;
......
PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/util/MCAST_Tester.cxx
View file @ c7243abd
......@@ -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 ) {
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment