From 4ebdf242e6cc0d9b6cf964a9ed2697c4dc072acf Mon Sep 17 00:00:00 2001
From: Tomas Dado <tomas.dado@cern.ch>
Date: Mon, 8 Feb 2021 17:46:52 +0100
Subject: [PATCH] First version
---
.../src/PixelConfigCondAlg.cxx | 44 +++-
.../PixelConditionsData/PixelHistoConverter.h | 64 ++++++
.../PixelConditionsData/PixelModuleData.h | 49 ++--
.../src/PixelHistoConverter.cxx | 182 +++++++++++++++
.../src/PixelModuleData.cxx | 38 +--
.../PixelDigitization/src/SensorSim3DTool.cxx | 143 ++++--------
.../PixelDigitization/src/SensorSim3DTool.h | 29 ++-
.../src/SensorSimPlanarTool.cxx | 216 ++++++++++--------
.../src/SensorSimPlanarTool.h | 18 +-
9 files changed, 524 insertions(+), 259 deletions(-)
create mode 100644 InnerDetector/InDetConditions/PixelConditionsData/PixelConditionsData/PixelHistoConverter.h
create mode 100644 InnerDetector/InDetConditions/PixelConditionsData/src/PixelHistoConverter.cxx
diff --git a/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/PixelConfigCondAlg.cxx b/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/PixelConfigCondAlg.cxx
index 7024b858a1a..8d266756d50 100644
--- a/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/PixelConfigCondAlg.cxx
+++ b/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/PixelConfigCondAlg.cxx
@@ -458,12 +458,17 @@ 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<TH3F*> ramoPotentialMapOld;
+ std::vector<TH2F*> lorentzMap_eOld;
+ std::vector<TH2F*> lorentzMap_hOld;
+ std::vector<TH2F*> distanceMap_eOld;
+ std::vector<TH2F*> distanceMap_hOld;
+
+ 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,13 +481,28 @@ 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"));
+ ramoPotentialMapOld.push_back(ramoPotentialMap_hold);
+
+ lorentzMap_eOld.push_back((TH2F*)mapsFile->Get("lorentz_map_e"));
+ lorentzMap_hOld.push_back((TH2F*)mapsFile->Get("lorentz_map_h"));
+ distanceMap_eOld.push_back((TH2F*)mapsFile->Get("edistance"));
+ distanceMap_hOld.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_eOld(lorentzMap_eOld);
+ writeCdo -> setLorentzMap_hOld(lorentzMap_hOld);
+ writeCdo -> setDistanceMap_eOld(distanceMap_eOld);
+ writeCdo -> setDistanceMap_hOld(distanceMap_hOld);
+ writeCdo -> setRamoPotentialMapOld(ramoPotentialMapOld);
writeCdo -> setLorentzMap_e(lorentzMap_e);
writeCdo -> setLorentzMap_h(lorentzMap_h);
writeCdo -> setDistanceMap_e(distanceMap_e);
diff --git a/InnerDetector/InDetConditions/PixelConditionsData/PixelConditionsData/PixelHistoConverter.h b/InnerDetector/InDetConditions/PixelConditionsData/PixelConditionsData/PixelHistoConverter.h
new file mode 100644
index 00000000000..484758f9d10
--- /dev/null
+++ b/InnerDetector/InDetConditions/PixelConditionsData/PixelConditionsData/PixelHistoConverter.h
@@ -0,0 +1,64 @@
+/*
+ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+ */
+/**
+ * @file PixelDigitization/PixelHistoConverter.h
+ * @author Tomas Dado <tomas.dado@cern.ch>
+ * @date February, 2020
+ * @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);
+
+ float GetContent(const std::size_t& x) const;
+ float GetContent(const std::size_t& x, const std::size_t& y) const;
+ float GetContent(const std::size_t& x, const std::size_t& y, const std::size_t& z) const;
+
+ bool IsOverflowZ(const float value) const;
+ 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<std::vector<std::vector<float> > > m_content;
+
+ bool SetAxis(Axis& axis, const TAxis* rootAxis);
+
+ std::size_t FindBin(const Axis& axis, const float value) const;
+
+};
+
+#endif
diff --git a/InnerDetector/InDetConditions/PixelConditionsData/PixelConditionsData/PixelModuleData.h b/InnerDetector/InDetConditions/PixelConditionsData/PixelConditionsData/PixelModuleData.h
index e1efb342b00..f612e1bc479 100644
--- a/InnerDetector/InDetConditions/PixelConditionsData/PixelConditionsData/PixelModuleData.h
+++ b/InnerDetector/InDetConditions/PixelConditionsData/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,31 @@ 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_eOld(std::vector<TH2F*> lorentzMap_e);
+ void setLorentzMap_hOld(std::vector<TH2F*> lorentzMap_h);
+ TH2F* getLorentzMap_eOld(int layer) const;
+ TH2F* getLorentzMap_hOld(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 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 setRamoPotentialMap(std::vector<TH3F*> ramoPotentialMap);
- TH3F* getRamoPotentialMap(int layer) const;
+ void setDistanceMap_eOld(std::vector<TH2F*> distanceMap_e);
+ void setDistanceMap_hOld(std::vector<TH2F*> distanceMap_h);
+ TH2F* getDistanceMap_eOld(int layer) const;
+ TH2F* getDistanceMap_hOld(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 setRamoPotentialMapOld(std::vector<TH3F*> ramoPotentialMap);
+ TH3F* getRamoPotentialMapOld(int layer) const;
+
+ void setRamoPotentialMap(std::vector<PixelHistoConverter> ramoPotentialMap);
+ const PixelHistoConverter& getRamoPotentialMap(int layer) const;
// Distortion parameters
void setDistortionInputSource(int distortionInputSource);
@@ -302,11 +316,16 @@ 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<TH2F*> m_lorentzMap_eOld;
+ std::vector<TH2F*> m_lorentzMap_hOld;
+ std::vector<TH2F*> m_distanceMap_eOld;
+ std::vector<TH2F*> m_distanceMap_hOld;
+ std::vector<TH3F*> m_ramoPotentialMapOld;
+ 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 00000000000..ad638026388
--- /dev/null
+++ b/InnerDetector/InDetConditions/PixelConditionsData/src/PixelHistoConverter.cxx
@@ -0,0 +1,182 @@
+/*
+ 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>
+
+//HistoConverter::HistoConverter() :
+//AthMessaging(Gaudi::svcLocator()->service< IMessageSvc >("MessageSvc"), "HistoConverter")
+PixelHistoConverter::PixelHistoConverter()
+{
+}
+
+StatusCode PixelHistoConverter::SetHisto1D(const TH1* histo) {
+
+ if (!histo) {
+ //ATH_MSG_ERROR("Nullptr passed");
+ return StatusCode::FAILURE;
+ }
+
+ if (!SetAxis(m_xAxis, histo->GetXaxis())) {
+ //ATH_MSG_ERROR("Failed to set the x axis");
+ return StatusCode::FAILURE;
+ }
+
+ /// fill the content
+ /// Add underflow and overflow bins
+ m_content.resize(m_xAxis.nBins + 2);
+ for (std::size_t x = 0; x <= m_xAxis.nBins + 1; ++x) {
+ m_content.at(x).resize(1);
+ m_content.at(x).at(0).resize(1);
+ m_content.at(x).at(0).at(0) = histo->GetBinContent(x);
+ }
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode PixelHistoConverter::SetHisto2D(const TH2* histo) {
+ if (!histo) {
+ //ATH_MSG_ERROR("Nullptr passed");
+ return StatusCode::FAILURE;
+ }
+
+ if (!SetAxis(m_xAxis, histo->GetXaxis())) {
+ //ATH_MSG_ERROR("Failed to set the x axis");
+ return StatusCode::FAILURE;
+ }
+ if (!SetAxis(m_yAxis, histo->GetYaxis())) {
+ //ATH_MSG_ERROR("Failed to set the y axis");
+ return StatusCode::FAILURE;
+ }
+
+ /// fill the content
+ /// Add underflow and overflow bins
+ m_content.resize(m_xAxis.nBins + 2);
+ for (std::size_t x = 0; x <= m_xAxis.nBins + 1; ++x) {
+ m_content.at(x).resize(m_yAxis.nBins + 2);
+ for (std::size_t y = 0; y <= m_yAxis.nBins + 1; ++y) {
+ m_content.at(x).at(y).resize(1);
+ m_content.at(x).at(y).at(0) = histo->GetBinContent(x,y);
+ }
+ }
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode PixelHistoConverter::SetHisto3D(const TH3* histo) {
+ if (!histo) {
+ //ATH_MSG_ERROR("Nullptr passed");
+ return StatusCode::FAILURE;
+ }
+
+ if (!SetAxis(m_xAxis, histo->GetXaxis())) {
+ //ATH_MSG_ERROR("Failed to set the x axis");
+ return StatusCode::FAILURE;
+ }
+ if (!SetAxis(m_yAxis, histo->GetYaxis())) {
+ //ATH_MSG_ERROR("Failed to set the y axis");
+ return StatusCode::FAILURE;
+ }
+ if (!SetAxis(m_zAxis, histo->GetZaxis())) {
+ //ATH_MSG_ERROR("Failed to set the Z axis");
+ return StatusCode::FAILURE;
+ }
+
+ /// fill the content
+ /// Add underflow and overflow bins
+ m_content.resize(m_xAxis.nBins + 2);
+ for (std::size_t x = 0; x <= m_xAxis.nBins + 1; ++x) {
+ m_content.at(x).resize(m_yAxis.nBins + 2);
+ for (std::size_t y = 0; y <= m_yAxis.nBins + 1; ++y) {
+ m_content.at(x).at(y).resize(m_zAxis.nBins + 2);
+ for (std::size_t z = 0; z <= m_zAxis.nBins + 1; ++z) {
+ m_content.at(x).at(y).at(z) = histo->GetBinContent(x,y,z);
+ }
+ }
+ }
+
+ return StatusCode::SUCCESS;
+}
+
+float PixelHistoConverter::GetContent(const std::size_t& x) const {
+
+ return m_content[x][0][0];
+}
+
+float PixelHistoConverter::GetContent(const std::size_t& x, const std::size_t& y) const {
+ return m_content[x][y][0];
+}
+
+float PixelHistoConverter::GetContent(const std::size_t& x, const std::size_t& y, const std::size_t& z) const {
+ return m_content[x][y][z];
+}
+
+bool PixelHistoConverter::IsOverflowZ(const float value) const {
+ return (value > m_zAxis.max) ? true : false;
+}
+
+bool PixelHistoConverter::IsFirstZ(const float value) const {
+ return (GetBinZ(value) == 1);
+}
+
+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) {
+ //ATH_MSG_ERROR("Nullptr passed");
+ return false;
+ }
+
+ axis.nBins = rootAxis->GetNbins();
+ axis.min = rootAxis->GetXmin();
+ axis.max = rootAxis->GetXmax();
+
+ if (axis.nBins < 1) {
+ //ATH_MSG_ERROR("Axis has less than 1 bin");
+ 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) {
+ //ATH_MSG_ERROR("Histogram does not have equidistant binning");
+ std::cerr << "Histogram does not have equidistant binning\n";
+
+ return false;
+ }
+ }
+
+ axis.width = 1./width;
+
+ return true;
+}
+
+std::size_t PixelHistoConverter::FindBin(const Axis& axis, const float value) const {
+ if (value < axis.min) return 0;
+ if (value > axis.max) return axis.nBins; // this is weird that it is not +1, but thats what we want
+
+ return ((value - axis.min) * axis.width) + 1;
+}
diff --git a/InnerDetector/InDetConditions/PixelConditionsData/src/PixelModuleData.cxx b/InnerDetector/InDetConditions/PixelConditionsData/src/PixelModuleData.cxx
index 16da803924f..3da1ca7d455 100644
--- a/InnerDetector/InDetConditions/PixelConditionsData/src/PixelModuleData.cxx
+++ b/InnerDetector/InDetConditions/PixelConditionsData/src/PixelModuleData.cxx
@@ -319,19 +319,31 @@ 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::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::setLorentzMap_eOld(std::vector<TH2F*> lorentzMap_e) { m_lorentzMap_eOld = lorentzMap_e; }
+void PixelModuleData::setLorentzMap_hOld(std::vector<TH2F*> lorentzMap_h) { m_lorentzMap_hOld = lorentzMap_h; }
+TH2F* PixelModuleData::getLorentzMap_eOld(int layer) const { return m_lorentzMap_eOld.at(layer); }
+TH2F* PixelModuleData::getLorentzMap_hOld(int layer) const { return m_lorentzMap_hOld.at(layer); }
+
+void PixelModuleData::setDistanceMap_eOld(std::vector<TH2F*> distanceMap_e) { m_distanceMap_eOld = distanceMap_e; }
+void PixelModuleData::setDistanceMap_hOld(std::vector<TH2F*> distanceMap_h) { m_distanceMap_hOld = distanceMap_h; }
+TH2F* PixelModuleData::getDistanceMap_eOld(int layer) const { return m_distanceMap_eOld.at(layer); }
+TH2F* PixelModuleData::getDistanceMap_hOld(int layer) const { return m_distanceMap_hOld.at(layer); }
+
+void PixelModuleData::setRamoPotentialMapOld(std::vector<TH3F*> ramoPotentialMap) { m_ramoPotentialMapOld = ramoPotentialMap; }
+TH3F* PixelModuleData::getRamoPotentialMapOld(int layer) const { return m_ramoPotentialMapOld.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<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 29a1fb0afb6..c499c033c2b 100644
--- a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.cxx
+++ b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.cxx
@@ -90,23 +90,17 @@ StatusCode SensorSim3DTool::initialize() {
//std::unique_ptr<TFile> mapsFile=std::make_unique<TFile>( (mapsPath_list.at(i)).c_str() ); //this is the ramo potential.
TFile* mapsFile = new TFile((mapsPath_list.at(i)).c_str()); //this is the ramo potential.
- 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 iayer.first=1
- Layer.second = i; //Layer: 0 = IBL Planar, 1=B-Layer, 2=Layer-1, 3=Layer-2
- //For 3D sensor, only possible index should be 0-0
- //May want to be changed in the future, since there's no point having an index with only one possible value
-
//Setup ramo weighting field map
- TH3F* ramoPotentialMap_hold;
+ TH2F* ramoPotentialMap_hold;
ramoPotentialMap_hold = 0;
- ramoPotentialMap_hold = (TH3F* ) mapsFile->Get("ramo");
+ ramoPotentialMap_hold = (TH2F* ) mapsFile->Get("ramo");
if (ramoPotentialMap_hold == 0) {
ATH_MSG_INFO("Did not find a Ramo potential map. Will use an approximate form.");
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().SetHisto2D(ramoPotentialMap_hold));
ATH_MSG_INFO("Using fluence " << m_fluence_layers.at(i));
@@ -118,7 +112,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 +135,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;
@@ -355,14 +346,8 @@ StatusCode SensorSim3DTool::induceCharge(const TimedHitPtr<SiHit>& phit, SiCharg
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));
-
+ double 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);
double xposFinal = getTrappingPositionY(yposDiff, xposDiff, std::min(driftTime, timeToElectrode), isHole);
@@ -375,6 +360,9 @@ StatusCode SensorSim3DTool::induceCharge(const TimedHitPtr<SiHit>& phit, SiCharg
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;
@@ -387,10 +375,8 @@ StatusCode SensorSim3DTool::induceCharge(const TimedHitPtr<SiHit>& phit, SiCharg
//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);
+ double ramoInit = m_ramoPotentialMap[index].GetContent(m_ramoPotentialMap[index].GetBinX(1000*(y_pix + 0.5*pixel_size_y - yNeighbor)), ramo_init_bin_y);
+ double 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);
@@ -504,8 +490,6 @@ StatusCode SensorSim3DTool::induceCharge(const TimedHitPtr<SiHit>& phit, SiCharg
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);
@@ -527,7 +511,7 @@ StatusCode SensorSim3DTool::induceCharge(const TimedHitPtr<SiHit>& phit, SiCharg
}
// 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);
@@ -563,7 +547,7 @@ 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) {
if (readout == "FEI4") {
for (std::multimap<std::pair<int, int>, double >::iterator it = m_probMapFEI4.begin(); it != m_probMapFEI4.end();
++it) {
@@ -586,7 +570,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") {
@@ -603,15 +587,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)
}
@@ -652,79 +629,41 @@ 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
+ std::size_t index = 0;
// Uses (x,y) position in um to return time to electrode in ns
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);
+ timeToElectrode = m_timeMap_e[index].GetContent(m_timeMap_e[index].GetBinX(1e3*x), m_timeMap_e[index].GetBinY(1e3*y));
}
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_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 1d4466900a0..e01087010f7 100644
--- a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.h
+++ b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.h
@@ -16,6 +16,7 @@
#include "GaudiKernel/ToolHandle.h"
#include "RadDamageUtil.h"
+#include "PixelConditionsData/PixelHistoConverter.h"
class SensorSim3DTool: public SensorSimTool {
public:
@@ -31,9 +32,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&);
+ double getProbMapEntry(const std::string&, int, int) const;
double getElectricField(double x, double y);
double getMobility(double electricField, bool isHoleBit);
@@ -41,7 +42,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();
@@ -50,19 +50,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 a5edb25dd32..3679cd6176f 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;
}
@@ -256,6 +259,14 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit, SiC
int numBins_weightingPotential_y = 0;
int numBins_weightingPotential_z = 0;
+
+ numBins_driftTime_e = moduleData->getDistanceMap_eOld(layer)->GetNbinsY();
+ numBins_driftTime_h = moduleData->getDistanceMap_hOld(layer)->GetNbinsY();
+ numBins_weightingPotential_x = moduleData->getRamoPotentialMapOld(layer)->GetNbinsX();
+ numBins_weightingPotential_y = moduleData->getRamoPotentialMapOld(layer)->GetNbinsY();
+ numBins_weightingPotential_z = moduleData->getRamoPotentialMapOld(layer)->GetNbinsZ();
+
+
if (m_doRadDamage && !(Module.isDBM()) && Module.isBarrel()) {
centreOfPixel_i = p_design.positionFromColumnRow(pixel_i.etaIndex(), pixel_i.phiIndex());
@@ -265,21 +276,6 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit, SiC
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
@@ -295,32 +291,28 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit, SiC
nontrappingProbability = exp(-dist_electrode / collectionDist);
}
+ double m_ramo_x_binMap;
+ double m_ramo_y_binMap;
+ double m_ramo_z_binMap;
+ const std::size_t depth_f_e_bin_x = m_doInterpolateEfield ? m_distanceMap_e[layer].GetBinX(dist_electrode) : moduleData->getDistanceMap_e(layer).GetBinX(dist_electrode);
+ const std::size_t depth_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);
+
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()));
+ (moduleData->getRamoPotentialMapOld(layer)->GetNbinsX() /
+ (moduleData->getRamoPotentialMapOld(layer)->GetXaxis()->GetXmax() -
+ moduleData->getRamoPotentialMapOld(layer)->GetXaxis()->GetXmin()));
m_ramo_y_binMap = 1000. *
- (moduleData->getRamoPotentialMap(layer)->GetNbinsY() /
- (moduleData->getRamoPotentialMap(layer)->GetYaxis()->GetXmax() -
- moduleData->getRamoPotentialMap(layer)->GetYaxis()->GetXmin()));
+ (moduleData->getRamoPotentialMapOld(layer)->GetNbinsY() /
+ (moduleData->getRamoPotentialMapOld(layer)->GetYaxis()->GetXmax() -
+ moduleData->getRamoPotentialMapOld(layer)->GetYaxis()->GetXmin()));
m_ramo_z_binMap = 1000. *
- (moduleData->getRamoPotentialMap(layer)->GetNbinsZ() /
- (moduleData->getRamoPotentialMap(layer)->GetZaxis()->GetXmax() -
- moduleData->getRamoPotentialMap(layer)->GetZaxis()->GetXmin()));
+ (moduleData->getRamoPotentialMapOld(layer)->GetNbinsZ() /
+ (moduleData->getRamoPotentialMapOld(layer)->GetZaxis()->GetXmax() -
+ moduleData->getRamoPotentialMapOld(layer)->GetZaxis()->GetXmin()));
}
for (int j = 0; j < ncharges; j++) {
@@ -335,49 +327,55 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit, SiC
double depth_f_h = 0.0;
double tanLorentz_e = 0.0;
double tanLorentz_h = 0.0;
+ double depth_f_eOld = 0.0;
+ double depth_f_hOld = 0.0;
+ double tanLorentz_eOld = 0.0;
+ double tanLorentz_hOld = 0.0;
//TODO: the holes map does not currently extend for a drift time long enough that, any hole will reach
//the corresponding electrode. This needs to be rectified by either (a) extrapolating the current map or
//(b) making a new map with a y-axis (drift time) that extends to at least 18 ns so all charge carriers reach
// 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(depth_f_e_bin_x, m_distanceMap_e[layer].GetBinY(drifttime_e));
+ depth_f_h = m_distanceMap_h[layer].GetContent(depth_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);
+ int nbin_z_e_xbin = moduleData->getDistanceMap_eOld(layer)->GetXaxis()->FindBin(dist_electrode);
+ int nbin_z_e_ybin = moduleData->getDistanceMap_eOld(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);
+ depth_f_eOld = moduleData->getDistanceMap_eOld(layer)->GetBinContent(nbin_z_e_xbin, nbin_z_e_ybin);
+ int nbin_z_h_xbin = moduleData->getDistanceMap_hOld(layer)->GetXaxis()->FindBin(dist_electrode);
+ int nbin_z_h_ybin = moduleData->getDistanceMap_hOld(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_hOld = moduleData->getDistanceMap_hOld(layer)->GetBinContent(nbin_z_h_xbin, nbin_z_h_ybin);
+ int nbin_Lorentz_e = moduleData->getLorentzMap_eOld(layer)->FindBin(dist_electrode, depth_f_eOld);
+ tanLorentz_eOld = moduleData->getLorentzMap_eOld(layer)->GetBinContent(nbin_Lorentz_e);
+ int nbin_Lorentz_h = moduleData->getLorentzMap_hOld(layer)->FindBin(dist_electrode, depth_f_hOld);
+ tanLorentz_hOld = moduleData->getLorentzMap_hOld(layer)->GetBinContent(nbin_Lorentz_h);
+ depth_f_e = moduleData->getDistanceMap_e(layer).GetContent(depth_f_e_bin_x, moduleData->getDistanceMap_e(layer).GetBinY(drifttime_e));
+ depth_f_h = moduleData->getDistanceMap_h(layer).GetContent(depth_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));
+ if (depth_f_e != depth_f_eOld || depth_f_h != depth_f_hOld) {
+ std::cout << "Depth E old: " << depth_f_eOld << ", new: " << depth_f_e << ", depth H old: " << depth_f_hOld << ", new: " << depth_f_h << std::endl;;
+ std::cout << "old depth E bin X: " << nbin_z_e_xbin << ", new: " << depth_f_e_bin_x << std::endl;
+ std::cout << "old depth E bin Y: " << nbin_z_e_ybin << ", new: " << moduleData->getDistanceMap_e(layer).GetBinY(drifttime_e) << std::endl;
+ std::cout << "old depth H bin X: " << nbin_z_h_xbin << ", new: " << depth_f_h_bin_x << std::endl;
+ std::cout << "old depth H bin Y: " << nbin_z_h_ybin << ", new: " << moduleData->getDistanceMap_h(layer).GetBinY(drifttime_h) << std::endl;
+ }
}
double dz_e = fabs(dist_electrode - depth_f_e);
double dz_h = fabs(depth_f_h - dist_electrode);
double coLorentz_e = std::sqrt(1.0 + std::pow(tanLorentz_e, 2));
+ double dz_eOld = fabs(dist_electrode - depth_f_eOld);
+ double dz_hOld = fabs(depth_f_hOld - dist_electrode);
+ double coLorentz_eOld = std::sqrt(1.0 + std::pow(tanLorentz_eOld, 2));
//Apply drift due to Lorentz force and diffusion
double phiRand = CLHEP::RandGaussZiggurat::shoot(rndmEngine);
@@ -394,7 +392,18 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit, SiC
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;
-
+
+ double coLorentz_hOld = std::sqrt(1.0 + std::pow(tanLorentz_hOld, 2));
+ double rdif_hOld = this->m_diffusionConstant * sqrt(fabs(dist_electrode - depth_f_hOld) * coLorentz_hOld / 0.3);
+ double phi_f_hOld = phi_i + dz_hOld * tanLorentz_hOld + rdif_hOld * phiRand;
+ double eta_f_hOld = eta_i + rdif_hOld * etaRand;
+
+ //if (depth_f_e != depth_f_eOld || depth_f_h != depth_f_hOld) {
+ // std::cout << "Depth E old: " << depth_f_eOld << ", new: " << depth_f_e << ", depth H old: " << depth_f_hOld << ", new: " << depth_f_h << "\n";
+ //}
+ //if (tanLorentz_e != tanLorentz_eOld || tanLorentz_h != tanLorentz_hOld) {
+ // std::cout << "TanLor E old: " << tanLorentz_eOld << ", new: " << tanLorentz_e << ", tanLor H old: " << tanLorentz_hOld << ", new: " << tanLorentz_e << "\n";
+ //}
// amount of energy to be converted into charges at current step
double energy_per_step = 1.0 * iHitRecord.second / 1.E+6 / ncharges;
@@ -409,8 +418,8 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit, SiC
// 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 );
+ nbin_ramo_f_e_z = int( minz + depth_f_eOld * m_ramo_z_binMap );
+ nbin_ramo_f_h_z = int( minz + depth_f_hOld * m_ramo_z_binMap );
// Check for overflow in ramo hists in z-direction
if (nbin_ramo_f_h_z > numBins_weightingPotential_z) {
@@ -419,6 +428,11 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit, SiC
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
@@ -465,12 +479,20 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit, SiC
double dPhi_f_h = pixelPhi_f_h - dPhi_nn_centre;
dEta_f_h *= scale_f;
+ double dEta_f_eOld = pixelEta_f_e - dEta_nn_centre;
+ double dPhi_f_eOld = pixelPhi_f_e - dPhi_nn_centre;
+ dEta_f_eOld *= scale_f;
+ double dEta_f_hOld = pixelEta_f_h - dEta_nn_centre;
+ double dPhi_f_hOld = pixelPhi_f_h - dPhi_nn_centre;
+ dEta_f_hOld *= scale_f;
//Boundary check on maps
double ramo_f_e = 0.0;
double ramo_f_h = 0.0;
+ double ramo_f_eOld = 0.0;
+ double ramo_f_hOld = 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 );
+ int nbin_ramo_f_e_x = int( 1 + std::abs(dPhi_f_eOld) * m_ramo_x_binMap );
+ int nbin_ramo_f_e_y = int( 1 + std::abs(dEta_f_eOld) * m_ramo_y_binMap );
// Check for overflow in ramo hists in x- and y-direction
if (nbin_ramo_f_e_x > numBins_weightingPotential_x) {
@@ -483,9 +505,10 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit, SiC
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 (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,
+ 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);
+ ramo_f_eOld = moduleData->getRamoPotentialMapOld(layer)->GetBinContent(nbin_ramo_f_e_x, nbin_ramo_f_e_y,
nbin_ramo_f_e_z);
}
}
@@ -504,34 +527,33 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit, SiC
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 (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,
+ 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);
+ ramo_f_hOld = moduleData->getRamoPotentialMapOld(layer)->GetBinContent(nbin_ramo_f_h_x, nbin_ramo_f_h_y,
nbin_ramo_f_h_z);
}
}
//Account for the imperfect binning that would cause charge to be double-counted
if (m_doInterpolateEfield) {
- if (m_ramoPotentialMap[layer]->GetZaxis()->FindBin(depth_f_h * 1000) ==
- m_ramoPotentialMap[layer]->GetNbinsZ() + 1) {
- ramo_f_h = 0.0;
+ } else {
+ if (moduleData->getRamoPotentialMapOld(layer)->GetZaxis()->FindBin(depth_f_h * 1000)
+ == moduleData->getRamoPotentialMapOld(layer)->GetNbinsZ() + 1) {
+ ramo_f_hOld = 0;
} //this means the hole has reached the back end
+ if (isOverflowZ_h) {
+ ramo_f_h = 0;
+ }
- if (m_ramoPotentialMap[layer]->GetZaxis()->FindBin(depth_f_e * 1000) == 1) {
+ if (moduleData->getRamoPotentialMapOld(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;
+ ramo_f_eOld = 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;
+ ramo_f_eOld = 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 (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) {
@@ -539,12 +561,21 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit, SiC
}
}
}
+
+ //if (ramo_f_e != ramo_f_eOld || ramo_f_h != ramo_f_hOld) {
+ // std::cout << "Ramo E old: " << ramo_f_eOld << ", new: " << ramo_f_e << ", Ramo H old: " << ramo_f_hOld << ", new: " << ramo_f_h << "\n";
+ //}
//Given final position of charge carrier, find induced charge. The difference in Ramo weighting potential
// gives the fraction of charge induced.
//The energy_per_step is transformed into charge with the eleholePair per Energy
+ double induced_chargeOld = (ramo_f_eOld - ramo_f_hOld) * energy_per_step * eleholePairEnergy;
double induced_charge = (ramo_f_e - ramo_f_h) * energy_per_step * eleholePairEnergy;
+ if (induced_chargeOld != induced_chargeOld) {
+ std::cout << "induced_chargeOld: " << induced_chargeOld << ", new: " << induced_charge << "\n";
+ }
+
//Collect charge in centre of each pixel, since location within pixel doesn't matter for record
SiLocalPosition chargePos = Module.hitLocalToLocal(centreOfPixel_nn.xEta(), centreOfPixel_nn.xPhi());
@@ -560,6 +591,9 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit, SiC
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 08aee096fd1..4a582b862f3 100644
--- a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.h
+++ b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.h
@@ -12,9 +12,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 {
public:
@@ -34,16 +35,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
{
--
GitLab