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

Introduce missing hits recovery in the SiStripHitEfficiency PCL workflow to be in synch with CalibTree based one #42175

Merged
Changes from 1 commit
Commits
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
Original file line number Diff line number Diff line change
@@ -401,8 +401,9 @@ void SiStripHitEfficiencyHarvester::dqmEndJob(DQMStore::IBooker& booker, DQMStor
const auto num = h_module_found->getValue(det);
const auto denom = h_module_total->getValue(det);
if (denom) {
assert(num <= denom); // can't have this happen
const auto eff = num / denom;
const auto eff_up = TEfficiency::Bayesian(denom, num, .99, 1, 1, true);
const auto eff_up = (num <= denom) ? TEfficiency::Bayesian(denom, num, .99, 1, 1, true) : 1.f;
mmusich marked this conversation as resolved.
Show resolved Hide resolved
mmusich marked this conversation as resolved.
Show resolved Hide resolved

if ((denom >= nModsMin_) && (eff_up < layer_min_eff)) {
//We have a bad module, put it in the list!
@@ -465,7 +466,7 @@ void SiStripHitEfficiencyHarvester::dqmEndJob(DQMStore::IBooker& booker, DQMStor
LOGPRINT << "Number of strips module " << det << " is " << nStrips;
badStripList.push_back(pQuality.encode(0, nStrips, 0));
//Now compact into a single bad module
LOGPRINT << "ID1 shoudl match list of modules above " << det;
LOGPRINT << "ID1 should match list of modules above " << det;
pQuality.compact(det, badStripList);
pQuality.put(det, SiStripQuality::Range(badStripList.begin(), badStripList.end()));
}
Original file line number Diff line number Diff line change
@@ -117,6 +117,7 @@ class SiStripHitEfficiencyWorker : public DQMEDAnalyzer {
bool useFirstMeas_;
bool useLastMeas_;
bool useAllHitsFromTracksWithMissingHits_;
bool doMissingHitsRecovery_;
unsigned int clusterMatchingMethod_;
float resXSig_;
float clusterTracjDist_;
@@ -130,6 +131,12 @@ class SiStripHitEfficiencyWorker : public DQMEDAnalyzer {
// output file
std::set<uint32_t> badModules_;

// for the missing hits recovery
std::vector<unsigned int> hitRecoveryCounters;
std::vector<unsigned int> hitTotalCounters;
int totalNbHits;
std::vector<int> missHitPerLayer;

struct EffME1 {
EffME1() : hTotal(nullptr), hFound(nullptr) {}
EffME1(MonitorElement* total, MonitorElement* found) : hTotal(total), hFound(found) {}
@@ -155,6 +162,14 @@ class SiStripHitEfficiencyWorker : public DQMEDAnalyzer {
}
}

bool check(uint32_t id) {
if (hTotal->getValue(id) < hFound->getValue(id)) {
return false;
} else {
return true;
}
}

std::unique_ptr<TkHistoMap> hTotal, hFound;
};

@@ -208,6 +223,7 @@ SiStripHitEfficiencyWorker::SiStripHitEfficiencyWorker(const edm::ParameterSet&
useFirstMeas_(conf.getParameter<bool>("useFirstMeas")),
useLastMeas_(conf.getParameter<bool>("useLastMeas")),
useAllHitsFromTracksWithMissingHits_(conf.getParameter<bool>("useAllHitsFromTracksWithMissingHits")),
doMissingHitsRecovery_(conf.getParameter<bool>("doMissingHitsRecovery")),
clusterMatchingMethod_(conf.getParameter<int>("ClusterMatchingMethod")),
resXSig_(conf.getParameter<double>("ResXSig")),
clusterTracjDist_(conf.getParameter<double>("ClusterTrajDist")),
@@ -216,6 +232,11 @@ SiStripHitEfficiencyWorker::SiStripHitEfficiencyWorker(const edm::ParameterSet&
bunchX_(conf.getUntrackedParameter<int>("BunchCrossing", 0)),
showRings_(conf.getUntrackedParameter<bool>("ShowRings", false)),
showTOB6TEC9_(conf.getUntrackedParameter<bool>("ShowTOB6TEC9", false)) {
hitRecoveryCounters.resize(k_END_OF_LAYERS, 0);
hitTotalCounters.resize(k_END_OF_LAYERS, 0);
missHitPerLayer.resize(k_END_OF_LAYERS, 0);
totalNbHits = 0;

nTEClayers_ = (showRings_ ? 7 : 9); // number of rings or wheels

const std::string badModulesFile = conf.getUntrackedParameter<std::string>("BadModulesFile", "");
@@ -477,13 +498,65 @@ void SiStripHitEfficiencyWorker::analyze(const edm::Event& e, const edm::EventSe

const bool highPurity = trajTrack.val->quality(reco::TrackBase::TrackQuality::highPurity);
auto TMeas = trajTrack.key->measurements();
totalNbHits += int(TMeas.size());

/*
const bool hasMissingHits = std::any_of(std::begin(TMeas), std::end(TMeas), [](const auto& tm) {
return tm.recHit()->getType() == TrackingRecHit::Type::missing;
});
*/

// Check whether the trajectory has some missing hits
bool hasMissingHits{false};
int previous_layer{999};
std::vector<unsigned int> missedLayers;

for (const auto& itm : TMeas) {
auto theHit = itm.recHit();
unsigned int iidd = theHit->geographicalId().rawId();
int layer = ::checkLayer(iidd, tTopo);
int missedLayer = layer + 1;
int diffPreviousLayer = (layer - previous_layer);
if (doMissingHitsRecovery_) {
//Layers from TIB + TOB
if (diffPreviousLayer == -2 && missedLayer > k_LayersStart && missedLayer < k_LayersAtTOBEnd) {
missHitPerLayer[missedLayer] += 1;
hasMissingHits = true;
}
//Layers from TID
else if (diffPreviousLayer == -2 && (missedLayer > k_LayersAtTOBEnd + 1 && missedLayer <= k_LayersAtTIDEnd)) {
missHitPerLayer[missedLayer] += 1;
hasMissingHits = true;
}
//Layers from TEC
else if (diffPreviousLayer == -2 && missedLayer > k_LayersAtTIDEnd && missedLayer <= k_LayersAtTECEnd) {
missHitPerLayer[missedLayer] += 1;
hasMissingHits = true;
}

//##### TID Layer 11 in transition TID -> TIB : layer is in TIB, previous layer = 12
if ((layer > k_LayersStart && layer <= k_LayersAtTIBEnd) && (previous_layer == 12)) {
missHitPerLayer[11] += 1;
hasMissingHits = true;
}

//##### TEC Layer 14 in transition TEC -> TOB : layer is in TOB, previous layer = 15
if ((layer > k_LayersAtTIBEnd && layer <= k_LayersAtTOBEnd) && (previous_layer == 15)) {
missHitPerLayer[14] += 1;
hasMissingHits = true;
}
}
if (theHit->getType() == TrackingRecHit::Type::missing)
hasMissingHits = true;

if (hasMissingHits)
missedLayers.push_back(layer);
previous_layer = layer;
}

// Loop on each measurement and take into consideration
//--------------------------------------------------------
unsigned int prev_TKlayers = 0;
for (auto itm = TMeas.cbegin(); itm != TMeas.cend(); ++itm) {
const auto theInHit = (*itm).recHit();

@@ -506,11 +579,19 @@ void SiStripHitEfficiencyWorker::analyze(const edm::Event& e, const edm::EventSe

// Test first and last points of the trajectory
// the list of measurements starts from outer layers !!! This could change -> should add a check
if ((!useFirstMeas_ && (itm == (TMeas.end() - 1))) || (!useLastMeas_ && (itm == (TMeas.begin()))) ||
// In case of missing hit in the track, check whether to use the other hits or not.
(!useAllHitsFromTracksWithMissingHits_ && hasMissingHits &&
theInHit->getType() != TrackingRecHit::Type::missing))
bool isFirstMeas = (itm == (TMeas.end() - 1));
bool isLastMeas = (itm == (TMeas.begin()));

if (!useFirstMeas_ && isFirstMeas)
continue;
if (!useLastMeas_ && isLastMeas)
continue;

// In case of missing hit in the track, check whether to use the other hits or not.
if (hasMissingHits && theInHit->getType() != TrackingRecHit::Type::missing &&
!useAllHitsFromTracksWithMissingHits_)
continue;

// If Trajectory measurement from TOB 6 or TEC 9, skip it because it's always valid they are filled later
if (TKlayers == bounds::k_LayersAtTOBEnd || TKlayers == bounds::k_LayersAtTECEnd) {
LogDebug("SiStripHitEfficiencyWorker") << "skipping original TM for TOB 6 or TEC 9";
@@ -541,6 +622,125 @@ void SiStripHitEfficiencyWorker::analyze(const edm::Event& e, const edm::EventSe
TMs.emplace_back(*itm, tTopo, tkgeom, propagator);
}

bool missingHitAdded{false};
std::vector<TrajectoryMeasurement> tmpTmeas;
unsigned int misLayer = TKlayers + 1;
//Use bool doMissingHitsRecovery to add possible missing hits based on actual/previous hit
if (doMissingHitsRecovery_) {
if (int(TKlayers) - int(prev_TKlayers) == -2) {
const DetLayer* detlayer = itm->layer();
const LayerMeasurements layerMeasurements{measTracker, *measurementTrackerEvent};
const TrajectoryStateOnSurface tsos = itm->updatedState();
std::vector<DetLayer::DetWithState> compatDets = detlayer->compatibleDets(tsos, prop, chi2Estimator);

if (misLayer > k_LayersAtTIDEnd && misLayer < k_LayersAtTECEnd) { //TEC
std::vector<ForwardDetLayer const*> negTECLayers = measTracker.geometricSearchTracker()->negTecLayers();
std::vector<ForwardDetLayer const*> posTECLayers = measTracker.geometricSearchTracker()->posTecLayers();
const DetLayer* tecLayerneg = negTECLayers[misLayer - k_LayersAtTIDEnd - 1];
const DetLayer* tecLayerpos = posTECLayers[misLayer - k_LayersAtTIDEnd - 1];
if (tTopo->tecSide(iidd) == 1) {
tmpTmeas = layerMeasurements.measurements(*tecLayerneg, tsos, prop, chi2Estimator);
} else if (tTopo->tecSide(iidd) == 2) {
tmpTmeas = layerMeasurements.measurements(*tecLayerpos, tsos, prop, chi2Estimator);
}
}

else if (misLayer == (k_LayersAtTIDEnd - 1) ||
misLayer == k_LayersAtTIDEnd) { // This is for TID layers 12 and 13
std::vector<ForwardDetLayer const*> negTIDLayers = measTracker.geometricSearchTracker()->negTidLayers();
std::vector<ForwardDetLayer const*> posTIDLayers = measTracker.geometricSearchTracker()->posTidLayers();
const DetLayer* tidLayerneg = negTIDLayers[misLayer - k_LayersAtTOBEnd - 1];
const DetLayer* tidLayerpos = posTIDLayers[misLayer - k_LayersAtTOBEnd - 1];

if (tTopo->tidSide(iidd) == 1) {
tmpTmeas = layerMeasurements.measurements(*tidLayerneg, tsos, prop, chi2Estimator);
} else if (tTopo->tidSide(iidd) == 2) {
tmpTmeas = layerMeasurements.measurements(*tidLayerpos, tsos, prop, chi2Estimator);
}
}

if (misLayer > k_LayersStart && misLayer < k_LayersAtTOBEnd) { // Barrel

std::vector<BarrelDetLayer const*> barrelTIBLayers = measTracker.geometricSearchTracker()->tibLayers();
std::vector<BarrelDetLayer const*> barrelTOBLayers = measTracker.geometricSearchTracker()->tobLayers();

if (misLayer > k_LayersStart && misLayer <= k_LayersAtTIBEnd) {
const DetLayer* tibLayer = barrelTIBLayers[misLayer - k_LayersStart - 1];
tmpTmeas = layerMeasurements.measurements(*tibLayer, tsos, prop, chi2Estimator);
} else if (misLayer > k_LayersAtTIBEnd && misLayer < k_LayersAtTOBEnd) {
const DetLayer* tobLayer = barrelTOBLayers[misLayer - k_LayersAtTIBEnd - 1];
tmpTmeas = layerMeasurements.measurements(*tobLayer, tsos, prop, chi2Estimator);
}
}
}
if ((int(TKlayers) > k_LayersStart && int(TKlayers) <= k_LayersAtTIBEnd) && int(prev_TKlayers) == 12) {
const DetLayer* detlayer = itm->layer();
const LayerMeasurements layerMeasurements{measTracker, *measurementTrackerEvent};
const TrajectoryStateOnSurface tsos = itm->updatedState();
std::vector<DetLayer::DetWithState> compatDets = detlayer->compatibleDets(tsos, prop, chi2Estimator);
std::vector<ForwardDetLayer const*> negTIDLayers = measTracker.geometricSearchTracker()->negTidLayers();
std::vector<ForwardDetLayer const*> posTIDLayers = measTracker.geometricSearchTracker()->posTidLayers();

const DetLayer* tidLayerneg = negTIDLayers[k_LayersStart];
const DetLayer* tidLayerpos = posTIDLayers[k_LayersStart];
if (tTopo->tidSide(iidd) == 1) {
tmpTmeas = layerMeasurements.measurements(*tidLayerneg, tsos, prop, chi2Estimator);
} else if (tTopo->tidSide(iidd) == 2) {
tmpTmeas = layerMeasurements.measurements(*tidLayerpos, tsos, prop, chi2Estimator);
}
}

if ((int(TKlayers) > k_LayersAtTIBEnd && int(TKlayers) <= k_LayersAtTOBEnd) && int(prev_TKlayers) == 15) {
const DetLayer* detlayer = itm->layer();
const LayerMeasurements layerMeasurements{measTracker, *measurementTrackerEvent};
const TrajectoryStateOnSurface tsos = itm->updatedState();
std::vector<DetLayer::DetWithState> compatDets = detlayer->compatibleDets(tsos, prop, chi2Estimator);

std::vector<ForwardDetLayer const*> negTECLayers = measTracker.geometricSearchTracker()->negTecLayers();
std::vector<ForwardDetLayer const*> posTECLayers = measTracker.geometricSearchTracker()->posTecLayers();

const DetLayer* tecLayerneg = negTECLayers[k_LayersStart];
const DetLayer* tecLayerpos = posTECLayers[k_LayersStart];
if (tTopo->tecSide(iidd) == 1) {
tmpTmeas = layerMeasurements.measurements(*tecLayerneg, tsos, prop, chi2Estimator);
} else if (tTopo->tecSide(iidd) == 2) {
tmpTmeas = layerMeasurements.measurements(*tecLayerpos, tsos, prop, chi2Estimator);
}
}

if (!tmpTmeas.empty()) {
TrajectoryMeasurement TM_tmp(tmpTmeas.back());
unsigned int iidd_tmp = TM_tmp.recHit()->geographicalId().rawId();
if (iidd_tmp != 0) {
LogDebug("SiStripHitEfficiency:HitEff") << " hit actually being added to TM vector";
if ((!useAllHitsFromTracksWithMissingHits_ || (!useFirstMeas_ && isFirstMeas)))
TMs.clear();
if (::isDoubleSided(iidd_tmp, tTopo)) {
TMs.push_back(TrajectoryAtInvalidHit(TM_tmp, tTopo, tkgeom, propagator, 1));
TMs.push_back(TrajectoryAtInvalidHit(TM_tmp, tTopo, tkgeom, propagator, 2));
} else
TMs.push_back(TrajectoryAtInvalidHit(TM_tmp, tTopo, tkgeom, propagator));
missingHitAdded = true;
hitRecoveryCounters[misLayer] += 1;
}
}
}

prev_TKlayers = TKlayers;
if (!useFirstMeas_ && isFirstMeas && !missingHitAdded)
continue;
if (!useLastMeas_ && isLastMeas)
continue;
bool hitsWithBias = false;
for (auto ilayer : missedLayers) {
if (ilayer < TKlayers)
hitsWithBias = true;
}
if (hasMissingHits && theInHit->getType() != TrackingRecHit::Type::missing && !missingHitAdded &&
hitsWithBias && !useAllHitsFromTracksWithMissingHits_) {
continue;
}

//////////////////////////////////////////////
//Now check for tracks at TOB6 and TEC9

@@ -987,6 +1187,7 @@ void SiStripHitEfficiencyWorker::fillForTraj(const TrajectoryAtInvalidHit& tm,
// efficiency without bad modules excluded
if (TKlayers) {
h_module.fill(iidd, !badflag);
assert(h_module.check(iidd));
}

/* Used in SiStripHitEffFromCalibTree:
@@ -1023,6 +1224,7 @@ void SiStripHitEfficiencyWorker::fillDescriptions(edm::ConfigurationDescriptions
desc.add<std::string>("dqmDir", "AlCaReco/SiStripHitEfficiency");
desc.add<bool>("UseOnlyHighPurityTracks", true);
desc.add<bool>("cutOnTracks", false);
desc.add<bool>("doMissingHitsRecovery", false);
desc.add<bool>("useAllHitsFromTracksWithMissingHits", false);
desc.add<bool>("useFirstMeas", false);
desc.add<bool>("useLastMeas", false);
Original file line number Diff line number Diff line change
@@ -67,11 +67,17 @@
useFirstMeas = False,
useLastMeas = False,
useAllHitsFromTracksWithMissingHits = False,
doMissingHitsRecovery = False,
## non-default settings
ClusterMatchingMethod = 4, # default 0 case0,1,2,3,4
ClusterTrajDist = 15, # default 64
)

from Configuration.Eras.Modifier_run3_common_cff import run3_common
run3_common.toModify(ALCARECOSiStripHitEff,
useAllHitsFromTracksWithMissingHits = True,
doMissingHitsRecovery = True)

# ----------------------------------------------------------------------------
MEtoEDMConvertSiStripHitEff = cms.EDProducer("MEtoEDMConverter",
Name = cms.untracked.string('MEtoEDMConverter'),