eHCALMatrix.cc

Go to the documentation of this file.

#include "Calibration/IsolatedParticles/interface/DebugInfo.h"
#include "Calibration/IsolatedParticles/interface/eHCALMatrix.h"
#include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
 
#include <algorithm>
#include <iostream>
 
//#define EDM_ML_DEBUG
 
namespace spr {
  double eHCALmatrix(const HcalTopology* topology,
                     const DetId& det0,
                     std::vector<PCaloHit>& hits,
                     int ieta,
                     int iphi,
                     bool includeHO,
                     double hbThr,
                     double heThr,
                     double hfThr,
                     double hoThr,
                     double tMin,
                     double tMax,
                     bool
#ifdef EDM_ML_DEBUG
                         debug
#endif
  ) {
 
    HcalDetId hcid0(det0.rawId());
    HcalDetId hcid(hcid0.subdet(), hcid0.ieta(), hcid0.iphi(), 1);
    DetId det(hcid.rawId());
#ifdef EDM_ML_DEBUG
    if (debug)
      std::cout << "Inside eHCALmatrix " << 2 * ieta + 1 << "X" << 2 * iphi + 1 << " Inclusion of HO Flag " << includeHO
                << std::endl;
#endif
    double energySum(0);
    std::vector<DetId> dets(1, det);
    std::vector<DetId> vdets = spr::matrixHCALIds(dets, topology, ieta, iphi, includeHO, false);
#ifdef EDM_ML_DEBUG
    if (debug) {
      std::cout << "matrixHCALIds::Total number of cells found is " << vdets.size() << std::endl;
      spr::debugHcalDets(0, vdets);
    }
#endif
    int khit(0);
    for (unsigned int i = 0; i < vdets.size(); i++) {
      std::vector<std::vector<PCaloHit>::const_iterator> hit = spr::findHit(hits, vdets[i]);
      double energy = 0;
      int subdet = ((HcalDetId)(vdets[i].rawId())).subdet();
      double eThr = spr::eHCALThreshold(subdet, hbThr, heThr, hfThr, hoThr);
      for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
        if (hit[ihit] != hits.end()) {
          khit++;
#ifdef EDM_ML_DEBUG
          if (debug)
            std::cout << "energyHCAL:: Hit " << khit << " " << (HcalDetId)vdets[i] << " E " << hit[ihit]->energy()
                      << " t " << hit[ihit]->time() << std::endl;
#endif
          if (hit[ihit]->time() > tMin && hit[ihit]->time() < tMax) {
            energy += hit[ihit]->energy();
          }
        }
      }
      if (energy > eThr)
        energySum += energy;
    }
 
#ifdef EDM_ML_DEBUG
    if (debug)
      std::cout << "eHCALmatrix::Total energy " << energySum << std::endl;
#endif
    return energySum;
  }
 
  double eHCALmatrix(const CaloGeometry* geo,
                     const HcalTopology* topology,
                     const DetId& det0,
                     std::vector<PCaloHit>& hits,
                     int ieta,
                     int iphi,
                     HcalDetId& hotCell,
                     bool includeHO,
                     bool debug) {
    HcalDetId hcid0(det0.rawId());
    HcalDetId hcid(hcid0.subdet(), hcid0.ieta(), hcid0.iphi(), 1);
    DetId det(hcid.rawId());
    std::vector<DetId> dets(1, det);
    std::vector<DetId> vdets = spr::matrixHCALIds(dets, topology, ieta, iphi, includeHO, debug);
    hotCell = hcid0;
 
    std::vector<std::vector<PCaloHit>::const_iterator> hitlist;
    for (unsigned int i = 0; i < vdets.size(); i++) {
      std::vector<std::vector<PCaloHit>::const_iterator> hit = spr::findHit(hits, vdets[i]);
      hitlist.insert(hitlist.end(), hit.begin(), hit.end());
    }
 
    double energySum(0);
    for (unsigned int ihit = 0; ihit < hitlist.size(); ihit++)
      energySum += hitlist[ihit]->energy();
 
    // Get hotCell ID
    dets.clear();
    std::vector<double> energies;
    for (unsigned int ihit = 0; ihit < hitlist.size(); ihit++) {
      double energy = hitlist[ihit]->energy();
      HcalDetId id0 = HcalDetId(hitlist[ihit]->id());
      if ((id0.subdet() != HcalOuter) || includeHO) {
        HcalDetId id1(id0.subdet(), id0.ieta(), id0.iphi(), 1);
        bool found(false);
        for (unsigned int idet = 0; idet < dets.size(); ++idet) {
          if (id1 == HcalDetId(dets[idet])) {
            energies[idet] += energy;
            found = true;
            break;
          }
        }
        if (!found) {
          dets.push_back(DetId(id1));
          energies.push_back(energy);
        }
      }
    }
    double energyMax(-99.);
    for (unsigned int ihit = 0; ihit < dets.size(); ihit++) {
      if (energies[ihit] > energyMax) {
        energyMax = energies[ihit];
        hotCell = HcalDetId(dets[ihit]);
      }
    }
    return energySum;
  }
 
  void energyHCALCell(HcalDetId detId,
                      std::vector<PCaloHit>& hits,
                      std::vector<std::pair<double, int> >& energyCell,
                      int maxDepth,
                      double hbThr,
                      double heThr,
                      double hfThr,
                      double hoThr,
                      double tMin,
                      double tMax,
                      int depthHE,
                      bool
#ifdef EDM_ML_DEBUG
                          debug
#endif
  ) {
 
    energyCell.clear();
    int subdet = detId.subdet();
    double eThr = spr::eHCALThreshold(subdet, hbThr, heThr, hfThr, hoThr);
    bool hbhe = (detId.ietaAbs() == 16);
#ifdef EDM_ML_DEBUG
    if (debug)
      std::cout << "energyHCALCell: input ID " << detId << " MaxDepth " << maxDepth << " Threshold (E) " << eThr
                << " (T) " << tMin << ":" << tMax << std::endl;
#endif
    for (int i = 0; i < maxDepth; i++) {
      HcalSubdetector subdet0 = (hbhe) ? ((i + 1 >= depthHE) ? HcalEndcap : HcalBarrel) : detId.subdet();
      HcalDetId hcid(subdet0, detId.ieta(), detId.iphi(), i + 1);
      DetId det(hcid.rawId());
      std::vector<std::vector<PCaloHit>::const_iterator> hit = spr::findHit(hits, det);
      double energy(0);
      for (unsigned int ihit = 0; ihit < hit.size(); ++ihit) {
        if (hit[ihit]->time() > tMin && hit[ihit]->time() < tMax)
          energy += hit[ihit]->energy();
#ifdef EDM_ML_DEBUG
        if (debug)
          std::cout << "energyHCALCell:: Hit[" << ihit << "] " << hcid << " E " << hit[ihit]->energy() << " t "
                    << hit[ihit]->time() << std::endl;
#endif
      }
#ifdef EDM_ML_DEBUG
      if (debug)
        std::cout << "energyHCALCell:: Cell " << hcid << " E " << energy << " from " << hit.size() << " threshold "
                  << eThr << std::endl;
#endif
      if (energy > eThr && !hit.empty()) {
        energyCell.push_back(std::pair<double, int>(energy, i + 1));
      }
    }
#ifdef EDM_ML_DEBUG
    if (debug) {
      std::cout << "energyHCALCell:: " << energyCell.size() << " entries from " << maxDepth << " depths:";
      for (unsigned int i = 0; i < energyCell.size(); ++i) {
        std::cout << " [" << i << "] (" << energyCell[i].first << ":" << energyCell[i].second << ")";
      }
      std::cout << std::endl;
    }
#endif
  }
 
  HcalDetId getHotCell(std::vector<HBHERecHitCollection::const_iterator>& hit, bool includeHO, int useRaw, bool) {
    std::vector<HcalDetId> dets;
    std::vector<double> energies;
    for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
      double energy = getRawEnergy(hit.at(ihit), useRaw);
      HcalDetId id0 = hit.at(ihit)->id();
      if ((id0.subdet() != HcalOuter) || includeHO) {
        HcalDetId id1(id0.subdet(), id0.ieta(), id0.iphi(), 1);
        bool found(false);
        for (unsigned int idet = 0; idet < dets.size(); ++idet) {
          if (id1 == dets[idet]) {
            energies[idet] += energy;
            found = true;
            break;
          }
        }
        if (!found) {
          dets.push_back(id1);
          energies.push_back(energy);
        }
      }
    }
    double energyMax(-99.);
    HcalDetId hotCell;
    for (unsigned int ihit = 0; ihit < dets.size(); ihit++) {
      if (energies[ihit] > energyMax) {
        energyMax = energies[ihit];
        hotCell = dets[ihit];
      }
    }
    return hotCell;
  }
 
  HcalDetId getHotCell(std::vector<std::vector<PCaloHit>::const_iterator>& hit, bool includeHO, int useRaw, bool) {
    std::vector<HcalDetId> dets;
    std::vector<double> energies;
    for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
      double energy = hit.at(ihit)->energy();
      HcalDetId id0 = getRawEnergy(hit.at(ihit), useRaw);
      if ((id0.subdet() != HcalOuter) || includeHO) {
        HcalDetId id1(id0.subdet(), id0.ieta(), id0.iphi(), 1);
        bool found(false);
        for (unsigned int idet = 0; idet < dets.size(); ++idet) {
          if (id1 == dets[idet]) {
            energies[idet] += energy;
            found = true;
            break;
          }
        }
        if (!found) {
          dets.push_back(id1);
          energies.push_back(energy);
        }
      }
    }
    double energyMax(-99.);
    HcalDetId hotCell;
    for (unsigned int ihit = 0; ihit < dets.size(); ihit++) {
      if (energies[ihit] > energyMax) {
        energyMax = energies[ihit];
        hotCell = dets[ihit];
      }
    }
    return hotCell;
  }
 
  double eHCALThreshold(int subdet, double hbThr, double heThr, double hfThr, double hoThr) {
    double eThr = hbThr;
    if (subdet == (int)(HcalEndcap))
      eThr = heThr;
    else if (subdet == (int)(HcalForward))
      eThr = hfThr;
    else if (subdet == (int)(HcalOuter))
      eThr = hoThr;
    return eThr;
  }
}  // namespace spr