CastorReco/src/CastorSimpleRecAlgo.cc

0001 #include "RecoLocalCalo/CastorReco/interface/CastorSimpleRecAlgo.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include "CalibCalorimetry/CastorCalib/interface/CastorTimeSlew.h"
0004 #include "DataFormats/METReco/interface/HcalCaloFlagLabels.h"
0005 #include <algorithm>  // for "max"
0006 #include <cmath>
0007
0008 constexpr double MaximumFractionalError = 0.0005;  // 0.05% error allowed from this source
0009
0010 CastorSimpleRecAlgo::CastorSimpleRecAlgo(
0011     int firstSample, int samplesToAdd, bool correctForTimeslew, bool correctForPulse, float phaseNS)
0012     : firstSample_(firstSample), samplesToAdd_(samplesToAdd), correctForTimeslew_(correctForTimeslew) {
0013   if (correctForPulse)
0014     pulseCorr_ = std::make_unique<CastorPulseContainmentCorrection>(samplesToAdd_, phaseNS, MaximumFractionalError);
0015 }
0016
0017 CastorSimpleRecAlgo::CastorSimpleRecAlgo(int firstSample, int samplesToAdd)
0018     : firstSample_(firstSample), samplesToAdd_(samplesToAdd), correctForTimeslew_(false) {}
0019
0020 ///Timeshift correction for HPDs based on the position of the peak ADC measurement.
0021 ///  Allows for an accurate determination of the relative phase of the pulse shape from
0022 ///  the HPD.  Calculated based on a weighted sum of the -1,0,+1 samples relative to the peak
0023 ///  as follows:  wpksamp = (0*sample[0] + 1*sample[1] + 2*sample[2]) / (sample[0] + sample[1] + sample[2])
0024 ///  where sample[1] is the maximum ADC sample value.
0025 //static float timeshift_ns_hbheho(float wpksamp);
0026
0027 ///Same as above, but for the HF PMTs.
0028 static float timeshift_ns_hf(float wpksamp);
0029
0030 namespace CastorSimpleRecAlgoImpl {
0031   template <class Digi, class RecHit>
0032   inline RecHit reco(const Digi& digi,
0033                      const CastorCoder& coder,
0034                      const CastorCalibrations& calibs,
0035                      int ifirst,
0036                      int n,
0037                      bool slewCorrect,
0038                      const CastorPulseContainmentCorrection* corr,
0039                      CastorTimeSlew::BiasSetting slewFlavor) {
0040     CaloSamples tool;
0041     coder.adc2fC(digi, tool);
0042
0043     double ampl = 0;
0044     int maxI = -1;
0045     double maxA = -1e10;
0046     float ta = 0;
0047     double fc_ampl = 0;
0048     for (int i = ifirst; i < tool.size() && i < n + ifirst; i++) {
0049       int capid = digi[i].capid();
0050       ta = (tool[i] - calibs.pedestal(capid));  // pedestal subtraction
0051       fc_ampl += ta;
0052       ta *= calibs.gain(capid);  // fC --> GeV
0053       ampl += ta;
0054       if (ta > maxA) {
0055         maxA = ta;
0056         maxI = i;
0057       }
0058     }
0059
0060     float time = -9999;
0061     ////Cannot calculate time value with max ADC sample at first or last position in window....
0062     if (maxI == 0 || maxI == (tool.size() - 1)) {
0063       // supress warning and use dummy time value for 2009 data
0064       /*
0065       edm::LogWarning("HCAL/CASTOR Pulse") << "CastorSimpleRecAlgo::reconstruct :"
0066                            << " Invalid max amplitude position, "
0067                            << " max Amplitude: "<< maxI
0068                            << " first: "<<ifirst
0069                            << " last: "<<(tool.size()-1)
0070                            << std::endl;
0071     */
0072     } else {
0073       maxA = fabs(maxA);
0074       int capid = digi[maxI - 1].capid();
0075       float t0 = fabs((tool[maxI - 1] - calibs.pedestal(capid)) * calibs.gain(capid));
0076       capid = digi[maxI + 1].capid();
0077       float t2 = fabs((tool[maxI + 1] - calibs.pedestal(capid)) * calibs.gain(capid));
0078       float wpksamp = (t0 + maxA + t2);
0079       if (wpksamp != 0)
0080         wpksamp = (maxA + 2.0 * t2) / wpksamp;
0081       time = (maxI - digi.presamples()) * 25.0 + timeshift_ns_hf(wpksamp);
0082
0083       if (corr != nullptr) {
0084         // Apply phase-based amplitude correction:
0085         ampl *= corr->getCorrection(fc_ampl);
0086         //      std::cout << fc_ampl << " --> " << corr->getCorrection(fc_ampl) << std::endl;
0087       }
0088
0089       if (slewCorrect)
0090         time -= CastorTimeSlew::delay(std::max(1.0, fc_ampl), slewFlavor);
0091     }
0092     return RecHit(digi.id(), ampl, time);
0093   }
0094
0095   // returns TRUE if ADC count is >= maxADCvalue
0096   template <class Digi>
0097   inline bool isSaturated(const Digi& digi, const int& maxADCvalue, int ifirst, int n) {
0098     for (int i = ifirst; i < digi.size() && i < n + ifirst; i++) {
0099       if (digi[i].adc() >= maxADCvalue)
0100         return true;
0101     }
0102
0103     return false;
0104   }
0105
0106   //++++ Saturation Correction +++++
0107   // returns TRUE when saturation correction was done
0108   template <class Digi, class RecHit>
0109   inline bool corrSaturation(RecHit& rechit,
0110                              const CastorCoder& coder,
0111                              const CastorCalibrations& calibs,
0112                              const Digi& digi,
0113                              const int& maxADCvalue,
0114                              const double& satCorrConst,
0115                              int ifirst,
0116                              int n) {
0117     // correction works only with 2 used TS
0118     if (n != 2)
0119       return false;
0120
0121     // to avoid segmentation violation
0122     if (ifirst + n > digi.size())
0123       return false;
0124
0125     // look if first TS is saturated not the second one
0126     // in case the second one is saturated do no correction
0127     if (digi[ifirst].adc() >= maxADCvalue && digi[ifirst + 1].adc() < maxADCvalue) {
0128       float ampl = 0;
0129
0130       CaloSamples tool;
0131       coder.adc2fC(digi, tool);
0132
0133       // get energy of first TS
0134       int capid = digi[ifirst].capid();
0135       float ta = tool[ifirst] - calibs.pedestal(capid);  // pedestal subtraction
0136       float ta1 = ta * calibs.gain(capid);               // fC --> GeV
0137
0138       // get energy of second TS
0139       capid = digi[ifirst + 1].capid();
0140       ta = tool[ifirst + 1] - calibs.pedestal(capid);  // pedestal subtraction
0141       float ta2 = ta * calibs.gain(capid);             // fC --> GeV
0142
0143       // use second TS to calc the first one
0144       // check if recalculated TS ampl is
0145       // realy greater than the saturated value
0146       if (ta2 / satCorrConst <= ta1)
0147         return false;
0148
0149       // ampl = TS5 + TS4 => TS4 = TS5 / TSratio
0150       // ampl = TS5 + TS5 / TSratio
0151       ampl = ta2 + ta2 / satCorrConst;
0152
0153       // reset energy of rechit
0154       rechit.setEnergy(ampl);
0155
0156       return true;
0157     }
0158
0159     return false;
0160   }
0161 }  // namespace CastorSimpleRecAlgoImpl
0162
0163 CastorRecHit CastorSimpleRecAlgo::reconstruct(const CastorDataFrame& digi,
0164                                               const CastorCoder& coder,
0165                                               const CastorCalibrations& calibs) const {
0166   return CastorSimpleRecAlgoImpl::reco<CastorDataFrame, CastorRecHit>(
0167       digi, coder, calibs, firstSample_, samplesToAdd_, false, nullptr, CastorTimeSlew::Fast);
0168 }
0169
0170 void CastorSimpleRecAlgo::checkADCSaturation(CastorRecHit& rechit,
0171                                              const CastorDataFrame& digi,
0172                                              const int& maxADCvalue) const {
0173   if (CastorSimpleRecAlgoImpl::isSaturated<CastorDataFrame>(digi, maxADCvalue, firstSample_, samplesToAdd_))
0174     rechit.setFlagField(1, HcalCaloFlagLabels::ADCSaturationBit);
0175 }
0176
0177 //++++ Saturation Correction +++++
0178 void CastorSimpleRecAlgo::recoverADCSaturation(CastorRecHit& rechit,
0179                                                const CastorCoder& coder,
0180                                                const CastorCalibrations& calibs,
0181                                                const CastorDataFrame& digi,
0182                                                const int& maxADCvalue,
0183                                                const double& satCorrConst) const {
0184   if (CastorSimpleRecAlgoImpl::corrSaturation<CastorDataFrame, CastorRecHit>(
0185           rechit, coder, calibs, digi, maxADCvalue, satCorrConst, firstSample_, samplesToAdd_))
0186     // use empty flag bit for recording saturation correction
0187     // see also CMSSW/DataFormats/METReco/interface/HcalCaloFlagLabels.h
0188     rechit.setFlagField(1, HcalCaloFlagLabels::UserDefinedBit0);
0189 }
0190
0191 // timeshift implementation
0192
0193 static const float wpksamp0_hf = 0.500635;
0194 static const float scale_hf = 0.999301;
0195 static const int num_bins_hf = 100;
0196
0197 static const float actual_ns_hf[num_bins_hf] = {
0198     0.00000,   // 0.000-0.010
0199     0.03750,   // 0.010-0.020
0200     0.07250,   // 0.020-0.030
0201     0.10750,   // 0.030-0.040
0202     0.14500,   // 0.040-0.050
0203     0.18000,   // 0.050-0.060
0204     0.21500,   // 0.060-0.070
0205     0.25000,   // 0.070-0.080
0206     0.28500,   // 0.080-0.090
0207     0.32000,   // 0.090-0.100
0208     0.35500,   // 0.100-0.110
0209     0.39000,   // 0.110-0.120
0210     0.42500,   // 0.120-0.130
0211     0.46000,   // 0.130-0.140
0212     0.49500,   // 0.140-0.150
0213     0.53000,   // 0.150-0.160
0214     0.56500,   // 0.160-0.170
0215     0.60000,   // 0.170-0.180
0216     0.63500,   // 0.180-0.190
0217     0.67000,   // 0.190-0.200
0218     0.70750,   // 0.200-0.210
0219     0.74250,   // 0.210-0.220
0220     0.78000,   // 0.220-0.230
0221     0.81500,   // 0.230-0.240
0222     0.85250,   // 0.240-0.250
0223     0.89000,   // 0.250-0.260
0224     0.92750,   // 0.260-0.270
0225     0.96500,   // 0.270-0.280
0226     1.00250,   // 0.280-0.290
0227     1.04250,   // 0.290-0.300
0228     1.08250,   // 0.300-0.310
0229     1.12250,   // 0.310-0.320
0230     1.16250,   // 0.320-0.330
0231     1.20500,   // 0.330-0.340
0232     1.24500,   // 0.340-0.350
0233     1.29000,   // 0.350-0.360
0234     1.33250,   // 0.360-0.370
0235     1.38000,   // 0.370-0.380
0236     1.42500,   // 0.380-0.390
0237     1.47500,   // 0.390-0.400
0238     1.52500,   // 0.400-0.410
0239     1.57750,   // 0.410-0.420
0240     1.63250,   // 0.420-0.430
0241     1.69000,   // 0.430-0.440
0242     1.75250,   // 0.440-0.450
0243     1.82000,   // 0.450-0.460
0244     1.89250,   // 0.460-0.470
0245     1.97500,   // 0.470-0.480
0246     2.07250,   // 0.480-0.490
0247     2.20000,   // 0.490-0.500
0248     19.13000,  // 0.500-0.510
0249     21.08750,  // 0.510-0.520
0250     21.57750,  // 0.520-0.530
0251     21.89000,  // 0.530-0.540
0252     22.12250,  // 0.540-0.550
0253     22.31000,  // 0.550-0.560
0254     22.47000,  // 0.560-0.570
0255     22.61000,  // 0.570-0.580
0256     22.73250,  // 0.580-0.590
0257     22.84500,  // 0.590-0.600
0258     22.94750,  // 0.600-0.610
0259     23.04250,  // 0.610-0.620
0260     23.13250,  // 0.620-0.630
0261     23.21500,  // 0.630-0.640
0262     23.29250,  // 0.640-0.650
0263     23.36750,  // 0.650-0.660
0264     23.43750,  // 0.660-0.670
0265     23.50500,  // 0.670-0.680
0266     23.57000,  // 0.680-0.690
0267     23.63250,  // 0.690-0.700
0268     23.69250,  // 0.700-0.710
0269     23.75000,  // 0.710-0.720
0270     23.80500,  // 0.720-0.730
0271     23.86000,  // 0.730-0.740
0272     23.91250,  // 0.740-0.750
0273     23.96500,  // 0.750-0.760
0274     24.01500,  // 0.760-0.770
0275     24.06500,  // 0.770-0.780
0276     24.11250,  // 0.780-0.790
0277     24.16000,  // 0.790-0.800
0278     24.20500,  // 0.800-0.810
0279     24.25000,  // 0.810-0.820
0280     24.29500,  // 0.820-0.830
0281     24.33750,  // 0.830-0.840
0282     24.38000,  // 0.840-0.850
0283     24.42250,  // 0.850-0.860
0284     24.46500,  // 0.860-0.870
0285     24.50500,  // 0.870-0.880
0286     24.54500,  // 0.880-0.890
0287     24.58500,  // 0.890-0.900
0288     24.62500,  // 0.900-0.910
0289     24.66500,  // 0.910-0.920
0290     24.70250,  // 0.920-0.930
0291     24.74000,  // 0.930-0.940
0292     24.77750,  // 0.940-0.950
0293     24.81500,  // 0.950-0.960
0294     24.85250,  // 0.960-0.970
0295     24.89000,  // 0.970-0.980
0296     24.92750,  // 0.980-0.990
0297     24.96250,  // 0.990-1.000
0298 };
0299
0300 float timeshift_ns_hf(float wpksamp) {
0301   float flx = (num_bins_hf * (wpksamp - wpksamp0_hf) / scale_hf);
0302   int index = (int)flx;
0303   float yval;
0304
0305   if (index < 0)
0306     return actual_ns_hf[0];
0307   else if (index >= num_bins_hf - 1)
0308     return actual_ns_hf[num_bins_hf - 1];
0309
0310   // else interpolate:
0311   float y1 = actual_ns_hf[index];
0312   float y2 = actual_ns_hf[index + 1];
0313
0314   // float delta_x  = 1/(float)num_bins_hf;
0315   // yval = y1 + (y2-y1)*(flx-(float)index)/delta_x;
0316
0317   yval = y1 + (y2 - y1) * (flx - (float)index);
0318   return yval;
0319 }