Skip to content

Commit

Permalink
Merge pull request #22121 from christopheralanwest/hcal-tp-all-correc…
Browse files Browse the repository at this point in the history 
…tions

Modifications to HCAL trigger LUTs
  • Loading branch information
@cmsbuild
cmsbuild committed Feb 17, 2018
2 parents 0b990b3 + bd9022a commit 5c51fbb
Show file tree
Hide file tree
Showing 10 changed files with 138 additions and 80 deletions.
1 change: 1 addition & 0 deletions CalibCalorimetry/HcalAlgos/interface/HcalPulseContainmentManager.h
View file
Expand Up @@ -13,6 +13,7 @@ class HcalPulseContainmentManager {
const HcalPulseContainmentCorrection * get(const HcalDetId & detId, int toAdd, float fixedphase_ns);

void beginRun(edm::EventSetup const & es);
void beginRun(const HcalTopology *topo, const edm::ESHandle<HcalTimeSlew>& delay, const edm::ESHandle<HcalMCParams>& mcParams, const edm::ESHandle<HcalRecoParams>& recoParams);
void endRun();

void setTimeSlew(const HcalTimeSlew* timeSlew) {
Expand Down
2 changes: 2 additions & 0 deletions CalibCalorimetry/HcalAlgos/interface/HcalPulseShapes.h
View file
Expand Up @@ -7,6 +7,7 @@
#include "CalibCalorimetry/HcalAlgos/interface/HcalPulseShape.h"
#include "DataFormats/HcalDetId/interface/HcalDetId.h"
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/ESHandle.h"

/** \class HcalPulseShapes
*
Expand All @@ -28,6 +29,7 @@ class HcalPulseShapes {
// only needed if you'll be getting shapes by DetId
void beginRun(edm::EventSetup const & es);
void endRun();
void beginRun(const HcalTopology* topo, const edm::ESHandle<HcalMCParams>& mcParams, const edm::ESHandle<HcalRecoParams>& recoParams);

const Shape& hbShape() const { return hpdShape_; }
const Shape& heShape() const { return hpdShape_; }
Expand Down
7 changes: 7 additions & 0 deletions CalibCalorimetry/HcalAlgos/src/HcalPulseContainmentManager.cc
View file
Expand Up @@ -26,6 +26,13 @@ void HcalPulseContainmentManager::endRun()
shapes_.endRun();
}

void HcalPulseContainmentManager::beginRun(const HcalTopology *topo, const edm::ESHandle<HcalTimeSlew>& delay, const edm::ESHandle<HcalMCParams>& mcParams, const edm::ESHandle<HcalRecoParams>& recoParams)
{
hcalTimeSlew_delay_ = &*delay;

shapes_.beginRun(topo, mcParams, recoParams);
}

double HcalPulseContainmentManager::correction(const HcalDetId & detId,
int toAdd, float fixedphase_ns, double fc_ampl)
{
Expand Down
9 changes: 9 additions & 0 deletions CalibCalorimetry/HcalAlgos/src/HcalPulseShapes.cc
View file
Expand Up @@ -127,6 +127,15 @@ void HcalPulseShapes::beginRun(edm::EventSetup const & es)
theRecoParams->setTopo(theTopology);
}

void HcalPulseShapes::beginRun(const HcalTopology* topo, const edm::ESHandle<HcalMCParams>& mcParams, const edm::ESHandle<HcalRecoParams>& recoParams)
{
theMCParams = new HcalMCParams(*mcParams.product());
theMCParams->setTopo(topo);

theRecoParams = new HcalRecoParams(*recoParams.product());
theRecoParams->setTopo(topo);
}


void HcalPulseShapes::endRun()
{
Expand Down
45 changes: 2 additions & 43 deletions CalibCalorimetry/HcalPlugins/src/HcalHardcodeCalibrations.cc
View file
Expand Up @@ -665,49 +665,8 @@ std::unique_ptr<HcalLutMetadata> HcalHardcodeCalibrations::produceLutMetadata (c
int granularity = 1;
int threshold = 0;

if (dbHardcode.useHEUpgrade() or dbHardcode.useHFUpgrade()) {
// Use values from 2016 as starting conditions for 2017+. These are
// averaged over the subdetectors, with the last two HE towers split
// off due to diverging correction values.
switch (cell.genericSubdet()) {
case HcalGenericDetId::HcalGenBarrel:
rcalib = 1.128;
break;
case HcalGenericDetId::HcalGenEndcap:
{
HcalDetId id(cell);
if (id.ietaAbs() >= 28)
rcalib = 1.188;
else
rcalib = 1.117;
// granularity is equal to 1 only for |ieta| == 17
if(id.ietaAbs() >= 18 && id.ietaAbs() <= 26) granularity = 2;
else if(id.ietaAbs() >=27 && id.ietaAbs() <= 29) granularity = 5;
}
break;
case HcalGenericDetId::HcalGenForward:
rcalib = 1.02;
break;
default:
break;
}

if (cell.isHcalTrigTowerDetId()) {
rcalib = 0.;
HcalTrigTowerDetId id(cell);
if(id.ietaAbs() <= 17) {
granularity = 1;
}
else if(id.ietaAbs() >= 18 && id.ietaAbs() <= 26) {
granularity = 2;
}
else if(id.ietaAbs() >= 27 && id.ietaAbs() <= 28) {
granularity = 5;
}
else {
granularity = 0;
}
}
if (cell.isHcalTrigTowerDetId()) {
rcalib = 0.;
}

HcalLutMetadatum item(cell.rawId(), rcalib, granularity, threshold);
Expand Down
12 changes: 11 additions & 1 deletion CalibCalorimetry/HcalTPGAlgos/interface/HcaluLUTTPGCoder.h
View file
Expand Up @@ -5,6 +5,7 @@
#include "CalibFormats/HcalObjects/interface/HcalNominalCoder.h"
#include "Geometry/CaloTopology/interface/HcalTopology.h"
#include "DataFormats/HcalDetId/interface/HcalDetId.h"
#include "CalibCalorimetry/HcalAlgos/interface/HcalPulseContainmentManager.h"

#include <bitset>
#include <vector>
Expand Down Expand Up @@ -32,7 +33,7 @@ class HcaluLUTTPGCoder : public HcalTPGCoder {
public:
static const float lsb_;

HcaluLUTTPGCoder(const HcalTopology* topo);
HcaluLUTTPGCoder(const HcalTopology* topo, const edm::ESHandle<HcalTimeSlew>& delay, const edm::ESHandle<HcalMCParams>& mcParams, const edm::ESHandle<HcalRecoParams>& recoParams);
~HcaluLUTTPGCoder() override;
void adc2Linear(const HBHEDataFrame& df, IntegerCaloSamples& ics) const override;
void adc2Linear(const HFDataFrame& df, IntegerCaloSamples& ics) const override;
Expand All @@ -44,6 +45,8 @@ class HcaluLUTTPGCoder : public HcalTPGCoder {
float getLUTGain(HcalDetId id) const override;
std::vector<unsigned short> getLinearizationLUT(HcalDetId id) const override;

double cosh_ieta(int ieta, int depth, HcalSubdetector subdet);
void make_cosh_ieta_map(void);
void update(const HcalDbService& conditions);
void update(const char* filename, bool appendMSB = false);
void updateXML(const char* filename);
Expand Down Expand Up @@ -80,6 +83,9 @@ class HcaluLUTTPGCoder : public HcalTPGCoder {

// member variables
const HcalTopology* topo_;
const edm::ESHandle<HcalTimeSlew>& delay_;
const edm::ESHandle<HcalMCParams>& mcParams_;
const edm::ESHandle<HcalRecoParams>& recoParams_;
bool LUTGenerationMode_;
unsigned int FG_HF_threshold_;
int bitToMask_;
Expand All @@ -89,8 +95,12 @@ class HcaluLUTTPGCoder : public HcalTPGCoder {
std::vector< Lut > inputLUT_;
std::vector<float> gain_;
std::vector<float> ped_;
std::vector<double> cosh_ieta_;
// 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_;
double linearLSB_QIE8_, linearLSB_QIE11_, linearLSB_QIE11Overlap_;
std::unique_ptr<HcalPulseContainmentManager> pulseCorr_;
};

#endif
104 changes: 78 additions & 26 deletions CalibCalorimetry/HcalTPGAlgos/src/HcaluLUTTPGCoder.cc
View file
Expand Up @@ -26,6 +26,7 @@
#include "CondFormats/HcalObjects/interface/HcalL1TriggerObject.h"
#include "CalibCalorimetry/HcalAlgos/interface/HcalDbASCIIIO.h"
#include "CalibCalorimetry/HcalAlgos/interface/HcalSiPMnonlinearity.h"
#include "CalibCalorimetry/HcalAlgos/interface/HcalPulseContainmentCorrection.h"
#include "CalibCalorimetry/HcalTPGAlgos/interface/XMLProcessor.h"
#include "CalibCalorimetry/HcalTPGAlgos/interface/LutXml.h"

Expand All @@ -35,8 +36,9 @@ const int HcaluLUTTPGCoder::QIE8_LUT_BITMASK;
const int HcaluLUTTPGCoder::QIE10_LUT_BITMASK;
const int HcaluLUTTPGCoder::QIE11_LUT_BITMASK;

constexpr double MaximumFractionalError = 0.002; // 0.2% error allowed from this source

HcaluLUTTPGCoder::HcaluLUTTPGCoder(const HcalTopology* top) : topo_(top), LUTGenerationMode_(true), bitToMask_(0), allLinear_(false), linearLSB_QIE8_(1.), linearLSB_QIE11_(1.) {
HcaluLUTTPGCoder::HcaluLUTTPGCoder(const HcalTopology* top, const edm::ESHandle<HcalTimeSlew>& delay, const edm::ESHandle<HcalMCParams>& mcParams, const edm::ESHandle<HcalRecoParams>& recoParams) : topo_(top), delay_(delay), mcParams_(mcParams), recoParams_(recoParams), LUTGenerationMode_(true), bitToMask_(0), allLinear_(false), linearLSB_QIE8_(1.), linearLSB_QIE11_(1.), pulseCorr_(std::make_unique<HcalPulseContainmentManager>(MaximumFractionalError)) {
firstHBEta_ = topo_->firstHBRing();
lastHBEta_ = topo_->lastHBRing();
nHBEta_ = (lastHBEta_-firstHBEta_+1);
Expand All @@ -56,6 +58,7 @@ HcaluLUTTPGCoder::HcaluLUTTPGCoder(const HcalTopology* top) : topo_(top), LUTGen
inputLUT_ = std::vector<HcaluLUTTPGCoder::Lut>(nluts);
gain_ = std::vector<float>(nluts, 0.);
ped_ = std::vector<float>(nluts, 0.);
make_cosh_ieta_map();
}

void HcaluLUTTPGCoder::compress(const IntegerCaloSamples& ics, const std::vector<bool>& featureBits, HcalTriggerPrimitiveDigi& tp) const {
Expand Down Expand Up @@ -224,26 +227,61 @@ void HcaluLUTTPGCoder::updateXML(const char* filename) {
XMLProcessor::getInstance()->terminate();
}

double HcaluLUTTPGCoder::cosh_ieta(int ieta, int depth, HcalSubdetector subdet) {
// ieta = 28 and 29 are both associated with trigger tower 28
// so special handling is required. HF ieta=29 channels included in TT30
// are already handled correctly in cosh_ieta_
if (abs(ieta) >= 28 && subdet == HcalEndcap && allLinear_) {
if (abs(ieta) == 29) return cosh_ieta_29_HE_;
if (abs(ieta) == 28) {
if (depth <= 3) return cosh_ieta_28_HE_low_depths_;
else return cosh_ieta_28_HE_high_depths_;
}
}

return cosh_ieta_[ieta];
}

void HcaluLUTTPGCoder::make_cosh_ieta_map(void) {

cosh_ieta_ = std::vector<double>(lastHFEta_ + 1, -1.0);

HcalTrigTowerGeometry triggeo(topo_);

for (int i = 1; i <= firstHFEta_; ++i) {
double eta_low = 0., eta_high = 0.;
triggeo.towerEtaBounds(i, 0, eta_low, eta_high);
cosh_ieta_[i] = cosh((eta_low + eta_high)/2.);
}
for (int i = firstHFEta_; i <= lastHFEta_; ++i){
std::pair<double,double> etas = topo_->etaRange(HcalForward,i);
double eta1 = etas.first;
double eta2 = etas.second;
cosh_ieta_[i] = cosh((eta1 + eta2)/2.);
}

// trigger tower 28 in HE has a more complicated geometry
std::pair<double, double> eta28 = topo_->etaRange(HcalEndcap, 28);
std::pair<double, double> eta29 = topo_->etaRange(HcalEndcap, 29);
cosh_ieta_29_HE_ = cosh((eta29.first + eta29.second)/2.);
cosh_ieta_28_HE_low_depths_ = cosh((eta28.first + eta28.second)/2.);
// for higher depths in ieta = 28, the trigger tower extends past
// the ieta = 29 channels
cosh_ieta_28_HE_high_depths_ = cosh((eta28.first + eta29.second)/2.);
}

void HcaluLUTTPGCoder::update(const HcalDbService& conditions) {

HcalCalibrations calibrations;
const HcalLutMetadata *metadata = conditions.getHcalLutMetadata();
assert(metadata !=nullptr);
float nominalgain_ = metadata->getNominalGain();

HcalTrigTowerGeometry triggeo(topo_);
std::map<int, float> cosh_ieta;
for (int i = 1; i <= firstHFEta_; ++i) {
double eta_low = 0., eta_high = 0.;
triggeo.towerEtaBounds(i, 0, eta_low, eta_high);
cosh_ieta[i] = fabs(cosh((eta_low + eta_high)/2.));
}
for (int i = firstHFEta_; i <= lastHFEta_; ++i){
std::pair<double,double> etas = topo_->etaRange(HcalForward,i);
double eta1 = etas.first;
double eta2 = etas.second;
cosh_ieta[i] = cosh((eta1 + eta2)/2.);
}
pulseCorr_->beginRun(topo_, delay_, mcParams_, recoParams_);
HcalRecoParams *theRecoParams = new HcalRecoParams(*recoParams_.product());
theRecoParams->setTopo(topo_);

make_cosh_ieta_map();

for (const auto& id: metadata->getAllChannels()) {

Expand All @@ -262,9 +300,9 @@ void HcaluLUTTPGCoder::update(const HcalDbService& conditions) {
unsigned int mipMax = 0;
unsigned int mipMin = 0;

if (topo_->triggerMode() >= HcalTopologyMode::TriggerMode_2018 or
topo_->triggerMode() == HcalTopologyMode::TriggerMode_2018legacy or
topo_->dddConstants()->isPlan1(cell)) {
bool is2018OrLater = topo_->triggerMode() >= HcalTopologyMode::TriggerMode_2018 or
topo_->triggerMode() == HcalTopologyMode::TriggerMode_2018legacy;
if (is2018OrLater or topo_->dddConstants()->isPlan1(cell)) {
const HcalTPChannelParameter *channelParameters = conditions.getHcalTPChannelParameter(cell);
mipMax = channelParameters->getFGBitInfo() >> 16;
mipMin = channelParameters->getFGBitInfo() & 0xFFFF;
Expand Down Expand Up @@ -325,21 +363,33 @@ void HcaluLUTTPGCoder::update(const HcalDbService& conditions) {

int granularity = meta->getLutGranularity();

double correctionPhaseNS = theRecoParams->getValues(cell)->correctionPhaseNS();
for (unsigned int adc = 0; adc < SIZE; ++adc) {
if (isMasked) lut[adc] = 0;
else {
double nonlinearityCorrection = 1.0;
if(qieType==QIE11) {
const HcalSiPMParameter& siPMParameter(*conditions.getHcalSiPMParameter(cell));
HcalSiPMnonlinearity corr(conditions.getHcalSiPMCharacteristics()->getNonLinearities(siPMParameter.getType()));
const double fcByPE = siPMParameter.getFCByPE();
const double effectivePixelsFired = adc2fC(adc)/fcByPE;
nonlinearityCorrection = corr.getRecoCorrectionFactor(effectivePixelsFired);
double containmentCorrection2TSCorrected = 1.0;
// SiPM nonlinearity was not corrected in 2017
// and containment corrections were not
// ET-dependent prior to 2018
if(is2018OrLater) {
double containmentCorrection1TS = pulseCorr_->correction(cell, 1, correctionPhaseNS, adc2fC(adc));
// 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);
if(qieType==QIE11) {
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);
}
}
if (allLinear_)
lut[adc] = (LutElement) std::min(std::max(0, int((adc2fC(adc) - ped) * gain * rcalib * nonlinearityCorrection / linearLSB / cosh_ieta[cell.ietaAbs()])), MASK);
lut[adc] = (LutElement) std::min(std::max(0, int((adc2fC(adc) - ped) * gain * rcalib * nonlinearityCorrection * containmentCorrection2TSCorrected / 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 / nominalgain_ / granularity)), MASK);
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) lut[adc] |= QIE11_LUT_MSB0;
Expand All @@ -352,12 +402,14 @@ void HcaluLUTTPGCoder::update(const HcalDbService& conditions) {
for (unsigned int adc = 0; adc < SIZE; ++adc) {
if (isMasked) lut[adc] = 0;
else {
lut[adc] = std::min(std::max(0,int((adc2fC(adc) - ped) * gain * rcalib / lsb_ / cosh_ieta[cell.ietaAbs()] )), MASK);
lut[adc] = std::min(std::max(0,int((adc2fC(adc) - ped) * gain * rcalib / lsb_ / cosh_ieta_[cell.ietaAbs()])), MASK);
if(adc>FG_HF_threshold_) lut[adc] |= QIE10_LUT_MSB;
}
}
}
}
pulseCorr_->endRun();
delete theRecoParams;
}

void HcaluLUTTPGCoder::adc2Linear(const HBHEDataFrame& df, IntegerCaloSamples& ics) const {
Expand Down
28 changes: 23 additions & 5 deletions CalibCalorimetry/HcalTPGEventSetup/src/HcalTPGCoderULUT.cc
View file
Expand Up @@ -47,7 +47,7 @@ class HcalTPGCoderULUT : public edm::ESProducer {

ReturnType produce(const HcalTPGRecord&);
private:
void buildCoder(const HcalTopology*);
void buildCoder(const HcalTopology*, const edm::ESHandle<HcalTimeSlew>&, const edm::ESHandle<HcalMCParams>&, const edm::ESHandle<HcalRecoParams>&);
// ----------member data ---------------------------
ReturnType coder_;
HcaluLUTTPGCoder* theCoder_;
Expand Down Expand Up @@ -96,9 +96,10 @@ HcalTPGCoderULUT::HcalTPGCoderULUT(const edm::ParameterSet& iConfig)
}


void HcalTPGCoderULUT::buildCoder(const HcalTopology* topo) {
void HcalTPGCoderULUT::buildCoder(const HcalTopology* topo, const edm::ESHandle<HcalTimeSlew>& delay,
const edm::ESHandle<HcalMCParams>& mcParams, const edm::ESHandle<HcalRecoParams>& recoParams) {
using namespace edm::es;
theCoder_ = new HcaluLUTTPGCoder(topo);
theCoder_ = new HcaluLUTTPGCoder(topo, delay, mcParams, recoParams);
if (read_Ascii_ || read_XML_){
edm::LogInfo("HCAL") << "Using ASCII/XML LUTs" << ifilename_.fullPath() << " for HcalTPGCoderULUT initialization";
if (read_Ascii_) {
Expand Down Expand Up @@ -139,7 +140,17 @@ HcalTPGCoderULUT::produce(const HcalTPGRecord& iRecord)
edm::ESHandle<HcalTopology> htopo;
iRecord.getRecord<HcalRecNumberingRecord>().get(htopo);
const HcalTopology* topo=&(*htopo);
buildCoder(topo);

edm::ESHandle<HcalTimeSlew> delay;
iRecord.getRecord<HcalDbRecord>().getRecord<HcalTimeSlewRecord>().get("HBHE", delay);

edm::ESHandle<HcalMCParams> mcParams;
iRecord.getRecord<HcalDbRecord>().getRecord<HcalMCParamsRcd>().get(mcParams);

edm::ESHandle<HcalRecoParams> recoParams;
iRecord.getRecord<HcalDbRecord>().getRecord<HcalRecoParamsRcd>().get(recoParams);

buildCoder(topo, delay, mcParams, recoParams);
}


Expand All @@ -153,7 +164,14 @@ void HcalTPGCoderULUT::dbRecordCallback(const HcalDbRecord& theRec) {
theRec.getRecord<HcalRecNumberingRecord>().get(htopo);
const HcalTopology* topo=&(*htopo);

buildCoder(topo);
edm::ESHandle<HcalTimeSlew> delay;
theRec.getRecord<HcalTimeSlewRecord>().get("HBHE", delay);
edm::ESHandle<HcalMCParams> mcParams;
theRec.getRecord<HcalMCParamsRcd>().get(mcParams);
edm::ESHandle<HcalRecoParams> recoParams;
theRec.getRecord<HcalRecoParamsRcd>().get(recoParams);

buildCoder(topo, delay, mcParams, recoParams);

theCoder_->update(*conditions);

Expand Down
4 changes: 2 additions & 2 deletions CalibFormats/HcalObjects/interface/HcalDbRecord.h
View file
Expand Up @@ -28,10 +28,10 @@
// class HcalDbRecord : public edm::eventsetup::EventSetupRecordImplementation<HcalDbRecord> {};

class HcalDbRecord : public edm::eventsetup::DependentRecordImplementation <HcalDbRecord,
boost::mpl::vector23<HcalRecNumberingRecord, IdealGeometryRecord, HcalPedestalsRcd, HcalPedestalWidthsRcd, HcalGainsRcd, HcalGainWidthsRcd,
boost::mpl::vector25<HcalRecNumberingRecord, IdealGeometryRecord, HcalPedestalsRcd, HcalPedestalWidthsRcd, HcalGainsRcd, HcalGainWidthsRcd,
HcalQIEDataRcd, HcalQIETypesRcd, HcalChannelQualityRcd, HcalZSThresholdsRcd, HcalRespCorrsRcd,
HcalL1TriggerObjectsRcd, HcalElectronicsMapRcd, HcalTimeCorrsRcd, HcalLUTCorrsRcd, HcalPFCorrsRcd,
HcalFrontEndMapRcd, HcalSiPMCharacteristicsRcd, HcalSiPMParametersRcd, HcalTPParametersRcd, HcalTPChannelParametersRcd,
HcalLutMetadataRcd, HcalMCParamsRcd > > {};
HcalLutMetadataRcd, HcalMCParamsRcd, HcalRecoParamsRcd, HcalTimeSlewRecord > > {};

#endif /* HCALDBPRODUCER_HCALDBRECORD_H */

0 comments on commit 5c51fbb

Please sign in to comment.