Calorimetry/src/HCALResponse.cc

0001 //updated by Reza Goldouzian
0002 //FastSimulation headers
0003 #include "FastSimulation/Calorimetry/interface/HCALResponse.h"
0004 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
0005 #include "FastSimulation/Utilities/interface/DoubleCrystalBallGenerator.h"
0006
0007 // CMSSW Headers
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010
0011 #include <iostream>
0012 #include <vector>
0013 #include <string>
0014 #include <cmath>
0015
0016 using namespace edm;
0017
0018 HCALResponse::HCALResponse(const edm::ParameterSet& pset) {
0019   //switches
0020   debug = pset.getParameter<bool>("debug");
0021   usemip = pset.getParameter<bool>("usemip");
0022
0023   //values for "old" response parameterizations
0024   //--------------------------------------------------------------------
0025   RespPar[HCAL][0][0] = pset.getParameter<double>("HadronBarrelResolution_Stochastic");
0026   RespPar[HCAL][0][1] = pset.getParameter<double>("HadronBarrelResolution_Constant");
0027   RespPar[HCAL][0][2] = pset.getParameter<double>("HadronBarrelResolution_Noise");
0028
0029   RespPar[HCAL][1][0] = pset.getParameter<double>("HadronEndcapResolution_Stochastic");
0030   RespPar[HCAL][1][1] = pset.getParameter<double>("HadronEndcapResolution_Constant");
0031   RespPar[HCAL][1][2] = pset.getParameter<double>("HadronEndcapResolution_Noise");
0032
0033   RespPar[VFCAL][0][0] = pset.getParameter<double>("HadronForwardResolution_Stochastic");
0034   RespPar[VFCAL][0][1] = pset.getParameter<double>("HadronForwardResolution_Constant");
0035   RespPar[VFCAL][0][2] = pset.getParameter<double>("HadronForwardResolution_Noise");
0036
0037   RespPar[VFCAL][1][0] = pset.getParameter<double>("ElectronForwardResolution_Stochastic");
0038   RespPar[VFCAL][1][1] = pset.getParameter<double>("ElectronForwardResolution_Constant");
0039   RespPar[VFCAL][1][2] = pset.getParameter<double>("ElectronForwardResolution_Noise");
0040
0041   eResponseScale[0] = pset.getParameter<double>("eResponseScaleHB");
0042   eResponseScale[1] = pset.getParameter<double>("eResponseScaleHE");
0043   eResponseScale[2] = pset.getParameter<double>("eResponseScaleHF");
0044
0045   eResponsePlateau[0] = pset.getParameter<double>("eResponsePlateauHB");
0046   eResponsePlateau[1] = pset.getParameter<double>("eResponsePlateauHE");
0047   eResponsePlateau[2] = pset.getParameter<double>("eResponsePlateauHF");
0048
0049   eResponseExponent = pset.getParameter<double>("eResponseExponent");
0050   eResponseCoefficient = pset.getParameter<double>("eResponseCoefficient");
0051
0052   //pion parameters
0053   //--------------------------------------------------------------------
0054   //energy values
0055   maxHDe[0] = pset.getParameter<int>("maxHBe");
0056   maxHDe[1] = pset.getParameter<int>("maxHEe");
0057   maxHDe[2] = pset.getParameter<int>("maxHFe");
0058   maxHDe[3] = pset.getParameter<int>("maxHFlowe");
0059
0060   eGridHD[0] = pset.getParameter<vec1>("eGridHB");
0061   eGridHD[1] = pset.getParameter<vec1>("eGridHE");
0062   eGridHD[2] = pset.getParameter<vec1>("eGridHF");
0063   eGridHD[3] = pset.getParameter<vec1>("loweGridHF");
0064
0065   //region eta indices calculated from eta values
0066   etaStep = pset.getParameter<double>("etaStep");
0067   //eta boundaries
0068   HDeta[0] = abs((int)(pset.getParameter<double>("HBeta") / etaStep));
0069   HDeta[1] = abs((int)(pset.getParameter<double>("HEeta") / etaStep));
0070   HDeta[2] = abs((int)(pset.getParameter<double>("HFeta") / etaStep));
0071   HDeta[3] = abs((int)(pset.getParameter<double>("maxHDeta") / etaStep));  //add 1 because this is the max index
0072   //eta ranges
0073   maxHDetas[0] = HDeta[1] - HDeta[0];
0074   maxHDetas[1] = HDeta[2] - HDeta[1];
0075   maxHDetas[2] = HDeta[3] - HDeta[2];
0076
0077   //parameter info
0078   nPar = pset.getParameter<int>("nPar");
0079   parNames = pset.getParameter<std::vector<std::string> >("parNames");
0080   std::string detNames[] = {"_HB", "_HE", "_HF"};
0081   std::string mipNames[] = {"_mip", "_nomip", ""};
0082   std::string fraction = "f";
0083   //setup parameters (5D vector)
0084   parameters = vec5(nPar, vec4(3, vec3(3)));
0085   for (int p = 0; p < nPar; p++) {   //loop over parameters
0086     for (int m = 0; m < 3; m++) {    //loop over mip, nomip, total
0087       for (int d = 0; d < 3; d++) {  //loop over dets: HB, HE, HF
0088         //get from python
0089         std::string pname = parNames[p] + detNames[d] + mipNames[m];
0090         vec1 tmp = pset.getParameter<vec1>(pname);
0091
0092         //resize vector for energy range of det d
0093         parameters[p][m][d].resize(maxHDe[d]);
0094
0095         for (int i = 0; i < maxHDe[d]; i++) {  //loop over energy for det d
0096           //resize vector for eta range of det d
0097           parameters[p][m][d][i].resize(maxHDetas[d]);
0098
0099           for (int j = 0; j < maxHDetas[d]; j++) {  //loop over eta for det d
0100             //fill in parameters vector from python
0101             parameters[p][m][d][i][j] = tmp[i * maxHDetas[d] + j];
0102           }
0103         }
0104       }
0105     }
0106   }
0107   //set up Poisson parameters for low energy Hadrons in HF
0108   //----------------------------------------------------------------------
0109   PoissonParameters = vec3(4);
0110   std::string PoissonParName[] = {"mean_overall", "shift_overall", "mean_between", "shift_between"};
0111   for (int d = 0; d < 4; d++) {  //loop over Poisson parameteres
0112     vec1 tmp1 = pset.getParameter<vec1>(PoissonParName[d]);
0113     for (int i = 0; i < maxHDe[3]; i++) {  //loop over energy for low HF energy points
0114       PoissonParameters[d].resize(maxHDe[3]);
0115       for (int j = 0; j < maxHDetas[2]; j++) {  //loop over HF eta points
0116         PoissonParameters[d][i].resize(maxHDetas[2]);
0117         PoissonParameters[d][i][j] = tmp1[i * maxHDetas[2] + j];
0118       }
0119     }
0120   }
0121
0122   //MIP fraction fill in 3d vector
0123   ////--------------------------------------------------------------------
0124   mipfraction = vec3(3);
0125   for (int d = 0; d < 3; d++) {  //loop over dets: HB, HE, HF
0126     //get from python
0127     std::string mipname = fraction + mipNames[0] + detNames[d];
0128     vec1 tmp1 = pset.getParameter<vec1>(mipname);
0129     mipfraction[d].resize(maxHDe[d]);
0130     for (int i = 0; i < maxHDe[d]; i++) {  //loop over energy for det d
0131       //resize vector for eta range of det d
0132       mipfraction[d][i].resize(maxHDetas[d]);
0133       for (int j = 0; j < maxHDetas[d]; j++) {  //loop over eta for det d
0134         //fill in parameters vector from python
0135         mipfraction[d][i][j] = tmp1[i * maxHDetas[d] + j];
0136       }
0137     }
0138   }
0139
0140   // MUON probability histos for bin size = 0.25 GeV (0-10 GeV, 40 bins)
0141   //--------------------------------------------------------------------
0142   muStep = pset.getParameter<double>("muStep");
0143   maxMUe = pset.getParameter<int>("maxMUe");
0144   maxMUeta = pset.getParameter<int>("maxMUeta");
0145   maxMUbin = pset.getParameter<int>("maxMUbin");
0146   eGridMU = pset.getParameter<vec1>("eGridMU");
0147   etaGridMU = pset.getParameter<vec1>("etaGridMU");
0148   vec1 _responseMU[2] = {pset.getParameter<vec1>("responseMUBarrel"), pset.getParameter<vec1>("responseMUEndcap")};
0149
0150   //get muon region eta indices from the eta grid
0151   double _barrelMUeta = pset.getParameter<double>("barrelMUeta");
0152   double _endcapMUeta = pset.getParameter<double>("endcapMUeta");
0153   barrelMUeta = endcapMUeta = maxMUeta;
0154   for (int i = 0; i < maxMUeta; i++) {
0155     if (fabs(_barrelMUeta) <= etaGridMU[i]) {
0156       barrelMUeta = i;
0157       break;
0158     }
0159   }
0160   for (int i = 0; i < maxMUeta; i++) {
0161     if (fabs(_endcapMUeta) <= etaGridMU[i]) {
0162       endcapMUeta = i;
0163       break;
0164     }
0165   }
0166   int maxMUetas[] = {endcapMUeta - barrelMUeta, maxMUeta - endcapMUeta};
0167
0168   //initialize 3D vector
0169   responseMU = vec3(maxMUe, vec2(maxMUeta, vec1(maxMUbin, 0)));
0170
0171   //fill in 3D vector
0172   //(complementary cumulative distribution functions, from normalized response distributions)
0173   int loc, eta_loc;
0174   loc = eta_loc = -1;
0175   for (int i = 0; i < maxMUe; i++) {
0176     for (int j = 0; j < maxMUeta; j++) {
0177       //check location - barrel, endcap, or forward
0178       if (j == barrelMUeta) {
0179         loc = 0;
0180         eta_loc = barrelMUeta;
0181       } else if (j == endcapMUeta) {
0182         loc = 1;
0183         eta_loc = endcapMUeta;
0184       }
0185
0186       for (int k = 0; k < maxMUbin; k++) {
0187         responseMU[i][j][k] = _responseMU[loc][i * maxMUetas[loc] * maxMUbin + (j - eta_loc) * maxMUbin + k];
0188
0189         if (debug) {
0190           //cout.width(6);
0191           LogInfo("FastCalorimetry") << " responseMU " << i << " " << j << " " << k << " = " << responseMU[i][j][k]
0192                                      << std::endl;
0193         }
0194       }
0195     }
0196   }
0197
0198   // values for EM response in HF
0199   //--------------------------------------------------------------------
0200   maxEMe = pset.getParameter<int>("maxEMe");
0201   maxEMeta = maxHDetas[2];
0202   respFactorEM = pset.getParameter<double>("respFactorEM");
0203   eGridEM = pset.getParameter<vec1>("eGridEM");
0204
0205   // e-gamma mean response and sigma in HF
0206   vec1 _meanEM = pset.getParameter<vec1>("meanEM");
0207   vec1 _sigmaEM = pset.getParameter<vec1>("sigmaEM");
0208
0209   //fill in 2D vectors (w/ correction factor applied)
0210   meanEM = vec2(maxEMe, vec1(maxEMeta, 0));
0211   sigmaEM = vec2(maxEMe, vec1(maxEMeta, 0));
0212   for (int i = 0; i < maxEMe; i++) {
0213     for (int j = 0; j < maxEMeta; j++) {
0214       meanEM[i][j] = respFactorEM * _meanEM[i * maxEMeta + j];
0215       sigmaEM[i][j] = respFactorEM * _sigmaEM[i * maxEMeta + j];
0216     }
0217   }
0218
0219   // HF correction for SL
0220   //---------------------
0221   maxEta = pset.getParameter<int>("maxEta");
0222   maxEne = pset.getParameter<int>("maxEne");
0223   energyHF = pset.getParameter<vec1>("energyHF");
0224   vec1 _corrHFgEm = pset.getParameter<vec1>("corrHFgEm");
0225   vec1 _corrHFgHad = pset.getParameter<vec1>("corrHFgHad");
0226   vec1 _corrHFhEm = pset.getParameter<vec1>("corrHFhEm");
0227   vec1 _corrHFhHad = pset.getParameter<vec1>("corrHFhHad");
0228   corrHFem = vec1(maxEta, 0);
0229   corrHFhad = vec1(maxEta, 0);
0230
0231   // initialize 2D vector corrHFgEm[energy][eta]
0232   corrHFgEm = vec2(maxEne, vec1(maxEta, 0));
0233   corrHFgHad = vec2(maxEne, vec1(maxEta, 0));
0234   corrHFhEm = vec2(maxEne, vec1(maxEta, 0));
0235   corrHFhHad = vec2(maxEne, vec1(maxEta, 0));
0236   // Fill
0237   for (int i = 0; i < maxEne; i++) {
0238     for (int j = 0; j < maxEta; j++) {
0239       corrHFgEm[i][j] = _corrHFgEm[i * maxEta + j];
0240       corrHFgHad[i][j] = _corrHFgHad[i * maxEta + j];
0241       corrHFhEm[i][j] = _corrHFhEm[i * maxEta + j];
0242       corrHFhHad[i][j] = _corrHFhHad[i * maxEta + j];
0243     }
0244   }
0245 }
0246
0247 double HCALResponse::getMIPfraction(double energy, double eta) {
0248   int ieta = abs((int)(eta / etaStep));
0249   int ie = -1;
0250   //check eta and det
0251   int det = getDet(ieta);
0252   int deta = ieta - HDeta[det];
0253   if (deta >= maxHDetas[det])
0254     deta = maxHDetas[det] - 1;
0255   else if (deta < 0)
0256     deta = 0;
0257   //find energy range
0258   for (int i = 0; i < maxHDe[det]; i++) {
0259     if (energy < eGridHD[det][i]) {
0260       if (i == 0)
0261         return mipfraction[det][0][deta];  // less than minimal - the first value is used instead of extrapolating
0262       else
0263         ie = i - 1;
0264       break;
0265     }
0266   }
0267   if (ie == -1)
0268     return mipfraction[det][maxHDe[det] - 1]
0269                       [deta];  // more than maximal - the last value is used instead of extrapolating
0270   double y1, y2;
0271   double x1 = eGridHD[det][ie];
0272   double x2 = eGridHD[det][ie + 1];
0273   y1 = mipfraction[det][ie][deta];
0274   y2 = mipfraction[det][ie + 1][deta];
0275   double mean = 0;
0276   mean = (y1 * (x2 - energy) + y2 * (energy - x1)) / (x2 - x1);
0277   return mean;
0278 }
0279
0280 double HCALResponse::responseHCAL(
0281     int _mip, double energy, double eta, int partype, RandomEngineAndDistribution const* random) {
0282   int ieta = abs((int)(eta / etaStep));
0283   int ie = -1;
0284
0285   int mip;
0286   if (usemip)
0287     mip = _mip;
0288   else
0289     mip = 2;  //ignore mip, use only overall (mip + nomip) parameters
0290
0291   double mean = 0;
0292
0293   // e/gamma in HF
0294   if (partype == 0) {
0295     //check eta
0296     ieta -= HDeta[2];  // HF starts at ieta=30 till ieta=51
0297     // but resp.vector from index=0 through 20
0298     if (ieta >= maxEMeta)
0299       ieta = maxEMeta - 1;
0300     else if (ieta < 0)
0301       ieta = 0;
0302
0303     //find energy range
0304     for (int i = 0; i < maxEMe; i++) {
0305       if (energy < eGridEM[i]) {
0306         if (i == 0)
0307           ie = 0;  // less than minimal - back extrapolation with the 1st interval
0308         else
0309           ie = i - 1;
0310         break;
0311       }
0312     }
0313     if (ie == -1)
0314       ie = maxEMe - 2;  // more than maximum - extrapolation with last interval
0315
0316     //do smearing
0317     mean = interEM(energy, ie, ieta, random);
0318   }
0319
0320   // hadrons
0321   else if (partype == 1) {
0322     //check eta and det
0323     int det = getDet(ieta);
0324     int deta = ieta - HDeta[det];
0325     if (deta >= maxHDetas[det])
0326       deta = maxHDetas[det] - 1;
0327     else if (deta < 0)
0328       deta = 0;
0329
0330     //find energy range
0331     for (int i = 0; i < maxHDe[det]; i++) {
0332       if (energy < eGridHD[det][i]) {
0333         if (i == 0)
0334           ie = 0;  // less than minimal - back extrapolation with the 1st interval
0335         else
0336           ie = i - 1;
0337         break;
0338       }
0339     }
0340     if (ie == -1)
0341       ie = maxHDe[det] - 2;  // more than maximum - extrapolation with last interval
0342
0343     //different energy smearing for low energy hadrons in HF
0344     if (det == 2 && energy < 20 && deta > 5) {
0345       for (int i = 0; i < maxHDe[3]; i++) {
0346         if (energy < eGridHD[3][i]) {
0347           if (i == 0)
0348             ie = 0;  // less than minimal - back extrapolation with the 1st interval
0349           else
0350             ie = i - 1;
0351           break;
0352         }
0353       }
0354     }
0355     //do smearing
0356     mean = interHD(mip, energy, ie, deta, det, random);
0357   }
0358
0359   // muons
0360   else if (partype == 2) {
0361     //check eta
0362     ieta = maxMUeta;
0363     for (int i = 0; i < maxMUeta; i++) {
0364       if (fabs(eta) < etaGridMU[i]) {
0365         ieta = i;
0366         break;
0367       }
0368     }
0369     if (ieta < 0)
0370       ieta = 0;
0371
0372     if (ieta < maxMUeta) {  // HB-HE
0373       //find energy range
0374       for (int i = 0; i < maxMUe; i++) {
0375         if (energy < eGridMU[i]) {
0376           if (i == 0)
0377             ie = 0;  // less than minimal - back extrapolation with the 1st interval
0378           else
0379             ie = i - 1;
0380           break;
0381         }
0382       }
0383       if (ie == -1)
0384         ie = maxMUe - 2;  // more than maximum - extrapolation with last interval
0385
0386       //do smearing
0387       mean = interMU(energy, ie, ieta, random);
0388       if (mean > energy)
0389         mean = energy;
0390     }
0391   }
0392
0393   // debugging
0394   if (debug) {
0395     //  cout.width(6);
0396     LogInfo("FastCalorimetry") << std::endl
0397                                << " HCALResponse::responseHCAL, partype = " << partype << " E, eta = " << energy << " "
0398                                << eta << "  mean = " << mean << std::endl;
0399   }
0400
0401   return mean;
0402 }
0403
0404 double HCALResponse::interMU(double e, int ie, int ieta, RandomEngineAndDistribution const* random) {
0405   double x = random->flatShoot();
0406
0407   int bin1 = maxMUbin;
0408   for (int i = 0; i < maxMUbin; i++) {
0409     if (x > responseMU[ie][ieta][i]) {
0410       bin1 = i - 1;
0411       break;
0412     }
0413   }
0414   int bin2 = maxMUbin;
0415   for (int i = 0; i < maxMUbin; i++) {
0416     if (x > responseMU[ie + 1][ieta][i]) {
0417       bin2 = i - 1;
0418       break;
0419     }
0420   }
0421
0422   double x1 = eGridMU[ie];
0423   double x2 = eGridMU[ie + 1];
0424   double y1 = (bin1 + random->flatShoot()) * muStep;
0425   double y2 = (bin2 + random->flatShoot()) * muStep;
0426
0427   if (debug) {
0428     //  cout.width(6);
0429     LogInfo("FastCalorimetry") << std::endl
0430                                << " HCALResponse::interMU  " << std::endl
0431                                << " x, x1-x2, y1-y2 = " << e << ", " << x1 << "-" << x2 << " " << y1 << "-" << y2
0432                                << std::endl;
0433   }
0434
0435   //linear interpolation
0436   double mean = (y1 * (x2 - e) + y2 * (e - x1)) / (x2 - x1);
0437
0438   if (debug) {
0439     //cout.width(6);
0440     LogInfo("FastCalorimetry") << std::endl
0441                                << " HCALResponse::interMU " << std::endl
0442                                << " e, ie, ieta = " << e << " " << ie << " " << ieta << std::endl
0443                                << " response  = " << mean << std::endl;
0444   }
0445
0446   return mean;
0447 }
0448
0449 double HCALResponse::interHD(int mip, double e, int ie, int ieta, int det, RandomEngineAndDistribution const* random) {
0450   double x1, x2;
0451   double y1, y2;
0452   if (det == 2)
0453     mip = 2;  //ignore mip status for HF
0454   double mean = 0;
0455   vec1 pars(nPar, 0);
0456
0457   // for ieta < 5 there is overlap between HE and HF, and measurement comes from HE
0458   if (det == 2 && ieta > 5 && e < 20) {
0459     for (int p = 0; p < 4; p++) {
0460       y1 = PoissonParameters[p][ie][ieta];
0461       y2 = PoissonParameters[p][ie + 1][ieta];
0462       if (e > 5) {
0463         x1 = eGridHD[det + 1][ie];
0464         x2 = eGridHD[det + 1][ie + 1];
0465         pars[p] = (y1 * (x2 - e) + y2 * (e - x1)) / (x2 - x1);
0466       } else
0467         pars[p] = y1;
0468     }
0469     mean =
0470         random->poissonShoot((int(PoissonShootNoNegative(pars[0], pars[1], random)) +
0471                               (int(PoissonShootNoNegative(pars[2], pars[3], random))) / 4 + random->flatShoot() / 4) *
0472                              6) /
0473         (0.3755 * 6);
0474   }
0475
0476   else {
0477     x1 = eGridHD[det][ie];
0478     x2 = eGridHD[det][ie + 1];
0479
0480     //calculate all parameters
0481     for (int p = 0; p < nPar; p++) {
0482       y1 = parameters[p][mip][det][ie][ieta];
0483       y2 = parameters[p][mip][det][ie + 1][ieta];
0484
0485       //par-specific checks
0486       double custom = 0;
0487       bool use_custom = false;
0488
0489       //do not let mu or sigma get extrapolated below zero for low energies
0490       //especially important for HF since extrapolation is used for E < 15 GeV
0491       if ((p == 0 || p == 1) && e < x1) {
0492         double tmp = (y1 * x2 - y2 * x1) / (x2 - x1);  //extrapolate down to e=0
0493         if (tmp < 0) {                                 //require mu,sigma > 0 for E > 0
0494           custom = y1 * e / x1;
0495           use_custom = true;
0496         }
0497       }
0498       //tail parameters have lower bounds - never extrapolate down
0499       else if ((p == 2 || p == 3 || p == 4 || p == 5)) {
0500         if (e < x1 && y1 < y2) {
0501           custom = y1;
0502           use_custom = true;
0503         } else if (e > x2 && y2 < y1) {
0504           custom = y2;
0505           use_custom = true;
0506         }
0507       }
0508
0509       //linear interpolation
0510       if (use_custom)
0511         pars[p] = custom;
0512       else
0513         pars[p] = (y1 * (x2 - e) + y2 * (e - x1)) / (x2 - x1);
0514     }
0515
0516     //random smearing
0517     if (nPar == 6)
0518       mean = cballShootNoNegative(pars[0], pars[1], pars[2], pars[3], pars[4], pars[5], random);
0519     else if (nPar == 2)
0520       mean = gaussShootNoNegative(pars[0], pars[1], random);  //gaussian fallback
0521   }
0522   return mean;
0523 }
0524
0525 double HCALResponse::interEM(double e, int ie, int ieta, RandomEngineAndDistribution const* random) {
0526   double y1 = meanEM[ie][ieta];
0527   double y2 = meanEM[ie + 1][ieta];
0528   double x1 = eGridEM[ie];
0529   double x2 = eGridEM[ie + 1];
0530
0531   if (debug) {
0532     //  cout.width(6);
0533     LogInfo("FastCalorimetry") << std::endl
0534                                << " HCALResponse::interEM mean " << std::endl
0535                                << " x, x1-x2, y1-y2 = " << e << ", " << x1 << "-" << x2 << " " << y1 << "-" << y2
0536                                << std::endl;
0537   }
0538
0539   //linear interpolation
0540   double mean = (y1 * (x2 - e) + y2 * (e - x1)) / (x2 - x1);
0541
0542   y1 = sigmaEM[ie][ieta];
0543   y2 = sigmaEM[ie + 1][ieta];
0544
0545   if (debug) {
0546     //  cout.width(6);
0547     LogInfo("FastCalorimetry") << std::endl
0548                                << " HCALResponse::interEM sigma" << std::endl
0549                                << " x, x1-x2, y1-y2 = " << e << ", " << x1 << "-" << x2 << " " << y1 << "-" << y2
0550                                << std::endl;
0551   }
0552
0553   //linear interpolation
0554   double sigma = (y1 * (x2 - e) + y2 * (e - x1)) / (x2 - x1);
0555
0556   //random smearing
0557   double rndm = gaussShootNoNegative(mean, sigma, random);
0558
0559   return rndm;
0560 }
0561
0562 // Old parameterization of the calo response to hadrons
0563 double HCALResponse::getHCALEnergyResponse(double e, int hit, RandomEngineAndDistribution const* random) {
0564   //response
0565   double s = eResponseScale[hit];
0566   double n = eResponseExponent;
0567   double p = eResponsePlateau[hit];
0568   double c = eResponseCoefficient;
0569
0570   double response = e * p / (1 + c * exp(n * log(s / e)));
0571
0572   if (response < 0.)
0573     response = 0.;
0574
0575   //resolution
0576   double resolution;
0577   if (hit == hcforward)
0578     resolution =
0579         e * sqrt(RespPar[VFCAL][1][0] * RespPar[VFCAL][1][0] / e + RespPar[VFCAL][1][1] * RespPar[VFCAL][1][1]);
0580   else
0581     resolution =
0582         e * sqrt(RespPar[HCAL][hit][0] * RespPar[HCAL][hit][0] / (e) + RespPar[HCAL][hit][1] * RespPar[HCAL][hit][1]);
0583
0584   //random smearing
0585   double rndm = gaussShootNoNegative(response, resolution, random);
0586
0587   return rndm;
0588 }
0589
0590 //find subdet and eta offset
0591 int HCALResponse::getDet(int ieta) {
0592   int d;
0593   for (d = 0; d < 2; d++) {
0594     if (ieta < HDeta[d + 1]) {
0595       break;
0596     }
0597   }
0598   return d;
0599 }
0600
0601 // Remove (most) hits with negative energies
0602 double HCALResponse::gaussShootNoNegative(double e, double sigma, RandomEngineAndDistribution const* random) {
0603   double out = random->gaussShoot(e, sigma);
0604   if (e >= 0.) {
0605     while (out < 0.)
0606       out = random->gaussShoot(e, sigma);
0607   }
0608   //else give up on re-trying, otherwise too much time can be lost before emeas comes out positive
0609
0610   return out;
0611 }
0612
0613 // Remove (most) hits with negative energies
0614 double HCALResponse::cballShootNoNegative(
0615     double mu, double sigma, double aL, double nL, double aR, double nR, RandomEngineAndDistribution const* random) {
0616   double out = cball.shoot(mu, sigma, aL, nL, aR, nR, random);
0617   if (mu >= 0.) {
0618     while (out < 0.)
0619       out = cball.shoot(mu, sigma, aL, nL, aR, nR, random);
0620   }
0621   //else give up on re-trying, otherwise too much time can be lost before emeas comes out positive
0622
0623   return out;
0624 }
0625 double HCALResponse::PoissonShootNoNegative(double e, double sigma, RandomEngineAndDistribution const* random) {
0626   double out = -1;
0627   while (out < 0.) {
0628     out = random->poissonShoot(e);
0629     out = out + sigma;
0630   }
0631   return out;
0632 }
0633
0634 void HCALResponse::correctHF(double ee, int type) {
0635   int jmin = 0;
0636   for (int i = 0; i < maxEne; i++) {
0637     if (ee >= energyHF[i])
0638       jmin = i;
0639   }
0640
0641   double x1, x2;
0642   double y1em, y2em;
0643   double y1had, y2had;
0644   for (int i = 0; i < maxEta; ++i) {
0645     if (ee < energyHF[0]) {
0646       if (abs(type) == 11 || abs(type) == 22) {
0647         corrHFem[i] = corrHFgEm[0][i];
0648         corrHFhad[i] = corrHFgHad[0][i];
0649       } else {
0650         corrHFem[i] = corrHFhEm[0][i];
0651         corrHFhad[i] = corrHFhHad[0][i];
0652       }
0653     } else if (jmin >= maxEne - 1) {
0654       if (abs(type) == 11 || abs(type) == 22) {
0655         corrHFem[i] = corrHFgEm[maxEta][i];
0656         corrHFhad[i] = corrHFgHad[maxEta][i];
0657       } else {
0658         corrHFem[i] = corrHFhEm[maxEta][i];
0659         corrHFhad[i] = corrHFhHad[maxEta][i];
0660       }
0661     } else {
0662       x1 = energyHF[jmin];
0663       x2 = energyHF[jmin + 1];
0664       if (abs(type) == 11 || abs(type) == 22) {
0665         y1em = corrHFgEm[jmin][i];
0666         y2em = corrHFgEm[jmin + 1][i];
0667         y1had = corrHFgHad[jmin][i];
0668         y2had = corrHFgHad[jmin + 1][i];
0669       } else {
0670         y1em = corrHFhEm[jmin][i];
0671         y2em = corrHFhEm[jmin + 1][i];
0672         y1had = corrHFhHad[jmin][i];
0673         y2had = corrHFhHad[jmin + 1][i];
0674       }
0675       corrHFem[i] = y1em + (ee - x1) * ((y2em - y1em) / (x2 - x1));
0676       corrHFhad[i] = y1had + (ee - x1) * ((y2had - y1had) / (x2 - x1));
0677     }
0678   }
0679 }