File indexing completed on 2024-10-08 05:12:11
0001 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0002 #include "Geometry/HcalCommonData/interface/HcalHitRelabeller.h"
0003 #include "Validation/HcalHits/interface/SimHitsValidationHcal.h"
0004
0005
0006
0007 SimHitsValidationHcal::SimHitsValidationHcal(const edm::ParameterSet &ps)
0008 : g4Label_(ps.getParameter<std::string>("ModuleLabel")),
0009 hcalHits_(ps.getParameter<std::string>("HitCollection")),
0010 verbose_(ps.getParameter<bool>("Verbose")),
0011 testNumber_(ps.getParameter<bool>("TestNumber")),
0012 tok_hits_(consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label_, hcalHits_))),
0013 tok_HRNDC_(esConsumes<HcalDDDRecConstants, HcalRecNumberingRecord, edm::Transition::BeginRun>()) {
0014 edm::LogVerbatim("HitsValidationHcal") << "Module Label: " << g4Label_ << " Hits: " << hcalHits_
0015 << " TestNumbering " << testNumber_;
0016 }
0017
0018 void SimHitsValidationHcal::bookHistograms(DQMStore::IBooker &ib, edm::Run const &run, edm::EventSetup const &es) {
0019 const auto &pHRNDC = es.getData(tok_HRNDC_);
0020 hcons = &pHRNDC;
0021 maxDepthHB_ = hcons->getMaxDepth(0);
0022 maxDepthHE_ = hcons->getMaxDepth(1);
0023 maxDepthHF_ = hcons->getMaxDepth(2);
0024 maxDepthHO_ = hcons->getMaxDepth(3);
0025
0026
0027
0028
0029 int NphiMax = hcons->getNPhi(0);
0030
0031 NphiMax = (hcons->getNPhi(1) > NphiMax ? hcons->getNPhi(1) : NphiMax);
0032 NphiMax = (hcons->getNPhi(2) > NphiMax ? hcons->getNPhi(2) : NphiMax);
0033 NphiMax = (hcons->getNPhi(3) > NphiMax ? hcons->getNPhi(3) : NphiMax);
0034
0035
0036 float iphi_min = 0.5;
0037 float iphi_max = NphiMax + 0.5;
0038 int iphi_bins = (int)(iphi_max - iphi_min);
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074 #ifdef EDM_ML_DEBUG
0075 edm::LogVerbatim("HitsValidationHcal") << " Maximum Depths HB:" << maxDepthHB_ << " HE:" << maxDepthHE_
0076 << " HO:" << maxDepthHO_ << " HF:" << maxDepthHF_;
0077 #endif
0078 std::vector<std::pair<std::string, std::string>> divisions = getHistogramTypes();
0079
0080 edm::LogVerbatim("HitsValidationHcal") << "Booking the Histograms";
0081 ib.setCurrentFolder("HcalHitsV/SimHitsValidationHcal");
0082
0083
0084
0085 std::string name, title;
0086 for (unsigned int i = 0; i < types.size(); ++i) {
0087 etaRange limit = getLimits(types[i]);
0088 name = "HcalHitEta" + divisions[i].first;
0089 title = "Hit energy as a function of eta tower index in " + divisions[i].second;
0090 meHcalHitEta_.push_back(ib.book1D(name, title, limit.bins, limit.low, limit.high));
0091
0092 name = "HcalHitTimeAEta" + divisions[i].first;
0093 title = "Hit time as a function of eta tower index in" + divisions[i].second;
0094 meHcalHitTimeEta_.push_back(ib.book1D(name, title, limit.bins, limit.low, limit.high));
0095
0096 name = "HcalHitE25" + divisions[i].first;
0097 title = "Energy in time window 0 to 25 for a tower in " + divisions[i].second;
0098 meHcalEnergyl25_.push_back(
0099 ib.book2D(name, title, limit.bins, limit.low, limit.high, iphi_bins, iphi_min, iphi_max));
0100
0101 name = "HcalHitE50" + divisions[i].first;
0102 title = "Energy in time window 0 to 50 for a tower in " + divisions[i].second;
0103 meHcalEnergyl50_.push_back(
0104 ib.book2D(name, title, limit.bins, limit.low, limit.high, iphi_bins, iphi_min, iphi_max));
0105
0106 name = "HcalHitE100" + divisions[i].first;
0107 title = "Energy in time window 0 to 100 for a tower in " + divisions[i].second;
0108 meHcalEnergyl100_.push_back(
0109 ib.book2D(name, title, limit.bins, limit.low, limit.high, iphi_bins, iphi_min, iphi_max));
0110
0111 name = "HcalHitE250" + divisions[i].first;
0112 title = "Energy in time window 0 to 250 for a tower in " + divisions[i].second;
0113 meHcalEnergyl250_.push_back(
0114 ib.book2D(name, title, limit.bins, limit.low, limit.high, iphi_bins, iphi_min, iphi_max));
0115 }
0116
0117 name = "Energy_HB";
0118 meEnergy_HB = ib.book1D(name, name, 100, 0, 1);
0119 name = "Energy_HE";
0120 meEnergy_HE = ib.book1D(name, name, 100, 0, 1);
0121 name = "Energy_HO";
0122 meEnergy_HO = ib.book1D(name, name, 100, 0, 1);
0123 name = "Energy_HF";
0124 meEnergy_HF = ib.book1D(name, name, 100, 0, 50);
0125
0126 name = "Time_HB";
0127 metime_HB = ib.book1D(name, name, 300, -150, 150);
0128 name = "Time_HE";
0129 metime_HE = ib.book1D(name, name, 300, -150, 150);
0130 name = "Time_HO";
0131 metime_HO = ib.book1D(name, name, 300, -150, 150);
0132 name = "Time_HF";
0133 metime_HF = ib.book1D(name, name, 300, -150, 150);
0134
0135 name = "Time_Enweighted_HB";
0136 metime_enweighted_HB = ib.book1D(name, name, 300, -150, 150);
0137 name = "Time_Enweighted_HE";
0138 metime_enweighted_HE = ib.book1D(name, name, 300, -150, 150);
0139 name = "Time_Enweighted_HO";
0140 metime_enweighted_HO = ib.book1D(name, name, 300, -150, 150);
0141 name = "Time_Enweighted_HF";
0142 metime_enweighted_HF = ib.book1D(name, name, 300, -150, 150);
0143 }
0144
0145 void SimHitsValidationHcal::analyze(const edm::Event &e, const edm::EventSetup &) {
0146 #ifdef EDM_ML_DEBUG
0147 edm::LogVerbatim("HitsValidationHcal") << "Run = " << e.id().run() << " Event = " << e.id().event();
0148 #endif
0149
0150 bool getHits = false;
0151 const edm::Handle<edm::PCaloHitContainer> &hitsHcal = e.getHandle(tok_hits_);
0152 if (hitsHcal.isValid())
0153 getHits = true;
0154 #ifdef EDM_ML_DEBUG
0155 edm::LogVerbatim("HitsValidationHcal") << "HitsValidationHcal.: Input flags Hits " << getHits;
0156 #endif
0157 if (getHits) {
0158 std::vector<PCaloHit> caloHits;
0159 caloHits.insert(caloHits.end(), hitsHcal->begin(), hitsHcal->end());
0160 #ifdef EDM_ML_DEBUG
0161 edm::LogVerbatim("HitsValidationHcal") << "testNumber_:" << testNumber_;
0162 #endif
0163 if (testNumber_) {
0164 for (unsigned int i = 0; i < caloHits.size(); ++i) {
0165 unsigned int id_ = caloHits[i].id();
0166 HcalDetId hid = HcalHitRelabeller::relabel(id_, hcons);
0167 caloHits[i].setID(hid.rawId());
0168 #ifdef EDM_ML_DEBUG
0169 edm::LogVerbatim("HitsValidationHcal") << "Hit[" << i << "] " << hid;
0170 #endif
0171 }
0172 }
0173 #ifdef EDM_ML_DEBUG
0174 edm::LogVerbatim("HitsValidationHcal") << "HitsValidationHcal: Hit buffer " << caloHits.size();
0175 #endif
0176 analyzeHits(caloHits);
0177 }
0178 }
0179
0180 void SimHitsValidationHcal::analyzeHits(std::vector<PCaloHit> &hits) {
0181 int nHit = hits.size();
0182 double entotHB = 0, entotHE = 0, entotHF = 0, entotHO = 0;
0183 double timetotHB = 0, timetotHE = 0, timetotHF = 0, timetotHO = 0;
0184
0185 std::map<std::pair<HcalDetId, unsigned int>, energysum> map_try;
0186 map_try.clear();
0187 std::map<std::pair<HcalDetId, unsigned int>, energysum>::iterator itr;
0188
0189 for (int i = 0; i < nHit; i++) {
0190 double energy = hits[i].energy();
0191 double time = hits[i].time();
0192 HcalDetId id = HcalDetId(hits[i].id());
0193 int itime = (int)(time);
0194 int subdet = id.subdet();
0195 int depth = id.depth();
0196 int eta = id.ieta();
0197 unsigned int dep = hits[i].depth();
0198
0199 std::pair<int, int> types = histId(subdet, eta, depth, dep);
0200 if (subdet == static_cast<int>(HcalBarrel)) {
0201 entotHB += energy;
0202 timetotHB += time;
0203 } else if (subdet == static_cast<int>(HcalEndcap)) {
0204 entotHE += energy;
0205 timetotHE += time;
0206 } else if (subdet == static_cast<int>(HcalOuter)) {
0207 entotHO += energy;
0208 timetotHO += time;
0209 } else if (subdet == static_cast<int>(HcalForward)) {
0210 entotHF += energy;
0211 timetotHF += time;
0212 }
0213
0214 std::pair<HcalDetId, unsigned int> id0(id, dep);
0215 energysum ensum;
0216 if (map_try.count(id0) != 0)
0217 ensum = map_try[id0];
0218 if (itime < 250) {
0219 ensum.e250 += energy;
0220 if (itime < 100) {
0221 ensum.e100 += energy;
0222 if (itime < 50) {
0223 ensum.e50 += energy;
0224 if (itime < 25)
0225 ensum.e25 += energy;
0226 }
0227 }
0228 }
0229 map_try[id0] = ensum;
0230
0231 #ifdef EDM_ML_DEBUG
0232 edm::LogVerbatim("HitsValidationHcal")
0233 << "Hit[" << i << "] ID " << std::dec << " " << id << std::dec << " Det " << id.det() << " Sub " << subdet
0234 << " depth " << depth << " depthX " << dep << " Eta " << eta << " Phi " << id.iphi() << " E " << energy
0235 << " time " << time << " type " << types.first << " " << types.second;
0236 #endif
0237
0238 double etax = eta - 0.5;
0239 if (eta < 0)
0240 etax += 1;
0241 if (types.first >= 0) {
0242 meHcalHitEta_[types.first]->Fill(etax, energy);
0243 meHcalHitTimeEta_[types.first]->Fill(etax, time);
0244 }
0245 if (types.second >= 0) {
0246 meHcalHitEta_[types.second]->Fill(etax, energy);
0247 meHcalHitTimeEta_[types.second]->Fill(etax, time);
0248 }
0249 }
0250
0251 meEnergy_HB->Fill(entotHB);
0252 meEnergy_HE->Fill(entotHE);
0253 meEnergy_HF->Fill(entotHF);
0254 meEnergy_HO->Fill(entotHO);
0255
0256 metime_HB->Fill(timetotHB);
0257 metime_HE->Fill(timetotHE);
0258 metime_HF->Fill(timetotHF);
0259 metime_HO->Fill(timetotHO);
0260
0261 metime_enweighted_HB->Fill(timetotHB, entotHB);
0262 metime_enweighted_HE->Fill(timetotHE, entotHE);
0263 metime_enweighted_HF->Fill(timetotHF, entotHF);
0264 metime_enweighted_HO->Fill(timetotHO, entotHO);
0265
0266 for (itr = map_try.begin(); itr != map_try.end(); ++itr) {
0267 HcalDetId id = (*itr).first.first;
0268 energysum ensum = (*itr).second;
0269 std::pair<int, int> types = histId((int)(id.subdet()), id.ieta(), id.depth(), (*itr).first.second);
0270 int eta = id.ieta();
0271 int phi = id.iphi();
0272 double etax = eta - 0.5;
0273 double phix = phi - 0.5;
0274 if (types.first >= 0) {
0275 meHcalEnergyl25_[types.first]->Fill(etax, phix, ensum.e25);
0276 meHcalEnergyl50_[types.first]->Fill(etax, phix, ensum.e50);
0277 meHcalEnergyl100_[types.first]->Fill(etax, phix, ensum.e100);
0278 meHcalEnergyl250_[types.first]->Fill(etax, phix, ensum.e250);
0279 }
0280 if (types.second >= 0) {
0281 meHcalEnergyl25_[types.second]->Fill(etax, phix, ensum.e25);
0282 meHcalEnergyl50_[types.second]->Fill(etax, phix, ensum.e50);
0283 meHcalEnergyl100_[types.second]->Fill(etax, phix, ensum.e100);
0284 meHcalEnergyl250_[types.second]->Fill(etax, phix, ensum.e250);
0285 }
0286
0287 #ifdef EDM_ML_DEBUG
0288 edm::LogVerbatim("HitsValidationHcal")
0289 << " energy of tower =" << (*itr).first.first << " in time 25ns is == " << (*itr).second.e25
0290 << " in time 25-50ns == " << (*itr).second.e50 << " in time 50-100ns == " << (*itr).second.e100
0291 << " in time 100-250 ns == " << (*itr).second.e250;
0292 #endif
0293 }
0294 }
0295
0296 SimHitsValidationHcal::etaRange SimHitsValidationHcal::getLimits(idType type) {
0297 int bins;
0298 std::pair<int, int> range;
0299 double low, high;
0300
0301 if (type.subdet == HcalBarrel) {
0302 range = hcons->getEtaRange(0);
0303 low = -range.second;
0304 high = range.second;
0305 bins = (high - low);
0306 } else if (type.subdet == HcalEndcap) {
0307 range = hcons->getEtaRange(1);
0308 bins = range.second - range.first;
0309 if (type.z == 1) {
0310 low = range.first;
0311 high = range.second;
0312 } else {
0313 low = -range.second;
0314 high = -range.first;
0315 }
0316 } else if (type.subdet == HcalOuter) {
0317 range = hcons->getEtaRange(3);
0318 low = -range.second;
0319 high = range.second;
0320 bins = high - low;
0321 } else if (type.subdet == HcalForward) {
0322 range = hcons->getEtaRange(2);
0323 bins = range.second - range.first;
0324 if (type.z == 1) {
0325 low = range.first;
0326 high = range.second;
0327 } else {
0328 low = -range.second;
0329 high = -range.first;
0330 }
0331 } else {
0332 bins = 82;
0333 low = -41;
0334 high = 41;
0335 }
0336 #ifdef EDM_ML_DEBUG
0337 edm::LogVerbatim("HitsValidationHcal") << "Subdetector:" << type.subdet << " z:" << type.z
0338 << " range.first:" << range.first << " and second:" << range.second;
0339 edm::LogVerbatim("HitsValidationHcal") << "bins: " << bins << " low:" << low << " high:" << high;
0340 #endif
0341 return SimHitsValidationHcal::etaRange(bins, low, high);
0342 }
0343
0344 std::pair<int, int> SimHitsValidationHcal::histId(int subdet, int eta, int depth, unsigned int dep) {
0345 int id1(-1), id2(-1);
0346 for (unsigned int k = 0; k < types.size(); ++k) {
0347 if (subdet == HcalForward) {
0348 if (subdet == (int)(types[k].subdet) && depth == types[k].depth1 && eta * types[k].z > 0 &&
0349 dep == (unsigned int)(types[k].depth2)) {
0350 id1 = k;
0351 break;
0352 }
0353 } else if (subdet == HcalEndcap) {
0354 if (subdet == (int)(types[k].subdet) && depth == types[k].depth1 && eta * types[k].z > 0) {
0355 id1 = k;
0356 break;
0357 }
0358 } else {
0359 if (subdet == (int)(types[k].subdet) && depth == types[k].depth1) {
0360 id1 = k;
0361 break;
0362 }
0363 }
0364 }
0365 if (subdet == HcalForward)
0366 depth += 2 * dep;
0367 for (unsigned int k = 0; k < types.size(); ++k) {
0368 if (types[k].subdet == HcalEmpty && types[k].depth1 == depth) {
0369 id2 = k;
0370 break;
0371 }
0372 }
0373 return std::pair<int, int>(id1, id2);
0374 }
0375
0376 std::vector<std::pair<std::string, std::string>> SimHitsValidationHcal::getHistogramTypes() {
0377 int maxDepth = std::max(maxDepthHB_, maxDepthHE_);
0378 maxDepth = std::max(maxDepth, maxDepthHF_);
0379 maxDepth = std::max(maxDepth, maxDepthHO_);
0380
0381 std::vector<std::pair<std::string, std::string>> divisions;
0382
0383 types.clear();
0384 std::pair<std::string, std::string> names;
0385 char name1[40], name2[40];
0386 SimHitsValidationHcal::idType type;
0387
0388 for (int depth = 0; depth < maxDepth; ++depth) {
0389 snprintf(name1, 40, "HC%d", depth);
0390 snprintf(name2, 40, "HCAL depth%d", depth + 1);
0391 names = std::pair<std::string, std::string>(std::string(name1), std::string(name2));
0392 type = SimHitsValidationHcal::idType(HcalEmpty, 0, depth + 1, depth + 1);
0393 divisions.push_back(names);
0394 types.push_back(type);
0395 }
0396
0397 for (int depth = 0; depth < maxDepthHB_; ++depth) {
0398 snprintf(name1, 40, "HB%d", depth);
0399 snprintf(name2, 40, "HB depth%d", depth + 1);
0400 names = std::pair<std::string, std::string>(std::string(name1), std::string(name2));
0401 type = SimHitsValidationHcal::idType(HcalBarrel, 0, depth + 1, depth + 1);
0402 divisions.push_back(names);
0403 types.push_back(type);
0404 }
0405
0406 for (int depth = 0; depth < maxDepthHE_; ++depth) {
0407 snprintf(name1, 40, "HE%d+z", depth);
0408 snprintf(name2, 40, "HE +z depth%d", depth + 1);
0409 names = std::pair<std::string, std::string>(std::string(name1), std::string(name2));
0410 type = SimHitsValidationHcal::idType(HcalEndcap, 1, depth + 1, depth + 1);
0411 divisions.push_back(names);
0412 types.push_back(type);
0413 snprintf(name1, 40, "HE%d-z", depth);
0414 snprintf(name2, 40, "HE -z depth%d", depth + 1);
0415 names = std::pair<std::string, std::string>(std::string(name1), std::string(name2));
0416 type = SimHitsValidationHcal::idType(HcalEndcap, -1, depth + 1, depth + 1);
0417 divisions.push_back(names);
0418 types.push_back(type);
0419 }
0420
0421 {
0422 int depth = maxDepthHO_;
0423 snprintf(name1, 40, "HO%d", depth);
0424 snprintf(name2, 40, "HO depth%d", depth);
0425 names = std::pair<std::string, std::string>(std::string(name1), std::string(name2));
0426 type = SimHitsValidationHcal::idType(HcalOuter, 0, depth, depth);
0427 divisions.push_back(names);
0428 types.push_back(type);
0429 }
0430
0431 std::string hfty1[4] = {"A", "W", "B", "J"};
0432 std::string hfty2[4] = {"Absorber", "Window", "Bundle", "Jungle"};
0433 int dept0[4] = {0, 1, 2, 3};
0434 for (int k = 0; k < 4; ++k) {
0435 for (int depth = 0; depth < maxDepthHF_; ++depth) {
0436 snprintf(name1, 40, "HF%s%d+z", hfty1[k].c_str(), depth);
0437 snprintf(name2, 40, "HF (%s) +z depth%d", hfty2[k].c_str(), depth + 1);
0438 names = std::pair<std::string, std::string>(std::string(name1), std::string(name2));
0439 type = SimHitsValidationHcal::idType(HcalForward, 1, depth + 1, dept0[k]);
0440 divisions.push_back(names);
0441 types.push_back(type);
0442 snprintf(name1, 40, "HF%s%d-z", hfty1[k].c_str(), depth);
0443 snprintf(name2, 40, "HF (%s) -z depth%d", hfty2[k].c_str(), depth + 1);
0444 names = std::pair<std::string, std::string>(std::string(name1), std::string(name2));
0445 type = SimHitsValidationHcal::idType(HcalForward, -1, depth + 1, dept0[k]);
0446 divisions.push_back(names);
0447 types.push_back(type);
0448 }
0449 }
0450
0451 return divisions;
0452 }
0453
0454 void SimHitsValidationHcal::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0455 edm::ParameterSetDescription desc;
0456 desc.add<std::string>("ModuleLabel", "g4SimHits");
0457 desc.add<std::string>("HitCollection", "HcalHits");
0458 desc.add<bool>("Verbose", false);
0459 desc.add<bool>("TestNumber", false);
0460
0461 descriptions.addDefault(desc);
0462 }
0463
0464 DEFINE_FWK_MODULE(SimHitsValidationHcal);