Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-03-11 02:41:14

0001 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0002 #include "Geometry/HcalCommonData/interface/HcalHitRelabeller.h"
0003 #include "Validation/HcalHits/interface/SimHitsValidationHcal.h"
0004 
0005 //#define EDM_ML_DEBUG
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   // Get Phi segmentation from geometry, use the max phi number so that all iphi
0027   // values are included.
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   // Center the iphi bins on the integers
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   // Retain classic behavior, all plots have same ieta range.
0046   // Comment out    code to allow each subdetector to have its on range
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   // Give an empty bin around the subdet ieta range to make it clear that all
0058   // ieta rings have been included float ieta_min_HB = -iEtaHBMax - 1.5; float
0059   // ieta_max_HB = iEtaHBMax + 1.5; int ieta_bins_HB = (int) (ieta_max_HB -
0060   // ieta_min_HB);
0061 
0062   // float ieta_min_HE = -iEtaHEMax - 1.5;
0063   // float ieta_max_HE = iEtaHEMax + 1.5;
0064   // int ieta_bins_HE = (int) (ieta_max_HE - ieta_min_HE);
0065 
0066   // float ieta_min_HF = -iEtaHFMax - 1.5;
0067   // float ieta_max_HF = iEtaHFMax + 1.5;
0068   // int ieta_bins_HF = (int) (ieta_max_HF - ieta_min_HF);
0069 
0070   // float ieta_min_HO = -iEtaHOMax - 1.5;
0071   // float ieta_max_HO = iEtaHOMax + 1.5;
0072   // int ieta_bins_HO = (int) (ieta_max_HO - ieta_min_HO);
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   // Histograms for Hits
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   int nHB = 0, nHE = 0, nHO = 0, nHF = 0;
0185 
0186   std::map<std::pair<HcalDetId, unsigned int>, energysum> map_try;
0187   map_try.clear();
0188   std::map<std::pair<HcalDetId, unsigned int>, energysum>::iterator itr;
0189 
0190   for (int i = 0; i < nHit; i++) {
0191     double energy = hits[i].energy();
0192     double time = hits[i].time();
0193     HcalDetId id = HcalDetId(hits[i].id());
0194     int itime = (int)(time);
0195     int subdet = id.subdet();
0196     int depth = id.depth();
0197     int eta = id.ieta();
0198     unsigned int dep = hits[i].depth();
0199 
0200     std::pair<int, int> types = histId(subdet, eta, depth, dep);
0201     if (subdet == static_cast<int>(HcalBarrel)) {
0202       entotHB += energy;
0203       timetotHB += time;
0204       nHB++;
0205     } else if (subdet == static_cast<int>(HcalEndcap)) {
0206       entotHE += energy;
0207       timetotHE += time;
0208       nHE++;
0209     } else if (subdet == static_cast<int>(HcalOuter)) {
0210       entotHO += energy;
0211       timetotHO += time;
0212       nHO++;
0213     } else if (subdet == static_cast<int>(HcalForward)) {
0214       entotHF += energy;
0215       timetotHF += time;
0216       nHF++;
0217     }
0218 
0219     std::pair<HcalDetId, unsigned int> id0(id, dep);
0220     energysum ensum;
0221     if (map_try.count(id0) != 0)
0222       ensum = map_try[id0];
0223     if (itime < 250) {
0224       ensum.e250 += energy;
0225       if (itime < 100) {
0226         ensum.e100 += energy;
0227         if (itime < 50) {
0228           ensum.e50 += energy;
0229           if (itime < 25)
0230             ensum.e25 += energy;
0231         }
0232       }
0233     }
0234     map_try[id0] = ensum;
0235 
0236 #ifdef EDM_ML_DEBUG
0237     edm::LogVerbatim("HitsValidationHcal")
0238         << "Hit[" << i << "] ID " << std::dec << " " << id << std::dec << " Det " << id.det() << " Sub " << subdet
0239         << " depth " << depth << " depthX " << dep << " Eta " << eta << " Phi " << id.iphi() << " E " << energy
0240         << " time " << time << " type " << types.first << " " << types.second;
0241 #endif
0242 
0243     double etax = eta - 0.5;
0244     if (eta < 0)
0245       etax += 1;
0246     if (types.first >= 0) {
0247       meHcalHitEta_[types.first]->Fill(etax, energy);
0248       meHcalHitTimeEta_[types.first]->Fill(etax, time);
0249     }
0250     if (types.second >= 0) {
0251       meHcalHitEta_[types.second]->Fill(etax, energy);
0252       meHcalHitTimeEta_[types.second]->Fill(etax, time);
0253     }
0254   }
0255 
0256   meEnergy_HB->Fill(entotHB);
0257   meEnergy_HE->Fill(entotHE);
0258   meEnergy_HF->Fill(entotHF);
0259   meEnergy_HO->Fill(entotHO);
0260 
0261   metime_HB->Fill(timetotHB);
0262   metime_HE->Fill(timetotHE);
0263   metime_HF->Fill(timetotHF);
0264   metime_HO->Fill(timetotHO);
0265 
0266   metime_enweighted_HB->Fill(timetotHB, entotHB);
0267   metime_enweighted_HE->Fill(timetotHE, entotHE);
0268   metime_enweighted_HF->Fill(timetotHF, entotHF);
0269   metime_enweighted_HO->Fill(timetotHO, entotHO);
0270 
0271   for (itr = map_try.begin(); itr != map_try.end(); ++itr) {
0272     HcalDetId id = (*itr).first.first;
0273     energysum ensum = (*itr).second;
0274     std::pair<int, int> types = histId((int)(id.subdet()), id.ieta(), id.depth(), (*itr).first.second);
0275     int eta = id.ieta();
0276     int phi = id.iphi();
0277     double etax = eta - 0.5;
0278     double phix = phi - 0.5;
0279     if (types.first >= 0) {
0280       meHcalEnergyl25_[types.first]->Fill(etax, phix, ensum.e25);
0281       meHcalEnergyl50_[types.first]->Fill(etax, phix, ensum.e50);
0282       meHcalEnergyl100_[types.first]->Fill(etax, phix, ensum.e100);
0283       meHcalEnergyl250_[types.first]->Fill(etax, phix, ensum.e250);
0284     }
0285     if (types.second >= 0) {
0286       meHcalEnergyl25_[types.second]->Fill(etax, phix, ensum.e25);
0287       meHcalEnergyl50_[types.second]->Fill(etax, phix, ensum.e50);
0288       meHcalEnergyl100_[types.second]->Fill(etax, phix, ensum.e100);
0289       meHcalEnergyl250_[types.second]->Fill(etax, phix, ensum.e250);
0290     }
0291 
0292 #ifdef EDM_ML_DEBUG
0293     edm::LogVerbatim("HitsValidationHcal")
0294         << " energy of tower =" << (*itr).first.first << " in time 25ns is == " << (*itr).second.e25
0295         << " in time 25-50ns == " << (*itr).second.e50 << " in time 50-100ns == " << (*itr).second.e100
0296         << " in time 100-250 ns == " << (*itr).second.e250;
0297 #endif
0298   }
0299 }
0300 
0301 SimHitsValidationHcal::etaRange SimHitsValidationHcal::getLimits(idType type) {
0302   int bins;
0303   std::pair<int, int> range;
0304   double low, high;
0305 
0306   if (type.subdet == HcalBarrel) {
0307     range = hcons->getEtaRange(0);
0308     low = -range.second;
0309     high = range.second;
0310     bins = (high - low);
0311   } else if (type.subdet == HcalEndcap) {
0312     range = hcons->getEtaRange(1);
0313     bins = range.second - range.first;
0314     if (type.z == 1) {
0315       low = range.first;
0316       high = range.second;
0317     } else {
0318       low = -range.second;
0319       high = -range.first;
0320     }
0321   } else if (type.subdet == HcalOuter) {
0322     range = hcons->getEtaRange(3);
0323     low = -range.second;
0324     high = range.second;
0325     bins = high - low;
0326   } else if (type.subdet == HcalForward) {
0327     range = hcons->getEtaRange(2);
0328     bins = range.second - range.first;
0329     if (type.z == 1) {
0330       low = range.first;
0331       high = range.second;
0332     } else {
0333       low = -range.second;
0334       high = -range.first;
0335     }
0336   } else {
0337     bins = 82;
0338     low = -41;
0339     high = 41;
0340   }
0341 #ifdef EDM_ML_DEBUG
0342   edm::LogVerbatim("HitsValidationHcal") << "Subdetector:" << type.subdet << " z:" << type.z
0343                                          << " range.first:" << range.first << " and second:" << range.second;
0344   edm::LogVerbatim("HitsValidationHcal") << "bins: " << bins << " low:" << low << " high:" << high;
0345 #endif
0346   return SimHitsValidationHcal::etaRange(bins, low, high);
0347 }
0348 
0349 std::pair<int, int> SimHitsValidationHcal::histId(int subdet, int eta, int depth, unsigned int dep) {
0350   int id1(-1), id2(-1);
0351   for (unsigned int k = 0; k < types.size(); ++k) {
0352     if (subdet == HcalForward) {
0353       if (subdet == (int)(types[k].subdet) && depth == types[k].depth1 && eta * types[k].z > 0 &&
0354           dep == (unsigned int)(types[k].depth2)) {
0355         id1 = k;
0356         break;
0357       }
0358     } else if (subdet == HcalEndcap) {
0359       if (subdet == (int)(types[k].subdet) && depth == types[k].depth1 && eta * types[k].z > 0) {
0360         id1 = k;
0361         break;
0362       }
0363     } else {
0364       if (subdet == (int)(types[k].subdet) && depth == types[k].depth1) {
0365         id1 = k;
0366         break;
0367       }
0368     }
0369   }
0370   if (subdet == HcalForward)
0371     depth += 2 * dep;
0372   for (unsigned int k = 0; k < types.size(); ++k) {
0373     if (types[k].subdet == HcalEmpty && types[k].depth1 == depth) {
0374       id2 = k;
0375       break;
0376     }
0377   }
0378   return std::pair<int, int>(id1, id2);
0379 }
0380 
0381 std::vector<std::pair<std::string, std::string>> SimHitsValidationHcal::getHistogramTypes() {
0382   int maxDepth = std::max(maxDepthHB_, maxDepthHE_);
0383   maxDepth = std::max(maxDepth, maxDepthHF_);
0384   maxDepth = std::max(maxDepth, maxDepthHO_);
0385 
0386   std::vector<std::pair<std::string, std::string>> divisions;
0387   // divisions and types need to be in sync
0388   types.clear();
0389   std::pair<std::string, std::string> names;
0390   char name1[40], name2[40];
0391   SimHitsValidationHcal::idType type;
0392   // first overall Hcal
0393   for (int depth = 0; depth < maxDepth; ++depth) {
0394     snprintf(name1, 40, "HC%d", depth);
0395     snprintf(name2, 40, "HCAL depth%d", depth + 1);
0396     names = std::pair<std::string, std::string>(std::string(name1), std::string(name2));
0397     type = SimHitsValidationHcal::idType(HcalEmpty, 0, depth + 1, depth + 1);
0398     divisions.push_back(names);
0399     types.push_back(type);
0400   }
0401   // HB
0402   for (int depth = 0; depth < maxDepthHB_; ++depth) {
0403     snprintf(name1, 40, "HB%d", depth);
0404     snprintf(name2, 40, "HB depth%d", depth + 1);
0405     names = std::pair<std::string, std::string>(std::string(name1), std::string(name2));
0406     type = SimHitsValidationHcal::idType(HcalBarrel, 0, depth + 1, depth + 1);
0407     divisions.push_back(names);
0408     types.push_back(type);
0409   }
0410   // HE
0411   for (int depth = 0; depth < maxDepthHE_; ++depth) {
0412     snprintf(name1, 40, "HE%d+z", depth);
0413     snprintf(name2, 40, "HE +z depth%d", depth + 1);
0414     names = std::pair<std::string, std::string>(std::string(name1), std::string(name2));
0415     type = SimHitsValidationHcal::idType(HcalEndcap, 1, depth + 1, depth + 1);
0416     divisions.push_back(names);
0417     types.push_back(type);
0418     snprintf(name1, 40, "HE%d-z", depth);
0419     snprintf(name2, 40, "HE -z depth%d", depth + 1);
0420     names = std::pair<std::string, std::string>(std::string(name1), std::string(name2));
0421     type = SimHitsValidationHcal::idType(HcalEndcap, -1, depth + 1, depth + 1);
0422     divisions.push_back(names);
0423     types.push_back(type);
0424   }
0425   // HO
0426   {
0427     int depth = maxDepthHO_;
0428     snprintf(name1, 40, "HO%d", depth);
0429     snprintf(name2, 40, "HO depth%d", depth);
0430     names = std::pair<std::string, std::string>(std::string(name1), std::string(name2));
0431     type = SimHitsValidationHcal::idType(HcalOuter, 0, depth, depth);
0432     divisions.push_back(names);
0433     types.push_back(type);
0434   }
0435   // HF (first absorber, then different types of abnormal hits)
0436   std::string hfty1[4] = {"A", "W", "B", "J"};
0437   std::string hfty2[4] = {"Absorber", "Window", "Bundle", "Jungle"};
0438   int dept0[4] = {0, 1, 2, 3};
0439   for (int k = 0; k < 4; ++k) {
0440     for (int depth = 0; depth < maxDepthHF_; ++depth) {
0441       snprintf(name1, 40, "HF%s%d+z", hfty1[k].c_str(), depth);
0442       snprintf(name2, 40, "HF (%s) +z depth%d", hfty2[k].c_str(), depth + 1);
0443       names = std::pair<std::string, std::string>(std::string(name1), std::string(name2));
0444       type = SimHitsValidationHcal::idType(HcalForward, 1, depth + 1, dept0[k]);
0445       divisions.push_back(names);
0446       types.push_back(type);
0447       snprintf(name1, 40, "HF%s%d-z", hfty1[k].c_str(), depth);
0448       snprintf(name2, 40, "HF (%s) -z depth%d", hfty2[k].c_str(), depth + 1);
0449       names = std::pair<std::string, std::string>(std::string(name1), std::string(name2));
0450       type = SimHitsValidationHcal::idType(HcalForward, -1, depth + 1, dept0[k]);
0451       divisions.push_back(names);
0452       types.push_back(type);
0453     }
0454   }
0455 
0456   return divisions;
0457 }
0458 
0459 void SimHitsValidationHcal::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0460   edm::ParameterSetDescription desc;
0461   desc.add<std::string>("ModuleLabel", "g4SimHits");
0462   desc.add<std::string>("HitCollection", "HcalHits");
0463   desc.add<bool>("Verbose", false);
0464   desc.add<bool>("TestNumber", false);
0465 
0466   descriptions.addDefault(desc);
0467 }
0468 
0469 DEFINE_FWK_MODULE(SimHitsValidationHcal);