Back to home page

Project CMSSW displayed by LXR
 
 

    

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 //#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 
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   // divisions and types need to be in sync
0383   types.clear();
0384   std::pair<std::string, std::string> names;
0385   char name1[40], name2[40];
0386   SimHitsValidationHcal::idType type;
0387   // first overall Hcal
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   // HB
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   // HE
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   // HO
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   // HF (first absorber, then different types of abnormal hits)
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);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!