26 #include <CLHEP/Units/GlobalPhysicalConstants.h>
27 #include <CLHEP/Units/GlobalSystemOfUnits.h>
31 #include <TProfile2D.h>
63 const std::vector<edm::EDGetTokenT<edm::PCaloHitContainer>>
toks_calo_;
82 : g4Label_(ps.getUntrackedParameter<std::
string>(
"ModuleLabel",
"g4SimHits")),
83 hitLab_(ps.getParameter<std::
vector<std::
string>>(
"HitCollection")),
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)),
90 print_(ps.getUntrackedParameter<bool>(
"PrintExcessEnergy",
true)),
92 edm::
InputTag(ps.getUntrackedParameter<std::
string>(
"SourceLabel",
"VtxSmeared")))),
105 std::vector<std::string>
labels = {
"EcalHitsEB",
"HcalHits"};
116 descriptions.
add(
"HOSimHitStudy", desc);
124 <<
"please add it to config file";
126 char name[60],
title[100];
128 sprintf(title,
"Incident Energy (GeV)");
130 eneInc_->GetXaxis()->SetTitle(title);
131 eneInc_->GetYaxis()->SetTitle(
"Events");
132 sprintf(title,
"Incident #eta");
134 etaInc_->GetXaxis()->SetTitle(title);
135 etaInc_->GetYaxis()->SetTitle(
"Events");
136 sprintf(title,
"Incident #phi");
138 phiInc_->GetXaxis()->SetTitle(title);
139 phiInc_->GetYaxis()->SetTitle(
"Events");
140 int itcut = (int)(
tcut_);
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());
145 hit_[
i]->GetXaxis()->SetTitle(title);
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());
150 time_[
i]->GetXaxis()->SetTitle(title);
151 time_[
i]->GetYaxis()->SetTitle(
"Hits");
156 sprintf(name,
"Edep%d",
i);
157 sprintf(title,
"Energy deposit (GeV) in %s", dets[
i].c_str());
159 edep_[
i]->GetXaxis()->SetTitle(title);
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);
164 edepT_[
i]->GetXaxis()->SetTitle(title);
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());
169 hitTow_[
i]->GetXaxis()->SetTitle(title);
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());
178 edepTW_[
i]->GetXaxis()->SetTitle(title);
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);
183 edepTWT_[
i]->GetXaxis()->SetTitle(title);
184 edepTWT_[
i]->GetYaxis()->SetTitle(
"Towers");
185 sprintf(name,
"EdepZone%d",
i);
186 sprintf(title,
"Energy deposit (GeV) in %s", dets[
i].c_str());
188 edepZon_[
i]->GetXaxis()->SetTitle(title);
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)");
198 eEB_->GetXaxis()->SetTitle(title);
199 eEB_->GetYaxis()->SetTitle(
"Events");
200 sprintf(title,
"Energy Measured in EB+HB (GeV)");
202 eEBHB_->GetXaxis()->SetTitle(title);
203 eEBHB_->GetYaxis()->SetTitle(
"Events");
204 sprintf(title,
"Energy Measured in EB+HB+HO (GeV)");
206 eEBHBHO_->GetXaxis()->SetTitle(title);
207 eEBHBHO_->GetYaxis()->SetTitle(
"Events");
208 sprintf(title,
"Energy Measured in EB (GeV) for t < %d ns", itcut);
210 eEBT_->GetXaxis()->SetTitle(title);
211 eEBT_->GetYaxis()->SetTitle(
"Events");
212 sprintf(title,
"Energy Measured in EB+HB (GeV) for t < %d ns", itcut);
214 eEBHBT_->GetXaxis()->SetTitle(title);
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");
222 eHO1_->GetXaxis()->SetTitle(title);
223 eHO1_->GetYaxis()->SetTitle(
"Events");
224 eHO2_ = tfile->
make<TProfile2D>(
"EHO2",
title, 30, -1.305, 1.305, 72, -3.1415926, 3.1415926);
225 eHO2_->GetXaxis()->SetTitle(title);
226 eHO2_->GetYaxis()->SetTitle(
"Events");
227 sprintf(title,
"SimHit energy in HO Layer 17");
229 eHO17_->GetXaxis()->SetTitle(title);
230 eHO17_->GetYaxis()->SetTitle(
"Events");
231 sprintf(title,
"SimHit energy in HO Layer 18");
233 eHO18_->GetXaxis()->SetTitle(title);
234 eHO18_->GetYaxis()->SetTitle(
"Events");
235 sprintf(title,
"SimHit energy in HO for t < %d ns", itcut);
237 eHO1T_->GetXaxis()->SetTitle(title);
238 eHO1T_->GetYaxis()->SetTitle(
"Events");
239 eHO2T_ = tfile->
make<TProfile2D>(
"EHO2T",
title, 30, -1.305, 1.305, 72, -3.1415926, 3.1415926);
240 eHO2T_->GetXaxis()->SetTitle(title);
241 eHO2T_->GetYaxis()->SetTitle(
"Events");
242 sprintf(title,
"SimHit energy in HO Layer 17 for t < %d ns", itcut);
244 eHO17T_->GetXaxis()->SetTitle(title);
245 eHO17T_->GetYaxis()->SetTitle(
"Events");
246 sprintf(title,
"SimHit energy in HO Layer 18 for t < %d ns", itcut);
248 eHO18T_->GetXaxis()->SetTitle(title);
249 eHO18T_->GetYaxis()->SetTitle(
"Events");
250 sprintf(title,
"Number of layers hit in HO");
252 nHO1_->GetXaxis()->SetTitle(title);
253 nHO1_->GetYaxis()->SetTitle(
"Events");
254 nHO2_ = tfile->
make<TProfile2D>(
"NHO2",
title, 30, -1.305, 1.305, 72, -3.1415926, 3.1415926);
255 nHO2_->GetXaxis()->SetTitle(title);
256 nHO2_->GetYaxis()->SetTitle(
"Events");
257 sprintf(title,
"Number of layers hit in HO for t < %d ns", itcut);
259 nHO1T_->GetXaxis()->SetTitle(title);
260 nHO1T_->GetYaxis()->SetTitle(
"Events");
261 nHO2T_ = tfile->
make<TProfile2D>(
"NHO2T",
title, 30, -1.305, 1.305, 72, -3.1415926, 3.1415926);
262 nHO2T_->GetXaxis()->SetTitle(title);
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);
268 eHOE_[
i]->GetXaxis()->SetTitle(title);
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);
273 eHOE17_[
i]->GetXaxis()->SetTitle(title);
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);
278 eHOE18_[
i]->GetXaxis()->SetTitle(title);
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);
283 eHOET_[
i]->GetXaxis()->SetTitle(title);
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);
288 eHOE17T_[
i]->GetXaxis()->SetTitle(title);
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);
293 eHOE18T_[
i]->GetXaxis()->SetTitle(title);
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);
298 eHOEta_[
i]->GetXaxis()->SetTitle(title);
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);
313 eHOEtaT_[
i]->GetXaxis()->SetTitle(title);
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);
328 nHOE1_[
i]->GetXaxis()->SetTitle(title);
329 nHOE1_[
i]->GetYaxis()->SetTitle(
"Events");
330 sprintf(name,
"NHOE2%d",
i + 1);
332 nHOE2_[
i]->GetXaxis()->SetTitle(title);
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);
337 nHOE1T_[
i]->GetXaxis()->SetTitle(title);
338 nHOE1T_[
i]->GetYaxis()->SetTitle(
"Events");
339 sprintf(name,
"NHOE2%dT",
i + 1);
341 nHOE2T_[
i]->GetXaxis()->SetTitle(title);
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);
346 nHOEta1_[
i]->GetXaxis()->SetTitle(title);
347 nHOEta1_[
i]->GetYaxis()->SetTitle(
"Events");
348 sprintf(name,
"NHOEta2%d",
i + 1);
350 nHOEta2_[
i]->GetXaxis()->SetTitle(title);
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);
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();
461 edm::LogVerbatim(
"HitStudy") <<
"HOSimHitStudy:: Hit " <<
k <<
" Subdet:" << subdet <<
" zside:" << zside
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;
490 time_[indx]->Fill(time);
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) {
610 if (ieta >= 0 && ieta < 15) {
611 eHOE_[ieta]->Fill(eHO17 + eHO18);
614 eHOET_[ieta]->Fill(eHO17T + eHO18T);
617 eHOEta_[ieta]->Fill(eHOE17[ieta] + eHOE18[ieta]);
620 nHOE1_[ieta]->Fill((
double)(nHO));
622 nHOE1T_[ieta]->Fill((
double)(nHOT));
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)
633 nHOEta1_[ieta]->Fill((
double)(nHOE));
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++) {
665 int ieta, iphi,
depth = 0;
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
EventNumber_t event() const
const edm::EventSetup & c
void beginRun(edm::Run const &, edm::EventSetup const &) override
ParameterDescriptionBase * addUntracked(U const &iLabel, T const &value)
void analyze(const edm::Event &e, const edm::EventSetup &c) override
std::ostream & print_(std::ostream &os, value_type const &hash)
#define DEFINE_FWK_MODULE(type)
T * make(const Args &...args) const
make new ROOT object
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_
Handle< PROD > getHandle(EDGetTokenT< PROD > token) const
Abs< T >::type abs(const T &t)
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
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_