File indexing completed on 2021-07-27 14:07:57
0001
0002
0003
0004
0005
0006 #include "DataFormats/Math/interface/FastMath.h"
0007 #include "DataFormats/ForwardDetId/interface/HFNoseDetId.h"
0008 #include "SimG4CMS/Calo/interface/HFNoseSD.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 <iostream>
0022 #include <fstream>
0023 #include <iomanip>
0024
0025
0026
0027 HFNoseSD::HFNoseSD(const std::string& name,
0028 const HGCalDDDConstants* hgc,
0029 const SensitiveDetectorCatalog& clg,
0030 edm::ParameterSet const& p,
0031 const SimTrackManager* manager)
0032 : CaloSD(name,
0033 clg,
0034 p,
0035 manager,
0036 (float)(p.getParameter<edm::ParameterSet>("HFNoseSD").getParameter<double>("TimeSliceUnit")),
0037 p.getParameter<edm::ParameterSet>("HFNoseSD").getParameter<bool>("IgnoreTrackID")),
0038 hgcons_(hgc),
0039 slopeMin_(0),
0040 levelT1_(99),
0041 levelT2_(99),
0042 tan30deg_(std::tan(30.0 * CLHEP::deg)) {
0043 numberingScheme_.reset(nullptr);
0044 mouseBite_.reset(nullptr);
0045
0046 edm::ParameterSet m_HFN = p.getParameter<edm::ParameterSet>("HFNoseSD");
0047 eminHit_ = m_HFN.getParameter<double>("EminHit") * CLHEP::MeV;
0048 fiducialCut_ = m_HFN.getParameter<bool>("FiducialCut");
0049 distanceFromEdge_ = m_HFN.getParameter<double>("DistanceFromEdge");
0050 storeAllG4Hits_ = m_HFN.getParameter<bool>("StoreAllG4Hits");
0051 rejectMB_ = m_HFN.getParameter<bool>("RejectMouseBite");
0052 waferRot_ = m_HFN.getParameter<bool>("RotatedWafer");
0053 cornerMinMask_ = m_HFN.getParameter<int>("CornerMinMask");
0054 angles_ = m_HFN.getUntrackedParameter<std::vector<double>>("WaferAngles");
0055
0056 nameX_ = ((name.find("HFNoseHits") != std::string::npos) ? "HGCalHFNoseSensitive" : "HFNoseSensitive");
0057
0058 if (storeAllG4Hits_) {
0059 setUseMap(false);
0060 setNumberCheckedHits(0);
0061 }
0062
0063 #ifdef EDM_ML_DEBUG
0064 edm::LogVerbatim("HFNSim") << "**************************************************"
0065 << "\n"
0066 << "* *"
0067 << "\n"
0068 << "* Constructing a HFNoseSD with name " << name << "\n"
0069 << "* *"
0070 << "\n"
0071 << "**************************************************";
0072 #endif
0073 edm::LogVerbatim("HFNSim") << "HFNoseSD:: Threshold for storing hits: " << eminHit_ << " for " << name;
0074 edm::LogVerbatim("HFNSim") << "Flag for storing individual Geant4 Hits " << storeAllG4Hits_;
0075 edm::LogVerbatim("HFNSim") << "Fiducial volume cut with cut from eta/phi "
0076 << "boundary " << fiducialCut_ << " at " << distanceFromEdge_;
0077 edm::LogVerbatim("HFNSim") << "Reject MosueBite Flag: " << rejectMB_ << " cuts along " << angles_.size()
0078 << " axes: " << angles_[0] << ", " << angles_[1];
0079 }
0080
0081 double HFNoseSD::getEnergyDeposit(const G4Step* aStep) {
0082 double r = aStep->GetPreStepPoint()->GetPosition().perp();
0083 double z = std::abs(aStep->GetPreStepPoint()->GetPosition().z());
0084 #ifdef EDM_ML_DEBUG
0085 G4int parCode = aStep->GetTrack()->GetDefinition()->GetPDGEncoding();
0086 G4String parName = aStep->GetTrack()->GetDefinition()->GetParticleName();
0087 G4LogicalVolume* lv = aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
0088 edm::LogVerbatim("HFNSim") << "HFNoseSD: Hit from standard path from " << lv->GetName() << " for Track "
0089 << aStep->GetTrack()->GetTrackID() << " (" << parCode << ":" << parName << ") R = " << r
0090 << " Z = " << z << " slope = " << r / z << ":" << slopeMin_;
0091 #endif
0092
0093 if (r < z * slopeMin_) {
0094 #ifdef EDM_ML_DEBUG
0095 edm::LogVerbatim("HFNSim") << "HFNoseSD: Fiducial Volume cut";
0096 #endif
0097 return 0.0;
0098 }
0099
0100 double wt1 = getResponseWt(aStep->GetTrack());
0101 double wt2 = aStep->GetTrack()->GetWeight();
0102 double destep = weight_ * wt1 * (aStep->GetTotalEnergyDeposit());
0103 if (wt2 > 0)
0104 destep *= wt2;
0105 #ifdef EDM_ML_DEBUG
0106 edm::LogVerbatim("HFNSim") << "HFNoseSD: weights= " << weight_ << ":" << wt1 << ":" << wt2 << " Total weight "
0107 << weight_ * wt1 * wt2 << " deStep: " << aStep->GetTotalEnergyDeposit() << ":" << destep;
0108 #endif
0109 return destep;
0110 }
0111
0112 uint32_t HFNoseSD::setDetUnitId(const G4Step* aStep) {
0113 const G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0114 const G4VTouchable* touch = preStepPoint->GetTouchable();
0115
0116
0117 G4ThreeVector hitPoint = preStepPoint->GetPosition();
0118 float globalZ = touch->GetTranslation(0).z();
0119 int iz(globalZ > 0 ? 1 : -1);
0120
0121 int layer, module, cell;
0122 if ((touch->GetHistoryDepth() == levelT1_) || (touch->GetHistoryDepth() == levelT2_)) {
0123 layer = touch->GetReplicaNumber(0);
0124 module = -1;
0125 cell = -1;
0126 #ifdef EDM_ML_DEBUG
0127 edm::LogVerbatim("HFNSim") << "DepthsTop: " << touch->GetHistoryDepth() << ":" << levelT1_ << ":" << levelT2_
0128 << " name " << touch->GetVolume(0)->GetName() << " layer:module:cell " << layer << ":"
0129 << module << ":" << cell;
0130 #endif
0131 } else {
0132 layer = touch->GetReplicaNumber(3);
0133 module = touch->GetReplicaNumber(2);
0134 cell = touch->GetReplicaNumber(1);
0135 #ifdef EDM_ML_DEBUG
0136 edm::LogVerbatim("HFNSim") << "DepthsInside: " << touch->GetHistoryDepth() << " name "
0137 << touch->GetVolume(0)->GetName() << " layer:module:cell " << layer << ":" << module
0138 << ":" << cell;
0139 #endif
0140 }
0141 #ifdef EDM_ML_DEBUG
0142 G4Material* mat = aStep->GetPreStepPoint()->GetMaterial();
0143 edm::LogVerbatim("HFNSim") << "Depths: " << touch->GetHistoryDepth() << " name " << touch->GetVolume(0)->GetName()
0144 << ":" << touch->GetReplicaNumber(0) << " " << touch->GetVolume(1)->GetName() << ":"
0145 << touch->GetReplicaNumber(1) << " " << touch->GetVolume(2)->GetName() << ":"
0146 << touch->GetReplicaNumber(2) << " " << touch->GetVolume(3)->GetName() << ":"
0147 << touch->GetReplicaNumber(3) << " " << touch->GetVolume(4)->GetName() << ":"
0148 << touch->GetReplicaNumber(4) << " "
0149 << " layer:module:cell " << layer << ":" << module << ":" << cell << " Material "
0150 << mat->GetName() << ":" << mat->GetRadlen();
0151 #endif
0152
0153 if (aStep->GetPreStepPoint()->GetMaterial()->GetRadlen() > 100000.)
0154 return 0;
0155
0156 uint32_t id = setDetUnitId(layer, module, cell, iz, hitPoint);
0157 if (rejectMB_ && id != 0) {
0158 auto uv = HFNoseDetId(id).waferUV();
0159 #ifdef EDM_ML_DEBUG
0160 edm::LogVerbatim("HFNSim") << "ID " << std::hex << id << std::dec << " " << HFNoseDetId(id);
0161 #endif
0162 if (mouseBite_->exclude(hitPoint, iz, uv.first, uv.second)) {
0163 id = 0;
0164 #ifdef EDM_ML_DEBUG
0165 edm::LogVerbatim("HFNSim") << "Rejected by mousebite cutoff *****";
0166 #endif
0167 }
0168 }
0169 return id;
0170 }
0171
0172 void HFNoseSD::update(const BeginOfJob* job) {
0173 if (hgcons_ != nullptr) {
0174 geom_mode_ = hgcons_->geomMode();
0175 slopeMin_ = hgcons_->minSlope();
0176 levelT1_ = hgcons_->levelTop(0);
0177 levelT2_ = hgcons_->levelTop(1);
0178 double waferSize = hgcons_->waferSize(false);
0179 double mouseBite = hgcons_->mouseBite(false);
0180 mouseBiteCut_ = waferSize * tan30deg_ - mouseBite;
0181 #ifdef EDM_ML_DEBUG
0182 edm::LogVerbatim("HFNSim") << "HFNoseSD::Initialized with mode " << geom_mode_ << " Slope cut " << slopeMin_
0183 << " top Level " << levelT1_ << ":" << levelT2_ << " wafer " << waferSize << ":"
0184 << mouseBite;
0185 #endif
0186
0187 numberingScheme_ = std::make_unique<HFNoseNumberingScheme>(*hgcons_);
0188 if (rejectMB_)
0189 mouseBite_ = std::make_unique<HGCMouseBite>(*hgcons_, angles_, mouseBiteCut_, waferRot_);
0190 } else {
0191 throw cms::Exception("Unknown", "HFNoseSD") << "Cannot find HGCalDDDConstants for " << nameX_ << "\n";
0192 }
0193 }
0194
0195 void HFNoseSD::initRun() {}
0196
0197 bool HFNoseSD::filterHit(CaloG4Hit* aHit, double time) {
0198 return ((time <= tmaxHit) && (aHit->getEnergyDeposit() > eminHit_));
0199 }
0200
0201 uint32_t HFNoseSD::setDetUnitId(int layer, int module, int cell, int iz, G4ThreeVector& pos) {
0202 uint32_t id = numberingScheme_ ? numberingScheme_->getUnitID(layer, module, cell, iz, pos, weight_) : 0;
0203 if (cornerMinMask_ > 2) {
0204 if (hgcons_->maskCell(DetId(id), cornerMinMask_))
0205 id = 0;
0206 }
0207 return id;
0208 }
0209
0210 bool HFNoseSD::isItinFidVolume(const G4ThreeVector& pos) {
0211 if (fiducialCut_) {
0212 return (hgcons_->distFromEdgeHex(pos.x(), pos.y(), pos.z()) > distanceFromEdge_);
0213 } else {
0214 return true;
0215 }
0216 }