File indexing completed on 2024-04-06 11:59:10
0001 #include "Calibration/IsolatedParticles/interface/DebugInfo.h"
0002 #include "Calibration/IsolatedParticles/interface/FindCaloHit.h"
0003 #include "Calibration/IsolatedParticles/interface/FindEtaPhi.h"
0004 #include "Calibration/IsolatedParticles/interface/FindDistCone.h"
0005 #include "Calibration/IsolatedParticles/interface/MatrixHCALDetIds.h"
0006 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008
0009 #include <algorithm>
0010 #include <iostream>
0011 #include <sstream>
0012
0013 namespace spr {
0014
0015
0016 template <typename T>
0017 double eHCALmatrix(const HcalTopology* topology,
0018 const DetId& det0,
0019 edm::Handle<T>& hits,
0020 int ieta,
0021 int iphi,
0022 bool includeHO,
0023 bool algoNew,
0024 double hbThr,
0025 double heThr,
0026 double hfThr,
0027 double hoThr,
0028 double tMin,
0029 double tMax,
0030 int useRaw,
0031 bool debug) {
0032 std::vector<typename T::const_iterator> hit;
0033 HcalDetId hcid0(det0.rawId());
0034 HcalDetId hcid(hcid0.subdet(), hcid0.ieta(), hcid0.iphi(), 1);
0035 DetId det(hcid.rawId());
0036 if (debug)
0037 edm::LogVerbatim("IsoTrack") << "Inside eHCALmatrix " << 2 * ieta + 1 << "X" << 2 * iphi + 1 << " AlgoNew "
0038 << algoNew << " Inclusion of HO Flag " << includeHO;
0039 double energySum = 0.0;
0040 if (algoNew) {
0041 energySum = spr::energyHCALmatrixNew(
0042 topology, det, hits, ieta, iphi, includeHO, hbThr, heThr, hfThr, hoThr, tMin, tMax, useRaw, debug);
0043 } else {
0044 edm::LogVerbatim("IsoTrack") << "eHCALmatrix:: Algorithm New = " << algoNew
0045 << " not supported - call directly spr::energyHCALmatrix ";
0046
0047 }
0048
0049 if (debug)
0050 edm::LogVerbatim("IsoTrack") << "eHCALmatrix::Total energy " << energySum;
0051 return energySum;
0052 }
0053
0054 template <typename T>
0055 double eHCALmatrix(const HcalTopology* topology,
0056 const DetId& det0,
0057 edm::Handle<T>& hits,
0058 int ietaE,
0059 int ietaW,
0060 int iphiN,
0061 int iphiS,
0062 bool includeHO,
0063 double hbThr,
0064 double heThr,
0065 double hfThr,
0066 double hoThr,
0067 double tMin,
0068 double tMax,
0069 int useRaw,
0070 bool debug) {
0071 HcalDetId hcid0(det0.rawId());
0072 HcalDetId hcid(hcid0.subdet(), hcid0.ieta(), hcid0.iphi(), 1);
0073 DetId det(hcid.rawId());
0074 if (debug)
0075 edm::LogVerbatim("IsoTrack") << "Inside eHCALmatrix " << ietaE + ietaW + 1 << "X" << iphiN + iphiS + 1
0076 << " Inclusion of HO Flag " << includeHO;
0077
0078 return spr::energyHCALmatrixTotal(topology,
0079 det,
0080 hits,
0081 ietaE,
0082 ietaW,
0083 iphiN,
0084 iphiS,
0085 includeHO,
0086 hbThr,
0087 heThr,
0088 hfThr,
0089 hoThr,
0090 tMin,
0091 tMax,
0092 useRaw,
0093 debug);
0094 }
0095
0096
0097 template <typename T>
0098 double eHCALmatrix(const HcalTopology* topology,
0099 const DetId& det0,
0100 edm::Handle<T>& hits,
0101 int ieta,
0102 int iphi,
0103 int& nRecHits,
0104 std::vector<int>& RH_ieta,
0105 std::vector<int>& RH_iphi,
0106 std::vector<double>& RH_ene,
0107 std::set<int>& uniqueIdset,
0108 int useRaw) {
0109 bool debug = false;
0110 bool includeHO = false;
0111 HcalDetId hcid0(det0.rawId());
0112 HcalDetId hcid(hcid0.subdet(), hcid0.ieta(), hcid0.iphi(), 1);
0113 DetId det(hcid.rawId());
0114
0115 std::vector<typename T::const_iterator> hit;
0116 spr::hitHCALmatrix(topology, det, hits, ieta, iphi, hit, includeHO, debug);
0117
0118 nRecHits = -99;
0119 double energySum = 0.0;
0120 for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
0121 energySum += spr::getEnergy(hit[ihit], useRaw);
0122 int eta(-99);
0123 int phi(-99);
0124
0125 spr::getEtaPhi(hit[ihit], eta, phi);
0126
0127 RH_ieta.push_back(hit[ihit]->id().ieta());
0128 RH_iphi.push_back(hit[ihit]->id().iphi());
0129 RH_ene.push_back(hit[ihit]->energy());
0130 int uniqueEtaPhiId = 100 * eta + phi;
0131 uniqueIdset.insert(uniqueEtaPhiId);
0132 }
0133 nRecHits = uniqueIdset.size();
0134 return energySum;
0135 }
0136
0137
0138
0139 template <typename T>
0140 double eHCALmatrix(const CaloGeometry* geo,
0141 const HcalTopology* topology,
0142 const DetId& det0,
0143 edm::Handle<T>& hits,
0144 int ieta,
0145 int iphi,
0146 int& nRecHits,
0147 std::vector<int>& RH_ieta,
0148 std::vector<int>& RH_iphi,
0149 std::vector<double>& RH_ene,
0150 GlobalPoint& gPosHotCell,
0151 int useRaw) {
0152 bool debug = false;
0153 bool includeHO = false;
0154 HcalDetId hcid0(det0.rawId());
0155 HcalDetId hcid(hcid0.subdet(), hcid0.ieta(), hcid0.iphi(), 1);
0156 DetId det(hcid.rawId());
0157
0158 std::vector<typename T::const_iterator> hit;
0159 spr::hitHCALmatrix(topology, det, hits, ieta, iphi, hit, includeHO, debug);
0160 double energyMax = -99.;
0161 int hottestIndex = -99;
0162
0163 nRecHits = -99;
0164 double energySum = 0.0;
0165 std::set<int> uniqueIdset;
0166
0167 for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
0168
0169 double energy = getEnergy(hit.at(ihit), useRaw);
0170 if (energy > energyMax) {
0171 energyMax = energy;
0172 hottestIndex = ihit;
0173 }
0174 energySum += energy;
0175 int eta(-99);
0176 int phi(-99);
0177
0178 spr::getEtaPhi(hit[ihit], eta, phi);
0179
0180 RH_ieta.push_back(hit[ihit]->id().ieta());
0181 RH_iphi.push_back(hit[ihit]->id().iphi());
0182 RH_ene.push_back(hit[ihit]->energy());
0183 int uniqueEtaPhiId = 100 * eta + phi;
0184 uniqueIdset.insert(uniqueEtaPhiId);
0185 }
0186
0187
0188 if (hottestIndex != -99) {
0189 gPosHotCell = spr::getGpos(geo, hit.at(hottestIndex));
0190 }
0191
0192 nRecHits = uniqueIdset.size();
0193 return energySum;
0194 }
0195
0196 template <typename T>
0197 double eHCALmatrix(const CaloGeometry* geo,
0198 const HcalTopology* topology,
0199 const DetId& det0,
0200 edm::Handle<T>& hits,
0201 int ieta,
0202 int iphi,
0203 HcalDetId& hotCell,
0204 bool includeHO,
0205 int useRaw,
0206 bool debug) {
0207 HcalDetId hcid0(det0.rawId());
0208 HcalDetId hcid(hcid0.subdet(), hcid0.ieta(), hcid0.iphi(), 1);
0209 DetId det(hcid.rawId());
0210
0211 std::vector<typename T::const_iterator> hit;
0212 spr::hitHCALmatrix(topology, det, hits, ieta, iphi, hit, includeHO, debug);
0213
0214 double energySum(0);
0215 for (unsigned int ihit = 0; ihit < hit.size(); ihit++)
0216 energySum += getEnergy(hit.at(ihit), useRaw);
0217
0218
0219 hotCell = getHotCell(hit, includeHO, useRaw, debug);
0220 return energySum;
0221 }
0222
0223 template <typename T>
0224 double energyHCALmatrixNew(const HcalTopology* topology,
0225 const DetId& det,
0226 edm::Handle<T>& hits,
0227 int ieta,
0228 int iphi,
0229 bool includeHO,
0230 double hbThr,
0231 double heThr,
0232 double hfThr,
0233 double hoThr,
0234 double tMin,
0235 double tMax,
0236 int useRaw,
0237 bool debug) {
0238 std::vector<DetId> dets(1, det);
0239 std::vector<DetId> vdets = spr::matrixHCALIds(dets, topology, ieta, iphi, includeHO, false);
0240 if (debug) {
0241 edm::LogVerbatim("IsoTrack") << "matrixHCALIds::Total number of cells found is " << vdets.size();
0242 spr::debugHcalDets(0, vdets);
0243 }
0244 return spr::energyHCAL(vdets, hits, hbThr, heThr, hfThr, hoThr, tMin, tMax, useRaw, debug);
0245 }
0246
0247 template <typename T>
0248 double energyHCALmatrixTotal(const HcalTopology* topology,
0249 const DetId& det,
0250 edm::Handle<T>& hits,
0251 int ietaE,
0252 int ietaW,
0253 int iphiN,
0254 int iphiS,
0255 bool includeHO,
0256 double hbThr,
0257 double heThr,
0258 double hfThr,
0259 double hoThr,
0260 double tMin,
0261 double tMax,
0262 int useRaw,
0263 bool debug) {
0264 std::vector<DetId> dets(1, det);
0265 std::vector<DetId> vdets = spr::matrixHCALIds(dets, topology, ietaE, ietaW, iphiN, iphiS, includeHO, debug);
0266 return spr::energyHCAL(vdets, hits, hbThr, heThr, hfThr, hoThr, tMin, tMax, useRaw, debug);
0267 }
0268
0269 template <typename T>
0270 void hitHCALmatrix(const HcalTopology* topology,
0271 const DetId& det,
0272 edm::Handle<T>& hits,
0273 int ieta,
0274 int iphi,
0275 std::vector<typename T::const_iterator>& hitlist,
0276 bool includeHO,
0277 bool debug) {
0278 std::vector<DetId> dets(1, det);
0279 std::vector<DetId> vdets = spr::matrixHCALIds(dets, topology, ieta, iphi, includeHO, debug);
0280 spr::hitsHCAL(vdets, hits, hitlist, debug);
0281 }
0282
0283 template <typename T>
0284 void hitHCALmatrixTotal(const HcalTopology* topology,
0285 const DetId& det,
0286 edm::Handle<T>& hits,
0287 int ietaE,
0288 int ietaW,
0289 int iphiN,
0290 int iphiS,
0291 std::vector<typename T::const_iterator>& hitlist,
0292 bool includeHO,
0293 bool debug) {
0294 std::vector<DetId> dets(1, det);
0295 std::vector<DetId> vdets = spr::matrixHCALIds(dets, topology, ietaE, ietaW, iphiN, iphiS, includeHO, debug);
0296 spr::hitsHCAL(vdets, hits, hitlist, debug);
0297 }
0298
0299 template <typename T>
0300 double energyHCAL(std::vector<DetId>& vdets,
0301 edm::Handle<T>& hits,
0302 double hbThr,
0303 double heThr,
0304 double hfThr,
0305 double hoThr,
0306 double tMin,
0307 double tMax,
0308 int useRaw,
0309 bool debug) {
0310 int khit = 0;
0311 double energySum = 0;
0312 unsigned int maxsize = vdets.size();
0313 for (unsigned int i = 0; i < maxsize; i++) {
0314 std::vector<typename T::const_iterator> hit = spr::findHit(hits, vdets[i]);
0315 double energy = 0;
0316 int subdet = ((HcalDetId)(vdets[i].rawId())).subdet();
0317 double eThr = spr::eHCALThreshold(subdet, hbThr, heThr, hfThr, hoThr);
0318 for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
0319 if (hit[ihit] != hits->end()) {
0320 khit++;
0321 if (debug)
0322 edm::LogVerbatim("IsoTrack") << "energyHCAL:: Hit " << khit << " " << (HcalDetId)vdets[i] << " E "
0323 << hit[ihit]->energy() << " t " << hit[ihit]->time();
0324 if (hit[ihit]->time() > tMin && hit[ihit]->time() < tMax) {
0325 energy += getRawEnergy(hit[ihit], useRaw);
0326 }
0327 }
0328 }
0329 if (energy > eThr)
0330 energySum += energy;
0331 }
0332
0333 if (debug)
0334 edm::LogVerbatim("IsoTrack") << "energyHCAL::Total Energy " << energySum << " from " << khit << " hits";
0335 return energySum;
0336 }
0337
0338 template <typename T>
0339 void energyHCALCell(HcalDetId detID,
0340 edm::Handle<T>& hits,
0341 std::vector<std::pair<double, int> >& energyCell,
0342 int maxDepth,
0343 double hbThr,
0344 double heThr,
0345 double hfThr,
0346 double hoThr,
0347 double tMin,
0348 double tMax,
0349 int useRaw,
0350 int depthHE,
0351 bool debug) {
0352 energyCell.clear();
0353 int subdet = detID.subdet();
0354 double eThr = spr::eHCALThreshold(subdet, hbThr, heThr, hfThr, hoThr);
0355 bool hbhe = (detID.ietaAbs() == 16);
0356 if (debug)
0357 edm::LogVerbatim("IsoTrack") << "energyHCALCell: input ID " << detID << " MaxDepth " << maxDepth
0358 << " Threshold (E) " << eThr << " (T) " << tMin << ":" << tMax;
0359 for (int i = 0; i < maxDepth; i++) {
0360 HcalSubdetector subdet0 = (hbhe) ? ((i + 1 >= depthHE) ? HcalEndcap : HcalBarrel) : detID.subdet();
0361 HcalDetId hcid(subdet0, detID.ieta(), detID.iphi(), i + 1);
0362 DetId det(hcid.rawId());
0363 std::vector<typename T::const_iterator> hit = spr::findHit(hits, det);
0364 double energy(0);
0365 for (unsigned int ihit = 0; ihit < hit.size(); ++ihit) {
0366 if (hit[ihit]->time() > tMin && hit[ihit]->time() < tMax)
0367 energy += getRawEnergy(hit[ihit], useRaw);
0368 if (debug)
0369 edm::LogVerbatim("IsoTrack") << "energyHCALCell:: Hit[" << ihit << "] " << hcid << " E "
0370 << hit[ihit]->energy() << " t " << hit[ihit]->time();
0371 }
0372 if (debug)
0373 edm::LogVerbatim("IsoTrack") << "energyHCALCell:: Cell " << hcid << " E " << energy << " from " << hit.size()
0374 << " threshold " << eThr;
0375 if (energy > eThr && !(hit.empty())) {
0376 energyCell.push_back(std::pair<double, int>(energy, i + 1));
0377 }
0378 }
0379 if (debug) {
0380 edm::LogVerbatim("IsoTrack") << "energyHCALCell:: " << energyCell.size() << " entries from " << maxDepth
0381 << " depths:";
0382 std::ostringstream st1;
0383 for (unsigned int i = 0; i < energyCell.size(); ++i) {
0384 st1 << " [" << i << "] (" << energyCell[i].first << ":" << energyCell[i].second << ")";
0385 }
0386 edm::LogVerbatim("IsoTrack") << st1.str();
0387 }
0388 }
0389
0390 template <typename T>
0391 void hitsHCAL(std::vector<DetId>& vdets,
0392 edm::Handle<T>& hits,
0393 std::vector<typename T::const_iterator>& hitlist,
0394 bool debug) {
0395 unsigned int maxsize = vdets.size();
0396 for (unsigned int i = 0; i < maxsize; i++) {
0397 std::vector<typename T::const_iterator> hit = spr::findHit(hits, vdets[i]);
0398 hitlist.insert(hitlist.end(), hit.begin(), hit.end());
0399 }
0400
0401 if (debug)
0402 edm::LogVerbatim("IsoTrack") << "hitsHCAL::Number of hits " << hitlist.size();
0403 }
0404 }