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

Calibration/IsolatedParticles/interface/eHCALMatrixExtra.icc
443 строки18 KB
Sunanda
Make use of MessageLogger wherever possible in Calibration/IsolatedParticles and HcalIsolatedTrackReco packages
21 июл 2022, 04:43
21 июл 2022, 04:43fe1c621
#include "Calibration/IsolatedParticles/interface/FindCaloHit.h"
#include "Calibration/IsolatedParticles/interface/FindEtaPhi.h"
#include "Calibration/IsolatedParticles/interface/MatrixHCALDetIds.h"
#include "DataFormats/HcalDetId/interface/HcalSubdetector.h"

#include <algorithm>
#include <iostream>

namespace spr {

  template <typename T>
  std::vector<std::pair<DetId, double> > eHCALmatrixCell(const HcalTopology* topology,
                                                         const DetId& det,
                                                         edm::Handle<T>& hits,
                                                         int ieta,
                                                         int iphi,
                                                         bool includeHO,
                                                         double hbThr,
                                                         double heThr,
                                                         double hfThr,
                                                         double hoThr,
                                                         bool debug) {
    std::vector<DetId> dets(1, det);
    std::vector<DetId> vdets = spr::matrixHCALIds(dets, topology, ieta, iphi, includeHO, debug);
    return spr::energyDetIdHCAL(vdets, hits, hbThr, heThr, hfThr, hoThr, debug);
  }

  //*** Returns energy in NxN Hcal towers
  template <typename T>
  std::pair<double, int> eHCALmatrixTotal(const HcalTopology* topology,
                                          const DetId& det0,
                                          edm::Handle<T>& hits,
                                          int ieta,
                                          int iphi,
                                          bool includeHO,
                                          double hbThr,
                                          double heThr,
                                          double hfThr,
                                          double hoThr,
                                          bool debug) {
    HcalDetId hcid0(det0.rawId());
    HcalDetId hcid(hcid0.subdet(), hcid0.ieta(), hcid0.iphi(), 1);
    DetId det(hcid.rawId());
    if (debug)
      edm::LogVerbatim("IsoTrack") << "Inside eHCALmatrixTotal " << ieta << "X" << iphi << " Inclusion of HO Flag "
                                   << includeHO;

    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++) {
      double energy = energyHCALmatrixTotal(topology,
                                            det,
                                            hits,
                                            etaphi.ietaE[itry],
                                            etaphi.ietaW[itry],
                                            etaphi.iphiN[itry],
                                            etaphi.iphiS[itry],
                                            includeHO,
                                            hbThr,
                                            heThr,
                                            hfThr,
                                            hoThr,
                                            debug);
      if (energy > energySum) {
        energySum = energy;
        itrym = itry;
      }
      if (debug)
        edm::LogVerbatim("IsoTrack") << "eHCALMatrix::Total energy " << energy << " Max " << energySum << "in trial # "
                                     << itrym;
    }

    return std::pair<double, int>(energySum, itrym);
  }

  //*** Returns vector of hits in 3x3 or 5x5 Hcal towers
  template <typename T>
  double energyHCALmatrix(const HcalTopology* topology,
                          const DetId& det0,
                          edm::Handle<T>& hits,
                          int ieta,
                          int iphi,
                          bool includeHO,
                          double hbThr,
                          double heThr,
                          double hfThr,
                          double hoThr,
                          bool debug) {
    HcalDetId hcid0(det0.rawId());
    HcalDetId hcid(hcid0.subdet(), hcid0.ieta(), hcid0.iphi(), 1);
    DetId det(hcid.rawId());
    if (debug)
      edm::LogVerbatim("IsoTrack") << "Inside energyHCALmatrix " << 2 * ieta + 1 << "X" << 2 * iphi + 1
                                   << " Inclusion of HO Flag " << includeHO;

    double energy = 0;
    if (ieta > 2) {
      edm::LogVerbatim("IsoTrack") << " no matrix > 5x5 is coded ! ";
      return energy;
    }
    // used dets during building the matrix
    std::vector<DetId> dets;
    dets.clear();

    // central tower
    std::vector<DetId> vNeighboursDetId(1, det);
    energy += energyHCAL(vNeighboursDetId, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug);

    if (ieta > 0) {
      // go to east from central tower
      vNeighboursDetId.clear();
      if (debug)
        edm::LogVerbatim("IsoTrack") << "hitHCALMatrix:: det " << (HcalDetId)det;
      vNeighboursDetId = topology->east(det);
      if (debug) {
        std::ostringstream st1;
        st1 << "Neighbour:: East";
        for (unsigned int i = 0; i < vNeighboursDetId.size(); i++)
          st1 << " " << (HcalDetId)vNeighboursDetId[i];
        edm::LogVerbatim("IsoTrack") << st1.str();
      }
      energy += energyHCAL(vNeighboursDetId, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug);
      if (ieta == 2) {
        for (int ii = 0; ii < (int)vNeighboursDetId.size(); ii++) {
          std::vector<DetId> vNeighboursDetIdc = topology->east(vNeighboursDetId[ii]);
          if (debug) {
            std::ostringstream st1;
            st1 << "Neighbour:: East";
            for (unsigned int i = 0; i < vNeighboursDetIdc.size(); i++)
              st1 << " " << (HcalDetId)vNeighboursDetIdc[i];
            edm::LogVerbatim("IsoTrack") << st1.str();
          }
          energy += energyHCAL(vNeighboursDetIdc, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug);
        }
      }

      // go to west from central tower
      vNeighboursDetId.clear();
      vNeighboursDetId = topology->west(det);
      if (debug) {
        std::ostringstream st1;
        st1 << "Neighbour:: West";
        for (unsigned int i = 0; i < vNeighboursDetId.size(); i++)
          st1 << " " << (HcalDetId)vNeighboursDetId[i];
        edm::LogVerbatim("IsoTrack") << st1.str();
      }
      energy += energyHCAL(vNeighboursDetId, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug);
      if (ieta == 2) {
        for (int ii = 0; ii < (int)vNeighboursDetId.size(); ii++) {
          std::vector<DetId> vNeighboursDetIdc = topology->west(vNeighboursDetId[ii]);
          if (debug) {
            std::ostringstream st1;
            st1 << "Neighbour:: West";
            for (unsigned int i = 0; i < vNeighboursDetIdc.size(); i++)
              st1 << " " << (HcalDetId)vNeighboursDetIdc[i];
            edm::LogVerbatim("IsoTrack") << st1.str();
          }
          energy += energyHCAL(vNeighboursDetIdc, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug);
        }
      }
    }

    // do steps in phi to north
    DetId detnorth = det;
    for (int inorth = 0; inorth < iphi; inorth++) {
      if (debug)
        edm::LogVerbatim("IsoTrack") << "hitHCALMatrix:: detNorth " << (HcalDetId)detnorth;
      std::vector<DetId> NorthDetId = topology->north(detnorth);
      if (debug) {
        std::ostringstream st1;
        st1 << "Neighbour:: North";
        for (unsigned int i = 0; i < NorthDetId.size(); i++)
          st1 << " " << (HcalDetId)NorthDetId[i];
        edm::LogVerbatim("IsoTrack") << st1.str();
      }
      if (!(NorthDetId.empty())) {
        energy += energyHCAL(NorthDetId, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug);
        if (ieta > 0) {
          // go to east
          vNeighboursDetId.clear();
          vNeighboursDetId = topology->east(NorthDetId[0]);
          if (debug) {
            std::ostringstream st1;
            st1 << "Neighbour:: East";
            for (unsigned int i = 0; i < vNeighboursDetId.size(); i++)
              st1 << " " << (HcalDetId)vNeighboursDetId[i];
            edm::LogVerbatim("IsoTrack") << st1.str();
          }
          energy += energyHCAL(vNeighboursDetId, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug);
          if (ieta == 2) {
            for (int ii = 0; ii < (int)vNeighboursDetId.size(); ii++) {
              std::vector<DetId> vNeighboursDetIdc = topology->east(vNeighboursDetId[ii]);
              if (debug) {
                std::ostringstream st1;
                st1 << "Neighbour:: East";
                for (unsigned int i = 0; i < vNeighboursDetIdc.size(); i++)
                  st1 << " " << (HcalDetId)vNeighboursDetIdc[i];
                edm::LogVerbatim("IsoTrack") << st1.str();
              }
              energy +=
                  energyHCAL(vNeighboursDetIdc, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug);
            }
          }

          // go to west
          vNeighboursDetId.clear();
          vNeighboursDetId = topology->west(NorthDetId[0]);
          if (debug) {
            std::ostringstream st1;
            st1 << "Neighbour:: West";
            for (unsigned int i = 0; i < vNeighboursDetId.size(); i++)
              st1 << " " << (HcalDetId)vNeighboursDetId[i];
            edm::LogVerbatim("IsoTrack") << st1.str();
          }
          energy += energyHCAL(vNeighboursDetId, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug);
          if (ieta == 2) {
            for (int ii = 0; ii < (int)vNeighboursDetId.size(); ii++) {
              std::vector<DetId> vNeighboursDetIdc = topology->west(vNeighboursDetId[ii]);
              if (debug) {
                std::ostringstream st1;
                st1 << "Neighbour:: West";
                for (unsigned int i = 0; i < vNeighboursDetIdc.size(); i++)
                  st1 << " " << (HcalDetId)vNeighboursDetIdc[i];
                edm::LogVerbatim("IsoTrack") << st1.str();
              }
              energy +=
                  energyHCAL(vNeighboursDetIdc, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug);
            }
          }
        }
        detnorth = NorthDetId[0];
      } else {
        break;
      }
    }

    // do steps in phi to south
    DetId detsouth = det;
    for (int isouth = 0; isouth < iphi; isouth++) {
      if (debug)
        edm::LogVerbatim("IsoTrack") << "hitHCALMatrix:: detSouth " << (HcalDetId)detsouth;
      vNeighboursDetId.clear();
      std::vector<DetId> SouthDetId = topology->south(detsouth);
      if (!(SouthDetId.empty())) {
        energy += energyHCAL(SouthDetId, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug);
        if (debug) {
          std::ostringstream st1;
          st1 << "Neighbour:: South";
          for (unsigned int i = 0; i < SouthDetId.size(); i++)
            st1 << " " << (HcalDetId)SouthDetId[i];
          edm::LogVerbatim("IsoTrack") << st1.str();
        }

        // go to east
        if (ieta > 0) {
          // go to east from central tower
          vNeighboursDetId.clear();
          vNeighboursDetId = topology->east(SouthDetId[0]);
          if (debug) {
            std::ostringstream st1;
            st1 << "Neighbour:: East";
            for (unsigned int i = 0; i < vNeighboursDetId.size(); i++)
              st1 << " " << (HcalDetId)vNeighboursDetId[i];
            edm::LogVerbatim("IsoTrack") << st1.str();
          }
          energy += energyHCAL(vNeighboursDetId, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug);
          if (ieta == 2) {
            for (int ii = 0; ii < (int)vNeighboursDetId.size(); ii++) {
              std::vector<DetId> vNeighboursDetIdc = topology->east(vNeighboursDetId[ii]);
              if (debug) {
                std::ostringstream st1;
                st1 << "Neighbour:: East";
                for (unsigned int i = 0; i < vNeighboursDetIdc.size(); i++)
                  st1 << " " << (HcalDetId)vNeighboursDetIdc[i];
                edm::LogVerbatim("IsoTrack") << st1.str();
              }
              energy +=
                  energyHCAL(vNeighboursDetIdc, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug);
            }
          }

          // go to west
          vNeighboursDetId.clear();
          vNeighboursDetId = topology->west(SouthDetId[0]);
          if (debug) {
            std::ostringstream st1;
            st1 << "Neighbour:: West";
            for (unsigned int i = 0; i < vNeighboursDetId.size(); i++)
              st1 << " " << (HcalDetId)vNeighboursDetId[i];
            edm::LogVerbatim("IsoTrack") << st1.str();
          }
          energy += energyHCAL(vNeighboursDetId, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug);
          if (ieta == 2) {
            for (int ii = 0; ii < (int)vNeighboursDetId.size(); ii++) {
              std::vector<DetId> vNeighboursDetIdc = topology->west(vNeighboursDetId[ii]);
              if (debug) {
                std::ostringstream st1;
                st1 << "Neighbour:: West";
                for (unsigned int i = 0; i < vNeighboursDetIdc.size(); i++)
                  st1 << " " << (HcalDetId)vNeighboursDetIdc[i];
                edm::LogVerbatim("IsoTrack") << st1.str();
              }
              energy +=
                  energyHCAL(vNeighboursDetIdc, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug);
            }
          }
        }
        detsouth = SouthDetId[0];
      } else {
        break;
      }
    }

    return energy;
  }

  template <typename T>
  double energyHCAL(std::vector<DetId>& vNeighboursDetId,
                    std::vector<DetId>& dets,
                    const HcalTopology* topology,
                    edm::Handle<T>& hits,
                    bool includeHO,
                    double hbThr,
                    double heThr,
                    double hfThr,
                    double hoThr,
                    bool debug) {
    double energySum = 0;
    int kHit = 0;
    for (int i = 0; i < (int)vNeighboursDetId.size(); i++) {
      int n = std::count(dets.begin(), dets.end(), vNeighboursDetId[i]);
      if (n != 0)
        continue;
      dets.push_back(vNeighboursDetId[i]);
      std::vector<typename T::const_iterator> hit = spr::findHit(hits, vNeighboursDetId[i]);

      double energy = 0;
      int subdet = ((HcalDetId)(vNeighboursDetId[i].rawId())).subdet();
      double eThr = hbThr;
      if (subdet == (int)(HcalEndcap))
        eThr = heThr;
      else if (subdet == (int)(HcalForward))
        eThr = hfThr;
      else if (subdet == (int)(HcalOuter))
        eThr = hoThr;
      for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
        if (hit[ihit] != hits->end()) {
          energy += hit[ihit]->energy();
          kHit++;
        }
      }
      if (energy > eThr)
        energySum += energy;

      HcalDetId depth = vNeighboursDetId[i];
      // max. three depths can be in endcap and we go to 2nd and 3rd from 1st where we are now
      for (int idepth = 0; idepth < 2; idepth++) {
        std::vector<DetId> vUpDetId = topology->up(depth);
        if (!(vUpDetId.empty())) {
          if (includeHO || vUpDetId[0].subdetId() != (int)(HcalOuter)) {
            int n = std::count(dets.begin(), dets.end(), vUpDetId[0]);
            if (n == 0) {
              dets.push_back(vUpDetId[0]);
              hit = spr::findHit(hits, vUpDetId[0]);
              energy = 0;
              subdet = ((HcalDetId)(vUpDetId[0].rawId())).subdet();
              eThr = hbThr;
              if (subdet == (int)(HcalEndcap))
                eThr = heThr;
              else if (subdet == (int)(HcalForward))
                eThr = hfThr;
              else if (subdet == (int)(HcalOuter))
                eThr = hoThr;
              for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
                if (hit[ihit] != hits->end()) {
                  energy += hit[ihit]->energy();
                  kHit++;
                }
              }
              if (energy > eThr)
                energySum += energy;
            }
          }
          depth = vUpDetId[0];
        }
      }
    }
    if (debug)
      edm::LogVerbatim("IsoTrack") << "energyHCAL:: Energy " << energySum << " from " << kHit << " hits";
    return energySum;
  }

  template <typename T>
  std::vector<std::pair<DetId, double> > energyDetIdHCAL(std::vector<DetId>& vdets,
                                                         edm::Handle<T>& hits,
                                                         double hbThr,
                                                         double heThr,
                                                         double hfThr,
                                                         double hoThr,
                                                         bool debug) {
    std::vector<std::pair<DetId, double> > energyDetId;
    int khit = 0;
    unsigned int maxsize = vdets.size();
    for (unsigned int i = 0; i < maxsize; i++) {
      std::vector<typename T::const_iterator> hit = spr::findHit(hits, vdets[i]);
      double energy = 0;
      bool ok = false;
      for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
        if (hit[ihit] != hits->end()) {
          energy += hit[ihit]->energy();
          khit++;
          ok = true;
          if (debug)
            edm::LogVerbatim("IsoTrack") << "energyDetIdHCALneighbours:: Hit " << khit << " " << (HcalDetId)vdets[i]
                                         << " energy " << hit[ihit]->energy();
        }
      }
      if (debug)
        edm::LogVerbatim("IsoTrack") << "energyDetIdHCALneighbours::Detector " << i << " " << (HcalDetId)vdets[i]
                                     << " energy " << energy;
      if (ok) {
        int subdet = ((HcalDetId)(vdets[i].rawId())).subdet();
        double eThr = hbThr;
        if (subdet == (int)(HcalEndcap))
          eThr = heThr;
        else if (subdet == (int)(HcalForward))
          eThr = hfThr;
        else if (subdet == (int)(HcalOuter))
          eThr = hoThr;
        if (energy > eThr)
          energyDetId.push_back(std::pair<DetId, double>(vdets[i], energy));
      }
    }

    if (debug)
      edm::LogVerbatim("IsoTrack") << "energyDetIdHCALneighbours::EnergyDetId buffer size " << energyDetId.size();
    return energyDetId;
  }

}  // namespace spr