HGCalTestBeam/src/HGCalTB16SD01.cc

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 //#define EDM_ML_DEBUG
0011
0012 HGCalTB16SD01::HGCalTB16SD01(const std::string& name,
0013                              const SensitiveDetectorCatalog& clg,
0014                              edm::ParameterSet const& p,
0015                              const SimTrackManager* manager)
0016     : CaloSD(name, clg, p, manager), initialize_(true) {
0017   // Values from NIM 80 (1970) 239-244: as implemented in Geant3
0018   edm::ParameterSet m_HC = p.getParameter<edm::ParameterSet>("HGCalTestBeamSD");
0019   matName_ = m_HC.getParameter<std::string>("Material");
0020   useBirk_ = m_HC.getParameter<bool>("UseBirkLaw");
0021   birk1_ = m_HC.getParameter<double>("BirkC1") * (g / (MeV * cm2));
0022   birk2_ = m_HC.getParameter<double>("BirkC2");
0023   birk3_ = m_HC.getParameter<double>("BirkC3");
0024   matScin_ = nullptr;
0025
0026   edm::LogVerbatim("HGCSim") << "HGCalTB16SD01:: Use of Birks law is set to " << useBirk_ << " for " << matName_
0027                              << " with three constants kB = " << birk1_ << ", C1 = " << birk2_ << ", C2 = " << birk3_;
0028 }
0029
0030 double HGCalTB16SD01::getEnergyDeposit(const G4Step* aStep) {
0031   auto const point = aStep->GetPreStepPoint();
0032   if (initialize_)
0033     initialize(point);
0034   double destep = aStep->GetTotalEnergyDeposit();
0035   double wt2 = aStep->GetTrack()->GetWeight();
0036   double weight = (wt2 > 0.0) ? wt2 : 1.0;
0037   if (useBirk_ && matScin_ == point->GetMaterial()) {
0038     weight *= getAttenuation(aStep, birk1_, birk2_, birk3_);
0039   }
0040 #ifdef EDM_ML_DEBUG
0041   edm::LogVerbatim("HGCSim") << "HGCalTB16SD01: Detector " << point->GetTouchable()->GetVolume()->GetName() << " with "
0042                              << point->GetMaterial()->GetName() << " weight " << weight << ":" << wt2;
0043 #endif
0044   return weight * destep;
0045 }
0046
0047 uint32_t HGCalTB16SD01::setDetUnitId(const G4Step* aStep) {
0048   const G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0049   const G4VTouchable* touch = preStepPoint->GetTouchable();
0050
0051   int det(1), x(0), y(0);
0052   int lay = (touch->GetReplicaNumber(0));
0053
0054   return packIndex(det, lay, x, y);
0055 }
0056
0057 uint32_t HGCalTB16SD01::packIndex(int det, int lay, int x, int y) {
0058   int ix = 0, ixx = x;
0059   if (x < 0) {
0060     ix = 1;
0061     ixx = -x;
0062   }
0063   int iy = 0, iyy = y;
0064   if (y < 0) {
0065     iy = 1;
0066     iyy = -y;
0067   }
0068   uint32_t idx = (det & 15) << 28;  // bits 28-31
0069   idx += (lay & 127) << 21;         // bits 21-27
0070   idx += (iy & 1) << 19;            // bit  19
0071   idx += (iyy & 511) << 10;         // bits 10-18
0072   idx += (ix & 1) << 9;             // bit   9
0073   idx += (ixx & 511);               // bits  0-8
0074
0075 #ifdef EDM_ML_DEBUG
0076   edm::LogVerbatim("HGCSim") << "HGCalTB16SD01: Detector " << det << " Layer " << lay << " x " << x << " " << ix << " "
0077                              << ixx << " y " << y << " " << iy << " " << iyy << " ID " << std::hex << idx << std::dec;
0078 #endif
0079   return idx;
0080 }
0081
0082 void HGCalTB16SD01::unpackIndex(const uint32_t& idx, int& det, int& lay, int& x, int& y) {
0083   det = (idx >> 28) & 15;
0084   lay = (idx >> 21) & 127;
0085   y = (idx >> 10) & 511;
0086   if (((idx >> 19) & 1) == 1)
0087     y = -y;
0088   x = (idx)&511;
0089   if (((idx >> 9) & 1) == 1)
0090     x = -x;
0091 }
0092
0093 void HGCalTB16SD01::initialize(const G4StepPoint* point) {
0094   if (matName_ == point->GetMaterial()->GetName()) {
0095     matScin_ = point->GetMaterial();
0096     initialize_ = false;
0097   }
0098 #ifdef EDM_ML_DEBUG
0099   edm::LogVerbatim("HGCSim") << "HGCalTB16SD01: Material pointer for " << matName_
0100                              << " is initialized to : " << matScin_;
0101 #endif
0102 }