diff --git a/Reconstruction/Jet/JetCalibTools/JetCalibTools/CalibrationMethods/InsituDataCorrection.h b/Reconstruction/Jet/JetCalibTools/JetCalibTools/CalibrationMethods/InsituDataCorrection.h
index b002886676842e114bf5ed066471c28b529e4a10..d1c9957741a4e463df11955d2adcb3401ca5af99 100644
--- a/Reconstruction/Jet/JetCalibTools/JetCalibTools/CalibrationMethods/InsituDataCorrection.h
+++ b/Reconstruction/Jet/JetCalibTools/JetCalibTools/CalibrationMethods/InsituDataCorrection.h
@@ -1,5 +1,5 @@
/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+ Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
*/
#ifndef JETCALIBTOOLS_INSITUDATACORRECTION_H
@@ -42,18 +42,21 @@ class InsituDataCorrection
private:
double getInsituCorr(double pt, double eta, std::string calibstep) const;
+ double getInsituCorr_JMS(double pt, double mass, double eta, std::string calibstep, bool isTAmass) const;
TH2D * combineCalibration(TH2D *h2d, TH1D *h);
+ TH2D * invertHistogram(TH2D *h2d);
private:
TEnv * m_config;
TString m_jetAlgo, m_calibAreaTag;
bool m_dev;
- TH2D * m_insituCorr;
- double m_insituEtaMax, m_insituPtMin, m_insituPtMax;
+ std::unique_ptr<TH2D> m_insituCorr;
+ std::unique_ptr<TH2D> m_insituCorr_JMS;
+ std::unique_ptr<TH2D> m_insituCorr_JMS_TA;
+ double m_insituEtaMax, m_insituPtMin, m_insituPtMax, m_insituEtaMax_JMS, m_insituPtMin_JMS, m_insituPtMax_JMS, m_insituMassMin_JMS, m_insituMassMax_JMS;
double m_relhistoPtMax, m_abshistoPtMax;
-
- TH2D * m_insituCorr_ResidualMCbased;
+ std::unique_ptr<TH2D> m_insituCorr_ResidualMCbased;
double m_insituEtaMax_ResidualMCbased, m_insituPtMin_ResidualMCbased, m_insituPtMax_ResidualMCbased;
bool m_applyRelativeandAbsoluteInsitu;
@@ -62,6 +65,9 @@ class InsituDataCorrection
bool m_applyResidualMCbasedInsitu;
bool m_applyEtaRestrictionResidualMCbased;
+ bool m_applyInsituCaloTAjets;
+ bool m_applyInsituJMS;
+
};
#endif
diff --git a/Reconstruction/Jet/JetCalibTools/JetCalibTools/CalibrationMethods/JMSCorrection.h b/Reconstruction/Jet/JetCalibTools/JetCalibTools/CalibrationMethods/JMSCorrection.h
index f9c7ffe5c9defc1d755f819a1b0570f8983d7231..dab61091198551af89aa43016d0f41b68f4c00d3 100644
--- a/Reconstruction/Jet/JetCalibTools/JetCalibTools/CalibrationMethods/JMSCorrection.h
+++ b/Reconstruction/Jet/JetCalibTools/JetCalibTools/CalibrationMethods/JMSCorrection.h
@@ -1,5 +1,5 @@
/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+ Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
*/
#ifndef JETCALIBTOOLS_JMSCORRECTION_H
@@ -81,6 +81,7 @@ class JMSCorrection
bool m_combination; // Mass Combination of calo mass with track-assisted mass
bool m_useCorrelatedWeights;
+ bool m_onlyCombination; //mass combination using insitu calibrated inputs
// Control the binning using a private class enum
enum class BinningParam { pt_mass_eta, e_LOGmOe_eta, e_LOGmOet_eta, e_LOGmOpt_eta, et_LOGmOet_eta };
diff --git a/Reconstruction/Jet/JetCalibTools/JetCalibTools/JetCalibrationTool.h b/Reconstruction/Jet/JetCalibTools/JetCalibTools/JetCalibrationTool.h
index 711c8f1972dac7fa127d018a679fa31fb90e2a76..eed0f4ecb474c5bc60a21a17cc5a9c158745b725 100644
--- a/Reconstruction/Jet/JetCalibTools/JetCalibTools/JetCalibrationTool.h
+++ b/Reconstruction/Jet/JetCalibTools/JetCalibTools/JetCalibrationTool.h
@@ -119,10 +119,13 @@ private:
std::vector<double> m_runBins;
bool m_doSetDetectorEta;
std::string m_vertexContainerName;
+ bool m_insituCombMassCalib;
+ std::vector<TString> m_insituCombMassConfig;
//TEnv to hold the global text config
TEnv * m_globalConfig;
std::vector<TEnv*> m_globalTimeDependentConfigs;
+ std::vector<TEnv*> m_globalInsituCombMassConfig;
//Bools/enums to avoid string comparisons at run time
jetScale m_jetScale;
@@ -145,6 +148,8 @@ private:
std::vector<JetCalibrationToolBase*> m_insituTimeDependentCorr;
JMSCorrection * m_jetMassCorr;
JetSmearingCorrection* m_jetSmearCorr;
+ JMSCorrection *m_insituCombMassCorr_tmp;
+ std::vector<JetCalibrationToolBase*> m_insituCombMassCorr;
};
diff --git a/Reconstruction/Jet/JetCalibTools/Root/InsituDataCorrection.cxx b/Reconstruction/Jet/JetCalibTools/Root/InsituDataCorrection.cxx
index 67cb11ddf3f7b7fb4aca666752fb8ef1935c4ee3..1f191a99527d4e0cce7c5875d6d91a211d9a8551 100644
--- a/Reconstruction/Jet/JetCalibTools/Root/InsituDataCorrection.cxx
+++ b/Reconstruction/Jet/JetCalibTools/Root/InsituDataCorrection.cxx
@@ -1,5 +1,5 @@
/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+ Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
*/
#include "JetCalibTools/CalibrationMethods/InsituDataCorrection.h"
@@ -7,24 +7,39 @@
InsituDataCorrection::InsituDataCorrection()
: JetCalibrationToolBase::JetCalibrationToolBase("InsituDataCorrection::InsituDataCorrection"),
- m_config(NULL), m_jetAlgo(""), m_calibAreaTag(""), m_dev(false), m_insituCorr(NULL), m_insituCorr_ResidualMCbased(NULL)
+ m_config(nullptr),
+ m_jetAlgo(""),
+ m_calibAreaTag(""),
+ m_dev(false),
+ m_insituCorr(nullptr),
+ m_insituCorr_JMS(nullptr),
+ m_insituCorr_ResidualMCbased(nullptr)
{ }
InsituDataCorrection::InsituDataCorrection(const std::string& name)
: JetCalibrationToolBase::JetCalibrationToolBase( name ),
- m_config(NULL), m_jetAlgo(""), m_calibAreaTag(""), m_dev(false), m_insituCorr(NULL), m_insituCorr_ResidualMCbased(NULL)
+ m_config(nullptr),
+ m_jetAlgo(""),
+ m_calibAreaTag(""),
+ m_dev(false),
+ m_insituCorr(nullptr),
+ m_insituCorr_JMS(nullptr),
+ m_insituCorr_ResidualMCbased(nullptr)
{ }
InsituDataCorrection::InsituDataCorrection(const std::string& name, TEnv * config, TString jetAlgo, TString calibAreaTag, bool dev)
: JetCalibrationToolBase::JetCalibrationToolBase( name ),
- m_config(config), m_jetAlgo(jetAlgo), m_calibAreaTag(calibAreaTag), m_dev(dev), m_insituCorr(NULL), m_insituCorr_ResidualMCbased(NULL)
+ m_config(config),
+ m_jetAlgo(jetAlgo),
+ m_calibAreaTag(calibAreaTag),
+ m_dev(dev),
+ m_insituCorr(nullptr),
+ m_insituCorr_JMS(nullptr),
+ m_insituCorr_ResidualMCbased(nullptr)
{ }
InsituDataCorrection::~InsituDataCorrection() {
- if(m_insituCorr) delete m_insituCorr;
- if(m_insituCorr_ResidualMCbased) delete m_insituCorr_ResidualMCbased;
-
}
StatusCode InsituDataCorrection::initializeTool(const std::string&) {
@@ -53,6 +68,10 @@ StatusCode InsituDataCorrection::initializeTool(const std::string&) {
double insitu_etarestriction_residualmcbased = m_config->GetValue("InsituEtaRestrictionResidualMCbased",0.8);
//Retrieve the Eta restriction on the Relative and Absolute insitu calibration
double insitu_etarestriction_relativeandabsolute = m_config->GetValue("InsituEtaRestrictionRelativeandAbsolute",0.8);
+ // Apply Insitu only to calo and TA mass calibrated jets (only for large jets)
+ m_applyInsituCaloTAjets = m_config->GetValue("ApplyInsituCaloTAJets", false);
+ // Apply in situ JMS:
+ m_applyInsituJMS = m_config->GetValue("ApplyInsituJMS", false);
//Find the absolute path to the insitu root file
if ( !insitu_filename.EqualTo("None") ){
@@ -82,14 +101,14 @@ StatusCode InsituDataCorrection::initializeTool(const std::string&) {
m_relhistoPtMax = rel_histo->GetXaxis()->GetBinLowEdge(rel_histo->GetNbinsX()+1);
m_abshistoPtMax = abs_histo->GetBinLowEdge(abs_histo->GetNbinsX()+1);
// combine in situ calibrations
- m_insituCorr = combineCalibration(rel_histo.get(),abs_histo.get());
+ m_insituCorr = std::unique_ptr<TH2D>(combineCalibration(rel_histo.get(),abs_histo.get()));
m_insituEtaMax = m_insituCorr->GetYaxis()->GetBinLowEdge(m_insituCorr->GetNbinsY()+1);
m_insituPtMin = m_insituCorr->GetXaxis()->GetBinLowEdge(1);
m_insituPtMax = m_insituCorr->GetXaxis()->GetBinLowEdge(m_insituCorr->GetNbinsX()+1);
if(m_applyEtaRestrictionRelativeandAbsolute) m_insituEtaMax = insitu_etarestriction_relativeandabsolute;
}
if(m_applyResidualMCbasedInsitu){
- m_insituCorr_ResidualMCbased = (TH2D*)JetCalibUtils::GetHisto2(insitu_file,residualmcbased_histoname);
+ m_insituCorr_ResidualMCbased = std::unique_ptr<TH2D>((TH2D*)JetCalibUtils::GetHisto2(insitu_file,residualmcbased_histoname));
if ( !m_insituCorr_ResidualMCbased ) {
ATH_MSG_FATAL( "\n Tool configured for the Residual MC based correction, but no residualmcbased in-situ histograms could be retrieved. Aborting..." );
return StatusCode::FAILURE;
@@ -102,12 +121,79 @@ StatusCode InsituDataCorrection::initializeTool(const std::string&) {
}
if(m_applyEtaRestrictionResidualMCbased) m_insituEtaMax_ResidualMCbased = insitu_etarestriction_residualmcbased;
}
- if(!m_applyRelativeandAbsoluteInsitu & !m_applyResidualMCbasedInsitu){
+ if(!m_applyRelativeandAbsoluteInsitu && !m_applyResidualMCbasedInsitu){
ATH_MSG_FATAL( "\n Tool configured for Insitu correction, but no in-situ histograms could be retrieved. Aborting..." );
return StatusCode::FAILURE;
}
- //insitu_file->Close();
+ //Large-R in situ JMS calibration
+ if(m_applyInsituJMS){
+ //Retrieve the name of root files containing the in-situ calibration histograms from the config
+ TString insituJMS_filename = m_config->GetValue("InsituCalibrationFile_JMS","None");
+ //Retrieve the name of the histogram for the absolute in-situ calibration
+ TString abs_histoname_JMS = m_config->GetValue("AbsoluteInsituCalibrationHistogram_JMS","");
+ TString abs_histoname_JMS_TA = m_config->GetValue("AbsoluteInsituCalibrationHistogram_JMS_TA","");
+ //Retrieve the eta range for the in-situ JMS calibration
+ double insitu_etarestriction_JMS = m_config->GetValue("InsituEtaRestriction_JMS",2.0);
+
+ //Find the absolute path to the insitu root file Low and High Mass
+ if ( !insituJMS_filename.EqualTo("None")){
+ if(m_dev){
+ insituJMS_filename.Remove(0,32);
+ insituJMS_filename.Insert(0,"JetCalibTools/");
+ }
+ else{
+ insituJMS_filename.Insert(14,m_calibAreaTag);
+ }
+ insituJMS_filename=PathResolverFindCalibFile(insituJMS_filename.Data());
+ }
+
+ TFile *insituJMS_file = TFile::Open(insituJMS_filename);
+ if ( !insituJMS_file ) { ATH_MSG_FATAL( "Cannot open InsituJMSCalibrationFile: " << insituJMS_filename ); return StatusCode::FAILURE; }
+
+ ATH_MSG_INFO("Reading In-situ JMS correction factors from: " << insituJMS_filename);
+
+ abs_histoname_JMS.ReplaceAll("JETALGO",m_jetAlgo);
+ if(m_applyInsituCaloTAjets){
+ abs_histoname_JMS_TA.ReplaceAll("JETALGO",m_jetAlgo);
+ }
+
+ if(m_applyRelativeandAbsoluteInsitu){
+ std::unique_ptr<TH2D> abs_histo_JMS(dynamic_cast<TH2D*>(JetCalibUtils::GetHisto2(insituJMS_file,abs_histoname_JMS)));
+ if ( !abs_histo_JMS ) {
+ ATH_MSG_FATAL( "\n Tool configured for data, but no in-situ JMS histogram could be retrieved. Aborting..." );
+ return StatusCode::FAILURE;
+ }
+ else {
+ gROOT->cd();
+ // combine in situ calibrations
+ m_insituCorr_JMS = std::unique_ptr<TH2D>(invertHistogram(abs_histo_JMS.get()));
+ m_insituEtaMax_JMS = insitu_etarestriction_JMS;
+ m_insituPtMin_JMS = m_insituCorr_JMS->GetXaxis()->GetBinLowEdge(1);
+ m_insituPtMax_JMS = m_insituCorr_JMS->GetXaxis()->GetBinLowEdge(m_insituCorr_JMS->GetNbinsX()+1);
+ m_insituMassMin_JMS = m_insituCorr_JMS->GetYaxis()->GetBinLowEdge(1);
+ m_insituMassMax_JMS = m_insituCorr_JMS->GetYaxis()->GetBinLowEdge((m_insituCorr_JMS->GetNbinsY()+1));
+
+ if(m_applyInsituCaloTAjets){
+
+ std::unique_ptr<TH2D> abs_histo_JMS_TA(dynamic_cast<TH2D*>(JetCalibUtils::GetHisto2(insituJMS_file,abs_histoname_JMS_TA)));
+
+ if ( !abs_histo_JMS_TA ){
+ ATH_MSG_FATAL( "\n Tool configured for data, but no in-situ JMS histogram for TA mass could be retrieved. Aborting..." );
+ return StatusCode::FAILURE;
+ }
+
+ gROOT->cd();
+ m_insituCorr_JMS_TA = std::unique_ptr<TH2D>(invertHistogram(abs_histo_JMS_TA.get()));
+ }
+ }
+ }
+ if(!m_applyRelativeandAbsoluteInsitu){
+ ATH_MSG_FATAL( "\n Tool configured for Insitu correction, but no in-situ histograms could be retrieved. Aborting..." );
+ return StatusCode::FAILURE;
+ }
+ }
+
ATH_MSG_INFO("Tool configured to calibrate data");
ATH_MSG_INFO("In-situ correction to be applied: " << insitu_desc);
return StatusCode::SUCCESS;
@@ -116,28 +202,117 @@ StatusCode InsituDataCorrection::initializeTool(const std::string&) {
StatusCode InsituDataCorrection::calibrateImpl(xAOD::Jet& jet, JetEventInfo&) const {
- xAOD::JetFourMom_t jetStartP4;
- ATH_CHECK( setStartP4(jet) );
- jetStartP4 = jet.jetP4();
-
float detectorEta = jet.getAttribute<float>("DetectorEta");
- xAOD::JetFourMom_t calibP4=jetStartP4;
+ // For small R jets or large-R jets without calo-TA combination
+ if(!m_applyInsituCaloTAjets){
+ xAOD::JetFourMom_t jetStartP4;
+ ATH_CHECK( setStartP4(jet) );
+ jetStartP4 = jet.jetP4();
+
+ xAOD::JetFourMom_t calibP4=jetStartP4;
- if(m_applyResidualMCbasedInsitu) calibP4=calibP4*getInsituCorr( calibP4.pt(), fabs(detectorEta), "ResidualMCbased" );
+ if(m_applyResidualMCbasedInsitu) calibP4=calibP4*getInsituCorr( calibP4.pt(), fabs(detectorEta), "ResidualMCbased" );
+
+ if(m_dev){
+ float insituFactor = getInsituCorr( jetStartP4.pt(), detectorEta, "RelativeAbs" );
+ jet.setAttribute<float>("JetRelativeAbsInsituCalibFactor",insituFactor);
+ }
- if(m_dev){
- float insituFactor = getInsituCorr( jetStartP4.pt(), detectorEta, "RelativeAbs" );
- jet.setAttribute<float>("JetRelativeAbsInsituCalibFactor",insituFactor);
+ if(m_applyRelativeandAbsoluteInsitu) calibP4=calibP4*getInsituCorr( jetStartP4.pt(), detectorEta, "RelativeAbs" );
+
+ // Only for large R jets with insitu JMS but no combination
+ if(m_applyInsituJMS){
+ xAOD::JetFourMom_t calibP4_JMS;
+ calibP4_JMS = calibP4;
+
+ calibP4_JMS=calibP4*getInsituCorr_JMS( calibP4.pt(), calibP4.M(), detectorEta, "RelativeAbs", false );
+
+ // pT doesn't change while applying in situ JMS
+ TLorentzVector TLVjet;
+ TLVjet.SetPtEtaPhiM( calibP4.pt(), calibP4.eta(), calibP4.phi(), calibP4_JMS.M() );
+ calibP4.SetPxPyPzE( TLVjet.Px(), TLVjet.Py(), TLVjet.Pz(), TLVjet.E() );
+ }
+
+ //Transfer calibrated jet properties to the Jet object
+ jet.setAttribute<xAOD::JetFourMom_t>("JetInsituScaleMomentum",calibP4);
+ jet.setJetP4( calibP4 );
}
- if(m_applyRelativeandAbsoluteInsitu) calibP4=calibP4*getInsituCorr( jetStartP4.pt(), detectorEta, "RelativeAbs" );
+ // For large-R jets: insitu needs to be applied to calo mass and TA mass (by default it's only applied to combined mass)
+ if(m_applyInsituCaloTAjets){
+
+ // calo mass calibrated jets
+ xAOD::JetFourMom_t jetStartP4_calo;
+ xAOD::JetFourMom_t calibP4_calo;
+ if(jet.getAttribute<xAOD::JetFourMom_t>("JetJMSScaleMomentumCalo",jetStartP4_calo)){
+ calibP4_calo=jetStartP4_calo;
+ }else{
+ ATH_MSG_FATAL( "Cannot retrieve JetJMSScaleMomentumCalo jets" );
+ return StatusCode::FAILURE;
+ }
+
+ if(m_applyResidualMCbasedInsitu) calibP4_calo=calibP4_calo*getInsituCorr( calibP4_calo.pt(), fabs(detectorEta), "ResidualMCbased" );
+
+ if(m_dev){
+ float insituFactor_calo = getInsituCorr( jetStartP4_calo.pt(), detectorEta, "RelativeAbs" );
+ jet.setAttribute<float>("JetRelativeAbsInsituCalibFactor_calo",insituFactor_calo);
+ }
+
+ if(m_applyRelativeandAbsoluteInsitu) calibP4_calo=calibP4_calo*getInsituCorr( jetStartP4_calo.pt(), detectorEta, "RelativeAbs" );
+
+ if(m_applyInsituJMS){
+ xAOD::JetFourMom_t calibP4_calo_JMS;
+ calibP4_calo_JMS = calibP4_calo;
+
+ calibP4_calo_JMS=calibP4_calo*getInsituCorr_JMS( calibP4_calo.pt(), calibP4_calo.M(), detectorEta, "RelativeAbs", false );
- // If the jet mass calibration was applied the calibrated mass needs to be used!
+ // pT doesn't change while applying in situ JMS
+ TLorentzVector TLVjet_calo;
+ TLVjet_calo.SetPtEtaPhiM( calibP4_calo.pt(), calibP4_calo.eta(), calibP4_calo.phi(), calibP4_calo_JMS.M() );
+ calibP4_calo.SetPxPyPzE( TLVjet_calo.Px(), TLVjet_calo.Py(), TLVjet_calo.Pz(), TLVjet_calo.E() );
+ }
+
+ //Transfer calibrated jet properties to the Jet object
+ jet.setAttribute<xAOD::JetFourMom_t>("JetInsituScaleMomentumCalo",calibP4_calo);
+ jet.setJetP4( calibP4_calo );
+
+ // TA mass calibrated jets
+ xAOD::JetFourMom_t jetStartP4_ta;
+ xAOD::JetFourMom_t calibP4_ta;
+ if(jet.getAttribute<xAOD::JetFourMom_t>("JetJMSScaleMomentumTA", jetStartP4_ta)){
+ calibP4_ta=jetStartP4_ta;
+ }else{
+ ATH_MSG_FATAL( "Cannot retrieve JetJMSScaleMomentumTA jets" );
+ return StatusCode::FAILURE;
+ }
- //Transfer calibrated jet properties to the Jet object
- jet.setAttribute<xAOD::JetFourMom_t>("JetInsituScaleMomentum",calibP4);
- jet.setJetP4( calibP4 );
+ if(m_applyResidualMCbasedInsitu) calibP4_ta=calibP4_ta*getInsituCorr( calibP4_ta.pt(), fabs(detectorEta), "ResidualMCbased" );
+
+ if(m_dev){
+ float insituFactor_ta = getInsituCorr( jetStartP4_ta.pt(), detectorEta, "RelativeAbs" );
+ jet.setAttribute<float>("JetRelativeAbsInsituCalibFactor_ta",insituFactor_ta);
+ }
+
+ if(m_applyRelativeandAbsoluteInsitu) calibP4_ta=calibP4_ta*getInsituCorr( jetStartP4_ta.pt(), detectorEta, "RelativeAbs" );
+
+ if(m_applyInsituJMS){
+ xAOD::JetFourMom_t calibP4_ta_JMS;
+ calibP4_ta_JMS = calibP4_ta;
+
+ calibP4_ta_JMS=calibP4_ta*getInsituCorr_JMS( calibP4_ta.pt(), calibP4_ta.M(), detectorEta, "RelativeAbs", true );
+
+ // pT doesn't change while applying in situ JMS
+ TLorentzVector TLVjet_ta;
+ TLVjet_ta.SetPtEtaPhiM( calibP4_ta.pt(), calibP4_ta.eta(), calibP4_ta.phi(), calibP4_ta_JMS.M() );
+ calibP4_ta.SetPxPyPzE( TLVjet_ta.Px(), TLVjet_ta.Py(), TLVjet_ta.Pz(), TLVjet_ta.E() );
+ }
+
+ //Transfer calibrated jet properties to the Jet object
+ jet.setAttribute<xAOD::JetFourMom_t>("JetInsituScaleMomentumTA",calibP4_ta);
+ jet.setJetP4( calibP4_ta );
+
+ }
return StatusCode::SUCCESS;
}
@@ -175,6 +350,49 @@ double InsituDataCorrection::getInsituCorr(double pt, double eta, std::string ca
return m_insituCorr->Interpolate(myPt,myEta);
}
+double InsituDataCorrection::getInsituCorr_JMS(double pt, double mass, double eta, std::string calibstep, bool isTAmass) const {
+
+ if(!isTAmass){
+ if (!m_insituCorr_JMS) return 1.0;
+ }
+ else{
+ if (!m_insituCorr_JMS_TA) return 1.0;
+ }
+
+ double myEta = eta, myPt = pt/m_GeV, myMass = mass/m_GeV;
+
+ //eta and pt ranges depends on the insitu calibration
+ double etaMax = m_insituEtaMax_JMS;
+ double ptMin = m_insituPtMin_JMS;
+ double ptMax = m_insituPtMax_JMS;
+ double massMin = m_insituMassMin_JMS;
+ double massMax = m_insituMassMax_JMS;
+
+ //protection against values outside the histogram range, snap back to the lowest/highest bin edge
+ if ( myPt <= ptMin ) myPt = ptMin + 1e-6;
+ else if ( myPt >= ptMax ) myPt = ptMax - 1e-6;
+ if (calibstep == "RelativeAbs" && m_applyEtaRestrictionRelativeandAbsolute) {
+ if(myEta>=etaMax) return 1.0;
+ else if(myEta<=-etaMax) return 1.0;
+ }
+ if (myEta <= -etaMax) myEta = 1e-6 - etaMax;
+ else if (myEta >= etaMax) myEta = etaMax - 1e-6;
+ if (myMass <= massMin ) myMass = massMin + 1e-6;
+ else if (myMass >= massMax ) myMass = massMax - 1e-6;
+
+ double calibFactor = 1.0;
+ if(!isTAmass){
+ calibFactor = m_insituCorr_JMS->Interpolate(myPt,myMass);
+ }
+ else{
+ calibFactor = m_insituCorr_JMS_TA->Interpolate(myPt,myMass);
+ }
+
+ return calibFactor;
+}
+
+
+
TH2D * InsituDataCorrection::combineCalibration(TH2D *h2d, TH1D *h) {
TH2D *prod = (TH2D*)h2d->Clone();
for (int xi=1;xi<=prod->GetNbinsX();xi++) {
@@ -189,3 +407,14 @@ TH2D * InsituDataCorrection::combineCalibration(TH2D *h2d, TH1D *h) {
}
return prod;
}
+
+TH2D * InsituDataCorrection::invertHistogram(TH2D *h2d){
+ TH2D *inv = (TH2D*)h2d->Clone();
+ for (int xi=1;xi<=inv->GetNbinsX();xi++) {
+ for (int yi=1;yi<=inv->GetNbinsY();yi++) {
+ inv->SetBinContent(xi,yi,1./h2d->GetBinContent(xi,yi));
+
+ }
+ }
+ return inv;
+}
diff --git a/Reconstruction/Jet/JetCalibTools/Root/JMSCorrection.cxx b/Reconstruction/Jet/JetCalibTools/Root/JMSCorrection.cxx
index 98144e15cbaac030b44f929efcff33d01279517e..2a5ee556e56eb5883199fc1e0c053bb45a5c0026 100644
--- a/Reconstruction/Jet/JetCalibTools/Root/JMSCorrection.cxx
+++ b/Reconstruction/Jet/JetCalibTools/Root/JMSCorrection.cxx
@@ -1,5 +1,5 @@
/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+ Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
*/
/*
@@ -66,7 +66,7 @@ StatusCode JMSCorrection::initializeTool(const std::string&) {
if ( !m_config ) { ATH_MSG_FATAL("Config file not specified. Aborting."); return StatusCode::FAILURE; }
- m_jetStartScale = m_config->GetValue("JMSStartingScale","JetGSCScaleMomentum");
+ m_jetStartScale = m_config->GetValue("JMSStartingScale","JetEtaJESScaleMomentum");
m_jetOutScale = m_config->GetValue("JMSOutScale","JetJMSScaleMomentum");
// Starting pT value to calibrate
m_pTMinCorr = m_config->GetValue("MinpT",180);
@@ -83,116 +83,116 @@ StatusCode JMSCorrection::initializeTool(const std::string&) {
// *Note* that it is assumed the calo and TA masses use the same histogram dimensionality
m_use3Dhisto = m_config->GetValue("MassCalibrationIs3D",false);
-
-
- //find the ROOT file containing response histograms, path comes from the config file.
- TString JMSFile = m_config->GetValue("MassCalibrationFile","empty");
- if ( JMSFile.EqualTo("empty") ) {
- ATH_MSG_FATAL("NO JMSFactorsFile specified. Aborting.");
- return StatusCode::FAILURE;
- }
- if(m_dev){
- JMSFile.Remove(0,33);
- JMSFile.Insert(0,"JetCalibTools/");
- }
- else{JMSFile.Insert(14,m_calibAreaTag);}
- TString fileName = PathResolverFindCalibFile(JMSFile.Data());
- TFile *inputFile = TFile::Open(fileName);
- if (!inputFile){
- ATH_MSG_FATAL("Cannot open JMS factors file" << fileName);
- return StatusCode::FAILURE;
- }
-
- //ATH_MSG_INFO(" for " << m_jetAlgo << " jets\n\n");
-
- if (!m_use3Dhisto) setMassEtaBins( JetCalibUtils::VectorizeD( m_config->GetValue("MassEtaBins","") ) );
-
- //Get a TList of TKeys pointing to the histograms contained in the ROOT file
- inputFile->cd();
- TList *keys = inputFile->GetListOfKeys();
- //fill the names of the TKeys into a vector of TStrings
- TIter ikeys(keys);
- while ( TKey *iterobj = (TKey*)ikeys() ) {
- TString histoName = iterobj->GetName();
- if ( histoName.Contains(m_jetAlgo) )
- {
- if (m_use3Dhisto)
- m_respFactorMass3D = dynamic_cast<TH3F*>(JetCalibUtils::GetHisto3(inputFile,histoName.Data()));
- else
- m_respFactorsMass.push_back( (TH2F*)JetCalibUtils::GetHisto2(inputFile,histoName.Data()) );
- }
- }
-
- //Make sure we put something in the vector of TH2Fs or we filled the TH3F
- if ( !m_use3Dhisto && m_respFactorsMass.size() < 3 ) {
- ATH_MSG_FATAL("Vector of mass correction histograms may be empty. Please check your mass calibration file: " << JMSFile);
- return StatusCode::FAILURE;
- }
- else if ( m_use3Dhisto && !m_respFactorMass3D)
- {
- ATH_MSG_FATAL("3D mass correction histogram may be missing. Please check your mass calibration file: " << JMSFile);
- return StatusCode::FAILURE;
- }
- else ATH_MSG_INFO("JMS Tool has been initialized with binning and eta fit factors from: " << fileName);
-
- // Shoulb be applied the mass combination?
+ // Should the mass combination be applied?
m_combination = m_config->GetValue("Combination",false); // true: turn on combination of calo mass with track-assisted mass
m_useCorrelatedWeights = m_config->GetValue("UseCorrelatedWeights",false); // true: turn on combination of calo mass with track-assisted mass
+ // Should we only apply the combination (e.g after in situ calibration)
+ m_onlyCombination = m_config->GetValue("OnlyCombination",false);
- // Track-Assisted Jet Mass correction
- m_trackAssistedJetMassCorr = m_config->GetValue("TrackAssistedJetMassCorr",false);
- TString JMS_TrackAssisted_File;
- TString file_trkAssisted_Name;
- TFile *inputFile_trkAssisted;
- if(m_trackAssistedJetMassCorr || m_combination){
- ATH_MSG_INFO("Track Assisted Jet Mass will be calibrated");
- JMS_TrackAssisted_File = m_config->GetValue("TrackAssistedMassCalibrationFile","empty");
- if ( JMS_TrackAssisted_File.EqualTo("empty") ) {
- ATH_MSG_FATAL("NO Track Assisted Mass Factors File specified. Aborting.");
+ //find the ROOT file containing response histograms, path comes from the config file.
+ TString JMSFile;
+
+ if(!m_onlyCombination){
+ JMSFile = m_config->GetValue("MassCalibrationFile","empty");
+ if ( JMSFile.EqualTo("empty") ) {
+ ATH_MSG_FATAL("NO JMSFactorsFile specified. Aborting.");
return StatusCode::FAILURE;
}
if(m_dev){
- JMS_TrackAssisted_File.Remove(0,33);
- JMS_TrackAssisted_File.Insert(0,"JetCalibTools/");
+ JMSFile.Remove(0,33);
+ JMSFile.Insert(0,"JetCalibTools/");
}
- else{JMS_TrackAssisted_File.Insert(14,m_calibAreaTag);}
- file_trkAssisted_Name = PathResolverFindCalibFile(JMS_TrackAssisted_File.Data());
- inputFile_trkAssisted = TFile::Open(file_trkAssisted_Name);
- if (!inputFile_trkAssisted){
- ATH_MSG_FATAL("Cannot open Track Assisted Mass factors file" << fileName);
+ else{JMSFile.Insert(14,m_calibAreaTag);}
+ TString fileName = PathResolverFindCalibFile(JMSFile.Data());
+ TFile *inputFile = TFile::Open(fileName);
+ if (!inputFile){
+ ATH_MSG_FATAL("Cannot open JMS factors file" << fileName);
return StatusCode::FAILURE;
}
-
- //setMassEtaBins( JetCalibUtils::VectorizeD( m_config->GetValue("MassEtaBins","") ) );
-
+
+ if (!m_use3Dhisto) setMassEtaBins( JetCalibUtils::VectorizeD( m_config->GetValue("MassEtaBins","") ) );
+
//Get a TList of TKeys pointing to the histograms contained in the ROOT file
- inputFile_trkAssisted->cd();
- TList *keys_trkAssisted = inputFile_trkAssisted->GetListOfKeys();
+ inputFile->cd();
+ TList *keys = inputFile->GetListOfKeys();
//fill the names of the TKeys into a vector of TStrings
- TIter ikeys_trkAssisted(keys_trkAssisted);
- while ( TKey *iterobj = (TKey*)ikeys_trkAssisted() ) {
+ TIter ikeys(keys);
+ while ( TKey *iterobj = (TKey*)ikeys() ) {
TString histoName = iterobj->GetName();
if ( histoName.Contains(m_jetAlgo) )
- {
- if (m_use3Dhisto)
- m_respFactorTrackAssistedMass3D = dynamic_cast<TH3F*>(JetCalibUtils::GetHisto3(inputFile_trkAssisted,histoName.Data()));
- else
- m_respFactorsTrackAssistedMass.push_back( (TH2F*)JetCalibUtils::GetHisto2(inputFile_trkAssisted,histoName.Data()) );
- }
+ {
+ if (m_use3Dhisto)
+ m_respFactorMass3D = dynamic_cast<TH3F*>(JetCalibUtils::GetHisto3(inputFile,histoName.Data()));
+ else
+ m_respFactorsMass.push_back( (TH2F*)JetCalibUtils::GetHisto2(inputFile,histoName.Data()) );
+ }
}
-
- //Make sure we put something in the vector of TH2Fs
- if ( !m_use3Dhisto && m_respFactorsTrackAssistedMass.size() < 3 ) {
- ATH_MSG_FATAL("Vector of track assisted mass correction histograms may be empty. Please check your track assisted mass calibration file: " << JMSFile);
+
+ //Make sure we put something in the vector of TH2Fs or we filled the TH3F
+ if ( !m_use3Dhisto && m_respFactorsMass.size() < 3 ) {
+ ATH_MSG_FATAL("Vector of mass correction histograms may be empty. Please check your mass calibration file: " << JMSFile);
return StatusCode::FAILURE;
}
- else if ( m_use3Dhisto && !m_respFactorTrackAssistedMass3D)
- {
- ATH_MSG_FATAL("3D track assisted mass correction histogram may be missing. Please check your mass calibration file: " << JMSFile);
- return StatusCode::FAILURE;
+ else if ( m_use3Dhisto && !m_respFactorMass3D)
+ {
+ ATH_MSG_FATAL("3D mass correction histogram may be missing. Please check your mass calibration file: " << JMSFile);
+ return StatusCode::FAILURE;
+ }
+ else ATH_MSG_INFO("JMS Tool has been initialized with binning and eta fit factors from: " << fileName);
+
+ // Track-Assisted Jet Mass correction
+ m_trackAssistedJetMassCorr = m_config->GetValue("TrackAssistedJetMassCorr",false);
+ TString JMS_TrackAssisted_File;
+ TString file_trkAssisted_Name;
+ TFile *inputFile_trkAssisted;
+ if(m_trackAssistedJetMassCorr){
+ ATH_MSG_INFO("Track Assisted Jet Mass will be calibrated");
+ JMS_TrackAssisted_File = m_config->GetValue("TrackAssistedMassCalibrationFile","empty");
+ if ( JMS_TrackAssisted_File.EqualTo("empty") ) {
+ ATH_MSG_FATAL("NO Track Assisted Mass Factors File specified. Aborting.");
+ return StatusCode::FAILURE;
+ }
+ if(m_dev){
+ JMS_TrackAssisted_File.Remove(0,33);
+ JMS_TrackAssisted_File.Insert(0,"JetCalibTools/");
+ }
+ else{JMS_TrackAssisted_File.Insert(14,m_calibAreaTag);}
+ file_trkAssisted_Name = PathResolverFindCalibFile(JMS_TrackAssisted_File.Data());
+ inputFile_trkAssisted = TFile::Open(file_trkAssisted_Name);
+ if (!inputFile_trkAssisted){
+ ATH_MSG_FATAL("Cannot open Track Assisted Mass factors file" << fileName);
+ return StatusCode::FAILURE;
+ }
+
+ //Get a TList of TKeys pointing to the histograms contained in the ROOT file
+ inputFile_trkAssisted->cd();
+ TList *keys_trkAssisted = inputFile_trkAssisted->GetListOfKeys();
+ //fill the names of the TKeys into a vector of TStrings
+ TIter ikeys_trkAssisted(keys_trkAssisted);
+ while ( TKey *iterobj = (TKey*)ikeys_trkAssisted() ) {
+ TString histoName = iterobj->GetName();
+ if ( histoName.Contains(m_jetAlgo) )
+ {
+ if (m_use3Dhisto)
+ m_respFactorTrackAssistedMass3D = dynamic_cast<TH3F*>(JetCalibUtils::GetHisto3(inputFile_trkAssisted,histoName.Data()));
+ else
+ m_respFactorsTrackAssistedMass.push_back( (TH2F*)JetCalibUtils::GetHisto2(inputFile_trkAssisted,histoName.Data()) );
+ }
+ }
+
+ //Make sure we put something in the vector of TH2Fs
+ if ( !m_use3Dhisto && m_respFactorsTrackAssistedMass.size() < 3 ) {
+ ATH_MSG_FATAL("Vector of track assisted mass correction histograms may be empty. Please check your track assisted mass calibration file: " << JMSFile);
+ return StatusCode::FAILURE;
+ }
+ else if ( m_use3Dhisto && !m_respFactorTrackAssistedMass3D)
+ {
+ ATH_MSG_FATAL("3D track assisted mass correction histogram may be missing. Please check your mass calibration file: " << JMSFile);
+ return StatusCode::FAILURE;
+ }
+ else ATH_MSG_INFO("JMS Tool has been initialized with binning and eta fit factors from: " << file_trkAssisted_Name);
}
- else ATH_MSG_INFO("JMS Tool has been initialized with binning and eta fit factors from: " << file_trkAssisted_Name);
- }
+ } //!m_onlyCombination
// Combination
TString Combination_File;
@@ -212,8 +212,8 @@ StatusCode JMSCorrection::initializeTool(const std::string&) {
else{Combination_File.Insert(14,m_calibAreaTag);}
file_combination_Name = PathResolverFindCalibFile(Combination_File.Data());
inputFile_combination = TFile::Open(file_combination_Name);
- if (!inputFile_combination){
- ATH_MSG_FATAL("Cannot open Mass Combination file" << fileName);
+ if (!inputFile_combination && !m_onlyCombination){
+ ATH_MSG_FATAL("Cannot open Mass Combination file" << inputFile_combination);
return StatusCode::FAILURE;
}
@@ -262,33 +262,36 @@ StatusCode JMSCorrection::initializeTool(const std::string&) {
}
//Make sure we put something in the vector of TH2Ds OR filled the TH3s
- if ( !m_use3Dhisto)
- {
- if ( m_caloResolutionMassCombination.size() < 1 ) {
- ATH_MSG_FATAL("Vector of mass combination histograms with calo factors may be empty. Please check your mass combination file: " << JMSFile);
- return StatusCode::FAILURE;
- }
- else if ( m_taResolutionMassCombination.size() < 1 ) {
- ATH_MSG_FATAL("Vector of mass combination histograms with trk-assisted factors may be empty. Please check your mass combination file: " << JMSFile);
- return StatusCode::FAILURE;
- }
- }
- else
- {
- if (!m_caloResolutionMassCombination3D)
- {
- ATH_MSG_FATAL("Mass combination 3D histogram with calo factors was not filled. Please check your mass combination file: " << JMSFile);
- return StatusCode::FAILURE;
- }
- else if (!m_taResolutionMassCombination3D)
- {
- ATH_MSG_FATAL("Mass combination 3D histogram with trk-assisted factors was not filled. Please check your mass combination file: " << JMSFile);
- return StatusCode::FAILURE;
- }
- }
-
+ if(!m_onlyCombination){
+ if ( !m_use3Dhisto)
+ {
+ if ( m_caloResolutionMassCombination.size() < 1 ) {
+ ATH_MSG_FATAL("Vector of mass combination histograms with calo factors may be empty. Please check your mass combination file: " << JMSFile);
+ return StatusCode::FAILURE;
+ }
+ else if ( m_taResolutionMassCombination.size() < 1 ) {
+ ATH_MSG_FATAL("Vector of mass combination histograms with trk-assisted factors may be empty. Please check your mass combination file: " << JMSFile);
+ return StatusCode::FAILURE;
+ }
+ }
+ else
+ {
+ if (!m_caloResolutionMassCombination3D)
+ {
+ ATH_MSG_FATAL("Mass combination 3D histogram with calo factors was not filled. Please check your mass combination file: " << JMSFile);
+ return StatusCode::FAILURE;
+ }
+ else if (!m_taResolutionMassCombination3D)
+ {
+ ATH_MSG_FATAL("Mass combination 3D histogram with trk-assisted factors was not filled. Please check your mass combination file: " << JMSFile);
+ return StatusCode::FAILURE;
+ }
+ }
+ } //m_onlyCombination
+
ATH_MSG_INFO("JMS Tool has been initialized with mass combination weights from: " << file_combination_Name);
- }
+
+ }//m_combination
// Determine the binning strategy
// History is to use pt_mass_eta, with many past config files that don't specify
@@ -547,260 +550,314 @@ StatusCode JMSCorrection::calibrateImpl(xAOD::Jet& jet, JetEventInfo&) const {
xAOD::JetFourMom_t calibP4 = jet.jetP4();
- // For combination
- float mass_ta; // saving calibrated trk-assisted mass
-
float mass_corr = jetStartP4.mass();
double pT_corr = jetStartP4.pt();
- // Determine mass eta bin to use (if using 2D histograms)
- int etabin=-99;
- if (!m_use3Dhisto && (m_massEtaBins.size()==0 || m_respFactorsMass.size() != m_massEtaBins.size()-1)){
- ATH_MSG_FATAL("Please check that the mass correction eta binning is properly set in your config file");
- return StatusCode::FAILURE;
- }
- else if (m_use3Dhisto && !m_respFactorMass3D)
- {
- ATH_MSG_FATAL("Please check that the mass correction 3D histogram is provided");
- return StatusCode::FAILURE;
- }
+ TLorentzVector caloCalibJet;
+ float mass_ta = 0;
- // Originally was jet.getConstituents().size() > 1
- // This essentially requires that the jet has a mass
- // However, constituents are not stored now in rel21 (LCOrigTopoClusters are transient)
- // Thus, getConstituents() breaks unless they are specifically written out
- // Instead, this has been changed to require a non-zero mass
- // Done by S. Schramm on Oct 21, 2017
+ if(!m_onlyCombination){
+ // Determine mass eta bin to use (if using 2D histograms)
+ int etabin=-99;
+ if (!m_use3Dhisto && (m_massEtaBins.size()==0 || m_respFactorsMass.size() != m_massEtaBins.size()-1)){
+ ATH_MSG_FATAL("Please check that the mass correction eta binning is properly set in your config file");
+ return StatusCode::FAILURE;
+ }
+ else if (m_use3Dhisto && !m_respFactorMass3D)
+ {
+ ATH_MSG_FATAL("Please check that the mass correction 3D histogram is provided");
+ return StatusCode::FAILURE;
+ }
- if ( ( ( !m_use3Dhisto && absdetectorEta < m_massEtaBins.back() ) ||
- ( m_use3Dhisto && absdetectorEta < m_respFactorMass3D->GetZaxis()->GetBinLowEdge(m_respFactorMass3D->GetNbinsZ()+1))
- ) && jetStartP4.mass() != 0 ) { // Fiducial cuts
- if (!m_use3Dhisto)
- {
- for (uint i=0; i<m_massEtaBins.size()-1; ++i) {
+ // Originally was jet.getConstituents().size() > 1
+ // This essentially requires that the jet has a mass
+ // However, constituents are not stored now in rel21 (LCOrigTopoClusters are transient)
+ // Thus, getConstituents() breaks unless they are specifically written out
+ // Instead, this has been changed to require a non-zero mass
+ // Done by S. Schramm on Oct 21, 2017
+
+ if ( ( ( !m_use3Dhisto && absdetectorEta < m_massEtaBins.back() ) ||
+ ( m_use3Dhisto && absdetectorEta < m_respFactorMass3D->GetZaxis()->GetBinLowEdge(m_respFactorMass3D->GetNbinsZ()+1))
+ ) && jetStartP4.mass() != 0 ) { // Fiducial cuts
+ if (!m_use3Dhisto)
+ {
+ for (uint i=0; i<m_massEtaBins.size()-1; ++i) {
if(absdetectorEta >= m_massEtaBins[i] && absdetectorEta < m_massEtaBins[i+1]) etabin = i;
- }
- if (etabin< 0){
- ATH_MSG_FATAL("There was a problem determining the eta bin to use for the mass correction");
- return StatusCode::FAILURE;
- }
- }
+ }
+ if (etabin< 0){
+ ATH_MSG_FATAL("There was a problem determining the eta bin to use for the mass correction");
+ return StatusCode::FAILURE;
+ }
+ }
- // Use the correct histogram binning parametrisation when reading the corrected mass
- double massFactor = 1;
- switch (m_binParam)
- {
- case BinningParam::pt_mass_eta:
- if (m_use3Dhisto)
+ // Use the correct histogram binning parametrisation when reading the corrected mass
+ double massFactor = 1;
+ switch (m_binParam)
+ {
+ case BinningParam::pt_mass_eta:
+ if (m_use3Dhisto)
massFactor = getMassCorr3D( jetStartP4.pt()/m_GeV, jetStartP4.mass()/m_GeV, absdetectorEta );
- else
+ else
massFactor = getMassCorr( jetStartP4.pt()/m_GeV, jetStartP4.mass()/m_GeV, etabin );
- break;
- case BinningParam::e_LOGmOe_eta:
- if (m_use3Dhisto)
+ break;
+ case BinningParam::e_LOGmOe_eta:
+ if (m_use3Dhisto)
massFactor = getMassCorr3D( jetStartP4.e()/m_GeV, log(jetStartP4.mass() / jetStartP4.e()), absdetectorEta);
- else
+ else
massFactor = getMassCorr( jetStartP4.e()/m_GeV, log(jetStartP4.mass() / jetStartP4.e()), etabin);
- break;
- case BinningParam::e_LOGmOet_eta:
- if (m_use3Dhisto)
+ break;
+ case BinningParam::e_LOGmOet_eta:
+ if (m_use3Dhisto)
massFactor = getMassCorr3D( jetStartP4.e()/m_GeV, log(jetStartP4.mass() / jetStartP4.Et()), absdetectorEta);
- else
+ else
massFactor = getMassCorr( jetStartP4.e()/m_GeV, log(jetStartP4.mass() / jetStartP4.Et()), etabin);
- break;
- case BinningParam::e_LOGmOpt_eta:
- if (m_use3Dhisto)
+ break;
+ case BinningParam::e_LOGmOpt_eta:
+ if (m_use3Dhisto)
massFactor = getMassCorr3D( jetStartP4.e()/m_GeV, log(jetStartP4.mass() / jetStartP4.pt()), absdetectorEta);
- else
+ else
massFactor = getMassCorr( jetStartP4.e()/m_GeV, log(jetStartP4.mass() / jetStartP4.pt()), etabin);
- break;
- case BinningParam::et_LOGmOet_eta:
- if (m_use3Dhisto)
+ break;
+ case BinningParam::et_LOGmOet_eta:
+ if (m_use3Dhisto)
massFactor = getMassCorr3D( jetStartP4.Et()/m_GeV, log(jetStartP4.mass() / jetStartP4.Et()), absdetectorEta);
- else
+ else
massFactor = getMassCorr( jetStartP4.Et()/m_GeV, log(jetStartP4.mass() / jetStartP4.Et()), etabin);
- break;
- default:
- ATH_MSG_FATAL("This should never be reached - if it happens, it's because a new BinningParam enum option was added, but how to handle it for the calo mass was not. Please contact the tool developer(s) to fix this.");
- return StatusCode::FAILURE;
- break;
+ break;
+ default:
+ ATH_MSG_FATAL("This should never be reached - if it happens, it's because a new BinningParam enum option was added, but how to handle it for the calo mass was not. Please contact the tool developer(s) to fix this.");
+ return StatusCode::FAILURE;
+ break;
+ }
+
+ mass_corr = jetStartP4.mass() / massFactor;
+
+ if(!m_pTfixed) pT_corr = sqrt(jetStartP4.e()*jetStartP4.e()-mass_corr*mass_corr)/cosh( jetStartP4.eta() );
}
- mass_corr = jetStartP4.mass() / massFactor;
+ TLorentzVector caloCalibJet;
+ caloCalibJet.SetPtEtaPhiM(pT_corr, jetStartP4.eta(), jetStartP4.phi(), mass_corr);
- if(!m_pTfixed) pT_corr = sqrt(jetStartP4.e()*jetStartP4.e()-mass_corr*mass_corr)/cosh( jetStartP4.eta() );
- }
-
- TLorentzVector caloCalibJet;
- caloCalibJet.SetPtEtaPhiM(pT_corr, jetStartP4.eta(), jetStartP4.phi(), mass_corr);
-
- if(!m_combination){
- calibP4.SetPxPyPzE( caloCalibJet.Px(), caloCalibJet.Py(), caloCalibJet.Pz(), caloCalibJet.E() );
+ if(!m_combination){
+ calibP4.SetPxPyPzE( caloCalibJet.Px(), caloCalibJet.Py(), caloCalibJet.Pz(), caloCalibJet.E() );
- //Transfer calibrated jet properties to the Jet object
- jet.setAttribute<xAOD::JetFourMom_t>("JetJMSScaleMomentum",calibP4);
- jet.setJetP4( calibP4 );
- }
-
- // Track Assisted Mass Correction
- if(m_trackAssistedJetMassCorr || m_combination ){
-
- double E_corr = jetStartP4.e();
-
- // Determine mass eta bin to use
- if (!m_use3Dhisto)
- {
- etabin=-99;
- if (m_massEtaBins.size()==0 || m_respFactorsTrackAssistedMass.size() != m_massEtaBins.size()-1){
- ATH_MSG_FATAL("Please check that the mass correction eta binning is properly set in your config file");
- if(m_combination) return StatusCode::FAILURE;
- }
- }
- else if (m_use3Dhisto && !m_respFactorTrackAssistedMass3D)
- {
- ATH_MSG_FATAL("Please check that the track assisted mass correction 3D histogram is provided");
- return StatusCode::FAILURE;
+ //Transfer calibrated jet properties to the Jet object
+ jet.setAttribute<xAOD::JetFourMom_t>("JetJMSScaleMomentum",calibP4);
+ jet.setJetP4( calibP4 );
}
- float trackSumMass;
- std::string TrackSumMassStr = "TrackSumMass";
- if(m_jetAlgo=="AntiKt4EMTopo" || m_jetAlgo=="AntiKt4LCTopo") TrackSumMassStr = "DFCommonJets_TrackSumMass";
- std::string TrackSumPtStr = "TrackSumPt";
- if(m_jetAlgo=="AntiKt4EMTopo" || m_jetAlgo=="AntiKt4LCTopo") TrackSumPtStr = "DFCommonJets_TrackSumPt";
- if( !jet.getAttribute<float>(TrackSumMassStr,trackSumMass) ) {
- if(!m_combination){
- //ATH_MSG_WARNING("Failed to retrieve TrackSumMass! Track Assisted Mass Correction will NOT be applied\n\n");
- if(m_warning_counter_mTACorr==0) ATH_MSG_WARNING("Failed to retrieve TrackSumMass! Track Assisted Mass Correction will NOT be applied");
- m_warning_counter_mTACorr++;
- return StatusCode::SUCCESS;
- } else{
- ATH_MSG_FATAL("Failed to retrieve TrackSumMass! Mass Combination can NOT be performed. Aborting.");
- return StatusCode::FAILURE;
+ // Track Assisted Mass Correction
+ if(m_trackAssistedJetMassCorr || m_combination){
+
+ double E_corr = jetStartP4.e();
+
+ // Determine mass eta bin to use
+ if (!m_use3Dhisto)
+ {
+ etabin=-99;
+ if (m_massEtaBins.size()==0 || m_respFactorsTrackAssistedMass.size() != m_massEtaBins.size()-1){
+ ATH_MSG_FATAL("Please check that the mass correction eta binning is properly set in your config file");
+ if(m_combination) return StatusCode::FAILURE;
+ }
+ }
+ else if (m_use3Dhisto && !m_respFactorTrackAssistedMass3D)
+ {
+ ATH_MSG_FATAL("Please check that the track assisted mass correction 3D histogram is provided");
+ return StatusCode::FAILURE;
+ }
+
+ float trackSumMass;
+ std::string TrackSumMassStr = "TrackSumMass";
+ if(m_jetAlgo=="AntiKt4EMTopo" || m_jetAlgo=="AntiKt4LCTopo") TrackSumMassStr = "DFCommonJets_TrackSumMass";
+ std::string TrackSumPtStr = "TrackSumPt";
+ if(m_jetAlgo=="AntiKt4EMTopo" || m_jetAlgo=="AntiKt4LCTopo") TrackSumPtStr = "DFCommonJets_TrackSumPt";
+ if( !jet.getAttribute<float>(TrackSumMassStr,trackSumMass) ) {
+ if(!m_combination){
+ //ATH_MSG_WARNING("Failed to retrieve TrackSumMass! Track Assisted Mass Correction will NOT be applied\n\n");
+ if(m_warning_counter_mTACorr==0) ATH_MSG_WARNING("Failed to retrieve TrackSumMass! Track Assisted Mass Correction will NOT be applied");
+ m_warning_counter_mTACorr++;
+ return StatusCode::SUCCESS;
+ } else{
+ ATH_MSG_FATAL("Failed to retrieve TrackSumMass! Mass Combination can NOT be performed. Aborting.");
+ return StatusCode::FAILURE;
+ }
}
- }
- float trackSumPt;
- if( !jet.getAttribute<float>(TrackSumPtStr,trackSumPt) ) {
- if(!m_combination){
- //ATH_MSG_WARNING("Failed to retrieve TrackSumPt! Track Assisted Mass Correction will NOT be applied\n\n");
- if(m_warning_counter_mTACorr==0) ATH_MSG_WARNING("Failed to retrieve TrackSumPt! Track Assisted Mass Correction will NOT be applied");
- m_warning_counter_mTACorr++;
- return StatusCode::SUCCESS;
- } else{
- ATH_MSG_FATAL("Failed to retrieve TrackSumPt! Mass Combination can NOT be performed. Aborting.");
- return StatusCode::FAILURE;
+ float trackSumPt;
+ if( !jet.getAttribute<float>(TrackSumPtStr,trackSumPt) ) {
+ if(!m_combination){
+ //ATH_MSG_WARNING("Failed to retrieve TrackSumPt! Track Assisted Mass Correction will NOT be applied\n\n");
+ if(m_warning_counter_mTACorr==0) ATH_MSG_WARNING("Failed to retrieve TrackSumPt! Track Assisted Mass Correction will NOT be applied");
+ m_warning_counter_mTACorr++;
+ return StatusCode::SUCCESS;
+ } else{
+ ATH_MSG_FATAL("Failed to retrieve TrackSumPt! Mass Combination can NOT be performed. Aborting.");
+ return StatusCode::FAILURE;
+ }
}
- }
- pT_corr = jetStartP4.pt();
- float mTA;
- if(trackSumPt==0) mTA = 0;
- else{mTA = (jetStartP4.pt()/trackSumPt)*trackSumMass;}
- if(mTA<0) mTA = 0;
- mass_corr = mTA;
-
- if ( ( ( !m_use3Dhisto && absdetectorEta < m_massEtaBins.back() ) ||
- ( m_use3Dhisto && absdetectorEta < m_respFactorMass3D->GetZaxis()->GetBinLowEdge(m_respFactorMass3D->GetNbinsZ()+1))
- ) && jetStartP4.mass() != 0 ) { // Fiducial cuts
- if (!m_use3Dhisto)
- {
- for (uint i=0; i<m_massEtaBins.size()-1; ++i) {
+ pT_corr = jetStartP4.pt();
+ float mTA;
+ if(trackSumPt==0) mTA = 0;
+ else{mTA = (jetStartP4.pt()/trackSumPt)*trackSumMass;}
+ if(mTA<0) mTA = 0;
+ mass_corr = mTA;
+
+ if ( ( ( !m_use3Dhisto && absdetectorEta < m_massEtaBins.back() ) ||
+ ( m_use3Dhisto && absdetectorEta < m_respFactorMass3D->GetZaxis()->GetBinLowEdge(m_respFactorMass3D->GetNbinsZ()+1))
+ ) && jetStartP4.mass() != 0 ) { // Fiducial cuts
+ if (!m_use3Dhisto)
+ {
+ for (uint i=0; i<m_massEtaBins.size()-1; ++i) {
if(absdetectorEta >= m_massEtaBins[i] && absdetectorEta < m_massEtaBins[i+1]) etabin = i;
- }
- if (etabin< 0){
- ATH_MSG_FATAL("There was a problem determining the eta bin to use for the track assisted mass correction");
- return StatusCode::FAILURE;
- }
- }
+ }
+ if (etabin< 0){
+ ATH_MSG_FATAL("There was a problem determining the eta bin to use for the track assisted mass correction");
+ return StatusCode::FAILURE;
+ }
+ }
- double mTAFactor = 1;
-
- if(mTA!=0){ // Read the calibration values from histograms only when this value is non-zero
- // Use the correct histogram binning parametrisation when reading the corrected mass
- switch (m_binParam)
- {
- case BinningParam::pt_mass_eta:
- if (m_use3Dhisto)
- mTAFactor = getTrackAssistedMassCorr3D( jetStartP4.pt()/m_GeV, mTA/m_GeV, absdetectorEta );
- else
- mTAFactor = getTrackAssistedMassCorr( jetStartP4.pt()/m_GeV, mTA/m_GeV, etabin );
- break;
- case BinningParam::e_LOGmOe_eta:
- if (m_use3Dhisto)
- mTAFactor = getTrackAssistedMassCorr3D( jetStartP4.e()/m_GeV, log(mTA / jetStartP4.e()), absdetectorEta);
- else
- mTAFactor = getTrackAssistedMassCorr( jetStartP4.e()/m_GeV, log(mTA / jetStartP4.e()), etabin);
- break;
- case BinningParam::e_LOGmOet_eta:
- if (m_use3Dhisto)
- mTAFactor = getTrackAssistedMassCorr3D( jetStartP4.e()/m_GeV, log(mTA / jetStartP4.Et()), absdetectorEta);
- else
- mTAFactor = getTrackAssistedMassCorr( jetStartP4.e()/m_GeV, log(mTA / jetStartP4.Et()), etabin);
- break;
- case BinningParam::e_LOGmOpt_eta:
- if (m_use3Dhisto)
- mTAFactor = getTrackAssistedMassCorr3D( jetStartP4.e()/m_GeV, log(mTA / jetStartP4.pt()), absdetectorEta);
- else
- mTAFactor = getTrackAssistedMassCorr( jetStartP4.e()/m_GeV, log(mTA / jetStartP4.pt()), etabin);
- break;
- case BinningParam::et_LOGmOet_eta:
- if (m_use3Dhisto)
- mTAFactor = getTrackAssistedMassCorr3D( jetStartP4.Et()/m_GeV, log(mTA / jetStartP4.Et()), absdetectorEta);
- else
- mTAFactor = getTrackAssistedMassCorr( jetStartP4.Et()/m_GeV, log(mTA / jetStartP4.Et()), etabin);
- break;
- default:
- ATH_MSG_FATAL("This should never be reached - if it happens, it's because a new BinningParam enum option was added, but how to handle it for the TA mass was not. Please contact the tool developer(s) to fix this.");
- return StatusCode::FAILURE;
- break;
- }
+ double mTAFactor = 1;
+
+ if(mTA!=0){ // Read the calibration values from histograms only when this value is non-zero
+ // Use the correct histogram binning parametrisation when reading the corrected mass
+ switch (m_binParam)
+ {
+ case BinningParam::pt_mass_eta:
+ if (m_use3Dhisto)
+ mTAFactor = getTrackAssistedMassCorr3D( jetStartP4.pt()/m_GeV, mTA/m_GeV, absdetectorEta );
+ else
+ mTAFactor = getTrackAssistedMassCorr( jetStartP4.pt()/m_GeV, mTA/m_GeV, etabin );
+ break;
+ case BinningParam::e_LOGmOe_eta:
+ if (m_use3Dhisto)
+ mTAFactor = getTrackAssistedMassCorr3D( jetStartP4.e()/m_GeV, log(mTA / jetStartP4.e()), absdetectorEta);
+ else
+ mTAFactor = getTrackAssistedMassCorr( jetStartP4.e()/m_GeV, log(mTA / jetStartP4.e()), etabin);
+ break;
+ case BinningParam::e_LOGmOet_eta:
+ if (m_use3Dhisto)
+ mTAFactor = getTrackAssistedMassCorr3D( jetStartP4.e()/m_GeV, log(mTA / jetStartP4.Et()), absdetectorEta);
+ else
+ mTAFactor = getTrackAssistedMassCorr( jetStartP4.e()/m_GeV, log(mTA / jetStartP4.Et()), etabin);
+ break;
+ case BinningParam::e_LOGmOpt_eta:
+ if (m_use3Dhisto)
+ mTAFactor = getTrackAssistedMassCorr3D( jetStartP4.e()/m_GeV, log(mTA / jetStartP4.pt()), absdetectorEta);
+ else
+ mTAFactor = getTrackAssistedMassCorr( jetStartP4.e()/m_GeV, log(mTA / jetStartP4.pt()), etabin);
+ break;
+ case BinningParam::et_LOGmOet_eta:
+ if (m_use3Dhisto)
+ mTAFactor = getTrackAssistedMassCorr3D( jetStartP4.Et()/m_GeV, log(mTA / jetStartP4.Et()), absdetectorEta);
+ else
+ mTAFactor = getTrackAssistedMassCorr( jetStartP4.Et()/m_GeV, log(mTA / jetStartP4.Et()), etabin);
+ break;
+ default:
+ ATH_MSG_FATAL("This should never be reached - if it happens, it's because a new BinningParam enum option was added, but how to handle it for the TA mass was not. Please contact the tool developer(s) to fix this.");
+ return StatusCode::FAILURE;
+ break;
+ }
+ }
+
+ if(mTAFactor!=0) mass_corr = mTA/mTAFactor;
+ else{
+ ATH_MSG_FATAL("The calibration histogram may have a bad filling bin that is causing mTAFactor to be zero. This value should be different from zero in order to take the ratio. Please contact the tool developer to fix this since the calibration histogram may be corrupted. ");
+ return StatusCode::FAILURE;
+ }
+
+ if(!m_pTfixed) pT_corr = sqrt(jetStartP4.e()*jetStartP4.e()-mass_corr*mass_corr)/cosh( jetStartP4.eta() );
+ else{E_corr = sqrt(jetStartP4.P()*jetStartP4.P()+mass_corr*mass_corr);}
}
-
- if(mTAFactor!=0) mass_corr = mTA/mTAFactor;
else{
- ATH_MSG_FATAL("The calibration histogram may have a bad filling bin that is causing mTAFactor to be zero. This value should be different from zero in order to take the ratio. Please contact the tool developer to fix this since the calibration histogram may be corrupted. ");
- return StatusCode::FAILURE;
+ mTA = 0;
+ mass_corr = 0;
+ if(!m_pTfixed) pT_corr = jetStartP4.e()/cosh( jetStartP4.eta() );
+ else{E_corr = jetStartP4.P();}
}
- if(!m_pTfixed) pT_corr = sqrt(jetStartP4.e()*jetStartP4.e()-mass_corr*mass_corr)/cosh( jetStartP4.eta() );
- else{E_corr = sqrt(jetStartP4.P()*jetStartP4.P()+mass_corr*mass_corr);}
- }
- else{
- mTA = 0;
- mass_corr = 0;
- if(!m_pTfixed) pT_corr = jetStartP4.e()/cosh( jetStartP4.eta() );
- else{E_corr = jetStartP4.P();}
- }
-
- // For combination
- mass_ta = mass_corr;
-
- if(!m_combination){
+ TLorentzVector TACalibJet;
+ xAOD::JetFourMom_t TACalibJet_pTfixed = jet.jetP4();
+ if(!m_pTfixed){
+ TACalibJet.SetPtEtaPhiM(pT_corr, jetStartP4.eta(), jetStartP4.phi(), mass_corr);
+ }else{
+ TACalibJet_pTfixed.SetPxPyPzE( jetStartP4.Px(), jetStartP4.Py(), jetStartP4.Pz(), E_corr );}
+
//Transfer calibrated track assisted mass property to the Jet object
jet.setAttribute<float>("JetTrackAssistedMassUnCalibrated",mTA);
jet.setAttribute<float>("JetTrackAssistedMassCalibrated",mass_corr);
if(!m_pTfixed) jet.setAttribute<float>("JetpTCorrByCalibratedTAMass",pT_corr);
else{jet.setAttribute<float>("JetECorrByCalibratedTAMass",E_corr);}
- }
- if(m_combination){
-
- //Transfer calibrated track assisted mass property to the Jet object
+ //float mass_ta;
+ mass_ta = mass_corr;
+
+ // Store calo and TA calibrated jets separetely to further apply insitu:
+ //Transfer calibrated calo mass property to the Jet object
+ xAOD::JetFourMom_t calibP4_calo = jet.jetP4();
+ calibP4_calo.SetCoordinates( caloCalibJet.Pt(), jetStartP4.eta(), jetStartP4.phi(), caloCalibJet.M() );
+ jet.setAttribute<xAOD::JetFourMom_t>("JetJMSScaleMomentumCalo",calibP4_calo);
+
+ //Transfer calibrated TA mass property to the Jet object
xAOD::JetFourMom_t calibP4_ta = jet.jetP4();
- TLorentzVector TLVjet_ta;
- TLVjet_ta.SetPtEtaPhiM( pT_corr, jetStartP4.eta(), jetStartP4.phi(), mass_corr );
- calibP4_ta.SetPxPyPzE( TLVjet_ta.Px(), TLVjet_ta.Py(), TLVjet_ta.Pz(), TLVjet_ta.E() );
+ if(!m_pTfixed){
+ calibP4_ta.SetCoordinates( TACalibJet.Pt(), jetStartP4.eta(), jetStartP4.phi(), TACalibJet.M() );
+ }else{
+ calibP4_ta.SetPxPyPzE( TACalibJet_pTfixed.Px(), TACalibJet_pTfixed.Py(), TACalibJet_pTfixed.Pz(), TACalibJet_pTfixed.E() );}
jet.setAttribute<xAOD::JetFourMom_t>("JetJMSScaleMomentumTA",calibP4_ta);
+ } //m_trackAssistedJetMassCorr
+ } //!m_onlyCombination
- float mass_comb = caloCalibJet.M(); // combined mass
- double pT_comb = caloCalibJet.Pt();
-
- // if one of the mass is null, use the other one
- if( (mass_comb==0) || (mass_corr==0) ) {
- mass_comb = mass_corr+mass_comb;
+ if(m_combination){
+ float mass_calo;
+ float Mass_comb = 0.;
+ double pT_calo;
+ double E_calo;
+ double Et_calo;
+
+ if(m_onlyCombination){
+ // Read input values (calo and TA insitu calibrated jets) for combination:
+
+ xAOD::JetFourMom_t jetInsituP4_calo;
+ xAOD::JetFourMom_t calibP4Insitu_calo;
+ if(jet.getAttribute<xAOD::JetFourMom_t>("JetInsituScaleMomentumCalo",jetInsituP4_calo)){
+ calibP4Insitu_calo=jetInsituP4_calo;
+ }else{
+ ATH_MSG_FATAL( "Cannot retrieve JetInsituScaleMomentumCalo jets" );
+ return StatusCode::FAILURE;
}
- else {
- // Determine mass combination eta bin to use
- int etabin=-99;
- if (!m_use3Dhisto)
+ TLorentzVector TLVCaloInsituCalib;
+ TLVCaloInsituCalib.SetPtEtaPhiM(calibP4Insitu_calo.pt(), calibP4Insitu_calo.eta(), calibP4Insitu_calo.phi(), calibP4Insitu_calo.mass());
+ mass_calo = TLVCaloInsituCalib.M();
+ pT_calo = TLVCaloInsituCalib.Pt();
+ E_calo = TLVCaloInsituCalib.E();
+ Et_calo = TLVCaloInsituCalib.Et();
+
+ xAOD::JetFourMom_t jetInsituP4_ta;
+ xAOD::JetFourMom_t calibP4Insitu_ta;
+ if(jet.getAttribute<xAOD::JetFourMom_t>("JetInsituScaleMomentumTA",jetInsituP4_ta)){
+ calibP4Insitu_ta=jetInsituP4_ta;
+ }else{
+ ATH_MSG_FATAL( "Cannot retrieve JetInsituScaleMomentumTA jets" );
+ return StatusCode::FAILURE;
+ }
+ TLorentzVector TLVTAInsituCalib;
+ TLVTAInsituCalib.SetPtEtaPhiM(calibP4Insitu_ta.pt(), calibP4Insitu_ta.eta(), calibP4Insitu_ta.phi(), calibP4Insitu_ta.mass());
+ mass_ta = TLVTAInsituCalib.M();
+ }else{
+ mass_calo = caloCalibJet.M(); // combined mass
+ pT_calo = caloCalibJet.Pt();
+ E_calo = caloCalibJet.E();
+ Et_calo = caloCalibJet.Et();
+ // mass_ta already defined above
+ }
+
+ // if one of the mass is null, use the other one
+ if( (mass_calo==0) || (mass_ta==0) ) {
+ Mass_comb = mass_ta+mass_calo;
+ }
+ else {
+ // Determine mass combination eta bin to use
+ int etabin=-99;
+ if (!m_use3Dhisto)
{
if (m_massCombinationEtaBins.size()==0 || m_caloResolutionMassCombination.size() != m_massCombinationEtaBins.size()-1){
ATH_MSG_FATAL("Please check that the mass combination eta binning is properly set in your config file");
@@ -811,23 +868,21 @@ StatusCode JMSCorrection::calibrateImpl(xAOD::Jet& jet, JetEventInfo&) const {
return StatusCode::FAILURE;
}
}
- else if (m_use3Dhisto && !m_caloResolutionMassCombination3D)
+ else if (m_use3Dhisto && !m_caloResolutionMassCombination3D)
{
ATH_MSG_FATAL("Please check that the mass resolution 3D histogram is provided");
return StatusCode::FAILURE;
}
- else if (m_use3Dhisto && !m_taResolutionMassCombination3D)
+ else if (m_use3Dhisto && !m_taResolutionMassCombination3D)
{
ATH_MSG_FATAL("Please check that the track assisted mass resolution 3D histogram is provided");
return StatusCode::FAILURE;
}
-
-
- if ( ( ( !m_use3Dhisto && absdetectorEta < m_massCombinationEtaBins.back() ) ||
- ( m_use3Dhisto && absdetectorEta < m_caloResolutionMassCombination3D->GetZaxis()->GetBinLowEdge(m_caloResolutionMassCombination3D->GetNbinsZ()+1)) ) ) {
-
- if (!m_use3Dhisto)
+ if ( ( ( !m_use3Dhisto && absdetectorEta < m_massCombinationEtaBins.back() ) ||
+ ( m_use3Dhisto && absdetectorEta < m_caloResolutionMassCombination3D->GetZaxis()->GetBinLowEdge(m_caloResolutionMassCombination3D->GetNbinsZ()+1)) ) ) {
+
+ if (!m_use3Dhisto)
{
for (uint i=0; i<m_massCombinationEtaBins.size()-1; ++i) {
if(absdetectorEta >= m_massCombinationEtaBins[i] && absdetectorEta < m_massCombinationEtaBins[i+1]) etabin = i;
@@ -838,130 +893,129 @@ StatusCode JMSCorrection::calibrateImpl(xAOD::Jet& jet, JetEventInfo&) const {
}
}
-
- // Use the correct histogram binning parametrisation when reading the combined mass weights
- double relCalo = 0;
- double relTA = 0;
- double rho = 0;
- switch (m_binParam)
+ // Use the correct histogram binning parametrisation when reading the combined mass weights
+ double relCalo = 0;
+ double relTA = 0;
+ double rho = 0;
+ switch (m_binParam)
{
- case BinningParam::pt_mass_eta:
- if (m_use3Dhisto)
+ case BinningParam::pt_mass_eta:
+ if (m_use3Dhisto)
{
- relCalo = getRelCalo3D( caloCalibJet.Pt()/m_GeV, caloCalibJet.M()/caloCalibJet.Pt(), absdetectorEta );
- relTA = getRelTA3D( caloCalibJet.Pt()/m_GeV, mass_ta /caloCalibJet.Pt(), absdetectorEta );
+ relCalo = getRelCalo3D( pT_calo/m_GeV, mass_calo/pT_calo, absdetectorEta );
+ relTA = getRelTA3D( pT_calo/m_GeV, mass_ta/pT_calo, absdetectorEta );
if (m_useCorrelatedWeights)
- rho = getRho3D( caloCalibJet.Pt()/m_GeV, caloCalibJet.M()/caloCalibJet.Pt(), absdetectorEta );
+ rho = getRho3D( pT_calo/m_GeV, mass_calo/pT_calo, absdetectorEta );
}
- else
+ else
{
- relCalo = getRelCalo( caloCalibJet.Pt()/m_GeV, caloCalibJet.M()/caloCalibJet.Pt(), etabin );
- relTA = getRelTA( caloCalibJet.Pt()/m_GeV, mass_ta /caloCalibJet.Pt(), etabin );
+ relCalo = getRelCalo( pT_calo/m_GeV, mass_calo/pT_calo, etabin );
+ relTA = getRelTA( pT_calo/m_GeV, mass_ta/pT_calo, etabin );
if (m_useCorrelatedWeights)
- rho = getRho( caloCalibJet.Pt()/m_GeV, caloCalibJet.M()/caloCalibJet.Pt(), etabin );
+ rho = getRho( pT_calo/m_GeV, mass_calo/pT_calo, etabin );
}
- break;
- case BinningParam::e_LOGmOe_eta:
- if (m_use3Dhisto)
+ break;
+ case BinningParam::e_LOGmOe_eta:
+ if (m_use3Dhisto)
{
- relCalo = getRelCalo3D( caloCalibJet.E()/m_GeV, log(caloCalibJet.M()/caloCalibJet.E()), absdetectorEta );
- relTA = getRelTA3D( caloCalibJet.E()/m_GeV, log(mass_ta /caloCalibJet.E()), absdetectorEta );
+ relCalo = getRelCalo3D( E_calo/m_GeV, log(mass_calo/E_calo), absdetectorEta );
+ relTA = getRelTA3D( E_calo/m_GeV, log(mass_ta/E_calo), absdetectorEta );
if (m_useCorrelatedWeights)
- rho = getRho3D( caloCalibJet.E()/m_GeV, log(caloCalibJet.M()/caloCalibJet.E()), absdetectorEta );
+ rho = getRho3D( E_calo/m_GeV, log(mass_calo/E_calo), absdetectorEta );
}
- else
+ else
{
- relCalo = getRelCalo( caloCalibJet.E()/m_GeV, log(caloCalibJet.M()/caloCalibJet.E()), etabin );
- relTA = getRelTA( caloCalibJet.E()/m_GeV, log(mass_ta /caloCalibJet.E()), etabin );
+ relCalo = getRelCalo( E_calo/m_GeV, log(mass_calo/E_calo), etabin );
+ relTA = getRelTA( E_calo/m_GeV, log(mass_ta/E_calo), etabin );
if (m_useCorrelatedWeights)
- rho = getRho( caloCalibJet.E()/m_GeV, log(caloCalibJet.M()/caloCalibJet.E()), etabin );
+ rho = getRho( E_calo/m_GeV, log(mass_calo/E_calo), etabin );
}
- break;
- case BinningParam::e_LOGmOet_eta:
- if (m_use3Dhisto)
+ break;
+ case BinningParam::e_LOGmOet_eta:
+ if (m_use3Dhisto)
{
- relCalo = getRelCalo3D( caloCalibJet.E()/m_GeV, log(caloCalibJet.M()/caloCalibJet.Et()), absdetectorEta );
- relTA = getRelTA3D( caloCalibJet.E()/m_GeV, log(mass_ta /caloCalibJet.Et()), absdetectorEta );
+ relCalo = getRelCalo3D( E_calo/m_GeV, log(mass_calo/Et_calo), absdetectorEta );
+ relTA = getRelTA3D( E_calo/m_GeV, log(mass_ta/Et_calo), absdetectorEta );
if (m_useCorrelatedWeights)
- rho = getRho3D( caloCalibJet.E()/m_GeV, log(caloCalibJet.M()/caloCalibJet.Et()), absdetectorEta );
+ rho = getRho3D( E_calo/m_GeV, log(mass_calo/Et_calo), absdetectorEta );
}
- else
+ else
{
- relCalo = getRelCalo( caloCalibJet.E()/m_GeV, log(caloCalibJet.M()/caloCalibJet.Et()), etabin );
- relTA = getRelTA( caloCalibJet.E()/m_GeV, log(mass_ta /caloCalibJet.Et()), etabin );
+ relCalo = getRelCalo( E_calo/m_GeV, log(mass_calo/Et_calo), etabin );
+ relTA = getRelTA( E_calo/m_GeV, log(mass_ta/Et_calo), etabin );
if (m_useCorrelatedWeights)
- rho = getRho( caloCalibJet.E()/m_GeV, log(caloCalibJet.M()/caloCalibJet.Et()), etabin );
+ rho = getRho( E_calo/m_GeV, log(mass_calo/Et_calo), etabin );
}
- break;
- case BinningParam::e_LOGmOpt_eta:
- if (m_use3Dhisto)
+ break;
+ case BinningParam::e_LOGmOpt_eta:
+ if (m_use3Dhisto)
{
- relCalo = getRelCalo3D( caloCalibJet.E()/m_GeV, log(caloCalibJet.M()/caloCalibJet.Pt()), absdetectorEta );
- relTA = getRelTA3D( caloCalibJet.E()/m_GeV, log(mass_ta /caloCalibJet.Pt()), absdetectorEta );
+ relCalo = getRelCalo3D( E_calo/m_GeV, log(mass_calo/pT_calo), absdetectorEta );
+ relTA = getRelTA3D( E_calo/m_GeV, log(mass_ta/pT_calo), absdetectorEta );
if (m_useCorrelatedWeights)
- rho = getRho3D( caloCalibJet.E()/m_GeV, log(caloCalibJet.M()/caloCalibJet.Pt()), absdetectorEta );
+ rho = getRho3D( E_calo/m_GeV, log(mass_calo/pT_calo), absdetectorEta );
}
- else
+ else
{
- relCalo = getRelCalo( caloCalibJet.E()/m_GeV, log(caloCalibJet.M()/caloCalibJet.Pt()), etabin );
- relTA = getRelTA( caloCalibJet.E()/m_GeV, log(mass_ta /caloCalibJet.Pt()), etabin );
+ relCalo = getRelCalo( E_calo/m_GeV, log(mass_calo/pT_calo), etabin );
+ relTA = getRelTA( E_calo/m_GeV, log(mass_ta/pT_calo), etabin );
if (m_useCorrelatedWeights)
- rho = getRho( caloCalibJet.E()/m_GeV, log(caloCalibJet.M()/caloCalibJet.Pt()), etabin );
+ rho = getRho( E_calo/m_GeV, log(mass_calo/pT_calo), etabin );
}
- break;
- case BinningParam::et_LOGmOet_eta:
- if (m_use3Dhisto)
+ break;
+ case BinningParam::et_LOGmOet_eta:
+ if (m_use3Dhisto)
{
- relCalo = getRelCalo3D( caloCalibJet.Et()/m_GeV, log(caloCalibJet.M()/caloCalibJet.Et()), absdetectorEta );
- relTA = getRelTA3D( caloCalibJet.Et()/m_GeV, log(mass_ta /caloCalibJet.Et()), absdetectorEta );
+ relCalo = getRelCalo3D( Et_calo/m_GeV, log(mass_calo/Et_calo), absdetectorEta );
+ relTA = getRelTA3D( Et_calo/m_GeV, log(mass_ta/Et_calo), absdetectorEta );
if (m_useCorrelatedWeights)
- rho = getRho3D( caloCalibJet.Et()/m_GeV, log(caloCalibJet.M()/caloCalibJet.Et()), absdetectorEta );
+ rho = getRho3D( Et_calo/m_GeV, log(mass_calo/Et_calo), absdetectorEta );
}
- else
+ else
{
- relCalo = getRelCalo( caloCalibJet.Et()/m_GeV, log(caloCalibJet.M()/caloCalibJet.Et()), etabin );
- relTA = getRelTA( caloCalibJet.Et()/m_GeV, log(mass_ta /caloCalibJet.Et()), etabin );
+ relCalo = getRelCalo( Et_calo/m_GeV, log(mass_calo/Et_calo), etabin );
+ relTA = getRelTA( Et_calo/m_GeV, log(mass_ta/Et_calo), etabin );
if (m_useCorrelatedWeights)
- rho = getRho( caloCalibJet.Et()/m_GeV, log(caloCalibJet.M()/caloCalibJet.Et()), etabin );
+ rho = getRho( Et_calo/m_GeV, log(mass_calo/Et_calo), etabin );
}
- break;
- default:
- ATH_MSG_FATAL("This should never be reached - if it happens, it's because a new BinningParam enum option was added, but how to handle it for the TA mass was not. Please contact the tool developer(s) to fix this.");
- return StatusCode::FAILURE;
- break;
+ break;
+ default:
+ ATH_MSG_FATAL("This should never be reached - if it happens, it's because a new BinningParam enum option was added, but how to handle it for the TA mass was not. Please contact the tool developer(s) to fix this.");
+ return StatusCode::FAILURE;
+ break;
}
-
-
- // Watch for division by zero
- if(m_useCorrelatedWeights && (relCalo*relCalo + relTA*relTA - 2 * rho* relCalo * relTA == 0)){
- ATH_MSG_ERROR("Encountered division by zero when calculating mass combination weight using correlated weights");
- return StatusCode::FAILURE;
- }
- const double Weight = ( relTA*relTA - rho *relCalo*relTA ) / ( relCalo*relCalo + relTA*relTA - 2 * rho* relCalo * relTA );
- mass_comb = ( caloCalibJet.M() * Weight ) + ( mass_ta * ( 1 - Weight) );
- // Protection
- if(mass_comb>jetStartP4.e()) mass_comb = caloCalibJet.M();
- else if(!m_pTfixed) pT_comb = sqrt(jetStartP4.e()*jetStartP4.e()-mass_comb*mass_comb)/cosh( jetStartP4.eta() );
- }
+ // Watch for division by zero
+ if(m_useCorrelatedWeights && (relCalo*relCalo + relTA*relTA - 2 * rho* relCalo * relTA == 0)){
+ ATH_MSG_ERROR("Encountered division by zero when calculating mass combination weight using correlated weights");
+ return StatusCode::FAILURE;
+ }
+ const double Weight = ( relTA*relTA - rho *relCalo*relTA ) / ( relCalo*relCalo + relTA*relTA - 2 * rho* relCalo * relTA );
+
+ // Zero should be only returned by resolution functions if jet mass is negative
+ if(relCalo == 0 && relTA == 0)
+ Mass_comb = 0;
+ else if(relCalo == 0)
+ Mass_comb = mass_ta;
+ else if(relTA == 0)
+ Mass_comb = mass_calo;
+ else
+ Mass_comb = ( mass_calo * Weight ) + ( mass_ta * ( 1 - Weight) );
+ // Protection
+ if(Mass_comb>jetStartP4.e()) Mass_comb = mass_calo;
+ else if(!m_pTfixed) pT_calo = sqrt(jetStartP4.e()*jetStartP4.e()-mass_calo*mass_calo)/cosh( jetStartP4.eta() );
}
-
-
- TLorentzVector TLVjet;
- TLVjet.SetPtEtaPhiM( pT_comb, jetStartP4.eta(), jetStartP4.phi(), mass_comb );
- calibP4.SetPxPyPzE( TLVjet.Px(), TLVjet.Py(), TLVjet.Pz(), TLVjet.E() );
+ }
- //Transfer calibrated jet properties to the Jet object
- jet.setAttribute<xAOD::JetFourMom_t>(m_jetOutScale.Data(),calibP4);
- jet.setJetP4( calibP4 );
-
- //Transfer calibrated calo mass property to the Jet object
- xAOD::JetFourMom_t calibP4_calo = jet.jetP4();
- calibP4_calo.SetCoordinates( caloCalibJet.Pt(), jetStartP4.eta(), jetStartP4.phi(), caloCalibJet.M() );
- jet.setAttribute<xAOD::JetFourMom_t>("JetJMSScaleMomentumCalo",calibP4_calo);
+ TLorentzVector TLVjet;
+ TLVjet.SetPtEtaPhiM( pT_calo, jetStartP4.eta(), jetStartP4.phi(), Mass_comb );
+ calibP4.SetPxPyPzE( TLVjet.Px(), TLVjet.Py(), TLVjet.Pz(), TLVjet.E() );
+
+ //Transfer calibrated jet properties to the Jet object
+ jet.setAttribute<xAOD::JetFourMom_t>(m_jetOutScale.Data(),calibP4);
+ jet.setJetP4( calibP4 );
- }
- }
+ } //m_combination
return StatusCode::SUCCESS;
diff --git a/Reconstruction/Jet/JetCalibTools/Root/JetCalibrationTool.cxx b/Reconstruction/Jet/JetCalibTools/Root/JetCalibrationTool.cxx
index 0fd60ef6b9f8237daab1c9c927901b513fdae429..3b96f6fba078a0440872480644531be874794ce1 100644
--- a/Reconstruction/Jet/JetCalibTools/Root/JetCalibrationTool.cxx
+++ b/Reconstruction/Jet/JetCalibTools/Root/JetCalibrationTool.cxx
@@ -24,7 +24,8 @@ JetCalibrationTool::JetCalibrationTool(const std::string& name)
m_rhkRhoKey(""),
m_jetAlgo(""), m_config(""), m_calibSeq(""), m_calibAreaTag(""), m_originScale(""), m_devMode(false), m_isData(true), m_timeDependentCalib(false), m_rhoKey("auto"), m_dir(""), m_eInfoName(""), m_globalConfig(NULL),
m_doJetArea(true), m_doResidual(true), m_doOrigin(true), m_doGSC(true), m_gscDepth("auto"),
- m_jetPileupCorr(NULL), m_etaJESCorr(NULL), m_globalSequentialCorr(NULL), m_insituDataCorr(NULL), m_jetMassCorr(NULL), m_jetSmearCorr(NULL)
+ m_jetPileupCorr(NULL), m_etaJESCorr(NULL), m_globalSequentialCorr(NULL), m_insituDataCorr(NULL),
+ m_jetMassCorr(NULL), m_jetSmearCorr(NULL), m_insituCombMassCorr_tmp(NULL)
{
declareProperty( "JetCollection", m_jetAlgo = "AntiKt4LCTopo" );
@@ -54,6 +55,7 @@ JetCalibrationTool::~JetCalibrationTool() {
if (m_insituDataCorr) delete m_insituDataCorr;
if (m_jetMassCorr) delete m_jetMassCorr;
if (m_jetSmearCorr) delete m_jetSmearCorr;
+ if (m_insituCombMassCorr_tmp) delete m_insituCombMassCorr_tmp;
}
@@ -215,6 +217,26 @@ StatusCode JetCalibrationTool::initializeTool(const std::string& name) {
}
}
+ //Combined Mass Calibration:
+ m_insituCombMassCalib = m_globalConfig->GetValue("InsituCombinedMassCorrection",false);
+ if(m_insituCombMassCalib && calibSeq.Contains("InsituCombinedMass")){ // Read Combination Config
+ m_insituCombMassConfig = JetCalibUtils::Vectorize( m_globalConfig->GetValue("InsituCombinedMassCorrectionFile","") );
+ if(m_insituCombMassConfig.size()==0) ATH_MSG_ERROR("Please check there is a combination config");
+ for(unsigned int i=0;i<m_insituCombMassConfig.size();++i){
+
+ std::string configPath_comb = dir+m_insituCombMassConfig.at(i).Data(); // Full path
+ TString fn_comb = PathResolverFindCalibFile(configPath_comb);
+
+ ATH_MSG_INFO("Reading combination settings from: " << m_insituCombMassConfig.at(i));
+ ATH_MSG_INFO("resolved in: " << fn_comb);
+
+ TEnv *globalInsituCombMass = new TEnv();
+ int status = globalInsituCombMass->ReadFile(fn_comb ,EEnvLevel(0));
+ if (status!=0) { ATH_MSG_FATAL("Cannot read config file " << fn_comb ); return StatusCode::FAILURE; }
+ m_globalInsituCombMassConfig.push_back(globalInsituCombMass);
+ }
+ }
+
//Loop over the request calib sequence
//Initialize derived classes for applying the requested calibrations and add them to a vector
std::vector<TString> vecCalibSeq = JetCalibUtils::Vectorize(calibSeq,"_");
@@ -303,6 +325,21 @@ StatusCode JetCalibrationTool::getCalibClass(const std::string&name, TString cal
m_calibClasses.push_back(m_jetMassCorr);
return StatusCode::SUCCESS;
}
+ } else if ( calibration.EqualTo("InsituCombinedMass") ){
+ ATH_MSG_INFO("Initializing Combined Mass Correction");
+ for(unsigned int i=0;i<m_insituCombMassConfig.size();++i){
+ suffix="_InsituCombinedMass"; suffix += "_"; suffix += std::to_string(i);
+ if(m_devMode) suffix+="_DEV";
+ m_insituCombMassCorr_tmp = new JMSCorrection(name+suffix,m_globalInsituCombMassConfig.at(i),jetAlgo,calibPath,m_devMode);
+ m_insituCombMassCorr_tmp->msg().setLevel( this->msg().level() );
+ if ( m_insituCombMassCorr_tmp->initializeTool(name+suffix).isFailure() ) {
+ ATH_MSG_FATAL("Couldn't initialize the Combined Mass correction. Aborting");
+ return StatusCode::FAILURE;
+ } else {
+ m_insituCombMassCorr.push_back(m_insituCombMassCorr_tmp);
+ return StatusCode::SUCCESS;
+ }
+ }
} else if ( calibration.EqualTo("Insitu") ) {
if(!m_timeDependentCalib){
ATH_MSG_INFO("Initializing Insitu correction.");
@@ -576,6 +613,17 @@ StatusCode JetCalibrationTool::calibrateImpl(xAOD::Jet& jet, JetEventInfo& jetEv
if(runNumber>m_runBins.at(i) && runNumber<=m_runBins.at(i+1)){ ATH_CHECK ( m_insituTimeDependentCorr.at(i)->calibrateImpl(jet,jetEventInfo) );}
}
}
+ if(CalibSeq.Contains("InsituCombinedMass") && m_insituCombMassCalib){
+ for(unsigned int i=0;i<m_insituCombMassConfig.size();++i){
+ //Retrive EventInfo container
+ const xAOD::EventInfo* eventInfo(nullptr);
+ if( evtStore()->retrieve(eventInfo,"EventInfo").isFailure() || !eventInfo ) {
+ ATH_MSG_ERROR(" JetCalibrationTool::calibrateImpl : Failed to retrieve EventInfo.");
+ }
+ ATH_CHECK ( m_insituCombMassCorr.at(i)->calibrateImpl(jet,jetEventInfo) );
+ }
+ }
+
return StatusCode::SUCCESS;
}