Calibration/IsolatedParticles/interface/eECALMatrix.icc - master - githubmirror/cmssw

Calibration/IsolatedParticles/interface/eECALMatrix.icc
508 строк19 KB
Shahzad Malik Muzaffar
[ALCA] Changes suggested by new llvm18 clang-format
05 сен 2024, 17:03
05 сен 2024, 17:03dfd2283
#include "Calibration/IsolatedParticles/interface/DebugInfo.h"
#include "Calibration/IsolatedParticles/interface/FindCaloHit.h"
#include "Calibration/IsolatedParticles/interface/FindEtaPhi.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include <sstream>

namespace spr {

  template <typename T>
  double eECALmatrix(const DetId& detId,
                     edm::Handle<T>& hitsEB,
                     edm::Handle<T>& hitsEE,
                     const CaloGeometry* geo,
                     const CaloTopology* caloTopology,
                     int ieta,
                     int iphi,
                     double ebThr,
                     double eeThr,
                     double tMin,
                     double tMax,
                     bool debug) {
    std::vector<DetId> vdets;
    spr::matrixECALIds(detId, ieta, iphi, geo, caloTopology, vdets, false);

    if (debug) {
      edm::LogVerbatim("IsoTrack") << "Inside eECALmatrix " << 2 * ieta + 1 << "X" << 2 * iphi + 1 << " nXtals "
                                   << vdets.size();
      spr::debugEcalDets(0, vdets);
    }

    return spr::energyECAL(vdets, hitsEB, hitsEE, ebThr, eeThr, tMin, tMax, debug);
  }

  template <typename T>
  double eECALmatrix(const DetId& detId,
                     edm::Handle<T>& hitsEB,
                     edm::Handle<T>& hitsEE,
                     const CaloGeometry* geo,
                     const CaloTopology* caloTopology,
                     const EcalTrigTowerConstituentsMap& ttMap,
                     int ieta,
                     int iphi,
                     double ebThr,
                     double eeThr,
                     double tMin,
                     double tMax,
                     bool debug) {
    std::vector<DetId> vdets;
    spr::matrixECALIds(detId, ieta, iphi, geo, caloTopology, vdets, debug);

    if (debug) {
      edm::LogVerbatim("IsoTrack") << "Inside eECALmatrix " << 2 * ieta + 1 << "X" << 2 * iphi + 1 << " nXtals "
                                   << vdets.size();
    }

    return spr::energyECAL(vdets, hitsEB, hitsEE, ttMap, ebThr, eeThr, tMin, tMax, debug);
  }

  template <typename T>
  double eECALmatrix(const DetId& detId,
                     edm::Handle<T>& hitsEB,
                     edm::Handle<T>& hitsEE,
                     const CaloGeometry* geo,
                     const CaloTopology* caloTopology,
                     int ietaE,
                     int ietaW,
                     int iphiN,
                     int iphiS,
                     double ebThr,
                     double eeThr,
                     double tMin,
                     double tMax,
                     bool debug) {
    std::vector<DetId> vdets;
    spr::matrixECALIds(detId, ietaE, ietaW, iphiN, iphiS, geo, caloTopology, vdets, debug);

    if (debug) {
      edm::LogVerbatim("IsoTrack") << "Inside eECALmatrix " << ietaE + ietaW + 1 << "X" << iphiN + iphiS + 1
                                   << " nXtals " << vdets.size();
    }

    return spr::energyECAL(vdets, hitsEB, hitsEE, ebThr, eeThr, tMin, tMax, debug);
  }

  template <typename T>
  void hitECALmatrix(CaloNavigator<DetId>& navigator,
                     edm::Handle<T>& hits,
                     int ieta,
                     int iphi,
                     std::vector<typename T::const_iterator>& hitlist,
                     bool debug) {
    DetId thisDet;

    for (int dx = -ieta; dx < ieta + 1; ++dx) {
      for (int dy = -iphi; dy < iphi + 1; ++dy) {
        // shift the navigator by dx/dy crystals in eta/phi
        thisDet = navigator.offsetBy(dx, dy);

        // Place the navigator back to the original position
        navigator.home();

        if (thisDet != DetId(0)) {
          std::vector<typename T::const_iterator> hit;
          spr::findHit(hits, thisDet, hit, debug);
          std::ostringstream st1;
          if (debug && !hit.empty()) {
            if (thisDet.subdetId() == EcalBarrel) {
              EBDetId id = thisDet;
              st1 << "hitECALmatrix::Cell 0x" << std::hex << thisDet() << std::dec << " " << id;
            } else if (thisDet.subdetId() == EcalEndcap) {
              EEDetId id = thisDet;
              st1 << "hitECALmatrix::Cell 0x" << std::hex << thisDet() << std::dec << " " << id;
            } else {
              st1 << "hitECALMatrix::Cell 0x" << std::hex << thisDet() << std::dec << " Unknown Type";
            }
          }
          for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
            if (hit[ihit] != hits->end()) {
              hitlist.push_back(hit[ihit]);
              if (debug)
                st1 << " hit " << ihit << " " << hit[ihit]->energy();
            }
          }
          if (debug && !hit.empty())
            edm::LogVerbatim("IsoTrack") << st1.str();
        }

      }  // iphi
    }  // ieta
  }

  template <typename T>
  double energyECAL(std::vector<DetId>& vdets,
                    edm::Handle<T>& hitsEB,
                    edm::Handle<T>& hitsEE,
                    double ebThr,
                    double eeThr,
                    double tMin,
                    double tMax,
                    bool debug) {
    double energySum = 0.0;
    for (unsigned int i1 = 0; i1 < vdets.size(); i1++) {
      if (vdets[i1] != DetId(0)) {
        std::vector<typename T::const_iterator> hit;
        if (vdets[i1].subdetId() == EcalBarrel) {
          spr::findHit(hitsEB, vdets[i1], hit, debug);
        } else if (vdets[i1].subdetId() == EcalEndcap) {
          spr::findHit(hitsEE, vdets[i1], hit, debug);
        }
        std::ostringstream st1;
        if (debug)
          st1 << "Crystal 0x" << std::hex << vdets[i1]() << std::dec;
        double ener = 0, ethr = ebThr;
        if (vdets[i1].subdetId() != EcalBarrel)
          ethr = eeThr;
        for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
          double en = 0, tt = 0;
          if (vdets[i1].subdetId() == EcalBarrel) {
            if (hit[ihit] != hitsEB->end()) {
              en = hit[ihit]->energy();
              tt = hit[ihit]->time();
            }
          } else if (vdets[i1].subdetId() == EcalEndcap) {
            if (hit[ihit] != hitsEE->end()) {
              en = hit[ihit]->energy();
              tt = hit[ihit]->time();
            }
          }
          if (debug)
            st1 << " " << ihit << " " << en << " Thr " << ethr;
          if (tt > tMin && tt < tMax)
            ener += en;
        }
        if (debug)
          edm::LogVerbatim("IsoTrack") << st1.str();
        if (ener > ethr)
          energySum += ener;
      }
    }
    if (debug)
      edm::LogVerbatim("IsoTrack") << "energyECAL: energySum = " << energySum;
    return energySum;
  }

  template <typename T>
  double energyECAL(std::vector<DetId>& vdets,
                    edm::Handle<T>& hitsEB,
                    edm::Handle<T>& hitsEE,
                    const EcalTrigTowerConstituentsMap& ttMap,
                    double ebThr,
                    double eeThr,
                    double tMin,
                    double tMax,
                    bool debug) {
    double energySum = 0.0;
    for (unsigned int i1 = 0; i1 < vdets.size(); i1++) {
      if (vdets[i1] != DetId(0)) {
        double eTower = spr::energyECALTower(vdets[i1], hitsEB, hitsEE, ttMap, debug);
        bool ok = true;
        if (vdets[i1].subdetId() == EcalBarrel)
          ok = (eTower > ebThr);
        else if (vdets[i1].subdetId() == EcalEndcap)
          ok = (eTower > eeThr);
        std::ostringstream st1;
        if (debug)
          st1 << "Crystal 0x" << std::hex << vdets[i1]() << std::dec << " Flag " << ok;
        if (ok) {
          std::vector<typename T::const_iterator> hit;
          if (vdets[i1].subdetId() == EcalBarrel) {
            spr::findHit(hitsEB, vdets[i1], hit, debug);
          } else if (vdets[i1].subdetId() == EcalEndcap) {
            spr::findHit(hitsEE, vdets[i1], hit, debug);
          }
          double ener = 0;
          for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
            double en = 0, tt = 0;
            if (vdets[i1].subdetId() == EcalBarrel) {
              if (hit[ihit] != hitsEB->end()) {
                en = hit[ihit]->energy();
                tt = hit[ihit]->time();
              }
            } else if (vdets[i1].subdetId() == EcalEndcap) {
              if (hit[ihit] != hitsEE->end()) {
                en = hit[ihit]->energy();
                tt = hit[ihit]->time();
              }
            }
            if (debug)
              st1 << " " << ihit << " E " << en << " time " << tt;
            if (tt > tMin && tt < tMax)
              ener += en;
          }
          energySum += ener;
        }
        if (debug)
          edm::LogVerbatim("IsoTrack") << st1.str();
      }
    }
    if (debug)
      edm::LogVerbatim("IsoTrack") << "energyECAL: energySum = " << energySum;
    return energySum;
  }

  template <typename T>
  double energyECALTower(const DetId& detId,
                         edm::Handle<T>& hitsEB,
                         edm::Handle<T>& hitsEE,
                         const EcalTrigTowerConstituentsMap& ttMap,
                         bool debug) {
    double ener = 0;
    EcalTrigTowerDetId trId = ttMap.towerOf(detId);
    std::vector<DetId> vdets = ttMap.constituentsOf(trId);
    if (debug) {
      std::ostringstream st1;
      st1 << "energyECALTower: ";
      if (detId.subdetId() == EcalBarrel) {
        EBDetId id = detId;
        st1 << "Cell 0x" << std::hex << detId() << std::dec << " " << id;
      } else if (detId.subdetId() == EcalEndcap) {
        EEDetId id = detId;
        st1 << "Cell 0x" << std::hex << detId() << std::dec << " " << id;
      } else {
        st1 << "Cell 0x" << std::hex << detId() << std::dec << " Unknown Type";
      }
      edm::LogVerbatim("IsoTrack") << st1.str() << " Tower " << trId << " with " << vdets.size() << " cells";
    }
    for (unsigned int i1 = 0; i1 < vdets.size(); i1++) {
      if (vdets[i1] != DetId(0)) {
        std::vector<typename T::const_iterator> hit;
        if (vdets[i1].subdetId() == EcalBarrel) {
          spr::findHit(hitsEB, vdets[i1], hit, debug);
        } else if (vdets[i1].subdetId() == EcalEndcap) {
          spr::findHit(hitsEE, vdets[i1], hit, debug);
        }
        std::ostringstream st1;
        if (debug)
          st1 << "Xtal 0x" << std::hex << vdets[i1]() << std::dec;
        double en = 0;
        for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
          if (vdets[i1].subdetId() == EcalBarrel) {
            if (hit[ihit] != hitsEB->end())
              en += hit[ihit]->energy();
          } else if (vdets[i1].subdetId() == EcalEndcap) {
            if (hit[ihit] != hitsEE->end())
              en += hit[ihit]->energy();
          }
        }
        if (debug)
          edm::LogVerbatim("IsoTrack") << st1.str() << " " << hit.size() << " E " << en;
        ener += en;
      }
    }
    if (debug)
      edm::LogVerbatim("IsoTrack") << "energyECALTower: Energy in the Tower = " << ener;
    return ener;
  }

  template <typename T>
  DetId hotCrystal(const DetId& detId,
                   edm::Handle<T>& hitsEB,
                   edm::Handle<T>& hitsEE,
                   const CaloGeometry* geo,
                   const CaloTopology* caloTopology,
                   int ieta,
                   int iphi,
                   double tMin,
                   double tMax,
                   bool debug) {
    std::vector<DetId> vdets;
    spr::matrixECALIds(detId, ieta, iphi, geo, caloTopology, vdets, false);

    if (debug) {
      edm::LogVerbatim("IsoTrack") << "Inside eECALmatrix " << 2 * ieta + 1 << "X" << 2 * iphi + 1 << " nXtals "
                                   << vdets.size();
      spr::debugEcalDets(0, vdets);
    }

    DetId det = spr::hotCrystal(vdets, hitsEB, hitsEE, tMin, tMax, debug);
    if (det == DetId(0))
      det = detId;
    return det;
  }

  template <typename T>
  DetId hotCrystal(
      std::vector<DetId>& vdets, edm::Handle<T>& hitsEB, edm::Handle<T>& hitsEE, double tMin, double tMax, bool debug) {
    DetId det = DetId(0);
    double eMax = -99999.;
    for (unsigned int i1 = 0; i1 < vdets.size(); i1++) {
      if (vdets[i1] != DetId(0)) {
        std::vector<typename T::const_iterator> hit;
        if (vdets[i1].subdetId() == EcalBarrel) {
          spr::findHit(hitsEB, vdets[i1], hit, debug);
        } else if (vdets[i1].subdetId() == EcalEndcap) {
          spr::findHit(hitsEE, vdets[i1], hit, debug);
        }
        std::ostringstream st1;
        if (debug)
          st1 << "Crystal 0x" << std::hex << vdets[i1]() << std::dec;
        double ener = 0;
        for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
          double en = 0, tt = 0;
          if (vdets[i1].subdetId() == EcalBarrel) {
            if (hit[ihit] != hitsEB->end()) {
              en = hit[ihit]->energy();
              tt = hit[ihit]->time();
            }
          } else if (vdets[i1].subdetId() == EcalEndcap) {
            if (hit[ihit] != hitsEE->end()) {
              en = hit[ihit]->energy();
              tt = hit[ihit]->time();
            }
          }
          if (debug)
            st1 << " " << ihit << " " << en << " " << tt;
          if (tt > tMin && tt < tMax)
            ener += en;
        }
        if (debug)
          edm::LogVerbatim("IsoTrack") << st1.str();
        if (ener > eMax) {
          det = vdets[i1];
          eMax = ener;
        }
      }
    }
    if (debug)
      edm::LogVerbatim("IsoTrack") << "energyECAL: maxEnegy = " << eMax;
    return det;
  }

  template <typename T>
  double eECALmatrix(CaloNavigator<DetId>& navigator,
                     edm::Handle<T>& hits,
                     int ieta,
                     int iphi,
                     const EcalSeverityLevelAlgo* sevlv,
                     bool debug) {
    std::vector<typename T::const_iterator> hit;
    spr::hitECALmatrix(navigator, hits, ieta, iphi, hit, debug);

    if (debug) {
      edm::LogVerbatim("IsoTrack") << "Inside eECALmatrix " << 2 * ieta + 1 << "X" << 2 * iphi + 1;
      std::ostringstream st1;
      st1 << "nXtals " << hit.size();
      for (unsigned int ihit = 0; ihit < hit.size(); ihit++)
        st1 << " ihit:" << ihit << " " << (unsigned int)hit[ihit]->id();
      edm::LogVerbatim("IsoTrack") << st1.str();
    }

    double energySum = 0.0;
    for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
      energySum += hit[ihit]->energy();
    }
    return energySum;
  }

  template <typename T>
  std::vector<std::pair<DetId, double> > eECALmatrixCell(const DetId& detId,
                                                         edm::Handle<T>& hitsEB,
                                                         edm::Handle<T>& hitsEE,
                                                         const CaloGeometry* geo,
                                                         const CaloTopology* caloTopology,
                                                         int ieta,
                                                         int iphi,
                                                         double ebThr,
                                                         double eeThr,
                                                         bool debug) {
    std::vector<DetId> vdets = spr::matrixECALIds(detId, ieta, iphi, geo, caloTopology, debug);

    if (debug) {
      edm::LogVerbatim("IsoTrack") << "Inside eECALmatrixCell " << 2 * ieta + 1 << "X" << 2 * iphi + 1 << " nXtals "
                                   << vdets.size();
    }

    return spr::energyECALCell(vdets, hitsEB, hitsEE, ebThr, eeThr, debug);
  }

  template <typename T>
  std::pair<double, int> eECALmatrixTotal(const DetId& detId,
                                          edm::Handle<T>& hitsEB,
                                          edm::Handle<T>& hitsEE,
                                          const CaloGeometry* geo,
                                          const CaloTopology* caloTopology,
                                          int ieta,
                                          int iphi,
                                          double ebThr,
                                          double eeThr,
                                          double tMin,
                                          double tMax,
                                          bool debug) {
    spr::EtaPhi etaphi = spr::getEtaPhi(ieta, iphi, debug);

    // Get maximum of all the trials
    double energySum = 0;
    int itrym = 0;
    for (int itry = 0; itry < etaphi.ntrys; itry++) {
      std::vector<DetId> vdets = spr::matrixECALIds(detId,
                                                    etaphi.ietaE[itry],
                                                    etaphi.ietaW[itry],
                                                    etaphi.iphiN[itry],
                                                    etaphi.iphiS[itry],
                                                    geo,
                                                    caloTopology,
                                                    debug);
      double energy = spr::energyECAL(vdets, hitsEB, hitsEE, ebThr, eeThr, tMin, tMax, debug);
      if (energy > energySum) {
        energySum = energy;
        itrym = itry;
      }
    }

    if (debug)
      edm::LogVerbatim("IsoTrack") << "eECALmatrixTotal:: energy deposit in " << ieta << "X" << iphi << " matrix is "
                                   << energySum << " for trial # " << itrym;
    return std::pair<double, int>(energySum, itrym);
  }

  template <typename T>
  std::vector<std::pair<DetId, double> > energyECALCell(std::vector<DetId>& vdets,
                                                        edm::Handle<T>& hitsEB,
                                                        edm::Handle<T>& hitsEE,
                                                        double ebThr,
                                                        double eeThr,
                                                        bool debug) {
    std::vector<std::pair<DetId, double> > energyDets;
    for (unsigned int i1 = 0; i1 < vdets.size(); i1++) {
      double energySum = 0.0;
      if (vdets[i1] != DetId(0)) {
        std::vector<typename T::const_iterator> hit;
        if (vdets[i1].subdetId() == EcalBarrel) {
          hit = spr::findHit(hitsEB, vdets[i1]);
        } else if (vdets[i1].subdetId() == EcalEndcap) {
          hit = spr::findHit(hitsEE, vdets[i1]);
        }
        std::ostringstream st1;
        if (debug)
          st1 << "Xtal 0x" << std::hex << vdets[i1]() << std::dec;
        bool ok = false;
        double ethr = ebThr;
        for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
          double en = 0;
          if (vdets[i1].subdetId() == EcalBarrel) {
            if (hit[ihit] != hitsEB->end()) {
              en = hit[ihit]->energy();
              ok = true;
            }
          } else if (vdets[i1].subdetId() == EcalEndcap) {
            if (hit[ihit] != hitsEE->end()) {
              en = hit[ihit]->energy();
              ok = true;
              ethr = eeThr;
            }
          }
          if (debug)
            st1 << " " << ihit << " " << en;
          energySum += en;
        }
        if (debug)
          edm::LogVerbatim("IsoTrack") << st1.str() << "\nenergyECALCell: energySum = " << energySum;
        if (ok && energySum > ethr)
          energyDets.push_back(std::pair<DetId, double>(vdets[i1], energySum));
      }
    }
    return energyDets;
  }

}  // namespace spr