CastorCalib/src/CastorPulseContainmentCorrection.cc

0001 #include "CalibCalorimetry/CastorCalib/interface/CastorPulseContainmentCorrection.h"
0002 #include "CalibCalorimetry/CastorCalib/interface/CastorTimeSlew.h"
0003 #include "CalibCalorimetry/CastorCalib/interface/CastorPulseShapes.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include <algorithm>  // for "max","min"
0006 #include <cmath>
0007 #include <iostream>
0008
0009 #include "CalibCalorimetry/HcalAlgos/interface/genlkupmap.h"
0010
0011 template <>
0012 RecoFCcorFactorAlgo<CastorPulseShapes::Shape>::RecoFCcorFactorAlgo(int num_samples, double fixedphase_ns) {
0013   fixedphasens_ = fixedphase_ns;
0014   CastorPulseShapes shapes;
0015   shape_ = shapes.castorShape();
0016   const int binsize_ns = 25;
0017
0018   // First set up controlling parameters for calculating the correction factor:
0019   // Integration window size...
0020   //
0021   integrationwindowns_ = (double)(binsize_ns * num_samples);
0022
0023   // First find the point at which time bin "1" exceeds time bin "0",
0024   // and call that point "time 0".
0025   //
0026   for (int shift_ns = 0; shift_ns < binsize_ns; shift_ns++) {
0027     // Digitize by integrating to find all time sample
0028     // bin values for this shift.
0029     //
0030     double tmin = -(double)shift_ns;
0031     double bin0val = (double)shape_.integrate(tmin, tmin + binsize_ns);
0032     double bin1val = (double)shape_.integrate(tmin + binsize_ns, tmin + 2 * binsize_ns);
0033
0034 #if 0
0035     char s[80];
0036     sprintf (s, "%7.3f %8.5f %8.5f\n", tmin, bin0val, bin1val);
0037     cout << s;
0038 #endif
0039
0040     if (bin1val > bin0val) {
0041       time0shiftns_ = shift_ns;
0042       break;
0043     }
0044   }
0045
0046 #if 0
0047   cout << "time0shiftns_ = " << time0shiftns_ << endl;
0048 #endif
0049 }
0050
0051 template <>
0052 std::pair<double, double> RecoFCcorFactorAlgo<CastorPulseShapes::Shape>::calcpair(double truefc) {
0053   double timeslew_ns = CastorTimeSlew::delay(std::max(0.0, (double)truefc), CastorTimeSlew::Medium);
0054   double shift_ns = fixedphasens_ - time0shiftns_ + timeslew_ns;
0055
0056   double tmin = -shift_ns;
0057   double tmax = tmin + integrationwindowns_;
0058
0059   double integral = shape_.integrate(tmin, tmax);
0060   double corfactor = 1.0 / integral;
0061   double recofc = (double)truefc * integral;
0062
0063 #if 0
0064   char s[80];
0065   sprintf (s, "%8.2f %8.4f %8.4f %8.5f %8.5f %8.5f ",
0066        truefc, tmin, tmax, integral, corfactor, recofc);
0067   cout << s;
0068 #endif
0069
0070   std::pair<double, double> thepair(recofc, corfactor);
0071   return thepair;
0072 }
0073
0074 ///Generate energy correction factors based on a predetermined phase of the hit + time slew
0075 //
0076 CastorPulseContainmentCorrection::CastorPulseContainmentCorrection(int num_samples,
0077                                                                    float fixedphase_ns,
0078                                                                    float max_fracerror) {
0079   RecoFCcorFactorAlgo<CastorPulseShapes::Shape> corFalgo(num_samples, (double)fixedphase_ns);
0080
0081   // Generate lookup map for the correction function, never exceeding
0082   // a maximum fractional error for lookups.
0083   //
0084   genlkupmap<RecoFCcorFactorAlgo<CastorPulseShapes::Shape> >(1.0,
0085                                                              5000.0f,        // generation domain
0086                                                              max_fracerror,  // maximum fractional error
0087                                                              1.0,            // min_xstep = minimum true fC increment
0088                                                              corFalgo,
0089                                                              mCorFactors_);  // return lookup map
0090 }
0091
0092 double CastorPulseContainmentCorrection::getCorrection(double fc_ampl) const {
0093   double correction;
0094
0095   std::map<double, double>::const_iterator fcupper, fclower;
0096
0097   fcupper = mCorFactors_.upper_bound(fc_ampl);
0098   fclower = fcupper;
0099   fclower--;
0100
0101   if (fcupper == mCorFactors_.end()) {
0102     correction = fclower->second;
0103   } else if (fcupper == mCorFactors_.begin()) {
0104     correction = fcupper->second;
0105   } else {
0106     if (fabs(fclower->first - fc_ampl) < fabs(fcupper->first - fc_ampl))
0107       correction = fclower->second;
0108     else
0109       correction = fcupper->second;
0110   }
0111
0112 #if 0
0113   char s[80];
0114   sprintf (s, "%7.1f (%8.5f %8.5f) (%8.5f %8.5f) %8.5f",
0115        fc_ampl,
0116        fclower->first, fclower->second,
0117        fcupper->first, fcupper->second,
0118        correction);
0119   cout << s << endl;
0120 #endif
0121
0122   return correction;
0123 }