First version of SCT_Digitization

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

Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_Amp.cxx

+/*
+  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;
+}

Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_Amp.h

+// -*- 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

Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_Digitization.cxx

+/*
+  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();
+}

Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_Digitization.h

+// -*- 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

Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_FrontEnd.cxx

+/*
+  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;
+  }
+}

Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_FrontEnd.h

+// -*- 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

Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_RandomDisabledCellGenerator.cxx

+/*
+  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);
+    }
+  }
+}

Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_RandomDisabledCellGenerator.h

+// -*- 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

Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_SurfaceChargesGenerator.cxx

+/*
+  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;
+// }

Tracker/TrackerDigitization/FaserSCT_Digitization/src/SCT_SurfaceChargesGenerator.h

+// -*- 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

Tracker/TrackerDigitization/FaserSCT_Digitization/src/components/FaserSCT_Digitization_entries.cxx

-// #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_SurfaceChargesGenerator.h"
+#include "../SCT_RandomDisabledCellGenerator.h"
+// Unused by ATLAS
 // #include "../SCT_DetailedSurfaceChargesGenerator.h" 
-// #include "../SCT_RandomDisabledCellGenerator.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 )
 

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};
 };