Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

HCAL Pulse Filtering for OOT PU Subtraction #31661

Merged
merged 11 commits into from
Nov 3, 2020
Merged

HCAL Pulse Filtering for OOT PU Subtraction #31661

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
741dd14
Initial commit of implementing weights method for OOT PU subtraction …
Apr 21, 2020
8621504
Actually pass the number of samples and presamples to the algo.
Apr 22, 2020
d00dbc0
Convert string to int when using as key
May 5, 2020
001bce1
Add ability to switch pulse containment correction from python as wel…
Apr 23, 2020
f22dd5e
Bring in ability to configure pulse containment things from python
Jun 8, 2020
f2b1fc8
modify how weights are used allowing mixed algos between HB and HE
Jun 30, 2020
577073e
expand to HB and HE versions of phase and contain
Jun 30, 2020
16b8d7a
refer to everything in ieta16 ring as HB
Jun 30, 2020
fe6f319
cleanup code formatting
Oct 5, 2020
90fe405
make more explicit when using 1TS vs 2TS containment
Oct 14, 2020
167641d
differentiate TP npresamples from npresamples used by FIR filter
Oct 16, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 6 additions & 0 deletions CalibCalorimetry/HcalTPGAlgos/interface/HcaluLUTTPGCoder.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,6 +65,10 @@ class HcaluLUTTPGCoder : public HcalTPGCoder {
linearLSB_QIE11_ = lsb11;
linearLSB_QIE11Overlap_ = lsb11overlap;
};
void set1TSContainHB(bool contain1TSHB) { contain1TSHB_ = contain1TSHB; }
void set1TSContainHE(bool contain1TSHE) { contain1TSHE_ = contain1TSHE; }
void setContainPhaseHB(double containPhaseNSHB) { containPhaseNSHB_ = containPhaseNSHB; }
void setContainPhaseHE(double containPhaseNSHE) { containPhaseNSHE_ = containPhaseNSHE; }
void lookupMSB(const HBHEDataFrame& df, std::vector<bool>& msb) const;
void lookupMSB(const QIE10DataFrame& df, std::vector<std::bitset<2>>& msb) const;
void lookupMSB(const QIE11DataFrame& df, std::vector<std::bitset<2>>& msb) const;
Expand Down Expand Up @@ -109,6 +113,8 @@ class HcaluLUTTPGCoder : public HcalTPGCoder {
// edge cases not covered by the cosh_ieta_ map
double cosh_ieta_28_HE_low_depths_, cosh_ieta_28_HE_high_depths_, cosh_ieta_29_HE_;
bool allLinear_;
bool contain1TSHB_, contain1TSHE_;
double containPhaseNSHB_, containPhaseNSHE_;
double linearLSB_QIE8_, linearLSB_QIE11_, linearLSB_QIE11Overlap_;
std::unique_ptr<HcalPulseContainmentManager> pulseCorr_;
};
Expand Down
42 changes: 32 additions & 10 deletions CalibCalorimetry/HcalTPGAlgos/src/HcaluLUTTPGCoder.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,6 +61,8 @@ HcaluLUTTPGCoder::HcaluLUTTPGCoder()
cosh_ieta_28_HE_high_depths_{},
cosh_ieta_29_HE_{},
allLinear_{},
contain1TSHB_{},
contain1TSHE_{},
linearLSB_QIE8_{},
linearLSB_QIE11_{},
linearLSB_QIE11Overlap_{} {}
Expand All @@ -74,6 +76,8 @@ void HcaluLUTTPGCoder::init(const HcalTopology* top, const HcalTimeSlew* delay)
FG_HF_thresholds_ = {0, 0};
bitToMask_ = 0;
allLinear_ = false;
contain1TSHB_ = false;
contain1TSHE_ = false;
linearLSB_QIE8_ = 1.;
linearLSB_QIE11_ = 1.;
linearLSB_QIE11Overlap_ = 1.;
Expand Down Expand Up @@ -414,12 +418,20 @@ void HcaluLUTTPGCoder::update(const HcalDbService& conditions) {
int granularity = meta->getLutGranularity();

double correctionPhaseNS = conditions.getHcalRecoParam(cell)->correctionPhaseNS();

// When containPhaseNS is not -999.0, and for QIE11 only, override from configuration
if (qieType == QIE11) {
if (containPhaseNSHB_ != -999.0 and cell.ietaAbs() <= topo_->lastHBRing())
correctionPhaseNS = containPhaseNSHB_;
else if (containPhaseNSHE_ != -999.0 and cell.ietaAbs() > topo_->lastHBRing())
correctionPhaseNS = containPhaseNSHE_;
}
for (unsigned int adc = 0; adc < SIZE; ++adc) {
if (isMasked)
lut[adc] = 0;
else {
double nonlinearityCorrection = 1.0;
double containmentCorrection2TSCorrected = 1.0;
double containmentCorrection = 1.0;
// SiPM nonlinearity was not corrected in 2017
// and containment corrections were not
// ET-dependent prior to 2018
Expand All @@ -428,28 +440,38 @@ void HcaluLUTTPGCoder::update(const HcalDbService& conditions) {
// Use the 1-TS containment correction to estimate the charge of the pulse
// from the individual samples
double correctedCharge = containmentCorrection1TS * adc2fC(adc);
containmentCorrection2TSCorrected = pulseCorr_->correction(cell, 2, correctionPhaseNS, correctedCharge);
double containmentCorrection2TSCorrected =
pulseCorr_->correction(cell, 2, correctionPhaseNS, correctedCharge);
if (qieType == QIE11) {
// When contain1TS_ is set, it should still only apply for QIE11-related things
if ((contain1TSHB_ and cell.ietaAbs() <= topo_->lastHBRing()) or
(contain1TSHE_ and cell.ietaAbs() > topo_->lastHBRing())) {
containmentCorrection = containmentCorrection1TS;
} else {
containmentCorrection = containmentCorrection2TSCorrected;
}

const HcalSiPMParameter& siPMParameter(*conditions.getHcalSiPMParameter(cell));
HcalSiPMnonlinearity corr(
conditions.getHcalSiPMCharacteristics()->getNonLinearities(siPMParameter.getType()));
const double fcByPE = siPMParameter.getFCByPE();
const double effectivePixelsFired = correctedCharge / fcByPE;
nonlinearityCorrection = corr.getRecoCorrectionFactor(effectivePixelsFired);
} else {
containmentCorrection = containmentCorrection2TSCorrected;
}
}
if (allLinear_)
lut[adc] = (LutElement)std::min(
std::max(0,
int((adc2fC(adc) - ped) * gain * rcalib * nonlinearityCorrection * containmentCorrection /
linearLSB / cosh_ieta(cell.ietaAbs(), cell.depth(), HcalEndcap))),
MASK);
else
lut[adc] = (LutElement)std::min(std::max(0,
int((adc2fC(adc) - ped) * gain * rcalib * nonlinearityCorrection *
containmentCorrection2TSCorrected / linearLSB /
cosh_ieta(cell.ietaAbs(), cell.depth(), HcalEndcap))),
containmentCorrection / nominalgain_ / granularity)),
MASK);
else
lut[adc] =
(LutElement)std::min(std::max(0,
int((adc2fC(adc) - ped) * gain * rcalib * nonlinearityCorrection *
containmentCorrection2TSCorrected / nominalgain_ / granularity)),
MASK);

if (qieType == QIE11) {
if (adc >= mipMin and adc < mipMax)
Expand Down
11 changes: 11 additions & 0 deletions CalibCalorimetry/HcalTPGEventSetup/src/HcalTPGCoderULUT.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -58,6 +58,8 @@ class HcalTPGCoderULUT : public edm::ESProducer {
edm::ESGetToken<HcalTimeSlew, HcalTimeSlewRecord> delayToken_;
edm::ESGetToken<HcalDbService, HcalDbRecord> serviceToken_;
bool read_FGLut_, read_Ascii_, read_XML_, LUTGenerationMode_, linearLUTs_;
bool contain1TSHB_, contain1TSHE_;
double containPhaseNSHB_, containPhaseNSHE_;
double linearLSB_QIE8_, linearLSB_QIE11Overlap_, linearLSB_QIE11_;
int maskBit_;
std::vector<uint32_t> FG_HF_thresholds_;
Expand All @@ -80,6 +82,10 @@ HcalTPGCoderULUT::HcalTPGCoderULUT(const edm::ParameterSet& iConfig) {
read_XML_ = iConfig.getParameter<bool>("read_XML_LUTs");
read_FGLut_ = iConfig.getParameter<bool>("read_FG_LUTs");
fgfile_ = iConfig.getParameter<edm::FileInPath>("FGLUTs");
contain1TSHB_ = iConfig.getParameter<bool>("contain1TSHB");
contain1TSHE_ = iConfig.getParameter<bool>("contain1TSHE");
containPhaseNSHB_ = iConfig.getParameter<double>("containPhaseNSHB");
containPhaseNSHE_ = iConfig.getParameter<double>("containPhaseNSHE");

//the following line is needed to tell the framework what
// data is being produced
Expand All @@ -105,6 +111,11 @@ HcalTPGCoderULUT::HcalTPGCoderULUT(const edm::ParameterSet& iConfig) {
void HcalTPGCoderULUT::buildCoder(const HcalTopology* topo, const HcalTimeSlew* delay, HcaluLUTTPGCoder* theCoder) {
using namespace edm::es;
theCoder->init(topo, delay);
theCoder->set1TSContainHB(contain1TSHB_);
theCoder->set1TSContainHE(contain1TSHE_);
theCoder->setContainPhaseHB(containPhaseNSHB_);
theCoder->setContainPhaseHE(containPhaseNSHE_);

if (read_Ascii_ || read_XML_) {
edm::LogInfo("HCAL") << "Using ASCII/XML LUTs" << ifilename_.fullPath() << " for HcalTPGCoderULUT initialization";
if (read_Ascii_) {
Expand Down
6 changes: 6 additions & 0 deletions SimCalorimetry/HcalTrigPrimAlgos/interface/HcalTriggerPrimitiveAlgo.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,8 @@ class HcalTriggerPrimitiveAlgo {
uint32_t ZS_threshold,
int numberOfSamples,
int numberOfPresamples,
int numberOfFilterPresamplesHBQIE11,
int numberOfFilterPresamplesHEQIE11,
int numberOfSamplesHF,
int numberOfPresamplesHF,
bool useTDCInMinBiasBits,
Expand Down Expand Up @@ -75,6 +77,7 @@ class HcalTriggerPrimitiveAlgo {
const HcalElectronicsMap* emap,
HcalTrigPrimDigiCollection& result);
void setPeakFinderAlgorithm(int algo);
void setWeightsQIE11(const edm::ParameterSet& weightsQIE11);
void setNCTScaleShift(int);
void setRCTScaleShift(int);

Expand Down Expand Up @@ -121,13 +124,16 @@ class HcalTriggerPrimitiveAlgo {
double theThreshold;
bool peakfind_;
std::vector<double> weights_;
std::map<int, std::vector<double>> weightsQIE11_;
int latency_;
uint32_t FG_threshold_;
std::vector<uint32_t> FG_HF_thresholds_;
uint32_t ZS_threshold_;
int ZS_threshold_I_;
int numberOfSamples_;
int numberOfPresamples_;
int numberOfFilterPresamplesHBQIE11_;
int numberOfFilterPresamplesHEQIE11_;
int numberOfSamplesHF_;
int numberOfPresamplesHF_;
bool useTDCInMinBiasBits_;
Expand Down
78 changes: 53 additions & 25 deletions SimCalorimetry/HcalTrigPrimAlgos/src/HcalTriggerPrimitiveAlgo.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,8 @@ HcalTriggerPrimitiveAlgo::HcalTriggerPrimitiveAlgo(bool pf,
uint32_t ZS_threshold,
int numberOfSamples,
int numberOfPresamples,
int numberOfFilterPresamplesHBQIE11,
int numberOfFilterPresamplesHEQIE11,
int numberOfSamplesHF,
int numberOfPresamplesHF,
bool useTDCInMinBiasBits,
Expand All @@ -44,6 +46,8 @@ HcalTriggerPrimitiveAlgo::HcalTriggerPrimitiveAlgo(bool pf,
ZS_threshold_(ZS_threshold),
numberOfSamples_(numberOfSamples),
numberOfPresamples_(numberOfPresamples),
numberOfFilterPresamplesHBQIE11_(numberOfFilterPresamplesHBQIE11),
numberOfFilterPresamplesHEQIE11_(numberOfFilterPresamplesHEQIE11),
numberOfSamplesHF_(numberOfSamplesHF),
numberOfPresamplesHF_(numberOfPresamplesHF),
useTDCInMinBiasBits_(useTDCInMinBiasBits),
Expand Down Expand Up @@ -369,26 +373,49 @@ void HcalTriggerPrimitiveAlgo::analyze(IntegerCaloSamples& samples, HcalTriggerP
void HcalTriggerPrimitiveAlgo::analyzeQIE11(IntegerCaloSamples& samples,
HcalTriggerPrimitiveDigi& result,
const HcalFinegrainBit& fg_algo) {
int shrink = weights_.size() - 1;
HcalDetId detId(samples.id());

// Get the |ieta| for current sample
unsigned int theIeta = detId.ietaAbs();

unsigned int dgSamples = samples.size();
unsigned int dgPresamples = samples.presamples();

unsigned int tpSamples = numberOfSamples_;
unsigned int tpPresamples = numberOfPresamples_;

unsigned int filterSamples = weightsQIE11_[theIeta].size();

unsigned int shift = dgPresamples - tpPresamples;

// shrink keeps the FIR filter from going off the end of the 8TS vector
unsigned int shrink = filterSamples - 1;

auto& msb = fgUpgradeMap_[samples.id()];
IntegerCaloSamples sum(samples.id(), samples.size());

HcalDetId detId(samples.id());
std::vector<HcalTrigTowerDetId> ids = theTrigTowerGeometry->towerIds(detId);
//slide algo window
for (int ibin = 0; ibin < int(samples.size()) - shrink; ++ibin) {
for (unsigned int ibin = 0; ibin < dgSamples - shrink; ++ibin) {
int algosumvalue = 0;
for (unsigned int i = 0; i < weights_.size(); i++) {
for (unsigned int i = 0; i < filterSamples; i++) {
//add up value * scale factor
// In addition, divide by two in the 10 degree phi segmentation region
// to mimic 5 degree segmentation for the trigger
unsigned int sample = samples[ibin + i];
if (sample > QIE11_MAX_LINEARIZATION_ET)
sample = QIE11_MAX_LINEARIZATION_ET;
if (ids.size() == 2)
algosumvalue += int(sample * 0.5 * weights_[i]);
else
algosumvalue += int(sample * weights_[i]);

// Usually use a segmentation factor of 1.0 but for ieta >= 21 use 0.5
double segmentationFactor = 1.0;
if (ids.size() == 2) {
segmentationFactor = 0.5;
}

// Based on the |ieta| of the sample, retrieve the correct region weight
double theWeight = weightsQIE11_[theIeta][i];

algosumvalue += int(sample * segmentationFactor * theWeight);
}
if (algosumvalue < 0)
sum[ibin] = 0; // low-side
Expand All @@ -398,30 +425,23 @@ void HcalTriggerPrimitiveAlgo::analyzeQIE11(IntegerCaloSamples& samples,
sum[ibin] = algosumvalue; //assign value to sum[]
}

// Align digis and TP
int dgPresamples = samples.presamples();
int tpPresamples = numberOfPresamples_;
int shift = dgPresamples - tpPresamples;
int dgSamples = samples.size();
int tpSamples = numberOfSamples_;

if ((shift < shrink) || (shift + tpSamples + shrink > dgSamples - (peak_finder_algorithm_ - 1))) {
edm::LogInfo("HcalTriggerPrimitiveAlgo::analyze") << "TP presample or size from the configuration file is out of "
"the accessible range. Using digi values from data instead...";
shift = shrink;
tpPresamples = dgPresamples - shrink;
tpSamples = dgSamples - (peak_finder_algorithm_ - 1) - shrink - shift;
}

std::vector<int> finegrain(tpSamples, false);

IntegerCaloSamples output(samples.id(), tpSamples);
output.setPresamples(tpPresamples);

for (int ibin = 0; ibin < tpSamples; ++ibin) {
for (unsigned int ibin = 0; ibin < tpSamples; ++ibin) {
// ibin - index for output TP
// idx - index for samples + shift
// idx - index for samples + shift - tpPresamples
// Subtract tpPresamples one more time to get SOI in the right position
int idx = ibin + shift;

// When idx is <= 0 peakfind would compare out-of-bounds of the vector. Avoid this ambiguity
if (idx <= 0) {
output[ibin] = 0;
continue;
}

bool isPeak = (sum[idx] > sum[idx - 1] && sum[idx] >= sum[idx + 1] && sum[idx] > theThreshold);

if (isPeak) {
Expand Down Expand Up @@ -834,6 +854,14 @@ void HcalTriggerPrimitiveAlgo::addUpgradeFG(const HcalTrigTowerDetId& id,
}
}

void HcalTriggerPrimitiveAlgo::setWeightsQIE11(const edm::ParameterSet& weightsQIE11) {
// Names are just abs(ieta) for HBHE
std::vector<std::string> ietaStrs = weightsQIE11.getParameterNames();
for (auto& ietaStr : ietaStrs) {
weightsQIE11_[std::stoi(ietaStr)] = weightsQIE11.getUntrackedParameter<std::vector<double>>(ietaStr);
}
}

void HcalTriggerPrimitiveAlgo::setPeakFinderAlgorithm(int algo) {
if (algo <= 0 || algo > 2)
throw cms::Exception("ERROR: Only algo 1 & 2 are supported.") << std::endl;
Expand Down
4 changes: 4 additions & 0 deletions SimCalorimetry/HcalTrigPrimProducers/python/hcaltpdigi_cff.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,10 @@
read_FG_LUTs = cms.bool(False),
LUTGenerationMode = cms.bool(True),
linearLUTs = cms.bool(False),
contain1TSHB = cms.bool(False),
contain1TSHE = cms.bool(False),
containPhaseNSHE = cms.double(-999.0),
containPhaseNSHB = cms.double(-999.0),
tpScales = tpScales,
MaskBit = cms.int32(0x8000),
FG_HF_thresholds = cms.vuint32(17, 255),
Expand Down
34 changes: 34 additions & 0 deletions SimCalorimetry/HcalTrigPrimProducers/python/hcaltpdigi_cfi.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,38 @@
simHcalTriggerPrimitiveDigis = cms.EDProducer("HcalTrigPrimDigiProducer",
peakFilter = cms.bool(True),
weights = cms.vdouble(1.0, 1.0), ##hardware algo
weightsQIE11 = cms.untracked.PSet(**dict([
("1", cms.untracked.vdouble(1.0, 1.0)),
("2", cms.untracked.vdouble(1.0, 1.0)),
("3", cms.untracked.vdouble(1.0, 1.0)),
("4", cms.untracked.vdouble(1.0, 1.0)),
("5", cms.untracked.vdouble(1.0, 1.0)),
("6", cms.untracked.vdouble(1.0, 1.0)),
("7", cms.untracked.vdouble(1.0, 1.0)),
("8", cms.untracked.vdouble(1.0, 1.0)),
("9", cms.untracked.vdouble(1.0, 1.0)),
("10", cms.untracked.vdouble(1.0, 1.0)),
("11", cms.untracked.vdouble(1.0, 1.0)),
("12", cms.untracked.vdouble(1.0, 1.0)),
("13", cms.untracked.vdouble(1.0, 1.0)),
("14", cms.untracked.vdouble(1.0, 1.0)),
("15", cms.untracked.vdouble(1.0, 1.0)),
("16", cms.untracked.vdouble(1.0, 1.0)),
("17", cms.untracked.vdouble(1.0, 1.0)),
("18", cms.untracked.vdouble(1.0, 1.0)),
("19", cms.untracked.vdouble(1.0, 1.0)),
("20", cms.untracked.vdouble(1.0, 1.0)),
("21", cms.untracked.vdouble(1.0, 1.0)),
("22", cms.untracked.vdouble(1.0, 1.0)),
("23", cms.untracked.vdouble(1.0, 1.0)),
("24", cms.untracked.vdouble(1.0, 1.0)),
("25", cms.untracked.vdouble(1.0, 1.0)),
("26", cms.untracked.vdouble(1.0, 1.0)),
("27", cms.untracked.vdouble(1.0, 1.0)),
("28", cms.untracked.vdouble(1.0, 1.0))
])
),

latency = cms.int32(1),
FG_threshold = cms.uint32(12), ## threshold for setting fine grain bit
FG_HF_thresholds = cms.vuint32(17, 255), ## thresholds for setting fine grain bit
Expand All @@ -24,6 +56,8 @@
numberOfPresamples = cms.int32(2),
numberOfSamplesHF = cms.int32(4),
numberOfPresamplesHF = cms.int32(2),
numberOfFilterPresamplesHBQIE11 = cms.int32(0),
numberOfFilterPresamplesHEQIE11 = cms.int32(0),
useTDCInMinBiasBits = cms.bool(False), # TDC information not used in MB fine grain bits
MinSignalThreshold = cms.uint32(0), # For HF PMT veto
PMTNoiseThreshold = cms.uint32(0), # For HF PMT veto
Expand Down
3 changes: 3 additions & 0 deletions SimCalorimetry/HcalTrigPrimProducers/src/HcalTrigPrimDigiProducer.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,8 @@ HcalTrigPrimDigiProducer::HcalTrigPrimDigiProducer(const edm::ParameterSet& ps)
ps.getParameter<uint32_t>("ZS_threshold"),
ps.getParameter<int>("numberOfSamples"),
ps.getParameter<int>("numberOfPresamples"),
ps.getParameter<int>("numberOfFilterPresamplesHBQIE11"),
ps.getParameter<int>("numberOfFilterPresamplesHEQIE11"),
ps.getParameter<int>("numberOfSamplesHF"),
ps.getParameter<int>("numberOfPresamplesHF"),
ps.getParameter<bool>("useTDCInMinBiasBits"),
Expand All @@ -43,6 +45,7 @@ HcalTrigPrimDigiProducer::HcalTrigPrimDigiProducer(const edm::ParameterSet& ps)
theAlgo_.overrideParameters(pset);
}
theAlgo_.setUpgradeFlags(upgrades[0], upgrades[1], upgrades[2]);
theAlgo_.setWeightsQIE11(ps.getUntrackedParameter<edm::ParameterSet>("weightsQIE11"));

HFEMB_ = false;
if (ps.exists("LSConfig")) {
Expand Down