diff --git a/CalibTracker/SiStripHitEfficiency/plugins/SiStripHitEfficiencyHarvester.cc b/CalibTracker/SiStripHitEfficiency/plugins/SiStripHitEfficiencyHarvester.cc index 6758675c647ea..0e8c43a1b8951 100644 --- a/CalibTracker/SiStripHitEfficiency/plugins/SiStripHitEfficiencyHarvester.cc +++ b/CalibTracker/SiStripHitEfficiency/plugins/SiStripHitEfficiencyHarvester.cc @@ -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! @@ -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())); } diff --git a/CalibTracker/SiStripHitEfficiency/plugins/SiStripHitEfficiencyWorker.cc b/CalibTracker/SiStripHitEfficiency/plugins/SiStripHitEfficiencyWorker.cc index a99b4f47f9411..8149c4a8275a4 100644 --- a/CalibTracker/SiStripHitEfficiency/plugins/SiStripHitEfficiencyWorker.cc +++ b/CalibTracker/SiStripHitEfficiency/plugins/SiStripHitEfficiencyWorker.cc @@ -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 badModules_; + // for the missing hits recovery + std::vector hitRecoveryCounters; + std::vector hitTotalCounters; + int totalNbHits; + std::vector 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 hTotal, hFound; }; @@ -208,6 +223,7 @@ SiStripHitEfficiencyWorker::SiStripHitEfficiencyWorker(const edm::ParameterSet& useFirstMeas_(conf.getParameter("useFirstMeas")), useLastMeas_(conf.getParameter("useLastMeas")), useAllHitsFromTracksWithMissingHits_(conf.getParameter("useAllHitsFromTracksWithMissingHits")), + doMissingHitsRecovery_(conf.getParameter("doMissingHitsRecovery")), clusterMatchingMethod_(conf.getParameter("ClusterMatchingMethod")), resXSig_(conf.getParameter("ResXSig")), clusterTracjDist_(conf.getParameter("ClusterTrajDist")), @@ -216,6 +232,11 @@ SiStripHitEfficiencyWorker::SiStripHitEfficiencyWorker(const edm::ParameterSet& bunchX_(conf.getUntrackedParameter("BunchCrossing", 0)), showRings_(conf.getUntrackedParameter("ShowRings", false)), showTOB6TEC9_(conf.getUntrackedParameter("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("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 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 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 compatDets = detlayer->compatibleDets(tsos, prop, chi2Estimator); + + if (misLayer > k_LayersAtTIDEnd && misLayer < k_LayersAtTECEnd) { //TEC + std::vector negTECLayers = measTracker.geometricSearchTracker()->negTecLayers(); + std::vector 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 negTIDLayers = measTracker.geometricSearchTracker()->negTidLayers(); + std::vector 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 barrelTIBLayers = measTracker.geometricSearchTracker()->tibLayers(); + std::vector 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 compatDets = detlayer->compatibleDets(tsos, prop, chi2Estimator); + std::vector negTIDLayers = measTracker.geometricSearchTracker()->negTidLayers(); + std::vector 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 compatDets = detlayer->compatibleDets(tsos, prop, chi2Estimator); + + std::vector negTECLayers = measTracker.geometricSearchTracker()->negTecLayers(); + std::vector 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("dqmDir", "AlCaReco/SiStripHitEfficiency"); desc.add("UseOnlyHighPurityTracks", true); desc.add("cutOnTracks", false); + desc.add("doMissingHitsRecovery", false); desc.add("useAllHitsFromTracksWithMissingHits", false); desc.add("useFirstMeas", false); desc.add("useLastMeas", false); diff --git a/Calibration/TkAlCaRecoProducers/python/ALCARECOPromptCalibProdSiStripHitEfficiency_cff.py b/Calibration/TkAlCaRecoProducers/python/ALCARECOPromptCalibProdSiStripHitEfficiency_cff.py index 8c515b1c1e400..d3af2bae4959d 100644 --- a/Calibration/TkAlCaRecoProducers/python/ALCARECOPromptCalibProdSiStripHitEfficiency_cff.py +++ b/Calibration/TkAlCaRecoProducers/python/ALCARECOPromptCalibProdSiStripHitEfficiency_cff.py @@ -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'),