MatrixFillMap.cc

Go to the documentation of this file.

#include "Calibration/EcalCalibAlgos/interface/MatrixFillMap.h"
#include "DataFormats/EcalDetId/interface/EEDetId.h"
#include "DataFormats/EcalDetId/interface/EBDetId.h"
#include "FWCore/Utilities/interface/isFinite.h"

MatrixFillMap::MatrixFillMap(int WindowX,
                             int WindowY,
                             const std::map<int, int>& xtalReg,
                             double minE,
                             double maxE,
                             const std::map<int, int>& IndexReg,
                             EcalIntercalibConstantMap* barrelMap,
                             EcalIntercalibConstantMap* endcapMap)
    :

      VFillMap(WindowX, WindowY, xtalReg, minE, maxE, IndexReg, barrelMap, endcapMap) {}

MatrixFillMap::~MatrixFillMap() {}

void MatrixFillMap::fillMap(const std::vector<std::pair<DetId, float> >& v1,
                            const DetId Max,
                            const EcalRecHitCollection* barrelHitsCollection,
                            const EcalRecHitCollection* endcapHitsCollection,
                            std::map<int, double>& xtlMap,
                            double& pSubtract) {
  if (Max.subdetId() == EcalBarrel) {
    EBDetId EBMax = Max;
    fillEBMap(EBMax, barrelHitsCollection, xtlMap, m_xtalRegionId[Max.rawId()], pSubtract);
  } else if (Max.subdetId() == EcalEndcap) {
    EEDetId EEMax = Max;
    fillEEMap(EEMax, endcapHitsCollection, xtlMap, m_xtalRegionId[Max.rawId()], pSubtract);
  }
}

void MatrixFillMap::fillEBMap(EBDetId EBmax,
                              const EcalRecHitCollection* barrelHitsCollection,
                              std::map<int, double>& EBRegionMap,
                              int EBNumberOfRegion,
                              double& pSubtract) {
  int curr_eta;
  int curr_phi;
  //reads the hits in a recoWindowSide^2 wide region around the MOX
  for (int ii = 0; ii < m_recoWindowSidex; ++ii)
    for (int ij = 0; ij < m_recoWindowSidey; ++ij) {
      curr_eta = EBmax.ieta() + ii - (m_recoWindowSidex / 2);
      curr_phi = EBmax.iphi() + ij - (m_recoWindowSidey / 2);
      //skips if the xtals matrix falls over the border
      if (abs(curr_eta) > 85)
        continue;
      //Couples with the zero gap in the barrel eta index
      if (curr_eta * EBmax.ieta() <= 0) {
        if (EBmax.ieta() > 0)
          curr_eta--;
        else
          curr_eta++;
      }  // JUMP over 0
      //The following 2 couples with the ciclicity of the phiIndex
      if (curr_phi < 1)
        curr_phi += 360;
      if (curr_phi >= 360)
        curr_phi -= 360;
      //checks if the detId is valid
      if (EBDetId::validDetId(curr_eta, curr_phi)) {
        EBDetId det = EBDetId(curr_eta, curr_phi, EBDetId::ETAPHIMODE);
        int ID = det.rawId();
        //finds the hit corresponding to the cell
        EcalRecHitCollection::const_iterator curr_recHit = barrelHitsCollection->find(det);
        double dummy = 0;
        dummy = curr_recHit->energy();
        //checks if the reading of the xtal is in a sensible range
        if (edm::isNotFinite(dummy)) {
          dummy = 0;
        }
        if (dummy < m_minEnergyPerCrystal)
          continue;
        if (dummy > m_maxEnergyPerCrystal)
          continue;
        //corrects the energy with the calibration coeff of the ring
        dummy *= (*m_barrelMap)[det];
        //sums the energy of the xtal to the appropiate ring
        if (m_xtalRegionId[ID] == EBNumberOfRegion)
          EBRegionMap[m_IndexInRegion[ID]] += dummy;
        //adds the reading to pSubtract when part of the matrix is outside the region
        else
          pSubtract += dummy;
      }
    }
}

void MatrixFillMap::fillEEMap(EEDetId EEmax,
                              const EcalRecHitCollection* endcapHitsCollection,
                              std::map<int, double>& EExtlMap,
                              int EENumberOfRegion,
                              double& pSubtract) {
  int curr_x;
  int curr_y;
  for (int ii = 0; ii < m_recoWindowSidex; ++ii)
    for (int ij = 0; ij < m_recoWindowSidey; ++ij) {
      //Works as fillEBMap
      curr_x = EEmax.ix() - m_recoWindowSidex / 2 + ii;
      curr_y = EEmax.iy() - m_recoWindowSidey / 2 + ij;
      if (EEDetId::validDetId(curr_x, curr_y, EEmax.zside())) {
        EEDetId det = EEDetId(curr_x, curr_y, EEmax.zside(), EEDetId::XYMODE);
        EcalRecHitCollection::const_iterator curr_recHit = endcapHitsCollection->find(det);
        double dummy = curr_recHit->energy();
        if (edm::isNotFinite(dummy)) {
          dummy = 0;
        }
        if (dummy < m_minEnergyPerCrystal)
          continue;
        if (dummy > m_maxEnergyPerCrystal) {
          continue;
        }
        dummy *= (*m_endcapMap)[det];
        int ID = det.rawId();
        if (m_xtalRegionId[ID] == EENumberOfRegion)
          EExtlMap[m_IndexInRegion[ID]] += dummy;
        else
          pSubtract += dummy;
      }
    }
}