File indexing completed on 2023-03-17 11:24:36
0001 #include "FWCore/Utilities/interface/Exception.h"
0002 #include "SimG4CMS/HGCalTestBeam/interface/HGCalTB16SD01.h"
0003
0004 #include "CLHEP/Units/GlobalSystemOfUnits.h"
0005 #include "G4LogicalVolumeStore.hh"
0006 #include "G4Track.hh"
0007
0008 #include <string>
0009
0010
0011
0012
0013 HGCalTB16SD01::HGCalTB16SD01(const std::string& name,
0014 const SensitiveDetectorCatalog& clg,
0015 edm::ParameterSet const& p,
0016 const SimTrackManager* manager)
0017 : CaloSD(name, clg, p, manager), initialize_(true) {
0018
0019 edm::ParameterSet m_HC = p.getParameter<edm::ParameterSet>("HGCalTestBeamSD");
0020 matName_ = m_HC.getParameter<std::string>("Material");
0021 useBirk_ = m_HC.getParameter<bool>("UseBirkLaw");
0022 birk1_ = m_HC.getParameter<double>("BirkC1") * (g / (MeV * cm2));
0023 birk2_ = m_HC.getParameter<double>("BirkC2");
0024 birk3_ = m_HC.getParameter<double>("BirkC3");
0025 matScin_ = nullptr;
0026
0027 edm::LogVerbatim("HGCSim") << "HGCalTB16SD01:: Use of Birks law is set to " << useBirk_ << " for " << matName_
0028 << " with three constants kB = " << birk1_ << ", C1 = " << birk2_ << ", C2 = " << birk3_;
0029 }
0030
0031 double HGCalTB16SD01::getEnergyDeposit(const G4Step* aStep) {
0032 auto const point = aStep->GetPreStepPoint();
0033 if (initialize_)
0034 initialize(point);
0035 double destep = aStep->GetTotalEnergyDeposit();
0036 double wt2 = aStep->GetTrack()->GetWeight();
0037 double weight = (wt2 > 0.0) ? wt2 : 1.0;
0038 if (useBirk_ && matScin_ == point->GetMaterial()) {
0039 weight *= getAttenuation(aStep, birk1_, birk2_, birk3_);
0040 }
0041 #ifdef EDM_ML_DEBUG
0042 edm::LogVerbatim("HGCSim") << "HGCalTB16SD01: Detector " << point->GetTouchable()->GetVolume()->GetName() << " with "
0043 << point->GetMaterial()->GetName() << " weight " << weight << ":" << wt2;
0044 #endif
0045 return weight * destep;
0046 }
0047
0048 uint32_t HGCalTB16SD01::setDetUnitId(const G4Step* aStep) {
0049 const G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0050 const G4VTouchable* touch = preStepPoint->GetTouchable();
0051
0052 int det(1), x(0), y(0);
0053 int lay = (touch->GetReplicaNumber(0));
0054
0055 return packIndex(det, lay, x, y);
0056 }
0057
0058 uint32_t HGCalTB16SD01::packIndex(int det, int lay, int x, int y) {
0059 int ix = 0, ixx = x;
0060 if (x < 0) {
0061 ix = 1;
0062 ixx = -x;
0063 }
0064 int iy = 0, iyy = y;
0065 if (y < 0) {
0066 iy = 1;
0067 iyy = -y;
0068 }
0069 uint32_t idx = (det & 15) << 28;
0070 idx += (lay & 127) << 21;
0071 idx += (iy & 1) << 19;
0072 idx += (iyy & 511) << 10;
0073 idx += (ix & 1) << 9;
0074 idx += (ixx & 511);
0075
0076 #ifdef EDM_ML_DEBUG
0077 edm::LogVerbatim("HGCSim") << "HGCalTB16SD01: Detector " << det << " Layer " << lay << " x " << x << " " << ix << " "
0078 << ixx << " y " << y << " " << iy << " " << iyy << " ID " << std::hex << idx << std::dec;
0079 #endif
0080 return idx;
0081 }
0082
0083 void HGCalTB16SD01::unpackIndex(const uint32_t& idx, int& det, int& lay, int& x, int& y) {
0084 det = (idx >> 28) & 15;
0085 lay = (idx >> 21) & 127;
0086 y = (idx >> 10) & 511;
0087 if (((idx >> 19) & 1) == 1)
0088 y = -y;
0089 x = (idx)&511;
0090 if (((idx >> 9) & 1) == 1)
0091 x = -x;
0092 }
0093
0094 void HGCalTB16SD01::initialize(const G4StepPoint* point) {
0095 if (matName_ == point->GetMaterial()->GetName()) {
0096 matScin_ = point->GetMaterial();
0097 initialize_ = false;
0098 }
0099 #ifdef EDM_ML_DEBUG
0100 edm::LogVerbatim("HGCSim") << "HGCalTB16SD01: Material pointer for " << matName_
0101 << " is initialized to : " << matScin_;
0102 #endif
0103 }