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.

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

code updates to add a noise correlation factor between adjacent time slices for HB and HE channels #31064

Merged
merged 9 commits into from Aug 13, 2020
Merged
Show file tree
Hide file tree
Changes from 4 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
7 changes: 7 additions & 0 deletions DataFormats/HcalRecHit/interface/HBHEChannelInfo.h
Expand Up @@ -28,6 +28,7 @@ class HBHEChannelInfo {
darkCurrent_{0},
fcByPE_{0},
lambda_{0},
noisecorr_{0},
riseTime_{0.f},
adc_{0},
dFcPerADC_{0.f},
Expand All @@ -51,6 +52,7 @@ class HBHEChannelInfo {
darkCurrent_{0},
fcByPE_{0},
lambda_{0},
noisecorr_{0},
riseTime_{0.f},
adc_{0},
dFcPerADC_{0.f},
Expand All @@ -73,6 +75,7 @@ class HBHEChannelInfo {
darkCurrent_ = 0;
fcByPE_ = 0;
lambda_ = 0, dropped_ = true;
noisecorr_ = 0;
hasLinkError_ = false;
hasCapidError_ = false;
}
Expand All @@ -85,6 +88,7 @@ class HBHEChannelInfo {
const double darkCurrent,
const double fcByPE,
const double lambda,
const double noisecorr,
const bool linkError,
const bool capidError,
const bool dropThisChannel) {
Expand All @@ -96,6 +100,7 @@ class HBHEChannelInfo {
darkCurrent_ = darkCurrent;
fcByPE_ = fcByPE;
lambda_ = lambda, dropped_ = dropThisChannel;
noisecorr_ = noisecorr;
hasLinkError_ = linkError;
hasCapidError_ = capidError;
}
Expand Down Expand Up @@ -136,6 +141,7 @@ class HBHEChannelInfo {
constexpr double darkCurrent() const { return darkCurrent_; }
constexpr double fcByPE() const { return fcByPE_; }
constexpr double lambda() const { return lambda_; }
constexpr double noisecorr() const { return noisecorr_; }
constexpr bool isDropped() const { return dropped_; }
constexpr bool hasLinkError() const { return hasLinkError_; }
constexpr bool hasCapidError() const { return hasCapidError_; }
Expand Down Expand Up @@ -254,6 +260,7 @@ class HBHEChannelInfo {
double darkCurrent_;
double fcByPE_;
double lambda_;
double noisecorr_;

// Signal rise time from TDC in ns (if provided)
float riseTime_[MAXSAMPLES];
Expand Down
3 changes: 2 additions & 1 deletion DataFormats/HcalRecHit/src/classes_def.xml
Expand Up @@ -21,7 +21,8 @@
<class name="HFPreRecHit" ClassVersion="3">
<version ClassVersion="3" checksum="224636588"/>
</class>
<class name="HBHEChannelInfo" ClassVersion="7">
<class name="HBHEChannelInfo" ClassVersion="8">
<version ClassVersion="8" checksum="158975327"/>
<version ClassVersion="7" checksum="530224705"/>
<version ClassVersion="3" checksum="184305963"/>
<version ClassVersion="4" checksum="1052948187"/>
Expand Down
6 changes: 4 additions & 2 deletions DataFormats/HcalRecHit/test/test_hcal_reco.cu
Expand Up @@ -23,7 +23,7 @@ __global__ void kernel_test_hcal_hfqie10info() { HFQIE10Info info; }

__global__ void kernel_test_hcal_hbhechinfo(HBHEChannelInfo *other) {
HBHEChannelInfo info{true, true};
info.setChannelInfo(HcalDetId{0}, 10, 10, 10, 1, 2.0, 2.0, 2.0, false, false, false);
info.setChannelInfo(HcalDetId{0}, 10, 10, 10, 1, 2.0, 2.0, 2.0, 0.0, false, false, false);
other->setChannelInfo(info.id(),
info.recoShape(),
info.nSamples(),
Expand All @@ -32,6 +32,7 @@ __global__ void kernel_test_hcal_hbhechinfo(HBHEChannelInfo *other) {
info.darkCurrent(),
info.fcByPE(),
info.lambda(),
info.noisecorr(),
info.hasLinkError(),
info.hasCapidError(),
info.isDropped());
Expand Down Expand Up @@ -85,7 +86,7 @@ void test_hcal_hbhechinfo() {
};

HBHEChannelInfo h_info, h_info_test{true, true};
h_info_test.setChannelInfo(HcalDetId{0}, 10, 10, 10, 1, 2.0, 2.0, 2.0, false, false, false);
h_info_test.setChannelInfo(HcalDetId{0}, 10, 10, 10, 1, 2.0, 2.0, 2.0, 0.0, false, false, false);
HBHEChannelInfo *d_info;

cudaMalloc((void **)&d_info, sizeof(HBHEChannelInfo));
Expand All @@ -103,6 +104,7 @@ void test_hcal_hbhechinfo() {
assert(h_info.darkCurrent() == h_info_test.darkCurrent());
assert(h_info.fcByPE() == h_info_test.fcByPE());
assert(h_info.lambda() == h_info_test.lambda());
assert(h_info.noisecorr() == h_info_test.noisecorr());
assert(h_info.hasLinkError() == h_info_test.hasLinkError());
assert(h_info.hasCapidError() == h_info_test.hasCapidError());

Expand Down
4 changes: 2 additions & 2 deletions RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h
Expand Up @@ -116,7 +116,7 @@ class MahiFit {

void phase1Debug(const HBHEChannelInfo& channelData, MahiDebugInfo& mdi) const;

void doFit(std::array<float, 3>& correctedOutput, const int nbx) const;
void doFit(std::array<float, 3>& correctedOutput, const int nbx, const double noisecorr) const;

void setPulseShapeTemplate(int pulseShapeId,
const HcalPulseShapes& ps,
Expand All @@ -130,7 +130,7 @@ class MahiFit {
private:
typedef std::pair<int, std::shared_ptr<FitterFuncs::PulseShapeFunctor> > ShapeWithId;

const float minimize() const;
const float minimize(const double noisecorr) const;
void onePulseMinimize() const;
void updateCov(const SampleMatrix& invCovMat) const;
void resetPulseShapeTemplate(int pulseShapeId, const HcalPulseShapes& ps, unsigned int nSamples);
Expand Down
22 changes: 15 additions & 7 deletions RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc
Expand Up @@ -78,10 +78,12 @@ void MahiFit::phase1Apply(const HBHEChannelInfo& channelData,
if (iTS == nnlsWork_.tsOffset)
tstrig += amplitude;
}

const double noisecorr = channelData.noisecorr();
tsTOT *= channelData.tsGain(0);
tstrig *= channelData.tsGain(0);

std::cout << noisecorr << std::endl;
TaeunKwon marked this conversation as resolved.
Show resolved Hide resolved

useTriple = false;
if (tstrig >= ts4Thresh_ && tsTOT > 0) {
//Average pedestal width (for covariance matrix constraint)
Expand All @@ -92,13 +94,13 @@ void MahiFit::phase1Apply(const HBHEChannelInfo& channelData,

// only do pre-fit with 1 pulse if chiSq threshold is positive
if (chiSqSwitch_ > 0) {
doFit(reconstructedVals, 1);
doFit(reconstructedVals, 1, noisecorr);
if (reconstructedVals[2] > chiSqSwitch_) {
doFit(reconstructedVals, 0); //nbx=0 means use configured BXs
doFit(reconstructedVals, 0, noisecorr); //nbx=0 means use configured BXs
useTriple = true;
}
} else {
doFit(reconstructedVals, 0);
doFit(reconstructedVals, 0, noisecorr);
useTriple = true;
}
} else {
Expand All @@ -112,7 +114,7 @@ void MahiFit::phase1Apply(const HBHEChannelInfo& channelData,
chi2 = reconstructedVals[2];
}

void MahiFit::doFit(std::array<float, 3>& correctedOutput, int nbx) const {
void MahiFit::doFit(std::array<float, 3>& correctedOutput, int nbx, const double noisecorr) const {
unsigned int bxSize = 1;

if (nbx == 1) {
Expand Down Expand Up @@ -175,7 +177,7 @@ void MahiFit::doFit(std::array<float, 3>& correctedOutput, int nbx) const {
}
}

const float chiSq = minimize();
const float chiSq = minimize(noisecorr);

bool foundintime = false;
unsigned int ipulseintime = 0;
Expand Down Expand Up @@ -203,14 +205,20 @@ void MahiFit::doFit(std::array<float, 3>& correctedOutput, int nbx) const {
}
}

const float MahiFit::minimize() const {
const float MahiFit::minimize(const double noisecorr) const {
nnlsWork_.invcovp.setZero(nnlsWork_.tsSize, nnlsWork_.nPulseTot);
nnlsWork_.ampVec.setZero(nnlsWork_.nPulseTot);

SampleMatrix invCovMat;
invCovMat.setConstant(nnlsWork_.tsSize, nnlsWork_.tsSize, nnlsWork_.pedVal);
invCovMat += nnlsWork_.noiseTerms.asDiagonal();

//Add off-Diagonal components up to first order
for (unsigned int i = 1; i < nnlsWork_.tsSize; ++i) {
TaeunKwon marked this conversation as resolved.
Show resolved Hide resolved
invCovMat(i - 1, i) += noisecorr * sqrt(nnlsWork_.noiseTerms.coeff(i - 1) * nnlsWork_.noiseTerms.coeff(i));
invCovMat(i, i - 1) += noisecorr * sqrt(nnlsWork_.noiseTerms.coeff(i) * nnlsWork_.noiseTerms.coeff(i - 1));
TaeunKwon marked this conversation as resolved.
Show resolved Hide resolved
}

float oldChiSq = 9999;
float chiSq = oldChiSq;

Expand Down
2 changes: 2 additions & 0 deletions RecoLocalCalo/HcalRecProducers/src/HBHEPhase1Reconstructor.cc
Expand Up @@ -492,6 +492,7 @@ void HBHEPhase1Reconstructor::processData(const Collection& coll,
const double fcByPE = properties.siPMParameter->getFCByPE();
double darkCurrent = 0.;
double lambda = 0.;
const double noisecorr = properties.siPMParameter->getauxi2();
if (!saveEffectivePeds || saveInfos_) {
// needed for the dark current in the M2 in alternative of the effectivePed
darkCurrent = properties.siPMParameter->getDarkCurrent();
Expand Down Expand Up @@ -571,6 +572,7 @@ void HBHEPhase1Reconstructor::processData(const Collection& coll,
darkCurrent,
fcByPE,
lambda,
noisecorr,
hwerr.first,
hwerr.second,
properties.taggedBadByDb || dropByZS);
Expand Down