File indexing completed on 2023-01-23 02:42:52
0001
0002
0003
0004
0005
0006
0007 #include "DataFormats/Math/interface/FastMath.h"
0008 #include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
0009 #include "SimG4CMS/Calo/interface/HGCScintSD.h"
0010 #include "SimG4CMS/Calo/interface/CaloSimUtils.h"
0011 #include "SimG4Core/Notification/interface/TrackInformation.h"
0012 #include "SimDataFormats/CaloTest/interface/HGCalTestNumbering.h"
0013 #include "FWCore/ParameterSet/interface/FileInPath.h"
0014 #include "FWCore/Utilities/interface/Exception.h"
0015 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0016 #include "Geometry/HGCalCommonData/interface/HGCalGeometryMode.h"
0017 #include "G4LogicalVolumeStore.hh"
0018 #include "G4LogicalVolume.hh"
0019 #include "G4Step.hh"
0020 #include "G4Track.hh"
0021 #include "G4ParticleTable.hh"
0022 #include "G4VProcess.hh"
0023 #include "G4Trap.hh"
0024
0025 #include <fstream>
0026 #include <iomanip>
0027 #include <iostream>
0028 #include <memory>
0029
0030
0031
0032 HGCScintSD::HGCScintSD(const std::string& name,
0033 const HGCalDDDConstants* hgc,
0034 const SensitiveDetectorCatalog& clg,
0035 edm::ParameterSet const& p,
0036 const SimTrackManager* manager)
0037 : CaloSD(name,
0038 clg,
0039 p,
0040 manager,
0041 (float)(p.getParameter<edm::ParameterSet>("HGCSD").getParameter<double>("TimeSliceUnit")),
0042 p.getParameter<edm::ParameterSet>("HGCSD").getParameter<bool>("IgnoreTrackID")),
0043 hgcons_(hgc),
0044 slopeMin_(0),
0045 levelT1_(99),
0046 levelT2_(99),
0047 firstLayer_(0) {
0048 numberingScheme_.reset(nullptr);
0049
0050 edm::ParameterSet m_HGC = p.getParameter<edm::ParameterSet>("HGCScintSD");
0051 eminHit_ = m_HGC.getParameter<double>("EminHit") * CLHEP::MeV;
0052 fiducialCut_ = m_HGC.getParameter<bool>("FiducialCut");
0053 distanceFromEdge_ = m_HGC.getParameter<double>("DistanceFromEdge");
0054 useBirk_ = m_HGC.getParameter<bool>("UseBirkLaw");
0055 birk1_ = m_HGC.getParameter<double>("BirkC1") * (CLHEP::g / (CLHEP::MeV * CLHEP::cm2));
0056 birk2_ = m_HGC.getParameter<double>("BirkC2");
0057 birk3_ = m_HGC.getParameter<double>("BirkC3");
0058 storeAllG4Hits_ = m_HGC.getParameter<bool>("StoreAllG4Hits");
0059 checkID_ = m_HGC.getUntrackedParameter<bool>("CheckID");
0060 fileName_ = m_HGC.getUntrackedParameter<std::string>("TileFileName");
0061
0062 if (storeAllG4Hits_) {
0063 setUseMap(false);
0064 setNumberCheckedHits(0);
0065 }
0066
0067
0068 G4String myName = name;
0069 mydet_ = DetId::Forward;
0070 nameX_ = "HGCal";
0071 if (myName.find("HitsHEback") != std::string::npos) {
0072 mydet_ = DetId::HGCalHSc;
0073 nameX_ = "HGCalHEScintillatorSensitive";
0074 }
0075
0076 #ifdef EDM_ML_DEBUG
0077 edm::LogVerbatim("HGCSim") << "**************************************************"
0078 << "\n"
0079 << "* *"
0080 << "\n"
0081 << "* Constructing a HGCScintSD with name " << name << "\n"
0082 << "* *"
0083 << "\n"
0084 << "**************************************************";
0085 #endif
0086 edm::LogVerbatim("HGCSim") << "HGCScintSD:: Threshold for storing hits: " << eminHit_ << " for " << nameX_
0087 << " detector " << mydet_ << " File " << fileName_;
0088 edm::LogVerbatim("HGCSim") << "Flag for storing individual Geant4 Hits " << storeAllG4Hits_;
0089 edm::LogVerbatim("HGCSim") << "Fiducial volume cut with cut from eta/phi "
0090 << "boundary " << fiducialCut_ << " at " << distanceFromEdge_;
0091 edm::LogVerbatim("HGCSim") << "Use of Birks law is set to " << useBirk_
0092 << " with three constants kB = " << birk1_ << ", C1 = " << birk2_ << ", C2 = " << birk3_;
0093
0094 if (!fileName_.empty()) {
0095 edm::FileInPath filetmp("SimG4CMS/Calo/data/" + fileName_);
0096 std::string fileName = filetmp.fullPath();
0097 std::ifstream fInput(fileName.c_str());
0098 if (!fInput.good()) {
0099 edm::LogVerbatim("HGCSim") << "Cannot open file " << fileName;
0100 } else {
0101 char buffer[80];
0102 while (fInput.getline(buffer, 80)) {
0103 std::vector<std::string> items = CaloSimUtils::splitString(std::string(buffer));
0104 if (items.size() > 2) {
0105 int layer = std::atoi(items[0].c_str());
0106 int ring = std::atoi(items[1].c_str());
0107 int phi = std::atoi(items[2].c_str());
0108 tiles_.emplace_back(HGCalTileIndex::tileIndex(layer, ring, phi));
0109 }
0110 }
0111 edm::LogVerbatim("HGCSim") << "Reads in " << tiles_.size() << " tile information from " << fileName_;
0112 fInput.close();
0113 }
0114 }
0115 }
0116
0117 double HGCScintSD::getEnergyDeposit(const G4Step* aStep) {
0118 double r = aStep->GetPreStepPoint()->GetPosition().perp();
0119 double z = std::abs(aStep->GetPreStepPoint()->GetPosition().z());
0120 #ifdef EDM_ML_DEBUG
0121 G4int parCode = aStep->GetTrack()->GetDefinition()->GetPDGEncoding();
0122 G4String parName = aStep->GetTrack()->GetDefinition()->GetParticleName();
0123 G4LogicalVolume* lv = aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
0124 edm::LogVerbatim("HGCSim") << "HGCScintSD: Hit from standard path from " << lv->GetName() << " for Track "
0125 << aStep->GetTrack()->GetTrackID() << " (" << parCode << ":" << parName << ") R = " << r
0126 << " Z = " << z << " slope = " << r / z << ":" << slopeMin_;
0127 #endif
0128
0129 if (r < z * slopeMin_) {
0130 #ifdef EDM_ML_DEBUG
0131 edm::LogVerbatim("HGCSim") << "HGCScintSD: Fiducial Volume cut";
0132 #endif
0133 return 0.0;
0134 }
0135
0136 double wt1 = getResponseWt(aStep->GetTrack());
0137 double wt2 = aStep->GetTrack()->GetWeight();
0138 double wt3 = (useBirk_ ? getAttenuation(aStep, birk1_, birk2_, birk3_) : 1.0);
0139 double destep = weight_ * wt1 * wt3 * (aStep->GetTotalEnergyDeposit());
0140 if (wt2 > 0)
0141 destep *= wt2;
0142 #ifdef EDM_ML_DEBUG
0143 edm::LogVerbatim("HGCalSim") << "HGCScintSD: weights= " << weight_ << ":" << wt1 << ":" << wt2 << ":" << wt3
0144 << " Total weight " << weight_ * wt1 * wt2 * wt3
0145 << " deStep: " << aStep->GetTotalEnergyDeposit() << ":" << destep;
0146 #endif
0147 return destep;
0148 }
0149
0150 uint32_t HGCScintSD::setDetUnitId(const G4Step* aStep) {
0151 const G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0152 const G4VTouchable* touch = preStepPoint->GetTouchable();
0153
0154 #ifdef EDM_ML_DEBUG
0155 edm::LogVerbatim("HGCSim") << "DepthsTop: " << touch->GetHistoryDepth() << ":" << levelT1_ << ":" << levelT2_;
0156 printDetectorLevels(touch);
0157 #endif
0158
0159 G4ThreeVector hitPoint = preStepPoint->GetPosition();
0160 float globalZ = touch->GetTranslation(0).z();
0161 int iz(globalZ > 0 ? 1 : -1);
0162
0163 int layer(0), module(-1), cell(-1);
0164 if ((geom_mode_ == HGCalGeometryMode::TrapezoidModule) || (geom_mode_ == HGCalGeometryMode::TrapezoidCassette)) {
0165 layer = touch->GetReplicaNumber(1);
0166 } else if ((touch->GetHistoryDepth() == levelT1_) || (touch->GetHistoryDepth() == levelT2_)) {
0167 layer = touch->GetReplicaNumber(0);
0168 #ifdef EDM_ML_DEBUG
0169 edm::LogVerbatim("HGCSim") << "DepthsTop: " << touch->GetHistoryDepth() << ":" << levelT1_ << ":" << levelT2_
0170 << " name " << touch->GetVolume(0)->GetName() << " layer:module:cell " << layer << ":"
0171 << module << ":" << cell;
0172 #endif
0173 } else {
0174 layer = touch->GetReplicaNumber(3);
0175 module = touch->GetReplicaNumber(2);
0176 cell = touch->GetReplicaNumber(1);
0177 #ifdef EDM_ML_DEBUG
0178 edm::LogVerbatim("HGCSim") << "DepthsInside: " << touch->GetHistoryDepth() << " name "
0179 << touch->GetVolume(0)->GetName() << " layer:module:cell " << layer << ":" << module
0180 << ":" << cell;
0181 #endif
0182 }
0183 #ifdef EDM_ML_DEBUG
0184 G4Material* mat = aStep->GetPreStepPoint()->GetMaterial();
0185 edm::LogVerbatim("HGCSim") << "Depths: " << touch->GetHistoryDepth() << " name " << touch->GetVolume(0)->GetName()
0186 << ":" << touch->GetReplicaNumber(0) << " " << touch->GetVolume(1)->GetName() << ":"
0187 << touch->GetReplicaNumber(1) << " " << touch->GetVolume(2)->GetName() << ":"
0188 << touch->GetReplicaNumber(2) << " " << touch->GetVolume(3)->GetName() << ":"
0189 << touch->GetReplicaNumber(3) << " " << touch->GetVolume(4)->GetName() << ":"
0190 << touch->GetReplicaNumber(4) << " "
0191 << " layer:module:cell " << layer << ":" << module << ":" << cell << " Material "
0192 << mat->GetName() << ":" << mat->GetRadlen();
0193 #endif
0194
0195 if (aStep->GetPreStepPoint()->GetMaterial()->GetRadlen() > 100000.)
0196 return 0;
0197
0198 uint32_t id = setDetUnitId(layer, module, cell, iz, hitPoint);
0199 bool debug(false);
0200 if (!tiles_.empty()) {
0201 HGCScintillatorDetId hid(id);
0202 int indx = HGCalTileIndex::tileIndex(firstLayer_ + hid.layer(), hid.ring(), hid.iphi());
0203 if (std::find(tiles_.begin(), tiles_.end(), indx) != tiles_.end())
0204 debug = true;
0205 }
0206 if (debug)
0207 edm::LogVerbatim("HGCSim") << "Layer:module:cell:iz " << layer << ":" << module << ":" << cell << ":" << iz
0208 << " Point (" << hitPoint.x() << ", " << hitPoint.y() << ", " << hitPoint.z() << ") "
0209 << HGCScintillatorDetId(id);
0210
0211 if (!isItinFidVolume(hitPoint)) {
0212 #ifdef EDM_ML_DEBUG
0213 edm::LogVerbatim("HGCSim") << "ID " << std::hex << id << std::dec << " " << HGCScintillatorDetId(id)
0214 << " is rejected by fiducilal volume cut";
0215 #endif
0216 id = 0;
0217 }
0218 if ((id != 0) && checkID_) {
0219 HGCScintillatorDetId hid1(id);
0220 std::string_view pid = ((hgcons_->cassetteShiftScintillator(hid1.layer(), hid1.iphi())) ? "HGCSim" : "HGCalSim");
0221 auto xy = hgcons_->locateCell(HGCScintillatorDetId(id), false);
0222 double dx = xy.first - (hitPoint.x() / CLHEP::cm);
0223 double dy = xy.second - (hitPoint.y() / CLHEP::cm);
0224 double diff = std::sqrt(dx * dx + dy * dy);
0225 constexpr double tol = 10.0;
0226 bool valid = hgcons_->isValidTrap(hid1.zside(), hid1.layer(), hid1.ring(), hid1.iphi());
0227 if ((diff > tol) || (!valid))
0228 pid = "HGCalError";
0229 edm::LogVerbatim(pid) << "CheckID " << HGCScintillatorDetId(id) << " input position: (" << hitPoint.x() / CLHEP::cm
0230 << ", " << hitPoint.y() / CLHEP::cm << "); position from ID (" << xy.first << ", "
0231 << xy.second << ") distance " << diff << " Valid " << valid
0232 << " Rho = " << hitPoint.perp() / CLHEP::cm;
0233 }
0234 return id;
0235 }
0236
0237 void HGCScintSD::update(const BeginOfJob* job) {
0238 if (hgcons_ != nullptr) {
0239 geom_mode_ = hgcons_->geomMode();
0240 slopeMin_ = hgcons_->minSlope();
0241 levelT1_ = hgcons_->levelTop(0);
0242 levelT2_ = hgcons_->levelTop(1);
0243 firstLayer_ = hgcons_->firstLayer() - 1;
0244 #ifdef EDM_ML_DEBUG
0245 edm::LogVerbatim("HGCSim") << "HGCScintSD::Initialized with mode " << geom_mode_ << " Slope cut " << slopeMin_
0246 << " top Level " << levelT1_ << ":" << levelT2_ << " FirstLayer " << firstLayer_;
0247 #endif
0248
0249 numberingScheme_ = std::make_unique<HGCalNumberingScheme>(*hgcons_, mydet_, nameX_, fileName_);
0250 } else {
0251 throw cms::Exception("Unknown", "HGCScintSD") << "Cannot find HGCalDDDConstants for " << nameX_ << "\n";
0252 }
0253 }
0254
0255 void HGCScintSD::initRun() {}
0256
0257 bool HGCScintSD::filterHit(CaloG4Hit* aHit, double time) {
0258 return ((time <= tmaxHit) && (aHit->getEnergyDeposit() > eminHit_));
0259 }
0260
0261 uint32_t HGCScintSD::setDetUnitId(int layer, int module, int cell, int iz, G4ThreeVector& pos) {
0262 uint32_t id = numberingScheme_ ? numberingScheme_->getUnitID(layer, module, cell, iz, pos, weight_) : 0;
0263 return id;
0264 }
0265
0266 bool HGCScintSD::isItinFidVolume(const G4ThreeVector& pos) {
0267 if (fiducialCut_) {
0268 return (hgcons_->distFromEdgeTrap(pos.x(), pos.y(), pos.z()) > distanceFromEdge_);
0269 } else {
0270 return true;
0271 }
0272 }