HcalRecAlgos/src/ZdcSimpleRecAlgo_Run3.cc

0001 #include "RecoLocalCalo/HcalRecAlgos/interface/ZdcSimpleRecAlgo_Run3.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include <algorithm>  // for "max"
0004 #include <iostream>
0005 #include <cmath>
0006
0007 #include <Eigen/Dense>  // for TemplateFit Method
0008
0009 // #include "DataFormats/METReco/interface/HcalCaloFlagLabels.h"
0010
0011 ZdcSimpleRecAlgo_Run3::ZdcSimpleRecAlgo_Run3(int recoMethod) : recoMethod_(recoMethod) {}
0012
0013 void ZdcSimpleRecAlgo_Run3::initCorrectionMethod(const int method, const int ZdcSection) {
0014   correctionMethod_[ZdcSection] = method;
0015 };
0016
0017 // Template fit is a linear combination of timeslices in a digi assuming there is a potential charge peak at each Bx
0018 // and the charges of the Ts follwing a peak are consistent with the chargeRatios
0019 void ZdcSimpleRecAlgo_Run3::initTemplateFit(const std::vector<unsigned int>& bxTs,
0020                                             const std::vector<double>& chargeRatios,
0021                                             const int nTs,
0022                                             const int ZdcSection) {
0023   int nRatios = chargeRatios.size();
0024   int nBx = bxTs.size();
0025   int nCols = nBx + 1;
0026   double val = 0;
0027   int index = 0;
0028   int timeslice = 0;
0029   nTs_ = nTs;
0030   Eigen::MatrixXf a(nTs, nCols);
0031   for (int j = 0; j < nBx; j++) {
0032     timeslice = bxTs.at(j);
0033     for (int i = 0; i < nTs; i++) {
0034       val = 0;
0035       index = i - timeslice;
0036       if (index >= 0 && index < nRatios)
0037         val = chargeRatios.at(index);
0038       a(i, j) = val;
0039     }
0040   }
0041   for (int i = 0; i < nTs; i++)
0042     a(i, nBx) = 1;
0043   Eigen::MatrixXf b = a.transpose() * a;
0044   if (std::fabs(b.determinant()) < 1E-8) {
0045     templateFitValid_[ZdcSection] = false;
0046     return;
0047   }
0048   templateFitValid_[ZdcSection] = true;
0049   Eigen::MatrixXf tfMatrix;
0050   tfMatrix = b.inverse() * a.transpose();
0051   for (int i = 0; i < nTs; i++)
0052     templateFitValues_[ZdcSection].push_back(tfMatrix.coeff(1, i));
0053 }
0054
0055 void ZdcSimpleRecAlgo_Run3::initRatioSubtraction(const float ratio, const float frac, const int ZdcSection) {
0056   ootpuRatio_[ZdcSection] = ratio;
0057   ootpuFrac_[ZdcSection] = frac;
0058 }
0059
0060 // helper functions for pedestal subtraction and noise calculations
0061 namespace zdchelper {
0062
0063   inline double subPedestal(const float charge, const float ped, const float width) {
0064     if (charge - ped > width)
0065       return (charge - ped);
0066     else
0067       return (0);
0068   }
0069
0070   // gets the noise with a check if the ratio of energy0 / energy1 > ootpuRatio
0071   // energy1 is noiseTS , energy0 is noiseTs-1
0072   inline double getNoiseOOTPURatio(const float energy0,
0073                                    const float energy1,
0074                                    const float ootpuRatio,
0075                                    const float ootpuFrac) {
0076     if (energy0 >= ootpuRatio * energy1 || ootpuRatio < 0)
0077       return (ootpuFrac * energy1);
0078     else
0079       return (energy1);
0080   }
0081
0082 }  // namespace zdchelper
0083
0084 ZDCRecHit ZdcSimpleRecAlgo_Run3::reco0(const QIE10DataFrame& digi,
0085                                        const HcalCoder& coder,
0086                                        const HcalCalibrations& calibs,
0087                                        const HcalPedestal& effPeds,
0088                                        const std::vector<unsigned int>& myNoiseTS,
0089                                        const std::vector<unsigned int>& mySignalTS) const {
0090   CaloSamples tool;
0091   coder.adc2fC(digi, tool);
0092   // Reads noiseTS and signalTS from database
0093   int ifirst = mySignalTS[0];
0094   double ampl = 0;
0095   int maxI = -1;
0096   double maxA = -1e10;
0097   double ta = 0;
0098   double energySOIp1 = 0;
0099   double ratioSOIp1 = -1.0;
0100   double chargeWeightedTime = -99.0;
0101
0102   double Allnoise = 0;
0103   int noiseslices = 0;
0104   int CurrentTS = 0;
0105   double noise = 0;
0106   int digi_size = digi.samples();
0107   HcalZDCDetId cell = digi.id();
0108   int zdcsection = cell.section();
0109
0110   // determining noise
0111   for (unsigned int iv = 0; iv < myNoiseTS.size(); ++iv) {
0112     CurrentTS = myNoiseTS[iv];
0113     int capid = digi[CurrentTS].capid();
0114     float ped = effPeds.getValue(capid);
0115     float width = effPeds.getWidth(capid);
0116     float gain = calibs.respcorrgain(capid);
0117     if (CurrentTS >= digi_size)
0118       continue;
0119     float energy1 = zdchelper::subPedestal(tool[CurrentTS], ped, width) * gain;
0120     float energy0 = 0;
0121     if (CurrentTS > 0) {
0122       capid = digi[CurrentTS - 1].capid();
0123       ped = effPeds.getValue(capid);
0124       width = effPeds.getWidth(capid);
0125       gain = calibs.respcorrgain(capid);
0126       energy0 = zdchelper::subPedestal(tool[CurrentTS - 1], ped, width) * gain;
0127     }
0128     Allnoise += zdchelper::getNoiseOOTPURatio(energy0, energy1, ootpuRatio_.at(zdcsection), ootpuFrac_.at(zdcsection));
0129     noiseslices++;
0130   }
0131   if (noiseslices != 0) {
0132     noise = (Allnoise) / double(noiseslices);
0133   }
0134
0135   // determining signal energy and max Ts
0136   for (unsigned int ivs = 0; ivs < mySignalTS.size(); ++ivs) {
0137     CurrentTS = mySignalTS[ivs];
0138     if (CurrentTS >= digi_size)
0139       continue;
0140     float energy1 = -1;
0141     int capid = digi[CurrentTS].capid();
0142     // float ped = calibs.pedestal(capid);
0143     float ped = effPeds.getValue(capid);
0144     float width = effPeds.getWidth(capid);
0145     float gain = calibs.respcorrgain(capid);
0146     float energy0 = std::max(0.0, zdchelper::subPedestal(tool[CurrentTS], ped, width)) * gain;
0147     if (CurrentTS < digi_size - 1) {
0148       capid = digi[CurrentTS].capid();
0149       ped = effPeds.getValue(capid);
0150       width = effPeds.getWidth(capid);
0151       gain = calibs.respcorrgain(capid);
0152       energy1 = std::max(0.0, zdchelper::subPedestal(tool[CurrentTS + 1], ped, width)) * gain;
0153     }
0154     ta = energy0 - noise;
0155     if (ta > 0)
0156       ampl += ta;
0157
0158     if (ta > maxA) {
0159       ratioSOIp1 = (energy0 > 0 && energy1 > 0) ? energy0 / energy1 : -1.0;
0160       maxA = ta;
0161       maxI = CurrentTS;
0162     }
0163   }
0164
0165   // determine energy if using Template Fit method
0166   if (correctionMethod_.at(zdcsection) == 1 && templateFitValid_.at(zdcsection)) {
0167     double energy = 0;
0168     for (int iv = 0; iv < nTs_; iv++) {
0169       int capid = digi[iv].capid();
0170       float ped = effPeds.getValue(capid);
0171       float width = effPeds.getWidth(capid);
0172       float gain = calibs.respcorrgain(capid);
0173       if (iv >= digi_size)
0174         continue;
0175       energy += zdchelper::subPedestal(tool[iv], ped, width) * (templateFitValues_.at(zdcsection)).at(iv) * gain;
0176     }
0177     ampl = std::max(0.0, energy);
0178   }
0179
0180   double time = -9999;
0181   // Time based on regular energy
0182   ////Cannot calculate time value with max ADC sample at first or last position in window....
0183   if (maxI == 0 || maxI == (tool.size() - 1)) {
0184     LogDebug("HCAL Pulse") << "ZdcSimpleRecAlgo::reco2 :"
0185                            << " Invalid max amplitude position, "
0186                            << " max Amplitude: " << maxI << " first: " << ifirst << " last: " << (tool.size() - 1)
0187                            << std::endl;
0188   } else {
0189     int capid = digi[maxI - 1].capid();
0190     double Energy0 = std::max(0.0, ((tool[maxI - 1]) * calibs.respcorrgain(capid)));
0191
0192     capid = digi[maxI].capid();
0193     double Energy1 = std::max(0.0, ((tool[maxI]) * calibs.respcorrgain(capid)));
0194     capid = digi[maxI + 1].capid();
0195     double Energy2 = std::max(0.0, ((tool[maxI + 1]) * calibs.respcorrgain(capid)));
0196
0197     double TSWeightEnergy = ((maxI - 1) * Energy0 + maxI * Energy1 + (maxI + 1) * Energy2);
0198     double EnergySum = Energy0 + Energy1 + Energy2;
0199     double AvgTSPos = 0.;
0200     if (EnergySum != 0)
0201       AvgTSPos = TSWeightEnergy / EnergySum;
0202     // If time is zero, set it to the "nonsensical" -99
0203     // Time should be between 75ns and 175ns (Timeslices 3-7)
0204     if (AvgTSPos == 0) {
0205       time = -99;
0206     } else {
0207       time = (AvgTSPos * 25.0);
0208     }
0209   }
0210
0211   // find energy for signal TS + 1
0212   for (unsigned int ivs = 0; ivs < mySignalTS.size(); ++ivs) {
0213     CurrentTS = mySignalTS[ivs] + 1;
0214     if (CurrentTS >= digi_size)
0215       continue;
0216     int capid = digi[CurrentTS].capid();
0217     // ta = tool[CurrentTS] - noise;
0218     ta = tool[CurrentTS];
0219     ta *= calibs.respcorrgain(capid);  // fC --> GeV
0220     if (ta > 0)
0221       energySOIp1 += ta;
0222   }
0223
0224   double tmp_energy = 0;
0225   double tmp_TSWeightedEnergy = 0;
0226   for (int ts = 0; ts < digi_size; ++ts) {
0227     int capid = digi[ts].capid();
0228
0229     // max sure there are no negative values in time calculation
0230     ta = std::max(0.0, tool[ts]);
0231     ta *= calibs.respcorrgain(capid);  // fC --> GeV
0232     if (ta > 0) {
0233       tmp_energy += ta;
0234       tmp_TSWeightedEnergy += (ts)*ta;
0235     }
0236   }
0237
0238   if (tmp_energy > 0)
0239     chargeWeightedTime = (tmp_TSWeightedEnergy / tmp_energy) * 25.0;
0240   auto rh = ZDCRecHit(digi.id(), ampl, time, -99);
0241   rh.setEnergySOIp1(energySOIp1);
0242
0243   if (maxI >= 0 && maxI < tool.size()) {
0244     float tmp_tdctime = 0;
0245     int le_tdc = digi[maxI].le_tdc();
0246     // TDC error codes will be 60=-1, 61 = -2, 62 = -3, 63 = -4
0247     if (le_tdc >= 60)
0248       tmp_tdctime = -1 * (le_tdc - 59);
0249     else
0250       tmp_tdctime = maxI * 25. + (le_tdc / 2.0);
0251     rh.setTDCtime(tmp_tdctime);
0252   }
0253
0254   rh.setChargeWeightedTime(chargeWeightedTime);
0255   rh.setRatioSOIp1(ratioSOIp1);
0256   return rh;
0257 }
0258
0259 ZDCRecHit ZdcSimpleRecAlgo_Run3::reconstruct(const QIE10DataFrame& digi,
0260                                              const std::vector<unsigned int>& myNoiseTS,
0261                                              const std::vector<unsigned int>& mySignalTS,
0262                                              const HcalCoder& coder,
0263                                              const HcalCalibrations& calibs,
0264                                              const HcalPedestal& effPeds) const {
0265   return ZdcSimpleRecAlgo_Run3::reco0(digi, coder, calibs, effPeds, myNoiseTS, mySignalTS);
0266
0267   edm::LogError("ZDCSimpleRecAlgoImpl::reconstruct, recoMethod was not declared");
0268   throw cms::Exception("ZDCSimpleRecoAlgos::exiting process");
0269 }