Skip to content

Commit

Permalink
Merge branch 'cms-sw:master' into dev_20220819
Browse files Browse the repository at this point in the history
  • Loading branch information
@davidwalter2
davidwalter2 committed Aug 21, 2022
2 parents bd620d4 + 3b570db commit c1dcaf9
Show file tree
Hide file tree
Showing 12 changed files with 246 additions and 140 deletions.
7 changes: 5 additions & 2 deletions CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,8 @@ class PPSTimingCalibrationPCLHarvester : public DQMEDHarvester {
const std::string dqmDir_;
const std::string formula_;
const unsigned int min_entries_;
static constexpr double resolution_ = 0.1;
static constexpr double offset_ = 0.;
TF1 interp_;
};

Expand Down Expand Up @@ -96,7 +98,7 @@ void PPSTimingCalibrationPCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQM
(int)detid.arm(), (int)detid.station(), (int)detid.plane(), (int)detid.channel()};

calib_params[key] = {0, 0, 0, 0};
calib_time[key] = std::make_pair(0.1, 0.);
calib_time[key] = std::make_pair(offset_, resolution_);

hists.leadingTime[chid] = iGetter.get(dqmDir_ + "/t_" + ch_name);
if (hists.leadingTime[chid] == nullptr) {
Expand Down Expand Up @@ -134,7 +136,8 @@ void PPSTimingCalibrationPCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQM
if (res == 0) {
calib_params[key] = {
interp_.GetParameter(0), interp_.GetParameter(1), interp_.GetParameter(2), interp_.GetParameter(3)};
calib_time[key] = std::make_pair(0.1, 0.); // hardcoded resolution/offset placeholder for the time being
calib_time[key] =
std::make_pair(offset_, resolution_); // hardcoded offset/resolution placeholder for the time being
// can possibly do something with interp_.GetChiSquare() in the near future
} else
edm::LogWarning("PPSTimingCalibrationPCLHarvester:dqmEndJob")
Expand Down
3 changes: 3 additions & 0 deletions CalibTracker/SiPixelLorentzAngle/interface/SiPixelLorentzAngleCalibrationStruct.h
View file
Original file line number Diff line number Diff line change
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
112 changes: 95 additions & 17 deletions CalibTracker/SiPixelLorentzAngle/src/SiPixelLorentzAnglePCLHarvester.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,8 @@
*/
namespace SiPixelLAHarvest {

enum covStatus { kNotCalculated = 0, kApproximated = 1, kMadePosDef = 2, kAccurate = 3 };
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 @@ -57,9 +58,29 @@ namespace SiPixelLAHarvest {
e0 = e1 = e2 = e3 = e4 = e5 = 0.;
chi2 = prob = dSq = redChi2 = -9999.;
tan_LA = error_LA = -9999.;
fitStatus = covMatrixStatus = ndf = -999;
fitStatus = covMatrixStatus = ndf = nentries = -999;
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 All @@ -68,12 +89,14 @@ namespace SiPixelLAHarvest {
double p5, e5;
double chi2;
int ndf;
int nentries;
double prob;
int fitStatus, covMatrixStatus;
double dSq;
double redChi2;
double tan_LA;
double error_LA;
std::bitset<4> quality; /* to store if passes cuts*/
};
} // namespace SiPixelLAHarvest

Expand All @@ -91,6 +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 isFitGood(SiPixelLAHarvest::fitResults& res);

// es tokens
edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geomEsToken_;
Expand Down Expand Up @@ -365,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 @@ -376,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 @@ -388,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 @@ -629,15 +669,20 @@ SiPixelLAHarvest::fitResults SiPixelLorentzAnglePCLHarvester::fitAndStore(
if (doChebyshevFit_) {
const int npar = order_ + 1;
auto cheb = std::make_unique<siPixelLACalibration::Chebyshev>(order_, theFitRange_.first, theFitRange_.second);
f1_ =
std::make_unique<TF1>("fChebyshev", cheb.release(), theFitRange_.first, theFitRange_.second, npar, "Chebyshev");
f1_ = std::make_unique<TF1>("f1", cheb.release(), theFitRange_.first, theFitRange_.second, npar, "Chebyshev");
} else {
f1_ = std::make_unique<TF1>("fPolynomial",
f1_ = std::make_unique<TF1>("f1",
"[0] + [1]*x + [2]*x*x + [3]*x*x*x + [4]*x*x*x*x + [5]*x*x*x*x*x",
theFitRange_.first,
theFitRange_.second);
}

// 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(strcmp(f1_->GetName(), "f1") == 0);

f1_->SetParName(0, "offset");
f1_->SetParName(1, "tan#theta_{LA}");
f1_->SetParName(2, "quad term");
Expand Down Expand Up @@ -665,20 +710,20 @@ 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::kNotCalculated;
res.covMatrixStatus = SiPixelLAHarvest::kUndefined;
}

// compute the error on the drift-at-half-width parameter (d^2)
float dSq{0.};
float cov00{0.}; // needed later for the error on the tan(theta_LA)
if (!doChebyshevFit_) {
if (res.covMatrixStatus != SiPixelLAHarvest::kNotCalculated) {
if (res.covMatrixStatus > SiPixelLAHarvest::kNotCalculated) {
for (int k = 0; k < order_; k++) {
for (int l = 0; l < order_; l++) {
dSq += (std::pow(half_width, k) * std::pow(half_width, l) * result->CovMatrix(k, l));
Expand Down Expand Up @@ -733,7 +778,7 @@ SiPixelLAHarvest::fitResults SiPixelLorentzAnglePCLHarvester::fitAndStore(
hists_.h_bySectDriftError_->setBinContent(i_index, res.dSq);
hists_.h_bySectDriftError_->setBinError(i_index, 0.);

int nentries = hists_.h_bySectOccupancy_->getBinContent(i_index); // number of on track hits in that sector
res.nentries = hists_.h_bySectOccupancy_->getBinContent(i_index); // number of on track hits in that sector

bool isNew = (i_index > hists_.nlay * hists_.nModules_[hists_.nlay - 1]);
int shiftIdx = i_index - hists_.nlay * hists_.nModules_[hists_.nlay - 1] - 1;
Expand All @@ -758,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::kAccurate);
// 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 && (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 @@ -800,6 +850,34 @@ SiPixelLAHarvest::fitResults SiPixelLorentzAnglePCLHarvester::fitAndStore(
return res;
}

// function to check the fit quality
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
const bool chi2Cut = doChebyshevFit_ ? (res.redChi2 < fitChi2Cut_) : true;
// check the covariance matrix status
const bool covMatrixStatusCut = (res.covMatrixStatus >= SiPixelLAHarvest::kMadePosDef);
// check that the number of entries is larger than the minimum amount
const bool nentriesCut = (res.nentries > minHitsCut_);

if (notZeroCut) {
res.setQualityCutBit(SiPixelLAHarvest::kZeroChi2);
}
if (chi2Cut) {
res.setQualityCutBit(SiPixelLAHarvest::kChi2Cut);
}
if (covMatrixStatusCut) {
res.setQualityCutBit(SiPixelLAHarvest::kCovStatus);
}
if (nentriesCut) {
res.setQualityCutBit(SiPixelLAHarvest::kNentries);
}

// check if all the bits are set
return (res.quality.all());
}

//------------------------------------------------------------------------------
void SiPixelLorentzAnglePCLHarvester::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
edm::ParameterSetDescription desc;
Expand Down
54 changes: 28 additions & 26 deletions Calibration/HcalAlCaRecoProducers/plugins/AlCaIsolatedBunchSelector.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,25 +27,25 @@

//#define EDM_ML_DEBUG

namespace AlCaIsolatedBunch {
namespace alCaIsolatedBunchSelector {
struct Counters {
Counters() : nAll_(0), nGood_(0) {}
mutable std::atomic<unsigned int> nAll_, nGood_;
};
} // namespace AlCaIsolatedBunch
} // namespace alCaIsolatedBunchSelector

class AlCaIsolatedBunchSelector : public edm::stream::EDFilter<edm::GlobalCache<AlCaIsolatedBunch::Counters> > {
class AlCaIsolatedBunchSelector : public edm::stream::EDFilter<edm::GlobalCache<alCaIsolatedBunchSelector::Counters> > {
public:
explicit AlCaIsolatedBunchSelector(edm::ParameterSet const&, const AlCaIsolatedBunch::Counters* count);
~AlCaIsolatedBunchSelector() override;
explicit AlCaIsolatedBunchSelector(edm::ParameterSet const&, const alCaIsolatedBunchSelector::Counters* count);
~AlCaIsolatedBunchSelector() override = default;

static std::unique_ptr<AlCaIsolatedBunch::Counters> initializeGlobalCache(edm::ParameterSet const& iConfig) {
return std::make_unique<AlCaIsolatedBunch::Counters>();
static std::unique_ptr<alCaIsolatedBunchSelector::Counters> initializeGlobalCache(edm::ParameterSet const& iConfig) {
return std::make_unique<alCaIsolatedBunchSelector::Counters>();
}

bool filter(edm::Event&, edm::EventSetup const&) override;
void endStream() override;
static void globalEndJob(const AlCaIsolatedBunch::Counters* counters);
static void globalEndJob(const alCaIsolatedBunchSelector::Counters* counters);
static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);

private:
Expand All @@ -65,7 +65,7 @@ class AlCaIsolatedBunchSelector : public edm::stream::EDFilter<edm::GlobalCache<
// constructors and destructor
//
AlCaIsolatedBunchSelector::AlCaIsolatedBunchSelector(const edm::ParameterSet& iConfig,
const AlCaIsolatedBunch::Counters* count)
const alCaIsolatedBunchSelector::Counters* count)
: nRun_(0),
nAll_(0),
nGood_(0),
Expand All @@ -78,8 +78,6 @@ AlCaIsolatedBunchSelector::AlCaIsolatedBunchSelector(const edm::ParameterSet& iC
<< std::endl;
}

AlCaIsolatedBunchSelector::~AlCaIsolatedBunchSelector() {}

//
// member functions
//
Expand All @@ -92,21 +90,25 @@ bool AlCaIsolatedBunchSelector::filter(edm::Event& iEvent, edm::EventSetup const
edm::LogVerbatim("AlCaIsoBunch") << "Run " << iEvent.id().run() << " Event " << iEvent.id().event() << " Luminosity "
<< iEvent.luminosityBlock() << " Bunch " << iEvent.bunchCrossing() << std::endl;
#endif
//Step1: Find if the event passes the chosen trigger
auto const& triggerResults = iEvent.getHandle(tok_trigRes_);
if (triggerResults.isValid()) {
const edm::TriggerNames& triggerNames = iEvent.triggerNames(*triggerResults);
const std::vector<std::string>& triggerNames_ = triggerNames.triggerNames();
for (unsigned int iHLT = 0; iHLT < triggerResults->size(); iHLT++) {
int hlt = triggerResults->accept(iHLT);
if (triggerNames_[iHLT].find(trigName_) != std::string::npos) {
if (hlt > 0) {
accept = true;
if (trigName_.empty()) {
accept = true;
} else {
//Step1: Find if the event passes the chosen trigger
auto const& triggerResults = iEvent.getHandle(tok_trigRes_);
if (triggerResults.isValid()) {
const edm::TriggerNames& triggerNames = iEvent.triggerNames(*triggerResults);
const std::vector<std::string>& triggerNames_ = triggerNames.triggerNames();
for (unsigned int iHLT = 0; iHLT < triggerResults->size(); iHLT++) {
int hlt = triggerResults->accept(iHLT);
if (triggerNames_[iHLT].find(trigName_) != std::string::npos) {
if (hlt > 0) {
accept = true;
#ifdef EDM_ML_DEBUG
edm::LogVerbatim("AlCaIsoBunch")
<< triggerNames_[iHLT] << " has got HLT flag " << hlt << ":" << accept << std::endl;
edm::LogVerbatim("AlCaIsoBunch")
<< triggerNames_[iHLT] << " has got HLT flag " << hlt << ":" << accept << std::endl;
#endif
break;
break;
}
}
}
}
Expand All @@ -124,7 +126,7 @@ void AlCaIsolatedBunchSelector::endStream() {
globalCache()->nGood_ += nGood_;
}

void AlCaIsolatedBunchSelector::globalEndJob(const AlCaIsolatedBunch::Counters* count) {
void AlCaIsolatedBunchSelector::globalEndJob(const alCaIsolatedBunchSelector::Counters* count) {
edm::LogVerbatim("AlCaIsoBunch") << "Selects " << count->nGood_ << " in " << count->nAll_ << " events" << std::endl;
}

Expand All @@ -145,7 +147,7 @@ void AlCaIsolatedBunchSelector::fillDescriptions(edm::ConfigurationDescriptions&
edm::ParameterSetDescription desc;
desc.add<edm::InputTag>("triggerResultLabel", edm::InputTag("TriggerResults", "", "HLT"));
desc.add<std::string>("processName", "HLT");
desc.add<std::string>("triggerName", "HLT_HcalIsolatedBunch");
desc.add<std::string>("triggerName", "");
descriptions.add("alcaIsolatedBunchSelector", desc);
}

Expand Down

0 comments on commit c1dcaf9

Please sign in to comment.