Skip to content

Commit

Permalink
Update benchmark calibration and correction (#78)
Browse files Browse the repository at this point in the history 
* update benchmark calibration

* update tests with benchmark correction

* Implement changes suggested by Juraj

* update tests

* add more safety checks and small improvements

* small modifications of safety checks
  • Loading branch information
@mmlynari
mmlynari authored May 17, 2024
1 parent 5e88ce2 commit 3a1b7f9
Show file tree
Hide file tree
Showing 8 changed files with 1,150 additions and 153 deletions.
281 changes: 164 additions & 117 deletions RecCalorimeter/src/components/CalibrateBenchmarkMethod.cpp
View file

Large diffs are not rendered by default.

31 changes: 20 additions & 11 deletions RecCalorimeter/src/components/CalibrateBenchmarkMethod.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,6 +64,8 @@ class CalibrateBenchmarkMethod : public GaudiAlgorithm {
ServiceHandle<ITHistSvc> m_histSvc;

double chiSquareFitBarrel(const Double_t *par) const;
void registerHistogram(const std::string& path, TH1F*& histogramName);
void runMinimization(int n_param, const std::vector<double>& variable, const std::vector<double>& steps, const std::vector<int>& fixedParameters) const;

/// Handle for electromagnetic barrel cells (input collection)
DataHandle<edm4hep::CalorimeterHitCollection> m_ecalBarrelCells{"ecalBarrelCells", Gaudi::DataHandle::Reader, this};
Expand All @@ -79,8 +81,8 @@ class CalibrateBenchmarkMethod : public GaudiAlgorithm {
TH1F* m_parameters;

/// vectors to store the energy in each ECal/HCal layer
std::vector<double> m_energyInECalLayer;
std::vector<double> m_energyInHCalLayer;
std::vector<double> m_energyInLayerECal;
std::vector<double> m_energyInLayerHCal;


/// Maximum energy for the x-axis range
Expand All @@ -90,21 +92,28 @@ class CalibrateBenchmarkMethod : public GaudiAlgorithm {
Gaudi::Property<size_t> m_numLayersECal{this, "numLayersECal", 12, "Number of ECal layers"};
Gaudi::Property<size_t> m_numLayersHCal{this, "numLayersHCal", 13, "Number of HCal layers"};

/// ID of the first HCal layer
/// ID of the first ECal/HCal layer
Gaudi::Property<uint> m_firstLayerECal{this, "firstLayerECal", 0, "ID of first ECal layer"};
Gaudi::Property<uint> m_firstLayerHCal{this, "firstLayerHCal", 0, "ID of first HCal layer"};

/// Names of the detector readouts, corresponding to system IDs
Gaudi::Property<std::vector<std::string>> m_readoutNames {this, "readoutNames", {"ECalBarrelReadout","HCalBarrelReadout"},"Names of the detector readout, corresponding to systemId"};
Gaudi::Property<std::vector<uint>> m_SystemID {this, "systemID", {4,8},"systemId"};
Gaudi::Property<uint> m_ECalSystemID{this, "ECalSystemID", 4, "ID of ECal system"};
Gaudi::Property<uint> m_HCalSystemID{this, "HCalSystemID", 8, "ID of the HCal system"};
Gaudi::Property<std::vector<uint>> m_systemID {this, "systemID", {4,8},"systemId"};
Gaudi::Property<uint> m_systemIDECal{this, "ECalSystemID", 4, "ID of ECal system"};
Gaudi::Property<uint> m_systemIDHCal{this, "HCalSystemID", 8, "ID of the HCal system"};

/// vectors containing the energy deposits to be used for minimization
std::vector<double> m_vecEgenerated;
std::vector<double> m_vecEinECaltotal;
std::vector<double> m_vecEinHCaltotal;
std::vector<double> m_vecEinHCalfirstLayer;
std::vector<double> m_vecEinECallastLayer;
std::vector<double> m_vecGeneratedEnergy;
std::vector<double> m_vecTotalEnergyinECal;
std::vector<double> m_vecTotalEnergyinHCal;
std::vector<double> m_vecEnergyInFirstLayerECal;
std::vector<double> m_vecEnergyInLastLayerECal;
std::vector<double> m_vecEnergyInFirstLayerHCal;

// benchmark parameters which should be fixed
// p[1] because HCal is already calibrated to HAD scale
// p[5] is the constant term and was not bringing large improvement, hence not minimized for the moment (might be tried again in the future)
Gaudi::Property<std::vector<int>> m_fixedParameters = {this, "fixedParameters", {1,5},"Fixed parameters that will not be minimized"};

};
#endif /* RECCALORIMETER_CALIBRATEBENCHMARKMETHOD_H */
259 changes: 242 additions & 17 deletions RecCalorimeter/src/components/CorrectCaloClusters.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,8 @@
// FCC Detectors
#include "detectorCommon/DetUtils_k4geo.h"
#include "detectorSegmentations/FCCSWGridPhiEta_k4geo.h"
#include "detectorSegmentations/FCCSWGridPhiTheta_k4geo.h"
#include "detectorSegmentations/FCCSWGridModuleThetaMerged_k4geo.h"

// DD4hep
#include "DD4hep/Detector.h"
Expand All @@ -24,6 +26,9 @@
// ROOT
#include "TF2.h"

// Include the <cmath> header for std::fabs
#include <cmath>

DECLARE_COMPONENT(CorrectCaloClusters)

CorrectCaloClusters::CorrectCaloClusters(const std::string& name,
Expand Down Expand Up @@ -72,21 +77,82 @@ StatusCode CorrectCaloClusters::initialize() {
error() << "Sizes of systemIDs vector and firstLayerIDs vector does not match, exiting!" << endmsg;
return StatusCode::FAILURE;
}
if (m_systemIDs.size() != m_upstreamFormulas.size()) {
error() << "Sizes of systemIDs vector and upstreamFormulas vector does not match, exiting!" << endmsg;
return StatusCode::FAILURE;
}
if (m_systemIDs.size() != m_upstreamParams.size()) {
error() << "Sizes of systemIDs vector and upstreamParams vector does not match, exiting!" << endmsg;
return StatusCode::FAILURE;
if (m_upstreamCorr) {
if (m_upstreamFormulas.empty()){
error() << "Upstream correction is turned on, but upstreamFormulas vector is empty, exiting!" << endmsg;
return StatusCode::FAILURE;
}
if (m_upstreamParams.empty()){
error() << "Upstream correction is turned on, but upstreamParams vector is empty, exiting!" << endmsg;
return StatusCode::FAILURE;
}
if (m_systemIDs.size() != m_upstreamFormulas.size()) {
error() << "Sizes of systemIDs vector and upstreamFormulas vector does not match, exiting!" << endmsg;
return StatusCode::FAILURE;
}
if (m_systemIDs.size() != m_upstreamParams.size()) {
error() << "Sizes of systemIDs vector and upstreamParams vector does not match, exiting!" << endmsg;
return StatusCode::FAILURE;
}
}
if (m_systemIDs.size() != m_downstreamFormulas.size()) {
error() << "Sizes of systemIDs vector and downstreamFormulas vector does not match, exiting!" << endmsg;
return StatusCode::FAILURE;
if (m_downstreamCorr) {
if (m_downstreamFormulas.empty()){
error() << "Downstream correction is turned on, but downstreamFormulas vector is empty, exiting!" << endmsg;
return StatusCode::FAILURE;
}
if (m_downstreamParams.empty()){
error() << "Downstream correction is turned on, but downstreamParams vector is empty, exiting!" << endmsg;
return StatusCode::FAILURE;
}
if (m_systemIDs.size() != m_downstreamFormulas.size()) {
error() << "Sizes of systemIDs vector and downstreamFormulas vector does not match, exiting!" << endmsg;
return StatusCode::FAILURE;
}
if (m_systemIDs.size() != m_downstreamParams.size()) {
error() << "Sizes of systemIDs vector and downstreamParams vector does not match, exiting!" << endmsg;
return StatusCode::FAILURE;
}
}
if (m_systemIDs.size() != m_downstreamParams.size()) {
error() << "Sizes of systemIDs vector and downstreamParams vector does not match, exiting!" << endmsg;
return StatusCode::FAILURE;
if (m_benchmarkCorr) {
// number of benchmark parameters is 6
size_t nBenchmarkParams = 6;
// check that all benchmark parameters and functions are not empty and have the correct size
if (m_benchmarkFormulas.empty()){
error() << "Benchmark correction is turned on, but benchmarkFormulas vector is empty, exiting!" << endmsg;
return StatusCode::FAILURE;
}
if (m_benchmarkParametrization.empty()){
error() << "Benchmark correction is turned on, but benchmarkParametrization vector is empty, exiting!" << endmsg;
return StatusCode::FAILURE;
}
if (m_benchmarkFormulas.size() != m_benchmarkParametrization.size()){
error() << "Sizes of benchmarkFormulas vector and benchmarkParametrization vector does not match, exiting! " << m_benchmarkFormulas.size() << endmsg;
return StatusCode::FAILURE;
}
if (m_benchmarkParamsApprox.size() != nBenchmarkParams){
error() << "The benchmarkParamsApprox vector does not have the proper size, the size should be " << nBenchmarkParams << endmsg;
return StatusCode::FAILURE;
}
if (m_benchmarkEneSwitch<0.){
if (m_benchmarkFormulas[0].size() != nBenchmarkParams){
error() << "A single formula parametrization for benchmark correction is turned on, but the benchmarkFormulas[0] vector does not have the proper size, the size should be " << nBenchmarkParams << endmsg;
return StatusCode::FAILURE;
}
}
if (m_benchmarkEneSwitch>0.){
if (m_benchmarkFormulas.size()<2){
error() << "The two formulas parametrization for benchmark correction is turned on, but the benchmarkFormulas vector does not have the proper size, the size should be 2." << endmsg;
return StatusCode::FAILURE;
}
if (m_benchmarkFormulas[0].size() != nBenchmarkParams){
error() << "The two formulas parametrization for benchmark correction is turned on, but the benchmarkFormulas[0] vector does not have the proper size, the size should be " << nBenchmarkParams << endmsg;
return StatusCode::FAILURE;
}
if (m_benchmarkFormulas[1].size() != nBenchmarkParams){
error() << "The two formulas parametrization for benchmark correction is turned on, but the benchmarkFormulas[1] vector does not have the proper size, the size should be " << nBenchmarkParams << endmsg;
return StatusCode::FAILURE;
}
}
}

// Prepare upstream and downstream correction functions
Expand All @@ -104,6 +170,13 @@ StatusCode CorrectCaloClusters::initialize() {
return sc;
}
}
{
StatusCode sc = initializeCorrFunctions(m_benchmarkFunctions, m_benchmarkFormulas, m_benchmarkParametrization, "benchmark");
if (sc.isFailure()) {
error() << "Initialization of benchmark correction functions not successful!" << endmsg;
return sc;
}
}

info() << "Initialized following upstream correction functions:" << endmsg;
for (size_t i = 0; i < m_upstreamFunctions.size(); ++i) {
Expand All @@ -127,6 +200,17 @@ StatusCode CorrectCaloClusters::initialize() {
}
}

info() << "Initialized following benchmark correction functions:" << endmsg;
for (size_t i = 0; i < m_benchmarkFunctions.size(); ++i) {
for (size_t j = 0; j < m_benchmarkFunctions[i].size(); ++j) {
auto func = m_benchmarkFunctions.at(i).at(j);
info() << " " << func->GetName() << ": " << func->GetExpFormula() << endmsg;
for (int k = 0; k < func->GetNpar(); ++k) {
info() << " " << func->GetParName(k) << ": " << func->GetParameter(k) << endmsg;
}
}
}

return StatusCode::SUCCESS;
}

Expand All @@ -149,19 +233,41 @@ StatusCode CorrectCaloClusters::execute() {
}

// Apply upstream correction
{
if (m_upstreamCorr) {
StatusCode sc = applyUpstreamCorr(inClusters, outClusters);
if (sc.isFailure()) {
return sc;
}
else{
verbose() << "Running the upstream correction." << endmsg;

}
}

// Apply downstream correction
{
if (m_downstreamCorr) {
StatusCode sc = applyDownstreamCorr(inClusters, outClusters);
if (sc.isFailure()) {
return sc;
}
else{
verbose() << "Running the downstream correction." << endmsg;
}
}

// Apply benchmark correction
if (m_benchmarkCorr) {
StatusCode sc = applyBenchmarkCorr(inClusters, outClusters);
if (sc.isFailure()) {
return sc;
}
else{
verbose() << "Running the benchmark correction." << endmsg;
}
}

if ((m_upstreamCorr && m_downstreamCorr && m_benchmarkCorr) || (m_upstreamCorr && m_benchmarkCorr) || (m_downstreamCorr && m_benchmarkCorr)){
warning() << "Too many corrections in the house, apply upstream and downstream on ECal standalone and benchmark on combined ECal and HCal simulation." << endmsg;
}

return StatusCode::SUCCESS;
Expand Down Expand Up @@ -233,7 +339,6 @@ StatusCode CorrectCaloClusters::initializeCorrFunctions(std::vector<std::vector<
func->SetParameter(k, paramVec.at(j));
std::string parName(1, 'a' + (char) j);
func->SetParName(k, parName.c_str());

j += 1;
}
}
Expand Down Expand Up @@ -318,6 +423,110 @@ StatusCode CorrectCaloClusters::applyDownstreamCorr(const edm4hep::ClusterCollec
return StatusCode::SUCCESS;
}

StatusCode CorrectCaloClusters::applyBenchmarkCorr(const edm4hep::ClusterCollection* inClusters,
edm4hep::ClusterCollection* outClusters) {

const size_t numReadoutNames = m_readoutNames.size();

int ecal_index = -1;
int hcal_index = -1;

// Identify ECal and HCal readout positions in the input parameters when running the calibration
for (size_t i=0; i<numReadoutNames; ++i){
if (m_systemIDs[i] == m_systemIDECal){
ecal_index = i;
}
else if (m_systemIDs[i] == m_systemIDHCal){
hcal_index = i;
}
}

for (size_t j = 0; j < inClusters->size(); ++j) {
double energyInLastLayerECal = getEnergyInLayer(inClusters->at(j),
m_readoutNames[ecal_index],
m_systemIDs[ecal_index],
m_lastLayerIDs[ecal_index]);

double energyInFirstLayerECal = getEnergyInLayer(inClusters->at(j),
m_readoutNames[ecal_index],
m_systemIDs[ecal_index],
m_firstLayerIDs[ecal_index]);

double energyInFirstLayerHCal = getEnergyInLayer(inClusters->at(j),
m_readoutNames[hcal_index],
m_systemIDs[hcal_index],
m_firstLayerIDs[hcal_index]);

double totalEnergyInECal = getTotalEnergy(inClusters->at(j),
m_readoutNames[ecal_index],
m_systemIDs[ecal_index]);

double totalEnergyInHCal = getTotalEnergy(inClusters->at(j),
m_readoutNames[hcal_index],
m_systemIDs[hcal_index]);


// calculate approximate benchmark energy using non energy dependent benchmark parameters
double approximateBenchmarkEnergy = m_benchmarkParamsApprox[0] * totalEnergyInECal +
m_benchmarkParamsApprox[1] * totalEnergyInHCal +
m_benchmarkParamsApprox[2] * sqrt(abs(energyInLastLayerECal * m_benchmarkParamsApprox[0] * energyInFirstLayerHCal * m_benchmarkParamsApprox[1])) +
m_benchmarkParamsApprox[3] * pow(totalEnergyInECal * m_benchmarkParamsApprox[0], 2) +
m_benchmarkParamsApprox[4] * energyInFirstLayerECal +
m_benchmarkParamsApprox[5];

// Calculate energy-dependent benchmark parameters p[0]-p[5]
auto benchmarkFormulasHighEne = m_benchmarkFunctions.at(0);
int nParam = benchmarkFormulasHighEne.size();
std::vector<double> benchmarkParameters(nParam,-1.);

if (m_benchmarkEneSwitch>0.)
{
auto benchmarkFormulasLowEne = m_benchmarkFunctions.at(1);
// ensure smooth transition between low- and high-energy formulas (parameter l controls how fast the transition is)
int l=2;
double transition = 0.5 * (1 + tanh(l * (approximateBenchmarkEnergy - m_benchmarkEneSwitch) / 2));
verbose() << "Using two formulas for benchmark calibration, the second formula provided will be used to correct energies below benchmarkEneSwitch threshold." << endmsg;
for (size_t k = 0; k < benchmarkFormulasHighEne.size(); ++k) {
auto func_low_ene = benchmarkFormulasLowEne.at(k);
auto func_high_ene = benchmarkFormulasHighEne.at(k);
benchmarkParameters[k] = (1 - transition) * func_low_ene->Eval(approximateBenchmarkEnergy) +
transition * func_high_ene->Eval(approximateBenchmarkEnergy);
}
}
else{
verbose() << "Using one formula for benchmark calibration." << endmsg;
for (size_t k = 0; k < benchmarkFormulasHighEne.size(); ++k) {
auto func = benchmarkFormulasHighEne.at(k);
benchmarkParameters[k] = func->Eval(approximateBenchmarkEnergy);
}
}

// Get final benchmark energy using the energy dependent benchmark parameters
double benchmarkEnergy = benchmarkParameters[0] * totalEnergyInECal +
benchmarkParameters[1] * totalEnergyInHCal +
benchmarkParameters[2] * std::sqrt(std::fabs(energyInLastLayerECal * benchmarkParameters[0] * energyInFirstLayerHCal * benchmarkParameters[1])) +
benchmarkParameters[3] * std::pow(totalEnergyInECal * benchmarkParameters[0], 2) +
benchmarkParameters[4] * energyInFirstLayerECal +
benchmarkParameters[5];

// Protection against negative energy (might be improved)
if (benchmarkEnergy < 0.0) {
outClusters->at(j).setEnergy(totalEnergyInECal * benchmarkParameters[0] + totalEnergyInHCal * benchmarkParameters[1]);
} else {
outClusters->at(j).setEnergy(benchmarkEnergy);
}

if (inClusters->at(j).getEnergy()>1.){
debug() << "********************************************************************" << endmsg;
debug() << "Cluster energy: " << inClusters->at(j).getEnergy() << endmsg;
debug() << "totalEnergyInECal+HCal from hits: " << totalEnergyInECal+totalEnergyInHCal << endmsg;
debug() << "********************************************************************" << endmsg;
debug() << "Corrected cluster energy benchmark: " << outClusters->at(j).getEnergy() << endmsg;
debug() << "********************************************************************" << endmsg;
}
}
return StatusCode::SUCCESS;
}

double CorrectCaloClusters::getEnergyInLayer(edm4hep::Cluster cluster,
const std::string& readoutName,
Expand All @@ -334,7 +543,6 @@ double CorrectCaloClusters::getEnergyInLayer(edm4hep::Cluster cluster,
if (decoder->get(cellID, "layer") != layerID) {
continue;
}

energy += cell->getEnergy();
}

Expand All @@ -349,3 +557,20 @@ double CorrectCaloClusters::getClusterTheta(edm4hep::Cluster cluster) {

return theta;
}


double CorrectCaloClusters::getTotalEnergy(edm4hep::Cluster cluster,
const std::string& readoutName,
int systemID) {
dd4hep::DDSegmentation::BitFieldCoder* decoder = m_geoSvc->getDetector()->readout(readoutName).idSpec().decoder();

double energy = 0;
for (auto cell = cluster.hits_begin(); cell != cluster.hits_end(); ++cell) {
dd4hep::DDSegmentation::CellID cellID = cell->getCellID();
if (decoder->get(cellID, "system") != systemID) {
continue;
}
energy += cell->getEnergy();
}
return energy;
}
Loading

0 comments on commit 3a1b7f9

Please sign in to comment.