File indexing completed on 2023-01-23 02:42:52
0001
0002
0003
0004
0005
0006 #include "DataFormats/Math/interface/FastMath.h"
0007 #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
0008 #include "SimG4CMS/Calo/interface/HGCalSD.h"
0009 #include "SimG4Core/Notification/interface/TrackInformation.h"
0010 #include "FWCore/Utilities/interface/Exception.h"
0011 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0012 #include "Geometry/HGCalCommonData/interface/HGCalGeometryMode.h"
0013 #include "G4LogicalVolumeStore.hh"
0014 #include "G4LogicalVolume.hh"
0015 #include "G4Step.hh"
0016 #include "G4Track.hh"
0017 #include "G4ParticleTable.hh"
0018 #include "G4VProcess.hh"
0019 #include "G4Trap.hh"
0020
0021 #include <fstream>
0022 #include <iomanip>
0023 #include <iostream>
0024 #include <memory>
0025
0026
0027
0028 HGCalSD::HGCalSD(const std::string& name,
0029 const HGCalDDDConstants* hgc,
0030 const SensitiveDetectorCatalog& clg,
0031 edm::ParameterSet const& p,
0032 const SimTrackManager* manager)
0033 : CaloSD(name,
0034 clg,
0035 p,
0036 manager,
0037 static_cast<float>(p.getParameter<edm::ParameterSet>("HGCSD").getParameter<double>("TimeSliceUnit")),
0038 p.getParameter<edm::ParameterSet>("HGCSD").getParameter<bool>("IgnoreTrackID")),
0039 hgcons_(hgc),
0040 slopeMin_(0),
0041 levelT1_(99),
0042 levelT2_(99),
0043 useSimWt_(0),
0044 tan30deg_(std::tan(30.0 * CLHEP::deg)) {
0045 numberingScheme_.reset(nullptr);
0046 mouseBite_.reset(nullptr);
0047
0048 edm::ParameterSet m_HGC = p.getParameter<edm::ParameterSet>("HGCSD");
0049 eminHit_ = m_HGC.getParameter<double>("EminHit") * CLHEP::MeV;
0050 fiducialCut_ = m_HGC.getParameter<bool>("FiducialCut");
0051 distanceFromEdge_ = m_HGC.getParameter<double>("DistanceFromEdge");
0052 storeAllG4Hits_ = m_HGC.getParameter<bool>("StoreAllG4Hits");
0053 rejectMB_ = m_HGC.getParameter<bool>("RejectMouseBite");
0054 waferRot_ = m_HGC.getParameter<bool>("RotatedWafer");
0055 cornerMinMask_ = m_HGC.getParameter<int>("CornerMinMask");
0056 angles_ = m_HGC.getUntrackedParameter<std::vector<double>>("WaferAngles");
0057 missingFile_ = m_HGC.getUntrackedParameter<std::string>("MissingWaferFile");
0058 checkID_ = m_HGC.getUntrackedParameter<bool>("CheckID");
0059
0060 if (storeAllG4Hits_) {
0061 setUseMap(false);
0062 setNumberCheckedHits(0);
0063 }
0064
0065
0066 G4String myName = name;
0067 mydet_ = DetId::Forward;
0068 nameX_ = "HGCal";
0069 if (myName.find("HitsEE") != std::string::npos) {
0070 mydet_ = DetId::HGCalEE;
0071 nameX_ = "HGCalEESensitive";
0072 } else if (myName.find("HitsHEfront") != std::string::npos) {
0073 mydet_ = DetId::HGCalHSi;
0074 nameX_ = "HGCalHESiliconSensitive";
0075 }
0076
0077 #ifdef EDM_ML_DEBUG
0078 edm::LogVerbatim("HGCSim") << "**************************************************"
0079 << "\n"
0080 << "* *"
0081 << "\n"
0082 << "* Constructing a HGCalSD with name " << name << "\n"
0083 << "* *"
0084 << "\n"
0085 << "**************************************************";
0086 #endif
0087 edm::LogVerbatim("HGCSim") << "HGCalSD:: Threshold for storing hits: " << eminHit_ << " for " << nameX_
0088 << " detector " << mydet_;
0089 edm::LogVerbatim("HGCSim") << "Flag for storing individual Geant4 Hits " << storeAllG4Hits_;
0090 edm::LogVerbatim("HGCSim") << "Fiducial volume cut with cut from eta/phi "
0091 << "boundary " << fiducialCut_ << " at " << distanceFromEdge_;
0092 edm::LogVerbatim("HGCSim") << "Reject MosueBite Flag: " << rejectMB_ << " cuts along " << angles_.size()
0093 << " axes: " << angles_[0] << ", " << angles_[1];
0094 }
0095
0096 double HGCalSD::getEnergyDeposit(const G4Step* aStep) {
0097 double r = aStep->GetPreStepPoint()->GetPosition().perp();
0098 double z = std::abs(aStep->GetPreStepPoint()->GetPosition().z());
0099 #ifdef EDM_ML_DEBUG
0100 G4int parCode = aStep->GetTrack()->GetDefinition()->GetPDGEncoding();
0101 G4String parName = aStep->GetTrack()->GetDefinition()->GetParticleName();
0102 G4LogicalVolume* lv = aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
0103 edm::LogVerbatim("HGCSim") << "HGCalSD: Hit from standard path from " << lv->GetName() << " for Track "
0104 << aStep->GetTrack()->GetTrackID() << " (" << parCode << ":" << parName << ") R = " << r
0105 << " Z = " << z << " slope = " << r / z << ":" << slopeMin_;
0106 #endif
0107
0108 if (r < z * slopeMin_) {
0109 #ifdef EDM_ML_DEBUG
0110 edm::LogVerbatim("HGCSim") << "HGCalSD: Fiducial Volume cut";
0111 #endif
0112 return 0.0;
0113 }
0114
0115 double wt1 = getResponseWt(aStep->GetTrack());
0116 double wt2 = aStep->GetTrack()->GetWeight();
0117 double destep = weight_ * wt1 * (aStep->GetTotalEnergyDeposit());
0118 if (wt2 > 0)
0119 destep *= wt2;
0120 #ifdef EDM_ML_DEBUG
0121 edm::LogVerbatim("HGCSim") << "HGCalSD: weights= " << weight_ << ":" << wt1 << ":" << wt2 << " Total weight "
0122 << weight_ * wt1 * wt2 << " deStep: " << aStep->GetTotalEnergyDeposit() << ":" << destep;
0123 #endif
0124 return destep;
0125 }
0126
0127 uint32_t HGCalSD::setDetUnitId(const G4Step* aStep) {
0128 const G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0129 const G4VTouchable* touch = preStepPoint->GetTouchable();
0130
0131 #ifdef EDM_ML_DEBUG
0132 edm::LogVerbatim("HGCSim") << "DepthsTop: " << touch->GetHistoryDepth() << ":" << levelT1_ << ":" << levelT2_;
0133 printDetectorLevels(touch);
0134 #endif
0135
0136 G4ThreeVector hitPoint = preStepPoint->GetPosition();
0137 float globalZ = touch->GetTranslation(0).z();
0138 int iz(globalZ > 0 ? 1 : -1);
0139
0140 int layer(0), module(-1), cell(-1);
0141 if ((geom_mode_ == HGCalGeometryMode::Hexagon8Module) || (geom_mode_ == HGCalGeometryMode::Hexagon8Cassette)) {
0142 if (useSimWt_ > 0) {
0143 layer = touch->GetReplicaNumber(2);
0144 module = touch->GetReplicaNumber(1);
0145 } else if (touch->GetHistoryDepth() > levelT2_) {
0146 layer = touch->GetReplicaNumber(4);
0147 module = touch->GetReplicaNumber(3);
0148 cell = touch->GetReplicaNumber(1);
0149 } else {
0150 layer = touch->GetReplicaNumber(3);
0151 module = touch->GetReplicaNumber(2);
0152 }
0153 #ifdef EDM_ML_DEBUG
0154 edm::LogVerbatim("HGCSim") << "DepthsTop: " << touch->GetHistoryDepth() << ":" << levelT1_ << ":" << levelT2_ << ":"
0155 << useSimWt_ << " name " << touch->GetVolume(0)->GetName() << " layer:module:cell "
0156 << layer << ":" << module << ":" << cell;
0157 printDetectorLevels(touch);
0158 #endif
0159 } else if ((touch->GetHistoryDepth() == levelT1_) || (touch->GetHistoryDepth() == levelT2_)) {
0160 layer = touch->GetReplicaNumber(0);
0161 #ifdef EDM_ML_DEBUG
0162 edm::LogVerbatim("HGCSim") << "DepthsTop: " << touch->GetHistoryDepth() << ":" << levelT1_ << ":" << levelT2_
0163 << " name " << touch->GetVolume(0)->GetName() << " layer:module:cell " << layer << ":"
0164 << module << ":" << cell;
0165 #endif
0166 } else {
0167 layer = touch->GetReplicaNumber(3);
0168 module = touch->GetReplicaNumber(2);
0169 cell = touch->GetReplicaNumber(1);
0170 #ifdef EDM_ML_DEBUG
0171 edm::LogVerbatim("HGCSim") << "DepthsInside: " << touch->GetHistoryDepth() << " name "
0172 << touch->GetVolume(0)->GetName() << " layer:module:cell " << layer << ":" << module
0173 << ":" << cell;
0174 #endif
0175 }
0176 #ifdef EDM_ML_DEBUG
0177 G4Material* mat = aStep->GetPreStepPoint()->GetMaterial();
0178 edm::LogVerbatim("HGCSim") << "Depths: " << touch->GetHistoryDepth() << " name " << touch->GetVolume(0)->GetName()
0179 << ":" << touch->GetReplicaNumber(0) << " " << touch->GetVolume(1)->GetName() << ":"
0180 << touch->GetReplicaNumber(1) << " " << touch->GetVolume(2)->GetName() << ":"
0181 << touch->GetReplicaNumber(2) << " " << touch->GetVolume(3)->GetName() << ":"
0182 << touch->GetReplicaNumber(3) << " " << touch->GetVolume(4)->GetName() << ":"
0183 << touch->GetReplicaNumber(4) << " "
0184 << " layer:module:cell " << layer << ":" << module << ":" << cell << " Material "
0185 << mat->GetName() << ":" << mat->GetRadlen();
0186 #endif
0187
0188 if (aStep->GetPreStepPoint()->GetMaterial()->GetRadlen() > 100000.)
0189 return 0;
0190
0191 uint32_t id = setDetUnitId(layer, module, cell, iz, hitPoint);
0192 if (rejectMB_ && id != 0) {
0193 auto uv = HGCSiliconDetId(id).waferUV();
0194 #ifdef EDM_ML_DEBUG
0195 edm::LogVerbatim("HGCSim") << "ID " << std::hex << id << std::dec << " " << HGCSiliconDetId(id);
0196 #endif
0197 if (mouseBite_->exclude(hitPoint, iz, uv.first, uv.second)) {
0198 id = 0;
0199 #ifdef EDM_ML_DEBUG
0200 edm::LogVerbatim("HGCSim") << "Rejected by mousebite cutoff *****";
0201 #endif
0202 }
0203 }
0204 #ifdef EDM_ML_DEBUG
0205 if (id != 0)
0206 edm::LogVerbatim("HGCSim") << HGCSiliconDetId(id);
0207 #endif
0208 if ((id != 0) && checkID_) {
0209 HGCSiliconDetId hid1(id);
0210 bool cshift = (hgcons_->cassetteShiftSilicon(hid1.layer(), hid1.waferU(), hid1.waferV()));
0211 std::string_view pid = (cshift ? "HGCSim" : "HGCalSim");
0212 auto xy = hgcons_->locateCell(hid1, false);
0213 double xx = (hid1.zside() > 0) ? xy.first : -xy.first;
0214 double dx = xx - (hitPoint.x() / CLHEP::cm);
0215 double dy = xy.second - (hitPoint.y() / CLHEP::cm);
0216 double diff = std::sqrt(dx * dx + dy * dy);
0217 constexpr double tol = 2.0;
0218 bool valid1 = hgcons_->isValidHex8(hid1.layer(), hid1.waferU(), hid1.waferV(), hid1.cellU(), hid1.cellV(), true);
0219 if ((diff > tol) || (!valid1))
0220 pid = "HGCalError";
0221 auto partn = hgcons_->waferTypeRotation(hid1.layer(), hid1.waferU(), hid1.waferV(), false, false);
0222 edm::LogVerbatim(pid) << "CheckID " << HGCSiliconDetId(id) << " input position: (" << hitPoint.x() / CLHEP::cm
0223 << ", " << hitPoint.y() / CLHEP::cm << "); position from ID (" << xx << ", " << xy.second
0224 << ") distance " << diff << " Valid " << valid1 << " Wafer type|rotation " << partn.first
0225 << ":" << partn.second << " CassetteShift " << cshift;
0226 }
0227 return id;
0228 }
0229
0230 void HGCalSD::update(const BeginOfJob* job) {
0231 if (hgcons_ != nullptr) {
0232 geom_mode_ = hgcons_->geomMode();
0233 slopeMin_ = hgcons_->minSlope();
0234 levelT1_ = hgcons_->levelTop(0);
0235 levelT2_ = hgcons_->levelTop(1);
0236 useSimWt_ = hgcons_->getParameter()->useSimWt_;
0237 double waferSize = hgcons_->waferSize(false);
0238 double mouseBite = hgcons_->mouseBite(false);
0239 mouseBiteCut_ = waferSize * tan30deg_ - mouseBite;
0240 #ifdef EDM_ML_DEBUG
0241 edm::LogVerbatim("HGCSim") << "HGCalSD::Initialized with mode " << geom_mode_ << " Slope cut " << slopeMin_
0242 << " top Level " << levelT1_ << ":" << levelT2_ << " useSimWt " << useSimWt_ << " wafer "
0243 << waferSize << ":" << mouseBite;
0244 #endif
0245
0246 numberingScheme_ = std::make_unique<HGCalNumberingScheme>(*hgcons_, mydet_, nameX_, missingFile_);
0247 if (rejectMB_)
0248 mouseBite_ = std::make_unique<HGCMouseBite>(*hgcons_, angles_, mouseBiteCut_, waferRot_);
0249 } else {
0250 throw cms::Exception("Unknown", "HGCalSD") << "Cannot find HGCalDDDConstants for " << nameX_ << "\n";
0251 }
0252 }
0253
0254 void HGCalSD::initRun() {}
0255
0256 bool HGCalSD::filterHit(CaloG4Hit* aHit, double time) {
0257 return ((time <= tmaxHit) && (aHit->getEnergyDeposit() > eminHit_));
0258 }
0259
0260 uint32_t HGCalSD::setDetUnitId(int layer, int module, int cell, int iz, G4ThreeVector& pos) {
0261 uint32_t id = numberingScheme_ ? numberingScheme_->getUnitID(layer, module, cell, iz, pos, weight_) : 0;
0262 if (cornerMinMask_ > 2) {
0263 if (hgcons_->maskCell(DetId(id), cornerMinMask_)) {
0264 id = 0;
0265 ignoreRejection();
0266 }
0267 }
0268 if (hgcons_->waferHexagon8File() || (id == 0))
0269 ignoreRejection();
0270 return id;
0271 }
0272
0273 bool HGCalSD::isItinFidVolume(const G4ThreeVector& pos) {
0274 if (fiducialCut_) {
0275 return (hgcons_->distFromEdgeHex(pos.x(), pos.y(), pos.z()) > distanceFromEdge_);
0276 } else {
0277 return true;
0278 }
0279 }