HCALProperties.cc

Go to the documentation of this file.

#include "FWCore/ParameterSet/interface/ParameterSet.h"

//This class header
#include "FastSimulation/CalorimeterProperties/interface/HCALProperties.h"
#include <cmath>
#include <iostream>
#include "FWCore/MessageLogger/interface/MessageLogger.h"


HCALProperties::HCALProperties(const edm::ParameterSet& fastDet) : CalorimeterProperties()
{

    edm::ParameterSet fastDetHCAL= fastDet.getParameter<edm::ParameterSet>("HadronicCalorimeterProperties");
  hOPi = fastDetHCAL.getParameter<double>("HCAL_PiOverE");
  spotFrac= fastDetHCAL.getParameter<double>("HCAL_Sampling");
  HCALAeff_= fastDetHCAL.getParameter<double>("HCALAeff");
  HCALZeff_= fastDetHCAL.getParameter<double>("HCALZeff");
  HCALrho_= fastDetHCAL.getParameter<double>("HCALrho");
  HCALradiationLengthIncm_= fastDetHCAL.getParameter< double>("HCALradiationLengthIncm");
  HCALradLenIngcm2_= fastDetHCAL.getParameter<double>("HCALradLenIngcm2");
  HCALmoliereRadius_= fastDetHCAL.getParameter<double>("HCALmoliereRadius");
  HCALcriticalEnergy_= fastDetHCAL.getParameter<double>("HCALcriticalEnergy");
  HCALinteractionLength_= fastDetHCAL.getParameter<double>("HCALinteractionLength");
  etatow_=fastDetHCAL.getParameter<std::vector<double>>("HCALetatow");
  hcalDepthLam_=fastDetHCAL.getParameter<std::vector<double>>("HCALDepthLam");

  // in principle this splitting into 42 bins may change with future detectors, but let's add a protection to make sure that differences are not typos in the configuration file:
  if (etatow_.size() != 42) std::cout << " HCALProperties::eta2ieta - WARNING: here we expect 42 entries instead of " << etatow_.size() << "; is the change intentional?" << std::endl;
  // splitting of  28-th tower is taken into account (2.65-2.853-3.0)
  if (hcalDepthLam_.size() != etatow_.size()-1) std::cout << " HCALProperties::eta2ieta - WARNING: the sizes of HCALetatow and HCALDepthLam should differ by 1 unit! HCALDepthLam has size " << hcalDepthLam_.size()<< " and HCALetatow has size " << etatow_.size() << std::endl;
  

}

double HCALProperties::getHcalDepth(double eta) const{

  int  ieta = eta2ieta(eta); 

  /* 
  std::cout << " HCALProperties::getHcalDepth for eta = " << eta 
        << "  returns lam.thickness = " << hcalDepthLam_[ieta] << std::endl;
  */

  return  hcalDepthLam_[ieta];

}



int HCALProperties::eta2ieta(double eta) const {
  // binary search in the array of towers eta edges
  int size = etatow_.size();

  double x = fabs(eta);
  int curr = size / 2;
  int step = size / 4;
  int iter;
  int prevdir = 0; 
  int actudir = 0; 

  for (iter = 0; iter < size ; iter++) {

    if( curr >= size || curr < 1 )
      std::cout <<  " HCALProperties::eta2ieta - wrong current index = "
        << curr << " !!!" << std::endl;

    if ((x <= etatow_[curr]) && (x > etatow_[curr-1])) break;
    prevdir = actudir;
    if(x > etatow_[curr]) {actudir =  1;}
    else                 {actudir = -1;}
    if(prevdir * actudir < 0) { if(step > 1) step /= 2;}
    curr += actudir * step;
    if(curr > size) curr = size;
    else { if(curr < 1) {curr = 1;}}

    /*
    std::cout << " HCALProperties::eta2ieta  end of iter." << iter 
          << " curr, etatow_[curr-1], etatow_[curr] = "
          << curr << " " << etatow_[curr-1] << " " << etatow_[curr] << std::endl;
    */
    
  }

  /*
  std::cout << " HCALProperties::eta2ieta  for input x = " << x 
        << "  found index = " << curr-1
        << std::endl;
  */
  
  return curr-1;
}