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)
 
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 ( const edm::ParameterSet &  pset)

Definition at line 18 of file HCALResponse.cc.

References funct::abs(), edmIntegrityCheck::d, debug, dedxEstimators_cff::fraction, edm::ParameterSet::getParameter(), HCAL, mps_fire::i, gen::k, create_public_lumi_plots::loc, funct::m, maxEta, AlCaHLTBitMon_ParallelJobs::p, metProducer_cfi::parameters, BPhysicsValidation_cfi::pname, AlCaHLTBitMon_QueryRunRegistry::string, tmp, 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 
123  //MIP fraction fill in 3d vector
125  mipfraction = vec3(3);
126  for(int d = 0; d < 3; d++){ //loop over dets: HB, HE, HF
127  //get from python
128  std::string mipname = fraction + mipNames[0] + detNames[d] ;
129  vec1 tmp1 = pset.getParameter<vec1>(mipname);
130  mipfraction[d].resize(maxHDe[d]);
131  for(int i = 0; i < maxHDe[d]; i++){ //loop over energy for det d
132  //resize vector for eta range of det d
133  mipfraction[d][i].resize(maxHDetas[d]);
134  for(int j = 0; j < maxHDetas[d]; j++){ //loop over eta for det d
135  //fill in parameters vector from python
136  mipfraction[d][i][j]= tmp1[i*maxHDetas[d] + j];
137  }
138  }
139  }
140 
141  // MUON probability histos for bin size = 0.25 GeV (0-10 GeV, 40 bins)
142  //--------------------------------------------------------------------
143  muStep = pset.getParameter<double>("muStep");
144  maxMUe = pset.getParameter<int>("maxMUe");
145  maxMUeta = pset.getParameter<int>("maxMUeta");
146  maxMUbin = pset.getParameter<int>("maxMUbin");
147  eGridMU = pset.getParameter<vec1>("eGridMU");
148  etaGridMU = pset.getParameter<vec1>("etaGridMU");
149  vec1 _responseMU[2] = {pset.getParameter<vec1>("responseMUBarrel"),pset.getParameter<vec1>("responseMUEndcap")};
150 
151  //get muon region eta indices from the eta grid
152  double _barrelMUeta = pset.getParameter<double>("barrelMUeta");
153  double _endcapMUeta = pset.getParameter<double>("endcapMUeta");
154  barrelMUeta = endcapMUeta = maxMUeta;
155  for(int i = 0; i < maxMUeta; i++) {
156  if(fabs(_barrelMUeta) <= etaGridMU[i]) { barrelMUeta = i; break; }
157  }
158  for(int i = 0; i < maxMUeta; i++) {
159  if(fabs(_endcapMUeta) <= etaGridMU[i]) { endcapMUeta = i; break; }
160  }
161  int maxMUetas[] = {endcapMUeta - barrelMUeta, maxMUeta - endcapMUeta};
162 
163  //initialize 3D vector
164  responseMU = vec3(maxMUe,vec2(maxMUeta,vec1(maxMUbin,0)));
165 
166  //fill in 3D vector
167  //(complementary cumulative distribution functions, from normalized response distributions)
168  int loc, eta_loc;
169  loc = eta_loc = -1;
170  for(int i = 0; i < maxMUe; i++){
171  for(int j = 0; j < maxMUeta; j++){
172  //check location - barrel, endcap, or forward
173  if(j==barrelMUeta) {loc = 0; eta_loc = barrelMUeta;}
174  else if(j==endcapMUeta) {loc = 1; eta_loc = endcapMUeta;}
175 
176  for(int k = 0; k < maxMUbin; k++){
177  responseMU[i][j][k] = _responseMU[loc][i*maxMUetas[loc]*maxMUbin + (j-eta_loc)*maxMUbin + k];
178 
179  if(debug) {
180  //cout.width(6);
181  LogInfo("FastCalorimetry") << " responseMU " << i << " " << j << " " << k << " = "
182  << responseMU[i][j][k] << std::endl;
183  }
184  }
185  }
186  }
187 
188  // values for EM response in HF
189  //--------------------------------------------------------------------
190  maxEMe = pset.getParameter<int>("maxEMe");
191  maxEMeta = maxHDetas[2];
192  respFactorEM = pset.getParameter<double>("respFactorEM");
193  eGridEM = pset.getParameter<vec1>("eGridEM");
194 
195  // e-gamma mean response and sigma in HF
196  vec1 _meanEM = pset.getParameter<vec1>("meanEM");
197  vec1 _sigmaEM = pset.getParameter<vec1>("sigmaEM");
198 
199  //fill in 2D vectors (w/ correction factor applied)
200  meanEM = vec2(maxEMe,vec1(maxEMeta,0));
201  sigmaEM = vec2(maxEMe,vec1(maxEMeta,0));
202  for(int i = 0; i < maxEMe; i++){
203  for(int j = 0; j < maxEMeta; j++){
204  meanEM[i][j] = respFactorEM * _meanEM[i*maxEMeta + j];
205  sigmaEM[i][j] = respFactorEM * _sigmaEM[i*maxEMeta + j];
206  }
207  }
208 
209  // HF correction for SL
210  //---------------------
211  maxEta = pset.getParameter<int>("maxEta");
212  maxEne = pset.getParameter<int>("maxEne");
213  energyHF = pset.getParameter<vec1>("energyHF");
214  vec1 _corrHFgEm = pset.getParameter<vec1>("corrHFgEm");
215  vec1 _corrHFgHad = pset.getParameter<vec1>("corrHFgHad");
216  vec1 _corrHFhEm = pset.getParameter<vec1>("corrHFhEm");
217  vec1 _corrHFhHad = pset.getParameter<vec1>("corrHFhHad");
218  corrHFem = vec1(maxEta,0);
219  corrHFhad = vec1(maxEta,0);
220 
221 // initialize 2D vector corrHFgEm[energy][eta]
222  corrHFgEm = vec2(maxEne,vec1(maxEta,0));
223  corrHFgHad = vec2(maxEne,vec1(maxEta,0));
224  corrHFhEm = vec2(maxEne,vec1(maxEta,0));
225  corrHFhHad = vec2(maxEne,vec1(maxEta,0));
226 // Fill
227  for(int i = 0; i < maxEne; i++){
228  for(int j = 0; j < maxEta; j++){
229  corrHFgEm[i][j] = _corrHFgEm[i*maxEta + j];
230  corrHFgHad[i][j] = _corrHFgHad[i*maxEta + j];
231  corrHFhEm[i][j] = _corrHFhEm[i*maxEta + j];
232  corrHFhHad[i][j] = _corrHFhHad[i*maxEta + j];
233  }
234  }
235 }
HCALResponse::etaStep
double etaStep
Definition: HCALResponse.h:86
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
HCALResponse::RespPar
double RespPar[3][2][3]
Definition: HCALResponse.h:74
mps_fire.i
i
Definition: mps_fire.py:156
HCALResponse::maxHDe
int maxHDe[4]
Definition: HCALResponse.h:84
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
create_public_lumi_plots.loc
loc
Definition: create_public_lumi_plots.py:1108
HCALResponse::HDeta
int HDeta[4]
Definition: HCALResponse.h:88
HCALResponse::eResponseCoefficient
double eResponseCoefficient
Definition: HCALResponse.h:80
HCALResponse::eResponseScale
double eResponseScale[3]
Definition: HCALResponse.h:77
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
HCALResponse::parNames
std::vector< std::string > parNames
Definition: HCALResponse.h:103
HCALResponse::eResponseExponent
double eResponseExponent
Definition: HCALResponse.h:79
vec5
std::vector< vec4 > vec5
Definition: HCALResponse.h:19
vec3
std::vector< vec2 > vec3
Definition: HCALResponse.h:17
HCALResponse::PoissonParameters
vec3 PoissonParameters
Definition: HCALResponse.h:114
vec1
std::vector< double > vec1
Definition: HCALResponse.h:15
HCALResponse::muStep
double muStep
Definition: HCALResponse.h:90
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
vec4
std::vector< vec3 > vec4
Definition: HCALResponse.h:18
BPhysicsValidation_cfi.pname
pname
Definition: BPhysicsValidation_cfi.py:6
gen::k
int k[5][pyjets_maxn]
Definition: Cascade2Hadronizer.cc:79
HCALResponse::eGridHD
vec1 eGridHD[4]
Definition: HCALResponse.h:95
HCALResponse::respFactorEM
double respFactorEM
Definition: HCALResponse.h:92
tmp
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
HCALResponse::maxHDetas
int maxHDetas[3]
Definition: HCALResponse.h:88
dedxEstimators_cff.fraction
fraction
Definition: dedxEstimators_cff.py:23
vec2
std::vector< vec1 > vec2
Definition: HCALResponse.h:16
HCALResponse::eResponsePlateau
double eResponsePlateau[3]
Definition: HCALResponse.h:78
HCALResponse::~HCALResponse ( )
inline

Definition at line 33 of file HCALResponse.h.

References MillePedeFileConverter_cfg::e, and PVValHelper::eta.

33 { }

Member Function Documentation

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

Definition at line 574 of file HCALResponse.cc.

References MillePedeFileConverter_cfg::out.

575  {
576  double out = cball.shoot(mu,sigma,aL,nL,aR,nR, random);
577  if (mu >= 0.) {
578  while (out < 0.) out = cball.shoot(mu,sigma,aL,nL,aR,nR, random);
579  }
580  //else give up on re-trying, otherwise too much time can be lost before emeas comes out positive
581 
582  return out;
583 }
DoubleCrystalBallGenerator::shoot
double shoot(double mu, double sigma, double aL, double nL, double aR, double nR, RandomEngineAndDistribution const *random)
Definition: DoubleCrystalBallGenerator.cc:11
random
TRandom random
Definition: MVATrainer.cc:138
RPCpg::mu
const int mu
Definition: Constants.h:22
MillePedeFileConverter_cfg.out
out
Definition: MillePedeFileConverter_cfg.py:31
HCALResponse::cball
DoubleCrystalBallGenerator cball
Definition: HCALResponse.h:117
void HCALResponse::correctHF ( double  e,
int  type 
)

Definition at line 594 of file HCALResponse.cc.

References funct::abs(), mps_fire::i, and maxEta.

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

594  {
595 
596  int jmin = 0;
597  for (int i = 0; i < maxEne; i++) {
598  if(ee >= energyHF[i]) jmin = i;
599  }
600 
601  double x1, x2;
602  double y1em, y2em;
603  double y1had, y2had;
604  for(int i=0; i<maxEta; ++i) {
605  if(ee < energyHF[0]) {
606  if(abs(type)==11 || abs(type)==22) {
607  corrHFem[i] = corrHFgEm[0][i];
608  corrHFhad[i] = corrHFgHad[0][i];
609  } else {
610  corrHFem[i] = corrHFhEm[0][i];
611  corrHFhad[i] = corrHFhHad[0][i];
612  }
613  } else if(jmin >= maxEne-1) {
614  if(abs(type)==11 || abs(type)==22) {
615  corrHFem[i] = corrHFgEm[maxEta][i];
616  corrHFhad[i] = corrHFgHad[maxEta][i];
617  } else {
618  corrHFem[i] = corrHFhEm[maxEta][i];
619  corrHFhad[i] = corrHFhHad[maxEta][i];
620  }
621  } else {
622  x1 = energyHF[jmin];
623  x2 = energyHF[jmin+1];
624  if(abs(type)==11 || abs(type)==22) {
625  y1em = corrHFgEm[jmin][i];
626  y2em = corrHFgEm[jmin+1][i];
627  y1had = corrHFgHad[jmin][i];
628  y2had = corrHFgHad[jmin+1][i];
629  } else {
630  y1em = corrHFhEm[jmin][i];
631  y2em = corrHFhEm[jmin+1][i];
632  y1had = corrHFhHad[jmin][i];
633  y2had = corrHFhHad[jmin+1][i];
634  }
635  corrHFem[i] = y1em + (ee-x1)*((y2em-y1em)/(x2-x1));
636  corrHFhad[i] = y1had + (ee-x1)*((y2had-y1had)/(x2-x1));
637  }
638  }
639 
640 }
type
type
Definition: HCALResponse.h:21
mps_fire.i
i
Definition: mps_fire.py:156
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
double HCALResponse::gaussShootNoNegative ( double  e,
double  sigma,
RandomEngineAndDistribution const *  random 
)
private

Definition at line 563 of file HCALResponse.cc.

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

563  {
564  double out = random->gaussShoot(e,sigma);
565  if (e >= 0.) {
566  while (out < 0.) out = random->gaussShoot(e,sigma);
567  }
568  //else give up on re-trying, otherwise too much time can be lost before emeas comes out positive
569 
570  return out;
571 }
random
TRandom random
Definition: MVATrainer.cc:138
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
MillePedeFileConverter_cfg.out
out
Definition: MillePedeFileConverter_cfg.py:31
vec1& HCALResponse::getCorrHFem ( )
inline

Definition at line 48 of file HCALResponse.h.

Referenced by CalorimetryManager::updateHCAL().

48 {return corrHFem;}
vec1& HCALResponse::getCorrHFhad ( )
inline

Definition at line 49 of file HCALResponse.h.

References RPCpg::mu.

Referenced by CalorimetryManager::updateHCAL().

49 {return corrHFhad;}
int HCALResponse::getDet ( int  ieta)
private

Definition at line 552 of file HCALResponse.cc.

References edmIntegrityCheck::d.

552  {
553  int d;
554  for(d = 0; d < 2; d++){
555  if(ieta < HDeta[d+1]){
556  break;
557  }
558  }
559  return d;
560 }
HCALResponse::HDeta
int HDeta[4]
Definition: HCALResponse.h:88
double HCALResponse::getHCALEnergyResponse ( double  e,
int  hit,
RandomEngineAndDistribution const *  random 
)

Definition at line 527 of file HCALResponse.cc.

References EnergyCorrector::c, JetChargeProducer_cfi::exp, hcforward, cmsBatch::log, gen::n, AlCaHLTBitMon_ParallelJobs::p, fftjetproducer_cfi::resolution, alignCSCRings::s, mathSSE::sqrt(), and VFCAL.

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

527  {
528  //response
529  double s = eResponseScale[hit];
530  double n = eResponseExponent;
531  double p = eResponsePlateau[hit];
532  double c = eResponseCoefficient;
533 
534  double response = e * p / (1+c*exp(n * log(s/e)));
535 
536  if(response<0.) response = 0.;
537 
538  //resolution
539  double resolution;
540  if(hit==hcforward)
541  resolution = e *sqrt( RespPar[VFCAL][1][0]*RespPar[VFCAL][1][0] / e + RespPar[VFCAL][1][1]*RespPar[VFCAL][1][1] );
542  else
543  resolution = e * sqrt( RespPar[HCAL][hit][0]*RespPar[HCAL][hit][0]/(e) + RespPar[HCAL][hit][1]*RespPar[HCAL][hit][1] );
544 
545  //random smearing
546  double rndm = gaussShootNoNegative(response, resolution, random);
547 
548  return rndm;
549 }
HCALResponse::RespPar
double RespPar[3][2][3]
Definition: HCALResponse.h:74
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
HCALResponse::eResponseCoefficient
double eResponseCoefficient
Definition: HCALResponse.h:80
HCALResponse::eResponseScale
double eResponseScale[3]
Definition: HCALResponse.h:77
random
TRandom random
Definition: MVATrainer.cc:138
HCALResponse::eResponseExponent
double eResponseExponent
Definition: HCALResponse.h:79
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
HCAL
Definition: HCAL.py:1
fftjetproducer_cfi.resolution
resolution
Definition: fftjetproducer_cfi.py:145
JetChargeProducer_cfi.exp
exp
Definition: JetChargeProducer_cfi.py:6
HCALResponse::eResponsePlateau
double eResponsePlateau[3]
Definition: HCALResponse.h:78
HCALResponse::gaussShootNoNegative
double gaussShootNoNegative(double e, double sigma, RandomEngineAndDistribution const *)
Definition: HCALResponse.cc:563
double HCALResponse::getMIPfraction ( double  energy,
double  eta 
)

Definition at line 237 of file HCALResponse.cc.

References funct::abs(), mps_fire::i, and RecoTauDiscriminantConfiguration::mean.

Referenced by CalorimetryManager::HDShowerSimulation().

237  {
238  int ieta = abs((int)(eta / etaStep)) ;
239  int ie = -1;
240  //check eta and det
241  int det = getDet(ieta);
242  int deta = ieta - HDeta[det];
243  if(deta >= maxHDetas[det]) deta = maxHDetas[det] - 1;
244  else if(deta < 0 ) deta = 0;
245  //find energy range
246  for (int i = 0; i < maxHDe[det]; i++) {
247  if(energy < eGridHD[det][i]) {
248  if(i == 0) return mipfraction [det][0][deta]; // less than minimal - the first value is used instead of extrapolating
249  else ie = i-1;
250  break;
251  }
252  }
253  if(ie == -1) return mipfraction [det][maxHDe[det]-1][deta]; // more than maximal - the last value is used instead of extrapolating
254  double y1, y2;
255  double x1 = eGridHD[det][ie];
256  double x2 = eGridHD[det][ie+1];
257  y1=mipfraction[det][ie][deta];
258  y2=mipfraction[det][ie+1][deta];
259  double mean = 0;
260  mean=(y1*(x2-energy) + y2*(energy-x1))/(x2-x1);
261  return mean;
262 }
HCALResponse::etaStep
double etaStep
Definition: HCALResponse.h:86
mps_fire.i
i
Definition: mps_fire.py:156
HCALResponse::maxHDe
int maxHDe[4]
Definition: HCALResponse.h:84
HCALResponse::HDeta
int HDeta[4]
Definition: HCALResponse.h:88
HCALResponse::getDet
int getDet(int ieta)
Definition: HCALResponse.cc:552
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
HCALResponse::eGridHD
vec1 eGridHD[4]
Definition: HCALResponse.h:95
HCALResponse::maxHDetas
int maxHDetas[3]
Definition: HCALResponse.h:88
double HCALResponse::interEM ( double  e,
int  ie,
int  ieta,
RandomEngineAndDistribution const *  random 
)
private

Definition at line 489 of file HCALResponse.cc.

References debug, MillePedeFileConverter_cfg::e, and RecoTauDiscriminantConfiguration::mean.

489  {
490  double y1 = meanEM[ie][ieta];
491  double y2 = meanEM[ie+1][ieta];
492  double x1 = eGridEM[ie];
493  double x2 = eGridEM[ie+1];
494 
495  if(debug) {
496  // cout.width(6);
497  LogInfo("FastCalorimetry") << std::endl
498  << " HCALResponse::interEM mean " << std::endl
499  << " x, x1-x2, y1-y2 = "
500  << e << ", " << x1 <<"-" << x2 << " " << y1 <<"-" << y2 << std::endl;
501  }
502 
503  //linear interpolation
504  double mean = (y1*(x2-e) + y2*(e-x1))/(x2-x1);
505 
506  y1 = sigmaEM[ie][ieta];
507  y2 = sigmaEM[ie+1][ieta];
508 
509  if(debug) {
510  // cout.width(6);
511  LogInfo("FastCalorimetry") << std::endl
512  << " HCALResponse::interEM sigma" << std::endl
513  << " x, x1-x2, y1-y2 = "
514  << e << ", " << x1 <<"-" << x2 << " " << y1 <<"-" << y2 << std::endl;
515  }
516 
517  //linear interpolation
518  double sigma = (y1*(x2-e) + y2*(e-x1))/(x2-x1);
519 
520  //random smearing
521  double rndm = gaussShootNoNegative(mean, sigma, random);
522 
523  return rndm;
524 }
random
TRandom random
Definition: MVATrainer.cc:138
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
HCALResponse::gaussShootNoNegative
double gaussShootNoNegative(double e, double sigma, RandomEngineAndDistribution const *)
Definition: HCALResponse.cc:563
double HCALResponse::interHD ( int  mip,
double  e,
int  ie,
int  ieta,
int  det,
RandomEngineAndDistribution const *  random 
)
private

Definition at line 418 of file HCALResponse.cc.

References MillePedeFileConverter_cfg::e, RandomEngineAndDistribution::flatShoot(), RecoTauDiscriminantConfiguration::mean, AlCaHLTBitMon_ParallelJobs::p, metProducer_cfi::parameters, RandomEngineAndDistribution::poissonShoot(), and tmp.

418  {
419  double x1, x2;
420  double y1, y2;
421  if(det==2) mip=2; //ignore mip status for HF
422  double mean = 0;
423  vec1 pars(nPar,0);
424 
425  // for ieta < 5 there is overlap between HE and HF, and measurement comes from HE
426  if (det==2 && ieta>5 && e<20){
427 
428  for(int p = 0; p < 4; p++){
429  y1=PoissonParameters[p][ie][ieta];
430  y2=PoissonParameters[p][ie+1][ieta];
431  if(e>5){
432  x1 = eGridHD[det+1][ie];
433  x2 = eGridHD[det+1][ie+1];
434  pars[p] = (y1*(x2-e) + y2*(e-x1))/(x2-x1);
435  }
436  else pars[p] = y1;
437  }
438  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);
439  }
440 
441  else{
442 
443  x1 = eGridHD[det][ie];
444  x2 = eGridHD[det][ie+1];
445 
446  //calculate all parameters
447  for(int p = 0; p < nPar; p++){
448  y1 = parameters[p][mip][det][ie][ieta];
449  y2 = parameters[p][mip][det][ie+1][ieta];
450 
451  //par-specific checks
452  double custom = 0;
453  bool use_custom = false;
454 
455  //do not let mu or sigma get extrapolated below zero for low energies
456  //especially important for HF since extrapolation is used for E < 15 GeV
457  if((p==0 || p==1) && e < x1){
458  double tmp = (y1*x2-y2*x1)/(x2-x1); //extrapolate down to e=0
459  if(tmp<0) { //require mu,sigma > 0 for E > 0
460  custom = y1*e/x1;
461  use_custom = true;
462  }
463  }
464  //tail parameters have lower bounds - never extrapolate down
465  else if((p==2 || p==3 || p==4 || p==5)){
466  if(e < x1 && y1 < y2){
467  custom = y1;
468  use_custom = true;
469  }
470  else if(e > x2 && y2 < y1){
471  custom = y2;
472  use_custom = true;
473  }
474  }
475 
476  //linear interpolation
477  if(use_custom) pars[p] = custom;
478  else pars[p] = (y1*(x2-e) + y2*(e-x1))/(x2-x1);
479  }
480 
481  //random smearing
482  if(nPar==6) mean = cballShootNoNegative(pars[0],pars[1],pars[2],pars[3],pars[4],pars[5], random);
483  else if(nPar==2) mean = gaussShootNoNegative(pars[0],pars[1], random); //gaussian fallback
484  }
485  return mean;
486 }
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
random
TRandom random
Definition: MVATrainer.cc:138
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
HCALResponse::PoissonShootNoNegative
double PoissonShootNoNegative(double e, double sigma, RandomEngineAndDistribution const *)
Definition: HCALResponse.cc:584
HCALResponse::PoissonParameters
vec3 PoissonParameters
Definition: HCALResponse.h:114
vec1
std::vector< double > vec1
Definition: HCALResponse.h:15
HCALResponse::eGridHD
vec1 eGridHD[4]
Definition: HCALResponse.h:95
tmp
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
HCALResponse::gaussShootNoNegative
double gaussShootNoNegative(double e, double sigma, RandomEngineAndDistribution const *)
Definition: HCALResponse.cc:563
HCALResponse::cballShootNoNegative
double cballShootNoNegative(double mu, double sigma, double aL, double nL, double aR, double nR, RandomEngineAndDistribution const *)
Definition: HCALResponse.cc:574
double HCALResponse::interMU ( double  e,
int  ie,
int  ieta,
RandomEngineAndDistribution const *  random 
)
private

Definition at line 372 of file HCALResponse.cc.

References debug, MillePedeFileConverter_cfg::e, RandomEngineAndDistribution::flatShoot(), mps_fire::i, RecoTauDiscriminantConfiguration::mean, and x().

372  {
373  double x = random->flatShoot();
374 
375  int bin1 = maxMUbin;
376  for(int i = 0; i < maxMUbin; i++) {
377  if(x > responseMU[ie][ieta][i]) {
378  bin1 = i-1;
379  break;
380  }
381  }
382  int bin2 = maxMUbin;
383  for(int i = 0; i < maxMUbin; i++) {
384  if(x > responseMU[ie+1][ieta][i]) {
385  bin2 = i-1;
386  break;
387  }
388  }
389 
390  double x1 = eGridMU[ie];
391  double x2 = eGridMU[ie+1];
392  double y1 = (bin1 + random->flatShoot()) * muStep;
393  double y2 = (bin2 + random->flatShoot()) * muStep;
394 
395  if(debug) {
396  // cout.width(6);
397  LogInfo("FastCalorimetry") << std::endl
398  << " HCALResponse::interMU " << std::endl
399  << " x, x1-x2, y1-y2 = "
400  << e << ", " << x1 <<"-" << x2 << " " << y1 <<"-" << y2 << std::endl;
401 
402  }
403 
404  //linear interpolation
405  double mean = (y1*(x2-e) + y2*(e-x1))/(x2-x1);
406 
407  if(debug) {
408  //cout.width(6);
409  LogInfo("FastCalorimetry") << std::endl
410  << " HCALResponse::interMU " << std::endl
411  << " e, ie, ieta = " << e << " " << ie << " " << ieta << std::endl
412  << " response = " << mean << std::endl;
413  }
414 
415  return mean;
416 }
mps_fire.i
i
Definition: mps_fire.py:156
random
TRandom random
Definition: MVATrainer.cc:138
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
HCALResponse::muStep
double muStep
Definition: HCALResponse.h:90
double HCALResponse::PoissonShootNoNegative ( double  e,
double  sigma,
RandomEngineAndDistribution const *  random 
)
private

Definition at line 584 of file HCALResponse.cc.

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

584  {
585  double out = -1;
586  while (out < 0.) {
587  out = random->poissonShoot(e);
588  out = out + sigma;
589  }
590  return out;
591 
592 }
random
TRandom random
Definition: MVATrainer.cc:138
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
MillePedeFileConverter_cfg.out
out
Definition: MillePedeFileConverter_cfg.py:31
double HCALResponse::responseHCAL ( int  _mip,
double  energy,
double  eta,
int  partype,
RandomEngineAndDistribution const *  random 
)

Definition at line 264 of file HCALResponse.cc.

References funct::abs(), debug, mps_fire::i, and RecoTauDiscriminantConfiguration::mean.

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

264  {
265 
266 
267  int ieta = abs((int)(eta / etaStep)) ;
268  int ie = -1;
269 
270  int mip;
271  if(usemip) mip = _mip;
272  else mip = 2; //ignore mip, use only overall (mip + nomip) parameters
273 
274  double mean = 0;
275 
276  // e/gamma in HF
277  if(partype == 0) {
278  //check eta
279  ieta -= HDeta[2]; // HF starts at ieta=30 till ieta=51
280  // but resp.vector from index=0 through 20
281  if(ieta >= maxEMeta ) ieta = maxEMeta-1;
282  else if(ieta < 0) ieta = 0;
283 
284  //find energy range
285  for (int i = 0; i < maxEMe; i++) {
286  if(energy < eGridEM[i]) {
287  if(i == 0) ie = 0; // less than minimal - back extrapolation with the 1st interval
288  else ie = i-1;
289  break;
290  }
291  }
292  if(ie == -1) ie = maxEMe - 2; // more than maximum - extrapolation with last interval
293 
294  //do smearing
295  mean = interEM(energy, ie, ieta, random);
296  }
297 
298  // hadrons
299  else if(partype == 1) {
300  //check eta and det
301  int det = getDet(ieta);
302  int deta = ieta - HDeta[det];
303  if(deta >= maxHDetas[det]) deta = maxHDetas[det] - 1;
304  else if(deta < 0 ) deta = 0;
305 
306  //find energy range
307  for (int i = 0; i < maxHDe[det]; i++) {
308  if(energy < eGridHD[det][i]) {
309  if(i == 0) ie = 0; // less than minimal - back extrapolation with the 1st interval
310  else ie = i-1;
311  break;
312  }
313  }
314  if(ie == -1) ie = maxHDe[det] - 2; // more than maximum - extrapolation with last interval
315 
316  //different energy smearing for low energy hadrons in HF
317  if(det==2 && energy <20 && deta>5){
318  for (int i = 0; i < maxHDe[3]; i++) {
319  if(energy < eGridHD[3][i]) {
320  if(i == 0) ie = 0; // less than minimal - back extrapolation with the 1st interval
321  else ie = i-1;
322  break;
323  }
324  }
325  }
326  //do smearing
327  mean = interHD(mip, energy, ie, deta, det, random);
328  }
329 
330 
331  // muons
332  else if(partype == 2) {
333  //check eta
334  ieta = maxMUeta;
335  for(int i = 0; i < maxMUeta; i++) {
336  if(fabs(eta) < etaGridMU[i]) {
337  ieta = i;
338  break;
339  }
340  }
341  if(ieta < 0) ieta = 0;
342 
343  if(ieta < maxMUeta) { // HB-HE
344  //find energy range
345  for (int i = 0; i < maxMUe; i++) {
346  if(energy < eGridMU[i]) {
347  if(i == 0) ie = 0; // less than minimal - back extrapolation with the 1st interval
348  else ie = i-1;
349  break;
350  }
351  }
352  if(ie == -1) ie = maxMUe - 2; // more than maximum - extrapolation with last interval
353 
354  //do smearing
355  mean = interMU(energy, ie, ieta, random);
356  if(mean > energy) mean = energy;
357  }
358  }
359 
360  // debugging
361  if(debug) {
362  // cout.width(6);
363  LogInfo("FastCalorimetry") << std::endl
364  << " HCALResponse::responseHCAL, partype = " << partype
365  << " E, eta = " << energy << " " << eta
366  << " mean = " << mean << std::endl;
367  }
368 
369  return mean;
370 }
HCALResponse::etaStep
double etaStep
Definition: HCALResponse.h:86
mps_fire.i
i
Definition: mps_fire.py:156
HCALResponse::maxHDe
int maxHDe[4]
Definition: HCALResponse.h:84
HCALResponse::HDeta
int HDeta[4]
Definition: HCALResponse.h:88
HCALResponse::getDet
int getDet(int ieta)
Definition: HCALResponse.cc:552
HCALResponse::interEM
double interEM(double e, int ie, int ieta, RandomEngineAndDistribution const *)
Definition: HCALResponse.cc:489
random
TRandom random
Definition: MVATrainer.cc:138
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
HCALResponse::interHD
double interHD(int mip, double e, int ie, int ieta, int det, RandomEngineAndDistribution const *)
Definition: HCALResponse.cc:418
HCALResponse::interMU
double interMU(double e, int ie, int ieta, RandomEngineAndDistribution const *)
Definition: HCALResponse.cc:372
HCALResponse::eGridHD
vec1 eGridHD[4]
Definition: HCALResponse.h:95
HCALResponse::maxHDetas
int maxHDetas[3]
Definition: HCALResponse.h:88

Member Data Documentation

int HCALResponse::barrelMUeta
private

Definition at line 88 of file HCALResponse.h.

DoubleCrystalBallGenerator HCALResponse::cball
private

Definition at line 117 of file HCALResponse.h.

vec1 HCALResponse::corrHFem
private

Definition at line 124 of file HCALResponse.h.

vec2 HCALResponse::corrHFgEm
private

Definition at line 122 of file HCALResponse.h.

vec2 HCALResponse::corrHFgHad
private

Definition at line 122 of file HCALResponse.h.

vec1 HCALResponse::corrHFhad
private

Definition at line 124 of file HCALResponse.h.

vec2 HCALResponse::corrHFhEm
private

Definition at line 123 of file HCALResponse.h.

vec2 HCALResponse::corrHFhHad
private

Definition at line 123 of file HCALResponse.h.

bool HCALResponse::debug
private

Definition at line 68 of file HCALResponse.h.

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

vec1 HCALResponse::eGridEM
private

Definition at line 96 of file HCALResponse.h.

vec1 HCALResponse::eGridHD[4]
private

Definition at line 95 of file HCALResponse.h.

vec1 HCALResponse::eGridMU
private

Definition at line 97 of file HCALResponse.h.

int HCALResponse::endcapMUeta
private

Definition at line 88 of file HCALResponse.h.

vec1 HCALResponse::energyHF
private

Definition at line 121 of file HCALResponse.h.

double HCALResponse::eResponseCoefficient
private

Definition at line 80 of file HCALResponse.h.

double HCALResponse::eResponseExponent
private

Definition at line 79 of file HCALResponse.h.

double HCALResponse::eResponsePlateau[3]
private

Definition at line 78 of file HCALResponse.h.

double HCALResponse::eResponseScale[3]
private

Definition at line 77 of file HCALResponse.h.

vec1 HCALResponse::etaGridMU
private

Definition at line 98 of file HCALResponse.h.

double HCALResponse::etaStep
private

Definition at line 86 of file HCALResponse.h.

int HCALResponse::HDeta[4]
private

Definition at line 88 of file HCALResponse.h.

int HCALResponse::maxEMe
private

Definition at line 83 of file HCALResponse.h.

int HCALResponse::maxEMeta
private

Definition at line 83 of file HCALResponse.h.

int HCALResponse::maxEne
private

Definition at line 120 of file HCALResponse.h.

int HCALResponse::maxEta
private

Definition at line 120 of file HCALResponse.h.

int HCALResponse::maxHDe[4]
private

Definition at line 84 of file HCALResponse.h.

int HCALResponse::maxHDetas[3]
private

Definition at line 88 of file HCALResponse.h.

int HCALResponse::maxMUbin
private

Definition at line 83 of file HCALResponse.h.

int HCALResponse::maxMUe
private

Definition at line 83 of file HCALResponse.h.

int HCALResponse::maxMUeta
private

Definition at line 83 of file HCALResponse.h.

vec2 HCALResponse::meanEM
private

Definition at line 108 of file HCALResponse.h.

vec3 HCALResponse::mipfraction
private

Definition at line 113 of file HCALResponse.h.

double HCALResponse::muStep
private

Definition at line 90 of file HCALResponse.h.

int HCALResponse::nPar
private

Definition at line 102 of file HCALResponse.h.

vec5 HCALResponse::parameters
private

Definition at line 104 of file HCALResponse.h.

Referenced by Vispa.Plugins.ConfigEditor.ConfigDataAccessor.ConfigDataAccessor::inputTags(), Vispa.Plugins.ConfigEditor.ConfigDataAccessor.ConfigDataAccessor::properties(), and Vispa.Plugins.ConfigEditor.ConfigDataAccessor.ConfigDataAccessor::recursePSetProperties().

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

Definition at line 103 of file HCALResponse.h.

vec3 HCALResponse::PoissonParameters
private

Definition at line 114 of file HCALResponse.h.

double HCALResponse::respFactorEM
private

Definition at line 92 of file HCALResponse.h.

vec3 HCALResponse::responseMU
private

Definition at line 112 of file HCALResponse.h.

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

Definition at line 74 of file HCALResponse.h.

vec2 HCALResponse::sigmaEM
private

Definition at line 108 of file HCALResponse.h.

bool HCALResponse::usemip
private

Definition at line 68 of file HCALResponse.h.

Generated for CMSSW Reference Manual by doxygen 1.8.11