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

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
 
bool debug
 
vec1 eGridEM
 
vec1 eGridHD [3]
 
vec1 eGridMU
 
int endcapMUeta
 
double eResponseCoefficient
 
double eResponseExponent
 
double eResponsePlateau [3]
 
double eResponseScale [3]
 
vec1 etaGridMU
 
double etaStep
 
int HDeta [4]
 
int maxEMe
 
int maxEMeta
 
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(), debug, edm::ParameterSet::getParameter(), HCAL, i, j, gen::k, m, 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 }
HCALResponse::etaStep
double etaStep
Definition: HCALResponse.h:81
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
HCALResponse::RespPar
double RespPar[3][2][3]
Definition: HCALResponse.h:69
i
int i
Definition: DBlmapReader.cc:9
HCALResponse::maxHDe
int maxHDe[4]
Definition: HCALResponse.h:79
HCALResponse::HDeta
int HDeta[4]
Definition: HCALResponse.h:83
HCALResponse::eResponseCoefficient
double eResponseCoefficient
Definition: HCALResponse.h:75
HCALResponse::eResponseScale
double eResponseScale[3]
Definition: HCALResponse.h:72
HCALResponse::parNames
std::vector< std::string > parNames
Definition: HCALResponse.h:98
HCALResponse::eResponseExponent
double eResponseExponent
Definition: HCALResponse.h:74
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
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:109
HCALResponse::eGridHD
vec1 eGridHD[3]
Definition: HCALResponse.h:90
vec1
std::vector< double > vec1
Definition: HCALResponse.h:15
HCALResponse::muStep
double muStep
Definition: HCALResponse.h:85
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
gen::k
int k[5][pyjets_maxn]
Definition: Cascade2Hadronizer.cc:79
HCALResponse::respFactorEM
double respFactorEM
Definition: HCALResponse.h:87
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:83
vec2
std::vector< vec1 > vec2
Definition: HCALResponse.h:16
HCALResponse::eResponsePlateau
double eResponsePlateau[3]
Definition: HCALResponse.h:73
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 543 of file HCALResponse.cc.

References dbtoconf::out.

544  {
545  double out = cball.shoot(mu,sigma,aL,nL,aR,nR, random);
546  if (mu >= 0.) {
547  while (out < 0.) out = cball.shoot(mu,sigma,aL,nL,aR,nR, random);
548  }
549  //else give up on re-trying, otherwise too much time can be lost before emeas comes out positive
550 
551  return out;
552 }
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:112
double HCALResponse::gaussShootNoNegative ( double  e,
double  sigma,
RandomEngineAndDistribution const *  random 
)
private

Definition at line 532 of file HCALResponse.cc.

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

532  {
533  double out = random->gaussShoot(e,sigma);
534  if (e >= 0.) {
535  while (out < 0.) out = random->gaussShoot(e,sigma);
536  }
537  //else give up on re-trying, otherwise too much time can be lost before emeas comes out positive
538 
539  return out;
540 }
random
TRandom random
Definition: MVATrainer.cc:138
dbtoconf.out
tuple out
Definition: dbtoconf.py:99
int HCALResponse::getDet ( int  ieta)
private

Definition at line 521 of file HCALResponse.cc.

521  {
522  int d;
523  for(d = 0; d < 2; d++){
524  if(ieta < HDeta[d+1]){
525  break;
526  }
527  }
528  return d;
529 }
HCALResponse::HDeta
int HDeta[4]
Definition: HCALResponse.h:83
double HCALResponse::getHCALEnergyResponse ( double  e,
int  hit,
RandomEngineAndDistribution const *  random 
)

Definition at line 496 of file HCALResponse.cc.

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

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

496  {
497  //response
498  double s = eResponseScale[hit];
499  double n = eResponseExponent;
500  double p = eResponsePlateau[hit];
501  double c = eResponseCoefficient;
502 
503  double response = e * p / (1+c*exp(n * log(s/e)));
504 
505  if(response<0.) response = 0.;
506 
507  //resolution
508  double resolution;
509  if(hit==hcforward)
510  resolution = e *sqrt( RespPar[VFCAL][1][0]*RespPar[VFCAL][1][0] / e + RespPar[VFCAL][1][1]*RespPar[VFCAL][1][1] );
511  else
512  resolution = e * sqrt( RespPar[HCAL][hit][0]*RespPar[HCAL][hit][0]/(e) + RespPar[HCAL][hit][1]*RespPar[HCAL][hit][1] );
513 
514  //random smearing
515  double rndm = gaussShootNoNegative(response, resolution, random);
516 
517  return rndm;
518 }
HCALResponse::RespPar
double RespPar[3][2][3]
Definition: HCALResponse.h:69
log
static std::vector< std::string > checklist log
Definition: LHERunInfoProduct.cc:208
create_public_lumi_plots.exp
tuple exp
Definition: create_public_lumi_plots.py:1087
HCALResponse::eResponseCoefficient
double eResponseCoefficient
Definition: HCALResponse.h:75
HCALResponse::eResponseScale
double eResponseScale[3]
Definition: HCALResponse.h:72
random
TRandom random
Definition: MVATrainer.cc:138
HCALResponse::eResponseExponent
double eResponseExponent
Definition: HCALResponse.h:74
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:73
HCALResponse::gaussShootNoNegative
double gaussShootNoNegative(double e, double sigma, RandomEngineAndDistribution const *)
Definition: HCALResponse.cc:532
double HCALResponse::getMIPfraction ( double  energy,
double  eta 
)

Definition at line 211 of file HCALResponse.cc.

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

Referenced by CalorimetryManager::HDShowerSimulation().

211  {
212  int ieta = abs((int)(eta / etaStep)) ;
213  int ie = -1;
214  //check eta and det
215  int det = getDet(ieta);
216  int deta = ieta - HDeta[det];
217  if(deta >= maxHDetas[det]) deta = maxHDetas[det] - 1;
218  else if(deta < 0 ) deta = 0;
219  //find energy range
220  for (int i = 0; i < maxHDe[det]; i++) {
221  if(energy < eGridHD[det][i]) {
222  if(i == 0) return mipfraction [det][0][deta]; // less than minimal - the first value is used instead of extrapolating
223  else ie = i-1;
224  break;
225  }
226  }
227  if(ie == -1) return mipfraction [det][maxHDe[det]-1][deta]; // more than maximal - the last value is used instead of extrapolating
228  double y1, y2;
229  double x1 = eGridHD[det][ie];
230  double x2 = eGridHD[det][ie+1];
231  y1=mipfraction[det][ie][deta];
232  y2=mipfraction[det][ie+1][deta];
233  double mean = 0;
234  mean=(y1*(x2-energy) + y2*(energy-x1))/(x2-x1);
235  return mean;
236  }
HCALResponse::etaStep
double etaStep
Definition: HCALResponse.h:81
i
int i
Definition: DBlmapReader.cc:9
HCALResponse::maxHDe
int maxHDe[4]
Definition: HCALResponse.h:79
HCALResponse::HDeta
int HDeta[4]
Definition: HCALResponse.h:83
HCALResponse::getDet
int getDet(int ieta)
Definition: HCALResponse.cc:521
eta
T eta() const
Definition: Basic3DVectorLD.h:177
HCALResponse::eGridHD
vec1 eGridHD[3]
Definition: HCALResponse.h:90
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
HCALResponse::maxHDetas
int maxHDetas[3]
Definition: HCALResponse.h:83
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 458 of file HCALResponse.cc.

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

458  {
459  double y1 = meanEM[ie][ieta];
460  double y2 = meanEM[ie+1][ieta];
461  double x1 = eGridEM[ie];
462  double x2 = eGridEM[ie+1];
463 
464  if(debug) {
465  // cout.width(6);
466  LogInfo("FastCalorimetry") << std::endl
467  << " HCALResponse::interEM mean " << std::endl
468  << " x, x1-x2, y1-y2 = "
469  << e << ", " << x1 <<"-" << x2 << " " << y1 <<"-" << y2 << std::endl;
470  }
471 
472  //linear interpolation
473  double mean = (y1*(x2-e) + y2*(e-x1))/(x2-x1);
474 
475  y1 = sigmaEM[ie][ieta];
476  y2 = sigmaEM[ie+1][ieta];
477 
478  if(debug) {
479  // cout.width(6);
480  LogInfo("FastCalorimetry") << std::endl
481  << " HCALResponse::interEM sigma" << std::endl
482  << " x, x1-x2, y1-y2 = "
483  << e << ", " << x1 <<"-" << x2 << " " << y1 <<"-" << y2 << std::endl;
484  }
485 
486  //linear interpolation
487  double sigma = (y1*(x2-e) + y2*(e-x1))/(x2-x1);
488 
489  //random smearing
490  double rndm = gaussShootNoNegative(mean, sigma, random);
491 
492  return rndm;
493 }
random
TRandom random
Definition: MVATrainer.cc:138
HCALResponse::gaussShootNoNegative
double gaussShootNoNegative(double e, double sigma, RandomEngineAndDistribution const *)
Definition: HCALResponse.cc:532
double HCALResponse::interHD ( int  mip,
double  e,
int  ie,
int  ieta,
int  det,
RandomEngineAndDistribution const *  random 
)
private

Definition at line 390 of file HCALResponse.cc.

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

390  {
391  double x1, x2;
392  double y1, y2;
393  if(det==2) mip=2; //ignore mip status for HF
394  double mean = 0;
395  vec1 pars(nPar,0);
396 
397  if (det==2 && ieta>5 && e<20){
398 
399  x1 = eGridHD[det+1][ie];
400  x2 = eGridHD[det+1][ie+1];
401  for(int p = 0; p < 4; p++){
402  y1=PoissonParameters[p][ie][ieta];
403  y2=PoissonParameters[p][ie+1][ieta];
404  if(e>5)pars[p] = (y1*(x2-e) + y2*(e-x1))/(x2-x1);
405  else pars[p] = y1;
406  }
407  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);
408  }
409 
410  else{
411 
412  x1 = eGridHD[det][ie];
413  x2 = eGridHD[det][ie+1];
414 
415  //calculate all parameters
416  for(int p = 0; p < nPar; p++){
417  y1 = parameters[p][mip][det][ie][ieta];
418  y2 = parameters[p][mip][det][ie+1][ieta];
419 
420  //par-specific checks
421  double custom = 0;
422  bool use_custom = false;
423 
424  //do not let mu or sigma get extrapolated below zero for low energies
425  //especially important for HF since extrapolation is used for E < 15 GeV
426  if((p==0 || p==1) && e < x1){
427  double tmp = (y1*x2-y2*x1)/(x2-x1); //extrapolate down to e=0
428  if(tmp<0) { //require mu,sigma > 0 for E > 0
429  custom = y1*e/x1;
430  use_custom = true;
431  }
432  }
433  //tail parameters have lower bounds - never extrapolate down
434  else if((p==2 || p==3 || p==4 || p==5)){
435  if(e < x1 && y1 < y2){
436  custom = y1;
437  use_custom = true;
438  }
439  else if(e > x2 && y2 < y1){
440  custom = y2;
441  use_custom = true;
442  }
443  }
444 
445  //linear interpolation
446  if(use_custom) pars[p] = custom;
447  else pars[p] = (y1*(x2-e) + y2*(e-x1))/(x2-x1);
448  }
449 
450  //random smearing
451  if(nPar==6) mean = cballShootNoNegative(pars[0],pars[1],pars[2],pars[3],pars[4],pars[5], random);
452  else if(nPar==2) mean = gaussShootNoNegative(pars[0],pars[1], random); //gaussian fallback
453  }
454  return mean;
455 }
random
TRandom random
Definition: MVATrainer.cc:138
HCALResponse::PoissonShootNoNegative
double PoissonShootNoNegative(double e, double sigma, RandomEngineAndDistribution const *)
Definition: HCALResponse.cc:553
HCALResponse::PoissonParameters
vec3 PoissonParameters
Definition: HCALResponse.h:109
HCALResponse::eGridHD
vec1 eGridHD[3]
Definition: HCALResponse.h:90
vec1
std::vector< double > vec1
Definition: HCALResponse.h:15
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:532
HCALResponse::cballShootNoNegative
double cballShootNoNegative(double mu, double sigma, double aL, double nL, double aR, double nR, RandomEngineAndDistribution const *)
Definition: HCALResponse.cc:543
double HCALResponse::interMU ( double  e,
int  ie,
int  ieta,
RandomEngineAndDistribution const *  random 
)
private

Definition at line 344 of file HCALResponse.cc.

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

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

Definition at line 553 of file HCALResponse.cc.

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

553  {
554  double out = -1;
555  while (out < 0.) {
556  out = random->poissonShoot(e);
557  out = out + sigma;
558  }
559  return out;
560 
561 }
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 238 of file HCALResponse.cc.

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

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

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

Member Data Documentation

int HCALResponse::barrelMUeta
private

Definition at line 83 of file HCALResponse.h.

DoubleCrystalBallGenerator HCALResponse::cball
private

Definition at line 112 of file HCALResponse.h.

bool HCALResponse::debug
private

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

vec1 HCALResponse::eGridHD[3]
private

Definition at line 90 of file HCALResponse.h.

vec1 HCALResponse::eGridMU
private

Definition at line 92 of file HCALResponse.h.

int HCALResponse::endcapMUeta
private

Definition at line 83 of file HCALResponse.h.

double HCALResponse::eResponseCoefficient
private

Definition at line 75 of file HCALResponse.h.

double HCALResponse::eResponseExponent
private

Definition at line 74 of file HCALResponse.h.

double HCALResponse::eResponsePlateau[3]
private

Definition at line 73 of file HCALResponse.h.

double HCALResponse::eResponseScale[3]
private

Definition at line 72 of file HCALResponse.h.

vec1 HCALResponse::etaGridMU
private

Definition at line 93 of file HCALResponse.h.

double HCALResponse::etaStep
private

Definition at line 81 of file HCALResponse.h.

int HCALResponse::HDeta[4]
private

Definition at line 83 of file HCALResponse.h.

int HCALResponse::maxEMe
private

Definition at line 78 of file HCALResponse.h.

int HCALResponse::maxEMeta
private

Definition at line 78 of file HCALResponse.h.

int HCALResponse::maxHDe[4]
private

Definition at line 79 of file HCALResponse.h.

int HCALResponse::maxHDetas[3]
private

Definition at line 83 of file HCALResponse.h.

int HCALResponse::maxMUbin
private

Definition at line 78 of file HCALResponse.h.

int HCALResponse::maxMUe
private

Definition at line 78 of file HCALResponse.h.

int HCALResponse::maxMUeta
private

Definition at line 78 of file HCALResponse.h.

vec2 HCALResponse::meanEM
private

Definition at line 103 of file HCALResponse.h.

vec3 HCALResponse::mipfraction
private

Definition at line 108 of file HCALResponse.h.

double HCALResponse::muStep
private

Definition at line 85 of file HCALResponse.h.

int HCALResponse::nPar
private

Definition at line 97 of file HCALResponse.h.

vec5 HCALResponse::parameters
private

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

vec3 HCALResponse::PoissonParameters
private

Definition at line 109 of file HCALResponse.h.

double HCALResponse::respFactorEM
private

Definition at line 87 of file HCALResponse.h.

vec3 HCALResponse::responseMU
private

Definition at line 107 of file HCALResponse.h.

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

Definition at line 69 of file HCALResponse.h.

vec2 HCALResponse::sigmaEM
private

Definition at line 103 of file HCALResponse.h.

bool HCALResponse::usemip
private

Definition at line 63 of file HCALResponse.h.

Generated for CMSSW Reference Manual by doxygen 1.8.5