File indexing completed on 2023-05-11 23:51:49
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 int iEtaHBMax = hcons->getEtaRange(0).second;
0041 int iEtaHEMax = std::max(hcons->getEtaRange(1).second, 1);
0042 int iEtaHFMax = hcons->getEtaRange(2).second;
0043 int iEtaHOMax = hcons->getEtaRange(3).second;
0044
0045
0046
0047
0048 int iEtaMax = (iEtaHBMax > iEtaHEMax ? iEtaHBMax : iEtaHEMax);
0049 iEtaMax = (iEtaMax > iEtaHFMax ? iEtaMax : iEtaHFMax);
0050 iEtaMax = (iEtaMax > iEtaHOMax ? iEtaMax : iEtaHOMax);
0051
0052 iEtaHBMax = iEtaMax;
0053 iEtaHEMax = iEtaMax;
0054 iEtaHFMax = iEtaMax;
0055 iEtaHOMax = iEtaMax;
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);