Manual sweep of MR !37385 to master

InnerDetector/InDetDigitization/SCT_Digitization/CMakeLists.txt

@@ -14,7 +14,7 @@ atlas_add_component( SCT_Digitization
                      src/*.cxx
                      src/components/*.cxx
                      INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${CLHEP_INCLUDE_DIRS}
-                     LINK_LIBRARIES ${ROOT_LIBRARIES} ${CLHEP_LIBRARIES} CommissionEvent AthenaBaseComps AthenaKernel CxxUtils PileUpToolsLib Identifier xAODEventInfo GaudiKernel SiDigitization InDetRawData InDetSimEvent HitManagement GeneratorObjects SiPropertiesToolLib InDetIdentifier InDetReadoutGeometry SCT_ReadoutGeometry InDetSimData InDetConditionsSummaryService PathResolver SCT_ConditionsToolsLib SCT_ModuleDistortionsLib )
+                     LINK_LIBRARIES ${ROOT_LIBRARIES} ${CLHEP_LIBRARIES} CommissionEvent AthenaBaseComps AthenaKernel CxxUtils PileUpToolsLib Identifier xAODEventInfo GaudiKernel SiDigitization InDetRawData InDetSimEvent HitManagement GeneratorObjects SiPropertiesToolLib InDetIdentifier InDetReadoutGeometry SCT_ReadoutGeometry InDetSimData InDetConditionsSummaryService PathResolver SCT_ConditionsToolsLib SCT_ModuleDistortionsLib MagFieldConditions )
 
 atlas_add_test( SCT_DigitizationMT_test
                 SCRIPT Digi_tf.py --inputHITSFile /cvmfs/atlas-nightlies.cern.ch/repo/data/data-art/DigitizationTests/HITS.04919495._001041.pool.root.1 --conditionsTag default:OFLCOND-RUN12-SDR-25 --digiSeedOffset1 170 --digiSeedOffset2 170 --geometryVersion ATLAS-R2-2015-03-01-00 --DataRunNumber 222525 --outputRDOFile mc15_2015_ttbar.RDO.pool.root --preInclude HITtoRDO:SimulationJobOptions/preInclude.SCTOnlyConfig.py,Digitization/ForceUseOfAlgorithms.py --postInclude Digitization/FixDataDependenciesForMT.py --skipEvents 0  --maxEvents 100 --athenaopts=--threads=10

InnerDetector/InDetDigitization/SCT_Digitization/src/SCT_InducedChargedModel.h → InnerDetector/InDetDigitization/SCT_Digitization/src/SCT_InducedChargeModel.h

@@ -15,15 +15,17 @@
 // Athena
 #include "AthenaKernel/IAtRndmGenSvc.h"
 #include "AthenaKernel/MsgStreamMember.h"
-#include "CxxUtils/checker_macros.h"
 #include "Identifier/IdentifierHash.h"
 #include "InDetConditionsSummaryService/ISiliconConditionsTool.h"
+#include "InDetReadoutGeometry/SiDetectorElement.h"
+#include "MagFieldConditions/AtlasFieldCacheCondObj.h"
 #include "PathResolver/PathResolver.h"
 #include "SiPropertiesTool/ISiPropertiesTool.h"
 
 // Gaudi
 #include "GaudiKernel/PhysicalConstants.h"
 #include "GaudiKernel/ToolHandle.h"
+#include "GaudiKernel/ServiceHandle.h"
 
 // Random number
 #include "CLHEP/Random/RandFlat.h"
@@ -32,51 +34,111 @@
 #include "CLHEP/Units/SystemOfUnits.h"
 
 // C++ Standard Library
+#include <array>
+#include <memory>
+#include <mutex>
 #include <string>
+#include <utility>
 #include <vector>
 
-class SCT_InducedChargedModel {
+class SCT_InducedChargeModel {
 
  public:
-  enum EFieldModel {UniformE=0, FlatDiodeModel=1, FEMsolutions=2};
+  enum EFieldModel {FlatDiodeModel=0, FEMsolutions=1, UniformE=2};
   enum TransportStep {NTransportSteps=100};
+  enum FEMNums {NRamoPoints=81, NEFieldPoints=17, NDepthPoints=115};
+  enum InducedStrips {StartStrip=-2, EndStrip=+2, Offset=-StartStrip, NStrips=EndStrip+Offset+1};
 
-  void init(const float vdepl, const float vbias, const IdentifierHash hashId,
-            ToolHandle<ISiliconConditionsTool>& siConditionsTool);
-  void setRandomEngine(CLHEP::HepRandomEngine* rndmEngine) { m_rndmEngine = rndmEngine; };
-  void setMagneticField(const double B) { m_B = - B*1000.; }; // m_B in [Tesla], B in kTesla
-
-  void transport(const bool isElectron,
+  struct SCT_InducedChargeModelData {
+    float m_VD; // full depletion voltage [Volt] negative for type-P
+    float m_VB; // applied bias voltage [Volt]
+    float m_T; // temperature
+    float m_depletion_depth;
+    const InDetDD::SiDetectorElement* m_element;
+    Amg::Vector3D m_magneticField;
+    EFieldModel m_EFieldModel;
+    CLHEP::HepRandomEngine* m_rndmEngine;
+    std::array<std::array<double, NDepthPoints>, NEFieldPoints> m_ExValue;
+    std::array<std::array<double, NDepthPoints>, NEFieldPoints> m_EyValue;
+
+    SCT_InducedChargeModelData(const float vdepl,
+                               const float vbias,
+                               const InDetDD::SiDetectorElement* element,
+                               const Amg::Vector3D& magneticField, // in kTesla
+                               const float bulk_depth,
+                               const EFieldModel model,
+                               const ToolHandle<ISiliconConditionsTool> siConditionsTool,
+                               CLHEP::HepRandomEngine* rndmEngine) {
+      m_VD = vdepl; // full depletion voltage [Volt] negative for type-P
+      m_VB = vbias; // applied bias voltage [Volt]
+      m_element = element;
+      m_magneticField = magneticField;
+
+      //------------ find delepletion deph for model=0 and 1 -------------
+      m_depletion_depth = bulk_depth;
+      // for type N (before type inversion)
+      if (m_VD >= 0.) {
+        if (m_VB < m_VD) m_depletion_depth = std::sqrt(m_VB/m_VD) * bulk_depth;
+      } else {
+        // for type P
+        if (m_VB <= std::abs(m_VD)) m_depletion_depth = 0.;
+      }
+
+      m_EFieldModel = model;
+      m_T = siConditionsTool->temperature(m_element->identifyHash()) + Gaudi::Units::STP_Temperature;
+      m_rndmEngine = rndmEngine;
+    }
+  };
+
+  SCT_InducedChargeModel(size_t maxHash, EFieldModel model=FEMsolutions);
+
+  SCT_InducedChargeModelData*
+    setWaferData(const float vdepl,
+                 const float vbias,
+                 const InDetDD::SiDetectorElement* element,
+                 const AtlasFieldCacheCondObj* fieldCondObj,
+                 const ToolHandle<ISiliconConditionsTool> siConditionsTool,
+                 CLHEP::HepRandomEngine* rndmEngine) const;
+
+  void setEField(SCT_InducedChargeModelData& data) const;
+
+  void transport(const SCT_InducedChargeModelData& data,
+                 const bool isElectron,
                  const double x0, const double y0,
                  double* Q_m2, double* Q_m1, double* Q_00, double* Q_p1, double* Q_p2,
                  const IdentifierHash hashId,
-                 ToolHandle<ISiPropertiesTool>& siPropertiesTool) const;
-  void holeTransport(const double x0, const double y0,
+                 const ToolHandle<ISiPropertiesTool> siPropertiesTool) const;
+  void holeTransport(const SCT_InducedChargeModelData& data,
+                     const double x0, const double y0,
                      double* Q_m2, double* Q_m1, double* Q_00, double* Q_p1, double* Q_p2,
                      const IdentifierHash hashId,
-                     ToolHandle<ISiPropertiesTool>& siPropertiesTool) const;
-  void electronTransport(const double x0, const double y0,
+                     const ToolHandle<ISiPropertiesTool> siPropertiesTool) const;
+  void electronTransport(const SCT_InducedChargeModelData& data,
+                         const double x0, const double y0,
                          double* Q_m2, double* Q_m1, double* Q_00, double* Q_p1, double* Q_p2,
                          const IdentifierHash hashId,
-                         ToolHandle<ISiPropertiesTool>& siPropertiesTool) const;
+                         const ToolHandle<ISiPropertiesTool> siPropertiesTool) const;
 
  private:
-  enum InducedStrips {StartStrip=-2, EndStrip=+2, Offset=-StartStrip, NStrips=EndStrip+Offset+1};
  
   void loadICMParameters();
 
-  bool getVxVyD(const bool isElectron,
+  bool getVxVyD(const SCT_InducedChargeModelData& data,
+                const bool isElectron,
                 const double x, const double y, double& vx, double& vy, double& D,
                 const IdentifierHash hashId,
-                ToolHandle<ISiPropertiesTool>& siPropertiesTool) const;
-  double induced(const int istrip, const double x, const double y) const;
-  void EField(const double x, const double y, double& Ex, double& Ey) const;
+                const ToolHandle<ISiPropertiesTool> siPropertiesTool) const;
+  double induced(const SCT_InducedChargeModelData& data,
+                 const int istrip, const double x, const double y) const;
+  void getEField(const SCT_InducedChargeModelData& data,
+                 const double x, const double y, double& Ex, double& Ey) const;
 
+  size_t getFEMIndex(SCT_InducedChargeModelData& data) const;
 
   /// Log a message using the Athena controlled logging system
   MsgStream& msg(MSG::Level lvl) const { return m_msg << lvl; }
   /// Check whether the logging system is active at the provided verbosity level
-  bool msgLvl(MSG::Level lvl) { return m_msg.get().level() <= lvl; }
+  bool msgLvl(MSG::Level lvl) const { return m_msg.get().level() <= lvl; }
 
   //-------- parameters for e, h transport --------------------------------
   static const double s_kB; // [m^2*kg/s^2/K]
@@ -85,21 +147,14 @@ class SCT_InducedChargedModel {
   // Private message stream member
   mutable Athena::MsgStreamMember m_msg ATLAS_THREAD_SAFE;
   
-  CLHEP::HepRandomEngine* m_rndmEngine; //!< Random number generation engine 
-
   //------parameters given externally by jobOptions ------------------
-  EFieldModel m_EFieldModel = FEMsolutions; // 0 (uniform E), 1 (flat diode model), 2 (FEM solusions)
-  double m_VD = -30.; // full depletion voltage [Volt] negative for type-P
-  double m_VB = 75.; // applied bias voltage [Volt]
-  double m_T;
-  double m_B;
+  EFieldModel m_EFieldModel; // 0 (flat diode model), 1 (FEM solusions), 2 (uniform E)
   double m_transportTimeStep = 0.50; // one step side in time [nsec]
   double m_transportTimeMax = m_transportTimeStep*NTransportSteps; // maximun tracing time [nsec]
 
   //------parameters mostly fixed but can be changed externally  ------------
   double m_bulk_depth =  0.0285; // in [cm]
   double m_strip_pitch = 0.0080; // in [cm]
-  double m_depletion_depth;
   double m_y_origin_min = 0.;
 
   //---------- arrays of FEM analysis -----------------------------------
@@ -110,11 +165,16 @@ class SCT_InducedChargedModel {
   // In sensor depth directions (for both potential and Electric field):
   //  114 divisions (115 points) with 2.5 nm intervals for 285 um.
 
-  enum FEMNums {NRamoPoints=81, NEFieldPoints=17, NDepthPoints=115};
+  std::vector<std::array<std::array<double, NDepthPoints>, NRamoPoints>> m_PotentialValue;
+  std::vector<std::array<std::array<double, NDepthPoints>, NEFieldPoints>> m_ExValue;
+  std::vector<std::array<std::array<double, NDepthPoints>, NEFieldPoints>> m_EyValue;
 
-  double m_PotentialValue[NRamoPoints][NDepthPoints];
-  double m_ExValue[NEFieldPoints][NDepthPoints];
-  double m_EyValue[NEFieldPoints][NDepthPoints];
+  // Cache of SCT_InducedChargeModelData for each wafer.
+  // Assuming wafer parameters do not change during a job.
+  // This should be almost always satisfied.
+  mutable std::vector<std::unique_ptr<SCT_InducedChargeModelData>> m_data ATLAS_THREAD_SAFE; // Guarded by m_mutex
+  // To protect concurrent access to m_data
+  mutable std::mutex m_mutex;
 
   static const std::vector<float> s_VFD0;
   static const std::vector<std::string> s_VFD0str;

InnerDetector/InDetDigitization/SCT_Digitization/src/SCT_SurfaceChargesGenerator.cxx

@@ -10,6 +10,7 @@
 
 // Athena
 #include "GeneratorObjects/HepMcParticleLink.h"
+#include "InDetIdentifier/SCT_ID.h"
 #include "InDetSimEvent/SiHit.h" // for SiHit, SiHit::::xDep, etc
 #include "HitManagement/TimedHitPtr.h" // for TimedHitPtr
 
@@ -49,6 +50,8 @@ StatusCode SCT_SurfaceChargesGenerator::initialize() {
   // Get ISiliconConditionsSvc
   ATH_CHECK(m_siConditionsTool.retrieve());
 
+  ATH_CHECK(m_fieldCacheCondObjInputKey.initialize(m_doInducedChargeModel));
+
   if (m_doTrapping) {
     // -- Get Radiation Damage Tool
     ATH_CHECK(m_radDamageTool.retrieve());
@@ -129,6 +132,13 @@ StatusCode SCT_SurfaceChargesGenerator::initialize() {
   }
   ///////////////////////////////////////////////////
 
+  // Induced Charge Module.
+  if (m_doInducedChargeModel) {
+    const SCT_ID* sct_id{nullptr};
+    ATH_CHECK(detStore()->retrieve(sct_id, "SCT_ID"));
+    m_InducedChargeModel = std::make_unique<SCT_InducedChargeModel>(sct_id->wafer_hash_max());
+  }
+
   // Surface drift time calculation Stuff
   m_tHalfwayDrift = m_tSurfaceDrift * 0.5;
   m_distHalfInterStrip = m_distInterStrip * 0.5;
@@ -304,7 +314,7 @@ void SCT_SurfaceChargesGenerator::processSiHit(const SiDetectorElement* element,
                                                const SiHit& phit,
                                                const ISiSurfaceChargesInserter& inserter,
                                                float p_eventTime,
-                                               unsigned short p_eventId, CLHEP::HepRandomEngine * rndmEngine) const {
+                                               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);
@@ -353,6 +363,24 @@ void SCT_SurfaceChargesGenerator::processSiHit(const SiDetectorElement* element,
   float yhit{xPhi};
   float zhit{xDep};
 
+  SCT_InducedChargeModel::SCT_InducedChargeModelData* data{nullptr};
+  if (m_doInducedChargeModel) { // Setting magnetic field for the ICM.
+    SG::ReadCondHandle<AtlasFieldCacheCondObj> readHandle{m_fieldCacheCondObjInputKey, Gaudi::Hive::currentContext()};
+    const AtlasFieldCacheCondObj* fieldCondObj{*readHandle};
+    float vdepl{m_vdepl};
+    float vbias{m_vbias};
+    if (m_useSiCondDB) {
+      vdepl = m_siConditionsTool->depletionVoltage(hashId);
+      vbias = m_siConditionsTool->biasVoltage(hashId);
+    }
+    data = m_InducedChargeModel->setWaferData(vdepl,
+                                              vbias,
+                                              element,
+                                              fieldCondObj,
+                                              m_siConditionsTool,
+                                              rndmEngine);
+  }
+
   if (m_doDistortions) {
     if (element->isBarrel()) {// Only apply disortions to barrel modules
       Amg::Vector2D BOW;
@@ -420,8 +448,12 @@ void SCT_SurfaceChargesGenerator::processSiHit(const SiDetectorElement* element,
     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)};
+
+      float sigma{0.};
+      if (not m_doInducedChargeModel) {
+        sigma = diffusionSigma(zReadout, element);
+        y1 += tanLorentz * zReadout; // !< Taking into account the magnetic field
+      } // These are treated in Induced Charge Model.
 
       for (int i{0}; i < m_numberOfCharges; ++i) {
         const float rx{CLHEP::RandGaussZiggurat::shoot(rndmEngine)};
@@ -491,21 +523,59 @@ void SCT_SurfaceChargesGenerator::processSiHit(const SiDetectorElement* element,
             } // 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);
+          if (m_doInducedChargeModel) { // Induced Charge Model
+            // Charges storages for 50 ns. 0.5 ns steps.
+            double Q_m2[SCT_InducedChargeModel::NTransportSteps]={0};
+            double Q_m1[SCT_InducedChargeModel::NTransportSteps]={0};
+            double Q_00[SCT_InducedChargeModel::NTransportSteps]={0};
+            double Q_p1[SCT_InducedChargeModel::NTransportSteps]={0};
+            double Q_p2[SCT_InducedChargeModel::NTransportSteps]={0};
+
+            const double mm2cm = 0.1; // For mm -> cm conversion
+            // Unit for y and z : mm -> cm in SCT_InducedChargeModel
+            m_InducedChargeModel->holeTransport(*data,
+                                                y0*mm2cm, z0*mm2cm,
+                                                Q_m2, Q_m1, Q_00, Q_p1, Q_p2,
+                                                hashId, m_siPropertiesTool);
+            m_InducedChargeModel->electronTransport(*data,
+                                                    y0*mm2cm, z0*mm2cm,
+                                                    Q_m2, Q_m1, Q_00, Q_p1, Q_p2,
+                                                    hashId, m_siPropertiesTool);
+
+            for (int it{0}; it<SCT_InducedChargeModel::NTransportSteps; it++) {
+              if (Q_00[it] == 0.0) continue;
+              double ICM_time{(it+0.5)*0.5 + timeOfFlight};
+              double Q_new[SCT_InducedChargeModel::NStrips]{
+                Q_m2[it], Q_m1[it], Q_00[it], Q_p1[it], Q_p2[it]
+              };
+              for (int strip{SCT_InducedChargeModel::StartStrip}; strip<=SCT_InducedChargeModel::EndStrip; strip++) {
+                double ystrip{y1 + strip * stripPitch};
+                SiLocalPosition position{element->hitLocalToLocal(x1, ystrip)};
+                if (design->inActiveArea(position)) {
+                  inserter(SiSurfaceCharge(position,
+                                           SiCharge(q1 * Q_new[strip+SCT_InducedChargeModel::Offset],
+                                                    ICM_time, hitproc, trklink)));
+                }
+              }
+            }
+          } else { // not m_doInducedChargeModel
+            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
-      }
+        }
+      } // end of loop on charges
     }
   }
   return;

InnerDetector/InDetDigitization/SCT_Digitization/src/SCT_SurfaceChargesGenerator.h

 // -*- C++ -*-
 
 /*
-  Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
 */
 
 /**
@@ -36,12 +36,16 @@
 #include "AthenaBaseComps/AthAlgTool.h"
 #include "SCT_Digitization/ISCT_SurfaceChargesGenerator.h"
 
+#include "SCT_InducedChargeModel.h"
+
 #include "Identifier/IdentifierHash.h"
 #include "InDetConditionsSummaryService/ISiliconConditionsTool.h"
 #include "InDetCondTools/ISiLorentzAngleTool.h"
+#include "MagFieldConditions/AtlasFieldCacheCondObj.h"
 #include "SCT_ConditionsTools/ISCT_RadDamageSummaryTool.h"
 #include "SCT_ModuleDistortions/ISCT_ModuleDistortionsTool.h"
 #include "SiPropertiesTool/ISiPropertiesTool.h"
+#include "StoreGate/ReadCondHandle.h"
 
 // Gaudi
 #include "GaudiKernel/ITHistSvc.h" // for ITHistSvc
@@ -123,6 +127,7 @@ class SCT_SurfaceChargesGenerator : public extends<AthAlgTool, ISCT_SurfaceCharg
   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"};
+  BooleanProperty m_doInducedChargeModel{this, "doInducedChargeModel", false, "Flag for Induced Charge Model"};
 
   //ToolHandles
   ToolHandle<ISCT_ModuleDistortionsTool> m_distortionsTool{this, "SCTDistortionsTool", "SCT_DistortionsTool", "Tool to retrieve SCT distortions"};
@@ -133,6 +138,9 @@ class SCT_SurfaceChargesGenerator : public extends<AthAlgTool, ISCT_SurfaceCharg
 
   ServiceHandle<ITHistSvc> m_thistSvc{this, "THistSvc", "THistSvc"};
 
+  SG::ReadCondHandleKey<AtlasFieldCacheCondObj> m_fieldCacheCondObjInputKey{
+    this, "AtlasFieldCacheCondObj", "fieldCondObj", "Name of the Magnetic Field conditions object key"};
+
   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
@@ -162,6 +170,9 @@ class SCT_SurfaceChargesGenerator : public extends<AthAlgTool, ISCT_SurfaceCharg
   TProfile* m_h_velocity_trap{nullptr};
   TProfile* m_h_mobility_trap{nullptr};
   TH1F* m_h_trap_pos{nullptr};
+
+  // For Induced Charge Module, M.Togawa
+  std::unique_ptr<SCT_InducedChargeModel> m_InducedChargeModel;
 };
 
 #endif // SCT_SURFACECHARGESGENERATOR_H