25 #include <CLHEP/Units/GlobalPhysicalConstants.h>
26 #include <CLHEP/Units/GlobalSystemOfUnits.h>
30 #include <TProfile2D.h>
88 for (
unsigned i = 0;
i != 2;
i++)
98 edm::LogVerbatim(
"HitStudy") <<
"HOSimHitStudy::Module Label: " << g4Label <<
" Hits: " << hitLab[0] <<
", "
99 << hitLab[1] <<
" MaxEnergy: " <<
maxEnergy <<
" Scale factor for EB " <<
scaleEB
116 descriptions.
add(
"HOSimHitStudy", desc);
124 <<
"please add it to config file";
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);
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();
463 edm::LogVerbatim(
"HitStudy") <<
"HOSimHitStudy:: Hit " <<
k <<
" Subdet:" << subdet <<
" zside:" << zside
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;
492 time_[indx]->Fill(time);
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) {
612 if (ieta >= 0 && ieta < 15) {
613 eHOE_[ieta]->Fill(eHO17 + eHO18);
616 eHOET_[ieta]->Fill(eHO17T + eHO18T);
619 eHOEta_[ieta]->Fill(eHOE17[ieta] + eHOE18[ieta]);
622 nHOE1_[ieta]->Fill((
double)(nHO));
624 nHOE1T_[ieta]->Fill((
double)(nHOT));
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)
635 nHOEta1_[ieta]->Fill((
double)(nHOE));
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++) {
667 int ieta, iphi,
depth = 0;
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;
static const std::string kSharedResource
Log< level::Info, true > LogVerbatim
EventNumber_t event() const
T getUntrackedParameter(std::string const &, T const &) 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
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
bool getByToken(EDGetToken token, Handle< PROD > &result) const
#define DEFINE_FWK_MODULE(type)
edm::EDGetTokenT< edm::HepMCProduct > tok_evt_
T * make(const Args &...args) const
make new ROOT object
void endRun(edm::Run const &, edm::EventSetup const &) override
Abs< T >::type abs(const T &t)
~HOSimHitStudy() override
void add(std::string const &label, ParameterSetDescription const &psetDescription)
edm::EDGetTokenT< edm::PCaloHitContainer > toks_calo_[2]
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)
HOSimHitStudy(const edm::ParameterSet &ps)