diff --git a/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/PixelConfigCondAlg.cxx b/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/PixelConfigCondAlg.cxx index 7024b858a1a6bd0a29b07e0465a0dfb83280c57a..2ba97d9cfcb5ad4a85ae66cddd3cfab82a24626d 100644 --- a/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/PixelConfigCondAlg.cxx +++ b/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/PixelConfigCondAlg.cxx @@ -1,5 +1,5 @@ /* - Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration */ #include "PixelConfigCondAlg.h" @@ -458,12 +458,11 @@ StatusCode PixelConfigCondAlg::execute(const EventContext& ctx) const { } } - std::vector<TH3F*> ramoPotentialMap; - std::vector<TH2F*> lorentzMap_e; - std::vector<TH2F*> lorentzMap_h; - std::vector<TH2F*> distanceMap_e; - std::vector<TH2F*> distanceMap_h; - + std::vector<PixelHistoConverter> ramoPotentialMap; + std::vector<PixelHistoConverter> lorentzMap_e; + std::vector<PixelHistoConverter> lorentzMap_h; + std::vector<PixelHistoConverter> distanceMap_e; + std::vector<PixelHistoConverter> distanceMap_h; for (unsigned int i=0; i<mapsPath_list.size(); i++) { ATH_MSG_INFO("Using maps located in: "<<mapsPath_list.at(i) << " for layer No." << i); TFile* mapsFile = new TFile((mapsPath_list.at(i)).c_str()); //this is the ramo potential. @@ -476,12 +475,17 @@ StatusCode PixelConfigCondAlg::execute(const EventContext& ctx) const { ATH_MSG_FATAL("Did not find a Ramo potential map and an approximate form is available yet. Exit..."); return StatusCode::FAILURE; } - ramoPotentialMap.push_back(ramoPotentialMap_hold); - lorentzMap_e.push_back((TH2F*)mapsFile->Get("lorentz_map_e")); - lorentzMap_h.push_back((TH2F*)mapsFile->Get("lorentz_map_h")); - distanceMap_e.push_back((TH2F*)mapsFile->Get("edistance")); - distanceMap_h.push_back((TH2F*)mapsFile->Get("hdistance")); + ramoPotentialMap.emplace_back(); + ATH_CHECK(ramoPotentialMap.back().setHisto3D(ramoPotentialMap_hold)); + lorentzMap_e.emplace_back(); + lorentzMap_h.emplace_back(); + distanceMap_e.emplace_back(); + distanceMap_h.emplace_back(); + ATH_CHECK(lorentzMap_e.back().setHisto2D((TH2F*)mapsFile->Get("lorentz_map_e"))); + ATH_CHECK(lorentzMap_h.back().setHisto2D((TH2F*)mapsFile->Get("lorentz_map_h"))); + ATH_CHECK(distanceMap_e.back().setHisto2D((TH2F*)mapsFile->Get("edistance"))); + ATH_CHECK(distanceMap_h.back().setHisto2D((TH2F*)mapsFile->Get("hdistance"))); } writeCdo -> setLorentzMap_e(lorentzMap_e); writeCdo -> setLorentzMap_h(lorentzMap_h); diff --git a/InnerDetector/InDetConditions/PixelConditionsData/PixelConditionsData/PixelHistoConverter.h b/InnerDetector/InDetConditions/PixelConditionsData/PixelConditionsData/PixelHistoConverter.h new file mode 100644 index 0000000000000000000000000000000000000000..13e03ce4813e87b0d7d8a9b0ef61c10622bb3a44 --- /dev/null +++ b/InnerDetector/InDetConditions/PixelConditionsData/PixelConditionsData/PixelHistoConverter.h @@ -0,0 +1,77 @@ +/* + Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration + */ +/** + * @file PixelDigitization/PixelHistoConverter.h + * @author Tomas Dado <tomas.dado@cern.ch> + * @date February, 2021 + * @brief A helper class to conver ROOT TH* objects into simple vectors for speed up + */ + +#ifndef PIXELCONDITIONSDATA_PIXELHISTOCONVRTER_H +#define PIXELCONDITIONSDATA_PIXELHISTOCONVRTER_H + +#include "AthenaBaseComps/AthMessaging.h" +#include "GaudiKernel/StatusCode.h" + +#include <vector> + +/// forward class declaration +class TAxis; +class TH1; +class TH2; +class TH3; + +class PixelHistoConverter{ +public: + PixelHistoConverter(); + virtual ~PixelHistoConverter() = default; + + StatusCode setHisto1D(const TH1* histo); + StatusCode setHisto2D(const TH2* histo); + StatusCode setHisto3D(const TH3* histo); + + inline float getContent(const std::size_t& x) const { + return m_content[x]; + } + inline float getContent(const std::size_t& x, const std::size_t& y) const { + return m_content[x + y*(m_xAxis.nBins)]; + } + inline float getContent(const std::size_t& x, const std::size_t& y, const std::size_t& z) const { + return m_content[x + m_xAxis.nBins*(y + (m_yAxis.nBins * z))]; + } + + inline bool isOverflowZ(const float value) const { + return (value >= m_zAxis.max); + } + bool isFirstZ(const float value) const; + float getBinX(const float value) const; + float getBinY(const float value) const; + float getBinZ(const float value) const; + +private: + struct Axis { + std::size_t nBins; + float min; + float max; + float width; + }; + + Axis m_xAxis; + Axis m_yAxis; + Axis m_zAxis; + + std::vector<float> m_content; + + bool setAxis(Axis& axis, const TAxis* rootAxis); + + inline std::size_t findBin(const Axis& axis, const float value) const { + if (value <= axis.min) return 0; + if (value >= axis.max) return (axis.nBins - 1); + + return ((value - axis.min) * axis.width); + } + +}; + +#endif diff --git a/InnerDetector/InDetConditions/PixelConditionsData/PixelConditionsData/PixelModuleData.h b/InnerDetector/InDetConditions/PixelConditionsData/PixelConditionsData/PixelModuleData.h index e1efb342b00c5038223ed7eb4917f80083297610..4d9fa5c2fad7ed4734295c0193027ecb82a1a6ee 100644 --- a/InnerDetector/InDetConditions/PixelConditionsData/PixelConditionsData/PixelModuleData.h +++ b/InnerDetector/InDetConditions/PixelConditionsData/PixelConditionsData/PixelModuleData.h @@ -1,5 +1,5 @@ /* - Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration */ /** * @file PixelConditionsData/PixelModuleData.h @@ -15,6 +15,7 @@ #include <map> #include "AthenaKernel/CondCont.h" +#include "PixelConditionsData/PixelHistoConverter.h" #include "TH1.h" #include "TH2.h" #include "TH3.h" @@ -167,18 +168,18 @@ class PixelModuleData { void setFluenceLayer(std::vector<double> fluenceLayer); double getFluenceLayer(int layer) const; - void setLorentzMap_e(std::vector<TH2F*> lorentzMap_e); - void setLorentzMap_h(std::vector<TH2F*> lorentzMap_h); - TH2F* getLorentzMap_e(int layer) const; - TH2F* getLorentzMap_h(int layer) const; + void setLorentzMap_e(std::vector<PixelHistoConverter> lorentzMap_e); + void setLorentzMap_h(std::vector<PixelHistoConverter> lorentzMap_h); + const PixelHistoConverter& getLorentzMap_e(int layer) const; + const PixelHistoConverter& getLorentzMap_h(int layer) const; - void setDistanceMap_e(std::vector<TH2F*> distanceMap_e); - void setDistanceMap_h(std::vector<TH2F*> distanceMap_h); - TH2F* getDistanceMap_e(int layer) const; - TH2F* getDistanceMap_h(int layer) const; + void setDistanceMap_e(std::vector<PixelHistoConverter> distanceMap_e); + void setDistanceMap_h(std::vector<PixelHistoConverter> distanceMap_h); + const PixelHistoConverter& getDistanceMap_e(int layer) const; + const PixelHistoConverter& getDistanceMap_h(int layer) const; - void setRamoPotentialMap(std::vector<TH3F*> ramoPotentialMap); - TH3F* getRamoPotentialMap(int layer) const; + void setRamoPotentialMap(std::vector<PixelHistoConverter> ramoPotentialMap); + const PixelHistoConverter& getRamoPotentialMap(int layer) const; // Distortion parameters void setDistortionInputSource(int distortionInputSource); @@ -302,11 +303,11 @@ class PixelModuleData { std::string m_cablingMapFileName; std::vector<double> m_fluenceLayer; - std::vector<TH2F*> m_lorentzMap_e; - std::vector<TH2F*> m_lorentzMap_h; - std::vector<TH2F*> m_distanceMap_e; - std::vector<TH2F*> m_distanceMap_h; - std::vector<TH3F*> m_ramoPotentialMap; + std::vector<PixelHistoConverter> m_lorentzMap_e; + std::vector<PixelHistoConverter> m_lorentzMap_h; + std::vector<PixelHistoConverter> m_distanceMap_e; + std::vector<PixelHistoConverter> m_distanceMap_h; + std::vector<PixelHistoConverter> m_ramoPotentialMap; int m_distortionInputSource; int m_distortionVersion; diff --git a/InnerDetector/InDetConditions/PixelConditionsData/src/PixelHistoConverter.cxx b/InnerDetector/InDetConditions/PixelConditionsData/src/PixelHistoConverter.cxx new file mode 100644 index 0000000000000000000000000000000000000000..70ffe84b4d151b5dedc9bb8aa604047f71c5b075 --- /dev/null +++ b/InnerDetector/InDetConditions/PixelConditionsData/src/PixelHistoConverter.cxx @@ -0,0 +1,145 @@ +/* + Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration + */ + +#include "PixelConditionsData/PixelHistoConverter.h" + +#include "StoreGate/StoreGateSvc.h" +#include "GaudiKernel/ISvcLocator.h" + +#include "TAxis.h" +#include "TH1.h" +#include "TH2.h" +#include "TH3.h" + +#include <iostream> + +PixelHistoConverter::PixelHistoConverter() +{ +} + +StatusCode PixelHistoConverter::setHisto1D(const TH1* histo) { + + if (!histo) { + return StatusCode::FAILURE; + } + + if (!setAxis(m_xAxis, histo->GetXaxis())) { + return StatusCode::FAILURE; + } + + /// fill the content + const std::size_t xSize = m_xAxis.nBins; + m_content.resize(xSize); + + for (std::size_t x = 0; x < xSize; ++x) { + m_content.at(x) = histo->GetBinContent(x+1); + } + + return StatusCode::SUCCESS; +} + +StatusCode PixelHistoConverter::setHisto2D(const TH2* histo) { + if (!histo) { + return StatusCode::FAILURE; + } + + if (!setAxis(m_xAxis, histo->GetXaxis())) { + return StatusCode::FAILURE; + } + if (!setAxis(m_yAxis, histo->GetYaxis())) { + return StatusCode::FAILURE; + } + + /// fill the content use linearized version for performance reasons + const std::size_t xSize = m_xAxis.nBins; + const std::size_t ySize = m_yAxis.nBins; + m_content.resize(xSize*ySize); + + for (std::size_t x = 0; x < xSize; ++x) { + for (std::size_t y = 0; y < ySize; ++y) { + const std::size_t position = x + y*xSize; + m_content.at(position) = histo->GetBinContent(x+1,y+1); + } + } + + return StatusCode::SUCCESS; +} + +StatusCode PixelHistoConverter::setHisto3D(const TH3* histo) { + if (!histo) { + return StatusCode::FAILURE; + } + + if (!setAxis(m_xAxis, histo->GetXaxis())) { + return StatusCode::FAILURE; + } + if (!setAxis(m_yAxis, histo->GetYaxis())) { + return StatusCode::FAILURE; + } + if (!setAxis(m_zAxis, histo->GetZaxis())) { + return StatusCode::FAILURE; + } + + /// fill the content use linearized version for performance reasons + const std::size_t xSize = m_xAxis.nBins; + const std::size_t ySize = m_yAxis.nBins; + const std::size_t zSize = m_zAxis.nBins; + m_content.resize(xSize*ySize*zSize); + + for (std::size_t x = 0; x < xSize; ++x) { + for (std::size_t y = 0; y < ySize; ++y) { + for (std::size_t z = 0; z < zSize; ++z) { + const std::size_t position = x + xSize*(y + (ySize * z)); + m_content.at(position) = histo->GetBinContent(x+1,y+1,z+1); + } + } + } + + return StatusCode::SUCCESS; +} + +bool PixelHistoConverter::isFirstZ(const float value) const { + return (getBinZ(value) == 0); +} + +float PixelHistoConverter::getBinX(const float value) const { + return findBin(m_xAxis, value); +} + +float PixelHistoConverter::getBinY(const float value) const { + return findBin(m_yAxis, value); +} + +float PixelHistoConverter::getBinZ(const float value) const { + return findBin(m_zAxis, value); +} + +bool PixelHistoConverter::setAxis(Axis& axis, const TAxis* rootAxis) { + + if (!rootAxis) { + return false; + } + + axis.nBins = rootAxis->GetNbins(); + axis.min = rootAxis->GetXmin(); + axis.max = rootAxis->GetXmax(); + + if (axis.nBins < 1) { + return false; + } + + /// check if the histogram has equidistant bins + const float width = rootAxis->GetBinWidth(1); + for (std::size_t ibin = 2; ibin <= axis.nBins; ++ibin) { + /// use a threshold for imperfect binning + if (std::abs(rootAxis->GetBinWidth(ibin) - width) > 0.01*width) { + return false; + } + } + + /// storing as 1/width to avoid (slow) division in retrieving + axis.width = 1.*axis.nBins/(axis.max - axis.min); + + return true; +} diff --git a/InnerDetector/InDetConditions/PixelConditionsData/src/PixelModuleData.cxx b/InnerDetector/InDetConditions/PixelConditionsData/src/PixelModuleData.cxx index 16da803924fcd45e2b570b9451e491b4fdcceaea..b815435be8486b168abe3aa5f2664ea6e0add5ac 100644 --- a/InnerDetector/InDetConditions/PixelConditionsData/src/PixelModuleData.cxx +++ b/InnerDetector/InDetConditions/PixelConditionsData/src/PixelModuleData.cxx @@ -1,5 +1,5 @@ /* - Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration */ #include "PixelConditionsData/PixelModuleData.h" @@ -319,19 +319,18 @@ std::string PixelModuleData::getCablingMapFileName() const { return m_cablingMap void PixelModuleData::setFluenceLayer(std::vector<double> fluenceLayer) { m_fluenceLayer = fluenceLayer; } double PixelModuleData::getFluenceLayer(int layer) const { return m_fluenceLayer.at(layer); } -void PixelModuleData::setLorentzMap_e(std::vector<TH2F*> lorentzMap_e) { m_lorentzMap_e = lorentzMap_e; } -void PixelModuleData::setLorentzMap_h(std::vector<TH2F*> lorentzMap_h) { m_lorentzMap_h = lorentzMap_h; } -TH2F* PixelModuleData::getLorentzMap_e(int layer) const { return m_lorentzMap_e.at(layer); } -TH2F* PixelModuleData::getLorentzMap_h(int layer) const { return m_lorentzMap_h.at(layer); } +void PixelModuleData::setLorentzMap_e(std::vector<PixelHistoConverter> lorentzMap_e) { m_lorentzMap_e = lorentzMap_e; } +void PixelModuleData::setLorentzMap_h(std::vector<PixelHistoConverter> lorentzMap_h) { m_lorentzMap_h = lorentzMap_h; } +const PixelHistoConverter& PixelModuleData::getLorentzMap_e(int layer) const { return m_lorentzMap_e.at(layer); } +const PixelHistoConverter& PixelModuleData::getLorentzMap_h(int layer) const { return m_lorentzMap_h.at(layer); } -void PixelModuleData::setDistanceMap_e(std::vector<TH2F*> distanceMap_e) { m_distanceMap_e = distanceMap_e; } -void PixelModuleData::setDistanceMap_h(std::vector<TH2F*> distanceMap_h) { m_distanceMap_h = distanceMap_h; } -TH2F* PixelModuleData::getDistanceMap_e(int layer) const { return m_distanceMap_e.at(layer); } -TH2F* PixelModuleData::getDistanceMap_h(int layer) const { return m_distanceMap_h.at(layer); } - -void PixelModuleData::setRamoPotentialMap(std::vector<TH3F*> ramoPotentialMap) { m_ramoPotentialMap = ramoPotentialMap; } -TH3F* PixelModuleData::getRamoPotentialMap(int layer) const { return m_ramoPotentialMap.at(layer); } +void PixelModuleData::setDistanceMap_e(std::vector<PixelHistoConverter> distanceMap_e) { m_distanceMap_e = distanceMap_e; } +void PixelModuleData::setDistanceMap_h(std::vector<PixelHistoConverter> distanceMap_h) { m_distanceMap_h = distanceMap_h; } +const PixelHistoConverter& PixelModuleData::getDistanceMap_e(int layer) const { return m_distanceMap_e.at(layer); } +const PixelHistoConverter& PixelModuleData::getDistanceMap_h(int layer) const { return m_distanceMap_h.at(layer); } +void PixelModuleData::setRamoPotentialMap(std::vector<PixelHistoConverter> ramoPotentialMap) { m_ramoPotentialMap = ramoPotentialMap; } +const PixelHistoConverter& PixelModuleData::getRamoPotentialMap(int layer) const { return m_ramoPotentialMap.at(layer); } // Distortion parameters void PixelModuleData::setDistortionInputSource(int distortionInputSource) { m_distortionInputSource = distortionInputSource; } int PixelModuleData::getDistortionInputSource() const { return m_distortionInputSource; } diff --git a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.cxx b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.cxx index 7e1b009b4afd5a985d56568952a6424b070cc9b8..6e36ec275af73a69a9bb84c73302ce886dcc12aa 100644 --- a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.cxx +++ b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.cxx @@ -1,5 +1,5 @@ /* - Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration */ #include "SensorSim3DTool.h" @@ -105,8 +105,8 @@ StatusCode SensorSim3DTool::initialize() { return StatusCode::FAILURE; //Obviously, remove this when gen. code is set up } //ramoPotentialMap.push_back(ramoPotentialMap_hold); - m_ramoPotentialMap[Layer] = ramoPotentialMap_hold; - m_fluence_layersMaps[Layer] = m_fluence_layers.at(i); + m_ramoPotentialMap.emplace_back(); + ATH_CHECK(m_ramoPotentialMap.back().setHisto3D(ramoPotentialMap_hold)); ATH_MSG_INFO("Using fluence " << m_fluence_layers.at(i)); @@ -118,7 +118,8 @@ StatusCode SensorSim3DTool::initialize() { ATH_MSG_INFO("Unable to load sensor e-field map, so generating one using approximations."); return StatusCode::FAILURE; //Obviously, remove this when gen. code is set up } - m_eFieldMap[Layer] = eFieldMap_hold; + m_eFieldMap.emplace_back(); + ATH_CHECK(m_eFieldMap.back().setHisto2D(eFieldMap_hold)); TH3F* xPositionMap_e_hold; TH3F* xPositionMap_h_hold; @@ -140,26 +141,22 @@ StatusCode SensorSim3DTool::initialize() { timeMap_e_hold = (TH2F* ) mapsFile->Get("etimes"); timeMap_h_hold = (TH2F* ) mapsFile->Get("htimes"); //Now, determine the time to reach the electrode and the trapping position. - if (xPositionMap_e_hold == 0 || xPositionMap_h_hold == 0 || yPositionMap_e_hold == 0 || yPositionMap_h_hold == 0 || - timeMap_e_hold == 0 || timeMap_h_hold == 0) { - ATH_MSG_INFO("Unable to load at least one of teh distance/time maps, so generating all using approximations."); - return StatusCode::FAILURE; //Obviously, remove this when gen. code is set up - } - m_xPositionMap_e[Layer] = xPositionMap_e_hold; - m_xPositionMap_h[Layer] = xPositionMap_h_hold; - m_yPositionMap_e[Layer] = yPositionMap_e_hold; - m_yPositionMap_h[Layer] = yPositionMap_h_hold; - m_timeMap_e[Layer] = timeMap_e_hold; - m_timeMap_h[Layer] = timeMap_h_hold; + m_xPositionMap_e.emplace_back(); + m_xPositionMap_h.emplace_back(); + m_yPositionMap_e.emplace_back(); + m_yPositionMap_h.emplace_back(); + ATH_CHECK(m_xPositionMap_e.back().setHisto3D(xPositionMap_e_hold)); + ATH_CHECK(m_xPositionMap_h.back().setHisto3D(xPositionMap_h_hold)); + ATH_CHECK(m_yPositionMap_e.back().setHisto3D(yPositionMap_e_hold)); + ATH_CHECK(m_yPositionMap_h.back().setHisto3D(yPositionMap_h_hold)); + m_timeMap_e.emplace_back(); + m_timeMap_h.emplace_back(); + ATH_CHECK(m_timeMap_e.back().setHisto2D(timeMap_e_hold)); + ATH_CHECK(m_timeMap_h.back().setHisto2D(timeMap_h_hold)); // Get average charge data (for charge chunk effect correction) - m_avgChargeMap_e = 0; - m_avgChargeMap_h = 0; - m_avgChargeMap_e = (TH2F*) mapsFile->Get("avgCharge_e"); - m_avgChargeMap_h = (TH2F*) mapsFile->Get("avgCharge_h"); - if (m_avgChargeMap_e == 0 || m_avgChargeMap_h == 0) { - ATH_MSG_INFO("Unsuccessful picking up histogram: m_avgChargeMap"); - } + ATH_CHECK(m_avgChargeMap_e.setHisto2D((TH2F*) mapsFile->Get("avgCharge_e"))); + ATH_CHECK(m_avgChargeMap_h.setHisto2D((TH2F*) mapsFile->Get("avgCharge_h"))); } return StatusCode::SUCCESS; @@ -359,13 +356,8 @@ StatusCode SensorSim3DTool::induceCharge(const TimedHitPtr<SiHit>& phit, ATH_MSG_DEBUG(" -- diffused position w.r.t. pixel edge = " << xposDiff << " " << yposDiff); - double average_charge = m_avgChargeMap_e->GetBinContent(m_avgChargeMap_e->GetYaxis()->FindBin( - y_pix * 1000), - m_avgChargeMap_e->GetXaxis()->FindBin(x_pix * - 1000)); - if (isHole) average_charge = - m_avgChargeMap_h->GetBinContent(m_avgChargeMap_h->GetYaxis()->FindBin( - y_pix * 1000), m_avgChargeMap_h->GetXaxis()->FindBin(x_pix * 1000)); + float average_charge = isHole ? m_avgChargeMap_h.getContent(m_avgChargeMap_h.getBinY(1e3*y_pix), m_avgChargeMap_h.getBinX(1e3*x_pix)) : + m_avgChargeMap_e.getContent(m_avgChargeMap_e.getBinY(1e3*y_pix), m_avgChargeMap_e.getBinX(1e3*x_pix)); ATH_MSG_DEBUG(" -- driftTime, timeToElectrode = " << driftTime << " " << timeToElectrode); @@ -379,6 +371,9 @@ StatusCode SensorSim3DTool::induceCharge(const TimedHitPtr<SiHit>& phit, for (int i = -1; i <= 1; i++) { double xNeighbor = i * pixel_size_x; // -- loop in the y-coordinate + const std::size_t index = 0; + const std::size_t ramo_init_bin_y = m_ramoPotentialMap[index].getBinY(1000*(x_pix + pixel_size_x * 3 - xNeighbor)); + const std::size_t ramo_final_bin_y = m_ramoPotentialMap[index].getBinY(1000*(xposFinal + pixel_size_x * 3 - xNeighbor)); for (int j = -1; j <= 1; j++) { double yNeighbor = j * pixel_size_y; @@ -391,10 +386,8 @@ StatusCode SensorSim3DTool::induceCharge(const TimedHitPtr<SiHit>& phit, //Ramo map over 500umx350um pixel area //Ramo init different if charge diffused into neighboring pixel -> change primary pixel!! - double ramoInit = getRamoPotential(y_pix + pixel_size_y * 1 / 2 - yNeighbor, - x_pix + pixel_size_x * 3 - xNeighbor); - double ramoFinal = getRamoPotential(yposFinal + pixel_size_y * 1 / 2 - yNeighbor, - xposFinal + pixel_size_x * 3 - xNeighbor); + float ramoInit = m_ramoPotentialMap[index].getContent(m_ramoPotentialMap[index].getBinX(1000*(y_pix + 0.5*pixel_size_y - yNeighbor)), ramo_init_bin_y); + float ramoFinal = m_ramoPotentialMap[index].getContent(m_ramoPotentialMap[index].getBinX(1000*(yposFinal + 0.5*pixel_size_y - yNeighbor)), ramo_final_bin_y); // Record deposit double eHitRamo = (1 - 2 * isHole) * eHit * (ramoFinal - ramoInit); @@ -508,8 +501,6 @@ StatusCode SensorSim3DTool::induceCharge(const TimedHitPtr<SiHit>& phit, double ed = es_current * eleholePairEnergy * ccprob_neighbor; // -- pixel coordinates --> module coordinates - //double x_mod = x_neighbor - half_pixel_size_x + pixel_size_x*nPixX -module_size_x/2.; - //double y_mod = y_neighbor - half_pixel_size_y + pixel_size_y*nPixY -module_size_y/2.; double x_mod = x_neighbor + pixel_size_x / 2 + pixel_size_x * nPixX - module_size_x / 2.; double y_mod = y_neighbor + pixel_size_y / 2 + pixel_size_y * nPixY - module_size_y / 2.; SiLocalPosition chargePos = Module.hitLocalToLocal(y_mod, x_mod); @@ -531,7 +522,7 @@ StatusCode SensorSim3DTool::induceCharge(const TimedHitPtr<SiHit>& phit, } // read the Charge Collection Prob Map from text file -StatusCode SensorSim3DTool::readProbMap(std::string fileE) { +StatusCode SensorSim3DTool::readProbMap(const std::string& fileE) { std::string line; const std::string fileName = fileE; std::string inputFile = PathResolverFindCalibFile(fileName); @@ -567,20 +558,18 @@ StatusCode SensorSim3DTool::readProbMap(std::string fileE) { } // -- Print out the Charge Collection Probability map (full map) -StatusCode SensorSim3DTool::printProbMap(std::string readout) { +StatusCode SensorSim3DTool::printProbMap(const std::string& readout) const { if (readout == "FEI4") { - for (std::multimap<std::pair<int, int>, double >::iterator it = m_probMapFEI4.begin(); it != m_probMapFEI4.end(); - ++it) { + for (const auto& it : m_probMapFEI4) { ATH_MSG_DEBUG( - "read full probMap FEI4 --- bin x " << it->first.first << " bin y " << it->first.second << " prob " << - it->second); + "read full probMap FEI4 --- bin x " << it.first.first << " bin y " << it.first.second << " prob " << + it.second); } } else if (readout == "FEI3") { - for (std::multimap<std::pair<int, int>, double >::iterator it = m_probMapFEI3.begin(); it != m_probMapFEI3.end(); - ++it) { + for (const auto& it : m_probMapFEI3) { ATH_MSG_DEBUG( - "read full probMap FEI3 --- bin x " << it->first.first << " bin y " << it->first.second << " prob " << - it->second); + "read full probMap FEI3 --- bin x " << it.first.first << " bin y " << it.first.second << " prob " << + it.second); } } else { ATH_MSG_ERROR("Error in printout Charge Coll Prob Maps! (readout should contain FEI3 or FEI4 strings) "); @@ -590,7 +579,7 @@ StatusCode SensorSim3DTool::printProbMap(std::string readout) { } // -- Returns the Charge Collection Probability at a given point (bin_x,bin_y) -double SensorSim3DTool::getProbMapEntry(std::string readout, int binx, int biny) { +double SensorSim3DTool::getProbMapEntry(const std::string& readout, int binx, int biny) const { std::pair<int, int> doublekey(binx, biny); double echarge; if (readout == "FEI4") { @@ -607,15 +596,8 @@ double SensorSim3DTool::getProbMapEntry(std::string readout, int binx, int biny) } double SensorSim3DTool::getElectricField(double x, double y) { - std::pair < int, int > Layer; // index for layer/end cap position - Layer.first = 0; //Barrel (0) or End Cap (1) - Now only for Barrel. If we want to add End Caps, put them at - // Layer.first=1 - Layer.second = 0; //Layer: 0 = IBL Planar, 1=B-Layer, 2=Layer-1, 3=Layer-2 - - int n_binx = m_eFieldMap[Layer]->GetXaxis()->FindBin(x * 1000); //position coordinates in um to return electric field - // in V/cm - int n_biny = m_eFieldMap[Layer]->GetYaxis()->FindBin(y * 1000); - double electricField = m_eFieldMap[Layer]->GetBinContent(n_binx, n_biny); + std::size_t index = 0; + double electricField = m_eFieldMap[index].getContent(m_eFieldMap[index].getBinX(1e3*x), m_eFieldMap[index].getBinY(1e3*y)); return electricField * 1.0E-7; //return efield in MV/mm (for mobility calculation) } @@ -656,79 +638,36 @@ double SensorSim3DTool::getDriftTime(bool isHoleBit) { } double SensorSim3DTool::getTimeToElectrode(double x, double y, bool isHoleBit) { - std::pair < int, int > Layer; // index for layer/end cap position - Layer.first = 0; //Barrel (0) or End Cap (1) - Now only for Barrel. If we want to add End Caps, put them at - // Layer.first=1 - Layer.second = 0; //Layer: 0 = IBL Planar, 1=B-Layer, 2=Layer-1, 3=Layer-2 - - // Uses (x,y) position in um to return time to electrode in ns + std::size_t index = 0; double timeToElectrode = 0; if (!isHoleBit) { - int n_binx = m_timeMap_e[Layer]->GetXaxis()->FindBin(x * 1000); //convert from mm to um - int n_biny = m_timeMap_e[Layer]->GetYaxis()->FindBin(y * 1000); - timeToElectrode = m_timeMap_e[Layer]->GetBinContent(n_binx, n_biny); - } - if (isHoleBit) { - int n_binx = m_timeMap_h[Layer]->GetXaxis()->FindBin(x * 1000); - int n_biny = m_timeMap_h[Layer]->GetYaxis()->FindBin(y * 1000); - timeToElectrode = m_timeMap_h[Layer]->GetBinContent(n_binx, n_biny); + timeToElectrode = m_timeMap_e[index].getContent(m_timeMap_e[index].getBinX(1e3*x), m_timeMap_e[index].getBinY(1e3*y)); + } else { + timeToElectrode = m_timeMap_h[index].getContent(m_timeMap_h[index].getBinX(1e3*x), m_timeMap_h[index].getBinY(1e3*y)); } return timeToElectrode; //[ns] } double SensorSim3DTool::getTrappingPositionX(double initX, double initY, double driftTime, bool isHoleBit) { - std::pair < int, int > Layer; // index for layer/end cap position - Layer.first = 0; //Barrel (0) or End Cap (1) - Now only for Barrel. If we want to add End Caps, put them at - // Layer.first=1 - Layer.second = 0; //Layer: 0 = IBL Planar, 1=B-Layer, 2=Layer-1, 3=Layer-2 - + std::size_t index = 0; double finalX = initX; if (!isHoleBit) { - int n_binx = m_xPositionMap_e[Layer]->GetXaxis()->FindBin(initX * 1000); - int n_biny = m_xPositionMap_e[Layer]->GetYaxis()->FindBin(initY * 1000); - int n_bint = m_xPositionMap_e[Layer]->GetZaxis()->FindBin(driftTime); - finalX = m_xPositionMap_e[Layer]->GetBinContent(n_binx, n_biny, n_bint); + finalX = m_xPositionMap_e[index].getContent(m_xPositionMap_e[index].getBinX(1e3*initX), m_xPositionMap_e[index].getBinY(1e3*initY), m_xPositionMap_e[index].getBinZ(driftTime)); } else { - int n_binx = m_xPositionMap_h[Layer]->GetXaxis()->FindBin(initX * 1000); - int n_biny = m_xPositionMap_h[Layer]->GetYaxis()->FindBin(initY * 1000); - int n_bint = m_xPositionMap_h[Layer]->GetZaxis()->FindBin(driftTime); - finalX = m_xPositionMap_h[Layer]->GetBinContent(n_binx, n_biny, n_bint); + finalX = m_xPositionMap_h[index].getContent(m_xPositionMap_h[index].getBinX(1e3*initX), m_xPositionMap_h[index].getBinY(1e3*initY), m_xPositionMap_h[index].getBinZ(driftTime)); } return finalX * 1e-3; //[mm] } double SensorSim3DTool::getTrappingPositionY(double initX, double initY, double driftTime, bool isHoleBit) { - std::pair < int, int > Layer; // index for layer/end cap position - Layer.first = 0; //Barrel (0) or End Cap (1) - Now only for Barrel. If we want to add End Caps, put them at - // Layer.first=1 - Layer.second = 0; //Layer: 0 = IBL Planar, 1=B-Layer, 2=Layer-1, 3=Layer-2 - + std::size_t index = 0; double finalY = initY; if (!isHoleBit) { - int n_binx = m_yPositionMap_e[Layer]->GetXaxis()->FindBin(initX * 1000); - int n_biny = m_yPositionMap_e[Layer]->GetYaxis()->FindBin(initY * 1000); - int n_bint = m_yPositionMap_e[Layer]->GetZaxis()->FindBin(driftTime); - finalY = m_yPositionMap_e[Layer]->GetBinContent(n_binx, n_biny, n_bint); + finalY = m_yPositionMap_e[index].getContent(m_yPositionMap_e[index].getBinX(1e3*initX), m_yPositionMap_e[index].getBinY(1e3*initY), m_yPositionMap_e[index].getBinZ(driftTime)); } else { - int n_binx = m_yPositionMap_h[Layer]->GetXaxis()->FindBin(initX * 1000); - int n_biny = m_yPositionMap_h[Layer]->GetYaxis()->FindBin(initY * 1000); - int n_bint = m_yPositionMap_h[Layer]->GetZaxis()->FindBin(driftTime); - finalY = m_yPositionMap_h[Layer]->GetBinContent(n_binx, n_biny, n_bint); + finalY = m_yPositionMap_h[index].getContent(m_yPositionMap_h[index].getBinX(1e3*initX), m_yPositionMap_h[index].getBinY(1e3*initY), m_yPositionMap_h[index].getBinZ(driftTime)); } return finalY * 1e-3; //[mm] } - -double SensorSim3DTool::getRamoPotential(double x, double y) { - std::pair < int, int > Layer; // index for layer/end cap position - Layer.first = 0; //Barrel (0) or End Cap (1) - Now only for Barrel. If we want to add End Caps, put them at - // Layer.first=1 - Layer.second = 0; //Layer: 0 = IBL Planar, 1=B-Layer, 2=Layer-1, 3=Layer-2 - - - int n_binx = m_ramoPotentialMap[Layer]->GetXaxis()->FindBin(x * 1000); - int n_biny = m_ramoPotentialMap[Layer]->GetYaxis()->FindBin(y * 1000); - double ramoPotential = m_ramoPotentialMap[Layer]->GetBinContent(n_binx, n_biny); - return ramoPotential; -} diff --git a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.h b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.h index 803e99eff0c7b17708923731e033f9b703d408b3..ea37734ebf756fcdfd8c4b0c37f0c7fc8b5a0ba7 100644 --- a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.h +++ b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.h @@ -17,6 +17,7 @@ #include "GaudiKernel/ToolHandle.h" #include "RadDamageUtil.h" +#include "PixelConditionsData/PixelHistoConverter.h" class SensorSim3DTool: public SensorSimTool { @@ -38,9 +39,9 @@ public: // 3D sensor simulation using probability density map (used in RUN-2 (no radiation damage) - StatusCode readProbMap(std::string); - StatusCode printProbMap(std::string); - double getProbMapEntry(std::string, int, int); + StatusCode readProbMap(const std::string&); + StatusCode printProbMap(const std::string&) const; + double getProbMapEntry(const std::string&, int, int) const; double getElectricField(double x, double y); double getMobility(double electricField, bool isHoleBit); @@ -48,7 +49,6 @@ public: double getTimeToElectrode(double x, double y, bool isHoleBit); double getTrappingPositionX(double initX, double initY, double driftTime, bool isHoleBit); double getTrappingPositionY(double initX, double initY, double driftTime, bool isHoleBit); - double getRamoPotential(double x, double y); private: SensorSim3DTool(); @@ -57,19 +57,18 @@ private: std::multimap<std::pair<int, int>, double> m_probMapFEI3; // Map for radiation damage simulation - std::map<std::pair<int, int>, TH3F*> m_ramoPotentialMap; - std::map<std::pair<int, int>, TH2F*> m_eFieldMap; - std::map<std::pair<int, int>, TH3F*> m_xPositionMap_e; - std::map<std::pair<int, int>, TH3F*> m_xPositionMap_h; - std::map<std::pair<int, int>, TH3F*> m_yPositionMap_e; - std::map<std::pair<int, int>, TH3F*> m_yPositionMap_h; - std::map<std::pair<int, int>, TH2F*> m_timeMap_e; - std::map<std::pair<int, int>, TH2F*> m_timeMap_h; - TH2F* m_avgChargeMap_e; - TH2F* m_avgChargeMap_h; + std::vector<PixelHistoConverter> m_ramoPotentialMap; + std::vector<PixelHistoConverter> m_eFieldMap; + std::vector<PixelHistoConverter> m_xPositionMap_e; + std::vector<PixelHistoConverter> m_xPositionMap_h; + std::vector<PixelHistoConverter> m_yPositionMap_e; + std::vector<PixelHistoConverter> m_yPositionMap_h; + std::vector<PixelHistoConverter> m_timeMap_e; + std::vector<PixelHistoConverter> m_timeMap_h; + PixelHistoConverter m_avgChargeMap_e; + PixelHistoConverter m_avgChargeMap_h; std::vector<double> m_fluence_layers; - std::map<std::pair<int, int>, double> m_fluence_layersMaps; Gaudi::Property<std::string> m_cc_prob_file_fei3 { diff --git a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.cxx b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.cxx index b6b379f47437aede0817ade2d01b1dce1c9933ac..8622bf6ec9f0bec78d10776b965c10a98b6eebd7 100644 --- a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.cxx +++ b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.cxx @@ -94,7 +94,8 @@ StatusCode SensorSimPlanarTool::initialize() { ATH_MSG_WARNING("Not implemented yet - exit"); return StatusCode::FAILURE; //Obviously, remove this when gen. code is set up } - m_ramoPotentialMap.push_back(ramoPotentialMap_hold); + m_ramoPotentialMap.emplace_back(); + ATH_CHECK(m_ramoPotentialMap.back().setHisto3D(ramoPotentialMap_hold)); //Now setup the E-field. TH1F* eFieldMap_hold; eFieldMap_hold = new TH1F(); @@ -135,13 +136,15 @@ StatusCode SensorSimPlanarTool::initialize() { ATH_MSG_INFO("Unable to load at least one of the distance/time/Lorentz angle maps."); return StatusCode::FAILURE;//Obviously, remove this when gen. code is set up } - m_distanceMap_e.push_back(distanceMap_e_hold); - m_distanceMap_h.push_back(distanceMap_h_hold); - m_lorentzMap_e.push_back(lorentzMap_e_hold); - m_lorentzMap_h.push_back(lorentzMap_h_hold); + m_distanceMap_e.emplace_back(); + m_distanceMap_h.emplace_back(); + ATH_CHECK(m_distanceMap_e.back().setHisto2D(distanceMap_e_hold)); + ATH_CHECK(m_distanceMap_h.back().setHisto2D(distanceMap_h_hold)); + m_lorentzMap_e.emplace_back(); + m_lorentzMap_h.emplace_back(); + ATH_CHECK(m_lorentzMap_e.back().setHisto2D(lorentzMap_e_hold)); + ATH_CHECK(m_lorentzMap_h.back().setHisto2D(lorentzMap_h_hold)); } - // Define necessary variables to describe hists stored in ramoPotentialMap - // to skip calling FindBin during induceCharge for speed reasons } return StatusCode::SUCCESS; } @@ -248,16 +251,10 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit, SiLocalPosition centreOfPixel_i; - int nnLoop_pixelEtaMax = 0; - int nnLoop_pixelEtaMin = 0; - int nnLoop_pixelPhiMax = 0; - int nnLoop_pixelPhiMin = 0; - - int numBins_driftTime_e = 0; - int numBins_driftTime_h = 0; - int numBins_weightingPotential_x = 0; - int numBins_weightingPotential_y = 0; - int numBins_weightingPotential_z = 0; + int nnLoop_pixelEtaMax(0); + int nnLoop_pixelEtaMin(0); + int nnLoop_pixelPhiMax(0); + int nnLoop_pixelPhiMin(0); if (m_doRadDamage && !(Module.isDBM()) && Module.isBarrel()) { centreOfPixel_i = p_design.positionFromColumnRow(pixel_i.etaIndex(), pixel_i.phiIndex()); @@ -268,21 +265,6 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit, nnLoop_pixelPhiMax = std::min(1, pixel_i.phiIndex()); nnLoop_pixelPhiMin = std::max(-1, pixel_i.phiIndex() + 1 - phiCells); - - //Setup values to check for overflow when using maps - if (m_doInterpolateEfield) { - numBins_driftTime_e = m_distanceMap_e[layer]->GetNbinsY(); //Returns nBins = totalBins - underflow - overflow - numBins_driftTime_h = m_distanceMap_h[layer]->GetNbinsY(); - numBins_weightingPotential_x = m_ramoPotentialMap[layer]->GetNbinsX(); - numBins_weightingPotential_y = m_ramoPotentialMap[layer]->GetNbinsY(); - numBins_weightingPotential_z = m_ramoPotentialMap[layer]->GetNbinsZ(); - } else { // use fluence value from conditions data - numBins_driftTime_e = moduleData->getDistanceMap_e(layer)->GetNbinsY(); - numBins_driftTime_h = moduleData->getDistanceMap_h(layer)->GetNbinsY(); - numBins_weightingPotential_x = moduleData->getRamoPotentialMap(layer)->GetNbinsX(); - numBins_weightingPotential_y = moduleData->getRamoPotentialMap(layer)->GetNbinsY(); - numBins_weightingPotential_z = moduleData->getRamoPotentialMap(layer)->GetNbinsZ(); - } } // Distance between charge and readout side. p_design->readoutSide() is @@ -298,33 +280,10 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit, nontrappingProbability = exp(-dist_electrode / collectionDist); } - if (m_doInterpolateEfield) { - m_ramo_x_binMap = 1000. * - (m_ramoPotentialMap[layer]->GetNbinsX() / - (m_ramoPotentialMap[layer]->GetXaxis()->GetXmax() - - m_ramoPotentialMap[layer]->GetXaxis()->GetXmin())); - m_ramo_y_binMap = 1000. * - (m_ramoPotentialMap[layer]->GetNbinsY() / - (m_ramoPotentialMap[layer]->GetYaxis()->GetXmax() - - m_ramoPotentialMap[layer]->GetYaxis()->GetXmin())); - m_ramo_z_binMap = 1000. * - (m_ramoPotentialMap[layer]->GetNbinsZ() / - (m_ramoPotentialMap[layer]->GetZaxis()->GetXmax() - - m_ramoPotentialMap[layer]->GetZaxis()->GetXmin())); - } else { - m_ramo_x_binMap = 1000. * - (moduleData->getRamoPotentialMap(layer)->GetNbinsX() / - (moduleData->getRamoPotentialMap(layer)->GetXaxis()->GetXmax() - - moduleData->getRamoPotentialMap(layer)->GetXaxis()->GetXmin())); - m_ramo_y_binMap = 1000. * - (moduleData->getRamoPotentialMap(layer)->GetNbinsY() / - (moduleData->getRamoPotentialMap(layer)->GetYaxis()->GetXmax() - - moduleData->getRamoPotentialMap(layer)->GetYaxis()->GetXmin())); - m_ramo_z_binMap = 1000. * - (moduleData->getRamoPotentialMap(layer)->GetNbinsZ() / - (moduleData->getRamoPotentialMap(layer)->GetZaxis()->GetXmax() - - moduleData->getRamoPotentialMap(layer)->GetZaxis()->GetXmin())); - } + const std::size_t distance_f_e_bin_x = m_doInterpolateEfield ? m_distanceMap_e[layer].getBinX(dist_electrode) : moduleData->getDistanceMap_e(layer).getBinX(dist_electrode); + const std::size_t distance_f_h_bin_x = m_doInterpolateEfield ? m_distanceMap_h[layer].getBinX(dist_electrode) : moduleData->getDistanceMap_h(layer).getBinX(dist_electrode); + const std::size_t tanLorentz_e_bin_x = m_doInterpolateEfield ? m_lorentzMap_e[layer].getBinX(dist_electrode) : moduleData->getLorentzMap_e(layer).getBinX(dist_electrode); + const std::size_t tanLorentz_h_bin_x = m_doInterpolateEfield ? m_lorentzMap_h[layer].getBinX(dist_electrode) : moduleData->getLorentzMap_h(layer).getBinX(dist_electrode); for (int j = 0; j < ncharges; j++) { if (m_doRadDamage && !(Module.isDBM()) && Module.isBarrel()) { @@ -344,39 +303,15 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit, // electrode. //However, if choose (b), will need to reduce granularity of map. if (m_doInterpolateEfield) { - int nbin_z_e_xbin = m_distanceMap_e[layer]->GetXaxis()->FindBin(dist_electrode); - int nbin_z_e_ybin = m_distanceMap_e[layer]->GetYaxis()->FindBin(drifttime_e); - if (nbin_z_e_ybin > numBins_driftTime_e) { - nbin_z_e_ybin = numBins_driftTime_e; - } - depth_f_e = m_distanceMap_e[layer]->GetBinContent(nbin_z_e_xbin, nbin_z_e_ybin); - int nbin_z_h_xbin = m_distanceMap_h[layer]->GetXaxis()->FindBin(dist_electrode); - int nbin_z_h_ybin = m_distanceMap_h[layer]->GetYaxis()->FindBin(drifttime_h); - if (nbin_z_h_ybin > numBins_driftTime_h) { - nbin_z_h_ybin = numBins_driftTime_h; - } - depth_f_h = m_distanceMap_h[layer]->GetBinContent(nbin_z_h_xbin, nbin_z_h_ybin); - int nbin_Lorentz_e = m_lorentzMap_e[layer]->FindBin(dist_electrode, depth_f_e); - tanLorentz_e = m_lorentzMap_e[layer]->GetBinContent(nbin_Lorentz_e); - int nbin_Lorentz_h = m_lorentzMap_h[layer]->FindBin(dist_electrode, depth_f_h); - tanLorentz_h = m_lorentzMap_h[layer]->GetBinContent(nbin_Lorentz_h); + depth_f_e = m_distanceMap_e[layer].getContent(distance_f_e_bin_x, m_distanceMap_e[layer].getBinY(drifttime_e)); + depth_f_h = m_distanceMap_h[layer].getContent(distance_f_h_bin_x, m_distanceMap_h[layer].getBinY(drifttime_h)); + tanLorentz_e = m_lorentzMap_e[layer].getContent(tanLorentz_e_bin_x, m_lorentzMap_e[layer].getBinY(depth_f_e)); + tanLorentz_h = m_lorentzMap_h[layer].getContent(tanLorentz_h_bin_x, m_lorentzMap_h[layer].getBinY(depth_f_h)); } else { // use fluence value from conditions data - int nbin_z_e_xbin = moduleData->getDistanceMap_e(layer)->GetXaxis()->FindBin(dist_electrode); - int nbin_z_e_ybin = moduleData->getDistanceMap_e(layer)->GetYaxis()->FindBin(drifttime_e); - if (nbin_z_e_ybin > numBins_driftTime_e) { - nbin_z_e_ybin = numBins_driftTime_e; - } - depth_f_e = moduleData->getDistanceMap_e(layer)->GetBinContent(nbin_z_e_xbin, nbin_z_e_ybin); - int nbin_z_h_xbin = moduleData->getDistanceMap_h(layer)->GetXaxis()->FindBin(dist_electrode); - int nbin_z_h_ybin = moduleData->getDistanceMap_h(layer)->GetYaxis()->FindBin(drifttime_h); - if (nbin_z_h_ybin > numBins_driftTime_h) { - nbin_z_h_ybin = numBins_driftTime_h; - } - depth_f_h = moduleData->getDistanceMap_h(layer)->GetBinContent(nbin_z_h_xbin, nbin_z_h_ybin); - int nbin_Lorentz_e = moduleData->getLorentzMap_e(layer)->FindBin(dist_electrode, depth_f_e); - tanLorentz_e = moduleData->getLorentzMap_e(layer)->GetBinContent(nbin_Lorentz_e); - int nbin_Lorentz_h = moduleData->getLorentzMap_h(layer)->FindBin(dist_electrode, depth_f_h); - tanLorentz_h = moduleData->getLorentzMap_h(layer)->GetBinContent(nbin_Lorentz_h); + depth_f_e = moduleData->getDistanceMap_e(layer).getContent(distance_f_e_bin_x, moduleData->getDistanceMap_e(layer).getBinY(drifttime_e)); + depth_f_h = moduleData->getDistanceMap_h(layer).getContent(distance_f_h_bin_x, moduleData->getDistanceMap_h(layer).getBinY(drifttime_h)); + tanLorentz_e = moduleData->getLorentzMap_e(layer).getContent(tanLorentz_e_bin_x, moduleData->getLorentzMap_e(layer).getBinY(depth_f_e)); + tanLorentz_h = moduleData->getLorentzMap_h(layer).getContent(tanLorentz_h_bin_x, moduleData->getLorentzMap_h(layer).getBinY(depth_f_h)); } double dz_e = fabs(dist_electrode - depth_f_e); double dz_h = fabs(depth_f_h - dist_electrode); @@ -397,7 +332,7 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit, double phi_f_h = phi_i + dz_h * tanLorentz_h + rdif_h * phiRand; etaRand = CLHEP::RandGaussZiggurat::shoot(rndmEngine); double eta_f_h = eta_i + rdif_h * etaRand; - + // amount of energy to be converted into charges at current step double energy_per_step = 1.0 * iHitRecord.second / 1.E+6 / ncharges; @@ -407,21 +342,11 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit, ATH_CHECK(applyIBLSlimEdges(energy_per_step, eta_f_h)); } - int nbin_ramo_f_e_z; - int nbin_ramo_f_h_z; - // distinction necessary because of min(z) = -0.5 - float minz = 1; - if (layer != 0) minz = 1.5; - nbin_ramo_f_e_z = int( minz + depth_f_e * m_ramo_z_binMap ); - nbin_ramo_f_h_z = int( minz + depth_f_h * m_ramo_z_binMap ); - - // Check for overflow in ramo hists in z-direction - if (nbin_ramo_f_h_z > numBins_weightingPotential_z) { - nbin_ramo_f_h_z = numBins_weightingPotential_z + 1; - } - if (nbin_ramo_f_e_z > numBins_weightingPotential_z) { - nbin_ramo_f_e_z = numBins_weightingPotential_z + 1; - } + const std::size_t ramo_f_e_bin_z = m_doInterpolateEfield ? m_ramoPotentialMap[layer].getBinZ(1e3*depth_f_e) : moduleData->getRamoPotentialMap(layer).getBinZ(1e3*depth_f_e); + const std::size_t ramo_f_h_bin_z = m_doInterpolateEfield ? m_ramoPotentialMap[layer].getBinZ(1e3*depth_f_h) : moduleData->getRamoPotentialMap(layer).getBinZ(1e3*depth_f_h); + + const bool isFirstZ_e = m_doInterpolateEfield ? m_ramoPotentialMap[layer].isFirstZ(1e3*depth_f_e) : moduleData->getRamoPotentialMap(layer).isFirstZ(1e3*depth_f_e); + const bool isOverflowZ_h = m_doInterpolateEfield ? m_ramoPotentialMap[layer].isOverflowZ(1e3*depth_f_h) : moduleData->getRamoPotentialMap(layer).isOverflowZ(1e3*depth_f_h); //Loop over nearest neighbours in x and y //We assume that the lateral diffusion is minimal @@ -472,74 +397,31 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit, double ramo_f_e = 0.0; double ramo_f_h = 0.0; - int nbin_ramo_f_e_x = int( 1 + std::abs(dPhi_f_e) * m_ramo_x_binMap ); - int nbin_ramo_f_e_y = int( 1 + std::abs(dEta_f_e) * m_ramo_y_binMap ); - - // Check for overflow in ramo hists in x- and y-direction - if (nbin_ramo_f_e_x > numBins_weightingPotential_x) { - nbin_ramo_f_e_x = numBins_weightingPotential_x + 1; - } - if (nbin_ramo_f_e_y > numBins_weightingPotential_y) { - nbin_ramo_f_e_y = numBins_weightingPotential_y + 1; - } - - if (nbin_ramo_f_e_x <= numBins_weightingPotential_x && nbin_ramo_f_e_y <= numBins_weightingPotential_y && - nbin_ramo_f_e_z <= numBins_weightingPotential_z) { + if (!isFirstZ_e) { if (m_doInterpolateEfield) { - ramo_f_e = m_ramoPotentialMap[layer]->GetBinContent(nbin_ramo_f_e_x, nbin_ramo_f_e_y, nbin_ramo_f_e_z); + ramo_f_e = m_ramoPotentialMap[layer].getContent(m_ramoPotentialMap[layer].getBinX(1e3*std::abs(dPhi_f_e)), m_ramoPotentialMap[layer].getBinY(1e3*std::abs(dEta_f_e)), ramo_f_e_bin_z); } else { - ramo_f_e = moduleData->getRamoPotentialMap(layer)->GetBinContent(nbin_ramo_f_e_x, nbin_ramo_f_e_y, - nbin_ramo_f_e_z); + ramo_f_e = moduleData->getRamoPotentialMap(layer).getContent(moduleData->getRamoPotentialMap(layer).getBinX(1e3*std::abs(dPhi_f_e)), moduleData->getRamoPotentialMap(layer).getBinY(1e3*std::abs(dEta_f_e)), ramo_f_e_bin_z); } } - int nbin_ramo_f_h_x = int( 1 + std::abs(dPhi_f_h) * m_ramo_x_binMap ); - int nbin_ramo_f_h_y = int( 1 + std::abs(dEta_f_h) * m_ramo_y_binMap ); - - // Check for overflow in ramo hists in x- and y-direction - if (nbin_ramo_f_h_x > numBins_weightingPotential_x) { - nbin_ramo_f_h_x = numBins_weightingPotential_x + 1; - } - if (nbin_ramo_f_h_y > numBins_weightingPotential_y) { - nbin_ramo_f_h_y = numBins_weightingPotential_y + 1; - } - - if (nbin_ramo_f_h_x <= numBins_weightingPotential_x && nbin_ramo_f_h_y <= numBins_weightingPotential_y && - nbin_ramo_f_h_z <= numBins_weightingPotential_z) { + if (!isOverflowZ_h) { if (m_doInterpolateEfield) { - ramo_f_h = m_ramoPotentialMap[layer]->GetBinContent(nbin_ramo_f_h_x, nbin_ramo_f_h_y, nbin_ramo_f_h_z); + ramo_f_h = m_ramoPotentialMap[layer].getContent(m_ramoPotentialMap[layer].getBinX(1e3*std::abs(dPhi_f_h)), m_ramoPotentialMap[layer].getBinY(1e3*std::abs(dEta_f_h)), ramo_f_h_bin_z); } else { - ramo_f_h = moduleData->getRamoPotentialMap(layer)->GetBinContent(nbin_ramo_f_h_x, nbin_ramo_f_h_y, - nbin_ramo_f_h_z); + ramo_f_h = moduleData->getRamoPotentialMap(layer).getContent(moduleData->getRamoPotentialMap(layer).getBinX(1e3*std::abs(dPhi_f_h)), moduleData->getRamoPotentialMap(layer).getBinY(1e3*std::abs(dEta_f_h)), ramo_f_h_bin_z); } } //Account for the imperfect binning that would cause charge to be double-counted + if (isOverflowZ_h) { + ramo_f_h = 0; + } - if (m_doInterpolateEfield) { - if (m_ramoPotentialMap[layer]->GetZaxis()->FindBin(depth_f_h * 1000) == - m_ramoPotentialMap[layer]->GetNbinsZ() + 1) { - ramo_f_h = 0.0; - } //this means the hole has reached the back end - - if (m_ramoPotentialMap[layer]->GetZaxis()->FindBin(depth_f_e * 1000) == 1) { - if (fabs(dEta_f_e) >= Module.etaPitch() / 2.0 || fabs(dPhi_f_e) >= Module.phiPitch() / 2.0) { - ramo_f_e = 0.0; - } else if (fabs(dEta_f_e) < Module.etaPitch() / 2.0 && fabs(dPhi_f_e) < Module.phiPitch() / 2.0) { - ramo_f_e = 1.0; - } - } - } else { - if (moduleData->getRamoPotentialMap(layer)->GetZaxis()->FindBin(depth_f_h * 1000) - == moduleData->getRamoPotentialMap(layer)->GetNbinsZ() + 1) { - ramo_f_h = 0; - } //this means the hole has reached the back end - - if (moduleData->getRamoPotentialMap(layer)->GetZaxis()->FindBin(depth_f_e * 1000) == 1) { - if (fabs(dEta_f_e) >= Module.etaPitch() / 2.0 || fabs(dPhi_f_e) >= Module.phiPitch() / 2.0) { - ramo_f_e = 0.0; - } else if (fabs(dEta_f_e) < Module.etaPitch() / 2.0 && fabs(dPhi_f_e) < Module.phiPitch() / 2.0) { - ramo_f_e = 1.0; - } + if (isFirstZ_e) { + if (fabs(dEta_f_e) >= Module.etaPitch() / 2.0 || fabs(dPhi_f_e) >= Module.phiPitch() / 2.0) { + ramo_f_e = 0.0; + } else if (fabs(dEta_f_e) < Module.etaPitch() / 2.0 && fabs(dPhi_f_e) < Module.phiPitch() / 2.0) { + ramo_f_e = 1.0; } } @@ -561,6 +443,7 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit, if (diode.isValid()) { chargedDiodes.add(diode, charge); } //IF + } //For q } //for p } else { //If no radDamage, run original diff --git a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.h b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.h index 085e28f748554d1e9625c5cdb10f833300c8188d..8dd2757762e66cd44bb101716e65329d98c76e0f 100644 --- a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.h +++ b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.h @@ -13,9 +13,10 @@ #define PIXELDIGITIZATION_SensorSimPlanarTool_H #include "AthenaBaseComps/AthAlgTool.h" -#include "SensorSimTool.h" #include "InDetCondTools/ISiLorentzAngleTool.h" +#include "SensorSimTool.h" #include "RadDamageUtil.h" +#include "PixelConditionsData/PixelHistoConverter.h" class SensorSimPlanarTool: public SensorSimTool { @@ -41,16 +42,11 @@ private: SensorSimPlanarTool(); // Map for radiation damage simulation - std::vector<TH3F*> m_ramoPotentialMap; - std::vector<TH2F*> m_distanceMap_e; - std::vector<TH2F*> m_distanceMap_h; - std::vector<TH2F*> m_lorentzMap_e; - std::vector<TH2F*> m_lorentzMap_h; - - // maps to directly get factor to calculate bin instead of calling FindBin - double m_ramo_x_binMap; - double m_ramo_y_binMap; - double m_ramo_z_binMap; + std::vector<PixelHistoConverter> m_ramoPotentialMap; + std::vector<PixelHistoConverter> m_distanceMap_e; + std::vector<PixelHistoConverter> m_distanceMap_h; + std::vector<PixelHistoConverter> m_lorentzMap_e; + std::vector<PixelHistoConverter> m_lorentzMap_h; Gaudi::Property<int> m_numberOfSteps {