HcalLayerDepthMap.cc

Go to the documentation of this file.

#include "Geometry/HcalCommonData/interface/HcalLayerDepthMap.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Utilities/interface/Exception.h"
#include<algorithm>

//#define EDM_ML_DEBUG

HcalLayerDepthMap::HcalLayerDepthMap() {
  subdet_  = 0;
  ietaMin_ = ietaMax_ = 0;
  depthMin_= 99;
  depthMax_=-1;
  dep29C_  = 2;
  wtl0C_   = 1.;
}

HcalLayerDepthMap::~HcalLayerDepthMap() {}

void HcalLayerDepthMap::initialize(const int subdet, const int ietaMax, 
                   const int dep16C, const int dep29C, 
                   const double wtl0C, 
                   std::vector<int> const& iphi,
                   std::vector<int> const& ieta,
                   std::vector<int> const& layer,
                   std::vector<int> const& depth) {
  
  subdet_  = subdet;
  ietaMin_ = ietaMax_ = ietaMax;
  dep16C_  = dep16C;
  dep29C_  = dep29C;
  wtl0C_   = wtl0C;
  iphi_.insert(iphi_.end(),iphi.begin(),iphi.end());
  layer2Depth_.clear(); depth2LayerF_.clear(); depth2LayerB_.clear();
  depthMin_= 99;
  depthMax_=-1;
  for (unsigned int k=0; k<ieta.size(); ++k) {
    if (ieta[k]  < ietaMin_)  ietaMin_  = ieta[k];
    if (depth[k] < depthMin_) depthMin_ = depth[k];
    if (depth[k] > depthMax_) depthMax_ = depth[k];
  }
  //Assume ieta, layer, depth are in increasing order of ieta and depth
  for (unsigned int k1=0; k1<ieta.size(); ++k1) {
    int ietaMin = ieta[k1];
    int ietaMax = ietaMax_;
    int layMin  = layer[k1];
    int layMax  = (k1+1 < ieta.size()) ? (layer[k1+1]-1) : maxLayers_;
    for (unsigned int k2=k1+1; k2<ieta.size(); ++k2) {
      if (ieta[k2] > ieta[k1]) {
    ietaMax = ieta[k2] - 1;
    if (k2 == k1+1) layMax = maxLayers_;
    break;
      }
    }
    for (int eta=ietaMin; eta<=ietaMax; ++eta) {
      depth2LayerF_[std::pair<int,int>(eta,depth[k1])] = layMin;
      depth2LayerB_[std::pair<int,int>(eta,depth[k1])] = layMax;
      for (int lay=layMin; lay<=layMax; ++lay)
    layer2Depth_[std::pair<int,int>(eta,lay)] = depth[k1];
    }
  }
  for (int eta=ietaMin_; eta<=ietaMax_; ++eta) {
    int dmin(99), dmax(-1);
    for (auto & itr : layer2Depth_) {
      if ((itr.first).first == eta) {
    if ((itr.second) < dmin) dmin = (itr.second);
    if ((itr.second) > dmax) dmax = (itr.second);
      }
    }
    if (subdet == 2) {
      if      (eta == ietaMin_) dmin = dep16C_;
      else if (eta == ietaMax_) dmax = dep29C_;
    }
    depthsEta_[eta] = std::pair<int,int>(dmin,dmax);
  }
#ifdef EDM_ML_DEBUG
  std::cout << "HcalLayerDepthMap: Subdet " << subdet_ << " iEta " << ietaMin_
        << ":" << ietaMax_ << " depth " << depthMin_ << ":" << depthMax_
        << "\nMaximum Depth for last HE towers " << dep29C_  
        << " Layer 0 Weight " << wtl0C_ << " iPhi";
  for (unsigned int k=0; k<iphi_.size(); ++k) std::cout << ":" << iphi_[k];
  std::cout << std::endl << "Layer2Depth_ with " << layer2Depth_.size()
        << " elements" << std::endl;
  for (std::map<std::pair<int,int>,int>::iterator itr=layer2Depth_.begin();
       itr != layer2Depth_.end(); ++itr)
    std::cout << "iEta " << (itr->first).first << " Layer " 
          << (itr->first).second << " Depth " << itr->second << std::endl;
  std::cout << "Depth2LayerFront with " << depth2LayerF_.size() 
        << " elemsts" << std::endl;
  for (std::map<std::pair<int,int>,int>::iterator itr=depth2LayerF_.begin();
       itr != depth2LayerF_.end(); ++itr)
    std::cout << "iEta " << (itr->first).first << " Depth " 
          << (itr->first).second << " Layer " << itr->second << std::endl;
  std::cout << "Depth2LayerBack with " << depth2LayerB_.size()
        << " elemets" << std::endl;
  for (std::map<std::pair<int,int>,int>::iterator itr=depth2LayerB_.begin();
       itr != depth2LayerB_.end(); ++itr)
    std::cout << "iEta " << (itr->first).first << " Depth " 
          << (itr->first).second << " Layer " << itr->second << std::endl;
  std::cout << "DepthsEta_ with " << depthsEta_.size() << " elements\n";
  for (std::map<int,std::pair<int,int> >::iterator itr=depthsEta_.begin();
       itr != depthsEta_.end(); ++itr)
    std::cout << "iEta " << itr->first << " Depths " 
          << (itr->second).first << ":" << (itr->second).second 
          << std::endl;
#endif
}

int HcalLayerDepthMap::getDepth(const int subdet, const int ieta, 
                const int iphi, const int zside,
                const int layer) const {
  int depth(-1);
  if (isValid(subdet,iphi,zside)) {
    std::map<std::pair<int,int>,int>::const_iterator itr = layer2Depth_.find(std::pair<int,int>(ieta,layer));
    if (itr != layer2Depth_.end()) depth = itr->second;
  }
#ifdef EDM_ML_DEBUG
  std::cout << "getDepth::Input " << subdet << ":" << ieta << ":" << iphi 
        << ":" << zside << ":" << layer << " Output " << depth <<std::endl;
#endif
  return depth;
}

int HcalLayerDepthMap::getDepth16(const int subdet, const int iphi, 
                  const int zside) const {
  int depth(-1);
  if (isValid(subdet,iphi,zside)) depth = dep16C_;
#ifdef EDM_ML_DEBUG
  std::cout << "getDepth16::Input " << subdet << ":" << iphi << ":" << zside
        << " Output "  << depth << std::endl;
#endif
  return depth;
}

int HcalLayerDepthMap::getDepthMin(const int subdet, const int iphi,
                   const int zside) const {
  int depth = (isValid(subdet,iphi,zside)) ? depthMin_ : -1;
#ifdef EDM_ML_DEBUG
  std::cout << "getDepthMin::Input " << subdet << ":" << iphi << ":" << zside
        << " Output " << depth << std::endl;
#endif
  return depth;
}

int HcalLayerDepthMap::getDepthMax(const int subdet, const int iphi,
                   const int zside) const {
  int depth = (isValid(subdet,iphi,zside)) ? depthMax_ : -1;
#ifdef EDM_ML_DEBUG
  std::cout << "getDepthMax::Input " << subdet << ":" << iphi << ":" << zside
        << " Output " << depth << std::endl;
#endif
  return depth;
}

int HcalLayerDepthMap::getDepthMax(const int subdet, const int ieta,
                   const int iphi, const int zside) const {
  int depth = (isValid(subdet,iphi,zside)) ? getDepth(subdet,ieta,iphi,zside,maxLayers_) : -1;
#ifdef EDM_ML_DEBUG
  std::cout << "getDepthMax::Input " << subdet << ":" << iphi << ":" << zside
        << " Output " << depth << std::endl;
#endif
  return depth;
}

std::pair<int,int> HcalLayerDepthMap::getDepths(const int eta) const {

  std::map<int,std::pair<int,int> >::const_iterator itr = depthsEta_.find(eta);
  if (itr == depthsEta_.end()) return std::pair<int,int>(-1,-1);
  else                         return itr->second;
}

int HcalLayerDepthMap::getLayerFront(const int subdet, const int ieta,
                     const int iphi, const int zside,
                     const int depth) const {
  int layer(-1);
  if (isValid(subdet,iphi,zside)) {
    std::map<std::pair<int,int>,int>::const_iterator itr = depth2LayerF_.find(std::pair<int,int>(ieta,depth));
    if (itr != depth2LayerF_.end()) layer = itr->second;
  }
#ifdef EDM_ML_DEBUG
  std::cout << "getLayerFront::Input " << subdet << ":" << ieta << ":" << iphi 
        << ":" << zside << ":" << depth << " Output " << layer <<std::endl;
#endif
  return layer;
}

int HcalLayerDepthMap::getLayerBack(const int subdet, const int ieta, 
                    const int iphi, const int zside,
                    const int depth) const {
  int layer(-1);
  if (isValid(subdet,iphi,zside)) {
    std::map<std::pair<int,int>,int>::const_iterator itr = depth2LayerB_.find(std::pair<int,int>(ieta,depth));
    if (itr != depth2LayerB_.end()) layer = itr->second;
  }
#ifdef EDM_ML_DEBUG
  std::cout << "getLayerBack::Input " << subdet << ":" << ieta << ":" << iphi 
        << ":" << zside << ":" << depth << " Output " << layer <<std::endl;
#endif
  return layer;
}

void HcalLayerDepthMap::getLayerDepth(const int subdet, const int eta, 
                      const int phi, const int zside,
                      std::map<int,int>& layers) const {
  layers.clear();
  if (isValid(subdet,phi,zside)) {
    for (const auto & itr : layer2Depth_) {
      if ((itr.first).first == eta) {
    layers[((itr.first).second)+1] = (itr.second);
      }
    }
  }
#ifdef EDM_ML_DEBUG
  std::cout << "getLayerDepth::Input " << subdet << ":" << eta << ":" << phi
        << ":" << zside << " Output " << layers.size() << " entries\n";
  for (std::map<int,int>::iterator itr=layers.begin(); itr != layers.end();
       ++itr) std::cout << " [" << itr->first << "] " << itr->second;
  std::cout << std::endl;
#endif
}

void HcalLayerDepthMap::getLayerDepth(const int eta, 
                      std::map<int,int>& layers) const {
  layers.clear();
  if (subdet_ > 0) {
    for (const auto & itr : layer2Depth_) {
      if ((itr.first).first == eta) {
    layers[((itr.first).second)+1] = (itr.second);
      }
    }
  }
#ifdef EDM_ML_DEBUG
  std::cout << "getLayerDepth::Input " << eta << " Output " << layers.size()
        << " entries\n";
  for (std::map<int,int>::iterator itr=layers.begin(); itr != layers.end();
       ++itr) std::cout << " [" << itr->first << "] " << itr->second;
  std::cout << std::endl;
#endif
}

int HcalLayerDepthMap::getMaxDepthLastHE(const int subdet, const int iphi,
                     const int zside) const {
  int depth = isValid(subdet,iphi,zside) ? dep29C_ : -1;
#ifdef EDM_ML_DEBUG
  std::cout << "getMaxDepthLastHE::Input " << subdet << ":" << iphi << ":"
        << zside << " Output " << depth << std::endl;
#endif
  return depth;
}

double HcalLayerDepthMap::getLayer0Wt(const int subdet, const int iphi,
                      const int zside) const {
  double wt = isValid(subdet,iphi,zside) ? wtl0C_ : -1.0;
#ifdef EDM_ML_DEBUG
  std::cout << "getLayer0Wt::Input " << subdet << ":" << iphi << ":" << zside 
        << " Output " << wt << std::endl;
#endif
  return wt;
}

bool HcalLayerDepthMap::isValid(const int subdet, const int iphi,
                const int zside) const {
  bool flag(false);
  int  kphi = (zside > 0) ? iphi : -iphi;
  if (subdet == subdet_)
    flag = (std::find(iphi_.begin(),iphi_.end(),kphi) != iphi_.end());
  return flag;
}

int HcalLayerDepthMap::validDet(std::vector<int>& phi) const {

  phi.insert(phi.end(),iphi_.begin(),iphi_.end());
  return subdet_;
}