diff --git a/Reconstruction/Jet/JetCalibTools/JetCalibTools/CalibrationMethods/JetSmearingCorrection.h b/Reconstruction/Jet/JetCalibTools/JetCalibTools/CalibrationMethods/JetSmearingCorrection.h
index 1c21a577f738ec7691661508e1edd04486a0d9f9..4e256fd01be1aea8ceceb786a75c0786940127c9 100644
--- a/Reconstruction/Jet/JetCalibTools/JetCalibTools/CalibrationMethods/JetSmearingCorrection.h
+++ b/Reconstruction/Jet/JetCalibTools/JetCalibTools/CalibrationMethods/JetSmearingCorrection.h
@@ -10,6 +10,7 @@
#include "JetCalibTools/JetCalibrationToolBase.h"
#include <memory>
+#include <vector>
#include "TRandom3.h"
@@ -35,10 +36,10 @@ class JetSmearingCorrection
private:
// Helper methods
StatusCode getSigmaSmear(xAOD::Jet& jet, double& sigmaSmear) const;
- StatusCode getNominalResolution(const xAOD::Jet& jet, TH1* histo, double& resolution) const;
+ StatusCode getNominalResolution(const xAOD::Jet& jet, TH1* histo, const std::vector< std::unique_ptr<TH1> >& projections, double& resolution) const;
StatusCode readHisto(double& returnValue, TH1* histo, double x) const;
- StatusCode readHisto(double& returnValue, TH1* histo, double x, double y) const;
- StatusCode readHisto(double& returnValue, TH1* histo, double x, double y, double z) const;
+ StatusCode readHisto(double& returnValue, TH1* histo, const std::vector< std::unique_ptr<TH1> >& projections, double x, double y) const;
+ StatusCode cacheProjections(TH1* fullHistogram, std::vector< std::unique_ptr<TH1> >& cacheLocation, const std::string& type);
// Private enums
enum class SmearType
@@ -79,6 +80,10 @@ class JetSmearingCorrection
InterpType m_interpType;
std::unique_ptr<TH1> m_smearResolutionMC;
std::unique_ptr<TH1> m_smearResolutionData;
+
+ // Variables to cache projections in case of 1-D interpolaton in 2-D or 3-D histograms
+ std::vector< std::unique_ptr<TH1> > m_cachedProjResMC;
+ std::vector< std::unique_ptr<TH1> > m_cachedProjResData;
};
diff --git a/Reconstruction/Jet/JetCalibTools/Root/JetSmearingCorrection.cxx b/Reconstruction/Jet/JetCalibTools/Root/JetSmearingCorrection.cxx
index 36bf9af353e804f3e8a2e3c75bec93cdf4c6bcd1..5e43783136a08e6c4292234d145d49c8d8e30343 100644
--- a/Reconstruction/Jet/JetCalibTools/Root/JetSmearingCorrection.cxx
+++ b/Reconstruction/Jet/JetCalibTools/Root/JetSmearingCorrection.cxx
@@ -20,6 +20,8 @@ JetSmearingCorrection::JetSmearingCorrection(const std::string name)
, m_interpType(InterpType::UNKNOWN)
, m_smearResolutionMC()
, m_smearResolutionData()
+ , m_cachedProjResMC()
+ , m_cachedProjResData()
{ }
JetSmearingCorrection::JetSmearingCorrection(const std::string& name, TEnv* config, TString jetAlgo, TString calibAreaTag, bool dev)
@@ -36,6 +38,8 @@ JetSmearingCorrection::JetSmearingCorrection(const std::string& name, TEnv* conf
, m_interpType(InterpType::UNKNOWN)
, m_smearResolutionMC()
, m_smearResolutionData()
+ , m_cachedProjResMC()
+ , m_cachedProjResData()
{ }
JetSmearingCorrection::~JetSmearingCorrection()
@@ -215,6 +219,15 @@ StatusCode JetSmearingCorrection::initializeTool(const std::string&)
return StatusCode::FAILURE;
}
+ // Pre-cache the histogram file in 1D projections if relevant (depends on InterpType)
+ if (m_interpType == InterpType::OnlyX || m_interpType == InterpType::OnlyY)
+ {
+ if (cacheProjections(m_smearResolutionMC.get(), m_cachedProjResMC,"mc").isFailure())
+ return StatusCode::FAILURE;
+ if (cacheProjections(m_smearResolutionData.get(),m_cachedProjResData,"data").isFailure())
+ return StatusCode::FAILURE;
+ }
+
return StatusCode::SUCCESS;
}
@@ -262,7 +275,7 @@ StatusCode JetSmearingCorrection::readHisto(double& returnValue, TH1* histo, dou
return StatusCode::SUCCESS;
}
-StatusCode JetSmearingCorrection::readHisto(double& returnValue, TH1* histo, double x, double y) const
+StatusCode JetSmearingCorrection::readHisto(double& returnValue, TH1* histo, const std::vector< std::unique_ptr<TH1> >& projections, double x, double y) const
{
// Ensure that the histogram exists
if (!histo)
@@ -293,10 +306,6 @@ StatusCode JetSmearingCorrection::readHisto(double& returnValue, TH1* histo, dou
y = minY + 1.e-6;
// Get the result, interpolating as appropriate
- // We need to use this as a 2D histogram to have access to the projection methods
- // We also need a histogram to project into
- TH2* localHist = NULL;
- TH1* projection = NULL;
switch (m_interpType)
{
case InterpType::Full:
@@ -308,27 +317,13 @@ StatusCode JetSmearingCorrection::readHisto(double& returnValue, TH1* histo, dou
break;
case InterpType::OnlyX:
- localHist = dynamic_cast<TH2*>(histo);
- if (!localHist)
- {
- ATH_MSG_ERROR("Could not convert the histogram to a TH2, please check inputs");
- return StatusCode::FAILURE;
- }
- projection = localHist->ProjectionX("projx_2D",histo->GetYaxis()->FindBin(y),histo->GetYaxis()->FindBin(y));
- returnValue = projection->Interpolate(x);
- delete projection;
+ // Determine the y-bin and use the cached projection to interpolate x
+ returnValue = projections.at(histo->GetYaxis()->FindBin(y))->Interpolate(x);
break;
case InterpType::OnlyY:
- localHist = dynamic_cast<TH2*>(histo);
- if (!localHist)
- {
- ATH_MSG_ERROR("Could not convert the histogram to a TH2, please check inputs");
- return StatusCode::FAILURE;
- }
- projection = localHist->ProjectionY("projy_2D",histo->GetXaxis()->FindBin(x),histo->GetXaxis()->FindBin(x));
- returnValue = projection->Interpolate(y);
- delete projection;
+ // Determine the x-bin and use the cached projection to interpolate y
+ returnValue = projections.at(histo->GetXaxis()->FindBin(x))->Interpolate(y);
break;
default:
@@ -339,88 +334,6 @@ StatusCode JetSmearingCorrection::readHisto(double& returnValue, TH1* histo, dou
return StatusCode::SUCCESS;
}
-StatusCode JetSmearingCorrection::readHisto(double& returnValue, TH1* histo, double x, double y, double z) const
-{
- // Ensure that the histogram exists
- if (!histo)
- {
- ATH_MSG_ERROR("Unable to read histogram - address is NULL");
- return StatusCode::FAILURE;
- }
-
- // Check dimensionality just to be safe
- if (histo->GetDimension() != 3)
- {
- ATH_MSG_ERROR("Blocking reading of a " << histo->GetDimension() << "D histogram as a 3D histogram");
- return StatusCode::FAILURE;
- }
-
- // Ensure we are within boundaries
- const double minX = histo->GetXaxis()->GetBinLowEdge(1);
- const double maxX = histo->GetXaxis()->GetBinLowEdge(histo->GetNbinsX()+1);
- if ( x >= maxX )
- x = maxX - 1.e-6;
- else if ( x <= minX )
- x = minX + 1.e-6;
- const double minY = histo->GetYaxis()->GetBinLowEdge(1);
- const double maxY = histo->GetYaxis()->GetBinLowEdge(histo->GetNbinsY()+1);
- if ( y >= maxY )
- y = maxY - 1.e-6;
- else if ( y <= minY )
- y = minY + 1.e-6;
- const double minZ = histo->GetZaxis()->GetBinLowEdge(1);
- const double maxZ = histo->GetZaxis()->GetBinLowEdge(histo->GetNbinsZ()+1);
- if ( z >= maxZ )
- z = maxZ - 1.e-6;
- else if ( z <= minZ )
- z = minZ + 1.e-6;
-
- // Get the result, interpolating as appropriate
- // We need to use this as a 3D histogram to have access to the projection methods
- // We also need a histogram to project into
- TH3* localHist = NULL;
- TH1* projection = NULL;
- switch (m_interpType)
- {
- case InterpType::Full:
- returnValue = histo->Interpolate(x,y,z);
- break;
-
- case InterpType::None:
- returnValue = histo->GetBinContent(histo->GetXaxis()->FindBin(x),histo->GetYaxis()->FindBin(y),histo->GetZaxis()->FindBin(z));
- break;
-
- case InterpType::OnlyX:
- localHist = dynamic_cast<TH3*>(histo);
- if (!localHist)
- {
- ATH_MSG_ERROR("Could not convert the histogram to a TH3, please check inputs");
- return StatusCode::FAILURE;
- }
- projection = localHist->ProjectionX("projx_3D",histo->GetYaxis()->FindBin(y),histo->GetYaxis()->FindBin(y),histo->GetZaxis()->FindBin(z),histo->GetZaxis()->FindBin(z));
- returnValue = projection->Interpolate(x);
- delete projection;
- break;
-
- case InterpType::OnlyY:
- localHist = dynamic_cast<TH3*>(histo);
- if (!localHist)
- {
- ATH_MSG_ERROR("Could not convert the histogram to a TH3, please check inputs");
- return StatusCode::FAILURE;
- }
- projection = localHist->ProjectionY("projy_3D",histo->GetXaxis()->FindBin(x),histo->GetXaxis()->FindBin(x),histo->GetZaxis()->FindBin(z),histo->GetZaxis()->FindBin(z));
- returnValue = projection->Interpolate(y);
- delete projection;
- break;
-
- default:
- ATH_MSG_ERROR("Unsupported interpolation type for a 3D histogram");
- return StatusCode::FAILURE;
- }
-
- return StatusCode::SUCCESS;
-}
StatusCode JetSmearingCorrection::getSigmaSmear(xAOD::Jet& jet, double& sigmaSmear) const
{
@@ -457,23 +370,23 @@ StatusCode JetSmearingCorrection::getSigmaSmear(xAOD::Jet& jet, double& sigmaSme
return StatusCode::SUCCESS;
}
-StatusCode JetSmearingCorrection::getNominalResolution(const xAOD::Jet& jet, TH1* histo, double& resolution) const
+StatusCode JetSmearingCorrection::getNominalResolution(const xAOD::Jet& jet, TH1* histo, const std::vector< std::unique_ptr<TH1> >& projections, double& resolution) const
{
double localRes = 0;
switch (m_histType)
{
case HistType::Pt:
- if (readHisto(localRes,histo,jet.pt()).isFailure())
+ if (readHisto(localRes,histo,jet.pt()/m_GeV).isFailure())
return StatusCode::FAILURE;
break;
case HistType::PtEta:
- if (readHisto(localRes,histo,jet.pt(),jet.eta()).isFailure())
+ if (readHisto(localRes,histo,projections,jet.pt()/m_GeV,jet.eta()).isFailure())
return StatusCode::FAILURE;
break;
case HistType::PtAbsEta:
- if (readHisto(localRes,histo,jet.pt(),fabs(jet.eta())).isFailure())
+ if (readHisto(localRes,histo,projections,jet.pt()/m_GeV,fabs(jet.eta())).isFailure())
return StatusCode::FAILURE;
break;
@@ -491,12 +404,12 @@ StatusCode JetSmearingCorrection::getNominalResolution(const xAOD::Jet& jet, TH1
StatusCode JetSmearingCorrection::getNominalResolutionData(const xAOD::Jet& jet, double& resolution) const
{
- return getNominalResolution(jet,m_smearResolutionData.get(),resolution);
+ return getNominalResolution(jet,m_smearResolutionData.get(),m_cachedProjResData,resolution);
}
StatusCode JetSmearingCorrection::getNominalResolutionMC(const xAOD::Jet& jet, double& resolution) const
{
- return getNominalResolution(jet,m_smearResolutionMC.get(),resolution);
+ return getNominalResolution(jet,m_smearResolutionMC.get(),m_cachedProjResMC,resolution);
}
@@ -546,4 +459,90 @@ StatusCode JetSmearingCorrection::calibrateImpl(xAOD::Jet& jet, JetEventInfo&) c
return StatusCode::SUCCESS;
}
+StatusCode JetSmearingCorrection::cacheProjections(TH1* fullHistogram, std::vector< std::unique_ptr<TH1> >& cacheLocation, const std::string& type)
+{
+ // Ensure the histogram exists
+ if (!fullHistogram)
+ {
+ ATH_MSG_FATAL("Cannot cache histogram as it doesn't exist: " << type);
+ return StatusCode::FAILURE;
+ }
+
+ // Ensure the number of dimensions is sane
+ if (fullHistogram->GetDimension() < 1 || fullHistogram->GetDimension() > 2)
+ {
+ ATH_MSG_FATAL("Unsupported histogram dimensionality for projection caching: " << fullHistogram->GetDimension());
+ return StatusCode::FAILURE;
+ }
+
+ // Protect vs InterpType
+ switch (m_interpType)
+ {
+ case InterpType::OnlyX:
+ // Simple case of 1D
+ if (fullHistogram->GetDimension() == 1)
+ return StatusCode::SUCCESS;
+ break;
+
+ case InterpType::OnlyY:
+ // Failure case of 1D
+ if (fullHistogram->GetDimension() == 1)
+ {
+ ATH_MSG_FATAL("Cannot project in Y for a 1D histogram: " << type);
+ return StatusCode::FAILURE;
+ }
+ break;
+
+ default:
+ ATH_MSG_FATAL("The interpolation type is not supported for caching: " << type);
+ return StatusCode::FAILURE;
+ }
+
+ // If we got here, then the request makes sense
+ // Start the projections
+ // Intentionally include underflow and overflow bins
+ // This keeps the same indexing scheme as root
+ // Avoids confusion and problems later at cost of a small amount of RAM
+ if (fullHistogram->GetDimension() == 2)
+ {
+ TH2* localHist = dynamic_cast<TH2*>(fullHistogram);
+ if (!localHist)
+ {
+ ATH_MSG_FATAL("Failed to convert histogram to a TH2, please check inputs: " << type);
+ return StatusCode::FAILURE;
+ }
+ if (m_interpType == InterpType::OnlyX)
+ {
+ for (Long64_t binY = 0; binY < localHist->GetNbinsY()+1; ++binY)
+ {
+ // Single bin of Y, interpolate across X
+ cacheLocation.push_back(std::unique_ptr<TH1>(localHist->ProjectionX(Form("projx_%s_%lld",type.c_str(),binY),binY,binY)));
+ }
+ }
+ else if (m_interpType == InterpType::OnlyY)
+ {
+ for (Long64_t binX = 0; binX < localHist->GetNbinsX()+1; ++binX)
+ {
+ // Single bin of X, interpolate across Y
+ cacheLocation.push_back(std::unique_ptr<TH1>(localHist->ProjectionY(Form("projy_%s_%lld",type.c_str(),binX),binX,binX)));
+ }
+ }
+ else
+ {
+ // We shouldn't make it here due to earlier checks
+ ATH_MSG_FATAL("Unexpected interpolation type, somehow escaped earlier checks: " << type);
+ return StatusCode::FAILURE;
+ }
+ }
+ else
+ {
+ // We shouldn't make it here due to earlier checks
+ ATH_MSG_FATAL("Unexpected dimensionality: " << fullHistogram->GetDimension());
+ return StatusCode::FAILURE;
+ }
+
+ // All done
+ return StatusCode::SUCCESS;
+}
+
diff --git a/Reconstruction/Jet/JetCalibTools/util/JetCalibTools_SmearingPlots.cxx b/Reconstruction/Jet/JetCalibTools/util/JetCalibTools_SmearingPlots.cxx
index 3c7cf7465ef7972d270b5fdac1c1063e1831df72..ace27098b2dfc86876397fa54f382694110faf5f 100644
--- a/Reconstruction/Jet/JetCalibTools/util/JetCalibTools_SmearingPlots.cxx
+++ b/Reconstruction/Jet/JetCalibTools/util/JetCalibTools_SmearingPlots.cxx
@@ -14,11 +14,13 @@
#include "AsgTools/AnaToolHandle.h"
#include <vector>
+#include <ctime>
#include "TString.h"
#include "TH1D.h"
#include "TF1.h"
#include "TCanvas.h"
#include "TLatex.h"
+#include "TRandom3.h"
namespace jet
@@ -41,9 +43,9 @@ namespace jet
int main (int argc, char* argv[])
{
// Check argument usage
- if (argc != 5 && argc != 6)
+ if (argc < 5 || argc > 7)
{
- std::cout << "USAGE: " << argv[0] << " <JetCollection> <ConfigFile> <OutputFile> <isData> (dev mode switch)" << std::endl;
+ std::cout << "USAGE: " << argv[0] << " <JetCollection> <ConfigFile> <OutputFile> <isData> (dev mode switch) (timing test switch)" << std::endl;
return 1;
}
@@ -53,6 +55,7 @@ int main (int argc, char* argv[])
const TString outFile = argv[3];
const bool isData = (TString(argv[4]) == "true");
const bool isDevMode = ( argc > 5 && (TString(argv[5]) == "true" || TString(argv[5]) == "dev") ) ? true : false;
+ const bool isTimeTest = ( argc > 6 && TString(argv[6]) == "true" ) ? true : false;
// Derived information
const bool outFileIsExtensible = outFile.EndsWith(".pdf") || outFile.EndsWith(".ps");
@@ -61,10 +64,11 @@ int main (int argc, char* argv[])
// Assumed constants
const TString calibSeq = "Smear"; // only want to apply the smearing correction here
const std::vector<float> ptVals = {20, 40, 60, 100, 400, 1000};
- const float eta = 0;
+ const float eta = 0.202;
const float phi = 0;
- const float mass = 100;
+ const float mass = 10;
const int maxNumIter = 1e5;
+ const int numTimeIter = 100;
// Accessor strings
const TString startingScaleString = "JetGSCScaleMomentum";
@@ -133,7 +137,46 @@ int main (int argc, char* argv[])
{
const TString baseName = Form("Smear_%.0f",pt);
hists_pt.push_back(new TH1D(baseName+"_pT",baseName+"_pT",100,0.25*pt,1.75*pt));
- hists_m.push_back(new TH1D(baseName+"_mass",baseName+"_pT",100,0.25*mass,1.75*mass));
+ hists_m.push_back(new TH1D(baseName+"_mass",baseName+"_mass",100,0.25*mass,1.75*mass));
+ }
+
+ // Run the timing test if specified
+ if (isTimeTest)
+ {
+ TRandom3 rand;
+ rand.SetSeed(0); // Deterministic random test
+
+ for (int numTest = 0; numTest < numTimeIter; ++numTest)
+ {
+ // Make a new calibration tool each time
+ JetCalibrationTool* jetCalibTool = new JetCalibrationTool(Form("mytool_%d",numTest));
+ if (jetCalibTool->setProperty("JetCollection",jetAlgo.Data()).isFailure())
+ exit(1);
+ if (jetCalibTool->setProperty("ConfigFile",config.Data()).isFailure())
+ exit(1);
+ if (jetCalibTool->setProperty("CalibSequence",calibSeq.Data()).isFailure())
+ exit(1);
+ if (jetCalibTool->setProperty("IsData",isData).isFailure())
+ exit(1);
+ if (isDevMode && jetCalibTool->setProperty("DEVmode",isDevMode).isFailure())
+ exit(1);
+ if (jetCalibTool->initialize().isFailure())
+ exit(1);
+
+ clock_t startTime = clock();
+ for (int numIter = 0; numIter < maxNumIter; ++numIter)
+ {
+ xAOD::JetFourMom_t fourvec(rand.Uniform(20.e3,1000.e3),rand.Uniform(-2,2),rand.Uniform(-3.14,3.14),10.e3);
+ jet->setJetP4(fourvec);
+ startingScale.setAttribute(*jet,fourvec);
+ xAOD::Jet* smearedJet = nullptr;
+ jetCalibTool->calibratedCopy(*jet,smearedJet);
+ //calibTool->calibratedCopy(*jet,smearedJet);
+ delete smearedJet;
+ }
+ delete jetCalibTool;
+ printf("Iteration %d: %f seconds\n",numTest+1,(clock()-startTime)/((double)CLOCKS_PER_SEC));
+ }
}
// Fill the histograms
@@ -147,12 +190,12 @@ int main (int argc, char* argv[])
hist_m->Sumw2();
- printf("Running for pT of %.0f GeV\n",pt);
+ printf("Running for pT of %.0f\n",pt);
for (int numIter = 0; numIter < maxNumIter; ++numIter)
{
// Set the jet four-vector
- xAOD::JetFourMom_t fourvec(pt,eta,phi,mass);
+ xAOD::JetFourMom_t fourvec(pt*1.e3,eta,phi,mass*1.e3);
jet->setJetP4(fourvec);
startingScale.setAttribute(*jet,fourvec);
@@ -173,8 +216,8 @@ int main (int argc, char* argv[])
}
// Fill the histograms
- hist_pt->Fill(smearedJet->pt());
- hist_m->Fill(smearedJet->m());
+ hist_pt->Fill(smearedJet->pt()/1.e3);
+ hist_m->Fill(smearedJet->m()/1.e3);
// Clean up
delete smearedJet;
@@ -190,7 +233,7 @@ int main (int argc, char* argv[])
{
// Set the jet four-vector
const float pt = nominalResData.GetXaxis()->GetBinCenter(binX);
- xAOD::JetFourMom_t fourvec(pt,eta,phi,mass);
+ xAOD::JetFourMom_t fourvec(pt*1.e3,eta,phi,mass*1.e3);
jet->setJetP4(fourvec);
// Jet kinematics set, now get the nominal resolutions