26 #include "CLHEP/Units/GlobalPhysicalConstants.h"
27 #include "CLHEP/Units/GlobalSystemOfUnits.h"
31 #include <TProfile2D.h>
89 for (
unsigned i = 0;
i != 2;
i++)
110 desc.addUntracked<
double>(
"MaxEnergy", 50.0);
111 desc.addUntracked<
double>(
"ScaleEB", 1.02);
112 desc.addUntracked<
double>(
"ScaleHB", 104.4);
113 desc.addUntracked<
double>(
"ScaleHO", 2.33);
114 desc.addUntracked<
double>(
"TimeCut", 100.0);
115 desc.addUntracked<
bool>(
"PrintExcessEnergy",
true);
116 desc.addUntracked<
bool>(
"TestNumbering",
false);
117 descriptions.
add(
"HOSimHitStudy",
desc);
123 if (!
tfile.isAvailable())
125 <<
"please add it to config file";
129 sprintf(
title,
"Incident Energy (GeV)");
132 eneInc_->GetYaxis()->SetTitle(
"Events");
133 sprintf(
title,
"Incident #eta");
136 etaInc_->GetYaxis()->SetTitle(
"Events");
137 sprintf(
title,
"Incident #phi");
140 phiInc_->GetYaxis()->SetTitle(
"Events");
142 for (
int i = 0;
i < 3;
i++) {
143 sprintf(
name,
"Hit%d",
i);
144 sprintf(
title,
"Number of hits in %s", dets[
i].c_str());
147 hit_[
i]->GetYaxis()->SetTitle(
"Events");
148 sprintf(
name,
"Time%d",
i);
149 sprintf(
title,
"Time of the hit (ns) in %s", dets[
i].c_str());
152 time_[
i]->GetYaxis()->SetTitle(
"Hits");
157 sprintf(
name,
"Edep%d",
i);
158 sprintf(
title,
"Energy deposit (GeV) in %s", dets[
i].c_str());
161 edep_[
i]->GetYaxis()->SetTitle(
"Hits");
162 sprintf(
name,
"EdepT%d",
i);
163 sprintf(
title,
"Energy deposit (GeV) in %s for t < %d ns", dets[
i].c_str(), itcut);
166 edepT_[
i]->GetYaxis()->SetTitle(
"Hits");
167 sprintf(
name,
"HitTow%d",
i);
168 sprintf(
title,
"Number of towers with hits in %s", dets[
i].c_str());
171 hitTow_[
i]->GetYaxis()->SetTitle(
"Events");
176 sprintf(
name,
"EdepTW%d",
i);
177 sprintf(
title,
"Energy deposit (GeV) in %s Tower", dets[
i].c_str());
180 edepTW_[
i]->GetYaxis()->SetTitle(
"Towers");
181 sprintf(
name,
"EdepTWT%d",
i);
182 sprintf(
title,
"Energy deposit (GeV) in %s Tower for t < %d ns", dets[
i].c_str(), itcut);
185 edepTWT_[
i]->GetYaxis()->SetTitle(
"Towers");
186 sprintf(
name,
"EdepZone%d",
i);
187 sprintf(
title,
"Energy deposit (GeV) in %s", dets[
i].c_str());
190 edepZon_[
i]->GetYaxis()->SetTitle(
"Events");
191 sprintf(
name,
"EdepZoneT%d",
i);
192 sprintf(
title,
"Energy deposit (GeV) in %s for t < %d ns", dets[
i].c_str(), itcut);
197 sprintf(
title,
"Energy Measured in EB (GeV)");
200 eEB_->GetYaxis()->SetTitle(
"Events");
201 sprintf(
title,
"Energy Measured in EB+HB (GeV)");
204 eEBHB_->GetYaxis()->SetTitle(
"Events");
205 sprintf(
title,
"Energy Measured in EB+HB+HO (GeV)");
208 eEBHBHO_->GetYaxis()->SetTitle(
"Events");
209 sprintf(
title,
"Energy Measured in EB (GeV) for t < %d ns", itcut);
212 eEBT_->GetYaxis()->SetTitle(
"Events");
213 sprintf(
title,
"Energy Measured in EB+HB (GeV) for t < %d ns", itcut);
216 eEBHBT_->GetYaxis()->SetTitle(
"Events");
217 sprintf(
title,
"Energy Measured in EB+HB+HO (GeV) for t < %d ns", itcut);
220 eEBHBHOT_->GetYaxis()->SetTitle(
"Events");
221 sprintf(
title,
"SimHit energy in HO");
224 eHO1_->GetYaxis()->SetTitle(
"Events");
225 eHO2_ =
tfile->make<TProfile2D>(
"EHO2",
title, 30, -1.305, 1.305, 72, -3.1415926, 3.1415926);
227 eHO2_->GetYaxis()->SetTitle(
"Events");
228 sprintf(
title,
"SimHit energy in HO Layer 17");
231 eHO17_->GetYaxis()->SetTitle(
"Events");
232 sprintf(
title,
"SimHit energy in HO Layer 18");
235 eHO18_->GetYaxis()->SetTitle(
"Events");
236 sprintf(
title,
"SimHit energy in HO for t < %d ns", itcut);
239 eHO1T_->GetYaxis()->SetTitle(
"Events");
240 eHO2T_ =
tfile->make<TProfile2D>(
"EHO2T",
title, 30, -1.305, 1.305, 72, -3.1415926, 3.1415926);
242 eHO2T_->GetYaxis()->SetTitle(
"Events");
243 sprintf(
title,
"SimHit energy in HO Layer 17 for t < %d ns", itcut);
246 eHO17T_->GetYaxis()->SetTitle(
"Events");
247 sprintf(
title,
"SimHit energy in HO Layer 18 for t < %d ns", itcut);
250 eHO18T_->GetYaxis()->SetTitle(
"Events");
251 sprintf(
title,
"Number of layers hit in HO");
254 nHO1_->GetYaxis()->SetTitle(
"Events");
255 nHO2_ =
tfile->make<TProfile2D>(
"NHO2",
title, 30, -1.305, 1.305, 72, -3.1415926, 3.1415926);
257 nHO2_->GetYaxis()->SetTitle(
"Events");
258 sprintf(
title,
"Number of layers hit in HO for t < %d ns", itcut);
261 nHO1T_->GetYaxis()->SetTitle(
"Events");
262 nHO2T_ =
tfile->make<TProfile2D>(
"NHO2T",
title, 30, -1.305, 1.305, 72, -3.1415926, 3.1415926);
264 nHO2T_->GetYaxis()->SetTitle(
"Events");
265 for (
int i = 0;
i < 15;
i++) {
266 sprintf(
name,
"EHOE%d",
i + 1);
267 sprintf(
title,
"SimHit energy in HO (Beam in #eta=%d bin)",
i + 1);
270 eHOE_[
i]->GetYaxis()->SetTitle(
"Events");
271 sprintf(
name,
"EHOE17%d",
i + 1);
272 sprintf(
title,
"SimHit energy in Layer 17 (Beam in #eta=%d bin)",
i + 1);
275 eHOE17_[
i]->GetYaxis()->SetTitle(
"Events");
276 sprintf(
name,
"EHOE18%d",
i + 1);
277 sprintf(
title,
"SimHit energy in Layer 18 (Beam in #eta=%d bin)",
i + 1);
280 eHOE18_[
i]->GetYaxis()->SetTitle(
"Events");
281 sprintf(
name,
"EHOE%dT",
i + 1);
282 sprintf(
title,
"SimHit energy in HO (Beam in #eta=%d bin, t < %d ns)",
i + 1, itcut);
285 eHOET_[
i]->GetYaxis()->SetTitle(
"Events");
286 sprintf(
name,
"EHOE17%dT",
i + 1);
287 sprintf(
title,
"SimHit energy in Layer 17 (Beam in #eta=%d bin, t < %d ns)",
i + 1, itcut);
290 eHOE17T_[
i]->GetYaxis()->SetTitle(
"Events");
291 sprintf(
name,
"EHOE18%dT",
i + 1);
292 sprintf(
title,
"SimHit energy in Layer 18 (Beam in #eta=%d bin, t < %d ns)",
i + 1, itcut);
295 eHOE18T_[
i]->GetYaxis()->SetTitle(
"Events");
296 sprintf(
name,
"EHOEta%d",
i + 1);
297 sprintf(
title,
"SimHit energy in HO #eta bin %d (Beam in #eta=%d bin)",
i + 1,
i + 1);
300 eHOEta_[
i]->GetYaxis()->SetTitle(
"Events");
301 sprintf(
name,
"EHOEta17%d",
i + 1);
302 sprintf(
title,
"SimHit energy in Layer 17 #eta bin %d (Beam in #eta=%d bin)",
i + 1,
i + 1);
306 sprintf(
name,
"EHOEta18%d",
i + 1);
307 sprintf(
title,
"SimHit energy in Layer 18 #eta bin %d (Beam in #eta=%d bin)",
i + 1,
i + 1);
311 sprintf(
name,
"EHOEta%dT",
i + 1);
312 sprintf(
title,
"SimHit energy in HO #eta bin %d (Beam in #eta=%d bin, t < %d ns)",
i + 1,
i + 1, itcut);
315 eHOEtaT_[
i]->GetYaxis()->SetTitle(
"Events");
316 sprintf(
name,
"EHOEta17%dT",
i + 1);
317 sprintf(
title,
"SimHit energy in Layer 17 #eta bin %d (Beam in #eta=%d bin, t < %d ns)",
i + 1,
i + 1, itcut);
321 sprintf(
name,
"EHOEta18%dT",
i + 1);
322 sprintf(
title,
"SimHit energy in Layer 18 #eta bin %d (Beam in #eta=%d bin, t < %d ns)",
i + 1,
i + 1, itcut);
326 sprintf(
name,
"NHOE1%d",
i + 1);
327 sprintf(
title,
"Number of layers hit in HO (Beam in #eta=%d bin)",
i + 1);
330 nHOE1_[
i]->GetYaxis()->SetTitle(
"Events");
331 sprintf(
name,
"NHOE2%d",
i + 1);
334 nHOE2_[
i]->GetYaxis()->SetTitle(
"Events");
335 sprintf(
name,
"NHOE1%dT",
i + 1);
336 sprintf(
title,
"Number of layers hit in HO (Beam in #eta=%d bin, t < %d ns)",
i + 1, itcut);
339 nHOE1T_[
i]->GetYaxis()->SetTitle(
"Events");
340 sprintf(
name,
"NHOE2%dT",
i + 1);
343 nHOE2T_[
i]->GetYaxis()->SetTitle(
"Events");
344 sprintf(
name,
"NHOEta1%d",
i + 1);
345 sprintf(
title,
"Number of layers hit in HO #eta bin %d (Beam in #eta=%d bin)",
i + 1,
i + 1);
348 nHOEta1_[
i]->GetYaxis()->SetTitle(
"Events");
349 sprintf(
name,
"NHOEta2%d",
i + 1);
352 nHOEta2_[
i]->GetYaxis()->SetTitle(
"Events");
353 sprintf(
name,
"NHOEta1%dT",
i + 1);
354 sprintf(
title,
"Number of layers hit in HO #eta bin %d (Beam in #eta=%d bin, t < %d ns)",
i + 1,
i + 1, itcut);
358 sprintf(
name,
"NHOEta2%dT",
i + 1);
366 edm::LogVerbatim(
"HitStudy") <<
"HOSimHitStudy::Run = " <<
e.id().run() <<
" Event = " <<
e.id().event();
373 HepMC::GenEvent::particle_const_iterator
p = myGenEvent->particles_begin();
374 if (
p != myGenEvent->particles_end()) {
375 eInc = (*p)->momentum().e();
376 etaInc = (*p)->momentum().eta();
377 phiInc = (*p)->momentum().phi();
381 <<
" Phi = " <<
phiInc / CLHEP::deg;
383 for (
int i = 0;
i < 2;
i++) {
384 bool getHits =
false;
413 double etot[3], etotT[3];
414 std::vector<unsigned int> ebID, hbID, hoID;
415 std::vector<double> ebEtow, hbEtow, hoEtow;
416 std::vector<double> ebEtowT, hbEtowT, hoEtowT;
417 for (
int k = 0;
k < 3;
k++) {
426 double eHO17 = 0, eHO18 = 0, eHO17T = 0, eHO18T = 0;
427 double eHOE17[15], eHOE18[15], eHOE17T[15], eHOE18T[15];
428 for (
int k = 0;
k < 15;
k++) {
429 eHOE17[
k] = eHOE18[
k] = eHOE17T[
k] = eHOE18T[
k] = 0;
432 for (
int k = 0;
k < 2;
k++) {
439 for (
int i = 0;
i < nHit;
i++) {
457 subdet =
id.subdet();
465 <<
" depth:" <<
depth <<
" layer:" << lay <<
" eta:" <<
eta <<
" phi:" <<
phi;
468 else if (subdet == static_cast<int>(
HcalOuter)) {
474 if (
eta >= 0 &&
eta < 15) {
477 eHOE17T[
eta] += edep;
483 if (
eta >= 0 &&
eta < 15) {
486 eHOE18T[
eta] += edep;
494 edep_[indx]->Fill(edep);
504 for (
unsigned int j = 0;
j < ebID.size();
j++) {
505 if (id_ == ebID[
j]) {
514 ebEtow.push_back(edep);
515 ebEtowT.push_back(edepT);
517 }
else if (indx == 1) {
519 for (
unsigned int j = 0;
j < hbID.size();
j++) {
520 if (id_ == hbID[
j]) {
529 hbEtow.push_back(edep);
530 hbEtowT.push_back(edepT);
534 for (
unsigned int j = 0;
j < hoID.size();
j++) {
535 if (id_ == hoID[
j]) {
544 hoEtow.push_back(edep);
545 hoEtowT.push_back(edepT);
553 for (
int k = 0;
k < 3;
k++) {
554 hit_[
k]->Fill(
double(nhit[
k]));
558 hitTow_[0]->Fill(
double(ebEtow.size()));
559 for (
unsigned int i = 0;
i < ebEtow.size();
i++) {
563 hitTow_[1]->Fill(
double(hbEtow.size()));
564 for (
unsigned int i = 0;
i < hbEtow.size();
i++) {
568 hitTow_[2]->Fill(
double(hoEtow.size()));
569 for (
unsigned int i = 0;
i < hoEtow.size();
i++) {
573 double eEB =
scaleEB * etot[0];
574 double eEBHB = eEB +
scaleHB * etot[1];
579 double eEBT =
scaleEB * etotT[0];
580 double eEBHBT = eEBT +
scaleHB * etotT[1];
593 int nHO = 0, nHOT = 0;
607 for (
int k = 0;
k < 15; ++
k) {
627 int nHOE = 0, nHOET = 0;
628 if (eHOE17[
ieta] > 0)
630 if (eHOE17T[
ieta] > 0)
632 if (eHOE18[
ieta] > 0)
634 if (eHOE18T[
ieta] > 0)
642 edm::LogVerbatim(
"HitStudy") <<
"HOSimHitStudy::analyzeHits: Hits in EB " << nhit[0] <<
" in " << ebEtow.size()
643 <<
" towers with total E " << etot[0] <<
"|" << etotT[0]
644 <<
"\n Hits in HB " << nhit[1] <<
" in " << hbEtow.size()
645 <<
" towers with total E " << etot[1] <<
"|" << etotT[1]
646 <<
"\n Hits in HO " << nhit[2] <<
" in " << hoEtow.size()
647 <<
" towers with total E " << etot[2] <<
"|" << etotT[2]
648 <<
"\n Energy in EB " << eEB <<
"|" << eEBT <<
" with HB "
649 << eEBHB <<
"|" << eEBHBT <<
" and with HO " << eEBHBHO <<
"|" << eEBHBHOT
650 <<
"\n E in HO layers " << eHO17 <<
"|" << eHO17T <<
" "
651 << eHO18 <<
"|" << eHO18T <<
" number of HO hits " << nHO <<
"|" << nHOT;
654 edm::LogVerbatim(
"HitStudy") <<
"HOSimHitStudy::Event with excessive energy: EB = " << eEB <<
" EB+HB = " << eEBHB
655 <<
" EB+HB+HO = " << eEBHBHO;
657 for (
int k = 0;
k < 2;
k++) {
664 for (
int i = 0;
i < nHit;
i++) {
683 int subdet =
id.subdet();
686 else if (subdet == static_cast<int>(
HcalOuter))
694 <<
" (" <<
ieta <<
"|" <<
iphi <<
"|" <<
depth <<
") " << std::setw(8) << edep
695 <<
" " << std::setw(8) <<
time;