Skip to content

Commit

Permalink
add fit quality matrix diagnostic ME
Browse files Browse the repository at this point in the history
  • Loading branch information
@mmusich
mmusich committed Aug 21, 2022
1 parent 2f16c61 commit 216a222
Show file tree
Hide file tree
Showing 2 changed files with 60 additions and 16 deletions.
3 changes: 3 additions & 0 deletions CalibTracker/SiPixelLorentzAngle/interface/SiPixelLorentzAngleCalibrationStruct.h
View file
Expand Up @@ -85,6 +85,9 @@ struct SiPixelLorentzAngleCalibrationHistograms {
dqm::reco::MonitorElement* h_bySectCovMatrixStatus_;
dqm::reco::MonitorElement* h_bySectDriftError_;

// for fit quality
dqm::reco::MonitorElement* h_bySectFitQuality_;

// ouput LA maps
std::vector<dqm::reco::MonitorElement*> h2_byLayerLA_;
std::vector<dqm::reco::MonitorElement*> h2_byLayerDiff_;
Expand Down
73 changes: 57 additions & 16 deletions CalibTracker/SiPixelLorentzAngle/src/SiPixelLorentzAnglePCLHarvester.cc
View file
Expand Up @@ -48,6 +48,7 @@
namespace SiPixelLAHarvest {

enum covStatus { kUndefined = -1, kNotCalculated = 0, kApproximated = 1, kMadePosDef = 2, kAccurate = 3 };
enum cuts { kZeroChi2 = 0, kChi2Cut = 1, kCovStatus = 2, kNentries = 3 };

struct fitResults {
public:
Expand All @@ -61,6 +62,25 @@ namespace SiPixelLAHarvest {
quality = {0b0000};
};

void setQualityCutBit(const SiPixelLAHarvest::cuts& cut) {
switch (cut) {
case kZeroChi2:
quality.set(0);
break;
case kChi2Cut:
quality.set(1);
break;
case kCovStatus:
quality.set(2);
break;
case kNentries:
quality.set(3);
break;
default:
throw cms::Exception("Inconsistent Data") << "Passed inexistent cut value";
}
}

double p0, e0;
double p1, e1;
double p2, e2;
Expand Down Expand Up @@ -94,7 +114,7 @@ class SiPixelLorentzAnglePCLHarvester : public DQMEDHarvester {
void endRun(const edm::Run&, const edm::EventSetup&) override;
void findMean(MonitorElement* h_drift_depth_adc_slice_, int i, int i_ring);
SiPixelLAHarvest::fitResults fitAndStore(std::shared_ptr<SiPixelLorentzAngle> theLA, int i_idx, int i_lay, int i_mod);
bool checkFitQuality(SiPixelLAHarvest::fitResults& res);
bool isFitGood(SiPixelLAHarvest::fitResults& res);

// es tokens
edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geomEsToken_;
Expand Down Expand Up @@ -369,7 +389,7 @@ void SiPixelLorentzAnglePCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQMS
iBooker.book1D("h_bySectorChi2", "Fit #chi^{2}/ndf by sector;pixel sector; fit #chi^{2}/ndf", maxSect, lo, hi);

hists_.h_bySectFitStatus_ =
iBooker.book1D("h_bySectFitStatus_", "Fit Status by sector;pixel sector; fit status", maxSect, lo, hi);
iBooker.book1D("h_bySectFitStatus", "Fit Status by sector;pixel sector; fit status", maxSect, lo, hi);

hists_.h_bySectCovMatrixStatus_ = iBooker.book1D(
"h_bySectorCovMatrixStatus", "Fit Covariance Matrix Status by sector;pixel sector; fit status", maxSect, lo, hi);
Expand All @@ -380,6 +400,17 @@ void SiPixelLorentzAnglePCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQMS
lo,
hi);

// set the bin labels for the fit quality dashboard histogram
std::vector<std::string> labels = {"#chi^{2}!=0", "#chi^{2} cut", "covStat>=2", "n. entries cut"};
hists_.h_bySectFitQuality_ = iBooker.book2D("h_bySectFitQuality",
"Fit quality by sector;pixel sector;quality cut",
maxSect,
lo,
hi,
labels.size(),
-0.5,
labels.size() - 0.5);

// copy the bin labels from the occupancy histogram
for (int bin = 1; bin <= maxSect; bin++) {
const auto& binName = hists_.h_bySectOccupancy_->getTH1()->GetXaxis()->GetBinLabel(bin);
Expand All @@ -392,6 +423,11 @@ void SiPixelLorentzAnglePCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQMS
hists_.h_bySectFitStatus_->setBinLabel(bin, binName);
hists_.h_bySectCovMatrixStatus_->setBinLabel(bin, binName);
hists_.h_bySectDriftError_->setBinLabel(bin, binName);
hists_.h_bySectFitQuality_->setBinLabel(bin, binName, 1);
}

for (unsigned int binY = 1; binY <= labels.size(); binY++) {
hists_.h_bySectFitQuality_->setBinLabel(binY, labels[binY - 1], 2);
}

// this will be booked in the Harvesting folder
Expand Down Expand Up @@ -644,7 +680,7 @@ SiPixelLAHarvest::fitResults SiPixelLorentzAnglePCLHarvester::fitAndStore(
// DO NOT REMOVE
// this is needed to ensure it stays synch with the render plugin:
// https://github.com/dmwm/deployment/blob/master/dqmgui/style/SiPixelLorentzAnglePCLRenderPlugin.cc#L199
//assert(std::string{f1_->GetName()}=="f1");
// assert(std::string{f1_->GetName()}=="f1");
assert(strcmp(f1_->GetName(), "f1") == 0);

f1_->SetParName(0, "offset");
Expand Down Expand Up @@ -674,10 +710,10 @@ SiPixelLAHarvest::fitResults SiPixelLorentzAnglePCLHarvester::fitAndStore(
edm::LogInfo("SiPixelLorentzAnglePCLHarvester") << "Fit status: " << res.fitStatus;
} else {
edm::LogError("SiPixelLorentzAnglePCLHarvester") << "Could not retrieve the fit status! Setting it to 0 by default";
res.fitStatus = 0;
res.fitStatus = -1;
}

if (res.fitStatus != 0) {
if (res.fitStatus >= 0) {
res.covMatrixStatus = result->CovMatrixStatus();
} else {
res.covMatrixStatus = SiPixelLAHarvest::kUndefined;
Expand Down Expand Up @@ -767,14 +803,19 @@ SiPixelLAHarvest::fitResults SiPixelLorentzAnglePCLHarvester::fitAndStore(

float LorentzAnglePerTesla_;
float currentLA = currentLorentzAngle_->getLorentzAngle(detIdsToFill.front());
// if the fit quality is OK

// check the result of the chi2 only if doing the chebyshev fit
const bool chi2Cut = doChebyshevFit_ ? (res.redChi2 < fitChi2Cut_) : true;
// check the covariance matrix status
const bool covMatrixStatusCut = (res.covMatrixStatus >= SiPixelLAHarvest::kMadePosDef);
// fill the fit status dashboard
bool goodFit = isFitGood(res);
for (unsigned int i_bin = 0; i_bin < res.quality.size(); i_bin++) {
if (res.quality[i_bin]) {
hists_.h_bySectFitQuality_->setBinContent(i_index, i_bin + 1, 1.);
} else {
hists_.h_bySectFitQuality_->setBinContent(i_index, i_bin + 1, 0.);
}
}

if ((res.redChi2 != 0.) && covMatrixStatusCut && chi2Cut && (res.nentries > minHitsCut_)) {
// if the fit quality is OK
if (goodFit) {
LorentzAnglePerTesla_ = res.tan_LA / theMagField_;
// fill the LA actually written to payload
hists_.h_bySectSetLA_->setBinContent(i_index, LorentzAnglePerTesla_);
Expand Down Expand Up @@ -810,7 +851,7 @@ SiPixelLAHarvest::fitResults SiPixelLorentzAnglePCLHarvester::fitAndStore(
}

// function to check the fit quality
bool SiPixelLorentzAnglePCLHarvester::checkFitQuality(SiPixelLAHarvest::fitResults& res) {
bool SiPixelLorentzAnglePCLHarvester::isFitGood(SiPixelLAHarvest::fitResults& res) {
// check if reduced chi2 is different from 0.
const bool notZeroCut = (res.redChi2 != 0.);
// check the result of the chi2 only if doing the chebyshev fit
Expand All @@ -821,16 +862,16 @@ bool SiPixelLorentzAnglePCLHarvester::checkFitQuality(SiPixelLAHarvest::fitResul
const bool nentriesCut = (res.nentries > minHitsCut_);

if (notZeroCut) {
res.quality.set(0);
res.setQualityCutBit(SiPixelLAHarvest::kZeroChi2);
}
if (chi2Cut) {
res.quality.set(1);
res.setQualityCutBit(SiPixelLAHarvest::kChi2Cut);
}
if (covMatrixStatusCut) {
res.quality.set(2);
res.setQualityCutBit(SiPixelLAHarvest::kCovStatus);
}
if (nentriesCut) {
res.quality.set(3);
res.setQualityCutBit(SiPixelLAHarvest::kNentries);
}

// check if all the bits are set
Expand Down

0 comments on commit 216a222

Please sign in to comment.