EcalRingCalibrationTools.cc

Go to the documentation of this file.

#include "Calibration/Tools/interface/EcalRingCalibrationTools.h"
#include "DataFormats/DetId/interface/DetId.h"
 
//Includes need to read from geometry
#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
 
#include "FWCore/MessageLogger/interface/MessageLogger.h"
 
#include <iostream>
 
//by default is not initialized, gets initialized at first call
std::atomic<bool> EcalRingCalibrationTools::isInitializedFromGeometry_(false);
short EcalRingCalibrationTools::endcapRingIndex_[EEDetId::IX_MAX][EEDetId::IY_MAX];
std::once_flag EcalRingCalibrationTools::once_;
 
constexpr short EcalRingCalibrationTools::N_RING_TOTAL;
constexpr short EcalRingCalibrationTools::N_RING_BARREL;
constexpr short EcalRingCalibrationTools::N_RING_ENDCAP;
 
constexpr short EcalRingCalibrationTools::N_MODULES_BARREL;
 
short EcalRingCalibrationTools::getRingIndex(DetId id) {
  if (id.det() != DetId::Ecal)
    return -1;
 
  if (id.subdetId() == EcalBarrel) {
    if (EBDetId(id).ieta() < 0)
      return EBDetId(id).ieta() + 85;
    else
      return EBDetId(id).ieta() + 84;
  }
  if (id.subdetId() == EcalEndcap) {
    //needed only for the EE, it can be replaced at some point with something smarter
    if (not isInitializedFromGeometry_)
      throw std::logic_error(
          "EcalRingCalibrationTools::initializeFromGeometry Ecal Endcap geometry is not initialized");
 
    EEDetId eid(id);
    short endcapRingIndex = endcapRingIndex_[eid.ix() - 1][eid.iy() - 1] + N_RING_BARREL;
    if (eid.zside() == 1)
      endcapRingIndex += N_RING_ENDCAP / 2;
    return endcapRingIndex;
  }
  return -1;
}
 
short EcalRingCalibrationTools::getModuleIndex(DetId id) {
  if (id.det() != DetId::Ecal)
    return -1;
 
  if (id.subdetId() == EcalBarrel) {
    short module = 4 * (EBDetId(id).ism() - 1) + EBDetId(id).im() - 1;
 
    //      std::cout<<"SM construction # : "<<EBDetId(id).ism() <<" SM installation # : "<<  installationSMNumber[ EBDetId(id).ism() -1 ]<< "Xtal is Module :"<<module<< std::endl;
    return module;
  }
  if (id.subdetId() == EcalEndcap) {
    return -1;
  }
 
  return -1;
}
 
std::vector<DetId> EcalRingCalibrationTools::getDetIdsInRing(short etaIndex) {
  std::vector<DetId> ringIds;
  if (etaIndex < 0)
    return ringIds;
 
  if (etaIndex < N_RING_BARREL) {
    int k = 0;
    if (etaIndex < 85)
      k = -85 + etaIndex;
    else
      k = etaIndex - 84;
 
    for (int iphi = EBDetId::MIN_IPHI; iphi <= EBDetId::MAX_IPHI; ++iphi)
 
      if (EBDetId::validDetId(k, iphi))
        ringIds.push_back(EBDetId(k, iphi));
  }
 
  else if (etaIndex < N_RING_TOTAL) {
    //needed only for the EE, it can be replaced at some point maybe with something smarter
    if (not isInitializedFromGeometry_)
      throw std::logic_error(
          "EcalRingCalibrationTools::initializeFromGeometry Ecal Endcap geometry is not initialized");
 
    int zside = (etaIndex < N_RING_BARREL + (N_RING_ENDCAP / 2)) ? -1 : 1;
    short eeEtaIndex = (etaIndex - N_RING_BARREL) % (N_RING_ENDCAP / 2);
 
    for (int ix = 0; ix < EEDetId::IX_MAX; ++ix)
      for (int iy = 0; iy < EEDetId::IY_MAX; ++iy)
        if (endcapRingIndex_[ix][iy] == eeEtaIndex)
          ringIds.push_back(EEDetId(ix + 1, iy + 1, zside));
  }
 
  return ringIds;
}
 
std::vector<DetId> EcalRingCalibrationTools::getDetIdsInECAL() {
  std::vector<DetId> ringIds;
 
  for (int ieta = -EBDetId::MAX_IETA; ieta <= EBDetId::MAX_IETA; ++ieta)
    for (int iphi = EBDetId::MIN_IPHI; iphi <= EBDetId::MAX_IPHI; ++iphi)
      if (EBDetId::validDetId(ieta, iphi))
        ringIds.push_back(EBDetId(ieta, iphi));
 
  //needed only for the EE, it can be replaced at some point maybe with something smarter
  if (not isInitializedFromGeometry_)
    throw std::logic_error("EcalRingCalibrationTools::initializeFromGeometry Ecal Endcap geometry is not initialized");
 
  for (int ix = 0; ix < EEDetId::IX_MAX; ++ix)
    for (int iy = 0; iy < EEDetId::IY_MAX; ++iy)
      for (int zside = -1; zside < 2; zside += 2)
        if (EEDetId::validDetId(ix + 1, iy + 1, zside))
          ringIds.push_back(EEDetId(ix + 1, iy + 1, zside));
 
  //  std::cout<<" [EcalRingCalibrationTools::getDetIdsInECAL()] DetId.size() is "<<ringIds.size()<<std::endl;
 
  return ringIds;
}
 
std::vector<DetId> EcalRingCalibrationTools::getDetIdsInModule(short moduleIndex) {
  std::vector<DetId> ringIds;
  if (moduleIndex < 0)
    return ringIds;
 
  short moduleBound[5] = {1, 26, 46, 66, 86};
  if (moduleIndex < EcalRingCalibrationTools::N_MODULES_BARREL) {
    short sm, moduleInSm, zsm;
 
    short minModuleiphi, maxModuleiphi, minModuleieta = 360, maxModuleieta = 0;
 
    //  if(moduleIndex%4 != 0 )
    sm = moduleIndex / 4 + 1;
    //  else
    //sm = moduleIndex/4;//i.e. module 8 belongs to sm=3, not sm=3
 
    //if(moduleIndex%4 != 0 )
    moduleInSm = moduleIndex % 4;
    //else
    //moduleInSm = 4;//moduleInSm is [1,2,3,4]
 
    if (moduleIndex > 71)
      zsm = -1;
    else
      zsm = 1;
 
    minModuleiphi = ((sm - 1) % 18 + 1) * 20 - 19;
    maxModuleiphi = ((sm - 1) % 18 + 1) * 20;
 
    if (zsm == 1) {
      minModuleieta = moduleBound[moduleInSm];
      maxModuleieta = moduleBound[moduleInSm + 1] - 1;
    } else if (zsm == -1) {
      minModuleieta = -moduleBound[moduleInSm + 1] + 1;
      maxModuleieta = -moduleBound[moduleInSm];
    }
 
    std::cout << "Called moduleIndex " << moduleIndex << std::endl;
    std::cout << "minModuleieta " << minModuleieta << " maxModuleieta " << maxModuleieta << " minModuleiphi "
              << minModuleiphi << " maxModuleiphi " << maxModuleiphi << std::endl;
 
    for (int ieta = minModuleieta; ieta <= maxModuleieta; ++ieta) {
      for (int iphi = minModuleiphi; iphi <= maxModuleiphi; ++iphi) {
        ringIds.push_back(EBDetId(ieta, iphi));
 
        //        std::cout<<"Putting Xtal with ieta: "<<ieta<<" iphi "<<iphi<<" of SM "<<sm<<" into Module "<<moduleIndex<<std::endl;
 
      }  //close loop on phi
    }    //close loop on eta
  }      //close if ( moduleInstallationNumber < 144)
 
  return ringIds;
}
 
void EcalRingCalibrationTools::setCaloGeometry(const CaloGeometry* geometry) {
  std::call_once(once_, EcalRingCalibrationTools::initializeFromGeometry, geometry);
}
 
void EcalRingCalibrationTools::initializeFromGeometry(CaloGeometry const* geometry) {
  if (not geometry)
    throw std::invalid_argument("EcalRingCalibrationTools::initializeFromGeometry called with a nullptr argument");
 
  float cellPosEta[EEDetId::IX_MAX][EEDetId::IY_MAX];
  for (int ix = 0; ix < EEDetId::IX_MAX; ++ix)
    for (int iy = 0; iy < EEDetId::IY_MAX; ++iy) {
      cellPosEta[ix][iy] = -1.;
      endcapRingIndex_[ix][iy] = -9;
    }
 
  CaloSubdetectorGeometry const* endcapGeometry = geometry->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
  if (not endcapGeometry)
    throw std::logic_error("EcalRingCalibrationTools::initializeFromGeometry Ecal Endcap geometry not found");
 
  std::vector<DetId> const& endcapCells = geometry->getValidDetIds(DetId::Ecal, EcalEndcap);
 
  for (std::vector<DetId>::const_iterator endcapIt = endcapCells.begin(); endcapIt != endcapCells.end(); ++endcapIt) {
    EEDetId ee(*endcapIt);
    if (ee.zside() == -1)
      continue;  //Just using +side to fill absEta x,y map
    auto cellGeometry = endcapGeometry->getGeometry(*endcapIt);
    int ics = ee.ix() - 1;
    int ips = ee.iy() - 1;
    cellPosEta[ics][ips] = fabs(cellGeometry->getPosition().eta());
    //std::cout<<"EE Xtal, |eta| is "<<fabs(cellGeometry->getPosition().eta())<<std::endl;
  }
 
  float eta_ring[N_RING_ENDCAP / 2];
  for (int ring = 0; ring < N_RING_ENDCAP / 2; ++ring)
    eta_ring[ring] = cellPosEta[ring][50];
 
  double etaBoundary[N_RING_ENDCAP / 2 + 1];
  etaBoundary[0] = 1.47;
  etaBoundary[N_RING_ENDCAP / 2] = 4.0;
 
  for (int ring = 1; ring < N_RING_ENDCAP / 2; ++ring)
    etaBoundary[ring] = (eta_ring[ring] + eta_ring[ring - 1]) / 2.;
 
  for (int ring = 0; ring < N_RING_ENDCAP / 2; ++ring) {
    // std::cout<<"***********************EE ring: "<<ring<<" eta "<<(etaBoundary[ring] + etaBoundary[ring+1])/2.<<std::endl;
    for (int ix = 0; ix < EEDetId::IX_MAX; ix++)
      for (int iy = 0; iy < EEDetId::IY_MAX; iy++)
        if (cellPosEta[ix][iy] > etaBoundary[ring] && cellPosEta[ix][iy] < etaBoundary[ring + 1]) {
          endcapRingIndex_[ix][iy] = ring;
          //std::cout<<"endcapRing_["<<ix+1<<"]["<<iy+1<<"] = "<<ring<<";"<<std::endl;
        }
  }
 
  std::vector<DetId> const& barrelCells = geometry->getValidDetIds(DetId::Ecal, EcalBarrel);
 
  for (std::vector<DetId>::const_iterator barrelIt = barrelCells.begin(); barrelIt != barrelCells.end(); ++barrelIt) {
    EBDetId eb(*barrelIt);
  }
 
  //EB
 
  isInitializedFromGeometry_ = true;
}