From 3a1b7f9444ee7cec4d1ad1f5e189f495536bb8ad Mon Sep 17 00:00:00 2001 From: mmlynari <104267016+mmlynari@users.noreply.github.com> Date: Fri, 17 May 2024 12:44:31 +0200 Subject: [PATCH] Update benchmark calibration and correction (#78) * 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 --- .../components/CalibrateBenchmarkMethod.cpp | 281 ++++---- .../src/components/CalibrateBenchmarkMethod.h | 31 +- .../src/components/CorrectCaloClusters.cpp | 259 ++++++- .../src/components/CorrectCaloClusters.h | 81 ++- RecFCCeeCalorimeter/CMakeLists.txt | 6 + .../tests/options/runCaloSim.py | 1 + .../tests/options/run_thetamodulemerged.py | 1 + .../run_thetamodulemerged_SW_benchmarkCorr.py | 643 ++++++++++++++++++ 8 files changed, 1150 insertions(+), 153 deletions(-) create mode 100644 RecFCCeeCalorimeter/tests/options/run_thetamodulemerged_SW_benchmarkCorr.py diff --git a/RecCalorimeter/src/components/CalibrateBenchmarkMethod.cpp b/RecCalorimeter/src/components/CalibrateBenchmarkMethod.cpp index 1347c028..52698d3e 100644 --- a/RecCalorimeter/src/components/CalibrateBenchmarkMethod.cpp +++ b/RecCalorimeter/src/components/CalibrateBenchmarkMethod.cpp @@ -19,6 +19,9 @@ #include "Math/Factory.h" #include "Math/Functor.h" +// Include the header for std::fabs +#include + DECLARE_COMPONENT(CalibrateBenchmarkMethod) @@ -69,74 +72,68 @@ StatusCode CalibrateBenchmarkMethod::initialize() { } } - // book histograms for checks - m_totalEnergyECal = new TH1F("totalEnergyECal", "Total deposited energy in ECal", 1000, 0, 1.5 * m_energy); - if (m_histSvc->regHist("/rec/ecal_total", m_totalEnergyECal).isFailure()) { - error() << "Couldn't register histogram" << endmsg; - return StatusCode::FAILURE; - } + // number of benchmark parameters + int n_param = 6; + // book histograms + m_totalEnergyECal = new TH1F("totalEnergyECal", "Total deposited energy in ECal", 1000, 0, 1.5 * m_energy); m_totalEnergyHCal = new TH1F("totalEnergyHCal", "Total deposited energy in HCal", 1000, 0, 1.5 * m_energy); - if (m_histSvc->regHist("/rec/hcal_total", m_totalEnergyHCal).isFailure()) { - error() << "Couldn't register histogram" << endmsg; - return StatusCode::FAILURE; - } - m_totalEnergyBoth = new TH1F("totalEnergyBoth", "Total deposited energy in ECal+HCal", 1000, 0, 1.5 * m_energy); - if (m_histSvc->regHist("/rec/both_total", m_totalEnergyBoth).isFailure()) { - error() << "Couldn't register histogram" << endmsg; - return StatusCode::FAILURE; - } + m_parameters = new TH1F("parameters","Benchmark parameters", n_param, 0, n_param); + + registerHistogram("/rec/ecal_total", m_totalEnergyECal); + registerHistogram("/rec/hcal_total", m_totalEnergyHCal); + registerHistogram("/rec/both_total", m_totalEnergyBoth); + registerHistogram("/rec/parameters", m_parameters); - m_parameters = new TH1F("parameters","", 4, 0,4); - if (m_histSvc->regHist("/rec/parameters", m_parameters).isFailure()) { - error() << "Couldn't register histogram" << endmsg; - return StatusCode::FAILURE; - } // clear vectors that will be later used for fitting - m_vecEgenerated.clear(); - m_vecEinECaltotal.clear(); - m_vecEinHCaltotal.clear(); - m_vecEinHCalfirstLayer.clear(); - m_vecEinECallastLayer.clear(); + m_vecGeneratedEnergy.clear(); + m_vecTotalEnergyinECal.clear(); + m_vecTotalEnergyinHCal.clear(); + m_vecEnergyInFirstLayerHCal.clear(); + m_vecEnergyInLastLayerECal.clear(); + m_vecEnergyInFirstLayerECal.clear(); - m_energyInECalLayer.resize(m_numLayersECal); - m_energyInHCalLayer.resize(m_numLayersHCal); + m_energyInLayerECal.resize(m_numLayersECal); + m_energyInLayerHCal.resize(m_numLayersHCal); return StatusCode::SUCCESS; } StatusCode CalibrateBenchmarkMethod::execute() { - double energyInFirstHCalLayer = 0; - double energyInLastECalLayer = 0; + double totalEnergyInECal = 0.; + double totalEnergyInHCal = 0.; + double energyInFirstLayerECal = 0; + double energyInLastLayerECal = 0; + double energyInFirstLayerHCal = 0; double energyInBoth = 0.; - int ecal_index = -1; - int hcal_index = -1; + int indexECal = -1; + int indexHCal = -1; - for (double& eneEcal : m_energyInECalLayer){ - eneEcal = 0.; + for (double& eneECal : m_energyInLayerECal){ + eneECal = 0.; } - for (double& eneHcal : m_energyInHCalLayer){ - eneHcal = 0.; + for (double& eneHCal : m_energyInLayerHCal){ + eneHCal = 0.; } // identify ECal and HCal readout position in the input parameters when running the calibration for (uint j=0; jgetDetector()->readout(m_readoutNames[ecal_index]).idSpec().decoder(); - auto decoder_HCal = m_geoSvc->getDetector()->readout(m_readoutNames[hcal_index]).idSpec().decoder(); + auto decoderECal = m_geoSvc->getDetector()->readout(m_readoutNames[indexECal]).idSpec().decoder(); + auto decoderHCal = m_geoSvc->getDetector()->readout(m_readoutNames[indexHCal]).idSpec().decoder(); std::vector cellIDs; @@ -150,72 +147,92 @@ StatusCode CalibrateBenchmarkMethod::execute() { // Loop over a collection of ECal cells and get energy in each ECal layer for (const auto& iCell : *ecalBarrelCells) { dd4hep::DDSegmentation::CellID cellID = iCell.getCellID(); - size_t layerIDecal = decoder_ECal->get(cellID, "layer"); - m_energyInECalLayer.at(layerIDecal) += iCell.getEnergy(); + size_t layerIDecal = decoderECal->get(cellID, "layer"); + m_energyInLayerECal.at(layerIDecal) += iCell.getEnergy(); } // Loop over a collection of HCal cells and get energy in each HCal layer for (const auto& iCell : *hcalBarrelCells) { dd4hep::DDSegmentation::CellID cellID = iCell.getCellID(); - size_t layerIDhcal = decoder_HCal->get(cellID, "layer"); - m_energyInHCalLayer.at(layerIDhcal) += iCell.getEnergy(); + size_t layerIDhcal = decoderHCal->get(cellID, "layer"); + m_energyInLayerHCal.at(layerIDhcal) += iCell.getEnergy(); } // Energy deposited in the whole ECal - const double energyInECal = std::accumulate(m_energyInECalLayer.begin(), m_energyInECalLayer.end(), 0); + for (size_t i = 0; i < m_energyInLayerECal.size(); ++i) { + totalEnergyInECal += m_energyInLayerECal.at(i); + } - // Energy deposited in the last ECal layer - energyInLastECalLayer = m_energyInECalLayer.at(m_numLayersECal-1); + // Energy deposited in the first/last ECal layer + energyInFirstLayerECal = m_energyInLayerECal.at(m_firstLayerECal); + energyInLastLayerECal = m_energyInLayerECal.at(m_numLayersECal-1); // Energy deposited in the whole HCal - const double energyInHCal = std::accumulate(m_energyInHCalLayer.begin(), m_energyInHCalLayer.end(), 0); + for (size_t i = 0; i < m_energyInLayerHCal.size(); ++i) { + totalEnergyInHCal += m_energyInLayerHCal.at(i); + } // Energy deposited in the first HCal layer - energyInFirstHCalLayer = m_energyInHCalLayer.at(m_firstLayerHCal); + energyInFirstLayerHCal = m_energyInLayerHCal.at(m_firstLayerHCal); // Total energy deposited in ECal and HCal - energyInBoth = energyInECal+energyInHCal; + energyInBoth = totalEnergyInECal+totalEnergyInHCal; // Fill histograms with ECal and HCal energy - m_totalEnergyECal->Fill(energyInECal); - m_totalEnergyHCal->Fill(energyInHCal); + m_totalEnergyECal->Fill(totalEnergyInECal); + m_totalEnergyHCal->Fill(totalEnergyInHCal); m_totalEnergyBoth->Fill(energyInBoth); // Fill vectors that will be later used for the energy fitting - length of the vector = number of events - m_vecEgenerated.push_back(m_energy); - m_vecEinECaltotal.push_back(energyInECal); - m_vecEinHCaltotal.push_back(energyInHCal); - m_vecEinHCalfirstLayer.push_back(energyInFirstHCalLayer); - m_vecEinECallastLayer.push_back(energyInLastECalLayer); + m_vecGeneratedEnergy.push_back(m_energy); + m_vecTotalEnergyinECal.push_back(totalEnergyInECal); + m_vecTotalEnergyinHCal.push_back(totalEnergyInHCal); + m_vecEnergyInFirstLayerECal.push_back(energyInFirstLayerECal); + m_vecEnergyInLastLayerECal.push_back(energyInLastLayerECal); + m_vecEnergyInFirstLayerHCal.push_back(energyInFirstLayerHCal); // Printouts for checks verbose() << "********************************************************************" << endmsg; - verbose() << "Energy in ECAL and HCAL: " << energyInECal << " GeV" << endmsg; - verbose() << "Energy in ECAL: " << energyInECal << " GeV" << endmsg; - verbose() << "Energy in HCAL: " << energyInHCal << " GeV" << endmsg; - verbose() << "Energy in ECAL last layer: " << energyInLastECalLayer << " GeV" << endmsg; - verbose() << "Energy in HCAL first layer: " << energyInFirstHCalLayer << " GeV" << endmsg; + verbose() << "Energy in ECAL and HCAL: " << energyInBoth << " GeV" << endmsg; + verbose() << "Energy in ECAL: " << totalEnergyInECal << " GeV" << endmsg; + verbose() << "Energy in HCAL: " << totalEnergyInHCal << " GeV" << endmsg; + verbose() << "Energy in ECAL last layer: " << energyInLastLayerECal << " GeV" << endmsg; + verbose() << "Energy in HCAL first layer: " << energyInFirstLayerHCal << " GeV" << endmsg; return StatusCode::SUCCESS; } // minimisation function for the benchmark method -double CalibrateBenchmarkMethod::chiSquareFitBarrel(const Double_t *par) const { +double CalibrateBenchmarkMethod::chiSquareFitBarrel(const Double_t *parameter) const { double fitvalue = 0.; // loop over all events, vector of a size of #evts is filled with the energies // ECal calibrated to EM scale, HCal calibrated to HAD scale - for(uint i=0; iSetMaxFunctionCalls(1000000); // for Minuit/Minuit2 - min->SetMaxIterations(10000); // for GSL - min->SetTolerance(0.001); - min->SetPrintLevel(1); - - // create funciton wrapper for minimizer - int n_param = 4; - ROOT::Math::Functor f(this, &CalibrateBenchmarkMethod::chiSquareFitBarrel,n_param); - - static const std::vector steps = {0.001, 0.001, 0.001, 0.001}; - static const std::vector variable = {1., 1., .5, .01}; - min->SetFunction(f); - - // Variables to be minimized - min->SetVariable(0,"p0", variable[0], steps[0] ); - min->SetVariable(1,"p1", variable[1], steps[1] ); - min->SetVariable(2,"p2", variable[2], steps[2] ); - min->SetVariable(3,"p3", variable[3], steps[3] ); - - // fix par1 because HCal is already calibrated to HAD scale - min->FixVariable(1); - // min->SetVariableLimits(2, 0, 1e6); - - // do the minimization - min->Minimize(); - - const Double_t *xs = min->X(); - const Double_t *ys = min->Errors(); - std::cout << "Minimum: f(" << xs[0] << "," << xs[1] << "," << xs[2] << "," << xs[3] << "): " << min->MinValue() << std::endl; - std::cout << "Minimum: #delta f(" << ys[0] << "," << ys[1] << "," << ys[2] << "," << ys[3] <SetBinContent(1, xs[0]); - m_parameters->SetBinContent(2, xs[1]); - m_parameters->SetBinContent(3, xs[2]); - m_parameters->SetBinContent(4, xs[3]); + // number of benchmark parameters + int n_param = 6; + // initial values for the minimization and steps while running the minimization + static const std::vector variable = {1., 1., .5, .01, 0.5, 0.}; + static const std::vector steps = {0.001, 0.001, 0.001, 0.001, 0.001, 0.001}; - m_parameters->SetBinError(1, ys[0]); - m_parameters->SetBinError(2, ys[1]); - m_parameters->SetBinError(3, ys[2]); - m_parameters->SetBinError(4, ys[3]); + runMinimization(n_param, variable, steps, m_fixedParameters); return GaudiAlgorithm::finalize(); } + +void CalibrateBenchmarkMethod::registerHistogram(const std::string& path, TH1F*& histogramName) { + if (m_histSvc->regHist(path, histogramName).isFailure()) { + error() << "Couldn't register histogram" << endmsg; + throw std::runtime_error("Histogram registration failure"); + } +} + +void CalibrateBenchmarkMethod::runMinimization(int n_param, const std::vector& variable, const std::vector& steps, const std::vector& fixedParameters) const { + // Set up the functor (the function to be minimized) + ROOT::Math::Functor f(this, &CalibrateBenchmarkMethod::chiSquareFitBarrel, n_param); + + // Initialize the minimizer + ROOT::Math::Minimizer* minimizer = ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad"); + minimizer->SetMaxFunctionCalls(1000000); + minimizer->SetMaxIterations(10000); + minimizer->SetTolerance(0.001); + minimizer->SetPrintLevel(1); + + // Set the function to minimize + minimizer->SetFunction(f); + + // Set parameters for the minimization + for (int i = 0; i < n_param; ++i) { + minimizer->SetVariable(i, "p" + std::to_string(i), variable[i], steps[i]); + + // Fix parameters listed in fixedParameters that will not be included in the fit minimization + if (std::find(fixedParameters.begin(), fixedParameters.end(), i) != fixedParameters.end()) { + minimizer->FixVariable(i); + } + } + + // Perform the minimization + minimizer->Minimize(); + + // Access the results if needed + const Double_t* xs = minimizer->X(); + const Double_t* ys = minimizer->Errors(); + + std::cout << "Minimum: f("; + for (int i = 0; i < n_param; ++i) { + std::cout << xs[i]; + if (i < n_param - 1) std::cout << ","; + } + std::cout << "): " << minimizer->MinValue() << std::endl; + + std::cout << "Minimum: #delta f("; + for (int i = 0; i < n_param; ++i) { + std::cout << ys[i]; + if (i < n_param - 1) std::cout << ","; + } + std::cout << std::endl; + + //histogram->SetBinContent(1, xs[0]); + //histogram->SetBinError(1, ys[0]); + + // each fitted parameter is stored in one bin of the output histogram + // bin 1 contains the value for par[0] etc + for (int i = 0; i < n_param; ++i) { + verbose() << "inside filling histos: " << xs[i] << " GeV" << endmsg; + m_parameters->SetBinContent(i+1, xs[i]); + m_parameters->SetBinError(i+1, ys[i]); + } +} + diff --git a/RecCalorimeter/src/components/CalibrateBenchmarkMethod.h b/RecCalorimeter/src/components/CalibrateBenchmarkMethod.h index 8d6083a1..8512efec 100644 --- a/RecCalorimeter/src/components/CalibrateBenchmarkMethod.h +++ b/RecCalorimeter/src/components/CalibrateBenchmarkMethod.h @@ -64,6 +64,8 @@ class CalibrateBenchmarkMethod : public GaudiAlgorithm { ServiceHandle 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& variable, const std::vector& steps, const std::vector& fixedParameters) const; /// Handle for electromagnetic barrel cells (input collection) DataHandle m_ecalBarrelCells{"ecalBarrelCells", Gaudi::DataHandle::Reader, this}; @@ -79,8 +81,8 @@ class CalibrateBenchmarkMethod : public GaudiAlgorithm { TH1F* m_parameters; /// vectors to store the energy in each ECal/HCal layer - std::vector m_energyInECalLayer; - std::vector m_energyInHCalLayer; + std::vector m_energyInLayerECal; + std::vector m_energyInLayerHCal; /// Maximum energy for the x-axis range @@ -90,21 +92,28 @@ class CalibrateBenchmarkMethod : public GaudiAlgorithm { Gaudi::Property m_numLayersECal{this, "numLayersECal", 12, "Number of ECal layers"}; Gaudi::Property m_numLayersHCal{this, "numLayersHCal", 13, "Number of HCal layers"}; - /// ID of the first HCal layer + /// ID of the first ECal/HCal layer + Gaudi::Property m_firstLayerECal{this, "firstLayerECal", 0, "ID of first ECal layer"}; Gaudi::Property m_firstLayerHCal{this, "firstLayerHCal", 0, "ID of first HCal layer"}; /// Names of the detector readouts, corresponding to system IDs Gaudi::Property> m_readoutNames {this, "readoutNames", {"ECalBarrelReadout","HCalBarrelReadout"},"Names of the detector readout, corresponding to systemId"}; - Gaudi::Property> m_SystemID {this, "systemID", {4,8},"systemId"}; - Gaudi::Property m_ECalSystemID{this, "ECalSystemID", 4, "ID of ECal system"}; - Gaudi::Property m_HCalSystemID{this, "HCalSystemID", 8, "ID of the HCal system"}; + Gaudi::Property> m_systemID {this, "systemID", {4,8},"systemId"}; + Gaudi::Property m_systemIDECal{this, "ECalSystemID", 4, "ID of ECal system"}; + Gaudi::Property m_systemIDHCal{this, "HCalSystemID", 8, "ID of the HCal system"}; /// vectors containing the energy deposits to be used for minimization - std::vector m_vecEgenerated; - std::vector m_vecEinECaltotal; - std::vector m_vecEinHCaltotal; - std::vector m_vecEinHCalfirstLayer; - std::vector m_vecEinECallastLayer; + std::vector m_vecGeneratedEnergy; + std::vector m_vecTotalEnergyinECal; + std::vector m_vecTotalEnergyinHCal; + std::vector m_vecEnergyInFirstLayerECal; + std::vector m_vecEnergyInLastLayerECal; + std::vector 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> m_fixedParameters = {this, "fixedParameters", {1,5},"Fixed parameters that will not be minimized"}; }; #endif /* RECCALORIMETER_CALIBRATEBENCHMARKMETHOD_H */ \ No newline at end of file diff --git a/RecCalorimeter/src/components/CorrectCaloClusters.cpp b/RecCalorimeter/src/components/CorrectCaloClusters.cpp index 939c7d99..5552e47c 100644 --- a/RecCalorimeter/src/components/CorrectCaloClusters.cpp +++ b/RecCalorimeter/src/components/CorrectCaloClusters.cpp @@ -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" @@ -24,6 +26,9 @@ // ROOT #include "TF2.h" +// Include the header for std::fabs +#include + DECLARE_COMPONENT(CorrectCaloClusters) CorrectCaloClusters::CorrectCaloClusters(const std::string& name, @@ -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 @@ -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) { @@ -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; } @@ -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; @@ -233,7 +339,6 @@ StatusCode CorrectCaloClusters::initializeCorrFunctions(std::vectorSetParameter(k, paramVec.at(j)); std::string parName(1, 'a' + (char) j); func->SetParName(k, parName.c_str()); - j += 1; } } @@ -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; isize(); ++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 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, @@ -334,7 +543,6 @@ double CorrectCaloClusters::getEnergyInLayer(edm4hep::Cluster cluster, if (decoder->get(cellID, "layer") != layerID) { continue; } - energy += cell->getEnergy(); } @@ -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; +} diff --git a/RecCalorimeter/src/components/CorrectCaloClusters.h b/RecCalorimeter/src/components/CorrectCaloClusters.h index 6696ce71..3d5168f1 100644 --- a/RecCalorimeter/src/components/CorrectCaloClusters.h +++ b/RecCalorimeter/src/components/CorrectCaloClusters.h @@ -30,6 +30,8 @@ namespace edm4hep { namespace dd4hep { namespace DDSegmentation { class FCCSWGridPhiEta_k4geo; + class FCCSWGridModuleThetaMerged_k4geo; + class FCCSWGridPhiTheta_k4geo; class MultiSegmentation; class BitFieldCoder; } @@ -44,10 +46,17 @@ namespace dd4hep { * * Downstream energy correction corrects for the energy lost due to punch trough. This energy is deposited in the * instrumentation behind the active volume of the calorimeter. It is parametrized in one (cluster energy) or two * variables (cluster energy, cluster angle) + * * Benchmark calibration should be used for the combined simulation of ECal and HCal when using charged pions. + * At the input level, the ECal should be calibrated to EM scale and HCal should be calibrated to HAD scale. + * The aim of the benchmark calibration is to bring ECal to HAD scale and also to take into account + * the energy loss between the ECal and HCal (e.g. in cryostat) - for this, the energy from the last ECal layer and the first HCal layer is used. + * While, downstream correction is part of the benchmark method (energy lost between ECal and HCal), as well as the upstream correction for ECal + * For standalone ECal, include upstream and downstream correction; for ECal+HCal simulation apply only benchmark correction should be applied + * To obtain the actual parameters run RecCalorimeter/tests/options/fcc_ee_caloBenchmarkCalibration.py which calls CalibrateBenchmarkMethod * * Based on similar corrections by Jana Faltova and Anna Zaborowska. * - * @author Juraj Smiesko + * @author Juraj Smiesko, benchmark calibration added by Michaela Mlynarikova */ class CorrectCaloClusters : public GaudiAlgorithm { @@ -79,7 +88,7 @@ class CorrectCaloClusters : public GaudiAlgorithm { StatusCode initializeCorrFunctions(std::vector>& functions, std::vector> formulas, std::vector> parameters, - const std::string& funcNameStem = "upDown"); + const std::string& funcNameStem = "upDownBenchmark"); /** * Apply upstream correction to the output clusters. @@ -103,12 +112,24 @@ class CorrectCaloClusters : public GaudiAlgorithm { StatusCode applyDownstreamCorr(const edm4hep::ClusterCollection* inClusters, edm4hep::ClusterCollection* outClusters); + /** + * Apply benchmark correction to the output clusters. + * + * @param[in] inClusters Pointer to the input cluster collection. + * @param[out] outClusters Pointer to the output cluster collection. + * + * @return Status code. + */ + StatusCode applyBenchmarkCorr(const edm4hep::ClusterCollection* inClusters, + edm4hep::ClusterCollection* outClusters); + /** * Get sum of energy from cells in specified layer. * This energy is not calibrated. * * @param[in] cluster Pointer to cluster of interest. * @param[in] readoutName Name of the readout. + * @param[in] systemID ID of the system. * @param[in] layerID ID of the layer of the readout. * * @return Energy in layer. @@ -118,6 +139,20 @@ class CorrectCaloClusters : public GaudiAlgorithm { int systemID, int layerID); + /** + * Get sum of energy from cells in the whole calorimeter. + * This energy is not calibrated. + * + * @param[in] cluster Pointer to cluster of interest. + * @param[in] readoutName Name of the readout. + * @param[in] systemID ID of the system. + * + * @return Total energy. + */ + double getTotalEnergy(edm4hep::Cluster cluster, + const std::string& readoutName, + int systemID); + /** * Get the theta angle of the specified cluster. * @@ -142,26 +177,30 @@ class CorrectCaloClusters : public GaudiAlgorithm { std::vector> m_upstreamFunctions; /// Pointers to downstream correction functions std::vector> m_downstreamFunctions; + /// Pointers to benchmark method correction functions + std::vector> m_benchmarkFunctions; /// IDs of the detectors Gaudi::Property> m_systemIDs { - this, "systemIDs", {4}, "IDs of systems" + this, "systemIDs", {4,8}, "IDs of systems" }; + Gaudi::Property m_systemIDECal{this, "systemIDECal", 4, "ID of ECal system"}; + Gaudi::Property m_systemIDHCal{this, "systemIDHCal", 8, "ID of the HCal system"}; /// Names of the detector readouts, corresponding to system IDs Gaudi::Property> m_readoutNames { - this, "readoutNames", {"ECalBarrelPhiEta"}, + this, "readoutNames", {"ECalBarrelModuleThetaMerged","BarHCal_Readout_phitheta"}, "Names of the detector readout, corresponding to systemID" }; /// Numbers of layers of the detectors Gaudi::Property> m_numLayers { - this, "numLayers", {12}, "Numbers of layers of the systems"}; + this, "numLayers", {12,13}, "Numbers of layers of the systems"}; /// IDs of the first layers of the detectors Gaudi::Property> m_firstLayerIDs { - this, "firstLayerIDs", {0}, "IDs of first layers in the systems" + this, "firstLayerIDs", {0,0}, "IDs of first layers in the systems" }; /// IDs of the last layers of the detectors Gaudi::Property> m_lastLayerIDs { - this, "lastLayerIDs", {7}, "IDs of last layers in the systems" + this, "lastLayerIDs", {11,12}, "IDs of last layers in the systems" }; /// Upstream correction parameters (a, b, c, ...) @@ -177,6 +216,32 @@ class CorrectCaloClusters : public GaudiAlgorithm { /// Downstream correction formulas in ROOT notation Gaudi::Property>> m_downstreamFormulas { this, "downstreamFormulas", {}, "Downstream formulas (in ROOT notation)"}; +/// Benchmark method correction parameters (a, b, c, ...) + Gaudi::Property>> m_benchmarkParametrization { + this, "benchmarkParametrization", {}, "Benchmark formula parameters"}; + /// Benchmark correction formulas in ROOT notation + Gaudi::Property>> m_benchmarkFormulas { + this, "benchmarkFormulas", {}, "Benchmark formulas (in ROOT notation)"}; + + /// Possibility to use different parametrization for very low energies (<10GeV) for benchmark method + /// For the energy (GeV) below this given threshold, the second energy formula in benchmarkFormulas is used + Gaudi::Property m_benchmarkEneSwitch { + this, "benchmarkEneSwitch", -1., + "Energy threshold in GeV to use the second formula. Set to a negative value to disable using the second formula." +}; + + /// benchmarkParamsApprox are used for the first estimate of the energy + /// which is then used to obtain energy dependent benchmark parameters + /// Parameters below were obtained by running a single benchmark calibration for 100 GeV pion + Gaudi::Property> m_benchmarkParamsApprox { + this, "benchmarkParamsApprox", {1.2281, 1, 1.07731, -0.00191622, 18.6772, 0.}, "Approximate benchmark parameters" + }; + /// Flag if upstream correction should be applied + Gaudi::Property m_upstreamCorr{this, "upstreamCorr", true}; + /// Flag if downstream correction should be applied + Gaudi::Property m_downstreamCorr{this, "downstreamCorr", true}; + /// Flag if benchmark correction should be applied + Gaudi::Property m_benchmarkCorr{this, "benchmarkCorr", false}; }; -#endif /* RECCALORIMETER_CORRECTCALOCLUSTERS_H */ +#endif /* RECCALORIMETER_CORRECTCALOCLUSTERS_H */ \ No newline at end of file diff --git a/RecFCCeeCalorimeter/CMakeLists.txt b/RecFCCeeCalorimeter/CMakeLists.txt index 31dd6889..c6d0beeb 100644 --- a/RecFCCeeCalorimeter/CMakeLists.txt +++ b/RecFCCeeCalorimeter/CMakeLists.txt @@ -74,3 +74,9 @@ add_test(NAME FCCeeLAr_xtalkNeighbours WORKING_DIRECTORY ${PROJECT_SOURCE_DIR} ) set_test_env(FCCeeLAr_xtalkNeighbours) + +add_test(NAME FCCeeLAr_benchmarkCorrection + COMMAND k4run RecFCCeeCalorimeter/tests/options/run_thetamodulemerged_SW_benchmarkCorr.py + WORKING_DIRECTORY ${PROJECT_SOURCE_DIR} +) +set_test_env(FCCeeLAr_benchmarkCorrection) diff --git a/RecFCCeeCalorimeter/tests/options/runCaloSim.py b/RecFCCeeCalorimeter/tests/options/runCaloSim.py index 22d585d2..8234999d 100644 --- a/RecFCCeeCalorimeter/tests/options/runCaloSim.py +++ b/RecFCCeeCalorimeter/tests/options/runCaloSim.py @@ -281,6 +281,7 @@ correctCaloClusters = CorrectCaloClusters("correctCaloClusters") correctCaloClusters.inClusters = createClusters.clusters.Path correctCaloClusters.outClusters = "Corrected"+createClusters.clusters.Path +correctCaloClusters.systemIDs = [4] correctCaloClusters.numLayers = [12] correctCaloClusters.firstLayerIDs = [0] correctCaloClusters.lastLayerIDs = [11] diff --git a/RecFCCeeCalorimeter/tests/options/run_thetamodulemerged.py b/RecFCCeeCalorimeter/tests/options/run_thetamodulemerged.py index b31539ab..812b61f8 100644 --- a/RecFCCeeCalorimeter/tests/options/run_thetamodulemerged.py +++ b/RecFCCeeCalorimeter/tests/options/run_thetamodulemerged.py @@ -520,6 +520,7 @@ correctCaloClusters = CorrectCaloClusters("correctCaloClusters", inClusters=createClusters.clusters.Path, outClusters="Corrected" + createClusters.clusters.Path, + systemIDs=[4], numLayers=[11], firstLayerIDs=[0], lastLayerIDs=[10], diff --git a/RecFCCeeCalorimeter/tests/options/run_thetamodulemerged_SW_benchmarkCorr.py b/RecFCCeeCalorimeter/tests/options/run_thetamodulemerged_SW_benchmarkCorr.py new file mode 100644 index 00000000..04f08712 --- /dev/null +++ b/RecFCCeeCalorimeter/tests/options/run_thetamodulemerged_SW_benchmarkCorr.py @@ -0,0 +1,643 @@ +from Configurables import ApplicationMgr +from Configurables import EventCounter +from Configurables import AuditorSvc, ChronoAuditor +from Configurables import PodioOutput +from Configurables import CaloTowerToolFCCee +from Configurables import CreateCaloClustersSlidingWindowFCCee +from Configurables import CorrectCaloClusters +from Configurables import CalibrateCaloClusters +from Configurables import CreateEmptyCaloCellsCollection +from Configurables import CreateCaloCellPositionsFCCee +from Configurables import CellPositionsECalBarrelModuleThetaSegTool +from Configurables import RedoSegmentation +from Configurables import CreateCaloCells +from Configurables import CalibrateCaloHitsTool +from Configurables import CalibrateInLayersTool +from Configurables import SimG4Alg +from Configurables import SimG4PrimariesFromEdmTool +from Configurables import SimG4SaveCalHits +from Configurables import SimG4ConstantMagneticFieldTool +from Configurables import SimG4Svc +from Configurables import SimG4FullSimActions +from Configurables import SimG4SaveParticleHistory +from Configurables import GeoSvc +from Configurables import HepMCToEDMConverter +from Configurables import GenAlg +from Configurables import FCCDataSvc +from Gaudi.Configuration import INFO +# from Gaudi.Configuration import * + +import os + +from GaudiKernel.SystemOfUnits import GeV, tesla, mm +from GaudiKernel.PhysicalConstants import pi, halfpi, twopi +from math import cos, sin, tan + +use_pythia = False +addNoise = False +dumpGDML = False +runHCal = True +# for big productions, save significant space removing hits and cells +# however, hits and cluster cells might be wanted for small productions for detailed event displays +# also, cluster cells are needed for the MVA training +saveHits = False +saveCells = False +saveClusterCells = False + +# cluster energy corrections +# simple parametrisations of up/downstream losses +applyUpDownstreamBenchmarkCorrections = True + +# Input for simulations (momentum is expected in GeV!) +# Parameters for the particle gun simulations, dummy if use_pythia is set +# to True +# theta from 80 to 100 degrees corresponds to -0.17 < eta < 0.17 +# reminder: cell granularity in theta = 0.5625 degrees +# (in strips: 0.5625/4=0.14) + +# Nevts = 20000 +Nevts = 10 +# Nevts = 1 +# Nevts=1000 + +# particle momentum and direction +# momentum = 100 # in GeV +momentum = 100 # in GeV +# momentum = 10 # in GeV +thetaMin = 69. # degrees +thetaMax = 69. # degrees +# thetaMin = 89 +# thetaMax = 91 +# thetaMin = 90 # degrees +# thetaMax = 90 # degrees +# phiMin = halfpi +# phiMax = halfpi +phiMin = 0 +phiMax = twopi + +# particle origin +# origR = 1000.0*mm +origR = 0.0 * mm +origTheta = halfpi +origPhi = 0.0 + +# particle type: 11 electron, 13 muon, 22 photon, 111 pi0, 211 pi+ +pdgCode = 211 +# pdgCode = 22 +# pdgCode = 111 +# pdgCode = 211 + +# Set to true if history from Geant4 decays is needed (e.g. to get the +# photons from pi0) +saveG4Hist = False +if (pdgCode == 111): + saveG4Hist = True + +magneticField = False + + +podioevent = FCCDataSvc("EventDataSvc") + +# Particle gun setup + +genAlg = GenAlg() +if use_pythia: + from Configurables import PythiaInterface + pythia8gentool = PythiaInterface() + pythia8gentool.pythiacard = os.path.join( + os.environ.get('PWD', ''), + "MCGeneration/ee_Zgamma_inclusive.cmd" + ) + # pythia8gentool.pythiacard = "MCGeneration/ee_Z_ee.cmd" + pythia8gentool.printPythiaStatistics = False + pythia8gentool.pythiaExtraSettings = [""] + genAlg.SignalProvider = pythia8gentool +else: + from Configurables import MomentumRangeParticleGun + pgun = MomentumRangeParticleGun("ParticleGun") + pgun.PdgCodes = [pdgCode] + pgun.MomentumMin = momentum * GeV + pgun.MomentumMax = momentum * GeV + pgun.PhiMin = phiMin + pgun.PhiMax = phiMax + pgun.ThetaMin = thetaMin * pi / 180. + pgun.ThetaMax = thetaMax * pi / 180. + genAlg.SignalProvider = pgun + +genAlg.hepmc.Path = "hepmc" + +# smear/shift vertex +if origR > 0.0: + origX = origR * cos(origPhi) + origY = origR * sin(origPhi) + origZ = origR / tan(origTheta) + print("Particle gun will be moved to %f %f %f" % (origX, origY, origZ)) + from Configurables import GaussSmearVertex + vertexSmearAndShiftTool = GaussSmearVertex() + vertexSmearAndShiftTool.xVertexSigma = 0. + vertexSmearAndShiftTool.yVertexSigma = 0. + vertexSmearAndShiftTool.tVertexSigma = 0. + vertexSmearAndShiftTool.xVertexMean = origX + vertexSmearAndShiftTool.yVertexMean = origY + vertexSmearAndShiftTool.zVertexMean = origZ + vertexSmearAndShiftTool.tVertexMean = 0. + genAlg.VertexSmearingTool = vertexSmearAndShiftTool + +# hepMC -> EDM converter +hepmc_converter = HepMCToEDMConverter() +hepmc_converter.hepmc.Path = "hepmc" +genParticlesOutputName = "genParticles" +hepmc_converter.GenParticles.Path = genParticlesOutputName +hepmc_converter.hepmcStatusList = [] +hepmc_converter.OutputLevel = INFO + +# Simulation setup +# Detector geometry +geoservice = GeoSvc("GeoSvc") +# if K4GEO is empty, this should use relative path to working directory +path_to_detector = os.environ.get("K4GEO", "") +print(path_to_detector) +detectors_to_use = [ + 'FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/ALLEGRO_o1_v02.xml' +] +# prefix all xmls with path_to_detector +geoservice.detectors = [ + os.path.join(path_to_detector, _det) for _det in detectors_to_use +] +geoservice.OutputLevel = INFO + +# Geant4 service +# Configures the Geant simulation: geometry, physics list and user actions +actions = SimG4FullSimActions() + +if saveG4Hist: + actions.enableHistory = True + actions.energyCut = 1.0 * GeV + saveHistTool = SimG4SaveParticleHistory("saveHistory") + +geantservice = SimG4Svc( + "SimG4Svc", + detector='SimG4DD4hepDetector', + physicslist="SimG4FtfpBert", + actions=actions +) + +# from Configurables import GeoToGdmlDumpSvc +if dumpGDML: + from Configurables import GeoToGdmlDumpSvc + gdmldumpservice = GeoToGdmlDumpSvc("GeoToGdmlDumpSvc") + +# Fixed seed to have reproducible results, change it for each job if you +# split one production into several jobs +# Mind that if you leave Gaudi handle random seed and some job start within +# the same second (very likely) you will have duplicates +geantservice.randomNumbersFromGaudi = False +geantservice.seedValue = 4242 + +# Range cut +geantservice.g4PreInitCommands += ["/run/setCut 0.1 mm"] + +# Magnetic field +if magneticField == 1: + field = SimG4ConstantMagneticFieldTool( + "SimG4ConstantMagneticFieldTool", + FieldComponentZ=-2 * tesla, + FieldOn=True, + IntegratorStepper="ClassicalRK4" + ) +else: + field = SimG4ConstantMagneticFieldTool( + "SimG4ConstantMagneticFieldTool", + FieldOn=False + ) + +# Geant4 algorithm +# Translates EDM to G4Event, passes the event to G4, writes out outputs +# via tools and a tool that saves the calorimeter hits + +# Detector readouts +# ECAL +ecalBarrelReadoutName = "ECalBarrelModuleThetaMerged" +ecalBarrelReadoutName2 = "ECalBarrelModuleThetaMerged2" +ecalEndcapReadoutName = "ECalEndcapPhiEta" +# HCAL +if runHCal: + hcalBarrelReadoutName = "HCalBarrelReadout" + hcalBarrelReadoutName2 = "BarHCal_Readout_phitheta" + hcalEndcapReadoutName = "HCalEndcapReadout" +else: + hcalBarrelReadoutName = "" + hcalBarrelReadoutName2 = "" + hcalEndcapReadoutName = "" + +# Configure saving of calorimeter hits +ecalBarrelHitsName = "ECalBarrelPositionedHits" +saveECalBarrelTool = SimG4SaveCalHits( + "saveECalBarrelHits", + readoutName=ecalBarrelReadoutName, + OutputLevel=INFO +) +saveECalBarrelTool.CaloHits.Path = ecalBarrelHitsName + +saveECalEndcapTool = SimG4SaveCalHits( + "saveECalEndcapHits", + readoutName=ecalEndcapReadoutName +) +saveECalEndcapTool.CaloHits.Path = "ECalEndcapHits" + +if runHCal: + hcalBarrelHitsName = "HCalBarrelPositionedHits" + saveHCalTool = SimG4SaveCalHits( + "saveHCalBarrelHits", + readoutName=hcalBarrelReadoutName + ) + saveHCalTool.CaloHits.Path = hcalBarrelHitsName + + # saveHCalEndcapTool = SimG4SaveCalHits( + # "saveHCalEndcapHits", + # readoutName = hcalEndcapReadoutName + # ) + # saveHCalEndcapTool.CaloHits.Path = "HCalEndcapHits" + +# next, create the G4 algorithm, giving the list of names of tools ("XX/YY") +particle_converter = SimG4PrimariesFromEdmTool("EdmConverter") +particle_converter.GenParticles.Path = genParticlesOutputName + +outputTools = [ + saveECalBarrelTool, + saveECalEndcapTool +] +if runHCal: + outputTools += [ + saveHCalTool, + # saveHCalEndcapTool + ] + +if saveG4Hist: + outputTools += [saveHistTool] + +geantsim = SimG4Alg("SimG4Alg", + outputs=outputTools, + eventProvider=particle_converter, + OutputLevel=INFO) + +# Digitization (Merging hits into cells, EM scale calibration) +# EM scale calibration (sampling fraction) +calibEcalBarrel = CalibrateInLayersTool("CalibrateECalBarrel", + samplingFraction=[0.3864252122990472] * 1 + [0.13597644835735828] * 1 + [0.14520427829645913] * 1 + [0.1510076084632846] * 1 + [0.1552347580991012] * 1 + [0.159694330729184] * 1 + [0.1632954482794191] * 1 + [0.16720711037339814] * 1 + [0.17047749048884808] * 1 + [0.17461698117974286] * 1 + [0.1798984163980135] * 1 + [0.17920355117405806] * 1, + readoutName=ecalBarrelReadoutName, + layerFieldName="layer") + +calibEcalEndcap = CalibrateCaloHitsTool( + "CalibrateECalEndcap", invSamplingFraction="4.27") +if runHCal: + calibHcells = CalibrateCaloHitsTool( + "CalibrateHCal", invSamplingFraction="30.4") + calibHcalEndcap = CalibrateCaloHitsTool( + "CalibrateHCalEndcap", invSamplingFraction="31.7") + +# Create cells in ECal barrel +# 1. step - merge hits into cells with theta and module segmentation +# (module is a 'physical' cell i.e. lead + LAr + PCB + LAr +lead) +# 2. step - rewrite the cellId using the merged theta-module segmentation +# (merging several modules and severla theta readout cells). +# Add noise at this step if you derived the noise already assuming merged cells + +# Step 1: merge hits into cells according to initial segmentation +ecalBarrelCellsName = "ECalBarrelCells" +createEcalBarrelCells = CreateCaloCells("CreateECalBarrelCells", + doCellCalibration=True, + calibTool=calibEcalBarrel, + addCellNoise=False, + filterCellNoise=False, + addPosition=True, + OutputLevel=INFO, + hits=ecalBarrelHitsName, + cells=ecalBarrelCellsName) + +# Step 2a: compute new cellID of cells based on new readout +# (merged module-theta segmentation with variable merging vs layer) +resegmentEcalBarrel = RedoSegmentation("ReSegmentationEcal", + # old bitfield (readout) + oldReadoutName=ecalBarrelReadoutName, + # specify which fields are going to be altered (deleted/rewritten) + oldSegmentationIds=["module", "theta"], + # new bitfield (readout), with new segmentation (merged modules and theta cells) + newReadoutName=ecalBarrelReadoutName2, + OutputLevel=INFO, + debugPrint=200, + inhits=ecalBarrelCellsName, + outhits="ECalBarrelCellsMerged") + +# Step 2b: merge new cells with same cellID together +# do not apply cell calibration again since cells were already +# calibrated in Step 1 +ecalBarrelCellsName2 = "ECalBarrelCells2" +createEcalBarrelCells2 = CreateCaloCells("CreateECalBarrelCells2", + doCellCalibration=False, + addCellNoise=False, + filterCellNoise=False, + OutputLevel=INFO, + hits="ECalBarrelCellsMerged", + cells=ecalBarrelCellsName2) + +# Add to Ecal barrel cells the position information +# (good for physics, all coordinates set properly) + +cellPositionEcalBarrelTool = CellPositionsECalBarrelModuleThetaSegTool( + "CellPositionsECalBarrel", + readoutName=ecalBarrelReadoutName, + OutputLevel=INFO +) +ecalBarrelPositionedCellsName = "ECalBarrelPositionedCells" +createEcalBarrelPositionedCells = CreateCaloCellPositionsFCCee( + "CreateECalBarrelPositionedCells", + OutputLevel=INFO +) +createEcalBarrelPositionedCells.positionsTool = cellPositionEcalBarrelTool +createEcalBarrelPositionedCells.hits.Path = ecalBarrelCellsName +createEcalBarrelPositionedCells.positionedHits.Path = ecalBarrelPositionedCellsName + +cellPositionEcalBarrelTool2 = CellPositionsECalBarrelModuleThetaSegTool( + "CellPositionsECalBarrel2", + readoutName=ecalBarrelReadoutName2, + OutputLevel=INFO +) +createEcalBarrelPositionedCells2 = CreateCaloCellPositionsFCCee( + "CreateECalBarrelPositionedCells2", + OutputLevel=INFO +) +createEcalBarrelPositionedCells2.positionsTool = cellPositionEcalBarrelTool2 +createEcalBarrelPositionedCells2.hits.Path = ecalBarrelCellsName2 +createEcalBarrelPositionedCells2.positionedHits.Path = "ECalBarrelPositionedCells2" + + +# Create cells in ECal endcap +createEcalEndcapCells = CreateCaloCells("CreateEcalEndcapCaloCells", + doCellCalibration=True, + calibTool=calibEcalEndcap, + addCellNoise=False, + filterCellNoise=False, + OutputLevel=INFO) +createEcalEndcapCells.hits.Path = "ECalEndcapHits" +createEcalEndcapCells.cells.Path = "ECalEndcapCells" + +if runHCal: + # Create cells in HCal + # 1 - merge hits into cells with the default readout + hcalBarrelCellsName = "HCalBarrelCells" + createHcalBarrelCells = CreateCaloCells("CreateHCalBarrelCells", + doCellCalibration=True, + calibTool=calibHcells, + addCellNoise=False, + filterCellNoise=False, + addPosition=True, + hits=hcalBarrelHitsName, + cells=hcalBarrelCellsName, + OutputLevel=INFO) + + # 2 - attach positions to the cells (cell positions needed for RedoSegmentation!) + from Configurables import CellPositionsHCalBarrelPhiThetaSegTool + cellPositionHcalBarrelTool = CellPositionsHCalBarrelPhiThetaSegTool( + "CellPositionsHCalBarrel", + readoutName=hcalBarrelReadoutName, + OutputLevel=INFO + ) + hcalBarrelPositionedCellsName = "HCalBarrelPositionedCells" + createHcalBarrelPositionedCells = CreateCaloCellPositionsFCCee( + "CreateHcalBarrelPositionedCells", + OutputLevel=INFO + ) + createHcalBarrelPositionedCells.positionsTool = cellPositionHcalBarrelTool + createHcalBarrelPositionedCells.hits.Path = hcalBarrelCellsName + createHcalBarrelPositionedCells.positionedHits.Path = hcalBarrelPositionedCellsName + + # 3 - compute new cellID of cells based on new readout - removing row information + # We use a RedoSegmentation. Using a RewriteBitField with removeIds=["row"], + # wont work because there are tiles with same layer/theta/phi but different row + # as a consequence there will be multiple cells with same cellID in the output collection + # and this will screw up the SW clustering + hcalBarrelCellsName2 = "HCalBarrelCells2" + + # first we create new hits with the readout without the row information + # and then merge them into new cells + rewriteHCalBarrel = RedoSegmentation("ReSegmentationHcal", + # old bitfield (readout) + oldReadoutName=hcalBarrelReadoutName, + # specify which fields are going to be altered (deleted/rewritten) + oldSegmentationIds=["row", "theta", "phi"], + # new bitfield (readout), with new segmentation (merged modules and theta cells) + newReadoutName=hcalBarrelReadoutName2, + OutputLevel=INFO, + debugPrint=200, + inhits=hcalBarrelPositionedCellsName, + outhits="HCalBarrelCellsWithoutRow") + + createHcalBarrelCells2 = CreateCaloCells("CreateHCalBarrelCells2", + doCellCalibration=False, + addCellNoise=False, + filterCellNoise=False, + OutputLevel=INFO, + hits=rewriteHCalBarrel.outhits.Path, + cells=hcalBarrelCellsName2) + + # 4 - attach positions to the new cells + from Configurables import CellPositionsHCalBarrelPhiThetaSegTool + hcalBarrelPositionedCellsName2 = "HCalBarrelPositionedCells2" + cellPositionHcalBarrelTool2 = CellPositionsHCalBarrelPhiThetaSegTool( + "CellPositionsHCalBarrel2", + readoutName=hcalBarrelReadoutName2, + OutputLevel=INFO + ) + createHcalBarrelPositionedCells2 = CreateCaloCellPositionsFCCee( + "CreateHCalBarrelPositionedCells2", + OutputLevel=INFO + ) + createHcalBarrelPositionedCells2.positionsTool = cellPositionHcalBarrelTool2 + createHcalBarrelPositionedCells2.hits.Path = hcalBarrelCellsName2 + createHcalBarrelPositionedCells2.positionedHits.Path = hcalBarrelPositionedCellsName2 + + + # createHcalEndcapCells = CreateCaloCells("CreateHcalEndcapCaloCells", + # doCellCalibration=True, + # calibTool=calibHcalEndcap, + # addCellNoise=False, + # filterCellNoise=False, + # OutputLevel=INFO) + # createHcalEndcapCells.hits.Path="HCalEndcapHits" + # createHcalEndcapCells.cells.Path="HCalEndcapCells" + +else: + hcalBarrelCellsName = "emptyCaloCells" + hcalBarrelPositionedCellsName = "emptyCaloCells" + hcalBarrelCellsName2 = "emptyCaloCells" + hcalBarrelPositionedCellsName2 = "emptyCaloCells" + cellPositionHcalBarrelTool = None + cellPositionHcalBarrelTool2 = None + +# Empty cells for parts of calorimeter not implemented yet +createemptycells = CreateEmptyCaloCellsCollection("CreateEmptyCaloCells") +createemptycells.cells.Path = "emptyCaloCells" + +# Produce sliding window clusters (ECAL+HCAL) +towers = CaloTowerToolFCCee("towers", + #deltaThetaTower=4 * 0.009817477 / 4, deltaPhiTower=2 * 2 * pi / 1536., + deltaThetaTower= 0.022180, + deltaPhiTower=2 * pi / 256., + max_layer=25, + ecalBarrelReadoutName=ecalBarrelReadoutName, + ecalEndcapReadoutName="", + ecalFwdReadoutName="", + hcalBarrelReadoutName=hcalBarrelReadoutName2, + hcalExtBarrelReadoutName="", + hcalEndcapReadoutName="", + hcalFwdReadoutName="", + OutputLevel=INFO) +towers.ecalBarrelCells.Path = ecalBarrelPositionedCellsName +towers.ecalEndcapCells.Path = "emptyCaloCells" +towers.ecalFwdCells.Path = "emptyCaloCells" + +towers.hcalBarrelCells.Path = hcalBarrelPositionedCellsName2 +towers.hcalExtBarrelCells.Path = "emptyCaloCells" +towers.hcalEndcapCells.Path = "emptyCaloCells" +towers.hcalFwdCells.Path = "emptyCaloCells" + +# Cluster variables (not optimized) +windT = 18 +windP = 34 +posT = 10 +posP = 22 +dupT = 14 +dupP = 26 +finT = 18 +finP = 34 +# Minimal energy to create a cluster in GeV (FCC-ee detectors have to reconstruct low energy particles) +threshold = 0.5 + +createClusters = CreateCaloClustersSlidingWindowFCCee("CreateClusters", + towerTool=towers, + nThetaWindow=windT, nPhiWindow=windP, + nThetaPosition=posT, nPhiPosition=posP, + nThetaDuplicates=dupT, nPhiDuplicates=dupP, + nThetaFinal=finT, nPhiFinal=finP, + energyThreshold=threshold, + energySharingCorrection=False, + attachCells=True, + OutputLevel=INFO + ) +createClusters.clusters.Path = "CaloClusters" +createClusters.clusterCells.Path = "CaloClusterCells" + +correctCaloClusters = CorrectCaloClusters("correctCaloClusters", + inClusters=createClusters.clusters.Path, + outClusters="Corrected" + createClusters.clusters.Path, + systemIDs=[4,8], + numLayers=[12,13], + firstLayerIDs=[0,0], + lastLayerIDs=[11,12], + readoutNames=[ecalBarrelReadoutName,hcalBarrelReadoutName2], + # do not split the following line or it will break scripts that update the values of the corrections + upstreamParameters = [[0.03900891447361534, -4.322941016402328, -139.1811369546787, 0.498342628339746, -3.3545078429754813, -13.99996971344221],[]], + upstreamFormulas=[['[0]+[1]/(x-[2])', '[0]+[1]/(x-[2])'],[]], + # do not split the following line or it will break scripts that update the values of the corrections + downstreamParameters = [[-0.0027661744480442195, 0.006059143775380306, 0.9788596364251927, -1.4951749409378743, -0.08491999337012696, 16.017621428757778],[]], + downstreamFormulas=[['[0]+[1]*x', '[0]+[1]/sqrt(x)', '[0]+[1]/x'],[]], + ## Approximate parameters for the first energy estimate obtained from the benchmark calib for 100 GeV pion + benchmarkParamsApprox = [1.22, 1, 1.04, -0.0019, 25.69, 0.], + ## below parameters and formulas for two functions with ene switch at 9~GeV, the constant term is set to zero + #benchmarkParametrization = [[1.0109, 2.4768, 1., 1.0837, -5.278, -1.9938, 0.0006, -0.2715, 45.29, -5674.31, 126.04, 0.],[2.0432, 0.2902, -0.0709, 0.004, 1., 1.2099, -432.2, 13.73, -0.0023, -0.2572, 1.99, -0.968, 0.259, -0.0152, 0.]], + #benchmarkFormulas = [['[0]+[1]/sqrt(x)', '[0]', '[0]+[1]/(x+[2])', '[0]+[1]/x', '[0]+[1]/(x+[2])', '[0]'],['[0]+[1]*x+[2]*x*x+[3]*x*x*x', '[0]', '[0]+[1]/(x+[2])**2', '[0]+[1]/x', '[0]+[1]*x+[2]*x*x+[3]*x*x*x', '[0]']], + ## below parameters and formulas for one function (redireved in Nov 2023) + benchmarkParametrization = [[2.04, 1.0626, 1., -3.5, 1.0182, -0.26, 0.0004, 45.9, -5906.11, -129.27, 0.],[]], + benchmarkFormulas = [['[0]/sqrt(x)+[1]', '[0]', '[0]/x+[1]', '[0]/x+[1]', '[0]+[1]/(x-[2])', '[0]'],[]], + benchmarkEneSwitch = -1., + upstreamCorr = False, + downstreamCorr = False, + benchmarkCorr = True, + OutputLevel=INFO + ) + +# Output +out = PodioOutput("out", + OutputLevel=INFO) + +if runHCal: + out.outputCommands = ["keep *", "drop emptyCaloCells"] +else: + out.outputCommands = ["keep *", "drop HCal*", "drop emptyCaloCells"] + +if not saveCells: + out.outputCommands.append("drop ECal*Cells*") +if not saveClusterCells: + out.outputCommands.append("drop *ClusterCells*") +if not saveHits: + out.outputCommands.append("drop ECal*Hits*") + out.outputCommands.append("drop HCal*Hits*") + +out.filename = "./output_evts_" + str(Nevts) + "_pdg_" + str(pdgCode) + "_pMin_" + str(momentum) + "_GeV" + "_ThetaMinMax_" + str(thetaMin) + "_" + str( + thetaMax) + "_PhiMinMax_" + str(phiMin) + "_" + str(phiMax) + "_MagneticField_" + str(magneticField) + "_Noise" + str(addNoise) + ".root" + +# CPU information +chra = ChronoAuditor() +audsvc = AuditorSvc() +audsvc.Auditors = [chra] +genAlg.AuditExecute = True +hepmc_converter.AuditExecute = True +geantsim.AuditExecute = True +createEcalBarrelCells.AuditExecute = True +createEcalBarrelPositionedCells.AuditExecute = True +if runHCal: + createHcalBarrelCells.AuditExecute = True +out.AuditExecute = True + +event_counter = EventCounter('event_counter') +event_counter.Frequency = 10 + +ExtSvc = [geoservice, podioevent, geantservice, audsvc] +if dumpGDML: + ExtSvc += [gdmldumpservice] + +TopAlg = [ + event_counter, + genAlg, + hepmc_converter, + geantsim, + createEcalBarrelCells, + createEcalBarrelPositionedCells, + resegmentEcalBarrel, + createEcalBarrelCells2, + createEcalBarrelPositionedCells2, + createEcalEndcapCells +] + +if runHCal: + TopAlg += [ + createHcalBarrelCells, + createHcalBarrelPositionedCells, + rewriteHCalBarrel, + createHcalBarrelCells2, + createHcalBarrelPositionedCells2, + # createHcalEndcapCells + ] + +TopAlg += [ + createemptycells, + createClusters, +] + +if applyUpDownstreamBenchmarkCorrections: + TopAlg += [ + correctCaloClusters, + ] + +TopAlg += [ + out +] + +ApplicationMgr( + TopAlg=TopAlg, + EvtSel='NONE', + EvtMax=Nevts, + ExtSvc=ExtSvc, + StopOnSignal=True, +)