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;

}