EcalCalibAlgos/src/MatrixFillMap.cc

0001 #include "Calibration/EcalCalibAlgos/interface/MatrixFillMap.h"
0002 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0003 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0004 #include "FWCore/Utilities/interface/isFinite.h"
0005
0006 MatrixFillMap::MatrixFillMap(int WindowX,
0007                              int WindowY,
0008                              const std::map<int, int>& xtalReg,
0009                              double minE,
0010                              double maxE,
0011                              const std::map<int, int>& IndexReg,
0012                              EcalIntercalibConstantMap* barrelMap,
0013                              EcalIntercalibConstantMap* endcapMap)
0014     :
0015
0016       VFillMap(WindowX, WindowY, xtalReg, minE, maxE, IndexReg, barrelMap, endcapMap) {}
0017
0018 MatrixFillMap::~MatrixFillMap() {}
0019
0020 void MatrixFillMap::fillMap(const std::vector<std::pair<DetId, float> >& v1,
0021                             const DetId Max,
0022                             const EcalRecHitCollection* barrelHitsCollection,
0023                             const EcalRecHitCollection* endcapHitsCollection,
0024                             std::map<int, double>& xtlMap,
0025                             double& pSubtract) {
0026   if (Max.subdetId() == EcalBarrel) {
0027     EBDetId EBMax = Max;
0028     fillEBMap(EBMax, barrelHitsCollection, xtlMap, m_xtalRegionId[Max.rawId()], pSubtract);
0029   } else if (Max.subdetId() == EcalEndcap) {
0030     EEDetId EEMax = Max;
0031     fillEEMap(EEMax, endcapHitsCollection, xtlMap, m_xtalRegionId[Max.rawId()], pSubtract);
0032   }
0033 }
0034
0035 void MatrixFillMap::fillEBMap(EBDetId EBmax,
0036                               const EcalRecHitCollection* barrelHitsCollection,
0037                               std::map<int, double>& EBRegionMap,
0038                               int EBNumberOfRegion,
0039                               double& pSubtract) {
0040   int curr_eta;
0041   int curr_phi;
0042   //reads the hits in a recoWindowSide^2 wide region around the MOX

0043   for (int ii = 0; ii < m_recoWindowSidex; ++ii)
0044     for (int ij = 0; ij < m_recoWindowSidey; ++ij) {
0045       curr_eta = EBmax.ieta() + ii - (m_recoWindowSidex / 2);
0046       curr_phi = EBmax.iphi() + ij - (m_recoWindowSidey / 2);
0047       //skips if the xtals matrix falls over the border

0048       if (abs(curr_eta) > 85)
0049         continue;
0050       //Couples with the zero gap in the barrel eta index

0051       if (curr_eta * EBmax.ieta() <= 0) {
0052         if (EBmax.ieta() > 0)
0053           curr_eta--;
0054         else
0055           curr_eta++;
0056       }  // JUMP over 0

0057       //The following 2 couples with the ciclicity of the phiIndex

0058       if (curr_phi < 1)
0059         curr_phi += 360;
0060       if (curr_phi >= 360)
0061         curr_phi -= 360;
0062       //checks if the detId is valid

0063       if (EBDetId::validDetId(curr_eta, curr_phi)) {
0064         EBDetId det = EBDetId(curr_eta, curr_phi, EBDetId::ETAPHIMODE);
0065         int ID = det.rawId();
0066         //finds the hit corresponding to the cell

0067         EcalRecHitCollection::const_iterator curr_recHit = barrelHitsCollection->find(det);
0068         double dummy = 0;
0069         dummy = curr_recHit->energy();
0070         //checks if the reading of the xtal is in a sensible range

0071         if (edm::isNotFinite(dummy)) {
0072           dummy = 0;
0073         }
0074         if (dummy < m_minEnergyPerCrystal)
0075           continue;
0076         if (dummy > m_maxEnergyPerCrystal)
0077           continue;
0078         //corrects the energy with the calibration coeff of the ring

0079         dummy *= (*m_barrelMap)[det];
0080         //sums the energy of the xtal to the appropiate ring

0081         if (m_xtalRegionId[ID] == EBNumberOfRegion)
0082           EBRegionMap[m_IndexInRegion[ID]] += dummy;
0083         //adds the reading to pSubtract when part of the matrix is outside the region

0084         else
0085           pSubtract += dummy;
0086       }
0087     }
0088 }
0089
0090 void MatrixFillMap::fillEEMap(EEDetId EEmax,
0091                               const EcalRecHitCollection* endcapHitsCollection,
0092                               std::map<int, double>& EExtlMap,
0093                               int EENumberOfRegion,
0094                               double& pSubtract) {
0095   int curr_x;
0096   int curr_y;
0097   for (int ii = 0; ii < m_recoWindowSidex; ++ii)
0098     for (int ij = 0; ij < m_recoWindowSidey; ++ij) {
0099       //Works as fillEBMap

0100       curr_x = EEmax.ix() - m_recoWindowSidex / 2 + ii;
0101       curr_y = EEmax.iy() - m_recoWindowSidey / 2 + ij;
0102       if (EEDetId::validDetId(curr_x, curr_y, EEmax.zside())) {
0103         EEDetId det = EEDetId(curr_x, curr_y, EEmax.zside(), EEDetId::XYMODE);
0104         EcalRecHitCollection::const_iterator curr_recHit = endcapHitsCollection->find(det);
0105         double dummy = curr_recHit->energy();
0106         if (edm::isNotFinite(dummy)) {
0107           dummy = 0;
0108         }
0109         if (dummy < m_minEnergyPerCrystal)
0110           continue;
0111         if (dummy > m_maxEnergyPerCrystal) {
0112           continue;
0113         }
0114         dummy *= (*m_endcapMap)[det];
0115         int ID = det.rawId();
0116         if (m_xtalRegionId[ID] == EENumberOfRegion)
0117           EExtlMap[m_IndexInRegion[ID]] += dummy;
0118         else
0119           pSubtract += dummy;
0120       }
0121     }
0122 }