From 622e95bd747dba67b5f453748fb66468d72aa730 Mon Sep 17 00:00:00 2001
From: Dave Casper <dcasper@localhost.localdomain>
Date: Wed, 15 Jan 2020 23:59:15 -0800
Subject: [PATCH] First version of SCT_Digitization
---
.../FaserSCT_Digitization/CMakeLists.txt | 4 +-
.../FaserSCT_Digitization/src/SCT_Amp.cxx | 143 +++
.../FaserSCT_Digitization/src/SCT_Amp.h | 69 ++
.../src/SCT_Digitization.cxx | 31 +
.../src/SCT_Digitization.h | 48 +
.../src/SCT_FrontEnd.cxx | 1044 +++++++++++++++++
.../FaserSCT_Digitization/src/SCT_FrontEnd.h | 140 +++
.../src/SCT_RandomDisabledCellGenerator.cxx | 40 +
.../src/SCT_RandomDisabledCellGenerator.h | 67 ++
.../src/SCT_SurfaceChargesGenerator.cxx | 578 +++++++++
.../src/SCT_SurfaceChargesGenerator.h | 164 +++
.../FaserSCT_Digitization_entries.cxx | 24 +-
.../TrackerSimEvent/FaserSiHit.h | 4 +-
13 files changed, 2341 insertions(+), 15 deletions(-)
create mode 100644 Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_Amp.cxx
create mode 100644 Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_Amp.h
create mode 100644 Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_Digitization.cxx
create mode 100644 Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_Digitization.h
create mode 100644 Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_FrontEnd.cxx
create mode 100644 Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_FrontEnd.h
create mode 100644 Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_RandomDisabledCellGenerator.cxx
create mode 100644 Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_RandomDisabledCellGenerator.h
create mode 100644 Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_SurfaceChargesGenerator.cxx
create mode 100644 Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_SurfaceChargesGenerator.h
diff --git a/Tracker/TrackerDigitization/FaserSCT_Digitization/CMakeLists.txt b/Tracker/TrackerDigitization/FaserSCT_Digitization/CMakeLists.txt
index e836c3e3..d0965e8d 100644
--- a/Tracker/TrackerDigitization/FaserSCT_Digitization/CMakeLists.txt
+++ b/Tracker/TrackerDigitization/FaserSCT_Digitization/CMakeLists.txt
@@ -24,7 +24,7 @@ atlas_depends_on_subdirs( PUBLIC
Generators/GeneratorObjects
#InnerDetector/InDetConditions/InDetConditionsSummaryService
#InnerDetector/InDetConditions/SCT_ConditionsTools
- #InnerDetector/InDetConditions/SiPropertiesTool
+ InnerDetector/InDetConditions/SiPropertiesTool
Tracker/TrackerDetDescr/TrackerIdentifier
Tracker/TrackerDetDescr/TrackerReadoutGeometry
#InnerDetector/InDetDetDescr/SCT_ModuleDistortions
@@ -40,7 +40,7 @@ atlas_add_component( FaserSCT_Digitization
src/*.cxx
src/components/*.cxx
INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${Boost_INCLUDE_DIRS} ${CLHEP_INCLUDE_DIRS}
- LINK_LIBRARIES ${ROOT_LIBRARIES} ${Boost_LIBRARIES} ${CLHEP_LIBRARIES} AthenaBaseComps AthenaKernel PileUpToolsLib Identifier xAODEventInfo GaudiKernel FaserSiDigitization InDetRawData TrackerSimEvent HitManagement GeneratorObjects TrackerIdentifier TrackerReadoutGeometry InDetSimData )
+ LINK_LIBRARIES ${ROOT_LIBRARIES} ${Boost_LIBRARIES} ${CLHEP_LIBRARIES} AthenaBaseComps AthenaKernel PileUpToolsLib Identifier xAODEventInfo GaudiKernel FaserSiDigitization InDetRawData TrackerSimEvent HitManagement GeneratorObjects SiPropertiesToolLib TrackerIdentifier TrackerReadoutGeometry InDetSimData )
# LINK_LIBRARIES ${ROOT_LIBRARIES} ${Boost_LIBRARIES} ${CLHEP_LIBRARIES} CommissionEvent AthenaBaseComps AthenaKernel PileUpToolsLib Identifier xAODEventInfo GaudiKernel SiDigitization InDetRawData InDetSimEvent HitManagement GeneratorObjects SiPropertiesToolLib InDetIdentifier InDetReadoutGeometry InDetSimData )
#atlas_add_test( SCT_DigitizationMT_test
diff --git a/Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_Amp.cxx b/Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_Amp.cxx
new file mode 100644
index 00000000..22352b26
--- /dev/null
+++ b/Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_Amp.cxx
@@ -0,0 +1,143 @@
+/*
+ Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "SCT_Amp.h"
+
+// CLHEP
+#include "CLHEP/Units/SystemOfUnits.h"
+
+//STD includes
+#include <cmath>
+#include <fstream>
+
+//#define SCT_DIG_DEBUG
+
+// constructor
+SCT_Amp::SCT_Amp(const std::string& type, const std::string& name, const IInterface* parent)
+ : base_class(type, name, parent)
+{
+}
+
+//----------------------------------------------------------------------
+// Initialize
+//----------------------------------------------------------------------
+StatusCode SCT_Amp::initialize() {
+
+ StatusCode sc{AthAlgTool::initialize()};
+ if (sc.isFailure()) {
+ ATH_MSG_FATAL("SCT_Amp::initialize() failed");
+ return sc;
+ }
+ ATH_MSG_DEBUG("SCT_Amp::initialize()");
+
+ /** CHLEP Units */
+ m_PeakTime.setValue(m_PeakTime.value() * CLHEP::ns);
+ m_dt.setValue(m_dt.value() * CLHEP::ns);
+ m_tmin.setValue(m_tmin.value() * CLHEP::ns);
+ m_tmax.setValue(m_tmax.value() * CLHEP::ns);
+
+ m_NormConstCentral = (exp(3.0)/27.0)*(1.0-m_CrossFactor2sides)*(1.0-m_CrossFactorBack);
+ m_NormConstNeigh = exp(3.0-sqrt(3.0))/(6*(2.0*sqrt(3.0)-3.0));
+ m_NormConstNeigh *= (m_CrossFactor2sides/2.0)*(1.0-m_CrossFactorBack);
+
+#ifdef SCT_DIG_DEBUG
+ ATH_MSG_INFO("\tAmp created, PeakTime = " << m_PeakTime);
+ ATH_MSG_INFO("\tResponse will be CR-RC^3 with tp = " << m_PeakTime/3.0);
+ ATH_MSG_INFO("\tCoupling to both neighbours = " << m_CrossFactor2sides);
+ ATH_MSG_INFO("\tCoupling to backplane = " << m_CrossFactorBack);
+ ATH_MSG_INFO("\tNormalization for central " << m_NormConstCentral);
+ ATH_MSG_INFO("\tNormalization for neighbour " << m_NormConstNeigh);
+#endif
+
+ return sc;
+}
+
+//----------------------------------------------------------------------
+// Finalize
+//----------------------------------------------------------------------
+StatusCode SCT_Amp::finalize() {
+ StatusCode sc{AthAlgTool::finalize()};
+ if (sc.isFailure()) {
+ ATH_MSG_FATAL("SCT_Amp::finalize() failed");
+ return sc;
+ }
+ ATH_MSG_DEBUG("SCT_Amp::finalize()");
+ return sc;
+}
+
+//----------------------------------------------------------------------
+// Electronique response is now CR-RC^3 of the charge diode
+//----------------------------------------------------------------------
+float SCT_Amp::response(const list_t& Charges, const float timeOfThreshold) const {
+ float resp{0.0};
+ float tp{static_cast<float>(m_PeakTime/3.0)}; // for CR-RC^3
+ for (const SiCharge& charge: Charges) {
+ float ch{static_cast<float>(charge.charge())};
+ float tC{static_cast<float>(timeOfThreshold - charge.time())};
+ if (tC > 0.0) {
+ tC/=tp; //to avoid doing it four times
+ resp += ch*tC*tC*tC*exp(-tC); //faster than pow
+ }
+ }
+ return resp*m_NormConstCentral;
+}
+
+void SCT_Amp::response(const list_t& Charges, const float timeOfThreshold, std::vector<float>& response) const {
+ short bin_max{static_cast<short>(response.size())};
+ std::fill(response.begin(), response.end(), 0.0);
+ float tp{static_cast<float>(m_PeakTime/3.0)}; // for CR-RC^3
+ for (const SiCharge& charge: Charges) {
+ float ch{static_cast<float>(charge.charge())};
+ float ch_time{static_cast<float>(charge.time())};
+ short bin_end{static_cast<short>(bin_max-1)};
+ for (short bin{-1}; bin<bin_end; ++bin) {
+ float bin_timeOfThreshold{timeOfThreshold + bin*25};//25, fix me
+ float tC{bin_timeOfThreshold - ch_time};
+ if (tC > 0.0) {
+ tC/=tp; //to avoid doing it four times
+ response[bin+1] += ch*tC*tC*tC*exp(-tC); //faster than pow
+ }
+ }
+ }
+ for (short bin{0}; bin<bin_max; ++bin) response[bin] = response[bin]*m_NormConstCentral;
+ return;
+}
+
+//----------------------------------------------------------------------
+// differenciated and scaled pulse on the neighbour strip!
+//----------------------------------------------------------------------
+float SCT_Amp::crosstalk(const list_t& Charges, const float timeOfThreshold) const {
+ float resp{0};
+ float tp{static_cast<float>(m_PeakTime/3.0)}; // for CR-RC^3
+ for (const SiCharge& charge: Charges) {
+ float ch{static_cast<float>(charge.charge())};
+ float tC{static_cast<float>(timeOfThreshold - charge.time())};
+ if (tC > 0.0) {
+ tC/=tp; //to avoid doing it four times
+ resp += ch*tC*tC*exp(-tC)*(3.0-tC); //faster than pow
+ }
+ }
+ return resp*m_NormConstNeigh;
+}
+
+void SCT_Amp::crosstalk(const list_t& Charges, const float timeOfThreshold, std::vector<float>& response) const {
+ short bin_max{static_cast<short>(response.size())};
+ std::fill(response.begin(), response.end(), 0.0);
+ float tp{static_cast<float>(m_PeakTime/3.0)}; // for CR-RC^3
+ for (const SiCharge& charge: Charges) {
+ float ch{static_cast<float>(charge.charge())};
+ float ch_time{static_cast<float>(charge.time())};
+ short bin_end{static_cast<short>(bin_max-1)};
+ for (short bin{-1}; bin<bin_end; ++bin) {
+ float bin_timeOfThreshold{timeOfThreshold + bin*25}; // 25, fix me
+ float tC{bin_timeOfThreshold - ch_time};
+ if (tC > 0.0) {
+ tC/=tp; //to avoid doing it four times
+ response[bin+1] += ch*tC*tC*exp(-tC)*(3.0-tC); //faster than pow
+ }
+ }
+ }
+ for (short bin{0}; bin<bin_max; ++bin) response[bin] = response[bin]*m_NormConstNeigh;
+ return;
+}
diff --git a/Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_Amp.h b/Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_Amp.h
new file mode 100644
index 00000000..d066e194
--- /dev/null
+++ b/Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_Amp.h
@@ -0,0 +1,69 @@
+// -*- C++ -*-
+
+/*
+ Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+*/
+
+/**
+ * Header file for SCT_Amp Class
+ * @brief A class to model an SCT amplifier and shaper. Gives a CRRC response to a
+ * list of charges with times. Also calculates average input and output for
+ * diagnostic purposes. Questions/comments to Szymon.Gadomski@cern.ch
+ * Name changed from SCTpreamp (misleading) on 09.05.01.
+ * @author Szymon.Gadomski@cern.ch, Awatif.Belymam@cern.ch, Davide.Costanzo@cern.ch,
+ * tgcornel@nikhef.nl, Grant.Gorfine@cern.ch, Paul.Bell@cern.ch, Jorgen.Dalmau@cern.ch,
+ * 23/08/2007 - Kondo.Gnanvo@cern.ch
+ * - Conversion of the SCT_Amp code AlgTool
+ */
+
+#ifndef FASERSCT_DIGITIZATION_SCTAMP_H
+#define FASERSCT_DIGITIZATION_SCTAMP_H
+
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "FaserSCT_Digitization/ISCT_Amp.h"
+
+#include "TrackerSimEvent/SiCharge.h"
+
+class SCT_Amp : public extends<AthAlgTool, ISCT_Amp> {
+ public:
+
+ /** constructor */
+ SCT_Amp(const std::string& type, const std::string& name, const IInterface* parent);
+ /** Destructor */
+ virtual ~SCT_Amp() = default;
+ /** AlgTool initialize */
+ virtual StatusCode initialize() override;
+ /** AlgTool finalize */
+ virtual StatusCode finalize() override;
+
+ /** main purpose: CR-RC^3 response to a list of charges with times */
+ virtual float response(const list_t& Charges, const float timeOverThreshold) const override;
+ virtual void response(const list_t& Charges, const float timeOverThreshold, std::vector<float>& resp) const override;
+
+ /** Neighbour strip cross talk response strip to a list of charges with times */
+ virtual float crosstalk(const list_t& Charges, const float timeOverThreshold) const override;
+ virtual void crosstalk(const list_t& Charges, const float timeOverThreshold, std::vector<float>& resp) const override;
+
+private:
+
+ /** signal peak time */
+ FloatProperty m_PeakTime{this, "PeakTime", 21., "Front End Electronics peaking time"};
+
+ /** Cross factor 2 side */
+ FloatProperty m_CrossFactor2sides{this, "CrossFactor2sides", 0.1, "Loss of charge to neighbour strip constant"};
+
+ /** cross factor */
+ FloatProperty m_CrossFactorBack{this, "CrossFactorBack", 0.07, "Loss of charge to back plane constant"};
+
+ FloatProperty m_tmin{this, "Tmin", -25.0};
+ FloatProperty m_tmax{this, "Tmax", 150.0};
+ FloatProperty m_dt{this, "deltaT", 1.0};
+
+ /** Normalisation factor for the signal response */
+ float m_NormConstCentral{0.};
+
+ /** Normalisation factor for the neighbour strip signal response */
+ float m_NormConstNeigh{0.};
+};
+
+#endif // FASERSCT_DIGITIZATION_SCTAMP_H
diff --git a/Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_Digitization.cxx b/Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_Digitization.cxx
new file mode 100644
index 00000000..f30bbd62
--- /dev/null
+++ b/Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_Digitization.cxx
@@ -0,0 +1,31 @@
+/*
+ Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "SCT_Digitization.h"
+#include "PileUpTools/IPileUpTool.h"
+
+//----------------------------------------------------------------------
+// Constructor with parameters:
+//----------------------------------------------------------------------
+SCT_Digitization::SCT_Digitization(const std::string& name, ISvcLocator* pSvcLocator) :
+ AthAlgorithm(name, pSvcLocator)
+{
+}
+
+//----------------------------------------------------------------------
+// Initialize method:
+//----------------------------------------------------------------------
+StatusCode SCT_Digitization::initialize() {
+ ATH_MSG_DEBUG("SCT_Digitization::initialize()");
+ return StatusCode::SUCCESS ;
+}
+
+//----------------------------------------------------------------------
+// Execute method:
+//----------------------------------------------------------------------
+
+StatusCode SCT_Digitization::execute() {
+ ATH_MSG_DEBUG("execute()");
+ return m_sctDigitizationTool->processAllSubEvents();
+}
diff --git a/Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_Digitization.h b/Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_Digitization.h
new file mode 100644
index 00000000..3f32cdf0
--- /dev/null
+++ b/Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_Digitization.h
@@ -0,0 +1,48 @@
+// -*- C++ -*-
+
+/*
+ Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+*/
+
+/** @file SCT_Digitization.h Header file for SCT_Digitization class.
+ * @author Szymon.Gadomski@cern.ch, Awatif.Belymam@cern.ch, Davide.Costanzo@cern.ch,
+ * tgcornel@nikhef.nl, Grant.Gorfine@cern.ch, Paul.Bell@cern.ch,
+ * Jorgen.Dalmau@cern.ch, Kondo.Gnanvo@cern.ch, and others
+ * Version 23/08/2007 Kondo.Gnanvo@cern.ch
+ * Conversion of the processors into AlgTool
+ */
+
+// Multiple inclusion protection
+#ifndef FASERSCT_DIGITIZATION_SCT_DIGITIZATION_H
+#define FASERSCT_DIGITIZATION_SCT_DIGITIZATION_H
+
+// Base class
+#include "AthenaBaseComps/AthAlgorithm.h"
+// Gaudi
+#include "GaudiKernel/ToolHandle.h"
+
+class IPileUpTool;
+
+/** Top algorithm class for SCT digitization */
+class SCT_Digitization : public AthAlgorithm {
+
+ public:
+
+ /** Constructor with parameters */
+ SCT_Digitization(const std::string& name, ISvcLocator* pSvcLocator);
+
+ /** Destructor */
+ virtual ~SCT_Digitization() = default;
+
+ /** Basic algorithm methods */
+ virtual StatusCode initialize() override final;
+ virtual StatusCode execute() override final;
+ virtual bool isClonable() const override final { return true; }
+
+ private:
+
+ ToolHandle<IPileUpTool> m_sctDigitizationTool{this, "DigitizationTool", "SCT_DigitizationTool", "SCT_DigitizationTool name"};
+
+};
+
+#endif // FASERSCT_DIGITIZATION_SCT_DIGITIZATION_H
diff --git a/Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_FrontEnd.cxx b/Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_FrontEnd.cxx
new file mode 100644
index 00000000..50d06c46
--- /dev/null
+++ b/Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_FrontEnd.cxx
@@ -0,0 +1,1044 @@
+/*
+ Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "SCT_FrontEnd.h"
+
+// Random number
+#include "CLHEP/Random/RandFlat.h"
+#include "CLHEP/Random/RandGaussZiggurat.h" // for RandGaussZiggurat
+#include "CLHEP/Random/RandPoisson.h"
+#include "CLHEP/Random/RandomEngine.h"
+
+#include "FaserSiDigitization/SiHelper.h"
+#include "FaserSCT_Digitization/ISCT_Amp.h"
+
+#include "TrackerIdentifier/FaserSCT_ID.h"
+
+#include "TrackerReadoutGeometry/SCT_DetectorManager.h"
+#include "TrackerReadoutGeometry/SCT_ModuleSideDesign.h"
+
+// C++ Standard Library
+#include <algorithm>
+#include <iostream>
+
+// #define SCT_DIG_DEBUG
+
+using namespace std;
+using namespace TrackerDD;
+
+// constructor
+SCT_FrontEnd::SCT_FrontEnd(const std::string& type, const std::string& name, const IInterface* parent)
+ : base_class(type, name, parent) {
+}
+
+// ----------------------------------------------------------------------
+// Initialize
+// ----------------------------------------------------------------------
+StatusCode SCT_FrontEnd::initialize() {
+ if (m_NoiseOn and (not m_analogueNoiseOn)) {
+ ATH_MSG_FATAL("AnalogueNoiseOn/m_analogueNoiseOn should be true if NoiseOn/m_NoiseOn is true.");
+ return StatusCode::FAILURE;
+ }
+
+ ATH_MSG_DEBUG("SCT_FrontEnd::initialize()");
+ // Get SCT helper
+ ATH_CHECK(detStore()->retrieve(m_sct_id, "FaserSCT_ID"));
+ // Get SCT detector manager
+ ATH_CHECK(detStore()->retrieve(m_SCTdetMgr, "SCT"));
+ // Get the amplifier tool
+ ATH_CHECK(m_sct_amplifier.retrieve());
+ ATH_MSG_DEBUG("SCT Amplifier tool located ");
+
+ // Get the SCT_ReadCaliDataSvc
+ if (m_useCalibData) {
+ // ATH_CHECK(m_ReadCalibChipDataTool.retrieve());
+ // ATH_MSG_DEBUG("CalibChipData Service located ");
+ ATH_MSG_FATAL("Use of calibration data not supported.");
+ } else {
+ // m_ReadCalibChipDataTool.disable();
+ }
+
+ // Get the maximum number of strips of any module
+ m_strip_max = m_SCTdetMgr->numerology().maxNumStrips();
+
+ constexpr float fC = 6242.2;
+ m_Threshold = m_Threshold * fC;
+
+#ifdef SCT_DIG_DEBUG
+ ATH_MSG_INFO("\tNoise factors:");
+ ATH_MSG_INFO("\tBarrel = " << m_NoiseBarrel << " Outer Barrel = " << m_NoiseBarrel3 <<
+ " EC, inners = " << m_NoiseInners << " EC, middles = " << m_NoiseMiddles <<
+ " EC, short middles = " << m_NoiseShortMiddles << " EC, outers = " << m_NoiseOuters);
+ ATH_MSG_INFO("\tThreshold=" << m_Threshold << " fC, time of Threshold=" << m_timeOfThreshold);
+#endif
+
+ ATH_MSG_INFO("m_Threshold (Threshold) = " << m_Threshold);
+ ATH_MSG_INFO("m_timeOfThreshold (TimeOfThreshold) = " << m_timeOfThreshold);
+ ATH_MSG_INFO("m_data_compression_mode (DataCompressionMode) = " << m_data_compression_mode);
+ ATH_MSG_INFO("m_data_readout_mode (DataReadOutMode) = " << m_data_readout_mode);
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode SCT_FrontEnd::finalize() {
+#ifdef SCT_DIG_DEBUG
+ ATH_MSG_INFO("SCT_FrontEnd::finalize()");
+#endif
+ return StatusCode::SUCCESS;
+}
+
+// ----------------------------------------------------------------------
+// Init the class variable vectors
+// ----------------------------------------------------------------------
+StatusCode SCT_FrontEnd::initVectors(int strips, SCT_FrontEndData& data) const {
+ data.m_Offset.assign(strips, 0.0);
+ data.m_GainFactor.assign(strips, 0.0);
+ data.m_NoiseFactor.assign(strips, 0.0);
+
+ if (m_data_compression_mode == 1 and m_data_readout_mode == 0) {
+ data.m_Analogue[0].reserve(1);
+ data.m_Analogue[1].reserve(strips);
+ } else if (m_data_compression_mode == 2 and m_data_readout_mode == 0) {
+ data.m_Analogue[0].reserve(strips);
+ data.m_Analogue[1].reserve(strips);
+ } else if (m_data_compression_mode == 3 or m_data_readout_mode == 1) {
+ data.m_Analogue[0].reserve(strips);
+ data.m_Analogue[1].reserve(strips);
+ data.m_Analogue[2].reserve(strips);
+ } else {
+ return StatusCode::FAILURE;
+ }
+
+ return StatusCode::SUCCESS;
+}
+
+// ----------------------------------------------------------------------
+// prepare gain and offset for the strips for a given module
+// ----------------------------------------------------------------------
+StatusCode SCT_FrontEnd::prepareGainAndOffset(SiChargedDiodeCollection& collection, /*const Identifier& moduleId,*/ CLHEP::HepRandomEngine * rndmEngine, SCT_FrontEndData& data) const {
+ // now we need to generate gain and offset channel by channel: some algebra
+ // for generation of partially correlated random numbers
+ float W = m_OGcorr * m_GainRMS * m_Ospread / (m_GainRMS * m_GainRMS - m_Ospread * m_Ospread);
+ float A = 4 * W * W + 1.0;
+ float x1 = (A - sqrt(A)) / (2.0 * A);
+ float sinfi = sqrt(x1);
+ float cosfi = sqrt(1.0 - x1);
+
+ sinfi = sinfi * m_OGcorr / fabs(m_OGcorr);
+ float S = m_GainRMS * m_GainRMS + m_Ospread * m_Ospread;
+ float D = (m_GainRMS * m_GainRMS - m_Ospread * m_Ospread) / (cosfi * cosfi - sinfi * sinfi);
+ float S1 = sqrt((S + D) * 0.5);
+ float S2 = sqrt((S - D) * 0.5);
+ float Noise = 0;
+ int mode = 1;
+
+ // To set noise values for different module types, barrel, EC, inners, middles, short middles, and outers
+ if (m_analogueNoiseOn) {
+ // if (m_sct_id->barrel_ec(moduleId) == 0) { // barrel_ec=0 corresponds to barrel
+ // if (m_sct_id->layer_disk(moduleId) == 3) { // outer barrel layer 10 degrees warmer
+ Noise = m_NoiseBarrel3;
+ // } else {
+ // Noise = m_NoiseBarrel;
+ // }
+ // } else {
+ // int moduleType = m_sct_id->eta_module(moduleId);
+ // switch (moduleType) { // eta = 0, 1, or 2 corresponds to outers, middles and inners?! (at least in the offline world)
+ // case 0: {
+ // Noise = m_NoiseOuters;
+ // break;
+ // }
+ // case 1: {
+ // if (m_sct_id->layer_disk(moduleId) == 7) {
+ // Noise = m_NoiseShortMiddles;
+ // } else {
+ // Noise = m_NoiseMiddles;
+ // }
+ // break;
+ // }
+ // case 2: {
+ // Noise = m_NoiseInners;
+ // break;
+ // }
+ // default: {
+ // Noise = m_NoiseBarrel;
+ // ATH_MSG_ERROR("moduleType(eta): " << moduleType << " unknown, using barrel");
+ // }
+ // }// end of switch structure
+ // }
+ }
+
+ // Loop over collection and setup gain/offset/noise for the hit and neighbouring strips
+ SiChargedDiodeIterator i_chargedDiode = collection.begin();
+ SiChargedDiodeIterator i_chargedDiode_end = collection.end();
+
+ for (; i_chargedDiode != i_chargedDiode_end; ++i_chargedDiode) {
+ SiChargedDiode diode = (*i_chargedDiode).second;
+ // should be const as we aren't trying to change it here - but getReadoutCell() is not a const method...
+ unsigned int flagmask = diode.flag() & 0xFE;
+ // Get the flag for this diode ( if flagmask = 1 If diode is disconnected/disabled skip it)
+ if (!flagmask) { // If the diode is OK (not flagged)
+ const SiReadoutCellId roCell = diode.getReadoutCell();
+ if (roCell.isValid()) {
+ int strip = roCell.strip();
+ int i = std::max(strip - 1, 0);
+ int i_end = std::min(strip + 2, m_strip_max.load());
+
+ // loop over strips
+ for (; i < i_end; i++) {
+ // Need to check if strip is already setup
+ if (data.m_Analogue[1][i] <= 0.0) {
+ float g = CLHEP::RandGaussZiggurat::shoot(rndmEngine, 0.0, S1);
+ float o = CLHEP::RandGaussZiggurat::shoot(rndmEngine, 0.0, S2);
+
+ data.m_GainFactor[i] = 1.0 + (cosfi * g + sinfi * o);
+ data.m_Offset[i] = (cosfi * o - sinfi * g);
+ data.m_NoiseFactor[i] = Noise * mode;
+
+ // Fill the noise and offset values into the Analogue
+ if (m_data_compression_mode == 1 and m_data_readout_mode == 0) { // level mode x1x
+ data.m_Analogue[1][i] = data.m_Offset[i] + data.m_NoiseFactor[i] * CLHEP::RandGaussZiggurat::shoot(rndmEngine);
+ } else if (m_data_compression_mode == 2 and m_data_readout_mode == 0) { // edge mode 01x
+ data.m_Analogue[0][i] = data.m_Offset[i] + data.m_NoiseFactor[i] * CLHEP::RandGaussZiggurat::shoot(rndmEngine);
+ data.m_Analogue[1][i] = data.m_Offset[i] + data.m_NoiseFactor[i] * CLHEP::RandGaussZiggurat::shoot(rndmEngine);
+ } else if (m_data_compression_mode == 3 or m_data_readout_mode == 1) { // any hit mode xxx or expanded read out mode
+ data.m_Analogue[0][i] = data.m_Offset[i] + data.m_NoiseFactor[i] * CLHEP::RandGaussZiggurat::shoot(rndmEngine);
+ data.m_Analogue[1][i] = data.m_Offset[i] + data.m_NoiseFactor[i] * CLHEP::RandGaussZiggurat::shoot(rndmEngine);
+ data.m_Analogue[2][i] = data.m_Offset[i] + data.m_NoiseFactor[i] * CLHEP::RandGaussZiggurat::shoot(rndmEngine);
+ }
+ }
+ }
+ }
+ }
+ }
+
+ return StatusCode::SUCCESS;
+}
+
+// ----------------------------------------------------------------------
+// prepare gain and offset for the strips for a given module using
+// Cond Db data to get the chip calibration data
+// ----------------------------------------------------------------------
+// StatusCode SCT_FrontEnd::prepareGainAndOffset(SiChargedDiodeCollection& collection, int side, const Identifier& moduleId, CLHEP::HepRandomEngine * rndmEngine, SCT_FrontEndData& data) const {
+// // Get chip data from calib DB
+// std::vector<float> gainByChipVect = m_ReadCalibChipDataTool->getNPtGainData(moduleId, side, "GainByChip");
+// std::vector<float> gainRMSByChipVect = m_ReadCalibChipDataTool->getNPtGainData(moduleId, side, "GainRMSByChip");
+// std::vector<float> offsetByChipVect = m_ReadCalibChipDataTool->getNPtGainData(moduleId, side, "OffsetByChip");
+// std::vector<float> offsetRMSByChipVect = m_ReadCalibChipDataTool->getNPtGainData(moduleId, side, "OffsetRMSByChip");
+// std::vector<float> noiseByChipVect(6, 0.0);
+
+// if (m_analogueNoiseOn) { // Check if noise should be on or off
+// noiseByChipVect = m_ReadCalibChipDataTool->getNPtGainData(moduleId, side, "NoiseByChip");
+// }
+
+// // Need to check if empty, most should have data, but a few old DEAD modules don't
+// if (gainByChipVect.empty() or noiseByChipVect.empty()) {
+// ATH_MSG_DEBUG("No calibration data in cond DB for module " << moduleId << " using JO values");
+// if (StatusCode::SUCCESS != prepareGainAndOffset(collection, moduleId, rndmEngine, data)) {
+// return StatusCode::FAILURE;
+// } else {
+// return StatusCode::SUCCESS;
+// }
+// }
+
+// // Don't really need to set up values for each chip...
+// float gainMeanValue = meanValue(gainByChipVect);
+// if (gainMeanValue < 0.0) {
+// ATH_MSG_DEBUG("All chip gain values are 0 for module " << moduleId << " using JO values");
+// if (StatusCode::SUCCESS != prepareGainAndOffset(collection, moduleId, rndmEngine, data)) {
+// return StatusCode::FAILURE;
+// } else {
+// return StatusCode::SUCCESS;
+// }
+// }
+
+// std::vector<float> gain(6, 0.0);
+// std::vector<float> offset(6, 0.0);
+// std::vector<float> S1(6, 0.0);
+// std::vector<float> S2(6, 0.0);
+// std::vector<float> sinfi(6, 0.0);
+// std::vector<float> cosfi(6, 0.0);
+// float gainRMS = 0.0;
+// float offsetRMS = 0.0;
+
+// for (int i = 0; i < 6; ++i) {
+// // Some very few chips have 0 values, dead/bypassed/etc, so check and use some fixed values instead
+// if (gainByChipVect[i] > 0.1) {
+// gain[i] = gainByChipVect[i] / gainMeanValue;
+// offset[i] = offsetByChipVect[i] / m_Threshold;
+// gainRMS = gainRMSByChipVect[i] / gainMeanValue;
+// offsetRMS = offsetRMSByChipVect[i] / m_Threshold;
+// } else {
+// gain[i] = 55.0 / gainMeanValue;
+// offset[i] = 42.0 / m_Threshold;
+// gainRMS = 1.3 / gainMeanValue;
+// offsetRMS = 2.0 / m_Threshold;
+// }
+
+// float W = m_OGcorr * gainRMS * offsetRMS / (gainRMS * gainRMS - offsetRMS * offsetRMS);
+// float A = 4 * W * W + 1.0;
+// float x1 = (A - sqrt(A)) / (2.0 * A);
+// sinfi[i] = sqrt(x1);
+// cosfi[i] = sqrt(1.0 - x1);
+// sinfi[i] = sinfi[i] * m_OGcorr / fabs(m_OGcorr);
+// float S = gainRMS * gainRMS + offsetRMS * offsetRMS;
+// float D = (gainRMS * gainRMS - offsetRMS * offsetRMS) / (cosfi[i] * cosfi[i] - sinfi[i] * sinfi[i]);
+// S1[i] = sqrt((S + D) / 2.0);
+// S2[i] = sqrt((S - D) / 2.0);
+// }
+
+// // Loop over collection and setup gain/offset/noise for the hit and neighbouring strips
+// SiChargedDiodeIterator i_chargedDiode = collection.begin();
+// SiChargedDiodeIterator i_chargedDiode_end = collection.end();
+
+// for (; i_chargedDiode != i_chargedDiode_end; ++i_chargedDiode) {
+// SiChargedDiode diode = (*i_chargedDiode).second;
+// // should be const as we aren't trying to change it here - but getReadoutCell() is not a const method...
+// unsigned int flagmask = diode.flag() & 0xFE;
+// // Get the flag for this diode ( if flagmask = 1 If diode is disconnected/disabled skip it)
+// if (!flagmask) { // If the diode is OK (not flagged)
+// const SiReadoutCellId roCell = diode.getReadoutCell();
+
+// if (roCell.isValid()) {
+// int strip = roCell.strip();
+// int i = std::max(strip - 1, 0);
+// int i_end = std::min(strip + 2, m_strip_max.load());
+
+// // loop over strips
+// for (; i < i_end; i++) {
+// // Need to check if strip is already setup
+// if (data.m_Analogue[1][i] <= 0.0) {
+// // Values depends on which chip the strip is on (complex when strip is on chip edge)
+// int chip = i / 128;
+// float g = CLHEP::RandGaussZiggurat::shoot(rndmEngine, 0.0, S1[chip]);
+// float o = CLHEP::RandGaussZiggurat::shoot(rndmEngine, 0.0, S2[chip]);
+
+// data.m_GainFactor[i] = gain[chip] + (cosfi[chip] * g + sinfi[chip] * o);
+// data.m_Offset[i] = offset[chip] + (cosfi[chip] * o - sinfi[chip] * g);
+// data.m_NoiseFactor[i] = noiseByChipVect[chip];
+
+// // Fill the noise and offset values into the Analogue
+// if (m_data_compression_mode == 1 and m_data_readout_mode == 0) { // level mode x1x
+// data.m_Analogue[1][i] = data.m_Offset[i] + data.m_NoiseFactor[i] * CLHEP::RandGaussZiggurat::shoot(rndmEngine);
+// } else if (m_data_compression_mode == 2 and m_data_readout_mode == 0) { // edge mode 01x
+// data.m_Analogue[0][i] = data.m_Offset[i] + data.m_NoiseFactor[i] * CLHEP::RandGaussZiggurat::shoot(rndmEngine);
+// data.m_Analogue[1][i] = data.m_Offset[i] + data.m_NoiseFactor[i] * CLHEP::RandGaussZiggurat::shoot(rndmEngine);
+// } else if (m_data_compression_mode == 3 or m_data_readout_mode == 1) { // any hit mode xxx or expanded read out mode
+// data.m_Analogue[0][i] = data.m_Offset[i] + data.m_NoiseFactor[i] * CLHEP::RandGaussZiggurat::shoot(rndmEngine);
+// data.m_Analogue[1][i] = data.m_Offset[i] + data.m_NoiseFactor[i] * CLHEP::RandGaussZiggurat::shoot(rndmEngine);
+// data.m_Analogue[2][i] = data.m_Offset[i] + data.m_NoiseFactor[i] * CLHEP::RandGaussZiggurat::shoot(rndmEngine);
+// }
+// }
+// }
+// }
+// }
+// }
+
+// return StatusCode::SUCCESS;
+// }
+
+// ----------------------------------------------------------------------
+//
+// ----------------------------------------------------------------------
+StatusCode SCT_FrontEnd::randomNoise(SiChargedDiodeCollection& collection, /*const Identifier& moduleId,*/ CLHEP::HepRandomEngine * rndmEngine, SCT_FrontEndData& data) const {
+ // Add random noise
+
+ double occupancy = 0.0;
+ double NoiseOccupancy = 0.0;
+ float Noise = 0.0;
+ int nNoisyStrips = 0;
+ double mode = 1.;
+
+ const bool noise_expanded_mode = (m_data_compression_mode == 3 and m_data_readout_mode == 1);
+
+ // Will give 3 times as much noise occupancy if running in any hit expanded mode
+ if (noise_expanded_mode) {
+ mode = 3.;
+ }
+
+ // Sets fixed noise occupancy values for different module types, barrel, EC,
+ // inners, middles
+ // short middles, and outers
+// if (m_sct_id->barrel_ec(moduleId) == 0) { // barrel_ec=0 corresponds to barrel
+// if (m_sct_id->layer_disk(moduleId) == 3) { // outer barrel layer 10 degrees warmer
+ NoiseOccupancy = m_NOBarrel3;
+ Noise = m_NoiseBarrel3;
+// } else {
+// NoiseOccupancy = m_NOBarrel;
+// Noise = m_NoiseBarrel;
+// }
+// } else {
+// int moduleType = m_sct_id->eta_module(moduleId);
+// switch (moduleType) { // eta = 0, 1, or 2 corresponds to outers, middles and inners?! (at least in the offline world)
+// case 0: {
+// NoiseOccupancy = m_NOOuters;
+// Noise = m_NoiseOuters;
+// break;
+// }
+// case 1: {
+// if (m_sct_id->layer_disk(moduleId) == 7) {
+// NoiseOccupancy = m_NOShortMiddles;
+// Noise = m_NoiseShortMiddles;
+// } else {
+// NoiseOccupancy = m_NOMiddles;
+// Noise = m_NoiseMiddles;
+// }
+// break;
+// }
+// case 2: {
+// NoiseOccupancy = m_NOInners;
+// Noise = m_NoiseInners;
+// break;
+// }
+// default: {
+// NoiseOccupancy = m_NOBarrel;
+// Noise = m_NoiseBarrel;
+// ATH_MSG_ERROR("moduleType(eta): " << moduleType << " unknown, using barrel");
+// }
+// }// end of switch structure
+// }
+
+ // Calculate the number of "free strips"
+ int nEmptyStrips = 0;
+ std::vector<int> emptyStrips;
+ emptyStrips.reserve(m_strip_max);
+ for (int i = 0; i < m_strip_max; i++) {
+ if (data.m_StripHitsOnWafer[i] == 0) {
+ emptyStrips.push_back(i);
+ ++nEmptyStrips;
+ }
+ }
+
+ if (nEmptyStrips != 0) {
+ // Should randomize the fixed NO values, so we get some differences per
+ // wafer
+ occupancy = CLHEP::RandGaussZiggurat::shoot(rndmEngine, NoiseOccupancy, NoiseOccupancy * 0.1);
+
+ // Modify the occupancy if threshold is not 1.0 fC
+ if (m_Threshold > 6242.3 or m_Threshold < 6242.1) {
+ const float fC = 6242.2;
+ occupancy = occupancy * exp(-(0.5 / (Noise * Noise) * (m_Threshold * m_Threshold - fC * fC)));
+ }
+ nNoisyStrips = CLHEP::RandPoisson::shoot(rndmEngine, m_strip_max * occupancy * mode);
+
+ // Check and adapt the number of noisy strips to the number of free strips
+ if (nEmptyStrips < nNoisyStrips) {
+ nNoisyStrips = nEmptyStrips;
+ }
+
+ // Find random strips to get noise hits
+ for (int i = 0; i < nNoisyStrips; i++) {
+ int index = CLHEP::RandFlat::shootInt(rndmEngine, nEmptyStrips - i); // strip == 10, 12 free strips
+ // have vector [10000100100200211001] 20 strips
+ int strip = emptyStrips.at(index);
+ emptyStrips.erase(emptyStrips.begin()+index); // Erase it not to use it again
+ if (data.m_StripHitsOnWafer[strip]!=0) {
+ ATH_MSG_ERROR(index << "-th empty strip, strip " << strip << " should be empty but is not empty! Something is wrong!");
+ }
+ data.m_StripHitsOnWafer[strip] = 3; // !< Random Noise hit
+ // Add tbin info to noise diode
+ if (noise_expanded_mode) { // !< if any hit mode, any time bin otherwise fixed tbin=2
+ int noise_tbin = CLHEP::RandFlat::shootInt(rndmEngine, 3);
+ // !< random number 0, 1 or 2
+ if (noise_tbin == 0) {
+ noise_tbin = 4; // !< now 1,2 or 4
+ }
+ if (StatusCode::SUCCESS != addNoiseDiode(collection, strip, noise_tbin)) {
+ ATH_MSG_ERROR("Can't add noise hit diode to collection");
+ }
+ } else {
+ if (StatusCode::SUCCESS != addNoiseDiode(collection, strip, 2)) {
+ ATH_MSG_ERROR("Can't add noise hit diode to collection");
+ }
+ }
+ }
+ }
+
+ return StatusCode::SUCCESS;
+}
+
+// ----------------------------------------------------------------------
+//
+// ----------------------------------------------------------------------
+// StatusCode SCT_FrontEnd::randomNoise(SiChargedDiodeCollection& collection, const Identifier& moduleId, int side, CLHEP::HepRandomEngine * rndmEngine, SCT_FrontEndData& data) const {
+// const int n_chips = 6;
+// const int chipStripmax = m_strip_max / n_chips;
+// std::vector<float> NOByChipVect(n_chips, 0.0);
+// std::vector<float> ENCByChipVect(n_chips, 0.0);
+// std::vector<int> nNoisyStrips(n_chips, 0);
+// double mode = 1.;
+
+// const bool noise_expanded_mode = (m_data_compression_mode == 3 and m_data_readout_mode == 1);
+
+// // Will give 3 times as much noise occupancy if running in any hit expanded mode
+// if (noise_expanded_mode) {
+// mode = 3.;
+// }
+
+// // Get chip data from calib DB
+// NOByChipVect = m_ReadCalibChipDataTool->getNoiseOccupancyData(moduleId, side, "OccupancyByChip");
+// ENCByChipVect = m_ReadCalibChipDataTool->getNPtGainData(moduleId, side, "NoiseByChip");
+
+// // Need to check if empty, most should have data, but a few old DEAD modules don't, and 9C...
+// if (NOByChipVect.empty()) {
+// ATH_MSG_DEBUG("No calibration data in cond DB for module " << moduleId << " using JO values");
+// if (StatusCode::SUCCESS != randomNoise(collection, moduleId, rndmEngine, data)) {
+// return StatusCode::FAILURE;
+// } else {
+// return StatusCode::SUCCESS;
+// }
+// } else {
+// for (int i = 0; i < n_chips; i++) {
+// // A 0 value can mean two things now, chip out of config for long time and no value was uploaded
+// // or its short middles and inners and the value is for all purposes 0! so ok.
+
+// // Modify the occupancy if threshold is not 1.0 fC
+// if (m_Threshold > 6242.3 or m_Threshold < 6242.1) {
+// constexpr float fC = 6242.2;
+// NOByChipVect[i] = NOByChipVect[i] * exp(-(0.5 / (ENCByChipVect[i]*ENCByChipVect[i]) * (m_Threshold*m_Threshold - fC*fC)));
+// }
+
+// nNoisyStrips[i] = CLHEP::RandPoisson::shoot(rndmEngine, chipStripmax * NOByChipVect[i] * mode);
+// }
+// }
+
+// // Loop over the chips on the wafer
+// for (int chip_index = 0; chip_index < n_chips; ++chip_index) {
+// int chip_strip_offset = chipStripmax * chip_index; // First strip number on chip
+
+// // Calculate the number of "free strips" on this chip
+// int nEmptyStripsOnChip = 0;
+// std::vector<int> emptyStripsOnChip;
+// emptyStripsOnChip.reserve(chipStripmax);
+// for (int i = 0; i < chipStripmax; i++) {
+// if (data.m_StripHitsOnWafer[i + chip_strip_offset] == 0) {
+// emptyStripsOnChip.push_back(i);
+// ++nEmptyStripsOnChip;
+// }
+// }
+
+// // if no empty strips on chip do nothing
+// if (nEmptyStripsOnChip != 0) {
+// // Check and adapt the number of noisy strips to the number of free strips
+// if (nEmptyStripsOnChip < nNoisyStrips[chip_index]) {
+// nNoisyStrips[chip_index] = nEmptyStripsOnChip;
+// }
+
+// // Find random strips to get noise hits
+// for (int i = 0; i < nNoisyStrips[chip_index]; i++) {
+// int index = CLHEP::RandFlat::shootInt(rndmEngine, nEmptyStripsOnChip - i);
+// int strip_on_chip = emptyStripsOnChip.at(index);
+// emptyStripsOnChip.erase(emptyStripsOnChip.begin()+index); // Erase it not to use it again
+// int strip = strip_on_chip + chip_strip_offset;
+// if (data.m_StripHitsOnWafer[strip]!=0) {
+// ATH_MSG_ERROR(index << "-th empty strip, strip " << strip << " should be empty but is not empty! Something is wrong!");
+// }
+// data.m_StripHitsOnWafer[strip] = 3; // !< Random Noise hit
+// // Add tbin info to noise diode
+// if (noise_expanded_mode) { // !< if any hit mode, any time bin
+// // !< otherwise fixed tbin=2
+// int noise_tbin = CLHEP::RandFlat::shootInt(rndmEngine, 3);
+// // !< random number 0, 1 or 2
+// if (noise_tbin == 0) {
+// noise_tbin = 4; // !< now 1, 2 or 4
+// }
+// if (StatusCode::SUCCESS != addNoiseDiode(collection, strip, noise_tbin)) {
+// ATH_MSG_ERROR("Can't add noise hit diode to collection");
+// }
+// } else {
+// if (StatusCode::SUCCESS != addNoiseDiode(collection, strip, 2)) {
+// ATH_MSG_ERROR("Can't add noise hit diode to collection");
+// }
+// }
+// }
+// }
+// }
+
+// return StatusCode::SUCCESS;
+// }
+
+// ----------------------------------------------------------------------
+// process the collection of pre digits this will need to go through
+// all single-strip pre-digits calculate the amplifier response add noise
+// (this could be moved elsewhere later) apply threshold do clustering
+// ----------------------------------------------------------------------
+void SCT_FrontEnd::process(SiChargedDiodeCollection& collection, CLHEP::HepRandomEngine * rndmEngine) const {
+ // get SCT module side design and check it
+ const SCT_ModuleSideDesign *p_design = dynamic_cast<const SCT_ModuleSideDesign*>(&(collection.design()));
+
+ if (p_design==nullptr) {
+ return;
+ }
+
+ SCT_FrontEndData data;
+
+ // Check number of strips in design and from manager(max number of strips on any module)
+ // The design value should always be equal or lower than the manager one
+ // However, no resising is now done in case of a lower value
+ const int strip_max = p_design->cells();
+ // Init vectors
+ if (StatusCode::SUCCESS != initVectors(strip_max, data)) {
+ ATH_MSG_ERROR("Can't resize class variable vectors");
+ return;
+ }
+
+ // Contains strip hit info, reset to 0 for each wafer processed
+ data.m_StripHitsOnWafer.assign(m_strip_max, 0);
+
+ // Containes the charge for each bin on each hit strip
+ if (m_data_compression_mode == 1 and m_data_readout_mode == 0) {
+ for (int i = 0; i < m_strip_max; ++i) {
+ data.m_Analogue[1][i] = 0.0;
+ }
+ } else if (m_data_compression_mode == 2 and m_data_readout_mode == 0) {
+ for (int i = 0; i < m_strip_max; ++i) {
+ data.m_Analogue[0][i] = 0.0;
+ data.m_Analogue[1][i] = 0.0;
+ }
+ } else if (m_data_compression_mode == 3 or m_data_readout_mode == 1) {
+ for (int i = 0; i < m_strip_max; ++i) {
+ data.m_Analogue[0][i] = 0.0;
+ data.m_Analogue[1][i] = 0.0;
+ data.m_Analogue[2][i] = 0.0;
+ }
+ }
+
+ // Get wafer, moduleId and side
+// Identifier waferId = collection.identify();
+// Identifier moduleId = m_sct_id->module_id(waferId);
+// const int side = m_sct_id->side(waferId);
+
+ // Check if collection empty
+ if (not collection.empty()) {
+ // Setup gain/offset/noise to the hit and neighbouring strips
+ if (m_useCalibData) { // Use calib cond DB data
+ ATH_MSG_FATAL("Use of calibration data not supported.");
+ // if (StatusCode::SUCCESS != prepareGainAndOffset(collection, side, moduleId, rndmEngine, data)) {
+ // ATH_MSG_ERROR("\tCan't prepare Gain and Offset");
+ // }
+ } else { // Use JO values
+ if (StatusCode::SUCCESS != prepareGainAndOffset(collection, /*moduleId,*/ rndmEngine, data)) {
+ ATH_MSG_ERROR("\tCan't prepare Gain and Offset");
+ }
+ }
+
+ if (StatusCode::SUCCESS != doSignalChargeForHits(collection, data)) {
+ ATH_MSG_ERROR("\tCan't doSignalChargeForHits");
+ }
+
+ if (StatusCode::SUCCESS != doThresholdCheckForRealHits(collection, data)) {
+ ATH_MSG_ERROR("\tCan't doThresholdCheckForRealHits");
+ }
+
+ if (StatusCode::SUCCESS != doThresholdCheckForCrosstalkHits(collection, data)) {
+ ATH_MSG_ERROR("\tCan't doThresholdCheckForCrosstalkHits");
+ }
+ }
+
+ if (m_NoiseOn) {
+ if (m_useCalibData) { // Check if using DB or not
+ ATH_MSG_FATAL("Use of calibration data not supported.");
+ // if (StatusCode::SUCCESS != randomNoise(collection, moduleId, side, rndmEngine, data)) {
+ // ATH_MSG_ERROR("\tCan't do random noise on wafer?!");
+ // }
+ } else { // Use JO fixed values
+ if (StatusCode::SUCCESS != randomNoise(collection, /*moduleId,*/ rndmEngine, data)) {
+ ATH_MSG_ERROR("\tCan't do random noise on wafer?!");
+ }
+ }
+ }
+
+ // Check for strips above threshold and do clustering
+ if (StatusCode::SUCCESS != doClustering(collection, data)) {
+ ATH_MSG_ERROR("\tCan't cluster the hits?!");
+ }
+}
+
+StatusCode SCT_FrontEnd::doSignalChargeForHits(SiChargedDiodeCollection& collection, SCT_FrontEndData& data) const {
+ typedef SiTotalCharge::list_t list_t;
+
+ // *****************************************************************************
+ // Loop over the diodes (strips ) and for each of them define the total signal
+ // *****************************************************************************
+
+ // set up number of needed bins depending on the compression mode
+ short bin_max = 0;
+ if (m_data_readout_mode == 0) {
+ bin_max = m_data_compression_mode;
+ } else {
+ bin_max = 3;
+ }
+
+ std::vector<float> response(bin_max);
+
+ SiChargedDiodeIterator i_chargedDiode = collection.begin();
+ SiChargedDiodeIterator i_chargedDiode_end = collection.end();
+ for (; i_chargedDiode != i_chargedDiode_end; ++i_chargedDiode) {
+ SiChargedDiode diode = (*i_chargedDiode).second;
+ // should be const as we aren't trying to change it here - but getReadoutCell() is not a const method...
+ unsigned int flagmask = diode.flag() & 0xFE;
+ // Get the flag for this diode ( if flagmask = 1 If diode is disconnected/disabled skip it)
+ if (!flagmask) { // If the diode is OK (not flagged)
+ const SiReadoutCellId roCell = diode.getReadoutCell();
+
+ if (roCell.isValid()) {
+ int strip = roCell.strip();
+
+ const list_t &ChargesOnStrip = diode.totalCharge().chargeComposition();
+
+ if (m_data_compression_mode == 1 and m_data_readout_mode == 0) { // level mode x1x
+ // Amplifier response
+ data.m_Analogue[1][strip] += data.m_GainFactor[strip] * m_sct_amplifier->response(ChargesOnStrip, m_timeOfThreshold);
+
+ // Add Crosstalk signal for neighboring strip
+ response[0] = m_sct_amplifier->crosstalk(ChargesOnStrip, m_timeOfThreshold);
+ if (strip + 1 < m_strip_max) {
+ data.m_Analogue[1][strip + 1] += data.m_GainFactor[strip + 1] * response[0];
+ }
+ if (strip > 0) {
+ data.m_Analogue[1][strip - 1] += data.m_GainFactor[strip - 1] * response[0];
+ }
+ } else if (m_data_compression_mode == 2 and m_data_readout_mode == 0) { // edge mode 01x
+ // Amplifier response
+ m_sct_amplifier->response(ChargesOnStrip, m_timeOfThreshold, response);
+ for (short bin = 0; bin < bin_max; ++bin) {
+ data.m_Analogue[bin][strip] += data.m_GainFactor[strip] * response[bin];
+ }
+ // Add Crosstalk signal for neighboring strip
+ m_sct_amplifier->crosstalk(ChargesOnStrip, m_timeOfThreshold, response);
+ for (short bin = 0; bin < bin_max; ++bin) {
+ if (strip + 1 < m_strip_max) {
+ data.m_Analogue[bin][strip + 1] += data.m_GainFactor[strip + 1] * response[bin];
+ }
+ if (strip > 0) {
+ data.m_Analogue[bin][strip - 1] += data.m_GainFactor[strip - 1] * response[bin];
+ }
+ }
+ } else if (m_data_compression_mode == 3 or m_data_readout_mode == 1) { // any hit mode xxx or expanded read out mode
+ // Amplifier response
+ m_sct_amplifier->response(ChargesOnStrip, m_timeOfThreshold, response);
+ for (short bin = 0; bin < bin_max; ++bin) {
+ data.m_Analogue[bin][strip] += data.m_GainFactor[strip] * response[bin];
+ }
+ // Add Crosstalk signal for neighboring strip
+ m_sct_amplifier->crosstalk(ChargesOnStrip, m_timeOfThreshold, response);
+ for (short bin = 0; bin < bin_max; ++bin) {
+ if (strip + 1 < m_strip_max) {
+ data.m_Analogue[bin][strip + 1] += data.m_GainFactor[strip + 1] * response[bin];
+ }
+ if (strip > 0) {
+ data.m_Analogue[bin][strip - 1] += data.m_GainFactor[strip - 1] * response[bin];
+ }
+ }
+ }
+ } else { // if roCell not valid
+ ATH_MSG_WARNING("\t Cannot get the cell ");
+ }
+ } else {// If diode is disconnected/disabled skip it
+ ATH_MSG_WARNING("\tDisabled or disconnected diode (strip)");
+ }
+ }
+ return StatusCode::SUCCESS;
+}
+
+StatusCode SCT_FrontEnd::doThresholdCheckForRealHits(SiChargedDiodeCollection& collection, SCT_FrontEndData& data) const {
+ // **********************************************************************************
+ // Flag strips below threshold and flag the threshold check into data.m_StripHitsOnWafer
+ // **********************************************************************************
+
+ SiChargedDiodeIterator i_chargedDiode = collection.begin();
+ SiChargedDiodeIterator i_chargedDiode_end = collection.end();
+
+ for (; i_chargedDiode != i_chargedDiode_end; ++i_chargedDiode) {
+ SiChargedDiode& diode = (*i_chargedDiode).second;
+ SiReadoutCellId roCell = diode.getReadoutCell();
+ if (roCell.isValid()) {
+ int strip = roCell.strip();
+ if (strip > -1 and strip < m_strip_max) {
+ if (m_data_compression_mode == 1 and m_data_readout_mode == 0) { // level mode x1x
+ if (data.m_Analogue[1][strip] < m_Threshold) {
+ SiHelper::belowThreshold(diode, true); // Below strip diode signal threshold
+ data.m_StripHitsOnWafer[strip] = -1;
+ } else if (((0x10 & diode.flag()) == 0x10) or ((0x4 & diode.flag()) == 0x4)) {
+ // previously a crazy strip number could have screwed things up here.
+ data.m_StripHitsOnWafer[strip] = -1;
+ } else {
+ data.m_StripHitsOnWafer[strip] = 1;
+ SiHelper::SetTimeBin(diode, 2, &msg()); // set timebin info
+ }
+ } else if (m_data_compression_mode == 2 and m_data_readout_mode == 0) { // edge mode 01x
+ if ((data.m_Analogue[0][strip] >= m_Threshold or data.m_Analogue[1][strip] < m_Threshold)) {
+ SiHelper::belowThreshold(diode, true); // Below strip diode signal threshold
+ data.m_StripHitsOnWafer[strip] = -1;
+ } else if (((0x10 & diode.flag()) == 0x10) or ((0x4 & diode.flag()) == 0x4)) {
+ // previously a crazy strip number could have screwed things up here.
+ data.m_StripHitsOnWafer[strip] = -1;
+ } else {
+ data.m_StripHitsOnWafer[strip] = 1;
+ SiHelper::SetTimeBin(diode, 2, &msg()); // set timebin info
+ }
+ } else if (m_data_compression_mode == 3 or m_data_readout_mode == 1) { // Check hit pattern
+ int have_hit_bin = 0;
+ if (data.m_Analogue[0][strip] >= m_Threshold) {
+ have_hit_bin = 4;
+ }
+ if (data.m_Analogue[1][strip] >= m_Threshold) {
+ have_hit_bin += 2;
+ }
+ if (data.m_Analogue[2][strip] >= m_Threshold) {
+ have_hit_bin += 1;
+ }
+ if (((0x10 & diode.flag()) == 0x10) || ((0x4 & diode.flag()) == 0x4)) {
+ // previously a crazy strip number could have screwed things up here.
+ data.m_StripHitsOnWafer[strip] = -1;
+ } else if (m_data_compression_mode == 1) { // !< level and expanded mode
+ if (have_hit_bin == 2 or have_hit_bin == 3 or have_hit_bin == 6 or have_hit_bin == 7) {
+ data.m_StripHitsOnWafer[strip] = 1;
+ SiHelper::SetTimeBin(diode, have_hit_bin, &msg());
+ } else {
+ SiHelper::belowThreshold(diode, true); // Below strip diode signal threshold
+ data.m_StripHitsOnWafer[strip] = -1;
+ }
+ } else if (m_data_compression_mode == 2) { // !< edge and expanded mode
+ if (have_hit_bin == 2 or have_hit_bin == 3) {
+ data.m_StripHitsOnWafer[strip] = 1;
+ SiHelper::SetTimeBin(diode, have_hit_bin, &msg());
+ } else {
+ SiHelper::belowThreshold(diode, true); // Below strip diode signal threshold
+ data.m_StripHitsOnWafer[strip] = -1;
+ }
+ } else if (m_data_compression_mode == 3) { // !< any hit mode
+ if (have_hit_bin == 0) {
+ SiHelper::belowThreshold(diode, true); // Below strip diode signal threshold
+ data.m_StripHitsOnWafer[strip] = -1;
+ } else {
+ data.m_StripHitsOnWafer[strip] = 1;
+ if (m_data_readout_mode == 1) { // !< check for exp mode or not
+ SiHelper::SetTimeBin(diode, have_hit_bin, &msg());
+ } else {
+ SiHelper::SetTimeBin(diode, 2, &msg());
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ return StatusCode::SUCCESS;
+}
+
+// ----------------------------------------------------------------------
+//
+// ----------------------------------------------------------------------
+StatusCode SCT_FrontEnd::doThresholdCheckForCrosstalkHits(SiChargedDiodeCollection& collection, SCT_FrontEndData& data) const {
+ // Check for noise+crosstalk strips above threshold
+ // data.m_StripHitsOnWafer: real hits above threshold == 1 or below/disconnected
+ // == -1
+ // =0 for free strips or strips with charge to be checked (data.m_Analogue[1]!=0)
+ // Set 2 for crosstalk noise hits and -2 for below ones
+
+ for (int strip = 0; strip < m_strip_max; strip++) {
+ // Find strips with data.m_StripHitsOnWafer[strip] == 0
+ if (data.m_StripHitsOnWafer[strip] != 0) { // real hits already checked
+ continue;
+ }
+ if (data.m_Analogue[1][strip] > 0) { // Better way of doing this?!
+ // set data.m_StripHitsOnWafer to x in prepareGainAndOffset
+ if (m_data_compression_mode == 1 and m_data_readout_mode == 0) { // level mode x1x
+ if (data.m_Analogue[1][strip] < m_Threshold) {
+ data.m_StripHitsOnWafer[strip] = -2; // Below threshold
+ } else {
+ data.m_StripHitsOnWafer[strip] = 2; // Crosstalk+Noise hit
+ if (StatusCode::SUCCESS != addNoiseDiode(collection, strip, 2)) {
+ ATH_MSG_ERROR("Can't add noise hit diode to collection");
+ }
+ }
+ } else if (m_data_compression_mode == 2 and m_data_readout_mode == 0) { // edge mode 01x
+ if ((data.m_Analogue[0][strip] >= m_Threshold or data.m_Analogue[1][strip] < m_Threshold)) {
+ data.m_StripHitsOnWafer[strip] = -2; // Below threshold
+ } else {
+ data.m_StripHitsOnWafer[strip] = 2; // Crosstalk+Noise hit
+ if (StatusCode::SUCCESS != addNoiseDiode(collection, strip, 2)) {
+ ATH_MSG_ERROR("Can't add noise hit diode to collection");
+ }
+ }
+ } else if (m_data_compression_mode == 3 or m_data_readout_mode == 1) {
+ int have_hit_bin = 0;
+ if (data.m_Analogue[0][strip] >= m_Threshold) {
+ have_hit_bin = 4;
+ }
+ if (data.m_Analogue[1][strip] >= m_Threshold) {
+ have_hit_bin += 2;
+ }
+ if (data.m_Analogue[2][strip] >= m_Threshold) {
+ have_hit_bin += 1;
+ }
+ if (m_data_compression_mode == 1) { // !< level and expanded mode
+ if (have_hit_bin == 2 or have_hit_bin == 3 or have_hit_bin == 6 or have_hit_bin == 7) {
+ data.m_StripHitsOnWafer[strip] = 2; // Crosstalk+Noise hit
+ if (StatusCode::SUCCESS != addNoiseDiode(collection, strip, have_hit_bin)) {
+ ATH_MSG_ERROR("Can't add noise hit diode to collection");
+ }
+ } else {
+ data.m_StripHitsOnWafer[strip] = -2; // Below threshold
+ }
+ } else if (m_data_compression_mode == 2) { // !< edge and expanded mode
+ if (have_hit_bin == 2 or have_hit_bin == 3) {
+ data.m_StripHitsOnWafer[strip] = 2; // Noise hit
+ if (StatusCode::SUCCESS != addNoiseDiode(collection, strip, have_hit_bin)) {
+ ATH_MSG_ERROR("Can't add noise hit diode to collection");
+ }
+ } else {
+ data.m_StripHitsOnWafer[strip] = -2; // Below threshold
+ }
+ } else if (m_data_compression_mode == 3) { // !< any hit mode
+ if (have_hit_bin == 0) {
+ data.m_StripHitsOnWafer[strip] = -2; // Below threshold
+ } else {
+ data.m_StripHitsOnWafer[strip] = 2; // !< Crosstalk+Noise hit
+ if (m_data_readout_mode == 1) { // !< check for exp mode or not
+ if (StatusCode::SUCCESS != addNoiseDiode(collection, strip, have_hit_bin)) {
+ ATH_MSG_ERROR("Can't add noise hit diode to collection");
+ }
+ } else {
+ if (StatusCode::SUCCESS != addNoiseDiode(collection, strip, 2)) {
+ ATH_MSG_ERROR("Can't add noise hit diode to collection");
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode SCT_FrontEnd::doClustering(SiChargedDiodeCollection& collection, SCT_FrontEndData& data) const {
+ // ********************************
+ // now do clustering
+ // ********************************
+ int strip = 0;
+ int clusterSize = 0;
+
+ const SCT_ModuleSideDesign& sctDesign = dynamic_cast<const SCT_ModuleSideDesign&>(collection.design());
+
+ Identifier hitStrip;
+
+ if (m_data_readout_mode == 0) {
+ do {
+ if (data.m_StripHitsOnWafer[strip] > 0) {
+ // ====== First step: Get the cluster size
+ // ===================================================
+ int clusterFirstStrip = strip;
+
+ // Find end of cluster. In multi-row sensors, cluster cannot span rows.
+ int row = sctDesign.row(strip);
+ if (row < 0) {
+ row = 0;
+ }
+
+ int lastStrip1DInRow = 0;
+ for (int i = 0; i < row + 1; ++i) {
+ lastStrip1DInRow += sctDesign.diodesInRow(i);
+ }
+
+ while (strip < lastStrip1DInRow-1 and data.m_StripHitsOnWafer[strip +1] > 0) {
+ ++strip; // !< Find first strip hit and add up the following strips
+ }
+ int clusterLastStrip = strip;
+
+ clusterSize = (clusterLastStrip - clusterFirstStrip) + 1;
+ hitStrip = m_sct_id->strip_id(collection.identify(), clusterFirstStrip);
+ SiChargedDiode& HitDiode = *(collection.find(hitStrip));
+ SiHelper::SetStripNum(HitDiode, clusterSize, &msg());
+
+ SiChargedDiode *PreviousHitDiode = &HitDiode;
+ for (int i = clusterFirstStrip+1; i <= clusterLastStrip; ++i) {
+ hitStrip = m_sct_id->strip_id(collection.identify(), i);
+ SiChargedDiode& HitDiode2 = *(collection.find(hitStrip));
+ SiHelper::ClusterUsed(HitDiode2, true);
+ PreviousHitDiode->setNextInCluster(&HitDiode2);
+ PreviousHitDiode = &HitDiode2;
+ }
+ }
+ ++strip; // !< This is the starting point of the next cluster within this collection
+ } while (strip < m_strip_max);
+ } else {
+ // Expanded read out mode, basically one RDO/strip per cluster
+ // But if consecutively fired strips have the same time bin, those are converted into one cluster.
+ do {
+ clusterSize = 1;
+ if (data.m_StripHitsOnWafer[strip] > 0) {
+ hitStrip = m_sct_id->strip_id(collection.identify(), strip);
+ SiChargedDiode& hitDiode = *(collection.find(hitStrip));
+ int timeBin = SiHelper::GetTimeBin(hitDiode);
+ SiChargedDiode* previousHitDiode = &hitDiode;
+ // Check if consecutively fired strips have the same time bin
+ for (int newStrip=strip+1; newStrip<m_strip_max; newStrip++) {
+ if (not (data.m_StripHitsOnWafer[newStrip]>0)) break;
+ Identifier newHitStrip = m_sct_id->strip_id(collection.identify(), newStrip);
+ SiChargedDiode& newHitDiode = *(collection.find(newHitStrip));
+ if (timeBin!=SiHelper::GetTimeBin(newHitDiode)) break;
+ SiHelper::ClusterUsed(newHitDiode, true);
+ previousHitDiode->setNextInCluster(&newHitDiode);
+ previousHitDiode = &newHitDiode;
+ clusterSize++;
+ }
+ SiHelper::SetStripNum(hitDiode, clusterSize, &msg());
+
+#ifdef SCT_DIG_DEBUG
+ ATH_MSG_DEBUG("RDO Strip = " << strip << ", tbin = " <<
+ SiHelper::GetTimeBin(hitDiode) <<
+ ", HitInfo(1=real, 2=crosstalk, 3=noise): " <<
+ data.m_StripHitsOnWafer[strip]);
+#endif
+ }
+ strip += clusterSize; // If more than one strip fires, those following strips are skipped.
+ } while (strip < m_strip_max);
+ }
+
+ // clusters below threshold, only from pre-digits that existed before no
+ // pure noise clusters below threshold will be made
+ // D. Costanzo. I don't see why we should cluster the strips below
+ // threshold. I'll pass on the info of strips below threshold
+ // to the SDOs. I'll leave the SCT experts the joy to change this if they
+ // don't like what I did or if this requires too much memory/disk
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode SCT_FrontEnd::addNoiseDiode(SiChargedDiodeCollection& collection, int strip, int tbin) const {
+ const SiCellId ndiode(strip); // !< create a new diode
+ const SiCharge noiseCharge(2 * m_Threshold, 0, SiCharge::noise); // !< add some noise to it
+ collection.add(ndiode, noiseCharge); // !< add it to the collection
+
+ // Get the strip back to check
+ const Identifier idstrip = m_sct_id->strip_id(collection.identify(), strip);
+ SiChargedDiode *NoiseDiode = (collection.find(idstrip));
+ if (NoiseDiode == nullptr) {
+ return StatusCode::FAILURE;
+ }
+ SiHelper::SetTimeBin(*NoiseDiode, tbin, &msg()); // set timebin info
+
+ return StatusCode::SUCCESS;
+}
+
+float SCT_FrontEnd::meanValue(std::vector<float>& calibDataVect) const {
+ float mean_value = 0.0;
+ int nData = 0;
+ const unsigned int vec_size = calibDataVect.size();
+
+ for (unsigned int i = 0; i < vec_size; ++i) {
+ if (calibDataVect[i] > 0.1) {
+ mean_value += calibDataVect[i];
+ ++nData;
+ }
+ }
+
+ if (nData == 0) {
+ return -1;
+ } else {
+ return mean_value / nData;
+ }
+}
diff --git a/Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_FrontEnd.h b/Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_FrontEnd.h
new file mode 100644
index 00000000..03f34b2c
--- /dev/null
+++ b/Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_FrontEnd.h
@@ -0,0 +1,140 @@
+// -*- C++ -*-
+
+/*
+ Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+*/
+
+/**
+ * Header file for SCT_FrontEnd Class
+ * @brief simulation of the SCT front-end electronics
+ * working as a SiPreDigitsProcessor
+ * models response of ABCD chip amplifiers to
+ * collected charges, also does cross-talk, offset
+ * variation and gain variation, in a correlated way
+ * Version 1.0 04.05.2001 Szymon Gadomski
+ * 1.1 07.06.2001 pass amplifier and design as pointers
+ * 1.2, 15/06/2001, compiles with framework and tested
+ * 31.07.2001 changes, constructor now without design
+ * design is only known to process method (from the collection)
+ * 10.08.2001 dedided not to use AlwaysSame flag
+ * Implementation for using it may be provided later
+ * @author Szymon.Gadomski@cern.ch, Awatif.Belymam@cern.ch, Davide.Costanzo@cern.ch,
+ * tgcornel@nikhef.nl, Grant.Gorfine@cern.ch, Paul.Bell@cern.ch, Jorgen.Dalmau@cern.ch,
+ * 23/08/2007 - Kondo.Gnanvo@cern.ch
+ * - Conversion of the SCT_Front code Al gTool
+ */
+
+#ifndef FASERSCT_DIGITIZATION_SCT_FRONTEND_H
+#define FASERSCT_DIGITIZATION_SCT_FRONTEND_H
+
+// Inheritance
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "FaserSCT_Digitization/ISCT_FrontEnd.h"
+
+// Athena
+#include "FaserSiDigitization/SiChargedDiodeCollection.h"
+// #include "SCT_ConditionsTools/ISCT_ReadCalibChipDataTool.h"
+
+// Gaudi
+#include "GaudiKernel/ToolHandle.h"
+
+// STL
+#include <atomic>
+#include <mutex>
+#include <vector>
+
+class ISCT_Amp;
+class FaserSCT_ID;
+
+namespace TrackerDD {
+ class SCT_DetectorManager;
+}
+
+namespace CLHEP {
+ class HepRandomEngine;
+}
+/**
+ * @brief simulation of the SCT front-end electronics
+ * working as a SiPreDigitsProcessor
+ * models response of ABCD chip amplifiers to
+ * collected charges, also does cross-talk, offset
+ * variation and gain variation, in a correlated way
+ */
+
+struct SCT_FrontEndData {
+ std::vector<float> m_Offset; //!< generate offset per channel
+ std::vector<float> m_GainFactor; //!< generate gain per channel (added to the gain per chip from calib data)
+ std::vector<float> m_NoiseFactor; //!< Kondo: 31/08/07 noise per channel (actually noise per chip from calib data)
+ std::vector<double> m_Analogue[3]; //!< To hold the noise and amplifier response
+ std::vector<int> m_StripHitsOnWafer; //!< Info about which strips are above threshold
+};
+
+class SCT_FrontEnd : public extends<AthAlgTool, ISCT_FrontEnd> {
+
+ public:
+
+ /** constructor */
+ SCT_FrontEnd(const std::string& type, const std::string& name, const IInterface* parent);
+ /** Destructor */
+ virtual ~SCT_FrontEnd() = default;
+ //PJ not needed after merging?!
+ /** AlgTool InterfaceID */
+ //static const InterfaceID& interfaceID();
+
+ /** AlgTool initialize */
+ virtual StatusCode initialize() override;
+ /** AlgTool finalize */
+ virtual StatusCode finalize() override;
+
+ /**
+ * process the collection of pre digits: needed to go through all single-strip pre-digits to calculate
+ * the amplifier response add noise (this could be moved elsewhere later) apply threshold do clustering
+ */
+ virtual void process(SiChargedDiodeCollection& collection, CLHEP::HepRandomEngine* rndmEngine) const override;
+ StatusCode doSignalChargeForHits(SiChargedDiodeCollection& collectione, SCT_FrontEndData& data) const;
+ StatusCode doThresholdCheckForRealHits(SiChargedDiodeCollection& collectione, SCT_FrontEndData& data) const;
+ StatusCode doThresholdCheckForCrosstalkHits(SiChargedDiodeCollection& collection, SCT_FrontEndData& data) const;
+ StatusCode doClustering(SiChargedDiodeCollection& collection, SCT_FrontEndData& data) const;
+ StatusCode prepareGainAndOffset(SiChargedDiodeCollection& collection, /*const Identifier& moduleId,*/ CLHEP::HepRandomEngine* rndmEngine, SCT_FrontEndData& data) const;
+// StatusCode prepareGainAndOffset(SiChargedDiodeCollection& collection, int side, const Identifier& moduleId, CLHEP::HepRandomEngine* rndmEngine, SCT_FrontEndData& data) const;
+ StatusCode randomNoise(SiChargedDiodeCollection& collection, /*const Identifier& moduleId,*/ CLHEP::HepRandomEngine* rndmEngine, SCT_FrontEndData& data) const;
+// StatusCode randomNoise(SiChargedDiodeCollection& collection, const Identifier& moduleId, int side, CLHEP::HepRandomEngine* rndmEngine, SCT_FrontEndData& data) const;
+ StatusCode addNoiseDiode(SiChargedDiodeCollection& collection, int strip, int tbin) const;
+ float meanValue(std::vector<float>& calibDataVect) const;
+ StatusCode initVectors(int strips, SCT_FrontEndData& data) const;
+
+ private:
+
+ FloatProperty m_NoiseBarrel{this, "NoiseBarrel", 1500.0, "Noise factor, Barrel (in the case of no use of calibration data)"};
+ FloatProperty m_NoiseBarrel3{this, "NoiseBarrel3", 1541.0, "Noise factor, Barrel3 (in the case of no use of calibration data)"};
+ FloatProperty m_NoiseInners{this, "NoiseInners", 1090.0, "Noise factor, EC Inners (in the case of no use of calibration data)"};
+ FloatProperty m_NoiseMiddles{this, "NoiseMiddles", 1557.0, "Noise factor, EC Middles (in the case of no use of calibration data)"};
+ FloatProperty m_NoiseShortMiddles{this, "NoiseShortMiddles", 940.0, "Noise factor, EC Short Middles (in the case of no use of calibration data)"};
+ FloatProperty m_NoiseOuters{this, "NoiseOuters", 1618.0, "Noise factor, Ec Outers (in the case of no use of calibration data)"};
+ DoubleProperty m_NOBarrel{this, "NOBarrel", 1.5e-5, "Noise factor, Barrel (in the case of no use of calibration data)"};
+ DoubleProperty m_NOBarrel3{this, "NOBarrel3", 2.1e-5, "Noise factor, Barrel3 (in the case of no use of calibration data)"};
+ DoubleProperty m_NOInners{this, "NOInners", 5.0e-9, "Noise Occupancy, EC Inners (in the case of no use of calibration data)"};
+ DoubleProperty m_NOMiddles{this, "NOMiddles", 2.7e-5, "Noise Occupancy, EC Middles (in the case of no use of calibration data)"};
+ DoubleProperty m_NOShortMiddles{this, "NOShortMiddles", 2.0e-9, "Noise Occupancy, EC Short Middles (in the case of no use of calibration data)"};
+ DoubleProperty m_NOOuters{this, "NOOuters", 3.5e-5, "Noise Occupancy, Ec Outers (in the case of no use of calibration data)"};
+ BooleanProperty m_NoiseOn{this, "NoiseOn", true, "To know if noise is on or off when using calibration data"};
+ BooleanProperty m_analogueNoiseOn{this, "AnalogueNoiseOn", true, "To know if analogue noise is on or off"};
+ FloatProperty m_GainRMS{this, "GainRMS", 0.031, "Gain spread parameter within the strips for a given Chip gain"};
+ FloatProperty m_Ospread{this, "Ospread", 0.0001, "offset spread within the strips for a given Chip offset"};
+ FloatProperty m_OGcorr{this, "OffsetGainCorrelation", 0.00001, "Gain/offset correlation for the strips"};
+ FloatProperty m_Threshold{this, "Threshold", 1.0, "Threshold"};
+ FloatProperty m_timeOfThreshold{this, "TimeOfThreshold", 30.0, "Threshold time"};
+ ShortProperty m_data_compression_mode{this, "DataCompressionMode", 1, "Front End Data Compression Mode"};
+ ShortProperty m_data_readout_mode{this, "DataReadOutMode", 0, "Front End Data Read out mode Mode"};
+ BooleanProperty m_useCalibData{this, "UseCalibData", false, "Flag to set the use of calibration data for noise, Gain,offset etc."}; // was true in ATLAS
+
+ ToolHandle<ISCT_Amp> m_sct_amplifier{this, "SCT_Amp", "SCT_Amp", "Handle the Amplifier tool"}; //!< Handle the Amplifier tool
+// ToolHandle<ISCT_ReadCalibChipDataTool> m_ReadCalibChipDataTool{this, "SCT_ReadCalibChipDataTool", "SCT_ReadCalibChipDataTool", "Tool to retrieve chip calibration information"}; //!< Handle to the Calibration ConditionsTool
+
+ const TrackerDD::SCT_DetectorManager* m_SCTdetMgr{nullptr}; //!< Handle to SCT detector manager
+ const FaserSCT_ID* m_sct_id{nullptr}; //!< Handle to SCT ID helper
+
+ mutable std::atomic_int m_strip_max{768}; //!< For SLHC studies
+};
+
+#endif //SCT_FRONTEND_H
diff --git a/Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_RandomDisabledCellGenerator.cxx b/Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_RandomDisabledCellGenerator.cxx
new file mode 100644
index 00000000..b170bd10
--- /dev/null
+++ b/Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_RandomDisabledCellGenerator.cxx
@@ -0,0 +1,40 @@
+/*
+ Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "SCT_RandomDisabledCellGenerator.h"
+
+#include "FaserSiDigitization/SiChargedDiodeCollection.h"
+
+#include "FaserSiDigitization/SiHelper.h"
+
+#include "CLHEP/Random/RandomEngine.h"
+#include "CLHEP/Random/RandFlat.h"
+
+// constructor
+SCT_RandomDisabledCellGenerator::SCT_RandomDisabledCellGenerator(const std::string& type, const std::string& name, const IInterface* parent)
+ : base_class(type, name, parent)
+{
+}
+
+StatusCode SCT_RandomDisabledCellGenerator::initialize() {
+ ATH_MSG_DEBUG("SCT_RandomDisabledCellGenerator::initialize()");
+ ATH_MSG_INFO("\tCreating missing bond generator with "<<m_disableProbability<<" probability");
+ return StatusCode::SUCCESS;
+}
+
+StatusCode SCT_RandomDisabledCellGenerator::finalize() {
+ ATH_MSG_INFO("SCT_RandomDisabledCellGenerator::finalize()");
+ return StatusCode::SUCCESS;
+}
+
+// process the collection
+void SCT_RandomDisabledCellGenerator::process(SiChargedDiodeCollection& collection, CLHEP::HepRandomEngine * rndmEngine) const {
+ // disabling is applied to all cells even unconnected or below threshold ones to be able to use these cells as well
+ // loop on all charged diodes
+ for (std::pair<const TrackerDD::SiCellId, SiChargedDiode>& chargedDiode: collection) {
+ if (CLHEP::RandFlat::shoot(rndmEngine)<m_disableProbability) {
+ SiHelper::disconnected(chargedDiode.second, true, false);
+ }
+ }
+}
diff --git a/Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_RandomDisabledCellGenerator.h b/Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_RandomDisabledCellGenerator.h
new file mode 100644
index 00000000..1ea27123
--- /dev/null
+++ b/Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_RandomDisabledCellGenerator.h
@@ -0,0 +1,67 @@
+// -*- C++ -*-
+
+/*
+ Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// SCT_RandomDisabledCellGenerator.h
+// Header file for class SCT_RandomDisabledCellGenerator
+///////////////////////////////////////////////////////////////////
+// Class to randomly disable cells for SCT
+///////////////////////////////////////////////////////////////////
+// Version 1.0 25/07/2007 Kondo Gnanvo
+// - Randomly select a cell and modify its flag by using the pointer to
+// a SiHelper function.
+// - Inherits from ISiChargedDiodeProcessor, as SiPreDiodeProcessor
+// has been discontinued
+// This can be used for disabling, disconnecting, and flagging bad_tot
+///////////////////////////////////////////////////////////////////
+
+#ifndef FASERSCT_DIGITIZATION_SCT_RANDOMDISABLEDCELLGENERATOR_H
+#define FASERSCT_DIGITIZATION_SCT_RANDOMDISABLEDCELLGENERATOR_H
+
+//Inheritance
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "FaserSCT_Digitization/ISCT_RandomDisabledCellGenerator.h"
+
+class SiChargedDiode;
+class SiChargedDiodeCollection;
+
+namespace CLHEP {
+ class HepRandomEngine;
+}
+
+class SCT_RandomDisabledCellGenerator : public extends<AthAlgTool, ISCT_RandomDisabledCellGenerator> {
+
+ ///////////////////////////////////////////////////////////////////
+ // Public methods:
+ ///////////////////////////////////////////////////////////////////
+
+ public:
+
+ /** constructor */
+ SCT_RandomDisabledCellGenerator(const std::string& type, const std::string& name, const IInterface* parent) ;
+
+ /** Destructor */
+ virtual ~SCT_RandomDisabledCellGenerator() = default;
+
+ /** AlgTool initialize */
+ virtual StatusCode initialize() override;
+ /** AlgTool finalize */
+ virtual StatusCode finalize() override;
+
+ virtual void process(SiChargedDiodeCollection& collection, CLHEP::HepRandomEngine * rndmEngine) const override;
+
+ private:
+
+ ///////////////////////////////////////////////////////////////////
+ // Private data:
+ ///////////////////////////////////////////////////////////////////
+ private:
+
+ FloatProperty m_disableProbability{this, "TotalBadChannels", 0.01, "probability that a cell is disabled"};
+
+};
+
+#endif //FASERSCT_DIGITIZATION_SCT_RANDOMDISABLEDCELLGENERATOR_H
diff --git a/Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_SurfaceChargesGenerator.cxx b/Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_SurfaceChargesGenerator.cxx
new file mode 100644
index 00000000..b05cd7b4
--- /dev/null
+++ b/Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_SurfaceChargesGenerator.cxx
@@ -0,0 +1,578 @@
+/*
+ Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "SCT_SurfaceChargesGenerator.h"
+
+// DD
+#include "TrackerReadoutGeometry/SiDetectorElement.h"
+#include "TrackerReadoutGeometry/SCT_ModuleSideDesign.h"
+
+// Athena
+#include "GeneratorObjects/HepMcParticleLink.h"
+#include "TrackerSimEvent/FaserSiHit.h" // for SiHit, SiHit::::xDep, etc
+#include "HitManagement/TimedHitPtr.h" // for TimedHitPtr
+
+// ROOT
+#include "TH1.h" // for TH1F
+#include "TH2.h" // for TH2F
+#include "TProfile.h" // for TProfile
+
+// CLHEP
+#include "CLHEP/Geometry/Point3D.h"
+#include "CLHEP/Random/RandomEngine.h"
+#include "CLHEP/Random/RandGaussZiggurat.h"
+#include "CLHEP/Units/SystemOfUnits.h"
+
+// STL
+#include <cmath>
+
+using TrackerDD::SiDetectorElement;
+using TrackerDD::SCT_ModuleSideDesign;
+using TrackerDD::SiLocalPosition;
+
+using namespace std;
+
+// constructor
+SCT_SurfaceChargesGenerator::SCT_SurfaceChargesGenerator(const std::string& type,
+ const std::string& name,
+ const IInterface* parent)
+ : base_class(type, name, parent) {
+}
+
+// ----------------------------------------------------------------------
+// Initialize
+// ----------------------------------------------------------------------
+StatusCode SCT_SurfaceChargesGenerator::initialize() {
+ ATH_MSG_DEBUG("SCT_SurfaceChargesGenerator::initialize()");
+
+ // Get ISiPropertiesTool
+ ATH_CHECK(m_siPropertiesTool.retrieve());
+
+ // Get ISiliconConditionsSvc
+// ATH_CHECK(m_siConditionsTool.retrieve());
+
+ if (m_doTrapping) {
+ // -- Get Radiation Damage Tool
+ ATH_MSG_FATAL("Radiation damage tool not supported.");
+// ATH_CHECK(m_radDamageTool.retrieve());
+// } else {
+// m_radDamageTool.disable();
+ }
+
+// if (m_doHistoTrap) {
+// // -- Get Histogram Service
+// ATH_CHECK(m_thistSvc.retrieve());
+
+// m_h_efieldz = new TProfile("efieldz", "", 50, 0., 0.4);
+// ATH_CHECK(m_thistSvc->regHist("/file1/efieldz", m_h_efieldz));
+
+// m_h_efield = new TH1F("efield", "", 100, 200., 800.);
+// ATH_CHECK(m_thistSvc->regHist("/file1/efield", m_h_efield));
+
+// m_h_spess = new TH1F("spess", "", 50, 0., 0.4);
+// ATH_CHECK(m_thistSvc->regHist("/file1/spess", m_h_spess));
+
+// m_h_depD = new TH1F("depD", "", 50, -0.3, 0.3);
+// ATH_CHECK(m_thistSvc->regHist("/file1/depD", m_h_depD));
+
+// m_h_drift_electrode = new TH2F("drift_electrode", "", 50, 0., 20., 50, 0., 20.);
+// ATH_CHECK(m_thistSvc->regHist("/file1/drift_electrode", m_h_drift_electrode));
+
+// m_h_drift_time = new TH1F("drift_time", "", 100, 0., 20.);
+// ATH_CHECK(m_thistSvc->regHist("/file1/drift_time", m_h_drift_time));
+
+// m_h_t_electrode = new TH1F("t_electrode", "", 100, 0., 20.);
+// ATH_CHECK(m_thistSvc->regHist("/file1/t_electrode", m_h_t_electrode));
+
+// m_h_ztrap = new TH1F("ztrap", "", 100, 0., 0.3);
+// ATH_CHECK(m_thistSvc->regHist("/file1/ztrap", m_h_ztrap));
+
+// // More histograms to check what's going on
+// m_h_zhit = new TH1F("zhit", "", 50, -0.2, 0.2);
+// ATH_CHECK(m_thistSvc->regHist("/file1/zhit", m_h_zhit));
+
+// m_h_ztrap_tot = new TH1F("ztrap_tot", "", 100, 0., 0.5);
+// ATH_CHECK(m_thistSvc->regHist("/file1/ztrap_tot", m_h_ztrap_tot));
+
+// m_h_no_ztrap = new TH1F("no_ztrap", "", 100, 0., 0.5);
+// ATH_CHECK(m_thistSvc->regHist("/file1/no_ztrap", m_h_no_ztrap));
+
+// m_h_trap_drift_t = new TH1F("trap_drift_t", "", 100, 0., 20.);
+// ATH_CHECK(m_thistSvc->regHist("/file1/trap_drift_t", m_h_trap_drift_t));
+
+// m_h_notrap_drift_t = new TH1F("notrap_drift_t", "", 100, 0., 20.);
+// ATH_CHECK(m_thistSvc->regHist("/file1/notrap_drift_t", m_h_notrap_drift_t));
+
+// m_h_mob_Char = new TProfile("mob_Char", "", 200, 100., 1000.);
+// ATH_CHECK(m_thistSvc->regHist("/file1/mob_Char", m_h_mob_Char));
+
+// m_h_vel = new TProfile("vel", "", 100, 100., 1000.);
+// ATH_CHECK(m_thistSvc->regHist("/file1/vel", m_h_vel));
+
+// m_h_drift1 = new TProfile("drift1", "", 50, 0., 0.3);
+// ATH_CHECK(m_thistSvc->regHist("/file1/drift1", m_h_drift1));
+
+// m_h_gen = new TProfile("gen", "", 50, 0., 0.3);
+// ATH_CHECK(m_thistSvc->regHist("/file1/gen", m_h_gen));
+
+// m_h_gen1 = new TProfile("gen1", "", 50, 0., 0.3);
+// ATH_CHECK(m_thistSvc->regHist("/file1/gen1", m_h_gen1));
+
+// m_h_gen2 = new TProfile("gen2", "", 50, 0., 0.3);
+// ATH_CHECK(m_thistSvc->regHist("/file1/gen2", m_h_gen2));
+
+// m_h_velocity_trap = new TProfile("velocity_trap", "", 50, 0., 1000.);
+// ATH_CHECK(m_thistSvc->regHist("/file1/velocity_trap", m_h_velocity_trap));
+
+// m_h_mobility_trap = new TProfile("mobility_trap", "", 50, 100., 1000.);
+// ATH_CHECK(m_thistSvc->regHist("/file1/mobility_trap", m_h_mobility_trap));
+
+// m_h_trap_pos = new TH1F("trap_pos", "", 100, 0., 0.3);
+// ATH_CHECK(m_thistSvc->regHist("/file1/trap_pos", m_h_trap_pos));
+// }
+ ///////////////////////////////////////////////////
+
+ m_smallStepLength.setValue(m_smallStepLength.value() * CLHEP::micrometer);
+ m_tSurfaceDrift.setValue(m_tSurfaceDrift.value() * CLHEP::ns);
+
+ // Surface drift time calculation Stuff
+ m_tHalfwayDrift = m_tSurfaceDrift * 0.5;
+ m_distHalfInterStrip = m_distInterStrip * 0.5;
+ if ((m_tSurfaceDrift > m_tHalfwayDrift) and (m_tHalfwayDrift >= 0.0) and
+ (m_distHalfInterStrip > 0.0) and (m_distHalfInterStrip < m_distInterStrip)) {
+ m_SurfaceDriftFlag = true;
+ } else {
+ ATH_MSG_INFO("\tsurface drift still not on, wrong params");
+ }
+
+ // Make sure these two flags are not set simultaneously
+ if (m_tfix>-998. and m_tsubtract>-998.) {
+ ATH_MSG_FATAL("\tCannot set both FixedTime and SubtractTime options!");
+ ATH_MSG_INFO("\tMake sure the two flags are not set simultaneously in jo");
+ return StatusCode::FAILURE;
+ }
+
+ if (m_doDistortions) {
+ ATH_MSG_FATAL("Distortions not supported.");
+// ATH_CHECK(m_distortionsTool.retrieve());
+// } else {
+// m_distortionsTool.disable();
+ }
+
+ ATH_CHECK(m_lorentzAngleTool.retrieve());
+
+ return StatusCode::SUCCESS;
+}
+
+// ----------------------------------------------------------------------
+// finalize
+// ----------------------------------------------------------------------
+StatusCode SCT_SurfaceChargesGenerator::finalize() {
+ ATH_MSG_DEBUG("SCT_SurfaceChargesGenerator::finalize()");
+ return StatusCode::SUCCESS;
+}
+
+// ----------------------------------------------------------------------
+// perpandicular Drift time calculation
+// ----------------------------------------------------------------------
+float SCT_SurfaceChargesGenerator::driftTime(float zhit, const SiDetectorElement* element) const {
+ if (element==nullptr) {
+ ATH_MSG_ERROR("SCT_SurfaceChargesGenerator::process element is nullptr");
+ return -2.0;
+ }
+ const SCT_ModuleSideDesign* design{dynamic_cast<const SCT_ModuleSideDesign*>(&(element->design()))};
+ if (design==nullptr) {
+ ATH_MSG_ERROR("SCT_SurfaceChargesGenerator::process can not get " << design);
+ return -2.0;
+ }
+ const double thickness{design->thickness()};
+ const IdentifierHash hashId{element->identifyHash()};
+
+ if ((zhit < 0.0) or (zhit > thickness)) {
+ ATH_MSG_DEBUG("driftTime: hit coordinate zhit=" << zhit / CLHEP::micrometer << " out of range");
+ return -2.0;
+ }
+
+ float depletionVoltage{0.};
+ float biasVoltage{0.};
+ if (m_useSiCondDB) {
+ ATH_MSG_FATAL("Conditions DB voltages not supported.");
+ // depletionVoltage = m_siConditionsTool->depletionVoltage(hashId) * CLHEP::volt;
+ // biasVoltage = m_siConditionsTool->biasVoltage(hashId) * CLHEP::volt;
+ } else {
+ depletionVoltage = m_vdepl * CLHEP::volt;
+ biasVoltage = m_vbias * CLHEP::volt;
+ }
+
+ const float denominator{static_cast<float>(depletionVoltage + biasVoltage - (2.0 * zhit * depletionVoltage / thickness))};
+ if (denominator <= 0.0) {
+ if (biasVoltage >= depletionVoltage) { // Should not happen
+ if(not m_isOverlay) {
+ ATH_MSG_ERROR("driftTime: negative argument X for log(X) " << zhit);
+ }
+ return -1.0;
+ } else {
+ // (m_biasVoltage<m_depletionVoltage) can happen with underdepleted sensors, lose charges in that volume
+ return -10.0;
+ }
+ }
+
+ float t_drift{log((depletionVoltage + biasVoltage) / denominator)};
+ t_drift *= thickness * thickness / (2.0 * m_siPropertiesTool->getSiProperties(hashId).holeDriftMobility() * depletionVoltage);
+ return t_drift;
+}
+
+// ----------------------------------------------------------------------
+// Sigma diffusion calculation
+// ----------------------------------------------------------------------
+float SCT_SurfaceChargesGenerator::diffusionSigma(float zhit, const SiDetectorElement* element) const {
+ if (element==nullptr) {
+ ATH_MSG_ERROR("SCT_SurfaceChargesGenerator::diffusionSigma element is nullptr");
+ return 0.0;
+ }
+ const IdentifierHash hashId{element->identifyHash()};
+ const float t{driftTime(zhit, element)}; // in ns
+
+ if (t > 0.0) {
+ const float sigma{static_cast<float>(sqrt(2. * m_siPropertiesTool->getSiProperties(hashId).holeDiffusionConstant() * t))}; // in mm
+ return sigma;
+ } else {
+ return 0.0;
+ }
+}
+
+// ----------------------------------------------------------------------
+// Maximum drift time
+// ----------------------------------------------------------------------
+float SCT_SurfaceChargesGenerator::maxDriftTime(const SiDetectorElement* element) const {
+ if (element) {
+ const float sensorThickness{static_cast<float>(element->thickness())};
+ return driftTime(sensorThickness, element);
+ } else {
+ ATH_MSG_INFO("Error: SiDetectorElement not set!");
+ return 0.;
+ }
+}
+
+// ----------------------------------------------------------------------
+// Maximum Sigma difusion
+// ----------------------------------------------------------------------
+float SCT_SurfaceChargesGenerator::maxDiffusionSigma(const SiDetectorElement* element) const {
+ if (element) {
+ const float sensorThickness{static_cast<float>(element->thickness())};
+ return diffusionSigma(sensorThickness, element);
+ } else {
+ ATH_MSG_INFO("Error: SiDetectorElement not set!");
+ return 0.;
+ }
+}
+
+// ----------------------------------------------------------------------
+// Calculating the surface drift time but I should confess that
+// I haven't found out yet where the calculation come from
+// ----------------------------------------------------------------------
+float SCT_SurfaceChargesGenerator::surfaceDriftTime(float ysurf) const {
+ if (m_SurfaceDriftFlag) {
+ if ((ysurf >= 0.0) and (ysurf <= m_distInterStrip)) {
+ float t_surfaceDrift{0.};
+ if (ysurf < m_distHalfInterStrip) {
+ const float y{ysurf / m_distHalfInterStrip};
+ t_surfaceDrift= m_tHalfwayDrift * y * y;
+ } else {
+ const float y{(m_distInterStrip - ysurf) / (m_distInterStrip - m_distHalfInterStrip)};
+ t_surfaceDrift = m_tSurfaceDrift + (m_tHalfwayDrift - m_tSurfaceDrift) * y * y;
+ }
+ return t_surfaceDrift;
+ } else {
+ ATH_MSG_INFO(" ysurf out of range " << ysurf);
+ return -1.0;
+ }
+ } else {
+ return 0.0;
+ }
+}
+
+// -------------------------------------------------------------------------------------------
+// create a list of surface charges from a hit - called from SCT_Digitization
+// AthAlgorithm
+// -------------------------------------------------------------------------------------------
+void SCT_SurfaceChargesGenerator::process(const SiDetectorElement* element,
+ const TimedHitPtr<FaserSiHit>& phit,
+ const ISiSurfaceChargesInserter& inserter, CLHEP::HepRandomEngine * rndmEngine) const {
+ ATH_MSG_VERBOSE("SCT_SurfaceChargesGenerator::process starts");
+ processSiHit(element, *phit, inserter, phit.eventTime(), phit.eventId(), rndmEngine);
+ return;
+}
+
+// -------------------------------------------------------------------------------------------
+// create a list of surface charges from a hit - called from both AthAlgorithm
+// and PileUpTool
+// -------------------------------------------------------------------------------------------
+void SCT_SurfaceChargesGenerator::processSiHit(const SiDetectorElement* element,
+ const FaserSiHit& phit,
+ const ISiSurfaceChargesInserter& inserter,
+ float p_eventTime,
+ unsigned short p_eventId, CLHEP::HepRandomEngine * rndmEngine) const {
+ const SCT_ModuleSideDesign* design{dynamic_cast<const SCT_ModuleSideDesign*>(&(element->design()))};
+ if (design==nullptr) {
+ ATH_MSG_ERROR("SCT_SurfaceChargesGenerator::process can not get " << design);
+ return;
+ }
+
+ const double thickness{design->thickness()};
+ const IdentifierHash hashId{element->identifyHash()};
+ const double tanLorentz{m_lorentzAngleTool->getTanLorentzAngle(hashId)};
+
+ // ---**************************************
+ // Time of Flight Calculation - separate method?
+ // ---**************************************
+ // --- Original calculation of Time of Flight of the particle Time needed by the particle to reach the sensor
+ float timeOfFlight{p_eventTime + hitTime(phit)};
+
+ // Kondo 19/09/2007: Use the coordinate of the center of the module to calculate the time of flight
+ timeOfFlight -= (element->center().mag()) / CLHEP::c_light;
+ // !< extract the distance to the origin of the module to Time of flight
+
+ // !< timing set from jo to adjust (subtract) the timing
+ if (m_tsubtract > -998.) {
+ timeOfFlight -= m_tsubtract;
+ }
+ // ---**************************************
+
+ const CLHEP::Hep3Vector pos{phit.localStartPosition()};
+ const float xEta{static_cast<float>(pos[FaserSiHit::xEta])};
+ const float xPhi{static_cast<float>(pos[FaserSiHit::xPhi])};
+ const float xDep{static_cast<float>(pos[FaserSiHit::xDep])};
+
+ const CLHEP::Hep3Vector endPos{phit.localEndPosition()};
+ const float cEta{static_cast<float>(endPos[FaserSiHit::xEta]) - xEta};
+ const float cPhi{static_cast<float>(endPos[FaserSiHit::xPhi]) - xPhi};
+ const float cDep{static_cast<float>(endPos[FaserSiHit::xDep]) - xDep};
+
+ const float LargeStep{sqrt(cEta*cEta + cPhi*cPhi + cDep*cDep)};
+ const int numberOfSteps{static_cast<int>(LargeStep / m_smallStepLength) + 1};
+ const float steps{static_cast<float>(m_numberOfCharges * numberOfSteps)};
+ const float e1{static_cast<float>(phit.energyLoss() / steps)};
+ const float q1{static_cast<float>(e1 * m_siPropertiesTool->getSiProperties(hashId).electronHolePairsPerEnergy())};
+
+ // in the following, to test the code, we will use the original coordinate
+ // system of the SCTtest3SurfaceChargesGenerator x is eta y is phi z is depth
+ float xhit{xEta};
+ float yhit{xPhi};
+ float zhit{xDep};
+
+ if (m_doDistortions) {
+ ATH_MSG_FATAL("Distortions not supported.");
+ // if (element->isBarrel()) {// Only apply disortions to barrel modules
+ // Amg::Vector2D BOW;
+ // BOW[0] = m_distortionsTool->correctSimulation(hashId, xhit, yhit, cEta, cPhi, cDep)[0];
+ // BOW[1] = m_distortionsTool->correctSimulation(hashId, xhit, yhit, cEta, cPhi, cDep)[1];
+ // xhit = BOW.x();
+ // yhit = BOW.y();
+ // }
+ }
+
+ const float StepX{cEta / numberOfSteps};
+ const float StepY{cPhi / numberOfSteps};
+ const float StepZ{cDep / numberOfSteps};
+
+ // check the status of truth information for this SiHit
+ // some Truth information is cut for pile up events
+ const HepMcParticleLink trklink{HepMcParticleLink(phit.trackNumber(), p_eventId)};
+ SiCharge::Process hitproc{SiCharge::track};
+ if (phit.trackNumber() != 0) {
+ if (not trklink.isValid()) {
+ hitproc = SiCharge::cut_track;
+ }
+ }
+
+ float dstep{-0.5};
+ for (int istep{0}; istep < numberOfSteps; ++istep) {
+ dstep += 1.0;
+ float z1{zhit + StepZ * dstep};
+
+ // Distance between charge and readout side.
+ // design->readoutSide() is +1 if readout side is in +ve depth axis direction and visa-versa.
+ float zReadout{static_cast<float>(0.5 * thickness - design->readoutSide() * z1)};
+ // const double spess{zReadout};
+
+ // if (m_doHistoTrap) {
+ // m_h_depD->Fill(z1);
+ // m_h_spess->Fill(spess);
+ // }
+
+ float t_drift{driftTime(zReadout, element)}; // !< t_drift: perpandicular drift time
+ if (t_drift>-2.0000002 and t_drift<-1.9999998) {
+ ATH_MSG_DEBUG("Checking for rounding errors in compression");
+ if ((fabs(z1) - 0.5 * thickness) < 0.000010) {
+ ATH_MSG_DEBUG("Rounding error found attempting to correct it. z1 = " << std::fixed << std::setprecision(8) << z1);
+ if (z1 < 0.0) {
+ z1 = 0.0000005 - 0.5 * thickness;
+ // set new coordinate to be 0.5nm inside wafer volume.
+ } else {
+ z1 = 0.5 * thickness - 0.0000005;
+ // set new coordinate to be 0.5nm inside wafer volume.
+ }
+ zReadout = 0.5 * thickness - design->readoutSide() * z1;
+ t_drift = driftTime(zReadout, element);
+ if (t_drift>-2.0000002 and t_drift<-1.9999998) {
+ ATH_MSG_WARNING("Attempt failed. Making no correction.");
+ } else {
+ ATH_MSG_DEBUG("Correction Successful! z1 = " << std::fixed << std::setprecision(8) << z1 << ", zReadout = " << zReadout << ", t_drift = " << t_drift);
+ }
+ } else {
+ ATH_MSG_DEBUG("No rounding error found. Making no correction.");
+ }
+ }
+ if (t_drift > 0.0) {
+ const float x1{xhit + StepX * dstep};
+ float y1{yhit + StepY * dstep};
+ y1 += tanLorentz * zReadout; // !< Taking into account the magnetic field
+ const float sigma{diffusionSigma(zReadout, element)};
+
+ for (int i{0}; i < m_numberOfCharges; ++i) {
+ const float rx{CLHEP::RandGaussZiggurat::shoot(rndmEngine)};
+ const float xd{x1 + sigma * rx};
+ const float ry{CLHEP::RandGaussZiggurat::shoot(rndmEngine)};
+ const float yd{y1 + sigma * ry};
+
+ // For charge trapping with Ramo potential
+ const double stripPitch{0.080}; // mm
+ double dstrip{y1 / stripPitch}; // mm
+ if (dstrip > 0.) {
+ dstrip -= static_cast<double>(static_cast<int>(dstrip));
+ } else {
+ dstrip -= static_cast<double>(static_cast<int>(dstrip)) + 1;
+ }
+
+ // now y will be x and z will be y ....just to make sure to confuse everebody
+ // double y0{dstrip * stripPitch}; // mm
+ // double z0{thickness - zReadout}; // mm
+
+ // -- Charge Trapping
+ // if (m_doTrapping) {
+ // if (m_doHistoTrap) {
+ // m_h_zhit->Fill(zhit);
+ // }
+ // double trap_pos{-999999.}, drift_time{-999999.}; // FIXME need better default values
+ // if (chargeIsTrapped(spess, element, trap_pos, drift_time)) {
+ // if (not m_doRamo) {
+ // break;
+ // } else { // if we want to take into account also Ramo Potential
+ // double Q_m2{0.}, Q_m1{0.}, Q_00{0.}, Q_p1{0.}, Q_p2{0.}; // Charges
+
+ // dstrip = y1 / stripPitch; // mm
+ // // need the distance from the nearest strips
+ // // edge not centre, xtaka = 1/2*stripPitch
+ // // centre of detector, y1=0, is in the middle of
+ // // an interstrip gap
+ // if (dstrip > 0.) {
+ // dstrip -= static_cast<double>(static_cast<int>(dstrip));
+ // } else {
+ // dstrip -= static_cast<double>(static_cast<int>(dstrip)) + 1;
+ // }
+
+ // // now y will be x and z will be y ....just to make sure to confuse everebody
+ // double yfin{dstrip * stripPitch}; // mm
+ // double zfin{thickness - trap_pos}; // mm
+
+ // m_radDamageTool->holeTransport(y0, z0, yfin, zfin, Q_m2, Q_m1, Q_00, Q_p1, Q_p2);
+ // for (int strip{-2}; strip<=2; strip++) {
+ // const double ystrip{yd + strip * stripPitch}; // mm
+ // const SiLocalPosition position(element->hitLocalToLocal(xd, ystrip));
+ // if (design->inActiveArea(position)) {
+ // double charge{0.};
+ // if (strip == -2) charge = Q_m2;
+ // if (strip == -1) charge = Q_m1;
+ // if (strip == 0) charge = Q_00;
+ // if (strip == 1) charge = Q_p1;
+ // if (strip == 2) charge = Q_p2;
+ // const double time{drift_time};
+ // if (charge != 0.) {
+ // inserter(SiSurfaceCharge(position, SiCharge(q1*charge, time, hitproc, trklink)));
+ // continue;
+ // }
+ // }
+ // }
+ // ATH_MSG_INFO("strip zero charge = " << Q_00); // debug
+ // } // m_doRamo==true
+ // } // chargeIsTrapped()
+ // } // m_doTrapping==true
+
+ if (not m_doRamo) {
+ const SiLocalPosition position{element->hitLocalToLocal(xd, yd)};
+ if (design->inActiveArea(position)) {
+ const float sdist{static_cast<float>(design->scaledDistanceToNearestDiode(position))}; // !< dist on the surface from the hit point to the nearest strip (diode)
+ const float t_surf{surfaceDriftTime(2.0 * sdist)}; // !< Surface drift time
+ const float totaltime{(m_tfix > -998.) ? m_tfix.value() : t_drift + timeOfFlight + t_surf}; // !< Total drift time
+ inserter(SiSurfaceCharge(position, SiCharge(q1, totaltime, hitproc, trklink)));
+ } else {
+ ATH_MSG_VERBOSE(std::fixed << std::setprecision(8) << "Local position (phi, eta, depth): ("
+ << position.xPhi() << ", " << position.xEta() << ", " << position.xDepth()
+ << ") of the element is out of active area, charge = " << q1);
+ }
+ } // end of loop on charges
+ }
+ }
+ }
+ return;
+}
+
+// bool SCT_SurfaceChargesGenerator::chargeIsTrapped(double spess,
+// const SiDetectorElement* element,
+// double& trap_pos,
+// double& drift_time) const {
+// if (element==nullptr) {
+// ATH_MSG_ERROR("SCT_SurfaceChargesGenerator::chargeIsTrapped element is nullptr");
+// return false;
+// }
+// bool isTrapped{false};
+// const IdentifierHash hashId{element->identifyHash()};
+// const SCT_ChargeTrappingCondData condData{m_radDamageTool->getCondData(hashId, spess)};
+// const double electric_field{condData.getElectricField()};
+
+// if (m_doHistoTrap) {
+// const double mobChar{condData.getHoleDriftMobility()};
+// m_h_efieldz->Fill(spess, electric_field);
+// m_h_efield->Fill(electric_field);
+// m_h_mob_Char->Fill(electric_field, mobChar);
+// m_h_vel->Fill(electric_field, electric_field * mobChar);
+// }
+// const double t_electrode{condData.getTimeToElectrode()};
+// drift_time = condData.getTrappingTime();
+// const double z_trap{condData.getTrappingPositionZ()};
+// trap_pos = spess - z_trap;
+// if (m_doHistoTrap) {
+// m_h_drift_time->Fill(drift_time);
+// m_h_t_electrode->Fill(t_electrode);
+// m_h_drift_electrode->Fill(drift_time, t_electrode);
+// m_h_ztrap_tot->Fill(z_trap);
+// }
+// // -- Calculate if the charge is trapped, and at which distance
+// // -- Charge gets trapped before arriving to the electrode
+// if (drift_time < t_electrode) {
+// isTrapped = true;
+// ATH_MSG_INFO("drift_time: " << drift_time << " t_electrode: " << t_electrode << " spess " << spess);
+// ATH_MSG_INFO("z_trap: " << z_trap);
+// if (m_doHistoTrap) {
+// m_h_ztrap->Fill(z_trap);
+// m_h_trap_drift_t->Fill(drift_time);
+// m_h_drift1->Fill(spess, drift_time / t_electrode);
+// m_h_gen->Fill(spess, drift_time);
+// m_h_gen1->Fill(spess, z_trap);
+// m_h_gen2->Fill(spess, z_trap / drift_time * t_electrode);
+// m_h_velocity_trap->Fill(electric_field, z_trap / drift_time);
+// m_h_mobility_trap->Fill(electric_field, z_trap / drift_time / electric_field);
+// m_h_trap_pos->Fill(trap_pos);
+// }
+// } else {
+// isTrapped = false;
+// if (m_doHistoTrap) {
+// const double z_trap{condData.getTrappingPositionZ()};
+// m_h_no_ztrap->Fill(z_trap);
+// m_h_notrap_drift_t->Fill(drift_time);
+// }
+// }
+// return isTrapped;
+// }
diff --git a/Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_SurfaceChargesGenerator.h b/Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_SurfaceChargesGenerator.h
new file mode 100644
index 00000000..f9ae5fcd
--- /dev/null
+++ b/Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_SurfaceChargesGenerator.h
@@ -0,0 +1,164 @@
+// -*- C++ -*-
+
+/*
+ Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+*/
+
+/**
+ * Header file for SCT_SurfaceChargesGenerator Class
+ * @brief Class to drift the charged hits to the sensor surface for SCT
+ * version 1.0a Szymon Gadomski 31.05.2001
+ * @author Szymon.Gadomski@cern.ch, Awatif.Belymam@cern.ch, Davide.Costanzo@cern.ch,
+ * tgcornel@nikhef.nl, Grant.Gorfine@cern.ch, Paul.Bell@cern.ch, Jorgen.Dalmau@cern.ch,
+ * This is based on the SCTtest3SurfaceChargesGenerator
+ * compared to "test3", changes are all related to interfaces
+ * this should work with the new SCT_BarrelModuleSideDesign
+ * Compared to "test2" the following have been added in "test3":
+ * - a possibility to work with long steps and to smear charge along
+ * the step in smaller steps
+ * - Lorentz angle
+ * 19/04/2001 - change of hit class to SiTrackerHit
+ * version 1.0 Szymon Gadomski 09/06/2001, compiles with framework
+ * 26/07/2001 - changed to work with SiDigitization-00-02-00, tested
+ * 20/10/2002 - Awatif Belymam
+ * 29/10/2002 - the thickness of sensors is obtained from the objects
+ * - of SCT_ModuleSideDesign type instead from its jobOptions file.
+ * - Changes are done as well in jobOptions files.
+ * 06/01/2004 - Everything is now in CLHEP units (mm, ns, MeV)
+ * 23/08/2007 - Kondo.Gnanvo@cern.ch
+ * - Conversion of the SCT_SurfaceChargesGenerator code Al gTool
+ */
+
+#ifndef FASERSCT_DIGITIZATION_SCT_SURFACECHARGESGENERATOR_H
+#define FASERSCT_DIGITIZATION_SCT_SURFACECHARGESGENERATOR_H
+
+//Inheritance
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "FaserSCT_Digitization/ISCT_SurfaceChargesGenerator.h"
+
+#include "Identifier/IdentifierHash.h"
+// #include "InDetConditionsSummaryService/ISiliconConditionsTool.h"
+#include "InDetCondTools/ISiLorentzAngleTool.h"
+// #include "SCT_ConditionsTools/ISCT_RadDamageSummaryTool.h"
+// #include "SCT_ModuleDistortions/ISCT_ModuleDistortionsTool.h"
+#include "SiPropertiesTool/ISiPropertiesTool.h"
+
+// Gaudi
+#include "GaudiKernel/ITHistSvc.h" // for ITHistSvc
+#include "GaudiKernel/ServiceHandle.h"
+#include "GaudiKernel/ToolHandle.h"
+
+#include <iostream>
+
+// Charges and hits
+class FaserSiHit;
+
+// -- do charge trapping histos
+class ITHistSvc;
+class TH1F;
+class TH2F;
+class TProfile;
+
+namespace TrackerDD {
+ class SiDetectorElement;
+ class SCT_ModuleSideDesign;
+}
+
+namespace CLHEP {
+ class HepRandomEngine;
+}
+
+template <class HIT> class TimedHitPtr;
+
+class SCT_SurfaceChargesGenerator : public extends<AthAlgTool, ISCT_SurfaceChargesGenerator> {
+ public:
+
+ /** constructor */
+ SCT_SurfaceChargesGenerator(const std::string& type, const std::string& name, const IInterface* parent);
+
+ /** Destructor */
+ virtual ~SCT_SurfaceChargesGenerator() = default;
+
+ /** AlgTool initialize */
+ virtual StatusCode initialize();
+
+ /** AlgTool finalize */
+ virtual StatusCode finalize();
+
+ private:
+
+ virtual void setFixedTime(float fixedTime) {m_tfix = fixedTime;}
+
+ /** create a list of surface charges from a hit */
+ virtual void process(const TrackerDD::SiDetectorElement* element, const TimedHitPtr<FaserSiHit>& phit, const ISiSurfaceChargesInserter& inserter, CLHEP::HepRandomEngine * rndmEngine) const;
+ void processSiHit(const TrackerDD::SiDetectorElement* element, const FaserSiHit& phit, const ISiSurfaceChargesInserter& inserter, float eventTime, unsigned short eventID, CLHEP::HepRandomEngine * rndmEngine) const;
+
+ // some diagnostics methods are needed here too
+ float driftTime(float zhit, const TrackerDD::SiDetectorElement* element) const; //!< calculate drift time perpandicular to the surface for a charge at distance zhit from mid gap
+ float diffusionSigma(float zhit, const TrackerDD::SiDetectorElement* element) const; //!< calculate diffusion sigma from a gaussian dist scattered charge
+ float surfaceDriftTime(float ysurf) const; //!< Calculate of the surface drift time
+ float maxDriftTime(const TrackerDD::SiDetectorElement* element) const; //!< max drift charge equivalent to the detector thickness
+ float maxDiffusionSigma(const TrackerDD::SiDetectorElement* element) const; //!< max sigma diffusion
+
+ // trap_pos and drift_time are updated based on spess.
+// bool chargeIsTrapped(double spess, const TrackerDD::SiDetectorElement* element, double& trap_pos, double& drift_time) const;
+
+ IntegerProperty m_numberOfCharges{this, "NumberOfCharges", 1, "number of charges"};
+ FloatProperty m_smallStepLength{this, "SmallStepLength", 5, "max internal step along the larger G4 step"};
+
+ /** related to the surface drift */
+ FloatProperty m_tSurfaceDrift{this, "SurfaceDriftTime", 10, "max surface drift time"};
+
+ FloatProperty m_tfix{this, "FixedTime", -999., "fixed time"};
+ FloatProperty m_tsubtract{this, "SubtractTime", -999., "subtract drift time from mid gap"};
+
+ BooleanProperty m_doDistortions{this, "doDistortions", false, "Simulation of module distortions"};
+ BooleanProperty m_useSiCondDB{this, "UseSiCondDB", false, "Usage of SiConditions DB values can be disabled to use setable ones"};
+ FloatProperty m_vdepl{this, "DepletionVoltage", 70., "depletion voltage, default 70V"};
+ FloatProperty m_vbias{this, "BiasVoltage", 150., "bias voltage, default 150V"};
+ BooleanProperty m_doTrapping{this, "doTrapping", false, "Flag to set Charge Trapping"};
+// BooleanProperty m_doHistoTrap{this, "doHistoTrap", false, "Histogram the charge trapping effect"};
+ BooleanProperty m_doRamo{this, "doRamo", false, "Ramo Potential for charge trapping effect"};
+ BooleanProperty m_isOverlay{this, "isOverlay", false, "flag for overlay"};
+
+ //ToolHandles
+// ToolHandle<ISCT_ModuleDistortionsTool> m_distortionsTool{this, "SCTDistortionsTool", "SCT_DistortionsTool", "Tool to retrieve SCT distortions"};
+ ToolHandle<ISiPropertiesTool> m_siPropertiesTool{this, "SiPropertiesTool", "SCT_SiPropertiesTool", "Tool to retrieve SCT silicon properties"};
+// ToolHandle<ISCT_RadDamageSummaryTool> m_radDamageTool{this, "RadDamageSummaryTool", "SCT_RadDamageSummaryTool", "Tool to retrieve SCT radiation damages"};
+// ToolHandle<ISiliconConditionsTool> m_siConditionsTool{this, "SiConditionsTool", "SCT_SiliconConditionsTool", "Tool to retrieve SCT silicon information"};
+ ToolHandle<ISiLorentzAngleTool> m_lorentzAngleTool{this, "LorentzAngleTool", "SiLorentzAngleTool/SCTLorentzAngleTool", "Tool to retreive Lorentz angle"};
+
+ ServiceHandle<ITHistSvc> m_thistSvc{this, "THistSvc", "THistSvc"};
+
+ float m_tHalfwayDrift{0.}; //!< Surface drift time
+ float m_distInterStrip{1.0}; //!< Inter strip distance normalized to 1
+ float m_distHalfInterStrip{0.}; //!< Half way distance inter strip
+
+ bool m_SurfaceDriftFlag{false}; //!< surface drift ON/OFF
+
+ // -- Histograms
+ TProfile* m_h_efieldz{nullptr};
+ TH1F* m_h_efield{nullptr};
+ TH1F* m_h_spess{nullptr};
+ TH1F* m_h_depD{nullptr};
+ TH2F* m_h_drift_electrode{nullptr};
+ TH1F* m_h_ztrap{nullptr};
+ TH1F* m_h_drift_time{nullptr};
+ TH1F* m_h_t_electrode{nullptr};
+ TH1F* m_h_zhit{nullptr};
+ TH1F* m_h_ztrap_tot{nullptr};
+ TH1F* m_h_no_ztrap{nullptr};
+ TH1F* m_h_trap_drift_t{nullptr};
+ TH1F* m_h_notrap_drift_t{nullptr};
+ TProfile* m_h_mob_Char{nullptr};
+ TProfile* m_h_vel{nullptr};
+ TProfile* m_h_drift1{nullptr};
+ TProfile* m_h_gen{nullptr};
+ TProfile* m_h_gen1{nullptr};
+ TProfile* m_h_gen2{nullptr};
+ TProfile* m_h_velocity_trap{nullptr};
+ TProfile* m_h_mobility_trap{nullptr};
+ TH1F* m_h_trap_pos{nullptr};
+};
+
+#endif // SCT_SURFACECHARGESGENERATOR_H
diff --git a/Tracker/TrackerDigitization/FaserSCT_Digitization/src/components/FaserSCT_Digitization_entries.cxx b/Tracker/TrackerDigitization/FaserSCT_Digitization/src/components/FaserSCT_Digitization_entries.cxx
index dc705d93..ea833080 100644
--- a/Tracker/TrackerDigitization/FaserSCT_Digitization/src/components/FaserSCT_Digitization_entries.cxx
+++ b/Tracker/TrackerDigitization/FaserSCT_Digitization/src/components/FaserSCT_Digitization_entries.cxx
@@ -1,16 +1,18 @@
-// #include "../SCT_Amp.h"
-// #include "../SCT_FrontEnd.h"
-// #include "../SCT_Digitization.h"
+#include "../SCT_Amp.h"
+#include "../SCT_FrontEnd.h"
+#include "../SCT_Digitization.h"
#include "../SCT_DigitizationTool.h"
-// #include "../SCT_SurfaceChargesGenerator.h"
-// #include "../SCT_DetailedSurfaceChargesGenerator.h"
-// #include "../SCT_RandomDisabledCellGenerator.h"
+#include "../SCT_SurfaceChargesGenerator.h"
+#include "../SCT_RandomDisabledCellGenerator.h"
+// Unused by ATLAS
+// #include "../SCT_DetailedSurfaceChargesGenerator.h"
-// DECLARE_COMPONENT( SCT_Digitization )
-// DECLARE_COMPONENT( SCT_Amp )
-// DECLARE_COMPONENT( SCT_FrontEnd )
+DECLARE_COMPONENT( SCT_Amp )
+DECLARE_COMPONENT( SCT_FrontEnd )
+DECLARE_COMPONENT( SCT_Digitization )
DECLARE_COMPONENT( SCT_DigitizationTool )
-// DECLARE_COMPONENT( SCT_SurfaceChargesGenerator )
+DECLARE_COMPONENT( SCT_SurfaceChargesGenerator )
+DECLARE_COMPONENT( SCT_RandomDisabledCellGenerator )
+// Unused by ATLAS
// DECLARE_COMPONENT( SCT_DetailedSurfaceChargesGenerator )
-// DECLARE_COMPONENT( SCT_RandomDisabledCellGenerator )
diff --git a/Tracker/TrackerSimEvent/TrackerSimEvent/FaserSiHit.h b/Tracker/TrackerSimEvent/TrackerSimEvent/FaserSiHit.h
index 8669b880..b17acf1d 100644
--- a/Tracker/TrackerSimEvent/TrackerSimEvent/FaserSiHit.h
+++ b/Tracker/TrackerSimEvent/TrackerSimEvent/FaserSiHit.h
@@ -141,8 +141,8 @@ private:
HepMcParticleLink m_partLink;
unsigned int m_ID;
public:
- // enum
- // { xDep = 2, xPhi = 0, xEta = 1};
+ enum
+ { xDep = 2, xPhi = 0, xEta = 1};
};
--
GitLab