File indexing completed on 2024-10-19 04:58:39
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
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
0018
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
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
0071
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 }
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
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
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
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
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
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
0182
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
0203
0204 if (AvgTSPos == 0) {
0205 time = -99;
0206 } else {
0207 time = (AvgTSPos * 25.0);
0208 }
0209 }
0210
0211
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
0218 ta = tool[CurrentTS];
0219 ta *= calibs.respcorrgain(capid);
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
0230 ta = std::max(0.0, tool[ts]);
0231 ta *= calibs.respcorrgain(capid);
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
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 }