HCALResponse Class Reference

#include <HCALResponse.h>

Public Member Functions
void	correctHF (double e, int type)
vec1 &	getCorrHFem ()
vec1 &	getCorrHFhad ()
double	getHCALEnergyResponse (double e, int hit, RandomEngineAndDistribution const *)
double	getMIPfraction (double energy, double eta)
	HCALResponse (const edm::ParameterSet &pset)
double	responseHCAL (int _mip, double energy, double eta, int partype, RandomEngineAndDistribution const *)
	~HCALResponse ()

Private Member Functions
double	cballShootNoNegative (double mu, double sigma, double aL, double nL, double aR, double nR, RandomEngineAndDistribution const *)
double	gaussShootNoNegative (double e, double sigma, RandomEngineAndDistribution const *)
int	getDet (int ieta)
double	interEM (double e, int ie, int ieta, RandomEngineAndDistribution const *)
double	interHD (int mip, double e, int ie, int ieta, int det, RandomEngineAndDistribution const *)
double	interMU (double e, int ie, int ieta, RandomEngineAndDistribution const *)
double	PoissonShootNoNegative (double e, double sigma, RandomEngineAndDistribution const *)

Private Attributes
int	barrelMUeta
DoubleCrystalBallGenerator	cball
vec1	corrHFem
vec2	corrHFgEm
vec2	corrHFgHad
vec1	corrHFhad
vec2	corrHFhEm
vec2	corrHFhHad
bool	debug
vec1	eGridEM
vec1	eGridHD [4]
vec1	eGridMU
int	endcapMUeta
vec1	energyHF
double	eResponseCoefficient
double	eResponseExponent
double	eResponsePlateau [3]
double	eResponseScale [3]
vec1	etaGridMU
double	etaStep
int	HDeta [4]
int	maxEMe
int	maxEMeta
int	maxEne
int	maxEta
int	maxHDe [4]
int	maxHDetas [3]
int	maxMUbin
int	maxMUe
int	maxMUeta
vec2	meanEM
vec3	mipfraction
double	muStep
int	nPar
vec5	parameters
std::vector< std::string >	parNames
vec3	PoissonParameters
double	respFactorEM
vec3	responseMU
double	RespPar [3][2][3]
vec2	sigmaEM
bool	usemip

Detailed Description

Definition at line 29 of file HCALResponse.h.

Constructor & Destructor Documentation

HCALResponse()

HCALResponse::HCALResponse

(

const edm::ParameterSet &

pset

)

Definition at line 18 of file HCALResponse.cc.

References funct::abs(), HcalResponse_cfi::barrelMUeta, HcalResponse_cfi::corrHFgEm, HcalResponse_cfi::corrHFgHad, HcalResponse_cfi::corrHFhEm, HcalResponse_cfi::corrHFhHad, ztail::d, debug, HcalResponse_cfi::eGridEM, HcalResponse_cfi::eGridMU, HcalResponse_cfi::endcapMUeta, HcalResponse_cfi::energyHF, HcalResponse_cfi::eResponseCoefficient, HcalResponse_cfi::eResponseExponent, HcalResponse_cfi::etaGridMU, HcalResponse_cfi::etaStep, HLT_2022v11_cff::fraction, HCAL, mps_fire::i, dqmiolumiharvest::j, dqmdumpme::k, visualization-live-secondInstance_cfg::m, HcalResponse_cfi::maxEMe, HcalResponse_cfi::maxEne, maxEta, HcalResponse_cfi::maxMUbin, HcalResponse_cfi::maxMUe, HcalResponse_cfi::maxMUeta, HcalResponse_cfi::meanEM, HcalResponse_cfi::muStep, HcalResponse_cfi::nPar, AlCaHLTBitMon_ParallelJobs::p, HcalResponse_cfi::parNames, unpackData-CaloStage2::pname, muonDTDigis_cfi::pset, HcalResponse_cfi::respFactorEM, HcalResponse_cfi::sigmaEM, AlCaHLTBitMon_QueryRunRegistry::string, createJobs::tmp, HcalResponse_cfi::usemip, and VFCAL.

                                                      {
  //switches
  debug = pset.getParameter<bool>("debug");
  usemip = pset.getParameter<bool>("usemip");

  //values for "old" response parameterizations
  //--------------------------------------------------------------------
  RespPar[HCAL][0][0] = pset.getParameter<double>("HadronBarrelResolution_Stochastic");
  RespPar[HCAL][0][1] = pset.getParameter<double>("HadronBarrelResolution_Constant");
  RespPar[HCAL][0][2] = pset.getParameter<double>("HadronBarrelResolution_Noise");

  RespPar[HCAL][1][0] = pset.getParameter<double>("HadronEndcapResolution_Stochastic");
  RespPar[HCAL][1][1] = pset.getParameter<double>("HadronEndcapResolution_Constant");
  RespPar[HCAL][1][2] = pset.getParameter<double>("HadronEndcapResolution_Noise");

  RespPar[VFCAL][0][0] = pset.getParameter<double>("HadronForwardResolution_Stochastic");
  RespPar[VFCAL][0][1] = pset.getParameter<double>("HadronForwardResolution_Constant");
  RespPar[VFCAL][0][2] = pset.getParameter<double>("HadronForwardResolution_Noise");

  RespPar[VFCAL][1][0] = pset.getParameter<double>("ElectronForwardResolution_Stochastic");
  RespPar[VFCAL][1][1] = pset.getParameter<double>("ElectronForwardResolution_Constant");
  RespPar[VFCAL][1][2] = pset.getParameter<double>("ElectronForwardResolution_Noise");

  eResponseScale[0] = pset.getParameter<double>("eResponseScaleHB");
  eResponseScale[1] = pset.getParameter<double>("eResponseScaleHE");
  eResponseScale[2] = pset.getParameter<double>("eResponseScaleHF");

  eResponsePlateau[0] = pset.getParameter<double>("eResponsePlateauHB");
  eResponsePlateau[1] = pset.getParameter<double>("eResponsePlateauHE");
  eResponsePlateau[2] = pset.getParameter<double>("eResponsePlateauHF");

  eResponseExponent = pset.getParameter<double>("eResponseExponent");
  eResponseCoefficient = pset.getParameter<double>("eResponseCoefficient");

  //pion parameters
  //--------------------------------------------------------------------
  //energy values
  maxHDe[0] = pset.getParameter<int>("maxHBe");
  maxHDe[1] = pset.getParameter<int>("maxHEe");
  maxHDe[2] = pset.getParameter<int>("maxHFe");
  maxHDe[3] = pset.getParameter<int>("maxHFlowe");

  eGridHD[0] = pset.getParameter<vec1>("eGridHB");
  eGridHD[1] = pset.getParameter<vec1>("eGridHE");
  eGridHD[2] = pset.getParameter<vec1>("eGridHF");
  eGridHD[3] = pset.getParameter<vec1>("loweGridHF");

  //region eta indices calculated from eta values
  etaStep = pset.getParameter<double>("etaStep");
  //eta boundaries
  HDeta[0] = abs((int)(pset.getParameter<double>("HBeta") / etaStep));
  HDeta[1] = abs((int)(pset.getParameter<double>("HEeta") / etaStep));
  HDeta[2] = abs((int)(pset.getParameter<double>("HFeta") / etaStep));
  HDeta[3] = abs((int)(pset.getParameter<double>("maxHDeta") / etaStep));  //add 1 because this is the max index
  //eta ranges
  maxHDetas[0] = HDeta[1] - HDeta[0];
  maxHDetas[1] = HDeta[2] - HDeta[1];
  maxHDetas[2] = HDeta[3] - HDeta[2];

  //parameter info
  nPar = pset.getParameter<int>("nPar");
  parNames = pset.getParameter<std::vector<std::string> >("parNames");
  std::string detNames[] = {"_HB", "_HE", "_HF"};
  std::string mipNames[] = {"_mip", "_nomip", ""};
  std::string fraction = "f";
  //setup parameters (5D vector)
  parameters = vec5(nPar, vec4(3, vec3(3)));
  for (int p = 0; p < nPar; p++) {   //loop over parameters
    for (int m = 0; m < 3; m++) {    //loop over mip, nomip, total
      for (int d = 0; d < 3; d++) {  //loop over dets: HB, HE, HF
        //get from python
        std::string pname = parNames[p] + detNames[d] + mipNames[m];
        vec1 tmp = pset.getParameter<vec1>(pname);

        //resize vector for energy range of det d
        parameters[p][m][d].resize(maxHDe[d]);

        for (int i = 0; i < maxHDe[d]; i++) {  //loop over energy for det d
          //resize vector for eta range of det d
          parameters[p][m][d][i].resize(maxHDetas[d]);

          for (int j = 0; j < maxHDetas[d]; j++) {  //loop over eta for det d
            //fill in parameters vector from python
            parameters[p][m][d][i][j] = tmp[i * maxHDetas[d] + j];
          }
        }
      }
    }
  }
  //set up Poisson parameters for low energy Hadrons in HF
  //----------------------------------------------------------------------
  PoissonParameters = vec3(4);
  std::string PoissonParName[] = {"mean_overall", "shift_overall", "mean_between", "shift_between"};
  for (int d = 0; d < 4; d++) {  //loop over Poisson parameteres
    vec1 tmp1 = pset.getParameter<vec1>(PoissonParName[d]);
    for (int i = 0; i < maxHDe[3]; i++) {  //loop over energy for low HF energy points
      PoissonParameters[d].resize(maxHDe[3]);
      for (int j = 0; j < maxHDetas[2]; j++) {  //loop over HF eta points
        PoissonParameters[d][i].resize(maxHDetas[2]);
        PoissonParameters[d][i][j] = tmp1[i * maxHDetas[2] + j];
      }
    }
  }

  //MIP fraction fill in 3d vector
  mipfraction = vec3(3);
  for (int d = 0; d < 3; d++) {  //loop over dets: HB, HE, HF
    //get from python
    std::string mipname = fraction + mipNames[0] + detNames[d];
    vec1 tmp1 = pset.getParameter<vec1>(mipname);
    mipfraction[d].resize(maxHDe[d]);
    for (int i = 0; i < maxHDe[d]; i++) {  //loop over energy for det d
      //resize vector for eta range of det d
      mipfraction[d][i].resize(maxHDetas[d]);
      for (int j = 0; j < maxHDetas[d]; j++) {  //loop over eta for det d
        //fill in parameters vector from python
        mipfraction[d][i][j] = tmp1[i * maxHDetas[d] + j];
      }
    }
  }

  // MUON probability histos for bin size = 0.25 GeV (0-10 GeV, 40 bins)
  //--------------------------------------------------------------------
  muStep = pset.getParameter<double>("muStep");
  maxMUe = pset.getParameter<int>("maxMUe");
  maxMUeta = pset.getParameter<int>("maxMUeta");
  maxMUbin = pset.getParameter<int>("maxMUbin");
  eGridMU = pset.getParameter<vec1>("eGridMU");
  etaGridMU = pset.getParameter<vec1>("etaGridMU");
  vec1 _responseMU[2] = {pset.getParameter<vec1>("responseMUBarrel"), pset.getParameter<vec1>("responseMUEndcap")};

  //get muon region eta indices from the eta grid
  double _barrelMUeta = pset.getParameter<double>("barrelMUeta");
  double _endcapMUeta = pset.getParameter<double>("endcapMUeta");
  barrelMUeta = endcapMUeta = maxMUeta;
  for (int i = 0; i < maxMUeta; i++) {
    if (fabs(_barrelMUeta) <= etaGridMU[i]) {
      barrelMUeta = i;
      break;
    }
  }
  for (int i = 0; i < maxMUeta; i++) {
    if (fabs(_endcapMUeta) <= etaGridMU[i]) {
      endcapMUeta = i;
      break;
    }
  }
  int maxMUetas[] = {endcapMUeta - barrelMUeta, maxMUeta - endcapMUeta};

  //initialize 3D vector
  responseMU = vec3(maxMUe, vec2(maxMUeta, vec1(maxMUbin, 0)));

  //fill in 3D vector
  //(complementary cumulative distribution functions, from normalized response distributions)
  int loc, eta_loc;
  loc = eta_loc = -1;
  for (int i = 0; i < maxMUe; i++) {
    for (int j = 0; j < maxMUeta; j++) {
      //check location - barrel, endcap, or forward
      if (j == barrelMUeta) {
        loc = 0;
        eta_loc = barrelMUeta;
      } else if (j == endcapMUeta) {
        loc = 1;
        eta_loc = endcapMUeta;
      }

      for (int k = 0; k < maxMUbin; k++) {
        responseMU[i][j][k] = _responseMU[loc][i * maxMUetas[loc] * maxMUbin + (j - eta_loc) * maxMUbin + k];

        if (debug) {
          //cout.width(6);
          LogInfo("FastCalorimetry") << " responseMU " << i << " " << j << " " << k << " = " << responseMU[i][j][k]
                                     << std::endl;
        }
      }
    }
  }

  // values for EM response in HF
  //--------------------------------------------------------------------
  maxEMe = pset.getParameter<int>("maxEMe");
  maxEMeta = maxHDetas[2];
  respFactorEM = pset.getParameter<double>("respFactorEM");
  eGridEM = pset.getParameter<vec1>("eGridEM");

  // e-gamma mean response and sigma in HF
  vec1 _meanEM = pset.getParameter<vec1>("meanEM");
  vec1 _sigmaEM = pset.getParameter<vec1>("sigmaEM");

  //fill in 2D vectors (w/ correction factor applied)
  meanEM = vec2(maxEMe, vec1(maxEMeta, 0));
  sigmaEM = vec2(maxEMe, vec1(maxEMeta, 0));
  for (int i = 0; i < maxEMe; i++) {
    for (int j = 0; j < maxEMeta; j++) {
      meanEM[i][j] = respFactorEM * _meanEM[i * maxEMeta + j];
      sigmaEM[i][j] = respFactorEM * _sigmaEM[i * maxEMeta + j];
    }
  }

  // HF correction for SL
  //---------------------
  maxEta = pset.getParameter<int>("maxEta");
  maxEne = pset.getParameter<int>("maxEne");
  energyHF = pset.getParameter<vec1>("energyHF");
  vec1 _corrHFgEm = pset.getParameter<vec1>("corrHFgEm");
  vec1 _corrHFgHad = pset.getParameter<vec1>("corrHFgHad");
  vec1 _corrHFhEm = pset.getParameter<vec1>("corrHFhEm");
  vec1 _corrHFhHad = pset.getParameter<vec1>("corrHFhHad");
  corrHFem = vec1(maxEta, 0);
  corrHFhad = vec1(maxEta, 0);

  // initialize 2D vector corrHFgEm[energy][eta]
  corrHFgEm = vec2(maxEne, vec1(maxEta, 0));
  corrHFgHad = vec2(maxEne, vec1(maxEta, 0));
  corrHFhEm = vec2(maxEne, vec1(maxEta, 0));
  corrHFhHad = vec2(maxEne, vec1(maxEta, 0));
  // Fill
  for (int i = 0; i < maxEne; i++) {
    for (int j = 0; j < maxEta; j++) {
      corrHFgEm[i][j] = _corrHFgEm[i * maxEta + j];
      corrHFgHad[i][j] = _corrHFgHad[i * maxEta + j];
      corrHFhEm[i][j] = _corrHFhEm[i * maxEta + j];
      corrHFhHad[i][j] = _corrHFhHad[i * maxEta + j];
    }
  }
}

~HCALResponse()

HCALResponse::~HCALResponse ( )

inline

Definition at line 32 of file HCALResponse.h.

{}

Member Function Documentation

cballShootNoNegative()

double HCALResponse::cballShootNoNegative ( double  mu, double  sigma, double  aL, double  nL, double  aR, double  nR, RandomEngineAndDistribution const *  random  )

private

Definition at line 614 of file HCALResponse.cc.

References amptDefaultParameters_cff::mu, and MillePedeFileConverter_cfg::out.

                                                                                                                    {
  double out = cball.shoot(mu, sigma, aL, nL, aR, nR, random);
  if (mu >= 0.) {
    while (out < 0.)
      out = cball.shoot(mu, sigma, aL, nL, aR, nR, random);
  }
  //else give up on re-trying, otherwise too much time can be lost before emeas comes out positive

  return out;
}

correctHF()

void HCALResponse::correctHF	(	double	e,
		int	type
	)

Definition at line 634 of file HCALResponse.cc.

References funct::abs(), HcalResponse_cfi::corrHFgEm, HcalResponse_cfi::corrHFgHad, HcalResponse_cfi::corrHFhEm, HcalResponse_cfi::corrHFhHad, HcalResponse_cfi::energyHF, mps_fire::i, HcalResponse_cfi::maxEne, maxEta, testProducerWithPsetDescEmpty_cfi::x1, and testProducerWithPsetDescEmpty_cfi::x2.

Referenced by CalorimetryManager::HDShowerSimulation(), and CalorimetryManager::reconstructTrack().

                                                {
  int jmin = 0;
  for (int i = 0; i < maxEne; i++) {
    if (ee >= energyHF[i])
      jmin = i;
  }

  double x1, x2;
  double y1em, y2em;
  double y1had, y2had;
  for (int i = 0; i < maxEta; ++i) {
    if (ee < energyHF[0]) {
      if (abs(type) == 11 || abs(type) == 22) {
        corrHFem[i] = corrHFgEm[0][i];
        corrHFhad[i] = corrHFgHad[0][i];
      } else {
        corrHFem[i] = corrHFhEm[0][i];
        corrHFhad[i] = corrHFhHad[0][i];
      }
    } else if (jmin >= maxEne - 1) {
      if (abs(type) == 11 || abs(type) == 22) {
        corrHFem[i] = corrHFgEm[maxEta][i];
        corrHFhad[i] = corrHFgHad[maxEta][i];
      } else {
        corrHFem[i] = corrHFhEm[maxEta][i];
        corrHFhad[i] = corrHFhHad[maxEta][i];
      }
    } else {
      x1 = energyHF[jmin];
      x2 = energyHF[jmin + 1];
      if (abs(type) == 11 || abs(type) == 22) {
        y1em = corrHFgEm[jmin][i];
        y2em = corrHFgEm[jmin + 1][i];
        y1had = corrHFgHad[jmin][i];
        y2had = corrHFgHad[jmin + 1][i];
      } else {
        y1em = corrHFhEm[jmin][i];
        y2em = corrHFhEm[jmin + 1][i];
        y1had = corrHFhHad[jmin][i];
        y2had = corrHFhHad[jmin + 1][i];
      }
      corrHFem[i] = y1em + (ee - x1) * ((y2em - y1em) / (x2 - x1));
      corrHFhad[i] = y1had + (ee - x1) * ((y2had - y1had) / (x2 - x1));
    }
  }
}

gaussShootNoNegative()

double HCALResponse::gaussShootNoNegative ( double  e, double  sigma, RandomEngineAndDistribution const *  random  )

private

Definition at line 602 of file HCALResponse.cc.

References MillePedeFileConverter_cfg::e, RandomEngineAndDistribution::gaussShoot(), and MillePedeFileConverter_cfg::out.

                                                                                                           {
  double out = random->gaussShoot(e, sigma);
  if (e >= 0.) {
    while (out < 0.)
      out = random->gaussShoot(e, sigma);
  }
  //else give up on re-trying, otherwise too much time can be lost before emeas comes out positive

  return out;
}

getCorrHFem()

vec1& HCALResponse::getCorrHFem ( )

inline

Definition at line 47 of file HCALResponse.h.

References corrHFem.

Referenced by CalorimetryManager::updateHCAL().

{ return corrHFem; }

getCorrHFhad()

vec1& HCALResponse::getCorrHFhad ( )

inline

Definition at line 48 of file HCALResponse.h.

References corrHFhad.

Referenced by CalorimetryManager::updateHCAL().

{ return corrHFhad; }

getDet()

int HCALResponse::getDet ( int  ieta )

private

Definition at line 591 of file HCALResponse.cc.

References ztail::d, and LEDCalibrationChannels::ieta.

                                 {
  int d;
  for (d = 0; d < 2; d++) {
    if (ieta < HDeta[d + 1]) {
      break;
    }
  }
  return d;
}

getHCALEnergyResponse()

double HCALResponse::getHCALEnergyResponse	(	double	e,
		int	hit,
		RandomEngineAndDistribution const *	random
	)

Definition at line 563 of file HCALResponse.cc.

References c, MillePedeFileConverter_cfg::e, HcalResponse_cfi::eResponseCoefficient, HcalResponse_cfi::eResponseExponent, JetChargeProducer_cfi::exp, hcforward, dqm-mbProfile::log, dqmiodumpmetadata::n, AlCaHLTBitMon_ParallelJobs::p, L1TObjectsTimingClient_cff::resolution, alignCSCRings::s, mathSSE::sqrt(), and VFCAL.

Referenced by CalorimetryManager::HDShowerSimulation(), and CalorimetryManager::reconstructHCAL().

                                                                                                       {
  //response
  double s = eResponseScale[hit];
  double n = eResponseExponent;
  double p = eResponsePlateau[hit];
  double c = eResponseCoefficient;

  double response = e * p / (1 + c * exp(n * log(s / e)));

  if (response < 0.)
    response = 0.;

  //resolution
  double resolution;
  if (hit == hcforward)
    resolution =
        e * sqrt(RespPar[VFCAL][1][0] * RespPar[VFCAL][1][0] / e + RespPar[VFCAL][1][1] * RespPar[VFCAL][1][1]);
  else
    resolution =
        e * sqrt(RespPar[HCAL][hit][0] * RespPar[HCAL][hit][0] / (e) + RespPar[HCAL][hit][1] * RespPar[HCAL][hit][1]);

  //random smearing
  double rndm = gaussShootNoNegative(response, resolution, random);

  return rndm;
}

getMIPfraction()

double HCALResponse::getMIPfraction	(	double	energy,
		double	eta
	)

Definition at line 247 of file HCALResponse.cc.

References funct::abs(), HCALHighEnergyHPDFilter_cfi::energy, PVValHelper::eta, HcalResponse_cfi::etaStep, mps_fire::i, LEDCalibrationChannels::ieta, SiStripPI::mean, testProducerWithPsetDescEmpty_cfi::x1, testProducerWithPsetDescEmpty_cfi::x2, testProducerWithPsetDescEmpty_cfi::y1, and testProducerWithPsetDescEmpty_cfi::y2.

Referenced by CalorimetryManager::HDShowerSimulation().

                                                             {
  int ieta = abs((int)(eta / etaStep));
  int ie = -1;
  //check eta and det
  int det = getDet(ieta);
  int deta = ieta - HDeta[det];
  if (deta >= maxHDetas[det])
    deta = maxHDetas[det] - 1;
  else if (deta < 0)
    deta = 0;
  //find energy range
  for (int i = 0; i < maxHDe[det]; i++) {
    if (energy < eGridHD[det][i]) {
      if (i == 0)
        return mipfraction[det][0][deta];  // less than minimal - the first value is used instead of extrapolating
      else
        ie = i - 1;
      break;
    }
  }
  if (ie == -1)
    return mipfraction[det][maxHDe[det] - 1]
                      [deta];  // more than maximal - the last value is used instead of extrapolating
  double y1, y2;
  double x1 = eGridHD[det][ie];
  double x2 = eGridHD[det][ie + 1];
  y1 = mipfraction[det][ie][deta];
  y2 = mipfraction[det][ie + 1][deta];
  double mean = 0;
  mean = (y1 * (x2 - energy) + y2 * (energy - x1)) / (x2 - x1);
  return mean;
}

interEM()

double HCALResponse::interEM ( double  e, int  ie, int  ieta, RandomEngineAndDistribution const *  random  )

private

Definition at line 525 of file HCALResponse.cc.

References debug, MillePedeFileConverter_cfg::e, HcalResponse_cfi::eGridEM, LEDCalibrationChannels::ieta, SiStripPI::mean, HcalResponse_cfi::meanEM, HcalResponse_cfi::sigmaEM, testProducerWithPsetDescEmpty_cfi::x1, testProducerWithPsetDescEmpty_cfi::x2, testProducerWithPsetDescEmpty_cfi::y1, and testProducerWithPsetDescEmpty_cfi::y2.

                                                                                                  {
  double y1 = meanEM[ie][ieta];
  double y2 = meanEM[ie + 1][ieta];
  double x1 = eGridEM[ie];
  double x2 = eGridEM[ie + 1];

  if (debug) {
    //  cout.width(6);
    LogInfo("FastCalorimetry") << std::endl
                               << " HCALResponse::interEM mean " << std::endl
                               << " x, x1-x2, y1-y2 = " << e << ", " << x1 << "-" << x2 << " " << y1 << "-" << y2
                               << std::endl;
  }

  //linear interpolation
  double mean = (y1 * (x2 - e) + y2 * (e - x1)) / (x2 - x1);

  y1 = sigmaEM[ie][ieta];
  y2 = sigmaEM[ie + 1][ieta];

  if (debug) {
    //  cout.width(6);
    LogInfo("FastCalorimetry") << std::endl
                               << " HCALResponse::interEM sigma" << std::endl
                               << " x, x1-x2, y1-y2 = " << e << ", " << x1 << "-" << x2 << " " << y1 << "-" << y2
                               << std::endl;
  }

  //linear interpolation
  double sigma = (y1 * (x2 - e) + y2 * (e - x1)) / (x2 - x1);

  //random smearing
  double rndm = gaussShootNoNegative(mean, sigma, random);

  return rndm;
}

interHD()

double HCALResponse::interHD ( int  mip, double  e, int  ie, int  ieta, int  det, RandomEngineAndDistribution const *  random  )

private

Definition at line 449 of file HCALResponse.cc.

References MillePedeFileConverter_cfg::e, RandomEngineAndDistribution::flatShoot(), LEDCalibrationChannels::ieta, SiStripPI::mean, HcalResponse_cfi::nPar, AlCaHLTBitMon_ParallelJobs::p, RandomEngineAndDistribution::poissonShoot(), createJobs::tmp, testProducerWithPsetDescEmpty_cfi::x1, testProducerWithPsetDescEmpty_cfi::x2, testProducerWithPsetDescEmpty_cfi::y1, and testProducerWithPsetDescEmpty_cfi::y2.

                                                                                                                    {
  double x1, x2;
  double y1, y2;
  if (det == 2)
    mip = 2;  //ignore mip status for HF
  double mean = 0;
  vec1 pars(nPar, 0);

  // for ieta < 5 there is overlap between HE and HF, and measurement comes from HE
  if (det == 2 && ieta > 5 && e < 20) {
    for (int p = 0; p < 4; p++) {
      y1 = PoissonParameters[p][ie][ieta];
      y2 = PoissonParameters[p][ie + 1][ieta];
      if (e > 5) {
        x1 = eGridHD[det + 1][ie];
        x2 = eGridHD[det + 1][ie + 1];
        pars[p] = (y1 * (x2 - e) + y2 * (e - x1)) / (x2 - x1);
      } else
        pars[p] = y1;
    }
    mean =
        random->poissonShoot((int(PoissonShootNoNegative(pars[0], pars[1], random)) +
                              (int(PoissonShootNoNegative(pars[2], pars[3], random))) / 4 + random->flatShoot() / 4) *
                             6) /
        (0.3755 * 6);
  }

  else {
    x1 = eGridHD[det][ie];
    x2 = eGridHD[det][ie + 1];

    //calculate all parameters
    for (int p = 0; p < nPar; p++) {
      y1 = parameters[p][mip][det][ie][ieta];
      y2 = parameters[p][mip][det][ie + 1][ieta];

      //par-specific checks
      double custom = 0;
      bool use_custom = false;

      //do not let mu or sigma get extrapolated below zero for low energies
      //especially important for HF since extrapolation is used for E < 15 GeV
      if ((p == 0 || p == 1) && e < x1) {
        double tmp = (y1 * x2 - y2 * x1) / (x2 - x1);  //extrapolate down to e=0
        if (tmp < 0) {                                 //require mu,sigma > 0 for E > 0
          custom = y1 * e / x1;
          use_custom = true;
        }
      }
      //tail parameters have lower bounds - never extrapolate down
      else if ((p == 2 || p == 3 || p == 4 || p == 5)) {
        if (e < x1 && y1 < y2) {
          custom = y1;
          use_custom = true;
        } else if (e > x2 && y2 < y1) {
          custom = y2;
          use_custom = true;
        }
      }

      //linear interpolation
      if (use_custom)
        pars[p] = custom;
      else
        pars[p] = (y1 * (x2 - e) + y2 * (e - x1)) / (x2 - x1);
    }

    //random smearing
    if (nPar == 6)
      mean = cballShootNoNegative(pars[0], pars[1], pars[2], pars[3], pars[4], pars[5], random);
    else if (nPar == 2)
      mean = gaussShootNoNegative(pars[0], pars[1], random);  //gaussian fallback
  }
  return mean;
}

interMU()

double HCALResponse::interMU ( double  e, int  ie, int  ieta, RandomEngineAndDistribution const *  random  )

private

Definition at line 404 of file HCALResponse.cc.

References debug, MillePedeFileConverter_cfg::e, HcalResponse_cfi::eGridMU, RandomEngineAndDistribution::flatShoot(), mps_fire::i, LEDCalibrationChannels::ieta, HcalResponse_cfi::maxMUbin, SiStripPI::mean, HcalResponse_cfi::muStep, x, testProducerWithPsetDescEmpty_cfi::x1, testProducerWithPsetDescEmpty_cfi::x2, testProducerWithPsetDescEmpty_cfi::y1, and testProducerWithPsetDescEmpty_cfi::y2.

                                                                                                  {
  double x = random->flatShoot();

  int bin1 = maxMUbin;
  for (int i = 0; i < maxMUbin; i++) {
    if (x > responseMU[ie][ieta][i]) {
      bin1 = i - 1;
      break;
    }
  }
  int bin2 = maxMUbin;
  for (int i = 0; i < maxMUbin; i++) {
    if (x > responseMU[ie + 1][ieta][i]) {
      bin2 = i - 1;
      break;
    }
  }

  double x1 = eGridMU[ie];
  double x2 = eGridMU[ie + 1];
  double y1 = (bin1 + random->flatShoot()) * muStep;
  double y2 = (bin2 + random->flatShoot()) * muStep;

  if (debug) {
    //  cout.width(6);
    LogInfo("FastCalorimetry") << std::endl
                               << " HCALResponse::interMU  " << std::endl
                               << " x, x1-x2, y1-y2 = " << e << ", " << x1 << "-" << x2 << " " << y1 << "-" << y2
                               << std::endl;
  }

  //linear interpolation
  double mean = (y1 * (x2 - e) + y2 * (e - x1)) / (x2 - x1);

  if (debug) {
    //cout.width(6);
    LogInfo("FastCalorimetry") << std::endl
                               << " HCALResponse::interMU " << std::endl
                               << " e, ie, ieta = " << e << " " << ie << " " << ieta << std::endl
                               << " response  = " << mean << std::endl;
  }

  return mean;
}

PoissonShootNoNegative()

double HCALResponse::PoissonShootNoNegative ( double  e, double  sigma, RandomEngineAndDistribution const *  random  )

private

Definition at line 625 of file HCALResponse.cc.

References MillePedeFileConverter_cfg::e, MillePedeFileConverter_cfg::out, and RandomEngineAndDistribution::poissonShoot().

                                                                                                             {
  double out = -1;
  while (out < 0.) {
    out = random->poissonShoot(e);
    out = out + sigma;
  }
  return out;
}

responseHCAL()

double HCALResponse::responseHCAL	(	int	_mip,
		double	energy,
		double	eta,
		int	partype,
		RandomEngineAndDistribution const *	random
	)

Definition at line 280 of file HCALResponse.cc.

References funct::abs(), debug, HcalResponse_cfi::eGridEM, HcalResponse_cfi::eGridMU, HCALHighEnergyHPDFilter_cfi::energy, PVValHelper::eta, HcalResponse_cfi::etaGridMU, HcalResponse_cfi::etaStep, mps_fire::i, LEDCalibrationChannels::ieta, HcalResponse_cfi::maxEMe, HcalResponse_cfi::maxMUe, HcalResponse_cfi::maxMUeta, SiStripPI::mean, and HcalResponse_cfi::usemip.

Referenced by CalorimetryManager::HDShowerSimulation(), and CalorimetryManager::reconstructHCAL().

                                                                                                 {
  int ieta = abs((int)(eta / etaStep));
  int ie = -1;

  int mip;
  if (usemip)
    mip = _mip;
  else
    mip = 2;  //ignore mip, use only overall (mip + nomip) parameters

  double mean = 0;

  // e/gamma in HF
  if (partype == 0) {
    //check eta
    ieta -= HDeta[2];  // HF starts at ieta=30 till ieta=51
    // but resp.vector from index=0 through 20
    if (ieta >= maxEMeta)
      ieta = maxEMeta - 1;
    else if (ieta < 0)
      ieta = 0;

    //find energy range
    for (int i = 0; i < maxEMe; i++) {
      if (energy < eGridEM[i]) {
        if (i == 0)
          ie = 0;  // less than minimal - back extrapolation with the 1st interval
        else
          ie = i - 1;
        break;
      }
    }
    if (ie == -1)
      ie = maxEMe - 2;  // more than maximum - extrapolation with last interval

    //do smearing
    mean = interEM(energy, ie, ieta, random);
  }

  // hadrons
  else if (partype == 1) {
    //check eta and det
    int det = getDet(ieta);
    int deta = ieta - HDeta[det];
    if (deta >= maxHDetas[det])
      deta = maxHDetas[det] - 1;
    else if (deta < 0)
      deta = 0;

    //find energy range
    for (int i = 0; i < maxHDe[det]; i++) {
      if (energy < eGridHD[det][i]) {
        if (i == 0)
          ie = 0;  // less than minimal - back extrapolation with the 1st interval
        else
          ie = i - 1;
        break;
      }
    }
    if (ie == -1)
      ie = maxHDe[det] - 2;  // more than maximum - extrapolation with last interval

    //different energy smearing for low energy hadrons in HF
    if (det == 2 && energy < 20 && deta > 5) {
      for (int i = 0; i < maxHDe[3]; i++) {
        if (energy < eGridHD[3][i]) {
          if (i == 0)
            ie = 0;  // less than minimal - back extrapolation with the 1st interval
          else
            ie = i - 1;
          break;
        }
      }
    }
    //do smearing
    mean = interHD(mip, energy, ie, deta, det, random);
  }

  // muons
  else if (partype == 2) {
    //check eta
    ieta = maxMUeta;
    for (int i = 0; i < maxMUeta; i++) {
      if (fabs(eta) < etaGridMU[i]) {
        ieta = i;
        break;
      }
    }
    if (ieta < 0)
      ieta = 0;

    if (ieta < maxMUeta) {  // HB-HE
      //find energy range
      for (int i = 0; i < maxMUe; i++) {
        if (energy < eGridMU[i]) {
          if (i == 0)
            ie = 0;  // less than minimal - back extrapolation with the 1st interval
          else
            ie = i - 1;
          break;
        }
      }
      if (ie == -1)
        ie = maxMUe - 2;  // more than maximum - extrapolation with last interval

      //do smearing
      mean = interMU(energy, ie, ieta, random);
      if (mean > energy)
        mean = energy;
    }
  }

  // debugging
  if (debug) {
    //  cout.width(6);
    LogInfo("FastCalorimetry") << std::endl
                               << " HCALResponse::responseHCAL, partype = " << partype << " E, eta = " << energy << " "
                               << eta << "  mean = " << mean << std::endl;
  }

  return mean;
}

Member Data Documentation

barrelMUeta

int HCALResponse::barrelMUeta

private

Definition at line 87 of file HCALResponse.h.

cball

DoubleCrystalBallGenerator HCALResponse::cball

private

Definition at line 116 of file HCALResponse.h.

corrHFem

vec1 HCALResponse::corrHFem

private

Definition at line 123 of file HCALResponse.h.

Referenced by getCorrHFem().

corrHFgEm

vec2 HCALResponse::corrHFgEm

private

Definition at line 121 of file HCALResponse.h.

corrHFgHad

vec2 HCALResponse::corrHFgHad

private

Definition at line 121 of file HCALResponse.h.

corrHFhad

vec1 HCALResponse::corrHFhad

private

Definition at line 123 of file HCALResponse.h.

Referenced by getCorrHFhad().

corrHFhEm

vec2 HCALResponse::corrHFhEm

private

Definition at line 122 of file HCALResponse.h.

corrHFhHad

vec2 HCALResponse::corrHFhHad

private

Definition at line 122 of file HCALResponse.h.

debug

bool HCALResponse::debug

private

Definition at line 67 of file HCALResponse.h.

Referenced by rrapi.RRApi::dprint(), pkg.AbstractPkg::generate(), rrapi.RRApi::get(), pkg.AbstractPkg::get_kwds(), runTauIdMVA.TauIDEmbedder::loadMVA_WPs_run2_2017(), runTauIdMVA.TauIDEmbedder::runTauID(), and pkg.AbstractPkg::write().

eGridEM

vec1 HCALResponse::eGridEM

private

Definition at line 95 of file HCALResponse.h.

eGridHD

vec1 HCALResponse::eGridHD[4]

private

Definition at line 94 of file HCALResponse.h.

eGridMU

vec1 HCALResponse::eGridMU

private

Definition at line 96 of file HCALResponse.h.

endcapMUeta

int HCALResponse::endcapMUeta

private

Definition at line 87 of file HCALResponse.h.

energyHF

vec1 HCALResponse::energyHF

private

Definition at line 120 of file HCALResponse.h.

eResponseCoefficient

double HCALResponse::eResponseCoefficient

private

Definition at line 79 of file HCALResponse.h.

eResponseExponent

double HCALResponse::eResponseExponent

private

Definition at line 78 of file HCALResponse.h.

eResponsePlateau

double HCALResponse::eResponsePlateau[3]

private

Definition at line 77 of file HCALResponse.h.

eResponseScale

double HCALResponse::eResponseScale[3]

private

Definition at line 76 of file HCALResponse.h.

etaGridMU

vec1 HCALResponse::etaGridMU

private

Definition at line 97 of file HCALResponse.h.

etaStep

double HCALResponse::etaStep

private

Definition at line 85 of file HCALResponse.h.

HDeta

int HCALResponse::HDeta[4]

private

Definition at line 87 of file HCALResponse.h.

maxEMe

int HCALResponse::maxEMe

private

Definition at line 82 of file HCALResponse.h.

maxEMeta

int HCALResponse::maxEMeta

private

Definition at line 82 of file HCALResponse.h.

maxEne

int HCALResponse::maxEne

private

Definition at line 119 of file HCALResponse.h.

maxEta

int HCALResponse::maxEta

private

Definition at line 119 of file HCALResponse.h.

maxHDe

int HCALResponse::maxHDe[4]

private

Definition at line 83 of file HCALResponse.h.

maxHDetas

int HCALResponse::maxHDetas[3]

private

Definition at line 87 of file HCALResponse.h.

maxMUbin

int HCALResponse::maxMUbin

private

Definition at line 82 of file HCALResponse.h.

maxMUe

int HCALResponse::maxMUe

private

Definition at line 82 of file HCALResponse.h.

maxMUeta

int HCALResponse::maxMUeta

private

Definition at line 82 of file HCALResponse.h.

meanEM

vec2 HCALResponse::meanEM

private

Definition at line 107 of file HCALResponse.h.

mipfraction

vec3 HCALResponse::mipfraction

private

Definition at line 112 of file HCALResponse.h.

muStep

double HCALResponse::muStep

private

Definition at line 89 of file HCALResponse.h.

nPar

int HCALResponse::nPar

private

Definition at line 101 of file HCALResponse.h.

parameters

vec5 HCALResponse::parameters

private

Definition at line 103 of file HCALResponse.h.

parNames

std::vector<std::string> HCALResponse::parNames

private

Definition at line 102 of file HCALResponse.h.

PoissonParameters

vec3 HCALResponse::PoissonParameters

private

Definition at line 113 of file HCALResponse.h.

respFactorEM

double HCALResponse::respFactorEM

private

Definition at line 91 of file HCALResponse.h.

responseMU

vec3 HCALResponse::responseMU

private

Definition at line 111 of file HCALResponse.h.

RespPar

double HCALResponse::RespPar[3][2][3]

private

Definition at line 73 of file HCALResponse.h.

sigmaEM

vec2 HCALResponse::sigmaEM

private

Definition at line 107 of file HCALResponse.h.

usemip

bool HCALResponse::usemip

private

Definition at line 67 of file HCALResponse.h.