CMS 3D CMS Logo

HGCSD.cc
Go to the documentation of this file.
1 // File: HGCSD.cc
3 // Description: Sensitive Detector class for Combined Forward Calorimeter
5 
7 
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 //#define plotDebug
31 
33  const HGCalDDDConstants* hgc,
34  const SensitiveDetectorCatalog& clg,
35  edm::ParameterSet const& p,
36  const SimTrackManager* manager)
37  : CaloSD(name,
38  clg,
39  p,
40  manager,
41  (float)(p.getParameter<edm::ParameterSet>("HGCSD").getParameter<double>("TimeSliceUnit")),
42  p.getParameter<edm::ParameterSet>("HGCSD").getParameter<bool>("IgnoreTrackID")),
43  hgcons_(hgc),
44  slopeMin_(0),
45  levelT_(99),
46  tree_(nullptr) {
47  numberingScheme_.reset(nullptr);
48  mouseBite_.reset(nullptr);
49 
50  edm::ParameterSet m_HGC = p.getParameter<edm::ParameterSet>("HGCSD");
51  eminHit_ = m_HGC.getParameter<double>("EminHit") * CLHEP::MeV;
52  storeAllG4Hits_ = m_HGC.getParameter<bool>("StoreAllG4Hits");
53  rejectMB_ = m_HGC.getParameter<bool>("RejectMouseBite");
54  waferRot_ = m_HGC.getParameter<bool>("RotatedWafer");
55  angles_ = m_HGC.getUntrackedParameter<std::vector<double>>("WaferAngles");
56  double waferSize = m_HGC.getUntrackedParameter<double>("WaferSize") * CLHEP::mm;
57  double mouseBite = m_HGC.getUntrackedParameter<double>("MouseBite") * CLHEP::mm;
58  mouseBiteCut_ = waferSize * tan(30.0 * CLHEP::deg) - mouseBite;
59 
60  if (storeAllG4Hits_) {
61  setUseMap(false);
63  }
64  //this is defined in the hgcsens.xml
65  G4String myName = name;
67  nameX_ = "HGCal";
68  if (myName.find("HitsEE") != std::string::npos) {
70  nameX_ = "HGCalEESensitive";
71  } else if (myName.find("HitsHEfront") != std::string::npos) {
73  nameX_ = "HGCalHESiliconSensitive";
74  } else if (myName.find("HitsHEback") != std::string::npos) {
76  nameX_ = "HGCalHEScintillatorSensitive";
77  }
78 
79 #ifdef EDM_ML_DEBUG
80  edm::LogVerbatim("HGCSim") << "**************************************************"
81  << "\n"
82  << "* *"
83  << "\n"
84  << "* Constructing a HGCSD with name " << name << "\n"
85  << "* *"
86  << "\n"
87  << "**************************************************";
88 #endif
89  edm::LogVerbatim("HGCSim") << "HGCSD:: Threshold for storing hits: " << eminHit << " for " << nameX_ << " subdet "
90  << myFwdSubdet_;
91  edm::LogVerbatim("HGCSim") << "Flag for storing individual Geant4 Hits " << storeAllG4Hits_;
92  edm::LogVerbatim("HGCSim") << "Reject MosueBite Flag: " << rejectMB_ << " Size of wafer " << waferSize
93  << " Mouse Bite " << mouseBite << ":" << mouseBiteCut_ << " along " << angles_.size()
94  << " axes";
95 }
96 
97 double HGCSD::getEnergyDeposit(const G4Step* aStep) {
98  double r = aStep->GetPreStepPoint()->GetPosition().perp();
99  double z = std::abs(aStep->GetPreStepPoint()->GetPosition().z());
100 
101 #ifdef EDM_ML_DEBUG
102  G4int parCode = aStep->GetTrack()->GetDefinition()->GetPDGEncoding();
103  G4LogicalVolume* lv = aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
104  edm::LogVerbatim("HGCSim") << "HGCSD: Hit from standard path from " << lv->GetName() << " for Track "
105  << aStep->GetTrack()->GetTrackID() << " ("
106  << aStep->GetTrack()->GetDefinition()->GetParticleName() << ":" << parCode << ") R = " << r
107  << " Z = " << z << " slope = " << r / z << ":" << slopeMin_;
108 #endif
109 
110  // Apply fiductial volume
111  if (r < z * slopeMin_) {
112  return 0.0;
113  }
114 
115  double wt1 = getResponseWt(aStep->GetTrack());
116  double wt2 = aStep->GetTrack()->GetWeight();
117  double destep = wt1 * aStep->GetTotalEnergyDeposit();
118  if (wt2 > 0)
119  destep *= wt2;
120 
121 #ifdef plotDebug
122  const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
123  G4double tmptrackE = aStep->GetTrack()->GetKineticEnergy();
124  G4int parCodex = aStep->GetTrack()->GetDefinition()->GetPDGEncoding();
125  G4double angle = (aStep->GetTrack()->GetMomentumDirection().theta()) / CLHEP::deg;
126  G4int layer = ((touch->GetHistoryDepth() == levelT_) ? touch->GetReplicaNumber(0) : touch->GetReplicaNumber(2));
127  G4int ilayer = (layer - 1) / 3;
128  if (aStep->GetTotalEnergyDeposit() > 0) {
129  t_Layer_.emplace_back(ilayer);
130  t_Parcode_.emplace_back(parCodex);
131  t_dEStep1_.emplace_back(aStep->GetTotalEnergyDeposit());
132  t_dEStep2_.emplace_back(destep);
133  t_TrackE_.emplace_back(tmptrackE);
134  t_Angle_.emplace_back(angle);
135  }
136 #endif
137 
138  return destep;
139 }
140 
141 uint32_t HGCSD::setDetUnitId(const G4Step* aStep) {
142  const G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
143  const G4VTouchable* touch = preStepPoint->GetTouchable();
144 
145  //determine the exact position in global coordinates in the mass geometry
146  G4ThreeVector hitPoint = preStepPoint->GetPosition();
147  float globalZ = touch->GetTranslation(0).z();
148  int iz(globalZ > 0 ? 1 : -1);
149 
150  //convert to local coordinates (=local to the current volume):
151  G4ThreeVector localpos = touch->GetHistory()->GetTopTransform().TransformPoint(hitPoint);
152 
153  //get the det unit id with
155 
156  int layer(-1), moduleLev(-1), module(-1), cell(-1);
157  if (touch->GetHistoryDepth() == levelT_) {
158  layer = touch->GetReplicaNumber(0);
159 #ifdef EDM_ML_DEBUG
160  edm::LogVerbatim("HGCSim") << "Depths: " << touch->GetHistoryDepth() << " name " << touch->GetVolume(0)->GetName()
161  << " layer:module:cell " << layer << ":" << moduleLev << ":" << module << ":" << cell;
162 #endif
163  } else {
164  layer = touch->GetReplicaNumber(2);
165  module = touch->GetReplicaNumber(1);
166  cell = touch->GetReplicaNumber(0);
167  moduleLev = 1;
168  }
169 #ifdef EDM_ML_DEBUG
170  const G4Material* mat = aStep->GetPreStepPoint()->GetMaterial();
171  edm::LogVerbatim("HGCSim") << "Depths: " << touch->GetHistoryDepth() << " name " << touch->GetVolume(0)->GetName()
172  << ":" << touch->GetReplicaNumber(0) << " " << touch->GetVolume(1)->GetName() << ":"
173  << touch->GetReplicaNumber(1) << " " << touch->GetVolume(2)->GetName() << ":"
174  << touch->GetReplicaNumber(2) << " layer:module:cell " << layer << ":" << moduleLev
175  << ":" << module << ":" << cell << " Material " << mat->GetName() << ":"
176  << mat->GetRadlen();
177  for (int k = 0; k < touch->GetHistoryDepth(); ++k)
178  edm::LogVerbatim("HGCSim") << "Level [" << k << "] " << touch->GetVolume(k)->GetName() << ":"
179  << touch->GetReplicaNumber(k);
180 #endif
181  // The following statement should be examined later before elimination
182  // VI: this is likely a check if media is vacuum - not needed
183  if (aStep->GetPreStepPoint()->GetMaterial()->GetRadlen() > 100000.)
184  return 0;
185 
186  uint32_t id = setDetUnitId(subdet, layer, module, cell, iz, localpos);
187  if (rejectMB_ && id != 0) {
188  int det, z, lay, wafer, type, ic;
189  HGCalTestNumbering::unpackHexagonIndex(id, det, z, lay, wafer, type, ic);
190 #ifdef EDM_ML_DEBUG
191  edm::LogVerbatim("HGCSim") << "ID " << std::hex << id << std::dec << " Input " << subdet << ":" << layer << ":"
192  << module << ":" << cell << ":" << iz << localpos.x() << ":" << localpos.y()
193  << " Decode " << det << ":" << z << ":" << lay << ":" << wafer << ":" << type << ":"
194  << ic;
195 #endif
196  G4ThreeVector local =
197  ((moduleLev >= 0) ? (touch->GetHistory()->GetTransform(moduleLev).TransformPoint(hitPoint)) : G4ThreeVector());
198  if (mouseBite_->exclude(local, z, layer, wafer, 0))
199  id = 0;
200  }
201  return id;
202 }
203 
204 void HGCSD::update(const BeginOfJob* job) {
205  if (hgcons_ != nullptr) {
208  levelT_ = hgcons_->levelTop();
209  numberingScheme_ = std::make_unique<HGCNumberingScheme>(*hgcons_, nameX_);
210  if (rejectMB_)
211  mouseBite_ = std::make_unique<HGCMouseBite>(*hgcons_, angles_, mouseBiteCut_, waferRot_);
212  } else {
213  edm::LogError("HGCSim") << "HGCSD : Cannot find HGCalDDDConstants for " << nameX_;
214  throw cms::Exception("Unknown", "HGCSD") << "Cannot find HGCalDDDConstants for " << nameX_ << "\n";
215  }
216 #ifdef EDM_ML_DEBUG
217  edm::LogVerbatim("HGCSim") << "HGCSD::Initialized with mode " << geom_mode_ << " Slope cut " << slopeMin_
218  << " top Level " << levelT_;
219 #endif
220 }
221 
223 #ifdef plotDebug
225  if (tfile.isAvailable()) {
226  tree_ = tfile->make<TTree>("TreeHGCSD", "TreeHGCSD");
227  tree_->Branch("EventID", &t_EventID_);
228  tree_->Branch("Layer", &t_Layer_);
229  tree_->Branch("ParticleCode", &t_Parcode_);
230  tree_->Branch("dEStepOriginal", &t_dEStep1_);
231  tree_->Branch("dEStepWeighted", &t_dEStep2_);
232  tree_->Branch("TrackEnergy", &t_TrackE_);
233  tree_->Branch("ThetaAngle", &t_Angle_);
234  }
235 #endif
236 }
237 
238 void HGCSD::initEvent(const BeginOfEvent* g4Event) {
239  const G4Event* evt = (*g4Event)();
240  t_EventID_ = evt->GetEventID();
241 #ifdef plotDebug
242  t_Layer_.clear();
243  t_Parcode_.clear();
244  t_dEStep1_.clear();
245  t_dEStep2_.clear();
246  t_TrackE_.clear();
247  t_Angle_.clear();
248 #endif
249 }
250 
252 #ifdef plotDebug
253  if (tree_)
254  tree_->Fill();
255 #endif
256 }
257 
258 bool HGCSD::filterHit(CaloG4Hit* aHit, double time) {
259  return ((time <= tmaxHit) && (aHit->getEnergyDeposit() > eminHit_));
260 }
261 
262 uint32_t HGCSD::setDetUnitId(ForwardSubdetector& subdet, int layer, int module, int cell, int iz, G4ThreeVector& pos) {
263  uint32_t id = numberingScheme_ ? numberingScheme_->getUnitID(subdet, layer, module, cell, iz, pos) : 0;
264  return id;
265 }
std::vector< double > t_TrackE_
Definition: HGCSD.h:62
Log< level::Info, true > LogVerbatim
double getEnergyDeposit(const G4Step *) override
Definition: HGCSD.cc:97
std::vector< double > t_dEStep1_
Definition: HGCSD.h:62
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
void setNumberCheckedHits(int val)
Definition: CaloSD.h:122
std::vector< double > t_dEStep2_
Definition: HGCSD.h:62
int levelT_
Definition: HGCSD.h:54
double mouseBiteCut_
Definition: HGCSD.h:56
Definition: CaloSD.h:40
std::vector< int > t_Layer_
Definition: HGCSD.h:61
std::string nameX_
Definition: HGCSD.h:47
std::vector< double > angles_
Definition: HGCSD.h:57
Log< level::Error, false > LogError
bool storeAllG4Hits_
Definition: HGCSD.h:55
HGCalGeometryMode::GeometryMode geomMode() const
ForwardSubdetector
bool waferRot_
Definition: HGCSD.h:55
void setUseMap(bool val)
Definition: CaloSD.h:111
double eminHit
Definition: CaloSD.h:144
T getUntrackedParameter(std::string const &, T const &) const
HGCalGeometryMode::GeometryMode geom_mode_
Definition: HGCSD.h:48
double slopeMin_
Definition: HGCSD.h:53
std::unique_ptr< HGCMouseBite > mouseBite_
Definition: HGCSD.h:50
Definition: tfile.py:1
const HGCalDDDConstants * hgcons_
Definition: HGCSD.h:46
void initRun() override
Definition: HGCSD.cc:222
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
std::vector< int > t_Parcode_
Definition: HGCSD.h:61
std::unique_ptr< HGCNumberingScheme > numberingScheme_
Definition: HGCSD.h:49
void initEvent(const BeginOfEvent *) override
Definition: HGCSD.cc:238
double tmaxHit
Definition: CaloSD.h:144
double eminHit_
Definition: HGCSD.h:51
std::vector< double > t_Angle_
Definition: HGCSD.h:63
bool rejectMB_
Definition: HGCSD.h:55
double minSlope() const
TTree * tree_
Definition: HGCSD.h:59
double getEnergyDeposit() const
Definition: CaloG4Hit.h:79
HGCSD(const std::string &, const HGCalDDDConstants *, const SensitiveDetectorCatalog &, edm::ParameterSet const &, const SimTrackManager *)
Definition: HGCSD.cc:32
HLT enums.
uint32_t setDetUnitId(const G4Step *step) override
Definition: HGCSD.cc:141
void update(const BeginOfJob *) override
This routine will be called when the appropriate signal arrives.
Definition: HGCSD.cc:204
bool filterHit(CaloG4Hit *, double) override
Definition: HGCSD.cc:258
double getResponseWt(const G4Track *)
Definition: CaloSD.cc:860
int levelTop(int ind=0) const
ForwardSubdetector myFwdSubdet_
Definition: HGCSD.h:52
void endEvent() override
Definition: HGCSD.cc:251
static void unpackHexagonIndex(const uint32_t &idx, int &subdet, int &z, int &lay, int &wafer, int &celltyp, int &cell)
uint32_t t_EventID_
Definition: HGCSD.h:60
T angle(T x1, T y1, T z1, T x2, T y2, T z2)
Definition: angle.h:11