diff --git a/Configuration/StandardSequences/python/Reconstruction_cff.py b/Configuration/StandardSequences/python/Reconstruction_cff.py index 021f5f090cbe1..b30aeb4617ca9 100644 --- a/Configuration/StandardSequences/python/Reconstruction_cff.py +++ b/Configuration/StandardSequences/python/Reconstruction_cff.py @@ -304,6 +304,7 @@ modulesToRemove.append(horeco) modulesToRemove.append(hcalnoise) modulesToRemove.append(zdcreco) +modulesToRemove.append(zdcrecoRun3) modulesToRemove.append(castorreco) ##it's OK according to Ronny modulesToRemove.append(CSCHaloData)#needs digis reconstruction_fromRECO = reconstruction.copyAndExclude(modulesToRemove+noTrackingAndDependent) diff --git a/DataFormats/HcalRecHit/interface/ZDCRecHit.h b/DataFormats/HcalRecHit/interface/ZDCRecHit.h index 7b63892b2f644..7b6793c19b2aa 100644 --- a/DataFormats/HcalRecHit/interface/ZDCRecHit.h +++ b/DataFormats/HcalRecHit/interface/ZDCRecHit.h @@ -19,8 +19,25 @@ class ZDCRecHit : public CaloRecHit { // follow EcalRecHit method of adding variable flagBits_ to CaloRecHit float lowGainEnergy() const { return lowGainEnergy_; }; + constexpr inline void setEnergySOIp1(const float en) { energySOIp1_ = en; }; + constexpr inline float energySOIp1() const { return energySOIp1_; }; // energy of Slice of Interest plus 1 + constexpr inline void setRatioSOIp1(const float ratio) { ratioSOIp1_ = ratio; }; + constexpr inline float ratioSOIp1() const { + return ratioSOIp1_; + }; // ratio of Energy of (Slice of Interest)/ (Slice of Interest plus 1) + constexpr inline void setTDCtime(const float time) { TDCtime_ = time; }; + constexpr inline float TDCtime() const { return TDCtime_; }; + constexpr inline void setChargeWeightedTime(const float time) { + chargeWeightedTime_ = time; + }; // time of activity determined by charged weighted average + constexpr inline float chargeWeightedTime() const { return chargeWeightedTime_; }; + private: float lowGainEnergy_; + float energySOIp1_; + float ratioSOIp1_; + float TDCtime_; + float chargeWeightedTime_; }; std::ostream& operator<<(std::ostream& s, const ZDCRecHit& hit); diff --git a/DataFormats/HcalRecHit/src/ZDCRecHit.cc b/DataFormats/HcalRecHit/src/ZDCRecHit.cc index 9d0b69dc42ed7..1c20349a073f4 100644 --- a/DataFormats/HcalRecHit/src/ZDCRecHit.cc +++ b/DataFormats/HcalRecHit/src/ZDCRecHit.cc @@ -1,9 +1,15 @@ #include "DataFormats/HcalRecHit/interface/ZDCRecHit.h" -ZDCRecHit::ZDCRecHit() : CaloRecHit(), lowGainEnergy_() {} +ZDCRecHit::ZDCRecHit() + : CaloRecHit(), lowGainEnergy_(), energySOIp1_(-99), ratioSOIp1_(-99), TDCtime_(-99), chargeWeightedTime_(-99) {} ZDCRecHit::ZDCRecHit(const HcalZDCDetId& id, float energy, float time, float lowGainEnergy) - : CaloRecHit(id, energy, time), lowGainEnergy_(lowGainEnergy) {} + : CaloRecHit(id, energy, time), + lowGainEnergy_(lowGainEnergy), + energySOIp1_(-99), + ratioSOIp1_(-99), + TDCtime_(-99), + chargeWeightedTime_(-99) {} std::ostream& operator<<(std::ostream& s, const ZDCRecHit& hit) { return s << hit.id() << ": " << hit.energy() << " GeV, " << hit.time() << " ns"; diff --git a/DataFormats/HcalRecHit/src/classes_def.xml b/DataFormats/HcalRecHit/src/classes_def.xml index 4cfb4c2be84a6..b5520847e8d86 100644 --- a/DataFormats/HcalRecHit/src/classes_def.xml +++ b/DataFormats/HcalRecHit/src/classes_def.xml @@ -35,7 +35,8 @@ - + + diff --git a/RecoHI/HiCentralityAlgos/python/HiCentrality_cfi.py b/RecoHI/HiCentralityAlgos/python/HiCentrality_cfi.py index 51503a73a0465..a1afbb4cdbf52 100644 --- a/RecoHI/HiCentralityAlgos/python/HiCentrality_cfi.py +++ b/RecoHI/HiCentralityAlgos/python/HiCentrality_cfi.py @@ -46,10 +46,14 @@ srcTracks = "generalTracks", srcVertex = "offlinePrimaryVertices" ) + +from Configuration.Eras.Modifier_run3_common_cff import run3_common +run3_common.toModify(hiCentrality, srcZDChits = "zdcrecoRun3",lowGainZDC = False) + from Configuration.ProcessModifiers.phase2_pp_on_AA_cff import phase2_pp_on_AA phase2_pp_on_AA.toModify(hiCentrality, isPhase2 = True, producePixelTracks = False, srcTracks = "generalTracks", srcVertex = "offlinePrimaryVertices" -) +) \ No newline at end of file diff --git a/RecoHI/HiCentralityAlgos/python/pACentrality_cfi.py b/RecoHI/HiCentralityAlgos/python/pACentrality_cfi.py index 15badf74fe346..ec2f847f5cabb 100644 --- a/RecoHI/HiCentralityAlgos/python/pACentrality_cfi.py +++ b/RecoHI/HiCentralityAlgos/python/pACentrality_cfi.py @@ -36,4 +36,6 @@ ) +from Configuration.Eras.Modifier_run3_common_cff import run3_common +run3_common.toModify(pACentrality, srcZDChits = "zdcrecoRun3",lowGainZDC =False) diff --git a/RecoHI/HiEvtPlaneAlgos/python/RecoHiEvtPlane_EventContent_cff.py b/RecoHI/HiEvtPlaneAlgos/python/RecoHiEvtPlane_EventContent_cff.py index 33e7b661b7c8f..d89eb94287faf 100644 --- a/RecoHI/HiEvtPlaneAlgos/python/RecoHiEvtPlane_EventContent_cff.py +++ b/RecoHI/HiEvtPlaneAlgos/python/RecoHiEvtPlane_EventContent_cff.py @@ -5,6 +5,7 @@ outputCommands = cms.untracked.vstring( 'keep recoEvtPlanes_hiEvtPlane_*_*', 'keep ZDCRecHitsSorted_zdcreco_*_*', + 'keep ZDCRecHitsSorted_zdcrecoRun3_*_*', 'keep ZDCDataFramesSorted_hcalDigis_*_*', 'keep HFRecHitsSorted_hfreco_*_*') ) diff --git a/RecoLocalCalo/Configuration/python/RecoLocalCalo_Cosmics_cff.py b/RecoLocalCalo/Configuration/python/RecoLocalCalo_Cosmics_cff.py index 25b8fed436ae2..ac5ac4b452f07 100644 --- a/RecoLocalCalo/Configuration/python/RecoLocalCalo_Cosmics_cff.py +++ b/RecoLocalCalo/Configuration/python/RecoLocalCalo_Cosmics_cff.py @@ -107,3 +107,11 @@ from RecoLocalCalo.Configuration.hcalLocalRecoNZS_cff import * calolocalrecoTaskCosmicsNZS = cms.Task(ecalLocalRecoTaskCosmics,hcalLocalRecoTask,hcalLocalRecoTaskNZS) calolocalrecoCosmicsNZS = cms.Sequence(calolocalrecoTaskCosmicsNZS) + +#--- for Run 3 and later +_run3_hcalLocalRecoTask = _phase1_hcalLocalRecoTask.copy() +from RecoLocalCalo.HcalRecProducers.zdcrecoRun3_cfi import zdcrecoRun3 +_run3_hcalLocalRecoTask.remove(zdcreco) +_run3_hcalLocalRecoTask.add(zdcrecoRun3) +from Configuration.Eras.Modifier_run3_common_cff import run3_common +run3_common.toReplaceWith(hcalLocalRecoTask, _run3_hcalLocalRecoTask) diff --git a/RecoLocalCalo/Configuration/python/RecoLocalCalo_EventContentCosmics_cff.py b/RecoLocalCalo/Configuration/python/RecoLocalCalo_EventContentCosmics_cff.py index f89863040ee6b..73dcdc1210797 100644 --- a/RecoLocalCalo/Configuration/python/RecoLocalCalo_EventContentCosmics_cff.py +++ b/RecoLocalCalo/Configuration/python/RecoLocalCalo_EventContentCosmics_cff.py @@ -24,6 +24,7 @@ 'keep ZDCDataFramesSorted_castorDigis_*_*', 'keep ZDCDataFramesSorted_simHcalUnsuppressedDigis_*_*', 'keep ZDCRecHitsSorted_zdcreco_*_*', + 'keep ZDCRecHitsSorted_zdcrecoRun3_*_*', 'keep HcalUnpackerReport_castorDigis_*_*', 'keep HcalUnpackerReport_hcalDigis_*_*') ) diff --git a/RecoLocalCalo/Configuration/python/RecoLocalCalo_EventContent_cff.py b/RecoLocalCalo/Configuration/python/RecoLocalCalo_EventContent_cff.py index 515b8374661da..056eab19a2fda 100644 --- a/RecoLocalCalo/Configuration/python/RecoLocalCalo_EventContent_cff.py +++ b/RecoLocalCalo/Configuration/python/RecoLocalCalo_EventContent_cff.py @@ -59,7 +59,8 @@ 'keep ZDCDataFramesSorted_hcalDigis_*_*', 'keep ZDCDataFramesSorted_castorDigis_*_*', 'keep QIE10DataFrameHcalDataFrameContainer_hcalDigis_ZDC_*', - 'keep ZDCRecHitsSorted_zdcreco_*_*') + 'keep ZDCRecHitsSorted_zdcreco_*_*', + 'keep ZDCRecHitsSorted_zdcrecoRun3_*_*') ) RecoLocalCaloRECO.outputCommands.extend(RecoLocalCaloAOD.outputCommands) RecoLocalCaloRECO.outputCommands.extend(ecalLocalRecoRECO.outputCommands) diff --git a/RecoLocalCalo/Configuration/python/hcalLocalReco_cff.py b/RecoLocalCalo/Configuration/python/hcalLocalReco_cff.py index feffe2e21b3fe..1b797e466f76e 100644 --- a/RecoLocalCalo/Configuration/python/hcalLocalReco_cff.py +++ b/RecoLocalCalo/Configuration/python/hcalLocalReco_cff.py @@ -60,8 +60,12 @@ #--- for Run 3 and later _run3_hcalLocalRecoTask = _phase1_hcalLocalRecoTask.copy() _run3_hcalLocalRecoTask.remove(hbheprereco) -from Configuration.Eras.Modifier_run3_HB_cff import run3_HB -run3_HB.toReplaceWith(hcalLocalRecoTask, _run3_hcalLocalRecoTask) + +from RecoLocalCalo.HcalRecProducers.zdcrecoRun3_cfi import zdcrecoRun3 +_run3_hcalLocalRecoTask.remove(zdcreco) +_run3_hcalLocalRecoTask.add(zdcrecoRun3) +from Configuration.Eras.Modifier_run3_common_cff import run3_common +run3_common.toReplaceWith(hcalLocalRecoTask, _run3_hcalLocalRecoTask) #--- for Run 3 on GPU from Configuration.ProcessModifiers.gpu_cff import gpu @@ -79,9 +83,10 @@ alpaka.toReplaceWith(hcalLocalRecoTask, _run3_hcalLocalRecoPortableTask) #--- HCAL-only workflow -hcalOnlyLocalRecoTask = hcalLocalRecoTask.copyAndExclude([zdcreco]) +hcalOnlyLocalRecoTask = hcalLocalRecoTask.copyAndExclude([zdcreco,zdcrecoRun3]) #--- HCAL-only workflow for Run 2 on GPU +from Configuration.Eras.Modifier_run3_HB_cff import run3_HB from RecoLocalCalo.HcalRecProducers.hcalCPURecHitsProducer_cfi import hcalCPURecHitsProducer as _hbheprerecoFromCUDA (gpu & ~run3_HB).toModify(hbheprereco, cuda = _hbheprerecoFromCUDA.clone( @@ -95,6 +100,6 @@ ) #--- for FastSim -_fastSim_hcalLocalRecoTask = hcalLocalRecoTask.copyAndExclude([zdcreco]) +_fastSim_hcalLocalRecoTask = hcalLocalRecoTask.copyAndExclude([zdcreco,zdcrecoRun3]) from Configuration.Eras.Modifier_fastSim_cff import fastSim fastSim.toReplaceWith( hcalLocalRecoTask, _fastSim_hcalLocalRecoTask ) diff --git a/RecoLocalCalo/HcalRecAlgos/interface/ZdcSimpleRecAlgo_Run3.h b/RecoLocalCalo/HcalRecAlgos/interface/ZdcSimpleRecAlgo_Run3.h new file mode 100644 index 0000000000000..39ccef9f9ea95 --- /dev/null +++ b/RecoLocalCalo/HcalRecAlgos/interface/ZdcSimpleRecAlgo_Run3.h @@ -0,0 +1,74 @@ +#ifndef ZDCSIMPLERECALGO_RUN3_H +#define ZDCSIMPLERECALGO_RUN3_H 1 + +#include "DataFormats/HcalDigi/interface/HBHEDataFrame.h" +#include "DataFormats/HcalDigi/interface/HFDataFrame.h" +#include "DataFormats/HcalDigi/interface/HODataFrame.h" +#include "DataFormats/HcalDigi/interface/ZDCDataFrame.h" +#include "DataFormats/HcalDigi/interface/HcalCalibDataFrame.h" +#include "DataFormats/HcalRecHit/interface/HBHERecHit.h" +#include "DataFormats/HcalRecHit/interface/HFRecHit.h" +#include "DataFormats/HcalRecHit/interface/HORecHit.h" +#include "DataFormats/HcalRecHit/interface/ZDCRecHit.h" +#include "DataFormats/HcalRecHit/interface/HcalCalibRecHit.h" +#include "CalibFormats/HcalObjects/interface/HcalCoder.h" +#include "CalibFormats/HcalObjects/interface/HcalCalibrations.h" +#include "CalibFormats/HcalObjects/interface/HcalCalibrationWidths.h" +#include "CondFormats/HcalObjects/interface/HcalPedestal.h" +#include "CalibCalorimetry/HcalAlgos/interface/HcalPulseContainmentCorrection.h" +#include + +#include "DataFormats/HcalDigi/interface/QIE10DataFrame.h" + +/** \class ZdcSimpleRecAlgo_Run3 + + This class reconstructs RecHits from Digis for ZDC by addition + of selected time samples, pedestal subtraction, and gain application. The + time of the hit is reconstructed using a weighted peak bin calculation + supplemented by precise time lookup table. A consumer of this class also + has the option of correcting the reconstructed time for energy-dependent + time slew associated with the QIE. + + A sencon method based on a based on a event by event substraction is also + implelented. signal = (S4 + S5 - 2*(S1+S2+S3 + S7+S8+S9+S10))*(ft-Gev constant) + where SN is the signal in the nth time slice + + \author E. Garcia CSU & J. Gomez UMD +*/ + +class ZdcSimpleRecAlgo_Run3 { +public: + /** Simple constructor for PMT-based detectors */ + ZdcSimpleRecAlgo_Run3(int recoMethod); + void initCorrectionMethod(const int method, const int ZdcSection); + void initTemplateFit(const std::vector& bxTs, + const std::vector& chargeRatios, + const int nTs, + const int ZdcSection); + void initRatioSubtraction(const float ratio, const float frac, const int ZdcSection); + + ZDCRecHit reco0(const QIE10DataFrame& digi, + const HcalCoder& coder, + const HcalCalibrations& calibs, + const HcalPedestal& effPeds, + const std::vector& myNoiseTS, + const std::vector& mySignalTS) const; + // reco method currently used to match L1 Trigger LUT energy formula + ZDCRecHit reconstruct(const QIE10DataFrame& digi, + const std::vector& myNoiseTS, + const std::vector& mySignalTS, + const HcalCoder& coder, + const HcalCalibrations& calibs, + const HcalPedestal& effPeds) const; + +private: + int recoMethod_; + int nTs_; + std::map> templateFitValues_; // Values[ZdcSection][Ts] + std::map templateFitValid_; // Values[ZdcSection] + std::map ootpuRatio_; // Values[ZdcSection] + std::map ootpuFrac_; // Values[ZdcSection] + std::map correctionMethod_; // Values[ZdcSection] +}; + +#endif diff --git a/RecoLocalCalo/HcalRecAlgos/src/ZdcSimpleRecAlgo_Run3.cc b/RecoLocalCalo/HcalRecAlgos/src/ZdcSimpleRecAlgo_Run3.cc new file mode 100644 index 0000000000000..cfd0bb26dcd23 --- /dev/null +++ b/RecoLocalCalo/HcalRecAlgos/src/ZdcSimpleRecAlgo_Run3.cc @@ -0,0 +1,270 @@ +#include "RecoLocalCalo/HcalRecAlgos/interface/ZdcSimpleRecAlgo_Run3.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include // for "max" +#include +#include + +#include // for TemplateFit Method + +// #include "DataFormats/METReco/interface/HcalCaloFlagLabels.h" + +ZdcSimpleRecAlgo_Run3::ZdcSimpleRecAlgo_Run3(int recoMethod) : recoMethod_(recoMethod) {} + +void ZdcSimpleRecAlgo_Run3::initCorrectionMethod(const int method, const int ZdcSection) { + correctionMethod_[ZdcSection] = method; +}; + +// Template fit is a linear combination of timeslices in a digi assuming there is a potential charge peak at each Bx +// and the charges of the Ts follwing a peak are consistent with the chargeRatios +void ZdcSimpleRecAlgo_Run3::initTemplateFit(const std::vector& bxTs, + const std::vector& chargeRatios, + const int nTs, + const int ZdcSection) { + int nRatios = chargeRatios.size(); + int nBx = bxTs.size(); + int nCols = nBx + 1; + double val = 0; + int index = 0; + int timeslice = 0; + nTs_ = nTs; + Eigen::MatrixXf a(nTs, nCols); + for (int j = 0; j < nBx; j++) { + timeslice = bxTs.at(j); + for (int i = 0; i < nTs; i++) { + val = 0; + index = i - timeslice; + if (index >= 0 && index < nRatios) + val = chargeRatios.at(index); + a(i, j) = val; + } + } + for (int i = 0; i < nTs; i++) + a(i, nBx) = 1; + Eigen::MatrixXf b = a.transpose() * a; + if (std::fabs(b.determinant()) < 1E-8) { + templateFitValid_[ZdcSection] = false; + return; + } + templateFitValid_[ZdcSection] = true; + Eigen::MatrixXf tfMatrix; + tfMatrix = b.inverse() * a.transpose(); + for (int i = 0; i < nTs; i++) + templateFitValues_[ZdcSection].push_back(tfMatrix.coeff(1, i)); +} + +void ZdcSimpleRecAlgo_Run3::initRatioSubtraction(const float ratio, const float frac, const int ZdcSection) { + ootpuRatio_[ZdcSection] = ratio; + ootpuFrac_[ZdcSection] = frac; +} + +// helper functions for pedestal subtraction and noise calculations +namespace zdchelper { + + inline double subPedestal(const float charge, const float ped, const float width) { + if (charge - ped > width) + return (charge - ped); + else + return (0); + } + + // gets the noise with a check if the ratio of energy0 / energy1 > ootpuRatio + // energy1 is noiseTS , energy0 is noiseTs-1 + inline double getNoiseOOTPURatio(const float energy0, + const float energy1, + const float ootpuRatio, + const float ootpuFrac) { + if (energy0 >= ootpuRatio * energy1 || ootpuRatio < 0) + return (ootpuFrac * energy1); + else + return (energy1); + } + +} // namespace zdchelper + +ZDCRecHit ZdcSimpleRecAlgo_Run3::reco0(const QIE10DataFrame& digi, + const HcalCoder& coder, + const HcalCalibrations& calibs, + const HcalPedestal& effPeds, + const std::vector& myNoiseTS, + const std::vector& mySignalTS) const { + CaloSamples tool; + coder.adc2fC(digi, tool); + // Reads noiseTS and signalTS from database + int ifirst = mySignalTS[0]; + double ampl = 0; + int maxI = -1; + double maxA = -1e10; + double ta = 0; + double energySOIp1 = 0; + double ratioSOIp1 = -1.0; + double chargeWeightedTime = 0; + + double Allnoise = 0; + int noiseslices = 0; + int CurrentTS = 0; + double noise = 0; + int digi_size = digi.size(); + HcalZDCDetId cell = digi.id(); + int zdcsection = cell.section(); + + // determining noise + for (unsigned int iv = 0; iv < myNoiseTS.size(); ++iv) { + CurrentTS = myNoiseTS[iv]; + int capid = digi[CurrentTS].capid(); + float ped = effPeds.getValue(capid); + float width = effPeds.getWidth(capid); + float gain = calibs.respcorrgain(capid); + if (CurrentTS >= digi_size) + continue; + float energy1 = zdchelper::subPedestal(tool[CurrentTS], ped, width) * gain; + float energy0 = 0; + if (CurrentTS > 0) { + capid = digi[CurrentTS - 1].capid(); + ped = effPeds.getValue(capid); + width = effPeds.getWidth(capid); + gain = calibs.respcorrgain(capid); + energy0 = zdchelper::subPedestal(tool[CurrentTS - 1], ped, width) * gain; + } + Allnoise += zdchelper::getNoiseOOTPURatio(energy0, energy1, ootpuRatio_.at(zdcsection), ootpuFrac_.at(zdcsection)); + noiseslices++; + } + if (noiseslices != 0) { + noise = (Allnoise) / double(noiseslices); + } + + // determining signal energy and max Ts + for (unsigned int ivs = 0; ivs < mySignalTS.size(); ++ivs) { + CurrentTS = mySignalTS[ivs]; + if (CurrentTS >= digi_size) + continue; + float energy1 = -1; + int capid = digi[CurrentTS].capid(); + // float ped = calibs.pedestal(capid); + float ped = effPeds.getValue(capid); + float width = effPeds.getWidth(capid); + float gain = calibs.respcorrgain(capid); + float energy0 = std::max(0.0, zdchelper::subPedestal(tool[CurrentTS], ped, width)) * gain; + if (CurrentTS < digi_size - 1) { + capid = digi[CurrentTS].capid(); + ped = effPeds.getValue(capid); + width = effPeds.getWidth(capid); + gain = calibs.respcorrgain(capid); + energy1 = std::max(0.0, zdchelper::subPedestal(tool[CurrentTS + 1], ped, width)) * gain; + } + ta = energy0 - noise; + if (ta > 0) + ampl += ta; + + if (ta > maxA) { + ratioSOIp1 = (energy0 > 0 && energy1 > 0) ? energy0 / energy1 : -1.0; + maxA = ta; + maxI = CurrentTS; + } + } + + // determine energy if using Template Fit method + if (correctionMethod_.at(zdcsection) == 1 && templateFitValid_.at(zdcsection)) { + double energy = 0; + int digi_size = digi.size(); + for (int iv = 0; iv < nTs_; iv++) { + int capid = digi[iv].capid(); + float ped = effPeds.getValue(capid); + float width = effPeds.getWidth(capid); + float gain = calibs.respcorrgain(capid); + if (iv >= digi_size) + continue; + energy += zdchelper::subPedestal(tool[iv], ped, width) * (templateFitValues_.at(zdcsection)).at(iv) * gain; + } + ampl = std::max(0.0, energy); + } + + double time = -9999; + // Time based on regular energy + ////Cannot calculate time value with max ADC sample at first or last position in window.... + if (maxI == 0 || maxI == (tool.size() - 1)) { + LogDebug("HCAL Pulse") << "ZdcSimpleRecAlgo::reco2 :" + << " Invalid max amplitude position, " + << " max Amplitude: " << maxI << " first: " << ifirst << " last: " << (tool.size() - 1) + << std::endl; + } else { + int capid = digi[maxI - 1].capid(); + double Energy0 = std::max(0.0, ((tool[maxI - 1]) * calibs.respcorrgain(capid))); + + capid = digi[maxI].capid(); + double Energy1 = std::max(0.0, ((tool[maxI]) * calibs.respcorrgain(capid))); + capid = digi[maxI + 1].capid(); + double Energy2 = std::max(0.0, ((tool[maxI + 1]) * calibs.respcorrgain(capid))); + + double TSWeightEnergy = ((maxI - 1) * Energy0 + maxI * Energy1 + (maxI + 1) * Energy2); + double EnergySum = Energy0 + Energy1 + Energy2; + double AvgTSPos = 0.; + if (EnergySum != 0) + AvgTSPos = TSWeightEnergy / EnergySum; + // If time is zero, set it to the "nonsensical" -99 + // Time should be between 75ns and 175ns (Timeslices 3-7) + if (AvgTSPos == 0) { + time = -99; + } else { + time = (AvgTSPos * 25.0); + } + } + + // find energy for signal TS + 1 + for (unsigned int ivs = 0; ivs < mySignalTS.size(); ++ivs) { + CurrentTS = mySignalTS[ivs] + 1; + if (CurrentTS >= digi_size) + continue; + int capid = digi[CurrentTS].capid(); + // ta = tool[CurrentTS] - noise; + ta = tool[CurrentTS]; + ta *= calibs.respcorrgain(capid); // fC --> GeV + if (ta > 0) + energySOIp1 += ta; + } + + double tmp_energy = 0; + double tmp_TSWeightedEnergy = 0; + for (int ts = 0; ts < nTs_; ++ts) { + if (CurrentTS >= digi_size) + continue; + int capid = digi[ts].capid(); + + // max sure there are no negative values in time calculation + ta = std::max(0.0, tool[ts]); + ta *= calibs.respcorrgain(capid); // fC --> GeV + if (ta > 0) { + tmp_energy += ta; + tmp_TSWeightedEnergy += (ts + 1) * ta; + } + } + + chargeWeightedTime = (tmp_TSWeightedEnergy / tmp_energy - 1) * 25.0; + auto rh = ZDCRecHit(digi.id(), ampl, time, -99); + rh.setEnergySOIp1(energySOIp1); + + if (maxI >= 0 && maxI < tool.size()) { + float tmp_tdctime = 0; + // TDC error codes will be 60=-1, 61 = -2, 62 = -3, 63 = -4 + if (digi[maxI].le_tdc() >= 60) + tmp_tdctime = -1 * (digi[maxI].le_tdc() - 59); + else + tmp_tdctime = maxI * 25. + (digi[maxI].le_tdc() / 2); + rh.setTDCtime(tmp_tdctime); + } + + rh.setChargeWeightedTime(chargeWeightedTime); + rh.setRatioSOIp1(ratioSOIp1); + return rh; +} + +ZDCRecHit ZdcSimpleRecAlgo_Run3::reconstruct(const QIE10DataFrame& digi, + const std::vector& myNoiseTS, + const std::vector& mySignalTS, + const HcalCoder& coder, + const HcalCalibrations& calibs, + const HcalPedestal& effPeds) const { + return ZdcSimpleRecAlgo_Run3::reco0(digi, coder, calibs, effPeds, myNoiseTS, mySignalTS); + + edm::LogError("ZDCSimpleRecAlgoImpl::reconstruct, recoMethod was not declared"); + throw cms::Exception("ZDCSimpleRecoAlgos::exiting process"); +} \ No newline at end of file diff --git a/RecoLocalCalo/HcalRecProducers/plugins/BuildFile.xml b/RecoLocalCalo/HcalRecProducers/plugins/BuildFile.xml index 351449b0d24e1..261576cdcc298 100644 --- a/RecoLocalCalo/HcalRecProducers/plugins/BuildFile.xml +++ b/RecoLocalCalo/HcalRecProducers/plugins/BuildFile.xml @@ -14,6 +14,10 @@ + + + + diff --git a/RecoLocalCalo/HcalRecProducers/plugins/ZdcHitReconstructor_Run3.cc b/RecoLocalCalo/HcalRecProducers/plugins/ZdcHitReconstructor_Run3.cc new file mode 100644 index 0000000000000..79bcbcd98871d --- /dev/null +++ b/RecoLocalCalo/HcalRecProducers/plugins/ZdcHitReconstructor_Run3.cc @@ -0,0 +1,242 @@ +#include "ZdcHitReconstructor_Run3.h" +#include "DataFormats/Common/interface/EDCollection.h" +#include "DataFormats/Common/interface/Handle.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" +#include "CalibFormats/HcalObjects/interface/HcalCoderDb.h" +#include "CalibFormats/HcalObjects/interface/HcalCalibrations.h" +#include "CalibFormats/HcalObjects/interface/HcalDbService.h" +#include "CalibFormats/HcalObjects/interface/HcalDbRecord.h" +#include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputer.h" +#include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputerRcd.h" +#include "Geometry/CaloTopology/interface/HcalTopology.h" + +#include + +#include + +#include +namespace zdchelper { + void setZDCSaturation(ZDCRecHit rh, QIE10DataFrame& digi, int maxValue) { + for (auto it = digi.begin(); it != digi.end(); it++) { + QIE10DataFrame::Sample sample = QIE10DataFrame::Sample(*it); + if (sample.adc() >= maxValue) { + rh.setFlagField(1, HcalCaloFlagLabels::ADCSaturationBit); + break; + } + } + } + +} // namespace zdchelper + +/* Zdc Hit reconstructor allows for CaloRecHits with status words */ +ZdcHitReconstructor_Run3::ZdcHitReconstructor_Run3(edm::ParameterSet const& conf) + + : reco_(conf.getParameter("recoMethod")), + saturationFlagSetter_(nullptr), + det_(DetId::Hcal), + correctionMethodEM_(conf.getParameter("correctionMethodEM")), + correctionMethodHAD_(conf.getParameter("correctionMethodHAD")), + correctionMethodRPD_(conf.getParameter("correctionMethodRPD")), + ootpuRatioEM_(conf.getParameter("ootpuRatioEM")), + ootpuRatioHAD_(conf.getParameter("ootpuRatioHAD")), + ootpuRatioRPD_(conf.getParameter("ootpuRatioRPD")), + ootpuFracEM_(conf.getParameter("ootpuFracEM")), + ootpuFracHAD_(conf.getParameter("ootpuFracHAD")), + ootpuFracRPD_(conf.getParameter("ootpuFracRPD")), + chargeRatiosEM_(conf.getParameter>("chargeRatiosEM")), + chargeRatiosHAD_(conf.getParameter>("chargeRatiosHAD")), + chargeRatiosRPD_(conf.getParameter>("chargeRatiosRPD")), + bxTs_(conf.getParameter>("bxTs")), + nTs_(conf.getParameter("nTs")), + forceSOI_(conf.getParameter("forceSOI")), + signalSOI_(conf.getParameter>("signalSOI")), + noiseSOI_(conf.getParameter>("noiseSOI")), + setSaturationFlags_(conf.getParameter("setSaturationFlags")), + dropZSmarkedPassed_(conf.getParameter("dropZSmarkedPassed")), + skipRPD_(conf.getParameter("skipRPD")) { + tok_input_QIE10 = consumes(conf.getParameter("digiLabelQIE10ZDC")); + + std::string subd = conf.getParameter("Subdetector"); + + if (setSaturationFlags_) { + const edm::ParameterSet& pssat = conf.getParameter("saturationParameters"); + maxADCvalue_ = pssat.getParameter("maxADCvalue"); + } + if (!strcasecmp(subd.c_str(), "ZDC")) { + det_ = DetId::Calo; + subdet_ = HcalZDCDetId::SubdetectorId; + produces(); + } else if (!strcasecmp(subd.c_str(), "CALIB")) { + subdet_ = HcalOther; + subdetOther_ = HcalCalibration; + produces(); + } else { + std::cout << "ZdcHitReconstructor_Run3 is not associated with a specific subdetector!" << std::endl; + } + reco_.initCorrectionMethod(correctionMethodEM_, 1); + reco_.initCorrectionMethod(correctionMethodHAD_, 2); + reco_.initCorrectionMethod(correctionMethodRPD_, 4); + reco_.initTemplateFit(bxTs_, chargeRatiosEM_, nTs_, 1); + reco_.initTemplateFit(bxTs_, chargeRatiosHAD_, nTs_, 2); + reco_.initTemplateFit(bxTs_, chargeRatiosRPD_, nTs_, 4); + reco_.initRatioSubtraction(ootpuRatioEM_, ootpuFracEM_, 1); + reco_.initRatioSubtraction(ootpuRatioHAD_, ootpuFracHAD_, 2); + reco_.initRatioSubtraction(ootpuRatioRPD_, ootpuFracRPD_, 4); + // ES tokens + htopoToken_ = esConsumes(); + paramsToken_ = esConsumes(); + conditionsToken_ = esConsumes(); + qualToken_ = esConsumes(edm::ESInputTag("", "withTopo")); + sevToken_ = esConsumes(); +} + +ZdcHitReconstructor_Run3::~ZdcHitReconstructor_Run3() { delete saturationFlagSetter_; } + +void ZdcHitReconstructor_Run3::beginRun(edm::Run const& r, edm::EventSetup const& es) { + const HcalTopology& htopo = es.getData(htopoToken_); + const HcalLongRecoParams& p = es.getData(paramsToken_); + longRecoParams_ = std::make_unique(p); + longRecoParams_->setTopo(&htopo); +} + +void ZdcHitReconstructor_Run3::endRun(edm::Run const& r, edm::EventSetup const& es) {} + +void ZdcHitReconstructor_Run3::produce(edm::Event& e, const edm::EventSetup& eventSetup) { + // get conditions + const HcalDbService* conditions = &eventSetup.getData(conditionsToken_); + const HcalChannelQuality* myqual = &eventSetup.getData(qualToken_); + const HcalSeverityLevelComputer* mySeverity = &eventSetup.getData(sevToken_); + + // define vectors to pass noiseTS and signalTS + std::vector mySignalTS; + std::vector myNoiseTS; + + if (det_ == DetId::Calo && subdet_ == HcalZDCDetId::SubdetectorId) { + edm::Handle digi; + e.getByToken(tok_input_QIE10, digi); + + // create empty output + auto rec = std::make_unique(); + rec->reserve(digi->size()); + + // testing QEI10 conditions + for (auto it = digi->begin(); it != digi->end(); it++) { + QIE10DataFrame QIE10_i = static_cast(*it); + HcalZDCDetId cell = QIE10_i.id(); + bool isRPD = cell.section() == 4; + if (isRPD && skipRPD_) + continue; + if (cell.section() == 1 && cell.channel() > 5) + continue; // ignore extra EM channels + + DetId detcell = (DetId)cell; + + // check on cells to be ignored and dropped: (rof,20.Feb.09) + const HcalChannelStatus* mydigistatus = myqual->getValues(detcell.rawId()); + if (mySeverity->dropChannel(mydigistatus->getValue())) + continue; + if (dropZSmarkedPassed_) + if (QIE10_i.zsMarkAndPass()) + continue; + + const HcalCalibrations& calibrations = conditions->getHcalCalibrations(cell); + const HcalQIECoder* channelCoder = conditions->getHcalCoder(cell); + const HcalQIEShape* shape = conditions->getHcalShape(channelCoder); + HcalCoderDb coder(*channelCoder, *shape); + + // pass the effective pedestals to rec hit since both ped value and width used in subtraction of pedestals + const HcalPedestal* effPeds = conditions->getEffectivePedestal(cell); + + if (forceSOI_) + rec->push_back(reco_.reconstruct(QIE10_i, noiseSOI_, signalSOI_, coder, calibrations, *effPeds)); + + else { + const HcalLongRecoParam* myParams = longRecoParams_->getValues(detcell); + mySignalTS.clear(); + myNoiseTS.clear(); + mySignalTS = myParams->signalTS(); + myNoiseTS = myParams->noiseTS(); + + rec->push_back(reco_.reconstruct(QIE10_i, myNoiseTS, mySignalTS, coder, calibrations, *effPeds)); + } + // saturationFlagSetter_ doesn't work with QIE10 + // created new function zdchelper::setZDCSaturation to work with QIE10 + (rec->back()).setFlags(0); + if (setSaturationFlags_) + zdchelper::setZDCSaturation(rec->back(), QIE10_i, maxADCvalue_); + } + // return result + e.put(std::move(rec)); + } // else if (det_==DetId::Calo...) + +} // void HcalHitReconstructor::produce(...) + +void ZdcHitReconstructor_Run3::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + // zdcreco + edm::ParameterSetDescription desc; + desc.add("digiLabelQIE10ZDC", edm::InputTag("hcalDigis", "ZDC")); + desc.add("Subdetector", "ZDC"); + desc.add("dropZSmarkedPassed", true); + desc.add("skipRPD", true); + desc.add("recoMethod", 1); + desc.add("correctionMethodEM", 1); + desc.add("correctionMethodHAD", 1); + desc.add("correctionMethodRPD", 0); + desc.add("ootpuRatioEM", 3.0); + desc.add("ootpuRatioHAD", 3.0); + desc.add("ootpuRatioRPD", -1.0); + desc.add("ootpuFracEM", 1.0); + desc.add("ootpuFracHAD", 1.0); + desc.add("ootpuFracRPD", 0.0); + desc.add>("chargeRatiosEM", + { + 1.0, + 0.23157, + 0.10477, + 0.06312, + }); + desc.add>("chargeRatiosHAD", + { + 1.0, + 0.23157, + 0.10477, + 0.06312, + }); + desc.add>("chargeRatiosRPD", + { + 1.0, + 0.23157, + 0.10477, + 0.06312, + }); + desc.add>("bxTs", + { + 0, + 2, + 4, + }); + desc.add("nTs", 6); + desc.add("forceSOI", false); + desc.add>("signalSOI", + { + 2, + }); + desc.add>("noiseSOI", + { + 1, + }); + desc.add("setSaturationFlags", false); + { + edm::ParameterSetDescription psd0; + psd0.add("maxADCvalue", 255); + desc.add("saturationParameters", psd0); + } + descriptions.add("zdcrecoRun3", desc); + // or use the following to generate the label from the module's C++ type + //descriptions.addWithDefaultLabel(desc); +} + +//define this as a plug-in +DEFINE_FWK_MODULE(ZdcHitReconstructor_Run3); \ No newline at end of file diff --git a/RecoLocalCalo/HcalRecProducers/plugins/ZdcHitReconstructor_Run3.h b/RecoLocalCalo/HcalRecProducers/plugins/ZdcHitReconstructor_Run3.h new file mode 100644 index 0000000000000..c46167e558c77 --- /dev/null +++ b/RecoLocalCalo/HcalRecProducers/plugins/ZdcHitReconstructor_Run3.h @@ -0,0 +1,99 @@ +#ifndef ZDCHITRECONSTRUCTOR_RUN3_H +#define ZDCHITRECONSTRUCTOR_RUN3_H 1 + +#include + +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/Framework/interface/Event.h" +#include "DataFormats/Common/interface/Handle.h" + +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Utilities/interface/ESGetToken.h" + +#include "RecoLocalCalo/HcalRecAlgos/interface/ZdcSimpleRecAlgo_Run3.h" +#include "RecoLocalCalo/HcalRecAlgos/interface/HcalHFStatusBitFromRecHits.h" +#include "RecoLocalCalo/HcalRecAlgos/interface/HcalHFStatusBitFromDigis.h" +#include "DataFormats/METReco/interface/HcalCaloFlagLabels.h" +#include "CondFormats/HcalObjects/interface/HcalChannelQuality.h" +#include "CondFormats/HcalObjects/interface/HcalChannelStatus.h" +#include "CondFormats/HcalObjects/interface/HcalLongRecoParams.h" +#include "CondFormats/HcalObjects/interface/HcalLongRecoParam.h" +#include "CondFormats/HcalObjects/interface/HcalPedestal.h" +#include "RecoLocalCalo/HcalRecAlgos/interface/HcalTimingCorrector.h" +#include "RecoLocalCalo/HcalRecAlgos/interface/HBHETimeProfileStatusBitSetter.h" +#include "RecoLocalCalo/HcalRecAlgos/interface/HBHETimingShapedFlag.h" +#include "RecoLocalCalo/HcalRecAlgos/interface/HcalADCSaturationFlag.h" +#include "RecoLocalCalo/HcalRecAlgos/interface/HFTimingTrustFlag.h" + +#include "DataFormats/HcalDigi/interface/HcalDigiCollections.h" +#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h" + +class HcalTopology; +class HcalRecNumberingRecord; +class HcalLongRecoParamsRcd; +class HcalDbService; +class HcalDbRecord; +class HcalChannelQuality; +class HcalChannelQualityRcd; +class HcalSeverityLevelComputer; +class HcalSeverityLevelComputerRcd; + +/** \class ZdcHitReconstructor_Run3 + + \author M. Nickel - Kansas + ** Based on ZDCHitReconstructor.h by E. Garcia + */ +class ZdcHitReconstructor_Run3 : public edm::stream::EDProducer<> { +public: + explicit ZdcHitReconstructor_Run3(const edm::ParameterSet& ps); + ~ZdcHitReconstructor_Run3() override; + void beginRun(edm::Run const& r, edm::EventSetup const& es) final; + void endRun(edm::Run const& r, edm::EventSetup const& es) final; + void produce(edm::Event& e, const edm::EventSetup& c) final; + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + +private: + ZdcSimpleRecAlgo_Run3 reco_; + HcalADCSaturationFlag* saturationFlagSetter_; + + DetId::Detector det_; + int subdet_; + HcalOtherSubdetector subdetOther_; + + edm::EDGetTokenT tok_input_QIE10; + + int correctionMethodEM_; + int correctionMethodHAD_; + int correctionMethodRPD_; + double ootpuRatioEM_; + double ootpuRatioHAD_; + double ootpuRatioRPD_; + double ootpuFracEM_; + double ootpuFracHAD_; + double ootpuFracRPD_; + std::vector chargeRatiosEM_; + std::vector chargeRatiosHAD_; + std::vector chargeRatiosRPD_; + std::vector bxTs_; + int nTs_; + bool forceSOI_; + std::vector signalSOI_; + std::vector noiseSOI_; + + bool setSaturationFlags_; // turn on/off flag indicating ADC saturation + bool dropZSmarkedPassed_; // turn on/off dropping of zero suppression marked and passed digis + bool skipRPD_; // ignore all channels but EM and HCAL if true + int maxADCvalue_; // max adc value for saturation Flag + + std::unique_ptr longRecoParams_; //noiseTS and signalTS from db + + // ES tokens + edm::ESGetToken htopoToken_; + edm::ESGetToken paramsToken_; + edm::ESGetToken conditionsToken_; + edm::ESGetToken qualToken_; + edm::ESGetToken sevToken_; +}; + +#endif