From 5baf95b9859e67be0295c409cfe8414349fedbbf Mon Sep 17 00:00:00 2001
From: Soshi Tsuno <soshi.tsuno@cern.ch>
Date: Mon, 5 Apr 2021 16:34:27 +0000
Subject: [PATCH] IBL 3D radiation map move into conditions data
(ATLASSIM-5148)
---
.../src/PixelConfigCondAlg.cxx | 136 ++++++++++-
.../src/PixelConfigCondAlg.h | 37 +++
.../PixelConditionsData/PixelModuleData.h | 49 +++-
.../src/PixelModuleData.cxx | 36 ++-
.../python/PixelDigitizationConfig.py | 2 -
.../python/PixelDigitizationConfigNew.py | 4 +-
.../PixelDigitization/src/SensorSim3DTool.cxx | 224 ++++--------------
.../PixelDigitization/src/SensorSim3DTool.h | 25 +-
8 files changed, 290 insertions(+), 223 deletions(-)
diff --git a/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/PixelConfigCondAlg.cxx b/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/PixelConfigCondAlg.cxx
index eb927777072..7d6faad1a2a 100644
--- a/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/PixelConfigCondAlg.cxx
+++ b/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/PixelConfigCondAlg.cxx
@@ -173,7 +173,7 @@ StatusCode PixelConfigCondAlg::execute(const EventContext& ctx) const {
// mapping files for radiation damage simulation
std::vector<std::string> mapsPath_list;
-
+ std::vector<std::string> mapsPath_list3D;
int currentRunNumber = ctx.eventID().run_number();
if (currentRunNumber<222222) {
@@ -214,6 +214,7 @@ StatusCode PixelConfigCondAlg::execute(const EventContext& ctx) const {
for (size_t i=0; i<m_BarrelFluenceMapRUN1.size(); i++) {
mapsPath_list.push_back(PathResolverFindCalibFile(m_BarrelFluenceMapRUN1[i]));
}
+
}
else if (currentRunNumber<240000) { // RUN2 (mc15)
writeCdo -> setBarrelToTThreshold(m_BarrelToTThreshold2016);
@@ -273,6 +274,13 @@ StatusCode PixelConfigCondAlg::execute(const EventContext& ctx) const {
for (size_t i=0; i<m_BarrelFluenceMap2016.size(); i++) {
mapsPath_list.push_back(PathResolverFindCalibFile(m_BarrelFluenceMap2016[i]));
}
+
+ // Radiation damage simulation for 3D sensor
+ writeCdo -> setFluenceLayer3D(m_3DFluence2016);
+ for (size_t i=0; i<m_3DFluenceMap2016.size(); i++) {
+ mapsPath_list3D.push_back(PathResolverFindCalibFile(m_3DFluenceMap2016[i]));
+ }
+
}
else if (currentRunNumber<250000) { // RUN4
writeCdo -> setBarrelToTThreshold(m_BarrelToTThresholdITK);
@@ -306,6 +314,12 @@ StatusCode PixelConfigCondAlg::execute(const EventContext& ctx) const {
mapsPath_list.push_back(PathResolverFindCalibFile(m_BarrelFluenceMapITK[i]));
}
+ // Radiation damage simulation for 3D sensor
+ writeCdo -> setFluenceLayer3D(m_3DFluenceITK);
+ for (size_t i=0; i<m_3DFluenceMapITK.size(); i++) {
+ mapsPath_list3D.push_back(PathResolverFindCalibFile(m_3DFluenceMapITK[i]));
+ }
+
}
else if (currentRunNumber<300000) { // RUN2 2015/2016 (mc16a)
writeCdo -> setBarrelToTThreshold(m_BarrelToTThreshold2016);
@@ -356,6 +370,13 @@ StatusCode PixelConfigCondAlg::execute(const EventContext& ctx) const {
for (size_t i=0; i<m_BarrelFluenceMap2016.size(); i++) {
mapsPath_list.push_back(PathResolverFindCalibFile(m_BarrelFluenceMap2016[i]));
}
+
+ // Radiation damage simulation for 3D sensor
+ writeCdo -> setFluenceLayer3D(m_3DFluence2016);
+ for (size_t i=0; i<m_3DFluenceMap2016.size(); i++) {
+ mapsPath_list3D.push_back(PathResolverFindCalibFile(m_3DFluenceMap2016[i]));
+ }
+
}
else if (currentRunNumber<310000) { // RUN2 2017 (mc16d)
writeCdo -> setBarrelToTThreshold(m_BarrelToTThreshold2017);
@@ -406,6 +427,13 @@ StatusCode PixelConfigCondAlg::execute(const EventContext& ctx) const {
for (size_t i=0; i<m_BarrelFluenceMap2017.size(); i++) {
mapsPath_list.push_back(PathResolverFindCalibFile(m_BarrelFluenceMap2017[i]));
}
+
+ // Radiation damage simulation for 3D sensor
+ writeCdo -> setFluenceLayer3D(m_3DFluence2017);
+ for (size_t i=0; i<m_3DFluenceMap2017.size(); i++) {
+ mapsPath_list3D.push_back(PathResolverFindCalibFile(m_3DFluenceMap2017[i]));
+ }
+
}
else { // RUN2 2018 (mc16e)
writeCdo -> setBarrelToTThreshold(m_BarrelToTThreshold2018);
@@ -456,8 +484,16 @@ StatusCode PixelConfigCondAlg::execute(const EventContext& ctx) const {
for (size_t i=0; i<m_BarrelFluenceMap2018.size(); i++) {
mapsPath_list.push_back(PathResolverFindCalibFile(m_BarrelFluenceMap2018[i]));
}
+
+ // Radiation damage simulation for 3D sensor
+ writeCdo -> setFluenceLayer3D(m_3DFluence2018);
+ for (size_t i=0; i<m_3DFluenceMap2018.size(); i++) {
+ mapsPath_list3D.push_back(PathResolverFindCalibFile(m_3DFluenceMap2018[i]));
+ }
+
}
+ // Create mapping file for radiation damage simulation
std::vector<PixelHistoConverter> ramoPotentialMap;
std::vector<PixelHistoConverter> lorentzMap_e;
std::vector<PixelHistoConverter> lorentzMap_h;
@@ -517,6 +553,104 @@ StatusCode PixelConfigCondAlg::execute(const EventContext& ctx) const {
writeCdo -> setDistanceMap_h(distanceMap_h);
writeCdo -> setRamoPotentialMap(ramoPotentialMap);
+ // Create mapping file for radiation damage simulation for 3D sensor
+ std::vector<PixelHistoConverter> ramoPotentialMap3D;
+ std::vector<PixelHistoConverter> eFieldMap3D;
+ std::vector<PixelHistoConverter> xPositionMap3D_e;
+ std::vector<PixelHistoConverter> xPositionMap3D_h;
+ std::vector<PixelHistoConverter> yPositionMap3D_e;
+ std::vector<PixelHistoConverter> yPositionMap3D_h;
+ std::vector<PixelHistoConverter> timeMap3D_e;
+ std::vector<PixelHistoConverter> timeMap3D_h;
+ PixelHistoConverter avgChargeMap3D_e;
+ PixelHistoConverter avgChargeMap3D_h;
+
+ for (unsigned int i=0; i<mapsPath_list3D.size(); i++) {
+ ATH_MSG_INFO("Using maps located in: "<<mapsPath_list3D.at(i) << " for 3D sensor layer No." << i);
+ std::unique_ptr<TFile> mapsFile3D(TFile::Open((mapsPath_list3D.at(i)).c_str(), "READ")); //this is the ramo potential
+
+ if (!mapsFile3D) {
+ ATH_MSG_FATAL("Cannot open file: " << mapsPath_list3D.at(i));
+ return StatusCode::FAILURE;
+ }
+
+ //Setup ramo weighting field map
+ std::unique_ptr<TH2F> ramoPotentialMap3D_hold(mapsFile3D->Get<TH2F>("ramo"));
+ std::unique_ptr<TH2F> eFieldMap3D_hold(mapsFile3D->Get<TH2F>("efield"));
+ if (!ramoPotentialMap3D_hold || !eFieldMap3D_hold) {
+ ATH_MSG_FATAL("Did not find a Ramo potential or e-field map for 3D and an approximate form is available yet. Exit...");
+ return StatusCode::FAILURE;
+ }
+ ramoPotentialMap3D_hold->SetDirectory(nullptr);
+ eFieldMap3D_hold->SetDirectory(nullptr);
+ ramoPotentialMap3D.emplace_back();
+ eFieldMap3D.emplace_back();
+ ATH_CHECK(ramoPotentialMap3D.back().setHisto2D(ramoPotentialMap3D_hold.get()));
+ ATH_CHECK(eFieldMap3D.back().setHisto2D(eFieldMap3D_hold.get()));
+
+ //Now setup the E-field.
+ std::unique_ptr<TH3F> xPositionMap3D_e_hold(mapsFile3D->Get<TH3F>("xPosition_e"));
+ std::unique_ptr<TH3F> xPositionMap3D_h_hold(mapsFile3D->Get<TH3F>("xPosition_h"));
+ std::unique_ptr<TH3F> yPositionMap3D_e_hold(mapsFile3D->Get<TH3F>("yPosition_e"));
+ std::unique_ptr<TH3F> yPositionMap3D_h_hold(mapsFile3D->Get<TH3F>("yPosition_h"));
+ std::unique_ptr<TH2F> timeMap3D_e_hold(mapsFile3D->Get<TH2F>("etimes"));
+ std::unique_ptr<TH2F> timeMap3D_h_hold(mapsFile3D->Get<TH2F>("htimes"));
+
+ if (!xPositionMap3D_e_hold || !xPositionMap3D_h_hold || !yPositionMap3D_e_hold || !yPositionMap3D_h_hold || !timeMap3D_e_hold || !timeMap3D_h_hold) {
+ ATH_MSG_FATAL("Cannot find one of the maps.");
+ return StatusCode::FAILURE;
+ }
+
+ xPositionMap3D_e_hold->SetDirectory(nullptr);
+ xPositionMap3D_h_hold->SetDirectory(nullptr);
+ yPositionMap3D_e_hold->SetDirectory(nullptr);
+ yPositionMap3D_h_hold->SetDirectory(nullptr);
+ timeMap3D_e_hold->SetDirectory(nullptr);
+ timeMap3D_h_hold->SetDirectory(nullptr);
+
+ //Now, determine the time to reach the electrode and the trapping position.
+ xPositionMap3D_e.emplace_back();
+ xPositionMap3D_h.emplace_back();
+ yPositionMap3D_e.emplace_back();
+ yPositionMap3D_h.emplace_back();
+ timeMap3D_e.emplace_back();
+ timeMap3D_h.emplace_back();
+ ATH_CHECK(xPositionMap3D_e.back().setHisto3D(xPositionMap3D_e_hold.get()));
+ ATH_CHECK(xPositionMap3D_h.back().setHisto3D(xPositionMap3D_h_hold.get()));
+ ATH_CHECK(yPositionMap3D_e.back().setHisto3D(yPositionMap3D_e_hold.get()));
+ ATH_CHECK(yPositionMap3D_h.back().setHisto3D(yPositionMap3D_h_hold.get()));
+ ATH_CHECK(timeMap3D_e.back().setHisto2D(timeMap3D_e_hold.get()));
+ ATH_CHECK(timeMap3D_h.back().setHisto2D(timeMap3D_h_hold.get()));
+
+ std::unique_ptr<TH2F> avgCharge3D_e_hold(mapsFile3D->Get<TH2F>("avgCharge_e"));
+ std::unique_ptr<TH2F> avgCharge3D_h_hold(mapsFile3D->Get<TH2F>("avgCharge_h"));
+
+ if (!avgCharge3D_e_hold || !avgCharge3D_h_hold) {
+ ATH_MSG_ERROR("Cannot find one of the charge maps.");
+ return StatusCode::FAILURE;
+ }
+
+ avgCharge3D_e_hold->SetDirectory(nullptr);
+ avgCharge3D_h_hold->SetDirectory(nullptr);
+
+ // Get average charge data (for charge chunk effect correction)
+ ATH_CHECK(avgChargeMap3D_e.setHisto2D(avgCharge3D_e_hold.get()));
+ ATH_CHECK(avgChargeMap3D_h.setHisto2D(avgCharge3D_h_hold.get()));
+
+ mapsFile3D->Close();
+ }
+
+ writeCdo -> setRamoPotentialMap3D(ramoPotentialMap3D);
+ writeCdo -> setEFieldMap3D(eFieldMap3D);
+ writeCdo -> setXPositionMap3D_e(xPositionMap3D_e);
+ writeCdo -> setXPositionMap3D_h(xPositionMap3D_h);
+ writeCdo -> setYPositionMap3D_e(yPositionMap3D_e);
+ writeCdo -> setYPositionMap3D_h(yPositionMap3D_h);
+ writeCdo -> setTimeMap3D_e(timeMap3D_e);
+ writeCdo -> setTimeMap3D_h(timeMap3D_h);
+ writeCdo -> setAvgChargeMap3D_e(avgChargeMap3D_e);
+ writeCdo -> setAvgChargeMap3D_h(avgChargeMap3D_h);
+
//=======================
// Combine time interval
//=======================
diff --git a/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/PixelConfigCondAlg.h b/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/PixelConfigCondAlg.h
index 7dacac1a898..c68e11312ba 100644
--- a/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/PixelConfigCondAlg.h
+++ b/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/PixelConfigCondAlg.h
@@ -125,6 +125,9 @@ class PixelConfigCondAlg : public AthReentrantAlgorithm {
// DisalbePix: [N/A] [9e-3,9e-3,9e-3] [9e-3,9e-3,9e-3] [9e-3,9e-3,9e-3]
// BiasVoltage: [N/A] [ 500, 500, 500] [ 500, 500, 500] [ 500, 500, 500]
//
+ // IBL 3D:
+ // Fluence(e14): [N/A] [0.50] [0.50] [50.0]
+ //
// See https://twiki.cern.ch/twiki/bin/view/Atlas/PixelConditionsRUN2
// for further details.
//
@@ -153,6 +156,8 @@ class PixelConfigCondAlg : public AthReentrantAlgorithm {
// BiasVoltage: [ 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150]
// Fluence(e14):[ n/a, n/a, n/a, n/a, n/a, n/a, n/a, n/a, n/a, n/a, n/a, n/a, n/a, n/a]
//
+ // ITK 3D:
+ // Fluence(e14):[ n/a]
//====================================================================================
//====================================================================================
@@ -258,6 +263,14 @@ class PixelConfigCondAlg : public AthReentrantAlgorithm {
Gaudi::Property<std::vector<double>> m_DBMDisableProbability2016
{this, "DBMDisableProbability2016", {9e-3,9e-3,9e-3}, "Disable probability of DBM pixel layers in 2015/2016"};
+ // IBL 3D RUN2 2015/2016
+ Gaudi::Property<std::vector<double>> m_3DFluence2016
+ {this, "Barrel3DFluence2016", {5.0e15}, "Barrel3D fluence for radiation damage in 2016"};
+
+ Gaudi::Property<std::vector<std::string>> m_3DFluenceMap2016
+ {this, "Barrel3DFluenceMap2016", {"PixelDigitization/TCAD_IBL_3Dsensors_efields/phi_5e15_160V.root"},
+ "Barrel3D fluence map for radiation damage in 2016"};
+
//====================================================================================
// Barrel RUN2 2017
Gaudi::Property<std::vector<int>> m_BarrelToTThreshold2017
@@ -357,6 +370,14 @@ class PixelConfigCondAlg : public AthReentrantAlgorithm {
Gaudi::Property<std::vector<double>> m_DBMDisableProbability2017
{this, "DBMDisableProbability2017", {9e-3,9e-3,9e-3}, "Disable probability of DBM pixel layers in 2017"};
+ // IBL 3D RUN2 2017
+ Gaudi::Property<std::vector<double>> m_3DFluence2017
+ {this, "Barrel3DFluence2017", {0.50e14}, "Barrel3D fluence for radiation damage in 2017"};
+
+ Gaudi::Property<std::vector<std::string>> m_3DFluenceMap2017
+ {this, "Barrel3DFluenceMap2017", {"PixelDigitization/TCAD_IBL_3Dsensors_efields/phi_5e15_160V.root"},
+ "Barrel3D fluence map for radiation damage in 2017"};
+
//====================================================================================
// Barrel RUN2 2018
Gaudi::Property<std::vector<int>> m_BarrelToTThreshold2018
@@ -456,6 +477,14 @@ class PixelConfigCondAlg : public AthReentrantAlgorithm {
Gaudi::Property<std::vector<double>> m_DBMDisableProbability2018
{this, "DBMDisableProbability2018", {9e-3,9e-3,9e-3}, "Disable probability of DBM pixel layers in 2018"};
+ // IBL 3D RUN2 2018
+ Gaudi::Property<std::vector<double>> m_3DFluence2018
+ {this, "Barrel3DFluence2018", {5.0e15}, "Barrel3D fluence for radiation damage in 2018"};
+
+ Gaudi::Property<std::vector<std::string>> m_3DFluenceMap2018
+ {this, "Barrel3DFluenceMap2018", {"PixelDigitization/TCAD_IBL_3Dsensors_efields/phi_5e15_160V.root"},
+ "Barrel3D fluence map for radiation damage in 2018"};
+
//====================================================================================
// Barrel RUN1
Gaudi::Property<std::vector<int>> m_BarrelToTThresholdRUN1
@@ -595,6 +624,14 @@ class PixelConfigCondAlg : public AthReentrantAlgorithm {
Gaudi::Property<std::vector<double>> m_EndcapLorentzAngleCorrITK
{this, "EndcapLorentzAngleCorrITK", {1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0}, "Scale factor for Lorentz angle of endcap pixel layers in ITK"};
+ // ITK 3D
+ Gaudi::Property<std::vector<double>> m_3DFluenceITK
+ {this, "Barrel3DFluenceITK", {0.0e14}, "Barrel3D fluence for radiation damage in ITK"};
+
+ Gaudi::Property<std::vector<std::string>> m_3DFluenceMapITK
+ {this, "Barrel3DFluenceMapITK", {"PixelDigitization/TCAD_IBL_3Dsensors_efields/phi_5e15_160V.root"},
+ "Barrel3D fluence map for radiation damage in ITK"};
+
//====================================================================================
// The following parameters are default values which will be overwritten by the one
// from the conditions DB. Otherwise the DB is not retrieved nor available, these
diff --git a/InnerDetector/InDetConditions/PixelConditionsData/PixelConditionsData/PixelModuleData.h b/InnerDetector/InDetConditions/PixelConditionsData/PixelConditionsData/PixelModuleData.h
index 4d9fa5c2fad..80dd8c1827e 100644
--- a/InnerDetector/InDetConditions/PixelConditionsData/PixelConditionsData/PixelModuleData.h
+++ b/InnerDetector/InDetConditions/PixelConditionsData/PixelConditionsData/PixelModuleData.h
@@ -170,17 +170,44 @@ class PixelModuleData {
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<PixelHistoConverter> distanceMap_e);
void setDistanceMap_h(std::vector<PixelHistoConverter> distanceMap_h);
+ void setRamoPotentialMap(std::vector<PixelHistoConverter> ramoPotentialMap);
+
+ const PixelHistoConverter& getLorentzMap_e(int layer) const;
+ const PixelHistoConverter& getLorentzMap_h(int layer) const;
const PixelHistoConverter& getDistanceMap_e(int layer) const;
const PixelHistoConverter& getDistanceMap_h(int layer) const;
-
- void setRamoPotentialMap(std::vector<PixelHistoConverter> ramoPotentialMap);
const PixelHistoConverter& getRamoPotentialMap(int layer) const;
+ // Map for radiation damage simulation for 3D sensor
+ // The implementation for 3D radiation damage is different from the one for planar sensor.
+ // Thus, define separately. In future, coherent treatment is preferrable.
+ void setFluenceLayer3D(std::vector<double> fluenceLayer);
+ double getFluenceLayer3D(int layer) const;
+
+ void setRamoPotentialMap3D(std::vector<PixelHistoConverter> ramoPotentialMap3D);
+ void setEFieldMap3D(std::vector<PixelHistoConverter> eFieldMap3D);
+ void setXPositionMap3D_e(std::vector<PixelHistoConverter> xPositionMap3D_e);
+ void setXPositionMap3D_h(std::vector<PixelHistoConverter> xPositionMap3D_h);
+ void setYPositionMap3D_e(std::vector<PixelHistoConverter> yPositionMap3D_e);
+ void setYPositionMap3D_h(std::vector<PixelHistoConverter> yPositionMap3D_h);
+ void setTimeMap3D_e(std::vector<PixelHistoConverter> timeMap3D_e);
+ void setTimeMap3D_h(std::vector<PixelHistoConverter> timeMap3D_h);
+ void setAvgChargeMap3D_e(PixelHistoConverter avgChargeMap3D_e);
+ void setAvgChargeMap3D_h(PixelHistoConverter avgChargeMap3D_h);
+
+ const PixelHistoConverter& getRamoPotentialMap3D(int layer) const;
+ const PixelHistoConverter& getEFieldMap3D(int layer) const;
+ const PixelHistoConverter& getXPositionMap3D_e(int layer) const;
+ const PixelHistoConverter& getXPositionMap3D_h(int layer) const;
+ const PixelHistoConverter& getYPositionMap3D_e(int layer) const;
+ const PixelHistoConverter& getYPositionMap3D_h(int layer) const;
+ const PixelHistoConverter& getTimeMap3D_e(int layer) const;
+ const PixelHistoConverter& getTimeMap3D_h(int layer) const;
+ const PixelHistoConverter& getAvgChargeMap3D_e() const;
+ const PixelHistoConverter& getAvgChargeMap3D_h() const;
+
// Distortion parameters
void setDistortionInputSource(int distortionInputSource);
int getDistortionInputSource() const;
@@ -309,6 +336,18 @@ class PixelModuleData {
std::vector<PixelHistoConverter> m_distanceMap_h;
std::vector<PixelHistoConverter> m_ramoPotentialMap;
+ std::vector<double> m_fluenceLayer3D;
+ std::vector<PixelHistoConverter> m_ramoPotentialMap3D;
+ std::vector<PixelHistoConverter> m_eFieldMap3D;
+ std::vector<PixelHistoConverter> m_xPositionMap3D_e;
+ std::vector<PixelHistoConverter> m_xPositionMap3D_h;
+ std::vector<PixelHistoConverter> m_yPositionMap3D_e;
+ std::vector<PixelHistoConverter> m_yPositionMap3D_h;
+ std::vector<PixelHistoConverter> m_timeMap3D_e;
+ std::vector<PixelHistoConverter> m_timeMap3D_h;
+ PixelHistoConverter m_avgChargeMap3D_e;
+ PixelHistoConverter m_avgChargeMap3D_h;
+
int m_distortionInputSource;
int m_distortionVersion;
double m_distortionR1;
diff --git a/InnerDetector/InDetConditions/PixelConditionsData/src/PixelModuleData.cxx b/InnerDetector/InDetConditions/PixelConditionsData/src/PixelModuleData.cxx
index b815435be84..d5fd5dc462e 100644
--- a/InnerDetector/InDetConditions/PixelConditionsData/src/PixelModuleData.cxx
+++ b/InnerDetector/InDetConditions/PixelConditionsData/src/PixelModuleData.cxx
@@ -321,16 +321,42 @@ double PixelModuleData::getFluenceLayer(int layer) const { return m_fluenceLayer
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; }
+void PixelModuleData::setRamoPotentialMap(std::vector<PixelHistoConverter> ramoPotentialMap) { m_ramoPotentialMap = ramoPotentialMap; }
+
+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); }
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); }
+
+// Map for radiation damage simulation for 3D sensor
+void PixelModuleData::setFluenceLayer3D(std::vector<double> fluenceLayer) { m_fluenceLayer3D = fluenceLayer; }
+double PixelModuleData::getFluenceLayer3D(int layer) const { return m_fluenceLayer3D.at(layer); }
+
+void PixelModuleData::setRamoPotentialMap3D(std::vector<PixelHistoConverter> ramoPotentialMap3D) { m_ramoPotentialMap3D = ramoPotentialMap3D; }
+void PixelModuleData::setEFieldMap3D(std::vector<PixelHistoConverter> eFieldMap3D) { m_eFieldMap3D = eFieldMap3D; }
+void PixelModuleData::setXPositionMap3D_e(std::vector<PixelHistoConverter> xPositionMap3D_e) { m_xPositionMap3D_e = xPositionMap3D_e; }
+void PixelModuleData::setXPositionMap3D_h(std::vector<PixelHistoConverter> xPositionMap3D_h) { m_xPositionMap3D_h = xPositionMap3D_h; }
+void PixelModuleData::setYPositionMap3D_e(std::vector<PixelHistoConverter> yPositionMap3D_e) { m_yPositionMap3D_e = yPositionMap3D_e; }
+void PixelModuleData::setYPositionMap3D_h(std::vector<PixelHistoConverter> yPositionMap3D_h) { m_yPositionMap3D_h = yPositionMap3D_h; }
+void PixelModuleData::setTimeMap3D_e(std::vector<PixelHistoConverter> timeMap3D_e) { m_timeMap3D_e = timeMap3D_e; }
+void PixelModuleData::setTimeMap3D_h(std::vector<PixelHistoConverter> timeMap3D_h) { m_timeMap3D_h = timeMap3D_h; }
+void PixelModuleData::setAvgChargeMap3D_e(PixelHistoConverter avgChargeMap3D_e) { m_avgChargeMap3D_e = avgChargeMap3D_e; }
+void PixelModuleData::setAvgChargeMap3D_h(PixelHistoConverter avgChargeMap3D_h) { m_avgChargeMap3D_h = avgChargeMap3D_h; }
+
+const PixelHistoConverter& PixelModuleData::getRamoPotentialMap3D(int layer) const { return m_ramoPotentialMap3D.at(layer); }
+const PixelHistoConverter& PixelModuleData::getEFieldMap3D(int layer) const { return m_eFieldMap3D.at(layer); }
+const PixelHistoConverter& PixelModuleData::getXPositionMap3D_e(int layer) const { return m_xPositionMap3D_e.at(layer); }
+const PixelHistoConverter& PixelModuleData::getXPositionMap3D_h(int layer) const { return m_xPositionMap3D_h.at(layer); }
+const PixelHistoConverter& PixelModuleData::getYPositionMap3D_e(int layer) const { return m_yPositionMap3D_e.at(layer); }
+const PixelHistoConverter& PixelModuleData::getYPositionMap3D_h(int layer) const { return m_yPositionMap3D_h.at(layer); }
+const PixelHistoConverter& PixelModuleData::getTimeMap3D_e(int layer) const { return m_timeMap3D_e.at(layer); }
+const PixelHistoConverter& PixelModuleData::getTimeMap3D_h(int layer) const { return m_timeMap3D_h.at(layer); }
+const PixelHistoConverter& PixelModuleData::getAvgChargeMap3D_e() const { return m_avgChargeMap3D_e; }
+const PixelHistoConverter& PixelModuleData::getAvgChargeMap3D_h() const { return m_avgChargeMap3D_h; }
+
// Distortion parameters
void PixelModuleData::setDistortionInputSource(int distortionInputSource) { m_distortionInputSource = distortionInputSource; }
int PixelModuleData::getDistortionInputSource() const { return m_distortionInputSource; }
diff --git a/InnerDetector/InDetDigitization/PixelDigitization/python/PixelDigitizationConfig.py b/InnerDetector/InDetDigitization/PixelDigitization/python/PixelDigitizationConfig.py
index 6b392f08141..4cbfa0f581c 100644
--- a/InnerDetector/InDetDigitization/PixelDigitization/python/PixelDigitizationConfig.py
+++ b/InnerDetector/InDetDigitization/PixelDigitization/python/PixelDigitizationConfig.py
@@ -60,8 +60,6 @@ def SensorSim3DTool(name="SensorSim3DTool", **kwargs):
from AthenaCommon.AppMgr import ToolSvc
kwargs.setdefault("SiPropertiesTool", ToolSvc.PixelSiPropertiesTool)
kwargs.setdefault("doRadDamage", digitizationFlags.doRadiationDamage.get_Value())
- # TODO This is 2018 fluence setting. These parameters should be controlled by the conditions data.
- kwargs.setdefault("fluence", 6)
return CfgMgr.SensorSim3DTool(name, **kwargs)
def SensorSimTool(name="SensorSimTool", **kwargs):
diff --git a/InnerDetector/InDetDigitization/PixelDigitization/python/PixelDigitizationConfigNew.py b/InnerDetector/InDetDigitization/PixelDigitization/python/PixelDigitizationConfigNew.py
index 3d4ce6823dd..81a251ddd5b 100644
--- a/InnerDetector/InDetDigitization/PixelDigitization/python/PixelDigitizationConfigNew.py
+++ b/InnerDetector/InDetDigitization/PixelDigitization/python/PixelDigitizationConfigNew.py
@@ -1,6 +1,6 @@
"""Define methods to construct configured Pixel Digitization tools and algorithms
-Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
"""
from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
from AthenaConfiguration.ComponentFactory import CompFactory
@@ -97,8 +97,6 @@ def SensorSim3DToolCfg(flags, name="SensorSim3DTool", **kwargs):
kwargs.setdefault("SiPropertiesTool", SiTool)
SensorSim3DTool = CompFactory.SensorSim3DTool
kwargs.setdefault("doRadDamage", flags.Digitization.DoRadiationDamage)
- # TODO This is 2018 fluence setting. These parameters should be controlled by the conditions data.
- kwargs.setdefault("fluence", 6)
acc.setPrivateTools(SensorSim3DTool(name, **kwargs))
return acc
diff --git a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.cxx b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.cxx
index ae287ed87c2..5702d134ad9 100644
--- a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.cxx
+++ b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.cxx
@@ -48,133 +48,6 @@ StatusCode SensorSim3DTool::initialize() {
ATH_CHECK(readProbMap(m_cc_prob_file_fei3));
ATH_CHECK(readProbMap(m_cc_prob_file_fei4));
- std::vector<std::string> mapsPath_list;
-
- if (!m_doRadDamage) {
- m_fluence = 0;
- }
-
- if (static_cast<int>(m_fluence) == 1) {
- mapsPath_list.push_back(PathResolverFindCalibFile("PixelDigitization/TCAD_IBL_3Dsensors_efields/phi_0_20V.root"));
- m_fluence_layers.push_back(1e-10);
- } else if (static_cast<int>(m_fluence) == 2) {
- mapsPath_list.push_back(PathResolverFindCalibFile("PixelDigitization/TCAD_IBL_3Dsensors_efields/phi_1e14_20V.root"));
- m_fluence_layers.push_back(1e14);
- } else if (static_cast<int>(m_fluence) == 3) {
- mapsPath_list.push_back(PathResolverFindCalibFile("PixelDigitization/TCAD_IBL_3Dsensors_efields/phi_2e14_30V.root"));
- m_fluence_layers.push_back(2e14);
- } else if (static_cast<int>(m_fluence) == 4) {
- mapsPath_list.push_back(PathResolverFindCalibFile("PixelDigitization/TCAD_IBL_3Dsensors_efields/phi_5e14_40V.root"));
- m_fluence_layers.push_back(5e14);
- } else if (static_cast<int>(m_fluence) == 5) {
- mapsPath_list.push_back(PathResolverFindCalibFile("PixelDigitization/TCAD_IBL_3Dsensors_efields/phi_1e15_50V.root"));
- m_fluence_layers.push_back(1e15);
- } else if (static_cast<int>(m_fluence) == 6) {
- mapsPath_list.push_back(PathResolverFindCalibFile("PixelDigitization/TCAD_IBL_3Dsensors_efields/phi_5e15_160V.root"));
- m_fluence_layers.push_back(5e15);
- } else if (static_cast<int>(m_fluence) == 7) {
- mapsPath_list.push_back(PathResolverFindCalibFile(
- "PixelDigitization/TCAD_IBL_3Dsensors_efields/phi_6e15_190V_new.root"));
- m_fluence_layers.push_back(6e15);
- } else if (static_cast<int>(m_fluence) == 8) {
- mapsPath_list.push_back(PathResolverFindCalibFile(
- "PixelDigitization/TCAD_IBL_3Dsensors_efields/phi_1e16_260V_new.root"));
- m_fluence_layers.push_back(1e16);
- }
-
- // *****************************
- // *** Setup Maps ****
- // *****************************
- //TODO This is only temporary until remotely stored maps and locally generated maps can be implemented
-
- for (unsigned int i = 0; i < mapsPath_list.size(); i++) {
- ATH_MSG_INFO("Using maps located in: " << mapsPath_list.at(i));
- std::unique_ptr<TFile> mapsFile(TFile::Open((mapsPath_list.at(i)).c_str(), "READ")); //this is the ramo potential.
- if (!mapsFile) {
- ATH_MSG_ERROR("Cannot open file: " << mapsPath_list.at(i));
- return StatusCode::FAILURE;
- }
-
- 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
- std::unique_ptr<TH2F> ramoPotentialMap_hold(mapsFile->Get<TH2F>("ramo"));
- if (!ramoPotentialMap_hold) {
- ATH_MSG_ERROR("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_hold->SetDirectory(nullptr);
- m_ramoPotentialMap.emplace_back();
- ATH_CHECK(m_ramoPotentialMap.back().setHisto2D(ramoPotentialMap_hold.get()));
-
- ATH_MSG_INFO("Using fluence " << m_fluence_layers.at(i));
-
- //Now setup the E-field.
- std::unique_ptr<TH2F> eFieldMap_hold(mapsFile->Get<TH2F>("efield"));
- if (!eFieldMap_hold) {
- ATH_MSG_ERROR("Unable to load sensor e-field map, so generating one using approximations.");
- return StatusCode::FAILURE; //Obviously, remove this when gen. code is set up
- }
- eFieldMap_hold->SetDirectory(nullptr);
- m_eFieldMap.emplace_back();
- ATH_CHECK(m_eFieldMap.back().setHisto2D(eFieldMap_hold.get()));
-
- std::unique_ptr<TH3F> xPositionMap_e_hold(mapsFile->Get<TH3F>("xPosition_e"));
- std::unique_ptr<TH3F> xPositionMap_h_hold(mapsFile->Get<TH3F>("xPosition_h"));
- std::unique_ptr<TH3F> yPositionMap_e_hold(mapsFile->Get<TH3F>("yPosition_e"));
- std::unique_ptr<TH3F> yPositionMap_h_hold(mapsFile->Get<TH3F>("yPosition_h"));
- std::unique_ptr<TH2F> timeMap_e_hold(mapsFile->Get<TH2F>("etimes"));
- std::unique_ptr<TH2F> timeMap_h_hold(mapsFile->Get<TH2F>("htimes"));
-
- if (!xPositionMap_e_hold || !xPositionMap_h_hold || !yPositionMap_e_hold || !yPositionMap_h_hold ||
- !timeMap_e_hold || !timeMap_h_hold) {
- ATH_MSG_ERROR("Cannot find one of the maps.");
- return StatusCode::FAILURE; //Obviously, remove this when gen. code is set up
- }
-
- xPositionMap_e_hold->SetDirectory(nullptr);
- xPositionMap_h_hold->SetDirectory(nullptr);
- yPositionMap_e_hold->SetDirectory(nullptr);
- yPositionMap_h_hold->SetDirectory(nullptr);
- timeMap_e_hold->SetDirectory(nullptr);
- timeMap_h_hold->SetDirectory(nullptr);
-
- //Now, determine the time to reach the electrode and the trapping position.
- 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.get()));
- ATH_CHECK(m_xPositionMap_h.back().setHisto3D(xPositionMap_h_hold.get()));
- ATH_CHECK(m_yPositionMap_e.back().setHisto3D(yPositionMap_e_hold.get()));
- ATH_CHECK(m_yPositionMap_h.back().setHisto3D(yPositionMap_h_hold.get()));
- m_timeMap_e.emplace_back();
- m_timeMap_h.emplace_back();
- ATH_CHECK(m_timeMap_e.back().setHisto2D(timeMap_e_hold.get()));
- ATH_CHECK(m_timeMap_h.back().setHisto2D(timeMap_h_hold.get()));
-
- std::unique_ptr<TH2F> avgCharge_e(mapsFile->Get<TH2F>("avgCharge_e"));
- std::unique_ptr<TH2F> avgCharge_h(mapsFile->Get<TH2F>("avgCharge_h"));
-
- if (!avgCharge_e || !avgCharge_h) {
- ATH_MSG_ERROR("Cannot find one of the charge maps.");
- return StatusCode::FAILURE; //Obviously, remove this when gen. code is set up
- }
-
- avgCharge_e->SetDirectory(nullptr);
- avgCharge_h->SetDirectory(nullptr);
-
- // Get average charge data (for charge chunk effect correction)
- ATH_CHECK(m_avgChargeMap_e.setHisto2D(avgCharge_e.get()));
- ATH_CHECK(m_avgChargeMap_h.setHisto2D(avgCharge_h.get()));
-
- mapsFile->Close();
- }
-
return StatusCode::SUCCESS;
}
@@ -193,7 +66,7 @@ StatusCode SensorSim3DTool::induceCharge(const TimedHitPtr<SiHit>& phit,
SiChargedDiodeCollection& chargedDiodes,
const InDetDD::SiDetectorElement& Module,
const InDetDD::PixelModuleDesign& p_design,
- const PixelModuleData */*moduleData*/,
+ const PixelModuleData *moduleData,
std::vector< std::pair<double, double> >& trfHitRecord,
std::vector<double>& initialConditions,
CLHEP::HepRandomEngine* rndmEngine,
@@ -216,10 +89,20 @@ StatusCode SensorSim3DTool::induceCharge(const TimedHitPtr<SiHit>& phit,
//Calculate trapping times based on fluence (already includes check for fluence=0)
if (m_doRadDamage) {
- std::pair < double, double > trappingTimes = m_radDamageUtil->getTrappingTimes(m_fluence_layers.at(0)); //0 = IBL
+ std::pair < double, double > trappingTimes = m_radDamageUtil->getTrappingTimes(moduleData->getFluenceLayer3D(0)); //0 = IBL
m_trappingTimeElectrons = trappingTimes.first;
m_trappingTimeHoles = trappingTimes.second;
}
+ const PixelHistoConverter& ramoPotentialMap = moduleData->getRamoPotentialMap3D(0);
+ const PixelHistoConverter& eFieldMap = moduleData->getEFieldMap3D(0);
+ const PixelHistoConverter& xPositionMap_e = moduleData->getXPositionMap3D_e(0);
+ const PixelHistoConverter& xPositionMap_h = moduleData->getXPositionMap3D_h(0);
+ const PixelHistoConverter& yPositionMap_e = moduleData->getYPositionMap3D_e(0);
+ const PixelHistoConverter& yPositionMap_h = moduleData->getYPositionMap3D_h(0);
+ const PixelHistoConverter& timeMap_e = moduleData->getTimeMap3D_e(0);
+ const PixelHistoConverter& timeMap_h = moduleData->getTimeMap3D_h(0);
+ const PixelHistoConverter& avgChargeMap_e = moduleData->getAvgChargeMap3D_e();
+ const PixelHistoConverter& avgChargeMap_h = moduleData->getAvgChargeMap3D_h();
double eta_0 = initialConditions[0];
double phi_0 = initialConditions[1];
@@ -250,7 +133,7 @@ StatusCode SensorSim3DTool::induceCharge(const TimedHitPtr<SiHit>& phit,
const EBC_EVCOLL evColl = EBC_MAINEVCOLL;
const HepMcParticleLink::PositionFlag idxFlag =
(phit.eventId() == 0) ? HepMcParticleLink::IS_POSITION : HepMcParticleLink::IS_INDEX;
- if (m_doRadDamage && m_fluence > 0) {
+ if (m_doRadDamage) {
//**************************************//
//*** Now diffuse charges to surface *** //
//**************************************//
@@ -288,7 +171,6 @@ StatusCode SensorSim3DTool::induceCharge(const TimedHitPtr<SiHit>& phit,
double x_new = chargepos.x() + 0.5*module_size_x;
double y_new = chargepos.y() + 0.5*module_size_y;
-
// -- change from module frame to pixel frame
int nPixX = int(x_new / pixel_size_x);
int nPixY = int(y_new / pixel_size_y);
@@ -308,7 +190,8 @@ StatusCode SensorSim3DTool::induceCharge(const TimedHitPtr<SiHit>& phit,
//only process hits which are not on the electrodes (E-field zero)
//all the maps have 250 as the x value, so need to invert x and y whenever reading maps
- double efield = getElectricField(y_pix, x_pix);
+ double efield = eFieldMap.getContent(eFieldMap.getBinX(1e3*y_pix),eFieldMap.getBinY(1e3*x_pix))*1.0E-7; //return efield in MV/mm (for mobility calculation)
+
if (efield == 0) {
ATH_MSG_DEBUG("Skipping since efield = 0 for x_pix = " << x_pix << " y_pix = " << y_pix);
continue;
@@ -335,7 +218,14 @@ StatusCode SensorSim3DTool::induceCharge(const TimedHitPtr<SiHit>& phit,
extraNPixX = nPixX;
extraNPixY = nPixY;
- const double timeToElectrode = getTimeToElectrode(y_pix, x_pix, isHole);
+ double timeToElectrode(0.0); //[ns]
+ if (!isHole) {
+ timeToElectrode = timeMap_e.getContent(timeMap_e.getBinX(1e3*y_pix),timeMap_e.getBinY(1e3*x_pix));
+ }
+ else {
+ timeToElectrode = timeMap_h.getContent(timeMap_h.getBinX(1e3*y_pix),timeMap_h.getBinY(1e3*x_pix));
+ }
+
const double driftTime = getDriftTime(isHole);
//Apply drift due to diffusion
@@ -366,28 +256,35 @@ StatusCode SensorSim3DTool::induceCharge(const TimedHitPtr<SiHit>& phit,
yposDiff = yposDiff + pixel_size_y;
}
-
- 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));
-
- double xposFinal = getTrappingPositionY(yposDiff, xposDiff, std::min(driftTime, timeToElectrode), isHole);
- double yposFinal = getTrappingPositionX(yposDiff, xposDiff, std::min(driftTime, timeToElectrode), isHole);
+ float average_charge = isHole ? avgChargeMap_h.getContent(avgChargeMap_h.getBinY(1e3*y_pix),avgChargeMap_h.getBinX(1e3*x_pix)) :
+ avgChargeMap_e.getContent(avgChargeMap_e.getBinY(1e3*y_pix),avgChargeMap_e.getBinX(1e3*x_pix));
+
+ double xposFinal(0); // Trapping position Y
+ double yposFinal(0); // Trapping position X
+ if (!isHole) {
+ xposFinal = yPositionMap_e.getContent(yPositionMap_e.getBinX(1e3*yposDiff),yPositionMap_e.getBinY(1e3*xposDiff),yPositionMap_e.getBinZ(std::min(driftTime,timeToElectrode)))*1e-3; //[mm]
+ yposFinal = xPositionMap_e.getContent(xPositionMap_e.getBinX(1e3*yposDiff),xPositionMap_e.getBinY(1e3*xposDiff),xPositionMap_e.getBinZ(std::min(driftTime,timeToElectrode)))*1e-3; //[mm]
+ }
+ else {
+ xposFinal = yPositionMap_h.getContent(yPositionMap_h.getBinX(1e3*yposDiff),yPositionMap_h.getBinY(1e3*xposDiff),yPositionMap_h.getBinZ(std::min(driftTime,timeToElectrode)))*1e-3; //[mm]
+ yposFinal = xPositionMap_h.getContent(xPositionMap_h.getBinX(1e3*yposDiff),xPositionMap_h.getBinY(1e3*xposDiff),xPositionMap_h.getBinZ(std::min(driftTime,timeToElectrode)))*1e-3; //[mm]
+ }
// -- Calculate signal in current pixel and in the neighboring ones
// -- loop in the x-coordinate
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));
+ const std::size_t ramo_init_bin_y = ramoPotentialMap.getBinY(1000*(x_pix + pixel_size_x * 3 - xNeighbor));
+ const std::size_t ramo_final_bin_y = ramoPotentialMap.getBinY(1000*(xposFinal + pixel_size_x * 3 - xNeighbor));
+
for (int j = -1; j <= 1; j++) {
double yNeighbor = j * pixel_size_y;
//Ramo map over 500umx350um pixel area
//Ramo init different if charge diffused into neighboring pixel -> change primary pixel!!
- 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);
+ float ramoInit = ramoPotentialMap.getContent(ramoPotentialMap.getBinX(1000*(y_pix + 0.5*pixel_size_y - yNeighbor)), ramo_init_bin_y);
+ float ramoFinal = ramoPotentialMap.getContent(ramoPotentialMap.getBinX(1000*(yposFinal + 0.5*pixel_size_y - yNeighbor)), ramo_final_bin_y);
// Record deposit
double eHitRamo = (isHole ? -1. : 1.) * eHit * (ramoFinal - ramoInit);
@@ -587,12 +484,6 @@ double SensorSim3DTool::getProbMapEntry(const SensorType& readout, int binx, int
return echarge;
}
-double SensorSim3DTool::getElectricField(double x, double y) {
- 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)
-}
-
double SensorSim3DTool::getMobility(double electricField, bool isHoleBit) {
//Not exactly the same as the getMobility function in RadDamageUtil, since we don't have a Hall effect for 3D sensors
// (B and E are parallel)
@@ -631,37 +522,4 @@ double SensorSim3DTool::getDriftTime(bool isHoleBit) {
return driftTime;
}
-double SensorSim3DTool::getTimeToElectrode(double x, double y, bool isHoleBit) {
- std::size_t index = 0;
- double timeToElectrode = 0;
- if (!isHoleBit) {
- 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::size_t index = 0;
- double finalX(0);
- if (!isHoleBit) {
- 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 {
- 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::size_t index = 0;
- double finalY(0);
- if (!isHoleBit) {
- 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 {
- 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]
-}
diff --git a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.h b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.h
index 91bda4541b9..42ad264bd48 100644
--- a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.h
+++ b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.h
@@ -42,12 +42,9 @@ public:
StatusCode readProbMap(const std::string&);
StatusCode printProbMap(const std::string&) const;
- double getElectricField(double x, double y);
double getMobility(double electricField, bool isHoleBit);
double getDriftTime(bool isHoleBit);
- 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);
+
private:
enum class SensorType {
@@ -63,20 +60,6 @@ private:
std::multimap<std::pair<int, int>, double> m_probMapFEI4;
std::multimap<std::pair<int, int>, double> m_probMapFEI3;
- // Map for radiation damage simulation
- 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;
-
Gaudi::Property<std::string> m_cc_prob_file_fei3
{
this, "CCProbMapFileFEI3", "PixelDigitization/3DFEI3-3E-problist-1um_v1.txt",
@@ -104,12 +87,6 @@ private:
this, "doChunkCorrection", false, "doChunkCorrection bool: should be flag"
};
- Gaudi::Property<double> m_fluence
- {
- this, "fluence", 0,
- "this is the fluence benchmark, 0-6. 0 is unirradiated, 1 is start of Run 2, 5 is end of 2018 and 6 is projected end of 2018"
- };
-
Gaudi::Property<double> m_trappingTimeElectrons
{
this, "trappingTimeElectrons", 0.0, "Characteristic time till electron is trapped [ns]"
--
GitLab