CMS 3D CMS Logo

SimG4HGCalValidation.cc
Go to the documentation of this file.
1 // File: SimG4HGCalValidation.cc
3 // Description: Main analysis class for HGCal Validation of G4 Hits
5 
12 
13 // to retreive hits
21 
27 
31 
35 
36 #include "G4Step.hh"
37 #include "G4Track.hh"
38 #include "G4NavigationHistory.hh"
39 #include "G4TouchableHistory.hh"
40 
41 #include "CLHEP/Units/GlobalSystemOfUnits.h"
42 #include "CLHEP/Units/GlobalPhysicalConstants.h"
43 
44 #include <cmath>
45 #include <iostream>
46 #include <iomanip>
47 #include <memory>
48 #include <vector>
49 #include <string>
50 
52  public Observer<const BeginOfJob*>,
53  public Observer<const BeginOfEvent*>,
54  public Observer<const G4Step*> {
55 public:
57  ~SimG4HGCalValidation() override;
58 
59  void produce(edm::Event&, const edm::EventSetup&) override;
60 
61 private:
62  SimG4HGCalValidation(const SimG4HGCalValidation&); // stop default
64 
65  void init();
66 
67  // observer classes
68  void update(const BeginOfJob* job) override;
69  void update(const BeginOfEvent* evt) override;
70  void update(const G4Step* step) override;
71 
72  // analysis related class
74  void clear();
75 
76 private:
77  //HGCal numbering scheme
78  std::vector<HGCNumberingScheme*> hgcNumbering_;
79  std::vector<HGCalNumberingScheme*> hgcalNumbering_;
80 
81  // to read from parameter set
82  std::vector<std::string> names_;
83  std::vector<int> types_, detTypes_, subdet_;
85 
86  // parameters from geometry
88 
89  // some private members for ananlysis
90  unsigned int count_;
93  std::vector<double> hgcEEedep_, hgcHEFedep_, hgcHEBedep_;
94  std::vector<unsigned int> dets_, hgchitDets_, hgchitIndex_;
95  std::vector<double> hgchitX_, hgchitY_, hgchitZ_;
96 };
97 
98 SimG4HGCalValidation::SimG4HGCalValidation(const edm::ParameterSet& p) : levelT1_(999), levelT2_(999), count_(0) {
99  edm::ParameterSet m_Anal = p.getParameter<edm::ParameterSet>("SimG4HGCalValidation");
100  names_ = m_Anal.getParameter<std::vector<std::string> >("Names");
101  types_ = m_Anal.getParameter<std::vector<int> >("Types");
102  detTypes_ = m_Anal.getParameter<std::vector<int> >("DetTypes");
103  labelLayer_ = m_Anal.getParameter<std::string>("LabelLayerInfo");
104  verbosity_ = m_Anal.getUntrackedParameter<int>("Verbosity", 0);
105 
106  produces<PHGCalValidInfo>(labelLayer_);
107 
108  if (verbosity_ > 0) {
109  edm::LogVerbatim("ValidHGCal") << "HGCalTestAnalysis:: Initialised as "
110  << "observer of begin events and of G4step "
111  << "with Parameter values: \n\tLabel : " << labelLayer_ << " and with "
112  << names_.size() << " detectors";
113  for (unsigned int k = 0; k < names_.size(); ++k)
114  edm::LogVerbatim("ValidHGCal") << " [" << k << "] " << names_[k] << " Type " << types_[k] << " DetType "
115  << detTypes_[k];
116  }
117 }
118 
120  for (auto number : hgcNumbering_)
121  delete number;
122  for (auto number : hgcalNumbering_)
123  delete number;
124 }
125 
127  std::unique_ptr<PHGCalValidInfo> productLayer(new PHGCalValidInfo);
128  layerAnalysis(*productLayer);
129  e.put(std::move(productLayer), labelLayer_);
130 }
131 
133  const edm::EventSetup* es = (*job)();
134  for (unsigned int type = 0; type < types_.size(); ++type) {
135  int layers(0);
136  int detType = detTypes_[type];
137  G4String nameX = "HGCal";
138  if (types_[type] <= 1) {
139  if (types_[type] == 0) {
141  if (detType == 0)
142  dets_.emplace_back((int)(DetId::HGCalEE));
143  else if (detType == 1)
144  dets_.emplace_back((int)(DetId::HGCalHSi));
145  else
146  dets_.emplace_back((int)(DetId::HGCalHSc));
147  } else {
148  dets_.push_back((unsigned int)(DetId::Forward));
149  if (detType == 0)
150  subdet_.emplace_back((int)(ForwardSubdetector::HGCEE));
151  else if (detType == 1)
152  subdet_.emplace_back((int)(ForwardSubdetector::HGCHEF));
153  else
154  subdet_.emplace_back((int)(ForwardSubdetector::HGCHEB));
155  }
156  if (detType == 0)
157  nameX = "HGCalEESensitive";
158  else if (detType == 1)
159  nameX = "HGCalHESiliconSensitive";
160  else
161  nameX = "HGCalHEScintillatorSensitive";
163  es->get<IdealGeometryRecord>().get(nameX, hdc);
164  if (hdc.isValid()) {
165  levelT1_ = hdc->levelTop(0);
166  levelT2_ = hdc->levelTop(1);
167  if (hdc->tileTrapezoid()) {
168  types_[type] = -1;
169  hgcalNumbering_.emplace_back(new HGCalNumberingScheme(*hdc, (DetId::Detector)(dets_[type]), nameX));
170  } else if (hdc->waferHexagon6()) {
171  types_[type] = 1;
172  hgcNumbering_.push_back(new HGCNumberingScheme(*hdc, nameX));
173  } else {
174  types_[type] = 0;
175  hgcalNumbering_.emplace_back(new HGCalNumberingScheme(*hdc, (DetId::Detector)(dets_[type]), nameX));
176  }
177  layers = hdc->layers(false);
178  } else {
179  edm::LogError("ValidHGCal") << "Cannot find HGCalDDDConstants for " << nameX;
180  throw cms::Exception("Unknown", "ValidHGCal") << "Cannot find HGCalDDDConstants for " << nameX << "\n";
181  }
182  } else {
183  edm::LogError("ValidHGCal") << "Wrong Type " << types_[type];
184  throw cms::Exception("Unknown", "ValidHGCal") << "Wrong Type " << types_[type] << "\n";
185  }
186  if (detType == 0) {
187  for (int i = 0; i < layers; ++i)
188  hgcEEedep_.push_back(0);
189  } else if (detType == 1) {
190  for (int i = 0; i < layers; ++i)
191  hgcHEFedep_.push_back(0);
192  } else {
193  for (int i = 0; i < layers; ++i)
194  hgcHEBedep_.push_back(0);
195  }
196  if (verbosity_ > 0)
197  edm::LogVerbatim("ValidHGCal") << "[" << type << "]: " << nameX << " det " << dets_[type] << " subdet "
198  << subdet_[type] << " with " << layers << " layers";
199  }
200 }
201 
202 //=================================================================== per EVENT
204  int iev = (*evt)()->GetEventID();
205  if (verbosity_ > 0)
206  edm::LogVerbatim("ValidHGCal") << "SimG4HGCalValidation: =====> Begin "
207  << "event = " << iev;
208 
209  ++count_;
210  edepEE_ = edepHEF_ = edepHEB_ = 0.;
211 
212  //HGCal variables
213  for (unsigned int i = 0; i < hgcEEedep_.size(); i++)
214  hgcEEedep_[i] = 0.;
215  for (unsigned int i = 0; i < hgcHEFedep_.size(); i++)
216  hgcHEFedep_[i] = 0.;
217  for (unsigned int i = 0; i < hgcHEBedep_.size(); i++)
218  hgcHEBedep_[i] = 0.;
219 
220  //Cache reset
221  clear();
222 }
223 
224 //=================================================================== each STEP
225 void SimG4HGCalValidation::update(const G4Step* aStep) {
226  if (aStep != nullptr) {
227  G4VPhysicalVolume* curPV = aStep->GetPreStepPoint()->GetPhysicalVolume();
228  const G4String& name = curPV->GetName();
229  G4VSensitiveDetector* curSD = aStep->GetPreStepPoint()->GetSensitiveDetector();
230  if (verbosity_ > 1)
231  edm::LogVerbatim("ValidHGCal") << "ValidHGCal::Step in " << name << " at "
232  << aStep->GetPreStepPoint()->GetPosition();
233 
234  // Only for Sensitive detector
235  if (curSD != nullptr) {
236  int type(-1);
237  for (unsigned int k = 0; k < names_.size(); ++k) {
238  if (name.find(names_[k].c_str()) != std::string::npos) {
239  type = k;
240  break;
241  }
242  }
243  int detType = (type >= 0) ? detTypes_[type] : -1;
244  if (verbosity_ > 0)
245  edm::LogVerbatim("ValidHGCal") << "ValidHGCal::In SD " << name << " type " << type << ":" << detType;
246 
247  // Right type of SD
248  if (type >= 0) {
249  //Get the 32-bit index of the hit cell
250  const G4TouchableHistory* touchable =
251  static_cast<const G4TouchableHistory*>(aStep->GetPreStepPoint()->GetTouchable());
252  unsigned int index(0);
253  int layer(0);
254  G4ThreeVector hitPoint = aStep->GetPreStepPoint()->GetPosition();
255  // HGCal
256  G4ThreeVector localpos = touchable->GetHistory()->GetTopTransform().TransformPoint(hitPoint);
257  float globalZ = touchable->GetTranslation(0).z();
258  int iz(globalZ > 0 ? 1 : -1);
259  int module(-1), cell(-1);
260  if (types_[type] == 1) {
261  if (touchable->GetHistoryDepth() == levelT1_) {
262  layer = touchable->GetReplicaNumber(0);
263  } else {
264  layer = touchable->GetReplicaNumber(2);
265  module = touchable->GetReplicaNumber(1);
266  cell = touchable->GetReplicaNumber(0);
267  }
268  index =
269  hgcNumbering_[type]->getUnitID((ForwardSubdetector)(subdet_[type]), layer, module, cell, iz, localpos);
270  } else {
271  if ((touchable->GetHistoryDepth() == levelT1_) || (touchable->GetHistoryDepth() == levelT2_)) {
272  layer = touchable->GetReplicaNumber(0);
273  } else {
274  layer = touchable->GetReplicaNumber(3);
275  module = touchable->GetReplicaNumber(2);
276  cell = touchable->GetReplicaNumber(1);
277  }
278  double weight(0);
279  index = hgcalNumbering_[type]->getUnitID(layer, module, cell, iz, hitPoint, weight);
280  }
281  if (verbosity_ > 1)
282  edm::LogVerbatim("ValidHGCal") << "HGCal: " << name << " Layer " << layer << " Module " << module << " Cell "
283  << cell;
284 
285  double edeposit = aStep->GetTotalEnergyDeposit();
286  if (verbosity_ > 0)
287  edm::LogVerbatim("ValidHGCal") << "Layer " << layer << " Index " << std::hex << index << std::dec << " Edep "
288  << edeposit << " hit at " << hitPoint;
289  if (detType == 0) {
290  edepEE_ += edeposit;
291  if (layer < (int)(hgcEEedep_.size()))
292  hgcEEedep_[layer] += edeposit;
293  } else if (detType == 1) {
294  edepHEF_ += edeposit;
295  if (layer < (int)(hgcHEFedep_.size()))
296  hgcHEFedep_[layer] += edeposit;
297  } else {
298  edepHEB_ += edeposit;
299  if (layer < (int)(hgcHEBedep_.size()))
300  hgcHEBedep_[layer] += edeposit;
301  }
302  G4String nextVolume("XXX");
303  if (aStep->GetTrack()->GetNextVolume() != nullptr)
304  nextVolume = aStep->GetTrack()->GetNextVolume()->GetName();
305 
306  if (nextVolume.c_str() != name.c_str()) { //save hit when it exits cell
307  if (std::find(hgchitIndex_.begin(), hgchitIndex_.end(), index) == hgchitIndex_.end()) {
308  hgchitDets_.push_back(dets_[type]);
309  hgchitIndex_.push_back(index);
310  hgchitX_.push_back(hitPoint.x());
311  hgchitY_.push_back(hitPoint.y());
312  hgchitZ_.push_back(hitPoint.z());
313  }
314  }
315  } // it is right type of SD
316  } // it is in a SD
317  } //if aStep!=NULL
318 } //end update aStep
319 
320 //================================================================ End of EVENT
321 
323  if (verbosity_ > 0)
324  edm::LogVerbatim("ValidHGCal") << "\n ===>>> SimG4HGCalValidation: Energy "
325  << "deposit\n at EE : " << std::setw(6) << edepEE_ / CLHEP::MeV
326  << "\n at HEF: " << std::setw(6) << edepHEF_ / CLHEP::MeV
327  << "\n at HEB: " << std::setw(6) << edepHEB_ / CLHEP::MeV;
328 
329  //Fill HGC Variables
332 }
333 
334 //---------------------------------------------------
336  hgchitDets_.erase(hgchitDets_.begin(), hgchitDets_.end());
337  hgchitIndex_.erase(hgchitIndex_.begin(), hgchitIndex_.end());
338  hgchitX_.erase(hgchitX_.begin(), hgchitX_.end());
339  hgchitY_.erase(hgchitY_.begin(), hgchitY_.end());
340  hgchitZ_.erase(hgchitZ_.begin(), hgchitZ_.end());
341 }
342 
SimG4HGCalValidation::levelT2_
int levelT2_
Definition: SimG4HGCalValidation.cc:87
SimG4HGCalValidation::hgcalNumbering_
std::vector< HGCalNumberingScheme * > hgcalNumbering_
Definition: SimG4HGCalValidation.cc:79
mps_fire.i
i
Definition: mps_fire.py:428
PHGCalValidInfo
Definition: PHGCalValidInfo.h:12
Observer
Definition: Observer.h:23
MessageLogger.h
SimG4HGCalValidation::edepEE_
double edepEE_
Definition: SimG4HGCalValidation.cc:92
iev
const HitContainer *__restrict__ const TkSoA *__restrict__ const Quality *__restrict__ const CAHitNtupletGeneratorKernelsGPU::HitToTuple *__restrict__ int32_t int iev
Definition: CAHitNtupletGeneratorKernelsImpl.h:862
ForwardEmpty
Definition: ForwardSubdetector.h:5
ESHandle.h
BeginOfJob.h
step
step
Definition: StallMonitor.cc:94
HGCNumberingScheme.h
ForwardSubdetector
ForwardSubdetector
Definition: ForwardSubdetector.h:4
SimProducer.h
mps_merge.weight
weight
Definition: mps_merge.py:88
HGCalGeometryMode.h
MeV
const double MeV
EndOfEvent.h
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
SimG4HGCalValidation::dets_
std::vector< unsigned int > dets_
Definition: SimG4HGCalValidation.cc:94
SimG4HGCalValidation::hgchitIndex_
std::vector< unsigned int > hgchitIndex_
Definition: SimG4HGCalValidation.cc:94
spr::find
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:19
SimG4HGCalValidation::names_
std::vector< std::string > names_
Definition: SimG4HGCalValidation.cc:82
SimG4HGCalValidation::init
void init()
ForwardSubdetector.h
Observer.h
SimG4HGCalValidation
Definition: SimG4HGCalValidation.cc:51
SimG4HGCalValidation::hgcNumbering_
std::vector< HGCNumberingScheme * > hgcNumbering_
Definition: SimG4HGCalValidation.cc:78
SimG4HGCalValidation::subdet_
std::vector< int > subdet_
Definition: SimG4HGCalValidation.cc:83
DEFINE_SIMWATCHER
#define DEFINE_SIMWATCHER(type)
Definition: SimWatcherFactory.h:13
SimG4HGCalValidation::hgchitY_
std::vector< double > hgchitY_
Definition: SimG4HGCalValidation.cc:95
SimG4HGCalValidation::hgchitX_
std::vector< double > hgchitX_
Definition: SimG4HGCalValidation.cc:95
DetId::HGCalHSi
Definition: DetId.h:33
DetId::HGCalEE
Definition: DetId.h:32
MakerMacros.h
HGCalNumberingScheme.h
HGCalDDDConstants::waferHexagon6
bool waferHexagon6() const
Definition: HGCalDDDConstants.h:140
BeginOfRun.h
edm::EventSetup::get
T get() const
Definition: EventSetup.h:87
SimG4HGCalValidation::labelLayer_
std::string labelLayer_
Definition: SimG4HGCalValidation.cc:84
contentValuesFiles.number
number
Definition: contentValuesFiles.py:53
SimG4HGCalValidation::produce
void produce(edm::Event &, const edm::EventSetup &) override
Definition: SimG4HGCalValidation.cc:126
SimG4HGCalValidation::types_
std::vector< int > types_
Definition: SimG4HGCalValidation.cc:83
SimG4HGCalValidation::hgcHEFedep_
std::vector< double > hgcHEFedep_
Definition: SimG4HGCalValidation.cc:93
edm::ESHandle< HGCalDDDConstants >
SimG4HGCalValidation::operator=
const SimG4HGCalValidation & operator=(const SimG4HGCalValidation &)
SimG4HGCalValidation::hgcEEedep_
std::vector< double > hgcEEedep_
Definition: SimG4HGCalValidation.cc:93
BeginOfJob
Definition: BeginOfJob.h:8
dqmdumpme.k
k
Definition: dqmdumpme.py:60
SimG4HGCalValidation::edepHEB_
double edepHEB_
Definition: SimG4HGCalValidation.cc:92
SimG4HGCalValidation::hgchitZ_
std::vector< double > hgchitZ_
Definition: SimG4HGCalValidation.cc:95
phase1PixelTopology::layer
constexpr std::array< uint8_t, layerIndexSize > layer
Definition: phase1PixelTopology.h:99
HGCalDDDConstants::tileTrapezoid
bool tileTrapezoid() const
Definition: HGCalDDDConstants.h:113
HGCEE
Definition: ForwardSubdetector.h:8
SimG4HGCalValidation::~SimG4HGCalValidation
~SimG4HGCalValidation() override
Definition: SimG4HGCalValidation.cc:119
edm::ParameterSet
Definition: ParameterSet.h:47
AlCaHLTBitMon_ParallelJobs.p
def p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
SimG4HGCalValidation::hgchitDets_
std::vector< unsigned int > hgchitDets_
Definition: SimG4HGCalValidation.cc:94
Event.h
SimProducer
Definition: SimProducer.h:64
type
type
Definition: SiPixelVCal_PayloadInspector.cc:39
SimWatcherFactory.h
gainCalibHelper::gainCalibPI::type
type
Definition: SiPixelGainCalibHelper.h:40
BeginOfEvent.h
HGCalDDDConstants::layers
unsigned int layers(bool reco) const
Definition: HGCalDDDConstants.cc:577
ModuleDef.h
BeginOfEvent
Definition: BeginOfEvent.h:6
IdealGeometryRecord.h
DetId::Detector
Detector
Definition: DetId.h:24
edm::EventSetup
Definition: EventSetup.h:58
edm::LogError
Log< level::Error, false > LogError
Definition: MessageLogger.h:123
get
#define get
edm::ESHandleBase::isValid
bool isValid() const
Definition: ESHandle.h:44
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
HGCNumberingScheme
Definition: HGCNumberingScheme.h:13
PHGCalValidInfo::fillhgcLayers
void fillhgcLayers(const double edepEE, const double edepHEF, const double edepHEB, const std::vector< double > &eedep, const std::vector< double > &hefdep, const std::vector< double > &hebdep)
Definition: PHGCalValidInfo.cc:21
callgraph.module
module
Definition: callgraph.py:61
SimG4HGCalValidation::edepHEF_
double edepHEF_
Definition: SimG4HGCalValidation.cc:92
eostools.move
def move(src, dest)
Definition: eostools.py:511
SimG4HGCalValidation::SimG4HGCalValidation
SimG4HGCalValidation(const edm::ParameterSet &p)
Definition: SimG4HGCalValidation.cc:98
DetId.h
SimG4HGCalValidation::detTypes_
std::vector< int > detTypes_
Definition: SimG4HGCalValidation.cc:83
SimG4HGCalValidation::count_
unsigned int count_
Definition: SimG4HGCalValidation.cc:90
edm::LogVerbatim
Log< level::Info, true > LogVerbatim
Definition: MessageLogger.h:128
DetId::HGCalHSc
Definition: DetId.h:34
hfnoseParametersInitialization_cfi.nameX
nameX
Definition: hfnoseParametersInitialization_cfi.py:12
SimG4HGCalValidation::clear
void clear()
Definition: SimG4HGCalValidation.cc:335
Exception
Definition: hltDiff.cc:245
Point3D.h
Skims_PA_cff.name
name
Definition: Skims_PA_cff.py:17
EventSetup.h
PHGCalValidInfo::fillhgcHits
void fillhgcHits(const std::vector< unsigned int > &hitdets, const std::vector< unsigned int > &hitindex, const std::vector< double > &hitvtxX, const std::vector< double > &hitvtxY, const std::vector< double > &hitvtxZ)
Definition: PHGCalValidInfo.cc:7
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
AlignmentPI::index
index
Definition: AlignmentPayloadInspectorHelper.h:46
HGCalDDDConstants::levelTop
int levelTop(int ind=0) const
Definition: HGCalDDDConstants.h:89
HGCHEF
Definition: ForwardSubdetector.h:9
HGCalDDDConstants.h
ParameterSet.h
SimG4HGCalValidation::hgcHEBedep_
std::vector< double > hgcHEBedep_
Definition: SimG4HGCalValidation.cc:93
SimG4HGCalValidation::verbosity_
int verbosity_
Definition: SimG4HGCalValidation.cc:91
DetId::Forward
Definition: DetId.h:30
HGCalNumberingScheme
Definition: HGCalNumberingScheme.h:16
edm::Event
Definition: Event.h:73
SimG4HGCalValidation::update
void update(const BeginOfJob *job) override
This routine will be called when the appropriate signal arrives.
Definition: SimG4HGCalValidation.cc:132
PHGCalValidInfo.h
edm::Log
Definition: MessageLogger.h:70
TauDecayModes.dec
dec
Definition: TauDecayModes.py:142
hgcalTopologyTester_cfi.layers
layers
Definition: hgcalTopologyTester_cfi.py:8
SimG4HGCalValidation::levelT1_
int levelT1_
Definition: SimG4HGCalValidation.cc:87
weight
Definition: weight.py:1
IdealGeometryRecord
Definition: IdealGeometryRecord.h:25
SimG4HGCalValidation::layerAnalysis
void layerAnalysis(PHGCalValidInfo &)
Definition: SimG4HGCalValidation.cc:322
HGCHEB
Definition: ForwardSubdetector.h:10
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
HGCalTestNumbering.h