Skip to content

Commit

Permalink
Merge pull request #42175 from mmusich/missingHitsRecoveryInSiStripHi…
Browse files Browse the repository at this point in the history 
…tEfficiencyWorkflow

Introduce missing hits recovery in the `SiStripHitEfficiency` PCL workflow to be in synch with `CalibTree` based one
  • Loading branch information
@cmsbuild
cmsbuild authored Jul 5, 2023
2 parents 6db9ecd + a30235c commit 4172a00
Show file tree
Hide file tree
Showing 3 changed files with 214 additions and 5 deletions.
3 changes: 2 additions & 1 deletion CalibTracker/SiStripHitEfficiency/plugins/SiStripHitEfficiencyHarvester.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -401,6 +401,7 @@ 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);

Expand Down Expand Up @@ -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()));
}
Expand Down
210 changes: 206 additions & 4 deletions CalibTracker/SiStripHitEfficiency/plugins/SiStripHitEfficiencyWorker.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -117,6 +117,7 @@ class SiStripHitEfficiencyWorker : public DQMEDAnalyzer {
bool useFirstMeas_;
bool useLastMeas_;
bool useAllHitsFromTracksWithMissingHits_;
bool doMissingHitsRecovery_;
unsigned int clusterMatchingMethod_;
float resXSig_;
float clusterTracjDist_;
Expand All @@ -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) {}
Expand All @@ -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;
};

Expand Down Expand Up @@ -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")),
Expand All @@ -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", "");
Expand Down Expand Up @@ -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();

Expand All @@ -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";
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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);
Expand Down
6 changes: 6 additions & 0 deletions Calibration/TkAlCaRecoProducers/python/ALCARECOPromptCalibProdSiStripHitEfficiency_cff.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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'),
Expand Down

0 comments on commit 4172a00

Please sign in to comment.