HcalRecAlgos/src/ZdcSimpleRecAlgo.cc

0001 #include "RecoLocalCalo/HcalRecAlgos/interface/ZdcSimpleRecAlgo.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include "CalibCalorimetry/HcalAlgos/interface/HcalTimeSlew.h"
0004 #include <algorithm>  // for "max"
0005 #include <iostream>
0006 #include <cmath>
0007
0008 constexpr double MaximumFractionalError = 0.0005;  // 0.05% error allowed from this source
0009
0010 ZdcSimpleRecAlgo::ZdcSimpleRecAlgo(
0011     bool correctForTimeslew, bool correctForPulse, float phaseNS, int recoMethod, int lowGainOffset, double lowGainFrac)
0012     : recoMethod_(recoMethod),
0013       correctForTimeslew_(correctForTimeslew),
0014       correctForPulse_(correctForPulse),
0015       phaseNS_(phaseNS),
0016       lowGainOffset_(lowGainOffset),
0017       lowGainFrac_(lowGainFrac) {}
0018
0019 ZdcSimpleRecAlgo::ZdcSimpleRecAlgo(int recoMethod) : recoMethod_(recoMethod), correctForTimeslew_(false) {}
0020
0021 void ZdcSimpleRecAlgo::initPulseCorr(int toadd, const HcalTimeSlew* hcalTimeSlew_delay) {
0022   if (correctForPulse_) {
0023     pulseCorr_ = std::make_unique<HcalPulseContainmentCorrection>(
0024         toadd, phaseNS_, false, MaximumFractionalError, hcalTimeSlew_delay);
0025   }
0026 }
0027 //static float timeshift_ns_zdc(float wpksamp);
0028
0029 namespace ZdcSimpleRecAlgoImpl {
0030   template <class Digi, class RecHit>
0031   inline RecHit reco1(const Digi& digi,
0032                       const HcalCoder& coder,
0033                       const HcalCalibrations& calibs,
0034                       const std::vector<unsigned int>& myNoiseTS,
0035                       const std::vector<unsigned int>& mySignalTS,
0036                       int lowGainOffset,
0037                       double lowGainFrac,
0038                       bool slewCorrect,
0039                       const HcalPulseContainmentCorrection* corr,
0040                       HcalTimeSlew::BiasSetting slewFlavor) {
0041     CaloSamples tool;
0042     coder.adc2fC(digi, tool);
0043     int ifirst = mySignalTS[0];
0044     int n = mySignalTS.size();
0045     double ampl = 0;
0046     int maxI = -1;
0047     double maxA = -1e10;
0048     double ta = 0;
0049     double fc_ampl = 0;
0050     double lowGEnergy = 0;
0051     double TempLGAmp = 0;
0052     // TS increment for regular energy to lowGainEnergy
0053     // Signal in higher TS (effective "low Gain") has a fraction of the whole signal
0054     // This constant for fC --> GeV is dervied from 2010 PbPb analysis of single neutrons
0055     // assumed similar fraction for EM and HAD sections
0056     // this variable converts from current assumed TestBeam values for fC--> GeV
0057     // to the lowGain TS region fraction value (based on 1N Had, assume EM same response)
0058     // regular energy
0059     for (int i = ifirst; i < tool.size() && i < n + ifirst; i++) {
0060       int capid = digi[i].capid();
0061       ta = (tool[i] - calibs.pedestal(capid));  // pedestal subtraction
0062       fc_ampl += ta;
0063       ta *= calibs.respcorrgain(capid);  // fC --> GeV
0064       ampl += ta;
0065       if (ta > maxA) {
0066         maxA = ta;
0067         maxI = i;
0068       }
0069     }
0070     // calculate low Gain Energy (in 2010 PbPb, signal TS 4,5,6, lowGain TS: 6,7,8)
0071     int topLowGain = 10;
0072     if ((n + ifirst + lowGainOffset) <= 10) {
0073       topLowGain = n + ifirst + lowGainOffset;
0074     } else {
0075       topLowGain = 10;
0076     }
0077     for (int iLG = (ifirst + lowGainOffset); iLG < tool.size() && iLG < topLowGain; iLG++) {
0078       int capid = digi[iLG].capid();
0079       TempLGAmp = (tool[iLG] - calibs.pedestal(capid));  // pedestal subtraction
0080       TempLGAmp *= calibs.respcorrgain(capid);           // fC --> GeV
0081       TempLGAmp *= lowGainFrac;                          // TS (signalRegion) --> TS (lowGainRegion)
0082       lowGEnergy += TempLGAmp;
0083     }
0084     double time = -9999;
0085     // Time based on regular energy (lowGainEnergy signal assumed to happen at same Time)
0086     ////Cannot calculate time value with max ADC sample at first or last position in window....
0087     if (maxI == 0 || maxI == (tool.size() - 1)) {
0088       LogDebug("HCAL Pulse") << "ZdcSimpleRecAlgo::reco1 :"
0089                              << " Invalid max amplitude position, "
0090                              << " max Amplitude: " << maxI << " first: " << ifirst << " last: " << (tool.size() - 1)
0091                              << std::endl;
0092     } else {
0093       int capid = digi[maxI - 1].capid();
0094       double Energy0 = ((tool[maxI - 1]) * calibs.respcorrgain(capid));
0095       // if any of the energies used in the weight are negative, make them 0 instead
0096       // these are actually QIE values, not energy
0097       if (Energy0 < 0) {
0098         Energy0 = 0.;
0099       }
0100       capid = digi[maxI].capid();
0101       double Energy1 = ((tool[maxI]) * calibs.respcorrgain(capid));
0102       if (Energy1 < 0) {
0103         Energy1 = 0.;
0104       }
0105       capid = digi[maxI + 1].capid();
0106       double Energy2 = ((tool[maxI + 1]) * calibs.respcorrgain(capid));
0107       if (Energy2 < 0) {
0108         Energy2 = 0.;
0109       }
0110       //
0111       double TSWeightEnergy = ((maxI - 1) * Energy0 + maxI * Energy1 + (maxI + 1) * Energy2);
0112       double EnergySum = Energy0 + Energy1 + Energy2;
0113       double AvgTSPos = 0.;
0114       if (EnergySum != 0)
0115         AvgTSPos = TSWeightEnergy / EnergySum;
0116       // If time is zero, set it to the "nonsensical" -99
0117       // Time should be between 75ns and 175ns (Timeslices 3-7)
0118       if (AvgTSPos == 0) {
0119         time = -99;
0120       } else {
0121         time = (AvgTSPos * 25.0);
0122       }
0123       if (corr != nullptr) {
0124         // Apply phase-based amplitude correction:
0125         ampl *= corr->getCorrection(fc_ampl);
0126       }
0127     }
0128     return RecHit(digi.id(), ampl, time, lowGEnergy);
0129   }
0130 }  // namespace ZdcSimpleRecAlgoImpl
0131
0132 namespace ZdcSimpleRecAlgoImpl {
0133   template <class Digi, class RecHit>
0134   inline RecHit reco2(const Digi& digi,
0135                       const HcalCoder& coder,
0136                       const HcalCalibrations& calibs,
0137                       const std::vector<unsigned int>& myNoiseTS,
0138                       const std::vector<unsigned int>& mySignalTS,
0139                       int lowGainOffset,
0140                       double lowGainFrac,
0141                       bool slewCorrect,
0142                       const HcalPulseContainmentCorrection* corr,
0143                       HcalTimeSlew::BiasSetting slewFlavor) {
0144     CaloSamples tool;
0145     coder.adc2fC(digi, tool);
0146     // Reads noiseTS and signalTS from database
0147     int ifirst = mySignalTS[0];
0148     //    int n = mySignalTS.size();
0149     double ampl = 0;
0150     int maxI = -1;
0151     double maxA = -1e10;
0152     double ta = 0;
0153     double fc_ampl = 0;
0154     double lowGEnergy = 0;
0155     double TempLGAmp = 0;
0156     //  TS increment for regular energy to lowGainEnergy
0157     // Signal in higher TS (effective "low Gain") has a fraction of the whole signal
0158     // This constant for fC --> GeV is dervied from 2010 PbPb analysis of single neutrons
0159     // assumed similar fraction for EM and HAD sections
0160     // this variable converts from current assumed TestBeam values for fC--> GeV
0161     // to the lowGain TS region fraction value (based on 1N Had, assume EM same response)
0162     double Allnoise = 0;
0163     int noiseslices = 0;
0164     int CurrentTS = 0;
0165     double noise = 0;
0166     // regular energy (both use same noise)
0167     for (unsigned int iv = 0; iv < myNoiseTS.size(); ++iv) {
0168       CurrentTS = myNoiseTS[iv];
0169       if (CurrentTS >= digi.size())
0170         continue;
0171       Allnoise += tool[CurrentTS];
0172       noiseslices++;
0173     }
0174     if (noiseslices != 0) {
0175       noise = (Allnoise) / double(noiseslices);
0176     } else {
0177       noise = 0;
0178     }
0179     for (unsigned int ivs = 0; ivs < mySignalTS.size(); ++ivs) {
0180       CurrentTS = mySignalTS[ivs];
0181       if (CurrentTS >= digi.size())
0182         continue;
0183       int capid = digi[CurrentTS].capid();
0184       //       if(noise<0){
0185       //       // flag hit as having negative noise, and don't subtract anything, because
0186       //       // it will falsely increase the energy
0187       //          noisefactor=0.;
0188       //       }
0189       ta = tool[CurrentTS] - noise;
0190       fc_ampl += ta;
0191       ta *= calibs.respcorrgain(capid);  // fC --> GeV
0192       ampl += ta;
0193       if (ta > maxA) {
0194         maxA = ta;
0195         maxI = CurrentTS;
0196       }
0197     }
0198     // calculate low Gain Energy (in 2010 PbPb, signal TS 4,5,6, lowGain TS: 6,7,8)
0199     for (unsigned int iLGvs = 0; iLGvs < mySignalTS.size(); ++iLGvs) {
0200       CurrentTS = mySignalTS[iLGvs] + lowGainOffset;
0201       if (CurrentTS >= digi.size())
0202         continue;
0203       int capid = digi[CurrentTS].capid();
0204       TempLGAmp = tool[CurrentTS] - noise;
0205       TempLGAmp *= calibs.respcorrgain(capid);  // fC --> GeV
0206       TempLGAmp *= lowGainFrac;                 // TS (signalRegion) --> TS (lowGainRegion)
0207       lowGEnergy += TempLGAmp;
0208     }
0209     //    if(ta<0){
0210     //      // flag hits that have negative energy
0211     //    }
0212
0213     double time = -9999;
0214     // Time based on regular energy (lowGainEnergy signal assumed to happen at same Time)
0215     ////Cannot calculate time value with max ADC sample at first or last position in window....
0216     if (maxI == 0 || maxI == (tool.size() - 1)) {
0217       LogDebug("HCAL Pulse") << "ZdcSimpleRecAlgo::reco2 :"
0218                              << " Invalid max amplitude position, "
0219                              << " max Amplitude: " << maxI << " first: " << ifirst << " last: " << (tool.size() - 1)
0220                              << std::endl;
0221     } else {
0222       int capid = digi[maxI - 1].capid();
0223       double Energy0 = ((tool[maxI - 1]) * calibs.respcorrgain(capid));
0224       // if any of the energies used in the weight are negative, make them 0 instead
0225       // these are actually QIE values, not energy
0226       if (Energy0 < 0) {
0227         Energy0 = 0.;
0228       }
0229       capid = digi[maxI].capid();
0230       double Energy1 = ((tool[maxI]) * calibs.respcorrgain(capid));
0231       if (Energy1 < 0) {
0232         Energy1 = 0.;
0233       }
0234       capid = digi[maxI + 1].capid();
0235       double Energy2 = ((tool[maxI + 1]) * calibs.respcorrgain(capid));
0236       if (Energy2 < 0) {
0237         Energy2 = 0.;
0238       }
0239       //
0240       double TSWeightEnergy = ((maxI - 1) * Energy0 + maxI * Energy1 + (maxI + 1) * Energy2);
0241       double EnergySum = Energy0 + Energy1 + Energy2;
0242       double AvgTSPos = 0.;
0243       if (EnergySum != 0)
0244         AvgTSPos = TSWeightEnergy / EnergySum;
0245       // If time is zero, set it to the "nonsensical" -99
0246       // Time should be between 75ns and 175ns (Timeslices 3-7)
0247       if (AvgTSPos == 0) {
0248         time = -99;
0249       } else {
0250         time = (AvgTSPos * 25.0);
0251       }
0252       if (corr != nullptr) {
0253         // Apply phase-based amplitude correction:
0254         ampl *= corr->getCorrection(fc_ampl);
0255       }
0256     }
0257     return RecHit(digi.id(), ampl, time, lowGEnergy);
0258   }
0259 }  // namespace ZdcSimpleRecAlgoImpl
0260
0261 ZDCRecHit ZdcSimpleRecAlgo::reconstruct(const ZDCDataFrame& digi,
0262                                         const std::vector<unsigned int>& myNoiseTS,
0263                                         const std::vector<unsigned int>& mySignalTS,
0264                                         const HcalCoder& coder,
0265                                         const HcalCalibrations& calibs) const {
0266   if (recoMethod_ == 1)
0267     return ZdcSimpleRecAlgoImpl::reco1<ZDCDataFrame, ZDCRecHit>(
0268         digi, coder, calibs, myNoiseTS, mySignalTS, lowGainOffset_, lowGainFrac_, false, nullptr, HcalTimeSlew::Fast);
0269   if (recoMethod_ == 2)
0270     return ZdcSimpleRecAlgoImpl::reco2<ZDCDataFrame, ZDCRecHit>(
0271         digi, coder, calibs, myNoiseTS, mySignalTS, lowGainOffset_, lowGainFrac_, false, nullptr, HcalTimeSlew::Fast);
0272
0273   edm::LogError("ZDCSimpleRecAlgoImpl::reconstruct, recoMethod was not declared");
0274   throw cms::Exception("ZDCSimpleRecoAlgos::exiting process");
0275 }