CMS 3D CMS Logo
HGCalSD.cc
Go to the documentation of this file.
1 // File: HGCalSD.cc
3 // Description: Sensitive Detector class for High Granularity Calorimeter
5 
6 #include "DataFormats/Math/interface/FastMath.h"
7 #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
8 #include "SimG4CMS/Calo/interface/HGCalSD.h"
9 #include "SimG4Core/Notification/interface/TrackInformation.h"
10 #include "FWCore/Utilities/interface/Exception.h"
11 #include "FWCore/Framework/interface/ESHandle.h"
12 #include "FWCore/Framework/interface/EventSetup.h"
13 #include "Geometry/Records/interface/IdealGeometryRecord.h"
14 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
15 #include "Geometry/HGCalCommonData/interface/HGCalGeometryMode.h"
16 #include "G4LogicalVolumeStore.hh"
17 #include "G4LogicalVolume.hh"
18 #include "G4Step.hh"
19 #include "G4Track.hh"
20 #include "G4ParticleTable.hh"
21 #include "G4VProcess.hh"
22 #include "G4Trap.hh"
23 
24 #include <fstream>
25 #include <iomanip>
26 #include <iostream>
27 #include <memory>
28 
29 //#define EDM_ML_DEBUG
30 
31 HGCalSD::HGCalSD(const std::string& name,
32  const edm::EventSetup& es,
33  const SensitiveDetectorCatalog& clg,
34  edm::ParameterSet const& p,
35  const SimTrackManager* manager)
36  : CaloSD(name,
37  es,
38  clg,
39  p,
40  manager,
41  static_cast<float>(p.getParameter<edm::ParameterSet>("HGCSD").getParameter<double>("TimeSliceUnit")),
42  p.getParameter<edm::ParameterSet>("HGCSD").getParameter<bool>("IgnoreTrackID")),
43  hgcons_(nullptr),
44  slopeMin_(0),
45  levelT1_(99),
46  levelT2_(99),
47  tan30deg_(std::tan(30.0 * CLHEP::deg)) {
48  numberingScheme_.reset(nullptr);
49  mouseBite_.reset(nullptr);
50 
51  edm::ParameterSet m_HGC = p.getParameter<edm::ParameterSet>("HGCSD");
52  eminHit_ = m_HGC.getParameter<double>("EminHit") * CLHEP::MeV;
53  fiducialCut_ = m_HGC.getParameter<bool>("FiducialCut");
54  distanceFromEdge_ = m_HGC.getParameter<double>("DistanceFromEdge");
55  storeAllG4Hits_ = m_HGC.getParameter<bool>("StoreAllG4Hits");
56  rejectMB_ = m_HGC.getParameter<bool>("RejectMouseBite");
57  waferRot_ = m_HGC.getParameter<bool>("RotatedWafer");
58  cornerMinMask_ = m_HGC.getParameter<int>("CornerMinMask");
59  angles_ = m_HGC.getUntrackedParameter<std::vector<double>>("WaferAngles");
60 
61  if (storeAllG4Hits_) {
62  setUseMap(false);
63  setNumberCheckedHits(0);
64  }
65 
66  //this is defined in the hgcsens.xml
67  G4String myName = name;
68  mydet_ = DetId::Forward;
69  nameX_ = "HGCal";
70  if (myName.find("HitsEE") != std::string::npos) {
71  mydet_ = DetId::HGCalEE;
72  nameX_ = "HGCalEESensitive";
73  } else if (myName.find("HitsHEfront") != std::string::npos) {
74  mydet_ = DetId::HGCalHSi;
75  nameX_ = "HGCalHESiliconSensitive";
76  }
77 
78 #ifdef EDM_ML_DEBUG
79  edm::LogVerbatim("HGCSim") << "**************************************************"
80  << "\n"
81  << "* *"
82  << "\n"
83  << "* Constructing a HGCalSD with name " << name << "\n"
84  << "* *"
85  << "\n"
86  << "**************************************************";
87 #endif
88  edm::LogVerbatim("HGCSim") << "HGCalSD:: Threshold for storing hits: " << eminHit_ << " for " << nameX_
89  << " detector " << mydet_;
90  edm::LogVerbatim("HGCSim") << "Flag for storing individual Geant4 Hits " << storeAllG4Hits_;
91  edm::LogVerbatim("HGCSim") << "Fiducial volume cut with cut from eta/phi "
92  << "boundary " << fiducialCut_ << " at " << distanceFromEdge_;
93  edm::LogVerbatim("HGCSim") << "Reject MosueBite Flag: " << rejectMB_ << " cuts along " << angles_.size()
94  << " axes: " << angles_[0] << ", " << angles_[1];
95 }
96 
97 double HGCalSD::getEnergyDeposit(const G4Step* aStep) {
98  double r = aStep->GetPreStepPoint()->GetPosition().perp();
99  double z = std::abs(aStep->GetPreStepPoint()->GetPosition().z());
100 #ifdef EDM_ML_DEBUG
101  G4int parCode = aStep->GetTrack()->GetDefinition()->GetPDGEncoding();
102  G4String parName = aStep->GetTrack()->GetDefinition()->GetParticleName();
103  G4LogicalVolume* lv = aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
104  edm::LogVerbatim("HGCSim") << "HGCalSD: Hit from standard path from " << lv->GetName() << " for Track "
105  << aStep->GetTrack()->GetTrackID() << " (" << parCode << ":" << parName << ") R = " << r
106  << " Z = " << z << " slope = " << r / z << ":" << slopeMin_;
107 #endif
108  // Apply fiducial cut
109  if (r < z * slopeMin_) {
110 #ifdef EDM_ML_DEBUG
111  edm::LogVerbatim("HGCSim") << "HGCalSD: Fiducial Volume cut";
112 #endif
113  return 0.0;
114  }
115 
116  double wt1 = getResponseWt(aStep->GetTrack());
117  double wt2 = aStep->GetTrack()->GetWeight();
118  double destep = weight_ * wt1 * (aStep->GetTotalEnergyDeposit());
119  if (wt2 > 0)
120  destep *= wt2;
121 #ifdef EDM_ML_DEBUG
122  edm::LogVerbatim("HGCSim") << "HGCalSD: weights= " << weight_ << ":" << wt1 << ":" << wt2 << " Total weight "
123  << weight_ * wt1 * wt2 << " deStep: " << aStep->GetTotalEnergyDeposit() << ":" << destep;
124 #endif
125  return destep;
126 }
127 
128 uint32_t HGCalSD::setDetUnitId(const G4Step* aStep) {
129  const G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
130  const G4VTouchable* touch = preStepPoint->GetTouchable();
131 
132 #ifdef EDM_ML_DEBUG
133  edm::LogVerbatim("HGCSim") << "DepthsTop: " << touch->GetHistoryDepth() << ":" << levelT1_ << ":" << levelT2_;
134  printDetectorLevels(touch);
135 #endif
136  //determine the exact position in global coordinates in the mass geometry
137  G4ThreeVector hitPoint = preStepPoint->GetPosition();
138  float globalZ = touch->GetTranslation(0).z();
139  int iz(globalZ > 0 ? 1 : -1);
140 
141  int layer(0), module(-1), cell(-1);
142  if (geom_mode_ == HGCalGeometryMode::Hexagon8Module) {
143  if (touch->GetHistoryDepth() > levelT2_) {
144  layer = touch->GetReplicaNumber(4);
145  module = touch->GetReplicaNumber(3);
146  cell = touch->GetReplicaNumber(1);
147  } else {
148  layer = touch->GetReplicaNumber(2);
149  module = touch->GetReplicaNumber(1);
150  }
151  } else if ((touch->GetHistoryDepth() == levelT1_) || (touch->GetHistoryDepth() == levelT2_)) {
152  layer = touch->GetReplicaNumber(0);
153 #ifdef EDM_ML_DEBUG
154  edm::LogVerbatim("HGCSim") << "DepthsTop: " << touch->GetHistoryDepth() << ":" << levelT1_ << ":" << levelT2_
155  << " name " << touch->GetVolume(0)->GetName() << " layer:module:cell " << layer << ":"
156  << module << ":" << cell;
157 #endif
158  } else {
159  layer = touch->GetReplicaNumber(3);
160  module = touch->GetReplicaNumber(2);
161  cell = touch->GetReplicaNumber(1);
162 #ifdef EDM_ML_DEBUG
163  edm::LogVerbatim("HGCSim") << "DepthsInside: " << touch->GetHistoryDepth() << " name "
164  << touch->GetVolume(0)->GetName() << " layer:module:cell " << layer << ":" << module
165  << ":" << cell;
166 #endif
167  }
168 #ifdef EDM_ML_DEBUG
169  G4Material* mat = aStep->GetPreStepPoint()->GetMaterial();
170  edm::LogVerbatim("HGCSim") << "Depths: " << touch->GetHistoryDepth() << " name " << touch->GetVolume(0)->GetName()
171  << ":" << touch->GetReplicaNumber(0) << " " << touch->GetVolume(1)->GetName() << ":"
172  << touch->GetReplicaNumber(1) << " " << touch->GetVolume(2)->GetName() << ":"
173  << touch->GetReplicaNumber(2) << " " << touch->GetVolume(3)->GetName() << ":"
174  << touch->GetReplicaNumber(3) << " " << touch->GetVolume(4)->GetName() << ":"
175  << touch->GetReplicaNumber(4) << " "
176  << " layer:module:cell " << layer << ":" << module << ":" << cell << " Material "
177  << mat->GetName() << ":" << mat->GetRadlen();
178 #endif
179  // The following statement should be examined later before elimination
180  if (aStep->GetPreStepPoint()->GetMaterial()->GetRadlen() > 100000.)
181  return 0;
182 
183  uint32_t id = setDetUnitId(layer, module, cell, iz, hitPoint);
184  if (rejectMB_ && id != 0) {
185  auto uv = HGCSiliconDetId(id).waferUV();
186 #ifdef EDM_ML_DEBUG
187  edm::LogVerbatim("HGCSim") << "ID " << std::hex << id << std::dec << " " << HGCSiliconDetId(id);
188 #endif
189  if (mouseBite_->exclude(hitPoint, iz, uv.first, uv.second)) {
190  id = 0;
191 #ifdef EDM_ML_DEBUG
192  edm::LogVerbatim("HGCSim") << "Rejected by mousebite cutoff *****";
193 #endif
194  }
195  }
196 #ifdef EDM_ML_DEBUG
197  if (id != 0)
198  edm::LogVerbatim("HGCSim") << HGCSiliconDetId(id);
199 #endif
200  return id;
201 }
202 
203 void HGCalSD::update(const BeginOfJob* job) {
204  const edm::EventSetup* es = (*job)();
205  edm::ESHandle<HGCalDDDConstants> hdc;
206  es->get<IdealGeometryRecord>().get(nameX_, hdc);
207  if (hdc.isValid()) {
208  hgcons_ = hdc.product();
209  geom_mode_ = hgcons_->geomMode();
210  slopeMin_ = hgcons_->minSlope();
211  levelT1_ = hgcons_->levelTop(0);
212  levelT2_ = hgcons_->levelTop(1);
213  double waferSize = hgcons_->waferSize(false);
214  double mouseBite = hgcons_->mouseBite(false);
215  mouseBiteCut_ = waferSize * tan30deg_ - mouseBite;
216 #ifdef EDM_ML_DEBUG
217  edm::LogVerbatim("HGCSim") << "HGCalSD::Initialized with mode " << geom_mode_ << " Slope cut " << slopeMin_
218  << " top Level " << levelT1_ << ":" << levelT2_ << " wafer " << waferSize << ":"
219  << mouseBite;
220 #endif
221 
222  numberingScheme_ = std::make_unique<HGCalNumberingScheme>(*hgcons_, mydet_, nameX_);
223  if (rejectMB_)
224  mouseBite_ = std::make_unique<HGCMouseBite>(*hgcons_, angles_, mouseBiteCut_, waferRot_);
225  } else {
226  throw cms::Exception("Unknown", "HGCalSD") << "Cannot find HGCalDDDConstants for " << nameX_ << "\n";
227  }
228 }
229 
230 void HGCalSD::initRun() {}
231 
232 bool HGCalSD::filterHit(CaloG4Hit* aHit, double time) {
233  return ((time <= tmaxHit) && (aHit->getEnergyDeposit() > eminHit_));
234 }
235 
236 uint32_t HGCalSD::setDetUnitId(int layer, int module, int cell, int iz, G4ThreeVector& pos) {
237  uint32_t id = numberingScheme_ ? numberingScheme_->getUnitID(layer, module, cell, iz, pos, weight_) : 0;
238  if (cornerMinMask_ > 2) {
239  if (hgcons_->maskCell(DetId(id), cornerMinMask_)) {
240  id = 0;
241  ignoreRejection();
242  }
243  }
244  if ((geom_mode_ == HGCalGeometryMode::Hexagon8File) || (geom_mode_ == HGCalGeometryMode::Hexagon8Module) || (id == 0))
245  ignoreRejection();
246  return id;
247 }
248 
249 bool HGCalSD::isItinFidVolume(const G4ThreeVector& pos) {
250  if (fiducialCut_) {
251  return (hgcons_->distFromEdgeHex(pos.x(), pos.y(), pos.z()) > distanceFromEdge_);
252  } else {
253  return true;
254  }
255 }
HGCalSD::storeAllG4Hits_
bool storeAllG4Hits_
Definition: HGCalSD.h:51
SimTrackManager
Definition: SimTrackManager.h:35
electrons_cff.bool
bool
Definition: electrons_cff.py:366
CaloSD::tmaxHit
double tmaxHit
Definition: CaloSD.h:141
dqmMemoryStats.float
float
Definition: dqmMemoryStats.py:127
ESHandle.h
HGCalSD::tan30deg_
const double tan30deg_
Definition: HGCalSD.h:53
CaloSD::ignoreRejection
void ignoreRejection()
Definition: CaloSD.h:105
HGCalSD::numberingScheme_
std::unique_ptr< HGCalNumberingScheme > numberingScheme_
Definition: HGCalSD.h:43
HGCalDDDConstants::geomMode
HGCalGeometryMode::GeometryMode geomMode() const
Definition: HGCalDDDConstants.h:54
edm
HLT enums.
Definition: AlignableModifier.h:19
HGCalGeometryMode.h
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
pos
Definition: PixelAliasList.h:18
HGCSiliconDetId.h
HGCSiliconDetId::waferUV
std::pair< int, int > waferUV() const
Definition: HGCSiliconDetId.h:78
HGCalSD::cornerMinMask_
int cornerMinMask_
Definition: HGCalSD.h:50
HGCalSD::geom_mode_
HGCalGeometryMode::GeometryMode geom_mode_
Definition: HGCalSD.h:47
HGCalSD::nameX_
std::string nameX_
Definition: HGCalSD.h:46
protons_cff.time
time
Definition: protons_cff.py:39
MeV
const double MeV
HGCalSD::fiducialCut_
bool fiducialCut_
Definition: HGCalSD.h:52
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
HGCalGeometryMode::Hexagon8Module
Definition: HGCalGeometryMode.h:34
HGCSiliconDetId
Definition: HGCSiliconDetId.h:22
HGCalSD.h
FastMath.h
CaloSD::getResponseWt
double getResponseWt(const G4Track *)
Definition: CaloSD.cc:876
DetId
Definition: DetId.h:17
DetId::HGCalHSi
Definition: DetId.h:33
DetId::HGCalEE
Definition: DetId.h:32
HGCalSD::HGCalSD
HGCalSD(const std::string &, const edm::EventSetup &, const SensitiveDetectorCatalog &, edm::ParameterSet const &, const SimTrackManager *)
Definition: HGCalSD.cc:31
edm::EventSetup::get
T get() const
Definition: EventSetup.h:87
CaloG4Hit::getEnergyDeposit
double getEnergyDeposit() const
Definition: CaloG4Hit.h:78
DDAxes::z
HGCalSD::distanceFromEdge_
double distanceFromEdge_
Definition: HGCalSD.h:48
edm::ESHandle
Definition: DTSurvey.h:22
SensitiveDetectorCatalog
Definition: SensitiveDetectorCatalog.h:10
HGCalSD::eminHit_
double eminHit_
Definition: HGCalSD.h:48
BeginOfJob
Definition: BeginOfJob.h:8
HGCalDDDConstants::mouseBite
double mouseBite(bool reco) const
Definition: HGCalDDDConstants.cc:880
phase1PixelTopology::layer
constexpr std::array< uint8_t, layerIndexSize > layer
Definition: phase1PixelTopology.h:99
CLHEP
Definition: CocoaGlobals.h:27
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
HGCalSD::setDetUnitId
uint32_t setDetUnitId(const G4Step *step) override
Definition: HGCalSD.cc:128
HGCalSD::mouseBite_
std::unique_ptr< HGCMouseBite > mouseBite_
Definition: HGCalSD.h:44
HGCalSD::update
void update(const BeginOfJob *) override
This routine will be called when the appropriate signal arrives.
Definition: HGCalSD.cc:203
edm::ParameterSet
Definition: ParameterSet.h:47
ParameterSet
Definition: Functions.h:16
HGCalSD::slopeMin_
double slopeMin_
Definition: HGCalSD.h:48
HGCalSD::levelT2_
int levelT2_
Definition: HGCalSD.h:50
funct::tan
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
CaloSD::printDetectorLevels
void printDetectorLevels(const G4VTouchable *) const
Definition: CaloSD.cc:1104
CaloSD::setNumberCheckedHits
void setNumberCheckedHits(int val)
Definition: CaloSD.h:119
IdealGeometryRecord.h
HGCalSD::getEnergyDeposit
double getEnergyDeposit(const G4Step *) override
Definition: HGCalSD.cc:97
edm::EventSetup
Definition: EventSetup.h:58
TrackInformation.h
CaloG4Hit
Definition: CaloG4Hit.h:32
get
#define get
HGCalSD::mouseBiteCut_
double mouseBiteCut_
Definition: HGCalSD.h:49
HGCalSD::isItinFidVolume
bool isItinFidVolume(const G4ThreeVector &)
Definition: HGCalSD.cc:249
alignCSCRings.r
r
Definition: alignCSCRings.py:93
HGCalSD::filterHit
bool filterHit(CaloG4Hit *, double) override
Definition: HGCalSD.cc:232
HGCalDDDConstants::minSlope
double minSlope() const
Definition: HGCalDDDConstants.h:96
std
Definition: JetResolutionObject.h:76
HGCalSD::mydet_
DetId::Detector mydet_
Definition: HGCalSD.h:45
HGCalSD::initRun
void initRun() override
Definition: HGCalSD.cc:230
edm::LogVerbatim
Log< level::Info, true > LogVerbatim
Definition: MessageLogger.h:128
triggerObjects_cff.id
id
Definition: triggerObjects_cff.py:29
HGCalSD::waferRot_
bool waferRot_
Definition: HGCalSD.h:52
Exception
Definition: hltDiff.cc:245
HGCalSD::hgcons_
const HGCalDDDConstants * hgcons_
Definition: HGCalSD.h:42
HGCalGeometryMode::Hexagon8File
Definition: HGCalGeometryMode.h:32
Skims_PA_cff.name
name
Definition: Skims_PA_cff.py:17
EventSetup.h
HGCalDDDConstants::distFromEdgeHex
double distFromEdgeHex(double x, double y, double z) const
Definition: HGCalDDDConstants.cc:279
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
Exception.h
HGCalSD::weight_
double weight_
Definition: HGCalSD.h:49
HGCalSD::rejectMB_
bool rejectMB_
Definition: HGCalSD.h:52
HGCalDDDConstants::levelTop
int levelTop(int ind=0) const
Definition: HGCalDDDConstants.h:88
CaloSD::setUseMap
void setUseMap(bool val)
Definition: CaloSD.h:108
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
HGCalDDDConstants.h
DetId::Forward
Definition: DetId.h:30
HGCalDDDConstants::waferSize
double waferSize(bool reco) const
Definition: HGCalDDDConstants.h:182
TauDecayModes.dec
dec
Definition: TauDecayModes.py:143
HGCalDDDConstants::maskCell
bool maskCell(const DetId &id, int corners) const
Definition: HGCalDDDConstants.cc:747
IdealGeometryRecord
Definition: IdealGeometryRecord.h:25
CaloSD
Definition: CaloSD.h:38
HGCalSD::angles_
std::vector< double > angles_
Definition: HGCalSD.h:54
HGCalSD::levelT1_
int levelT1_
Definition: HGCalSD.h:50
Generated for CMSSW Reference Manual by doxygen 1.8.16