25 #include "CLHEP/Units/GlobalPhysicalConstants.h"
26 #include "CLHEP/Units/GlobalSystemOfUnits.h"
30 #include <TProfile2D.h>
88 for (
unsigned i = 0;
i != 2;
i++)
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();
372 HepMC::GenEvent::particle_const_iterator
p = myGenEvent->particles_begin();
373 if (
p != myGenEvent->particles_end()) {
374 eInc = (*p)->momentum().e();
375 etaInc = (*p)->momentum().eta();
376 phiInc = (*p)->momentum().phi();
380 <<
" Phi = " <<
phiInc / CLHEP::deg;
382 for (
int i = 0;
i < 2;
i++) {
383 bool getHits =
false;
412 double etot[3], etotT[3];
413 std::vector<unsigned int> ebID, hbID, hoID;
414 std::vector<double> ebEtow, hbEtow, hoEtow;
415 std::vector<double> ebEtowT, hbEtowT, hoEtowT;
416 for (
int k = 0;
k < 3;
k++) {
425 double eHO17 = 0, eHO18 = 0, eHO17T = 0, eHO18T = 0;
426 double eHOE17[15], eHOE18[15], eHOE17T[15], eHOE18T[15];
427 for (
int k = 0;
k < 15;
k++) {
428 eHOE17[
k] = eHOE18[
k] = eHOE17T[
k] = eHOE18T[
k] = 0;
431 for (
int k = 0;
k < 2;
k++) {
438 for (
int i = 0;
i < nHit;
i++) {
456 subdet =
id.subdet();
464 <<
" depth:" <<
depth <<
" layer:" << lay <<
" eta:" <<
eta <<
" phi:" <<
phi;
467 else if (subdet == static_cast<int>(
HcalOuter)) {
473 if (
eta >= 0 &&
eta < 15) {
476 eHOE17T[
eta] += edep;
482 if (
eta >= 0 &&
eta < 15) {
485 eHOE18T[
eta] += edep;
493 edep_[indx]->Fill(edep);
503 for (
unsigned int j = 0;
j < ebID.size();
j++) {
504 if (id_ == ebID[
j]) {
513 ebEtow.push_back(edep);
514 ebEtowT.push_back(edepT);
516 }
else if (indx == 1) {
518 for (
unsigned int j = 0;
j < hbID.size();
j++) {
519 if (id_ == hbID[
j]) {
528 hbEtow.push_back(edep);
529 hbEtowT.push_back(edepT);
533 for (
unsigned int j = 0;
j < hoID.size();
j++) {
534 if (id_ == hoID[
j]) {
543 hoEtow.push_back(edep);
544 hoEtowT.push_back(edepT);
552 for (
int k = 0;
k < 3;
k++) {
553 hit_[
k]->Fill(
double(nhit[
k]));
557 hitTow_[0]->Fill(
double(ebEtow.size()));
558 for (
unsigned int i = 0;
i < ebEtow.size();
i++) {
562 hitTow_[1]->Fill(
double(hbEtow.size()));
563 for (
unsigned int i = 0;
i < hbEtow.size();
i++) {
567 hitTow_[2]->Fill(
double(hoEtow.size()));
568 for (
unsigned int i = 0;
i < hoEtow.size();
i++) {
572 double eEB =
scaleEB * etot[0];
573 double eEBHB = eEB +
scaleHB * etot[1];
578 double eEBT =
scaleEB * etotT[0];
579 double eEBHBT = eEBT +
scaleHB * etotT[1];
592 int nHO = 0, nHOT = 0;
606 for (
int k = 0;
k < 15; ++
k) {
626 int nHOE = 0, nHOET = 0;
627 if (eHOE17[
ieta] > 0)
629 if (eHOE17T[
ieta] > 0)
631 if (eHOE18[
ieta] > 0)
633 if (eHOE18T[
ieta] > 0)
641 edm::LogVerbatim(
"HitStudy") <<
"HOSimHitStudy::analyzeHits: Hits in EB " << nhit[0] <<
" in " << ebEtow.size()
642 <<
" towers with total E " << etot[0] <<
"|" << etotT[0]
643 <<
"\n Hits in HB " << nhit[1] <<
" in " << hbEtow.size()
644 <<
" towers with total E " << etot[1] <<
"|" << etotT[1]
645 <<
"\n Hits in HO " << nhit[2] <<
" in " << hoEtow.size()
646 <<
" towers with total E " << etot[2] <<
"|" << etotT[2]
647 <<
"\n Energy in EB " << eEB <<
"|" << eEBT <<
" with HB "
648 << eEBHB <<
"|" << eEBHBT <<
" and with HO " << eEBHBHO <<
"|" << eEBHBHOT
649 <<
"\n E in HO layers " << eHO17 <<
"|" << eHO17T <<
" "
650 << eHO18 <<
"|" << eHO18T <<
" number of HO hits " << nHO <<
"|" << nHOT;
653 edm::LogVerbatim(
"HitStudy") <<
"HOSimHitStudy::Event with excessive energy: EB = " << eEB <<
" EB+HB = " << eEBHB
654 <<
" EB+HB+HO = " << eEBHBHO;
656 for (
int k = 0;
k < 2;
k++) {
663 for (
int i = 0;
i < nHit;
i++) {
682 int subdet =
id.subdet();
685 else if (subdet == static_cast<int>(
HcalOuter))
693 <<
" (" <<
ieta <<
"|" <<
iphi <<
"|" <<
depth <<
") " << std::setw(8) << edep
694 <<
" " << std::setw(8) <<
time;