CMS 3D CMS Logo

List of all members | Public Member Functions | Private Member Functions | Private Attributes
HCALResponse Class Reference

#include <HCALResponse.h>

Public Member Functions

void correctHF (double e, int type)
 
vec1getCorrHFem ()
 
vec1getCorrHFhad ()
 
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_2024v14_cff::fraction, HCAL, mps_fire::i, dqmiolumiharvest::j, dqmdumpme::k, visualization-live-secondInstance_cfg::m, HcalResponse_cfi::maxEMe, HcalResponse_cfi::maxEne, razorScouting_cff::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.

18  {
19  //switches
20  debug = pset.getParameter<bool>("debug");
21  usemip = pset.getParameter<bool>("usemip");
22 
23  //values for "old" response parameterizations
24  //--------------------------------------------------------------------
25  RespPar[HCAL][0][0] = pset.getParameter<double>("HadronBarrelResolution_Stochastic");
26  RespPar[HCAL][0][1] = pset.getParameter<double>("HadronBarrelResolution_Constant");
27  RespPar[HCAL][0][2] = pset.getParameter<double>("HadronBarrelResolution_Noise");
28 
29  RespPar[HCAL][1][0] = pset.getParameter<double>("HadronEndcapResolution_Stochastic");
30  RespPar[HCAL][1][1] = pset.getParameter<double>("HadronEndcapResolution_Constant");
31  RespPar[HCAL][1][2] = pset.getParameter<double>("HadronEndcapResolution_Noise");
32 
33  RespPar[VFCAL][0][0] = pset.getParameter<double>("HadronForwardResolution_Stochastic");
34  RespPar[VFCAL][0][1] = pset.getParameter<double>("HadronForwardResolution_Constant");
35  RespPar[VFCAL][0][2] = pset.getParameter<double>("HadronForwardResolution_Noise");
36 
37  RespPar[VFCAL][1][0] = pset.getParameter<double>("ElectronForwardResolution_Stochastic");
38  RespPar[VFCAL][1][1] = pset.getParameter<double>("ElectronForwardResolution_Constant");
39  RespPar[VFCAL][1][2] = pset.getParameter<double>("ElectronForwardResolution_Noise");
40 
41  eResponseScale[0] = pset.getParameter<double>("eResponseScaleHB");
42  eResponseScale[1] = pset.getParameter<double>("eResponseScaleHE");
43  eResponseScale[2] = pset.getParameter<double>("eResponseScaleHF");
44 
45  eResponsePlateau[0] = pset.getParameter<double>("eResponsePlateauHB");
46  eResponsePlateau[1] = pset.getParameter<double>("eResponsePlateauHE");
47  eResponsePlateau[2] = pset.getParameter<double>("eResponsePlateauHF");
48 
49  eResponseExponent = pset.getParameter<double>("eResponseExponent");
50  eResponseCoefficient = pset.getParameter<double>("eResponseCoefficient");
51 
52  //pion parameters
53  //--------------------------------------------------------------------
54  //energy values
55  maxHDe[0] = pset.getParameter<int>("maxHBe");
56  maxHDe[1] = pset.getParameter<int>("maxHEe");
57  maxHDe[2] = pset.getParameter<int>("maxHFe");
58  maxHDe[3] = pset.getParameter<int>("maxHFlowe");
59 
60  eGridHD[0] = pset.getParameter<vec1>("eGridHB");
61  eGridHD[1] = pset.getParameter<vec1>("eGridHE");
62  eGridHD[2] = pset.getParameter<vec1>("eGridHF");
63  eGridHD[3] = pset.getParameter<vec1>("loweGridHF");
64 
65  //region eta indices calculated from eta values
66  etaStep = pset.getParameter<double>("etaStep");
67  //eta boundaries
68  HDeta[0] = abs((int)(pset.getParameter<double>("HBeta") / etaStep));
69  HDeta[1] = abs((int)(pset.getParameter<double>("HEeta") / etaStep));
70  HDeta[2] = abs((int)(pset.getParameter<double>("HFeta") / etaStep));
71  HDeta[3] = abs((int)(pset.getParameter<double>("maxHDeta") / etaStep)); //add 1 because this is the max index
72  //eta ranges
73  maxHDetas[0] = HDeta[1] - HDeta[0];
74  maxHDetas[1] = HDeta[2] - HDeta[1];
75  maxHDetas[2] = HDeta[3] - HDeta[2];
76 
77  //parameter info
78  nPar = pset.getParameter<int>("nPar");
79  parNames = pset.getParameter<std::vector<std::string> >("parNames");
80  std::string detNames[] = {"_HB", "_HE", "_HF"};
81  std::string mipNames[] = {"_mip", "_nomip", ""};
82  std::string fraction = "f";
83  //setup parameters (5D vector)
84  parameters = vec5(nPar, vec4(3, vec3(3)));
85  for (int p = 0; p < nPar; p++) { //loop over parameters
86  for (int m = 0; m < 3; m++) { //loop over mip, nomip, total
87  for (int d = 0; d < 3; d++) { //loop over dets: HB, HE, HF
88  //get from python
89  std::string pname = parNames[p] + detNames[d] + mipNames[m];
90  vec1 tmp = pset.getParameter<vec1>(pname);
91 
92  //resize vector for energy range of det d
93  parameters[p][m][d].resize(maxHDe[d]);
94 
95  for (int i = 0; i < maxHDe[d]; i++) { //loop over energy for det d
96  //resize vector for eta range of det d
97  parameters[p][m][d][i].resize(maxHDetas[d]);
98 
99  for (int j = 0; j < maxHDetas[d]; j++) { //loop over eta for det d
100  //fill in parameters vector from python
101  parameters[p][m][d][i][j] = tmp[i * maxHDetas[d] + j];
102  }
103  }
104  }
105  }
106  }
107  //set up Poisson parameters for low energy Hadrons in HF
108  //----------------------------------------------------------------------
109  PoissonParameters = vec3(4);
110  std::string PoissonParName[] = {"mean_overall", "shift_overall", "mean_between", "shift_between"};
111  for (int d = 0; d < 4; d++) { //loop over Poisson parameteres
112  vec1 tmp1 = pset.getParameter<vec1>(PoissonParName[d]);
113  for (int i = 0; i < maxHDe[3]; i++) { //loop over energy for low HF energy points
114  PoissonParameters[d].resize(maxHDe[3]);
115  for (int j = 0; j < maxHDetas[2]; j++) { //loop over HF eta points
116  PoissonParameters[d][i].resize(maxHDetas[2]);
117  PoissonParameters[d][i][j] = tmp1[i * maxHDetas[2] + j];
118  }
119  }
120  }
121 
122  //MIP fraction fill in 3d vector
124  mipfraction = vec3(3);
125  for (int d = 0; d < 3; d++) { //loop over dets: HB, HE, HF
126  //get from python
127  std::string mipname = fraction + mipNames[0] + detNames[d];
128  vec1 tmp1 = pset.getParameter<vec1>(mipname);
129  mipfraction[d].resize(maxHDe[d]);
130  for (int i = 0; i < maxHDe[d]; i++) { //loop over energy for det d
131  //resize vector for eta range of det d
132  mipfraction[d][i].resize(maxHDetas[d]);
133  for (int j = 0; j < maxHDetas[d]; j++) { //loop over eta for det d
134  //fill in parameters vector from python
135  mipfraction[d][i][j] = tmp1[i * maxHDetas[d] + j];
136  }
137  }
138  }
139 
140  // MUON probability histos for bin size = 0.25 GeV (0-10 GeV, 40 bins)
141  //--------------------------------------------------------------------
142  muStep = pset.getParameter<double>("muStep");
143  maxMUe = pset.getParameter<int>("maxMUe");
144  maxMUeta = pset.getParameter<int>("maxMUeta");
145  maxMUbin = pset.getParameter<int>("maxMUbin");
146  eGridMU = pset.getParameter<vec1>("eGridMU");
147  etaGridMU = pset.getParameter<vec1>("etaGridMU");
148  vec1 _responseMU[2] = {pset.getParameter<vec1>("responseMUBarrel"), pset.getParameter<vec1>("responseMUEndcap")};
149 
150  //get muon region eta indices from the eta grid
151  double _barrelMUeta = pset.getParameter<double>("barrelMUeta");
152  double _endcapMUeta = pset.getParameter<double>("endcapMUeta");
154  for (int i = 0; i < maxMUeta; i++) {
155  if (fabs(_barrelMUeta) <= etaGridMU[i]) {
156  barrelMUeta = i;
157  break;
158  }
159  }
160  for (int i = 0; i < maxMUeta; i++) {
161  if (fabs(_endcapMUeta) <= etaGridMU[i]) {
162  endcapMUeta = i;
163  break;
164  }
165  }
166  int maxMUetas[] = {endcapMUeta - barrelMUeta, maxMUeta - endcapMUeta};
167 
168  //initialize 3D vector
170 
171  //fill in 3D vector
172  //(complementary cumulative distribution functions, from normalized response distributions)
173  int loc, eta_loc;
174  loc = eta_loc = -1;
175  for (int i = 0; i < maxMUe; i++) {
176  for (int j = 0; j < maxMUeta; j++) {
177  //check location - barrel, endcap, or forward
178  if (j == barrelMUeta) {
179  loc = 0;
180  eta_loc = barrelMUeta;
181  } else if (j == endcapMUeta) {
182  loc = 1;
183  eta_loc = endcapMUeta;
184  }
185 
186  for (int k = 0; k < maxMUbin; k++) {
187  responseMU[i][j][k] = _responseMU[loc][i * maxMUetas[loc] * maxMUbin + (j - eta_loc) * maxMUbin + k];
188 
189  if (debug) {
190  //cout.width(6);
191  LogInfo("FastCalorimetry") << " responseMU " << i << " " << j << " " << k << " = " << responseMU[i][j][k]
192  << std::endl;
193  }
194  }
195  }
196  }
197 
198  // values for EM response in HF
199  //--------------------------------------------------------------------
200  maxEMe = pset.getParameter<int>("maxEMe");
201  maxEMeta = maxHDetas[2];
202  respFactorEM = pset.getParameter<double>("respFactorEM");
203  eGridEM = pset.getParameter<vec1>("eGridEM");
204 
205  // e-gamma mean response and sigma in HF
206  vec1 _meanEM = pset.getParameter<vec1>("meanEM");
207  vec1 _sigmaEM = pset.getParameter<vec1>("sigmaEM");
208 
209  //fill in 2D vectors (w/ correction factor applied)
210  meanEM = vec2(maxEMe, vec1(maxEMeta, 0));
211  sigmaEM = vec2(maxEMe, vec1(maxEMeta, 0));
212  for (int i = 0; i < maxEMe; i++) {
213  for (int j = 0; j < maxEMeta; j++) {
214  meanEM[i][j] = respFactorEM * _meanEM[i * maxEMeta + j];
215  sigmaEM[i][j] = respFactorEM * _sigmaEM[i * maxEMeta + j];
216  }
217  }
218 
219  // HF correction for SL
220  //---------------------
221  maxEta = pset.getParameter<int>("maxEta");
222  maxEne = pset.getParameter<int>("maxEne");
223  energyHF = pset.getParameter<vec1>("energyHF");
224  vec1 _corrHFgEm = pset.getParameter<vec1>("corrHFgEm");
225  vec1 _corrHFgHad = pset.getParameter<vec1>("corrHFgHad");
226  vec1 _corrHFhEm = pset.getParameter<vec1>("corrHFhEm");
227  vec1 _corrHFhHad = pset.getParameter<vec1>("corrHFhHad");
228  corrHFem = vec1(maxEta, 0);
229  corrHFhad = vec1(maxEta, 0);
230 
231  // initialize 2D vector corrHFgEm[energy][eta]
232  corrHFgEm = vec2(maxEne, vec1(maxEta, 0));
233  corrHFgHad = vec2(maxEne, vec1(maxEta, 0));
234  corrHFhEm = vec2(maxEne, vec1(maxEta, 0));
235  corrHFhHad = vec2(maxEne, vec1(maxEta, 0));
236  // Fill
237  for (int i = 0; i < maxEne; i++) {
238  for (int j = 0; j < maxEta; j++) {
239  corrHFgEm[i][j] = _corrHFgEm[i * maxEta + j];
240  corrHFgHad[i][j] = _corrHFgHad[i * maxEta + j];
241  corrHFhEm[i][j] = _corrHFhEm[i * maxEta + j];
242  corrHFhHad[i][j] = _corrHFhHad[i * maxEta + j];
243  }
244  }
245 }
double etaStep
Definition: HCALResponse.h:85
double RespPar[3][2][3]
Definition: HCALResponse.h:73
int maxHDe[4]
Definition: HCALResponse.h:83
int HDeta[4]
Definition: HCALResponse.h:87
double eResponseCoefficient
Definition: HCALResponse.h:79
double eResponseScale[3]
Definition: HCALResponse.h:76
std::vector< std::string > parNames
Definition: HCALResponse.h:102
double eResponseExponent
Definition: HCALResponse.h:78
std::vector< vec4 > vec5
Definition: HCALResponse.h:19
std::vector< vec2 > vec3
Definition: HCALResponse.h:17
vec3 PoissonParameters
Definition: HCALResponse.h:113
std::vector< double > vec1
Definition: HCALResponse.h:15
double muStep
Definition: HCALResponse.h:89
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
std::vector< vec3 > vec4
Definition: HCALResponse.h:18
d
Definition: ztail.py:151
Log< level::Info, false > LogInfo
vec1 eGridHD[4]
Definition: HCALResponse.h:94
double respFactorEM
Definition: HCALResponse.h:91
int maxHDetas[3]
Definition: HCALResponse.h:87
tmp
align.sh
Definition: createJobs.py:716
std::vector< vec1 > vec2
Definition: HCALResponse.h:16
double eResponsePlateau[3]
Definition: HCALResponse.h:77

◆ ~HCALResponse()

HCALResponse::~HCALResponse ( )
inline

Definition at line 32 of file HCALResponse.h.

32 {}

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.

615  {
616  double out = cball.shoot(mu, sigma, aL, nL, aR, nR, random);
617  if (mu >= 0.) {
618  while (out < 0.)
619  out = cball.shoot(mu, sigma, aL, nL, aR, nR, random);
620  }
621  //else give up on re-trying, otherwise too much time can be lost before emeas comes out positive
622 
623  return out;
624 }
double shoot(double mu, double sigma, double aL, double nL, double aR, double nR, RandomEngineAndDistribution const *random)
DoubleCrystalBallGenerator cball
Definition: HCALResponse.h:116

◆ 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, and razorScouting_cff::maxEta.

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

634  {
635  int jmin = 0;
636  for (int i = 0; i < maxEne; i++) {
637  if (ee >= energyHF[i])
638  jmin = i;
639  }
640 
641  double x1, x2;
642  double y1em, y2em;
643  double y1had, y2had;
644  for (int i = 0; i < maxEta; ++i) {
645  if (ee < energyHF[0]) {
646  if (abs(type) == 11 || abs(type) == 22) {
647  corrHFem[i] = corrHFgEm[0][i];
648  corrHFhad[i] = corrHFgHad[0][i];
649  } else {
650  corrHFem[i] = corrHFhEm[0][i];
651  corrHFhad[i] = corrHFhHad[0][i];
652  }
653  } else if (jmin >= maxEne - 1) {
654  if (abs(type) == 11 || abs(type) == 22) {
655  corrHFem[i] = corrHFgEm[maxEta][i];
657  } else {
658  corrHFem[i] = corrHFhEm[maxEta][i];
660  }
661  } else {
662  x1 = energyHF[jmin];
663  x2 = energyHF[jmin + 1];
664  if (abs(type) == 11 || abs(type) == 22) {
665  y1em = corrHFgEm[jmin][i];
666  y2em = corrHFgEm[jmin + 1][i];
667  y1had = corrHFgHad[jmin][i];
668  y2had = corrHFgHad[jmin + 1][i];
669  } else {
670  y1em = corrHFhEm[jmin][i];
671  y2em = corrHFhEm[jmin + 1][i];
672  y1had = corrHFhHad[jmin][i];
673  y2had = corrHFhHad[jmin + 1][i];
674  }
675  corrHFem[i] = y1em + (ee - x1) * ((y2em - y1em) / (x2 - x1));
676  corrHFhad[i] = y1had + (ee - x1) * ((y2had - y1had) / (x2 - x1));
677  }
678  }
679 }
Abs< T >::type abs(const T &t)
Definition: Abs.h:22

◆ 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.

602  {
603  double out = random->gaussShoot(e, sigma);
604  if (e >= 0.) {
605  while (out < 0.)
606  out = random->gaussShoot(e, sigma);
607  }
608  //else give up on re-trying, otherwise too much time can be lost before emeas comes out positive
609 
610  return out;
611 }

◆ getCorrHFem()

vec1& HCALResponse::getCorrHFem ( )
inline

Definition at line 47 of file HCALResponse.h.

References corrHFem.

Referenced by CalorimetryManager::updateHCAL().

47 { return corrHFem; }

◆ getCorrHFhad()

vec1& HCALResponse::getCorrHFhad ( )
inline

Definition at line 48 of file HCALResponse.h.

References corrHFhad.

Referenced by CalorimetryManager::updateHCAL().

48 { return corrHFhad; }

◆ getDet()

int HCALResponse::getDet ( int  ieta)
private

Definition at line 591 of file HCALResponse.cc.

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

591  {
592  int d;
593  for (d = 0; d < 2; d++) {
594  if (ieta < HDeta[d + 1]) {
595  break;
596  }
597  }
598  return d;
599 }
int HDeta[4]
Definition: HCALResponse.h:87
d
Definition: ztail.py:151

◆ getHCALEnergyResponse()

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

Definition at line 563 of file HCALResponse.cc.

References DummyCfis::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, testHGCalSimSingleMuonPt100_cfg::rndm, alignCSCRings::s, mathSSE::sqrt(), and VFCAL.

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

563  {
564  //response
565  double s = eResponseScale[hit];
566  double n = eResponseExponent;
567  double p = eResponsePlateau[hit];
568  double c = eResponseCoefficient;
569 
570  double response = e * p / (1 + c * exp(n * log(s / e)));
571 
572  if (response < 0.)
573  response = 0.;
574 
575  //resolution
576  double resolution;
577  if (hit == hcforward)
578  resolution =
579  e * sqrt(RespPar[VFCAL][1][0] * RespPar[VFCAL][1][0] / e + RespPar[VFCAL][1][1] * RespPar[VFCAL][1][1]);
580  else
581  resolution =
582  e * sqrt(RespPar[HCAL][hit][0] * RespPar[HCAL][hit][0] / (e) + RespPar[HCAL][hit][1] * RespPar[HCAL][hit][1]);
583 
584  //random smearing
585  double rndm = gaussShootNoNegative(response, resolution, random);
586 
587  return rndm;
588 }
double RespPar[3][2][3]
Definition: HCALResponse.h:73
double eResponseCoefficient
Definition: HCALResponse.h:79
double eResponseScale[3]
Definition: HCALResponse.h:76
double eResponseExponent
Definition: HCALResponse.h:78
T sqrt(T t)
Definition: SSEVec.h:23
Definition: HCAL.py:1
double eResponsePlateau[3]
Definition: HCALResponse.h:77
double gaussShootNoNegative(double e, double sigma, RandomEngineAndDistribution const *)

◆ getMIPfraction()

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

Definition at line 247 of file HCALResponse.cc.

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

Referenced by CalorimetryManager::HDShowerSimulation().

247  {
248  int ieta = abs((int)(eta / etaStep));
249  int ie = -1;
250  //check eta and det
251  int det = getDet(ieta);
252  int deta = ieta - HDeta[det];
253  if (deta >= maxHDetas[det])
254  deta = maxHDetas[det] - 1;
255  else if (deta < 0)
256  deta = 0;
257  //find energy range
258  for (int i = 0; i < maxHDe[det]; i++) {
259  if (energy < eGridHD[det][i]) {
260  if (i == 0)
261  return mipfraction[det][0][deta]; // less than minimal - the first value is used instead of extrapolating
262  else
263  ie = i - 1;
264  break;
265  }
266  }
267  if (ie == -1)
268  return mipfraction[det][maxHDe[det] - 1]
269  [deta]; // more than maximal - the last value is used instead of extrapolating
270  double y1, y2;
271  double x1 = eGridHD[det][ie];
272  double x2 = eGridHD[det][ie + 1];
273  y1 = mipfraction[det][ie][deta];
274  y2 = mipfraction[det][ie + 1][deta];
275  double mean = 0;
276  mean = (y1 * (x2 - energy) + y2 * (energy - x1)) / (x2 - x1);
277  return mean;
278 }
double etaStep
Definition: HCALResponse.h:85
int maxHDe[4]
Definition: HCALResponse.h:83
int HDeta[4]
Definition: HCALResponse.h:87
int getDet(int ieta)
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
vec1 eGridHD[4]
Definition: HCALResponse.h:94
int maxHDetas[3]
Definition: HCALResponse.h:87

◆ 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, hcalRecHitTable_cff::ieta, SiStripPI::mean, HcalResponse_cfi::meanEM, testHGCalSimSingleMuonPt100_cfg::rndm, HcalResponse_cfi::sigmaEM, testProducerWithPsetDescEmpty_cfi::y1, and testProducerWithPsetDescEmpty_cfi::y2.

525  {
526  double y1 = meanEM[ie][ieta];
527  double y2 = meanEM[ie + 1][ieta];
528  double x1 = eGridEM[ie];
529  double x2 = eGridEM[ie + 1];
530 
531  if (debug) {
532  // cout.width(6);
533  LogInfo("FastCalorimetry") << std::endl
534  << " HCALResponse::interEM mean " << std::endl
535  << " x, x1-x2, y1-y2 = " << e << ", " << x1 << "-" << x2 << " " << y1 << "-" << y2
536  << std::endl;
537  }
538 
539  //linear interpolation
540  double mean = (y1 * (x2 - e) + y2 * (e - x1)) / (x2 - x1);
541 
542  y1 = sigmaEM[ie][ieta];
543  y2 = sigmaEM[ie + 1][ieta];
544 
545  if (debug) {
546  // cout.width(6);
547  LogInfo("FastCalorimetry") << std::endl
548  << " HCALResponse::interEM sigma" << std::endl
549  << " x, x1-x2, y1-y2 = " << e << ", " << x1 << "-" << x2 << " " << y1 << "-" << y2
550  << std::endl;
551  }
552 
553  //linear interpolation
554  double sigma = (y1 * (x2 - e) + y2 * (e - x1)) / (x2 - x1);
555 
556  //random smearing
557  double rndm = gaussShootNoNegative(mean, sigma, random);
558 
559  return rndm;
560 }
Log< level::Info, false > LogInfo
double gaussShootNoNegative(double e, double sigma, RandomEngineAndDistribution const *)

◆ 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(), hcalRecHitTable_cff::ieta, SiStripPI::mean, HcalResponse_cfi::nPar, AlCaHLTBitMon_ParallelJobs::p, RandomEngineAndDistribution::poissonShoot(), createJobs::tmp, testProducerWithPsetDescEmpty_cfi::y1, and testProducerWithPsetDescEmpty_cfi::y2.

449  {
450  double x1, x2;
451  double y1, y2;
452  if (det == 2)
453  mip = 2; //ignore mip status for HF
454  double mean = 0;
455  vec1 pars(nPar, 0);
456 
457  // for ieta < 5 there is overlap between HE and HF, and measurement comes from HE
458  if (det == 2 && ieta > 5 && e < 20) {
459  for (int p = 0; p < 4; p++) {
460  y1 = PoissonParameters[p][ie][ieta];
461  y2 = PoissonParameters[p][ie + 1][ieta];
462  if (e > 5) {
463  x1 = eGridHD[det + 1][ie];
464  x2 = eGridHD[det + 1][ie + 1];
465  pars[p] = (y1 * (x2 - e) + y2 * (e - x1)) / (x2 - x1);
466  } else
467  pars[p] = y1;
468  }
469  mean =
470  random->poissonShoot((int(PoissonShootNoNegative(pars[0], pars[1], random)) +
471  (int(PoissonShootNoNegative(pars[2], pars[3], random))) / 4 + random->flatShoot() / 4) *
472  6) /
473  (0.3755 * 6);
474  }
475 
476  else {
477  x1 = eGridHD[det][ie];
478  x2 = eGridHD[det][ie + 1];
479 
480  //calculate all parameters
481  for (int p = 0; p < nPar; p++) {
482  y1 = parameters[p][mip][det][ie][ieta];
483  y2 = parameters[p][mip][det][ie + 1][ieta];
484 
485  //par-specific checks
486  double custom = 0;
487  bool use_custom = false;
488 
489  //do not let mu or sigma get extrapolated below zero for low energies
490  //especially important for HF since extrapolation is used for E < 15 GeV
491  if ((p == 0 || p == 1) && e < x1) {
492  double tmp = (y1 * x2 - y2 * x1) / (x2 - x1); //extrapolate down to e=0
493  if (tmp < 0) { //require mu,sigma > 0 for E > 0
494  custom = y1 * e / x1;
495  use_custom = true;
496  }
497  }
498  //tail parameters have lower bounds - never extrapolate down
499  else if ((p == 2 || p == 3 || p == 4 || p == 5)) {
500  if (e < x1 && y1 < y2) {
501  custom = y1;
502  use_custom = true;
503  } else if (e > x2 && y2 < y1) {
504  custom = y2;
505  use_custom = true;
506  }
507  }
508 
509  //linear interpolation
510  if (use_custom)
511  pars[p] = custom;
512  else
513  pars[p] = (y1 * (x2 - e) + y2 * (e - x1)) / (x2 - x1);
514  }
515 
516  //random smearing
517  if (nPar == 6)
518  mean = cballShootNoNegative(pars[0], pars[1], pars[2], pars[3], pars[4], pars[5], random);
519  else if (nPar == 2)
520  mean = gaussShootNoNegative(pars[0], pars[1], random); //gaussian fallback
521  }
522  return mean;
523 }
double PoissonShootNoNegative(double e, double sigma, RandomEngineAndDistribution const *)
vec3 PoissonParameters
Definition: HCALResponse.h:113
std::vector< double > vec1
Definition: HCALResponse.h:15
vec1 eGridHD[4]
Definition: HCALResponse.h:94
tmp
align.sh
Definition: createJobs.py:716
double gaussShootNoNegative(double e, double sigma, RandomEngineAndDistribution const *)
double cballShootNoNegative(double mu, double sigma, double aL, double nL, double aR, double nR, RandomEngineAndDistribution const *)

◆ 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, hcalRecHitTable_cff::ieta, HcalResponse_cfi::maxMUbin, SiStripPI::mean, HcalResponse_cfi::muStep, x, testProducerWithPsetDescEmpty_cfi::y1, and testProducerWithPsetDescEmpty_cfi::y2.

404  {
405  double x = random->flatShoot();
406 
407  int bin1 = maxMUbin;
408  for (int i = 0; i < maxMUbin; i++) {
409  if (x > responseMU[ie][ieta][i]) {
410  bin1 = i - 1;
411  break;
412  }
413  }
414  int bin2 = maxMUbin;
415  for (int i = 0; i < maxMUbin; i++) {
416  if (x > responseMU[ie + 1][ieta][i]) {
417  bin2 = i - 1;
418  break;
419  }
420  }
421 
422  double x1 = eGridMU[ie];
423  double x2 = eGridMU[ie + 1];
424  double y1 = (bin1 + random->flatShoot()) * muStep;
425  double y2 = (bin2 + random->flatShoot()) * muStep;
426 
427  if (debug) {
428  // cout.width(6);
429  LogInfo("FastCalorimetry") << std::endl
430  << " HCALResponse::interMU " << std::endl
431  << " x, x1-x2, y1-y2 = " << e << ", " << x1 << "-" << x2 << " " << y1 << "-" << y2
432  << std::endl;
433  }
434 
435  //linear interpolation
436  double mean = (y1 * (x2 - e) + y2 * (e - x1)) / (x2 - x1);
437 
438  if (debug) {
439  //cout.width(6);
440  LogInfo("FastCalorimetry") << std::endl
441  << " HCALResponse::interMU " << std::endl
442  << " e, ie, ieta = " << e << " " << ie << " " << ieta << std::endl
443  << " response = " << mean << std::endl;
444  }
445 
446  return mean;
447 }
double muStep
Definition: HCALResponse.h:89
Log< level::Info, false > LogInfo

◆ 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().

625  {
626  double out = -1;
627  while (out < 0.) {
628  out = random->poissonShoot(e);
629  out = out + sigma;
630  }
631  return out;
632 }

◆ 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, hcalRecHitTable_cff::energy, PVValHelper::eta, HcalResponse_cfi::etaGridMU, HcalResponse_cfi::etaStep, mps_fire::i, hcalRecHitTable_cff::ieta, HcalResponse_cfi::maxEMe, HcalResponse_cfi::maxMUe, HcalResponse_cfi::maxMUeta, SiStripPI::mean, and HcalResponse_cfi::usemip.

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

281  {
282  int ieta = abs((int)(eta / etaStep));
283  int ie = -1;
284 
285  int mip;
286  if (usemip)
287  mip = _mip;
288  else
289  mip = 2; //ignore mip, use only overall (mip + nomip) parameters
290 
291  double mean = 0;
292 
293  // e/gamma in HF
294  if (partype == 0) {
295  //check eta
296  ieta -= HDeta[2]; // HF starts at ieta=30 till ieta=51
297  // but resp.vector from index=0 through 20
298  if (ieta >= maxEMeta)
299  ieta = maxEMeta - 1;
300  else if (ieta < 0)
301  ieta = 0;
302 
303  //find energy range
304  for (int i = 0; i < maxEMe; i++) {
305  if (energy < eGridEM[i]) {
306  if (i == 0)
307  ie = 0; // less than minimal - back extrapolation with the 1st interval
308  else
309  ie = i - 1;
310  break;
311  }
312  }
313  if (ie == -1)
314  ie = maxEMe - 2; // more than maximum - extrapolation with last interval
315 
316  //do smearing
317  mean = interEM(energy, ie, ieta, random);
318  }
319 
320  // hadrons
321  else if (partype == 1) {
322  //check eta and det
323  int det = getDet(ieta);
324  int deta = ieta - HDeta[det];
325  if (deta >= maxHDetas[det])
326  deta = maxHDetas[det] - 1;
327  else if (deta < 0)
328  deta = 0;
329 
330  //find energy range
331  for (int i = 0; i < maxHDe[det]; i++) {
332  if (energy < eGridHD[det][i]) {
333  if (i == 0)
334  ie = 0; // less than minimal - back extrapolation with the 1st interval
335  else
336  ie = i - 1;
337  break;
338  }
339  }
340  if (ie == -1)
341  ie = maxHDe[det] - 2; // more than maximum - extrapolation with last interval
342 
343  //different energy smearing for low energy hadrons in HF
344  if (det == 2 && energy < 20 && deta > 5) {
345  for (int i = 0; i < maxHDe[3]; i++) {
346  if (energy < eGridHD[3][i]) {
347  if (i == 0)
348  ie = 0; // less than minimal - back extrapolation with the 1st interval
349  else
350  ie = i - 1;
351  break;
352  }
353  }
354  }
355  //do smearing
356  mean = interHD(mip, energy, ie, deta, det, random);
357  }
358 
359  // muons
360  else if (partype == 2) {
361  //check eta
362  ieta = maxMUeta;
363  for (int i = 0; i < maxMUeta; i++) {
364  if (fabs(eta) < etaGridMU[i]) {
365  ieta = i;
366  break;
367  }
368  }
369  if (ieta < 0)
370  ieta = 0;
371 
372  if (ieta < maxMUeta) { // HB-HE
373  //find energy range
374  for (int i = 0; i < maxMUe; i++) {
375  if (energy < eGridMU[i]) {
376  if (i == 0)
377  ie = 0; // less than minimal - back extrapolation with the 1st interval
378  else
379  ie = i - 1;
380  break;
381  }
382  }
383  if (ie == -1)
384  ie = maxMUe - 2; // more than maximum - extrapolation with last interval
385 
386  //do smearing
387  mean = interMU(energy, ie, ieta, random);
388  if (mean > energy)
389  mean = energy;
390  }
391  }
392 
393  // debugging
394  if (debug) {
395  // cout.width(6);
396  LogInfo("FastCalorimetry") << std::endl
397  << " HCALResponse::responseHCAL, partype = " << partype << " E, eta = " << energy << " "
398  << eta << " mean = " << mean << std::endl;
399  }
400 
401  return mean;
402 }
double etaStep
Definition: HCALResponse.h:85
int maxHDe[4]
Definition: HCALResponse.h:83
int HDeta[4]
Definition: HCALResponse.h:87
int getDet(int ieta)
double interEM(double e, int ie, int ieta, RandomEngineAndDistribution const *)
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
double interHD(int mip, double e, int ie, int ieta, int det, RandomEngineAndDistribution const *)
Log< level::Info, false > LogInfo
double interMU(double e, int ie, int ieta, RandomEngineAndDistribution const *)
vec1 eGridHD[4]
Definition: HCALResponse.h:94
int maxHDetas[3]
Definition: HCALResponse.h:87

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

◆ 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.