File indexing completed on 2021-07-27 14:07:58
0001
0002
0003
0004
0005
0006 #include "DataFormats/Math/interface/FastMath.h"
0007
0008 #include "SimG4CMS/Calo/interface/HGCSD.h"
0009 #include "SimG4Core/Notification/interface/TrackInformation.h"
0010 #include "SimDataFormats/CaloTest/interface/HGCalTestNumbering.h"
0011 #include "FWCore/Utilities/interface/Exception.h"
0012 #include "FWCore/ServiceRegistry/interface/Service.h"
0013 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0014 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0015 #include "Geometry/HGCalCommonData/interface/HGCalGeometryMode.h"
0016 #include "G4LogicalVolumeStore.hh"
0017 #include "G4LogicalVolume.hh"
0018 #include "G4Step.hh"
0019 #include "G4Track.hh"
0020 #include "G4ParticleTable.hh"
0021 #include "G4VProcess.hh"
0022 #include "G4Trap.hh"
0023
0024 #include <fstream>
0025 #include <iomanip>
0026 #include <iostream>
0027 #include <memory>
0028
0029
0030
0031
0032 HGCSD::HGCSD(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 levelT_(99),
0046 tree_(nullptr) {
0047 numberingScheme_.reset(nullptr);
0048 mouseBite_.reset(nullptr);
0049
0050 edm::ParameterSet m_HGC = p.getParameter<edm::ParameterSet>("HGCSD");
0051 eminHit_ = m_HGC.getParameter<double>("EminHit") * CLHEP::MeV;
0052 storeAllG4Hits_ = m_HGC.getParameter<bool>("StoreAllG4Hits");
0053 rejectMB_ = m_HGC.getParameter<bool>("RejectMouseBite");
0054 waferRot_ = m_HGC.getParameter<bool>("RotatedWafer");
0055 angles_ = m_HGC.getUntrackedParameter<std::vector<double>>("WaferAngles");
0056 double waferSize = m_HGC.getUntrackedParameter<double>("WaferSize") * CLHEP::mm;
0057 double mouseBite = m_HGC.getUntrackedParameter<double>("MouseBite") * CLHEP::mm;
0058 mouseBiteCut_ = waferSize * tan(30.0 * CLHEP::deg) - mouseBite;
0059
0060 if (storeAllG4Hits_) {
0061 setUseMap(false);
0062 setNumberCheckedHits(0);
0063 }
0064
0065 G4String myName = name;
0066 myFwdSubdet_ = ForwardSubdetector::ForwardEmpty;
0067 nameX_ = "HGCal";
0068 if (myName.find("HitsEE") != std::string::npos) {
0069 myFwdSubdet_ = ForwardSubdetector::HGCEE;
0070 nameX_ = "HGCalEESensitive";
0071 } else if (myName.find("HitsHEfront") != std::string::npos) {
0072 myFwdSubdet_ = ForwardSubdetector::HGCHEF;
0073 nameX_ = "HGCalHESiliconSensitive";
0074 } else if (myName.find("HitsHEback") != std::string::npos) {
0075 myFwdSubdet_ = ForwardSubdetector::HGCHEB;
0076 nameX_ = "HGCalHEScintillatorSensitive";
0077 }
0078
0079 #ifdef EDM_ML_DEBUG
0080 edm::LogVerbatim("HGCSim") << "**************************************************"
0081 << "\n"
0082 << "* *"
0083 << "\n"
0084 << "* Constructing a HGCSD with name " << name << "\n"
0085 << "* *"
0086 << "\n"
0087 << "**************************************************";
0088 #endif
0089 edm::LogVerbatim("HGCSim") << "HGCSD:: Threshold for storing hits: " << eminHit << " for " << nameX_ << " subdet "
0090 << myFwdSubdet_;
0091 edm::LogVerbatim("HGCSim") << "Flag for storing individual Geant4 Hits " << storeAllG4Hits_;
0092 edm::LogVerbatim("HGCSim") << "Reject MosueBite Flag: " << rejectMB_ << " Size of wafer " << waferSize
0093 << " Mouse Bite " << mouseBite << ":" << mouseBiteCut_ << " along " << angles_.size()
0094 << " axes";
0095 }
0096
0097 double HGCSD::getEnergyDeposit(const G4Step* aStep) {
0098 double r = aStep->GetPreStepPoint()->GetPosition().perp();
0099 double z = std::abs(aStep->GetPreStepPoint()->GetPosition().z());
0100
0101 #ifdef EDM_ML_DEBUG
0102 G4int parCode = aStep->GetTrack()->GetDefinition()->GetPDGEncoding();
0103 G4LogicalVolume* lv = aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
0104 edm::LogVerbatim("HGCSim") << "HGCSD: Hit from standard path from " << lv->GetName() << " for Track "
0105 << aStep->GetTrack()->GetTrackID() << " ("
0106 << aStep->GetTrack()->GetDefinition()->GetParticleName() << ":" << parCode << ") R = " << r
0107 << " Z = " << z << " slope = " << r / z << ":" << slopeMin_;
0108 #endif
0109
0110
0111 if (r < z * slopeMin_) {
0112 return 0.0;
0113 }
0114
0115 double wt1 = getResponseWt(aStep->GetTrack());
0116 double wt2 = aStep->GetTrack()->GetWeight();
0117 double destep = wt1 * aStep->GetTotalEnergyDeposit();
0118 if (wt2 > 0)
0119 destep *= wt2;
0120
0121 #ifdef plotDebug
0122 const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
0123 G4double tmptrackE = aStep->GetTrack()->GetKineticEnergy();
0124 G4int parCode = aStep->GetTrack()->GetDefinition()->GetPDGEncoding();
0125 G4double angle = (aStep->GetTrack()->GetMomentumDirection().theta()) / CLHEP::deg;
0126 G4int layer = ((touch->GetHistoryDepth() == levelT_) ? touch->GetReplicaNumber(0) : touch->GetReplicaNumber(2));
0127 G4int ilayer = (layer - 1) / 3;
0128 if (aStep->GetTotalEnergyDeposit() > 0) {
0129 t_Layer_.emplace_back(ilayer);
0130 t_Parcode_.emplace_back(parCode);
0131 t_dEStep1_.emplace_back(aStep->GetTotalEnergyDeposit());
0132 t_dEStep2_.emplace_back(destep);
0133 t_TrackE_.emplace_back(tmptrackE);
0134 t_Angle_.emplace_back(angle);
0135 }
0136 #endif
0137
0138 return destep;
0139 }
0140
0141 uint32_t HGCSD::setDetUnitId(const G4Step* aStep) {
0142 const G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0143 const G4VTouchable* touch = preStepPoint->GetTouchable();
0144
0145
0146 G4ThreeVector hitPoint = preStepPoint->GetPosition();
0147 float globalZ = touch->GetTranslation(0).z();
0148 int iz(globalZ > 0 ? 1 : -1);
0149
0150
0151 G4ThreeVector localpos = touch->GetHistory()->GetTopTransform().TransformPoint(hitPoint);
0152
0153
0154 ForwardSubdetector subdet = myFwdSubdet_;
0155
0156 int layer, module, cell;
0157 if (touch->GetHistoryDepth() == levelT_) {
0158 layer = touch->GetReplicaNumber(0);
0159 module = -1;
0160 cell = -1;
0161 #ifdef EDM_ML_DEBUG
0162 edm::LogVerbatim("HGCSim") << "Depths: " << touch->GetHistoryDepth() << " name " << touch->GetVolume(0)->GetName()
0163 << " layer:module:cell " << layer << ":" << module << ":" << cell;
0164 #endif
0165 } else {
0166 layer = touch->GetReplicaNumber(2);
0167 module = touch->GetReplicaNumber(1);
0168 cell = touch->GetReplicaNumber(0);
0169 }
0170 #ifdef EDM_ML_DEBUG
0171 const G4Material* mat = aStep->GetPreStepPoint()->GetMaterial();
0172 edm::LogVerbatim("HGCSim") << "Depths: " << touch->GetHistoryDepth() << " name " << touch->GetVolume(0)->GetName()
0173 << ":" << touch->GetReplicaNumber(0) << " " << touch->GetVolume(1)->GetName() << ":"
0174 << touch->GetReplicaNumber(1) << " " << touch->GetVolume(2)->GetName() << ":"
0175 << touch->GetReplicaNumber(2) << " layer:module:cell " << layer << ":" << module << ":"
0176 << cell << " Material " << mat->GetName() << ":" << mat->GetRadlen();
0177
0178
0179 #endif
0180
0181
0182 if (aStep->GetPreStepPoint()->GetMaterial()->GetRadlen() > 100000.)
0183 return 0;
0184
0185 uint32_t id = setDetUnitId(subdet, layer, module, cell, iz, localpos);
0186 if (rejectMB_ && id != 0) {
0187 int det, z, lay, wafer, type, ic;
0188 HGCalTestNumbering::unpackHexagonIndex(id, det, z, lay, wafer, type, ic);
0189 #ifdef EDM_ML_DEBUG
0190 edm::LogVerbatim("HGCSim") << "ID " << std::hex << id << std::dec << " Input " << subdet << ":" << layer << ":"
0191 << module << ":" << cell << ":" << iz << localpos.x() << ":" << localpos.y()
0192 << " Decode " << det << ":" << z << ":" << lay << ":" << wafer << ":" << type << ":"
0193 << ic;
0194 #endif
0195 if (mouseBite_->exclude(hitPoint, z, wafer, 0))
0196 id = 0;
0197 }
0198 return id;
0199 }
0200
0201 void HGCSD::update(const BeginOfJob* job) {
0202 if (hgcons_ != nullptr) {
0203 geom_mode_ = hgcons_->geomMode();
0204 slopeMin_ = hgcons_->minSlope();
0205 levelT_ = hgcons_->levelTop();
0206 numberingScheme_ = std::make_unique<HGCNumberingScheme>(*hgcons_, nameX_);
0207 if (rejectMB_)
0208 mouseBite_ = std::make_unique<HGCMouseBite>(*hgcons_, angles_, mouseBiteCut_, waferRot_);
0209 } else {
0210 edm::LogError("HGCSim") << "HGCSD : Cannot find HGCalDDDConstants for " << nameX_;
0211 throw cms::Exception("Unknown", "HGCSD") << "Cannot find HGCalDDDConstants for " << nameX_ << "\n";
0212 }
0213 #ifdef EDM_ML_DEBUG
0214 edm::LogVerbatim("HGCSim") << "HGCSD::Initialized with mode " << geom_mode_ << " Slope cut " << slopeMin_
0215 << " top Level " << levelT_;
0216 #endif
0217 }
0218
0219 void HGCSD::initRun() {
0220 #ifdef plotDebug
0221 edm::Service<TFileService> tfile;
0222 if (tfile.isAvailable()) {
0223 tree_ = tfile->make<TTree>("TreeHGCSD", "TreeHGCSD");
0224 tree_->Branch("EventID", &t_EventID_);
0225 tree_->Branch("Layer", &t_Layer_);
0226 tree_->Branch("ParticleCode", &t_Parcode_);
0227 tree_->Branch("dEStepOriginal", &t_dEStep1_);
0228 tree_->Branch("dEStepWeighted", &t_dEStep2_);
0229 tree_->Branch("TrackEnergy", &t_TrackE_);
0230 tree_->Branch("ThetaAngle", &t_Angle_);
0231 }
0232 #endif
0233 }
0234
0235 void HGCSD::initEvent(const BeginOfEvent* g4Event) {
0236 const G4Event* evt = (*g4Event)();
0237 t_EventID_ = evt->GetEventID();
0238 #ifdef plotDebug
0239 t_Layer_.clear();
0240 t_Parcode_.clear();
0241 t_dEStep1_.clear();
0242 t_dEStep2_.clear();
0243 t_TrackE_.clear();
0244 t_Angle_.clear();
0245 #endif
0246 }
0247
0248 void HGCSD::endEvent() {
0249 #ifdef plotDebug
0250 if (tree_)
0251 tree_->Fill();
0252 #endif
0253 }
0254
0255 bool HGCSD::filterHit(CaloG4Hit* aHit, double time) {
0256 return ((time <= tmaxHit) && (aHit->getEnergyDeposit() > eminHit_));
0257 }
0258
0259 uint32_t HGCSD::setDetUnitId(ForwardSubdetector& subdet, int layer, int module, int cell, int iz, G4ThreeVector& pos) {
0260 uint32_t id = numberingScheme_ ? numberingScheme_->getUnitID(subdet, layer, module, cell, iz, pos) : 0;
0261 return id;
0262 }