File indexing completed on 2023-10-25 10:04:15
0001 #include "SimG4CMS/HGCalTestBeam/interface/AHCalSD.h"
0002 #include "SimG4CMS/HGCalTestBeam/interface/AHCalDetId.h"
0003 #include "SimG4Core/Notification/interface/TrackInformation.h"
0004 #include "Geometry/HGCalCommonData/interface/AHCalParameters.h"
0005
0006 #include "DataFormats/DetId/interface/DetId.h"
0007 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0008
0009 #include "G4LogicalVolume.hh"
0010 #include "G4LogicalVolumeStore.hh"
0011 #include "G4ParticleTable.hh"
0012 #include "G4Track.hh"
0013 #include "G4VProcess.hh"
0014
0015 #include "G4PhysicalConstants.hh"
0016 #include "G4SystemOfUnits.hh"
0017
0018 #include <iomanip>
0019 #include <map>
0020
0021
0022
0023 AHCalSD::AHCalSD(const std::string& name,
0024 const SensitiveDetectorCatalog& clg,
0025 edm::ParameterSet const& p,
0026 const SimTrackManager* manager)
0027 : CaloSD(name,
0028 clg,
0029 p,
0030 manager,
0031 (float)(p.getParameter<edm::ParameterSet>("AHCalSD").getParameter<double>("TimeSliceUnit")),
0032 p.getParameter<edm::ParameterSet>("AHCalSD").getParameter<bool>("IgnoreTrackID")) {
0033 edm::ParameterSet m_HC = p.getParameter<edm::ParameterSet>("AHCalSD");
0034 useBirk = m_HC.getParameter<bool>("UseBirkLaw");
0035 birk1 = m_HC.getParameter<double>("BirkC1") * (CLHEP::g / (CLHEP::MeV * CLHEP::cm2));
0036 birk2 = m_HC.getParameter<double>("BirkC2");
0037 birk3 = m_HC.getParameter<double>("BirkC3");
0038 eminHit = m_HC.getParameter<double>("EminHit") * CLHEP::MeV;
0039
0040 edm::LogVerbatim("HcalSim") << "AHCalSD:: Use of Birks law is set to " << useBirk
0041 << " with three constants kB = " << birk1 << ", C1 = " << birk2 << ", C2 = " << birk3
0042 << "\nAHCalSD:: Threshold for storing"
0043 << " hits: " << eminHit;
0044 }
0045
0046 double AHCalSD::getEnergyDeposit(const G4Step* aStep) {
0047 double destep = aStep->GetTotalEnergyDeposit();
0048 double wt2 = aStep->GetTrack()->GetWeight();
0049 double weight = (wt2 > 0.0) ? wt2 : 1.0;
0050 #ifdef EDM_ML_DEBUG
0051 double weight0 = weight;
0052 #endif
0053 if (useBirk)
0054 weight *= getAttenuation(aStep, birk1, birk2, birk3);
0055 #ifdef EDM_ML_DEBUG
0056 edm::LogVerbatim("HcalSim") << "AHCalSD: weight " << weight0 << " " << weight;
0057 #endif
0058 double edep = weight * destep;
0059 return edep;
0060 }
0061
0062 uint32_t AHCalSD::setDetUnitId(const G4Step* aStep) {
0063 const G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0064 const G4VTouchable* touch = preStepPoint->GetTouchable();
0065
0066 int depth = (touch->GetReplicaNumber(1));
0067 int incol = ((touch->GetReplicaNumber(0)) % AHCalParameters::kColumn_);
0068 int inrow = ((touch->GetReplicaNumber(0)) / AHCalParameters::kColumn_) % AHCalParameters::kRow_;
0069 int jncol = ((touch->GetReplicaNumber(0)) / AHCalParameters::kRowColumn_) % AHCalParameters::kSign_;
0070 int jnrow = ((touch->GetReplicaNumber(0)) / AHCalParameters::kSignRowColumn_) % AHCalParameters::kSign_;
0071 int col = (jncol == 0) ? incol : -incol;
0072 int row = (jnrow == 0) ? inrow : -inrow;
0073 uint32_t index = AHCalDetId(row, col, depth).rawId();
0074 #ifdef EDM_ML_DEBUG
0075 edm::LogVerbatim("HcalSim") << "AHCalSD: det = " << HcalOther << " depth = " << depth << " row = " << row
0076 << " column = " << col << " packed index = 0x" << std::hex << index << std::dec;
0077 bool flag = unpackIndex(index, row, col, depth);
0078 edm::LogVerbatim("HcalSim") << "Results from unpacker for 0x" << std::hex << index << std::dec << " Flag " << flag
0079 << " Row " << row << " Col " << col << " Depth " << depth;
0080 #endif
0081 return index;
0082 }
0083
0084 bool AHCalSD::unpackIndex(const uint32_t& idx, int& row, int& col, int& depth) {
0085 DetId gen(idx);
0086 HcalSubdetector subdet = (HcalSubdetector(gen.subdetId()));
0087 bool rcode = (gen.det() == DetId::Hcal && subdet != HcalOther);
0088 row = col = depth = 0;
0089 if (rcode) {
0090 row = AHCalDetId(idx).irow();
0091 col = AHCalDetId(idx).icol();
0092 depth = AHCalDetId(idx).depth();
0093 }
0094 #ifdef EDM_ML_DEBUG
0095 edm::LogVerbatim("HcalSim") << "AHCalSD: packed index = 0x" << std::hex << idx << std::dec << " Row " << row
0096 << " Column " << col << " Depth " << depth << " OK " << rcode;
0097 #endif
0098 return rcode;
0099 }
0100
0101 bool AHCalSD::filterHit(CaloG4Hit* aHit, double time) {
0102 return ((time <= tmaxHit) && (aHit->getEnergyDeposit() > eminHit));
0103 }
