CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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 & getCorrHF ()
 
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 corrHF
 
vec1 corrHFg
 
vec1 corrHFh
 
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(), ztail::d, debug, HLT_25ns14e33_v1_cff::fraction, edm::ParameterSet::getParameter(), HCAL, i, j, relval_steps::k, visualization-live-secondInstance_cfg::m, maxEta, AlCaHLTBitMon_ParallelJobs::p, Parameters::parameters, 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  corrHFg = pset.getParameter<vec1>("corrHFg");
215  corrHFh = pset.getParameter<vec1>("corrHFh");
216  corrHF = vec1(maxEta,0);
217 }
HCALResponse::etaStep
double etaStep
Definition: HCALResponse.h:85
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
HCALResponse::RespPar
double RespPar[3][2][3]
Definition: HCALResponse.h:73
i
int i
Definition: DBlmapReader.cc:9
HLT_25ns14e33_v1_cff.fraction
tuple fraction
Definition: HLT_25ns14e33_v1_cff.py:28703
HCALResponse::maxHDe
int maxHDe[4]
Definition: HCALResponse.h:83
HCALResponse::HDeta
int HDeta[4]
Definition: HCALResponse.h:87
HCALResponse::eResponseCoefficient
double eResponseCoefficient
Definition: HCALResponse.h:79
HCALResponse::eResponseScale
double eResponseScale[3]
Definition: HCALResponse.h:76
HCALResponse::parNames
std::vector< std::string > parNames
Definition: HCALResponse.h:102
HCALResponse::eResponseExponent
double eResponseExponent
Definition: HCALResponse.h:78
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
ztail.d
tuple d
Definition: ztail.py:151
visualization-live-secondInstance_cfg.m
tuple m
Definition: visualization-live-secondInstance_cfg.py:45
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:113
vec1
std::vector< double > vec1
Definition: HCALResponse.h:15
HCALResponse::muStep
double muStep
Definition: HCALResponse.h:89
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
j
int j
Definition: DBlmapReader.cc:9
vec4
std::vector< vec3 > vec4
Definition: HCALResponse.h:18
HCALResponse::eGridHD
vec1 eGridHD[4]
Definition: HCALResponse.h:94
HCALResponse::respFactorEM
double respFactorEM
Definition: HCALResponse.h:91
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
tmp
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
HCALResponse::maxHDetas
int maxHDetas[3]
Definition: HCALResponse.h:87
relval_steps.k
list k
Definition: relval_steps.py:1523
vec2
std::vector< vec1 > vec2
Definition: HCALResponse.h:16
HCALResponse::eResponsePlateau
double eResponsePlateau[3]
Definition: HCALResponse.h:77
HCALResponse::~HCALResponse ( )
inline

Definition at line 33 of file HCALResponse.h.

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 556 of file HCALResponse.cc.

References dbtoconf::out.

557  {
558  double out = cball.shoot(mu,sigma,aL,nL,aR,nR, random);
559  if (mu >= 0.) {
560  while (out < 0.) out = cball.shoot(mu,sigma,aL,nL,aR,nR, random);
561  }
562  //else give up on re-trying, otherwise too much time can be lost before emeas comes out positive
563 
564  return out;
565 }
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
dbtoconf.out
tuple out
Definition: dbtoconf.py:99
HCALResponse::cball
DoubleCrystalBallGenerator cball
Definition: HCALResponse.h:116
void HCALResponse::correctHF ( double  e,
int  type 
)

Definition at line 576 of file HCALResponse.cc.

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

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

576  {
577 
578  int jmin = 0;
579  for (int i = 0; i < maxEne; i++) {
580  if(ee >= energyHF[i]) jmin = i;
581  }
582 
583  double x1, x2, y1, y2;
584  for(int i=0; i<maxEta; ++i) {
585  if(ee < energyHF[0]) {
586  if(abs(type)==11 || abs(type)==22) corrHF[i] = corrHFg[i];
587  else corrHF[i] = corrHFh[i];
588  } else if(jmin >= maxEne-1) {
589  if(abs(type)==11 || abs(type)==22) corrHF[i] = corrHFg[maxEta*jmin+i];
590  else corrHF[i] = corrHFh[maxEta*jmin+i];
591  } else {
592  x1 = energyHF[jmin];
593  x2 = energyHF[jmin+1];
594  if(abs(type)==11 || abs(type)==22) {
595  y1 = corrHFg[maxEta*jmin+i];
596  y2 = corrHFg[maxEta*(jmin+1)+i];
597  } else {
598  y1 = corrHFh[maxEta*jmin+i];
599  y2 = corrHFh[maxEta*(jmin+1)+i];
600  }
601  corrHF[i] = y1 + (ee-x1)*((y2-y1)/(x2-x1));
602  }
603  }
604 
605 }
type
type
Definition: HCALResponse.h:21
i
int i
Definition: DBlmapReader.cc:9
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 545 of file HCALResponse.cc.

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

545  {
546  double out = random->gaussShoot(e,sigma);
547  if (e >= 0.) {
548  while (out < 0.) out = random->gaussShoot(e,sigma);
549  }
550  //else give up on re-trying, otherwise too much time can be lost before emeas comes out positive
551 
552  return out;
553 }
random
TRandom random
Definition: MVATrainer.cc:138
dbtoconf.out
tuple out
Definition: dbtoconf.py:99
vec1& HCALResponse::getCorrHF ( )
inline

Definition at line 48 of file HCALResponse.h.

References corrHF.

Referenced by CalorimetryManager::updateHCAL().

48 {return corrHF;}
int HCALResponse::getDet ( int  ieta)
private

Definition at line 534 of file HCALResponse.cc.

References ztail::d.

534  {
535  int d;
536  for(d = 0; d < 2; d++){
537  if(ieta < HDeta[d+1]){
538  break;
539  }
540  }
541  return d;
542 }
HCALResponse::HDeta
int HDeta[4]
Definition: HCALResponse.h:87
ztail.d
tuple d
Definition: ztail.py:151
double HCALResponse::getHCALEnergyResponse ( double  e,
int  hit,
RandomEngineAndDistribution const *  random 
)

Definition at line 509 of file HCALResponse.cc.

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

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

509  {
510  //response
511  double s = eResponseScale[hit];
512  double n = eResponseExponent;
513  double p = eResponsePlateau[hit];
514  double c = eResponseCoefficient;
515 
516  double response = e * p / (1+c*exp(n * log(s/e)));
517 
518  if(response<0.) response = 0.;
519 
520  //resolution
521  double resolution;
522  if(hit==hcforward)
523  resolution = e *sqrt( RespPar[VFCAL][1][0]*RespPar[VFCAL][1][0] / e + RespPar[VFCAL][1][1]*RespPar[VFCAL][1][1] );
524  else
525  resolution = e * sqrt( RespPar[HCAL][hit][0]*RespPar[HCAL][hit][0]/(e) + RespPar[HCAL][hit][1]*RespPar[HCAL][hit][1] );
526 
527  //random smearing
528  double rndm = gaussShootNoNegative(response, resolution, random);
529 
530  return rndm;
531 }
HCALResponse::RespPar
double RespPar[3][2][3]
Definition: HCALResponse.h:73
create_public_lumi_plots.exp
tuple exp
Definition: create_public_lumi_plots.py:1087
HCALResponse::eResponseCoefficient
double eResponseCoefficient
Definition: HCALResponse.h:79
HCALResponse::eResponseScale
double eResponseScale[3]
Definition: HCALResponse.h:76
random
TRandom random
Definition: MVATrainer.cc:138
HCALResponse::eResponseExponent
double eResponseExponent
Definition: HCALResponse.h:78
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:48
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
dtDQMClient_cfg.resolution
tuple resolution
Definition: dtDQMClient_cfg.py:21
HCALResponse::eResponsePlateau
double eResponsePlateau[3]
Definition: HCALResponse.h:77
HCALResponse::gaussShootNoNegative
double gaussShootNoNegative(double e, double sigma, RandomEngineAndDistribution const *)
Definition: HCALResponse.cc:545
cmsBatch.log
tuple log
Definition: cmsBatch.py:347
double HCALResponse::getMIPfraction ( double  energy,
double  eta 
)

Definition at line 219 of file HCALResponse.cc.

References funct::abs(), relval_parameters_module::energy, i, and timingPdfMaker::mean.

Referenced by CalorimetryManager::HDShowerSimulation().

219  {
220  int ieta = abs((int)(eta / etaStep)) ;
221  int ie = -1;
222  //check eta and det
223  int det = getDet(ieta);
224  int deta = ieta - HDeta[det];
225  if(deta >= maxHDetas[det]) deta = maxHDetas[det] - 1;
226  else if(deta < 0 ) deta = 0;
227  //find energy range
228  for (int i = 0; i < maxHDe[det]; i++) {
229  if(energy < eGridHD[det][i]) {
230  if(i == 0) return mipfraction [det][0][deta]; // less than minimal - the first value is used instead of extrapolating
231  else ie = i-1;
232  break;
233  }
234  }
235  if(ie == -1) return mipfraction [det][maxHDe[det]-1][deta]; // more than maximal - the last value is used instead of extrapolating
236  double y1, y2;
237  double x1 = eGridHD[det][ie];
238  double x2 = eGridHD[det][ie+1];
239  y1=mipfraction[det][ie][deta];
240  y2=mipfraction[det][ie+1][deta];
241  double mean = 0;
242  mean=(y1*(x2-energy) + y2*(energy-x1))/(x2-x1);
243  return mean;
244  }
HCALResponse::etaStep
double etaStep
Definition: HCALResponse.h:85
i
int i
Definition: DBlmapReader.cc:9
HCALResponse::maxHDe
int maxHDe[4]
Definition: HCALResponse.h:83
HCALResponse::HDeta
int HDeta[4]
Definition: HCALResponse.h:87
HCALResponse::getDet
int getDet(int ieta)
Definition: HCALResponse.cc:534
eta
T eta() const
Definition: Basic3DVectorLD.h:177
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
HCALResponse::eGridHD
vec1 eGridHD[4]
Definition: HCALResponse.h:94
HCALResponse::maxHDetas
int maxHDetas[3]
Definition: HCALResponse.h:87
relval_parameters_module.energy
string energy
Definition: relval_parameters_module.py:29
double HCALResponse::interEM ( double  e,
int  ie,
int  ieta,
RandomEngineAndDistribution const *  random 
)
private

Definition at line 471 of file HCALResponse.cc.

References debug, alignCSCRings::e, and timingPdfMaker::mean.

471  {
472  double y1 = meanEM[ie][ieta];
473  double y2 = meanEM[ie+1][ieta];
474  double x1 = eGridEM[ie];
475  double x2 = eGridEM[ie+1];
476 
477  if(debug) {
478  // cout.width(6);
479  LogInfo("FastCalorimetry") << std::endl
480  << " HCALResponse::interEM mean " << std::endl
481  << " x, x1-x2, y1-y2 = "
482  << e << ", " << x1 <<"-" << x2 << " " << y1 <<"-" << y2 << std::endl;
483  }
484 
485  //linear interpolation
486  double mean = (y1*(x2-e) + y2*(e-x1))/(x2-x1);
487 
488  y1 = sigmaEM[ie][ieta];
489  y2 = sigmaEM[ie+1][ieta];
490 
491  if(debug) {
492  // cout.width(6);
493  LogInfo("FastCalorimetry") << std::endl
494  << " HCALResponse::interEM sigma" << std::endl
495  << " x, x1-x2, y1-y2 = "
496  << e << ", " << x1 <<"-" << x2 << " " << y1 <<"-" << y2 << std::endl;
497  }
498 
499  //linear interpolation
500  double sigma = (y1*(x2-e) + y2*(e-x1))/(x2-x1);
501 
502  //random smearing
503  double rndm = gaussShootNoNegative(mean, sigma, random);
504 
505  return rndm;
506 }
random
TRandom random
Definition: MVATrainer.cc:138
HCALResponse::gaussShootNoNegative
double gaussShootNoNegative(double e, double sigma, RandomEngineAndDistribution const *)
Definition: HCALResponse.cc:545
double HCALResponse::interHD ( int  mip,
double  e,
int  ie,
int  ieta,
int  det,
RandomEngineAndDistribution const *  random 
)
private

Definition at line 400 of file HCALResponse.cc.

References alignCSCRings::e, RandomEngineAndDistribution::flatShoot(), timingPdfMaker::mean, AlCaHLTBitMon_ParallelJobs::p, Parameters::parameters, RandomEngineAndDistribution::poissonShoot(), and tmp.

400  {
401  double x1, x2;
402  double y1, y2;
403  if(det==2) mip=2; //ignore mip status for HF
404  double mean = 0;
405  vec1 pars(nPar,0);
406 
407  // for ieta < 5 there is overlap between HE and HF, and measurement comes from HE
408  if (det==2 && ieta>5 && e<20){
409 
410  for(int p = 0; p < 4; p++){
411  y1=PoissonParameters[p][ie][ieta];
412  y2=PoissonParameters[p][ie+1][ieta];
413  if(e>5){
414  x1 = eGridHD[det+1][ie];
415  x2 = eGridHD[det+1][ie+1];
416  pars[p] = (y1*(x2-e) + y2*(e-x1))/(x2-x1);
417  }
418  else pars[p] = y1;
419  }
420  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);
421  }
422 
423  else{
424 
425  x1 = eGridHD[det][ie];
426  x2 = eGridHD[det][ie+1];
427 
428  //calculate all parameters
429  for(int p = 0; p < nPar; p++){
430  y1 = parameters[p][mip][det][ie][ieta];
431  y2 = parameters[p][mip][det][ie+1][ieta];
432 
433  //par-specific checks
434  double custom = 0;
435  bool use_custom = false;
436 
437  //do not let mu or sigma get extrapolated below zero for low energies
438  //especially important for HF since extrapolation is used for E < 15 GeV
439  if((p==0 || p==1) && e < x1){
440  double tmp = (y1*x2-y2*x1)/(x2-x1); //extrapolate down to e=0
441  if(tmp<0) { //require mu,sigma > 0 for E > 0
442  custom = y1*e/x1;
443  use_custom = true;
444  }
445  }
446  //tail parameters have lower bounds - never extrapolate down
447  else if((p==2 || p==3 || p==4 || p==5)){
448  if(e < x1 && y1 < y2){
449  custom = y1;
450  use_custom = true;
451  }
452  else if(e > x2 && y2 < y1){
453  custom = y2;
454  use_custom = true;
455  }
456  }
457 
458  //linear interpolation
459  if(use_custom) pars[p] = custom;
460  else pars[p] = (y1*(x2-e) + y2*(e-x1))/(x2-x1);
461  }
462 
463  //random smearing
464  if(nPar==6) mean = cballShootNoNegative(pars[0],pars[1],pars[2],pars[3],pars[4],pars[5], random);
465  else if(nPar==2) mean = gaussShootNoNegative(pars[0],pars[1], random); //gaussian fallback
466  }
467  return mean;
468 }
random
TRandom random
Definition: MVATrainer.cc:138
HCALResponse::PoissonShootNoNegative
double PoissonShootNoNegative(double e, double sigma, RandomEngineAndDistribution const *)
Definition: HCALResponse.cc:566
HCALResponse::PoissonParameters
vec3 PoissonParameters
Definition: HCALResponse.h:113
vec1
std::vector< double > vec1
Definition: HCALResponse.h:15
HCALResponse::eGridHD
vec1 eGridHD[4]
Definition: HCALResponse.h:94
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
tmp
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
HCALResponse::gaussShootNoNegative
double gaussShootNoNegative(double e, double sigma, RandomEngineAndDistribution const *)
Definition: HCALResponse.cc:545
HCALResponse::cballShootNoNegative
double cballShootNoNegative(double mu, double sigma, double aL, double nL, double aR, double nR, RandomEngineAndDistribution const *)
Definition: HCALResponse.cc:556
double HCALResponse::interMU ( double  e,
int  ie,
int  ieta,
RandomEngineAndDistribution const *  random 
)
private

Definition at line 354 of file HCALResponse.cc.

References debug, alignCSCRings::e, RandomEngineAndDistribution::flatShoot(), i, timingPdfMaker::mean, and x.

354  {
355  double x = random->flatShoot();
356 
357  int bin1 = maxMUbin;
358  for(int i = 0; i < maxMUbin; i++) {
359  if(x > responseMU[ie][ieta][i]) {
360  bin1 = i-1;
361  break;
362  }
363  }
364  int bin2 = maxMUbin;
365  for(int i = 0; i < maxMUbin; i++) {
366  if(x > responseMU[ie+1][ieta][i]) {
367  bin2 = i-1;
368  break;
369  }
370  }
371 
372  double x1 = eGridMU[ie];
373  double x2 = eGridMU[ie+1];
374  double y1 = (bin1 + random->flatShoot()) * muStep;
375  double y2 = (bin2 + random->flatShoot()) * muStep;
376 
377  if(debug) {
378  // cout.width(6);
379  LogInfo("FastCalorimetry") << std::endl
380  << " HCALResponse::interMU " << std::endl
381  << " x, x1-x2, y1-y2 = "
382  << e << ", " << x1 <<"-" << x2 << " " << y1 <<"-" << y2 << std::endl;
383 
384  }
385 
386  //linear interpolation
387  double mean = (y1*(x2-e) + y2*(e-x1))/(x2-x1);
388 
389  if(debug) {
390  //cout.width(6);
391  LogInfo("FastCalorimetry") << std::endl
392  << " HCALResponse::interMU " << std::endl
393  << " e, ie, ieta = " << e << " " << ie << " " << ieta << std::endl
394  << " response = " << mean << std::endl;
395  }
396 
397  return mean;
398 }
i
int i
Definition: DBlmapReader.cc:9
random
TRandom random
Definition: MVATrainer.cc:138
HCALResponse::muStep
double muStep
Definition: HCALResponse.h:89
x
Definition: DDAxes.h:10
double HCALResponse::PoissonShootNoNegative ( double  e,
double  sigma,
RandomEngineAndDistribution const *  random 
)
private

Definition at line 566 of file HCALResponse.cc.

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

566  {
567  double out = -1;
568  while (out < 0.) {
569  out = random->poissonShoot(e);
570  out = out + sigma;
571  }
572  return out;
573 
574 }
random
TRandom random
Definition: MVATrainer.cc:138
dbtoconf.out
tuple out
Definition: dbtoconf.py:99
double HCALResponse::responseHCAL ( int  _mip,
double  energy,
double  eta,
int  partype,
RandomEngineAndDistribution const *  random 
)

Definition at line 246 of file HCALResponse.cc.

References funct::abs(), debug, relval_parameters_module::energy, i, and timingPdfMaker::mean.

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

246  {
247 
248 
249  int ieta = abs((int)(eta / etaStep)) ;
250  int ie = -1;
251 
252  int mip;
253  if(usemip) mip = _mip;
254  else mip = 2; //ignore mip, use only overall (mip + nomip) parameters
255 
256  double mean = 0;
257 
258  // e/gamma in HF
259  if(partype == 0) {
260  //check eta
261  ieta -= HDeta[2]; // HF starts at ieta=30 till ieta=51
262  // but resp.vector from index=0 through 20
263  if(ieta >= maxEMeta ) ieta = maxEMeta-1;
264  else if(ieta < 0) ieta = 0;
265 
266  //find energy range
267  for (int i = 0; i < maxEMe; i++) {
268  if(energy < eGridEM[i]) {
269  if(i == 0) ie = 0; // less than minimal - back extrapolation with the 1st interval
270  else ie = i-1;
271  break;
272  }
273  }
274  if(ie == -1) ie = maxEMe - 2; // more than maximum - extrapolation with last interval
275 
276  //do smearing
277  mean = interEM(energy, ie, ieta, random);
278  }
279 
280  // hadrons
281  else if(partype == 1) {
282  //check eta and det
283  int det = getDet(ieta);
284  int deta = ieta - HDeta[det];
285  if(deta >= maxHDetas[det]) deta = maxHDetas[det] - 1;
286  else if(deta < 0 ) deta = 0;
287 
288  //find energy range
289  for (int i = 0; i < maxHDe[det]; i++) {
290  if(energy < eGridHD[det][i]) {
291  if(i == 0) ie = 0; // less than minimal - back extrapolation with the 1st interval
292  else ie = i-1;
293  break;
294  }
295  }
296  if(ie == -1) ie = maxHDe[det] - 2; // more than maximum - extrapolation with last interval
297 
298 //different energy smearing for low energy hadrons in HF
299  if(det==2 && energy <20 && deta>5){
300  for (int i = 0; i < maxHDe[3]; i++) {
301  if(energy < eGridHD[3][i]) {
302  if(i == 0) ie = 0; // less than minimal - back extrapolation with the 1st interval
303  else ie = i-1;
304  break;
305  }
306  }
307  }
308  //do smearing
309  mean = interHD(mip, energy, ie, deta, det, random);
310  }
311 
312 
313  // muons
314  else if(partype == 2) {
315  //check eta
316  ieta = maxMUeta;
317  for(int i = 0; i < maxMUeta; i++) {
318  if(fabs(eta) < etaGridMU[i]) {
319  ieta = i;
320  break;
321  }
322  }
323  if(ieta < 0) ieta = 0;
324 
325  if(ieta < maxMUeta) { // HB-HE
326  //find energy range
327  for (int i = 0; i < maxMUe; i++) {
328  if(energy < eGridMU[i]) {
329  if(i == 0) ie = 0; // less than minimal - back extrapolation with the 1st interval
330  else ie = i-1;
331  break;
332  }
333  }
334  if(ie == -1) ie = maxMUe - 2; // more than maximum - extrapolation with last interval
335 
336  //do smearing
337  mean = interMU(energy, ie, ieta, random);
338  if(mean > energy) mean = energy;
339  }
340  }
341 
342  // debugging
343  if(debug) {
344  // cout.width(6);
345  LogInfo("FastCalorimetry") << std::endl
346  << " HCALResponse::responseHCAL, partype = " << partype
347  << " E, eta = " << energy << " " << eta
348  << " mean = " << mean << std::endl;
349  }
350 
351  return mean;
352 }
HCALResponse::etaStep
double etaStep
Definition: HCALResponse.h:85
i
int i
Definition: DBlmapReader.cc:9
HCALResponse::maxHDe
int maxHDe[4]
Definition: HCALResponse.h:83
HCALResponse::HDeta
int HDeta[4]
Definition: HCALResponse.h:87
HCALResponse::getDet
int getDet(int ieta)
Definition: HCALResponse.cc:534
HCALResponse::interEM
double interEM(double e, int ie, int ieta, RandomEngineAndDistribution const *)
Definition: HCALResponse.cc:471
random
TRandom random
Definition: MVATrainer.cc:138
eta
T eta() const
Definition: Basic3DVectorLD.h:177
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:400
HCALResponse::interMU
double interMU(double e, int ie, int ieta, RandomEngineAndDistribution const *)
Definition: HCALResponse.cc:354
HCALResponse::eGridHD
vec1 eGridHD[4]
Definition: HCALResponse.h:94
HCALResponse::maxHDetas
int maxHDetas[3]
Definition: HCALResponse.h:87
relval_parameters_module.energy
string energy
Definition: relval_parameters_module.py:29

Member Data Documentation

int HCALResponse::barrelMUeta
private

Definition at line 87 of file HCALResponse.h.

DoubleCrystalBallGenerator HCALResponse::cball
private

Definition at line 116 of file HCALResponse.h.

vec1 HCALResponse::corrHF
private

Definition at line 122 of file HCALResponse.h.

Referenced by getCorrHF().

vec1 HCALResponse::corrHFg
private

Definition at line 121 of file HCALResponse.h.

vec1 HCALResponse::corrHFh
private

Definition at line 121 of file HCALResponse.h.

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(), and pkg.AbstractPkg::write().

vec1 HCALResponse::eGridEM
private

Definition at line 95 of file HCALResponse.h.

vec1 HCALResponse::eGridHD[4]
private

Definition at line 94 of file HCALResponse.h.

vec1 HCALResponse::eGridMU
private

Definition at line 96 of file HCALResponse.h.

int HCALResponse::endcapMUeta
private

Definition at line 87 of file HCALResponse.h.

vec1 HCALResponse::energyHF
private

Definition at line 120 of file HCALResponse.h.

double HCALResponse::eResponseCoefficient
private

Definition at line 79 of file HCALResponse.h.

double HCALResponse::eResponseExponent
private

Definition at line 78 of file HCALResponse.h.

double HCALResponse::eResponsePlateau[3]
private

Definition at line 77 of file HCALResponse.h.

double HCALResponse::eResponseScale[3]
private

Definition at line 76 of file HCALResponse.h.

vec1 HCALResponse::etaGridMU
private

Definition at line 97 of file HCALResponse.h.

double HCALResponse::etaStep
private

Definition at line 85 of file HCALResponse.h.

int HCALResponse::HDeta[4]
private

Definition at line 87 of file HCALResponse.h.

int HCALResponse::maxEMe
private

Definition at line 82 of file HCALResponse.h.

int HCALResponse::maxEMeta
private

Definition at line 82 of file HCALResponse.h.

int HCALResponse::maxEne
private

Definition at line 119 of file HCALResponse.h.

int HCALResponse::maxEta
private

Definition at line 119 of file HCALResponse.h.

int HCALResponse::maxHDe[4]
private

Definition at line 83 of file HCALResponse.h.

int HCALResponse::maxHDetas[3]
private

Definition at line 87 of file HCALResponse.h.

int HCALResponse::maxMUbin
private

Definition at line 82 of file HCALResponse.h.

int HCALResponse::maxMUe
private

Definition at line 82 of file HCALResponse.h.

int HCALResponse::maxMUeta
private

Definition at line 82 of file HCALResponse.h.

vec2 HCALResponse::meanEM
private

Definition at line 107 of file HCALResponse.h.

vec3 HCALResponse::mipfraction
private

Definition at line 112 of file HCALResponse.h.

double HCALResponse::muStep
private

Definition at line 89 of file HCALResponse.h.

int HCALResponse::nPar
private

Definition at line 101 of file HCALResponse.h.

vec5 HCALResponse::parameters
private

Definition at line 103 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 102 of file HCALResponse.h.

vec3 HCALResponse::PoissonParameters
private

Definition at line 113 of file HCALResponse.h.

double HCALResponse::respFactorEM
private

Definition at line 91 of file HCALResponse.h.

vec3 HCALResponse::responseMU
private

Definition at line 111 of file HCALResponse.h.

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

Definition at line 73 of file HCALResponse.h.

vec2 HCALResponse::sigmaEM
private

Definition at line 107 of file HCALResponse.h.

bool HCALResponse::usemip
private

Definition at line 67 of file HCALResponse.h.

Generated for CMSSW Reference Manual by doxygen 1.8.5