Skip to content

Commit

Permalink
introduce missing hits recovery in the SiStripHitEfficiency PCL workf…
Browse files Browse the repository at this point in the history 
…low to be in synch with CalibTree based one
  • Loading branch information
@mmusich
mmusich committed Jul 3, 2023
1 parent 248f510 commit c0ec2d5
Show file tree
Hide file tree
Showing 3 changed files with 215 additions and 6 deletions.
5 changes: 3 additions & 2 deletions CalibTracker/SiStripHitEfficiency/plugins/SiStripHitEfficiencyHarvester.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;

if ((denom >= nModsMin_) && (eff_up < layer_min_eff)) {
//We have a bad module, put it in the list!
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 c0ec2d5

Please sign in to comment.