Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-08-05 11:15:00

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