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
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   // Values from NIM 80 (1970) 239-244: as implemented in Geant3
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;  // bits 28-31
0070   idx += (lay & 127) << 21;         // bits 21-27
0071   idx += (iy & 1) << 19;            // bit  19
0072   idx += (iyy & 511) << 10;         // bits 10-18
0073   idx += (ix & 1) << 9;             // bit   9
0074   idx += (ixx & 511);               // bits  0-8
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 }