26 #include <CLHEP/Units/PhysicalConstants.h> 27 #include <CLHEP/Units/SystemOfUnits.h> 31 #include <TProfile2D.h> 63 const std::vector<edm::EDGetTokenT<edm::PCaloHitContainer>>
toks_calo_;
82 : g4Label_(ps.getUntrackedParameter<
std::
string>(
"ModuleLabel",
"g4SimHits")),
84 maxEnergy_(ps.getUntrackedParameter<double>(
"MaxEnergy", 200.0)),
85 scaleEB_(ps.getUntrackedParameter<double>(
"ScaleEB", 1.0)),
86 scaleHB_(ps.getUntrackedParameter<double>(
"ScaleHB", 100.0)),
87 scaleHO_(ps.getUntrackedParameter<double>(
"ScaleHO", 2.0)),
88 tcut_(ps.getUntrackedParameter<double>(
"TimeCut", 100.0)),
89 scheme_(ps.getUntrackedParameter<
bool>(
"TestNumbering",
false)),
98 edm::LogVerbatim(
"HitStudy") <<
"HOSimHitStudy::Module Label: " << g4Label_ <<
" Hits: " << hitLab_[0] <<
", " 99 << hitLab_[1] <<
" MaxEnergy: " << maxEnergy_ <<
" Scale factor for EB " << scaleEB_
100 <<
", for HB " << scaleHB_ <<
" and for HO " << scaleHO_ <<
" time Cut " << tcut_;
105 std::vector<std::string>
labels = {
"EcalHitsEB",
"HcalHits"};
108 desc.addUntracked<std::vector<std::string>>(
"HitCollection",
labels);
109 desc.addUntracked<
double>(
"MaxEnergy", 50.0);
110 desc.addUntracked<
double>(
"ScaleEB", 1.02);
111 desc.addUntracked<
double>(
"ScaleHB", 104.4);
112 desc.addUntracked<
double>(
"ScaleHO", 2.33);
113 desc.addUntracked<
double>(
"TimeCut", 100.0);
114 desc.addUntracked<
bool>(
"PrintExcessEnergy",
true);
115 desc.addUntracked<
bool>(
"TestNumbering",
false);
116 descriptions.
add(
"HOSimHitStudy",
desc);
122 if (!
tfile.isAvailable())
124 <<
"please add it to config file";
128 sprintf(
title,
"Incident Energy (GeV)");
131 eneInc_->GetYaxis()->SetTitle(
"Events");
132 sprintf(
title,
"Incident #eta");
135 etaInc_->GetYaxis()->SetTitle(
"Events");
136 sprintf(
title,
"Incident #phi");
139 phiInc_->GetYaxis()->SetTitle(
"Events");
141 for (
int i = 0;
i < 3;
i++) {
142 sprintf(
name,
"Hit%d",
i);
143 sprintf(
title,
"Number of hits in %s", dets[
i].c_str());
146 hit_[
i]->GetYaxis()->SetTitle(
"Events");
147 sprintf(
name,
"Time%d",
i);
148 sprintf(
title,
"Time of the hit (ns) in %s", dets[
i].c_str());
151 time_[
i]->GetYaxis()->SetTitle(
"Hits");
156 sprintf(
name,
"Edep%d",
i);
157 sprintf(
title,
"Energy deposit (GeV) in %s", dets[
i].c_str());
160 edep_[
i]->GetYaxis()->SetTitle(
"Hits");
161 sprintf(
name,
"EdepT%d",
i);
162 sprintf(
title,
"Energy deposit (GeV) in %s for t < %d ns", dets[
i].c_str(), itcut);
165 edepT_[
i]->GetYaxis()->SetTitle(
"Hits");
166 sprintf(
name,
"HitTow%d",
i);
167 sprintf(
title,
"Number of towers with hits in %s", dets[
i].c_str());
170 hitTow_[
i]->GetYaxis()->SetTitle(
"Events");
175 sprintf(
name,
"EdepTW%d",
i);
176 sprintf(
title,
"Energy deposit (GeV) in %s Tower", dets[
i].c_str());
179 edepTW_[
i]->GetYaxis()->SetTitle(
"Towers");
180 sprintf(
name,
"EdepTWT%d",
i);
181 sprintf(
title,
"Energy deposit (GeV) in %s Tower for t < %d ns", dets[
i].c_str(), itcut);
184 edepTWT_[
i]->GetYaxis()->SetTitle(
"Towers");
185 sprintf(
name,
"EdepZone%d",
i);
186 sprintf(
title,
"Energy deposit (GeV) in %s", dets[
i].c_str());
189 edepZon_[
i]->GetYaxis()->SetTitle(
"Events");
190 sprintf(
name,
"EdepZoneT%d",
i);
191 sprintf(
title,
"Energy deposit (GeV) in %s for t < %d ns", dets[
i].c_str(), itcut);
196 sprintf(
title,
"Energy Measured in EB (GeV)");
199 eEB_->GetYaxis()->SetTitle(
"Events");
200 sprintf(
title,
"Energy Measured in EB+HB (GeV)");
203 eEBHB_->GetYaxis()->SetTitle(
"Events");
204 sprintf(
title,
"Energy Measured in EB+HB+HO (GeV)");
207 eEBHBHO_->GetYaxis()->SetTitle(
"Events");
208 sprintf(
title,
"Energy Measured in EB (GeV) for t < %d ns", itcut);
211 eEBT_->GetYaxis()->SetTitle(
"Events");
212 sprintf(
title,
"Energy Measured in EB+HB (GeV) for t < %d ns", itcut);
215 eEBHBT_->GetYaxis()->SetTitle(
"Events");
216 sprintf(
title,
"Energy Measured in EB+HB+HO (GeV) for t < %d ns", itcut);
219 eEBHBHOT_->GetYaxis()->SetTitle(
"Events");
220 sprintf(
title,
"SimHit energy in HO");
223 eHO1_->GetYaxis()->SetTitle(
"Events");
224 eHO2_ =
tfile->make<TProfile2D>(
"EHO2",
title, 30, -1.305, 1.305, 72, -3.1415926, 3.1415926);
226 eHO2_->GetYaxis()->SetTitle(
"Events");
227 sprintf(
title,
"SimHit energy in HO Layer 17");
230 eHO17_->GetYaxis()->SetTitle(
"Events");
231 sprintf(
title,
"SimHit energy in HO Layer 18");
234 eHO18_->GetYaxis()->SetTitle(
"Events");
235 sprintf(
title,
"SimHit energy in HO for t < %d ns", itcut);
238 eHO1T_->GetYaxis()->SetTitle(
"Events");
239 eHO2T_ =
tfile->make<TProfile2D>(
"EHO2T",
title, 30, -1.305, 1.305, 72, -3.1415926, 3.1415926);
241 eHO2T_->GetYaxis()->SetTitle(
"Events");
242 sprintf(
title,
"SimHit energy in HO Layer 17 for t < %d ns", itcut);
245 eHO17T_->GetYaxis()->SetTitle(
"Events");
246 sprintf(
title,
"SimHit energy in HO Layer 18 for t < %d ns", itcut);
249 eHO18T_->GetYaxis()->SetTitle(
"Events");
250 sprintf(
title,
"Number of layers hit in HO");
253 nHO1_->GetYaxis()->SetTitle(
"Events");
254 nHO2_ =
tfile->make<TProfile2D>(
"NHO2",
title, 30, -1.305, 1.305, 72, -3.1415926, 3.1415926);
256 nHO2_->GetYaxis()->SetTitle(
"Events");
257 sprintf(
title,
"Number of layers hit in HO for t < %d ns", itcut);
260 nHO1T_->GetYaxis()->SetTitle(
"Events");
261 nHO2T_ =
tfile->make<TProfile2D>(
"NHO2T",
title, 30, -1.305, 1.305, 72, -3.1415926, 3.1415926);
263 nHO2T_->GetYaxis()->SetTitle(
"Events");
264 for (
int i = 0;
i < 15;
i++) {
265 sprintf(
name,
"EHOE%d",
i + 1);
266 sprintf(
title,
"SimHit energy in HO (Beam in #eta=%d bin)",
i + 1);
269 eHOE_[
i]->GetYaxis()->SetTitle(
"Events");
270 sprintf(
name,
"EHOE17%d",
i + 1);
271 sprintf(
title,
"SimHit energy in Layer 17 (Beam in #eta=%d bin)",
i + 1);
274 eHOE17_[
i]->GetYaxis()->SetTitle(
"Events");
275 sprintf(
name,
"EHOE18%d",
i + 1);
276 sprintf(
title,
"SimHit energy in Layer 18 (Beam in #eta=%d bin)",
i + 1);
279 eHOE18_[
i]->GetYaxis()->SetTitle(
"Events");
280 sprintf(
name,
"EHOE%dT",
i + 1);
281 sprintf(
title,
"SimHit energy in HO (Beam in #eta=%d bin, t < %d ns)",
i + 1, itcut);
284 eHOET_[
i]->GetYaxis()->SetTitle(
"Events");
285 sprintf(
name,
"EHOE17%dT",
i + 1);
286 sprintf(
title,
"SimHit energy in Layer 17 (Beam in #eta=%d bin, t < %d ns)",
i + 1, itcut);
289 eHOE17T_[
i]->GetYaxis()->SetTitle(
"Events");
290 sprintf(
name,
"EHOE18%dT",
i + 1);
291 sprintf(
title,
"SimHit energy in Layer 18 (Beam in #eta=%d bin, t < %d ns)",
i + 1, itcut);
294 eHOE18T_[
i]->GetYaxis()->SetTitle(
"Events");
295 sprintf(
name,
"EHOEta%d",
i + 1);
296 sprintf(
title,
"SimHit energy in HO #eta bin %d (Beam in #eta=%d bin)",
i + 1,
i + 1);
299 eHOEta_[
i]->GetYaxis()->SetTitle(
"Events");
300 sprintf(
name,
"EHOEta17%d",
i + 1);
301 sprintf(
title,
"SimHit energy in Layer 17 #eta bin %d (Beam in #eta=%d bin)",
i + 1,
i + 1);
305 sprintf(
name,
"EHOEta18%d",
i + 1);
306 sprintf(
title,
"SimHit energy in Layer 18 #eta bin %d (Beam in #eta=%d bin)",
i + 1,
i + 1);
310 sprintf(
name,
"EHOEta%dT",
i + 1);
311 sprintf(
title,
"SimHit energy in HO #eta bin %d (Beam in #eta=%d bin, t < %d ns)",
i + 1,
i + 1, itcut);
314 eHOEtaT_[
i]->GetYaxis()->SetTitle(
"Events");
315 sprintf(
name,
"EHOEta17%dT",
i + 1);
316 sprintf(
title,
"SimHit energy in Layer 17 #eta bin %d (Beam in #eta=%d bin, t < %d ns)",
i + 1,
i + 1, itcut);
320 sprintf(
name,
"EHOEta18%dT",
i + 1);
321 sprintf(
title,
"SimHit energy in Layer 18 #eta bin %d (Beam in #eta=%d bin, t < %d ns)",
i + 1,
i + 1, itcut);
325 sprintf(
name,
"NHOE1%d",
i + 1);
326 sprintf(
title,
"Number of layers hit in HO (Beam in #eta=%d bin)",
i + 1);
329 nHOE1_[
i]->GetYaxis()->SetTitle(
"Events");
330 sprintf(
name,
"NHOE2%d",
i + 1);
333 nHOE2_[
i]->GetYaxis()->SetTitle(
"Events");
334 sprintf(
name,
"NHOE1%dT",
i + 1);
335 sprintf(
title,
"Number of layers hit in HO (Beam in #eta=%d bin, t < %d ns)",
i + 1, itcut);
338 nHOE1T_[
i]->GetYaxis()->SetTitle(
"Events");
339 sprintf(
name,
"NHOE2%dT",
i + 1);
342 nHOE2T_[
i]->GetYaxis()->SetTitle(
"Events");
343 sprintf(
name,
"NHOEta1%d",
i + 1);
344 sprintf(
title,
"Number of layers hit in HO #eta bin %d (Beam in #eta=%d bin)",
i + 1,
i + 1);
347 nHOEta1_[
i]->GetYaxis()->SetTitle(
"Events");
348 sprintf(
name,
"NHOEta2%d",
i + 1);
351 nHOEta2_[
i]->GetYaxis()->SetTitle(
"Events");
352 sprintf(
name,
"NHOEta1%dT",
i + 1);
353 sprintf(
title,
"Number of layers hit in HO #eta bin %d (Beam in #eta=%d bin, t < %d ns)",
i + 1,
i + 1, itcut);
357 sprintf(
name,
"NHOEta2%dT",
i + 1);
365 edm::LogVerbatim(
"HitStudy") <<
"HOSimHitStudy::Run = " <<
e.id().run() <<
" Event = " <<
e.id().event();
371 HepMC::GenEvent::particle_const_iterator
p = myGenEvent->particles_begin();
372 if (
p != myGenEvent->particles_end()) {
373 eInc = (*p)->momentum().e();
374 etaInc = (*p)->momentum().eta();
375 phiInc = (*p)->momentum().phi();
379 <<
" Phi = " <<
phiInc / CLHEP::deg;
381 for (
int i = 0;
i < 2;
i++) {
382 bool getHits =
false;
410 double etot[3], etotT[3];
411 std::vector<unsigned int> ebID, hbID, hoID;
412 std::vector<double> ebEtow, hbEtow, hoEtow;
413 std::vector<double> ebEtowT, hbEtowT, hoEtowT;
414 for (
int k = 0;
k < 3;
k++) {
423 double eHO17 = 0, eHO18 = 0, eHO17T = 0, eHO18T = 0;
424 double eHOE17[15], eHOE18[15], eHOE17T[15], eHOE18T[15];
425 for (
int k = 0;
k < 15;
k++) {
426 eHOE17[
k] = eHOE18[
k] = eHOE17T[
k] = eHOE18T[
k] = 0;
429 for (
int k = 0;
k < 2;
k++) {
436 for (
int i = 0;
i < nHit;
i++) {
454 subdet =
id.subdet();
462 <<
" depth:" <<
depth <<
" layer:" << lay <<
" eta:" <<
eta <<
" phi:" <<
phi;
465 else if (subdet == static_cast<int>(
HcalOuter)) {
471 if (
eta >= 0 &&
eta < 15) {
474 eHOE17T[
eta] += edep;
480 if (
eta >= 0 &&
eta < 15) {
483 eHOE18T[
eta] += edep;
491 edep_[indx]->Fill(edep);
501 for (
unsigned int j = 0;
j < ebID.size();
j++) {
502 if (id_ == ebID[
j]) {
511 ebEtow.push_back(edep);
512 ebEtowT.push_back(edepT);
514 }
else if (indx == 1) {
516 for (
unsigned int j = 0;
j < hbID.size();
j++) {
517 if (id_ == hbID[
j]) {
526 hbEtow.push_back(edep);
527 hbEtowT.push_back(edepT);
531 for (
unsigned int j = 0;
j < hoID.size();
j++) {
532 if (id_ == hoID[
j]) {
541 hoEtow.push_back(edep);
542 hoEtowT.push_back(edepT);
550 for (
int k = 0;
k < 3;
k++) {
551 hit_[
k]->Fill(
double(nhit[
k]));
555 hitTow_[0]->Fill(
double(ebEtow.size()));
556 for (
unsigned int i = 0;
i < ebEtow.size();
i++) {
560 hitTow_[1]->Fill(
double(hbEtow.size()));
561 for (
unsigned int i = 0;
i < hbEtow.size();
i++) {
565 hitTow_[2]->Fill(
double(hoEtow.size()));
566 for (
unsigned int i = 0;
i < hoEtow.size();
i++) {
571 double eEBHB = eEB +
scaleHB_ * etot[1];
577 double eEBHBT = eEBT +
scaleHB_ * etotT[1];
590 int nHO = 0, nHOT = 0;
604 for (
int k = 0;
k < 15; ++
k) {
624 int nHOE = 0, nHOET = 0;
625 if (eHOE17[
ieta] > 0)
627 if (eHOE17T[
ieta] > 0)
629 if (eHOE18[
ieta] > 0)
631 if (eHOE18T[
ieta] > 0)
639 edm::LogVerbatim(
"HitStudy") <<
"HOSimHitStudy::analyzeHits: Hits in EB " << nhit[0] <<
" in " << ebEtow.size()
640 <<
" towers with total E " << etot[0] <<
"|" << etotT[0]
641 <<
"\n Hits in HB " << nhit[1] <<
" in " << hbEtow.size()
642 <<
" towers with total E " << etot[1] <<
"|" << etotT[1]
643 <<
"\n Hits in HO " << nhit[2] <<
" in " << hoEtow.size()
644 <<
" towers with total E " << etot[2] <<
"|" << etotT[2]
645 <<
"\n Energy in EB " << eEB <<
"|" << eEBT <<
" with HB " 646 << eEBHB <<
"|" << eEBHBT <<
" and with HO " << eEBHBHO <<
"|" << eEBHBHOT
647 <<
"\n E in HO layers " << eHO17 <<
"|" << eHO17T <<
" " 648 << eHO18 <<
"|" << eHO18T <<
" number of HO hits " << nHO <<
"|" << nHOT;
651 edm::LogVerbatim(
"HitStudy") <<
"HOSimHitStudy::Event with excessive energy: EB = " << eEB <<
" EB+HB = " << eEBHB
652 <<
" EB+HB+HO = " << eEBHBHO;
654 for (
int k = 0;
k < 2;
k++) {
661 for (
int i = 0;
i < nHit;
i++) {
680 int subdet =
id.subdet();
683 else if (subdet == static_cast<int>(
HcalOuter))
691 <<
" (" <<
ieta <<
"|" <<
iphi <<
"|" <<
depth <<
") " << std::setw(8) << edep
692 <<
" " << std::setw(8) <<
time;
static const std::string kSharedResource
Log< level::Info, true > LogVerbatim
void beginRun(edm::Run const &, edm::EventSetup const &) override
void analyze(const edm::Event &e, const edm::EventSetup &c) override
std::ostream & print_(std::ostream &os, value_type const &hash)
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
void endRun(edm::Run const &, edm::EventSetup const &) override
auto vector_transform(std::vector< InputType > const &input, Function predicate) -> std::vector< typename std::remove_cv< typename std::remove_reference< decltype(predicate(input.front()))>::type >::type >
const std::vector< std::string > hitLab_
Abs< T >::type abs(const T &t)
#define DEFINE_FWK_MODULE(type)
const HepMC::GenEvent * GetEvent() const
const std::vector< edm::EDGetTokenT< edm::PCaloHitContainer > > toks_calo_
~HOSimHitStudy() override=default
void add(std::string const &label, ParameterSetDescription const &psetDescription)
std::vector< PCaloHit > hcalHits
std::vector< PCaloHit > ecalHits
static void unpackHcalIndex(const uint32_t &idx, int &det, int &z, int &depth, int &eta, int &phi, int &lay)
const std::string g4Label_
HOSimHitStudy(const edm::ParameterSet &ps)
const edm::EDGetTokenT< edm::HepMCProduct > tok_evt_