408 LogDebug(
"SiStripHitEfficiencyWorker") <<
"beginning analyze from HitEff";
424 }
else if (metaData.
isValid()) {
428 edm::LogWarning(
"SiStripHitEfficiencyWorker") <<
"could not find a source for the Luminosity and PU";
460 if (not fedErrorIdsCol_h.isValid() and not fedErrorIdsVec_h.isValid()) {
462 <<
"no valid product for SiStrip DetIds with FED errors (see parameter \"siStripDigis\"), " 463 "neither for new format (DetIdVector) nor old format (DetIdCollection)";
465 auto const& fedErrorIds = fedErrorIdsVec_h.isValid() ? *fedErrorIdsVec_h : fedErrorIdsCol_h->as_vector();
468 for (
const auto& fedErr : fedErrorIds) {
492 LogDebug(
"SiStripHitEfficiencyWorker") <<
"number ckf tracks found = " << tracksCKF->size();
502 if (!tracksCKF->empty()) {
506 LogDebug(
"SiStripHitEfficiencyWorker")
512 for (
const auto& trajTrack : *trajTrackAssociationHandle) {
516 auto TMeas = trajTrack.key->measurements();
526 bool hasMissingHits{
false};
527 int previous_layer{999};
528 std::vector<unsigned int> missedLayers;
530 for (
const auto& itm : TMeas) {
531 auto theHit = itm.recHit();
532 unsigned int iidd = theHit->geographicalId().rawId();
533 int layer = ::checkLayer(iidd, tTopo);
534 int missedLayer =
layer + 1;
535 int diffPreviousLayer = (
layer - previous_layer);
538 if (diffPreviousLayer == -2 && missedLayer > k_LayersStart && missedLayer < k_LayersAtTOBEnd) {
540 hasMissingHits =
true;
543 else if (diffPreviousLayer == -2 && (missedLayer > k_LayersAtTOBEnd + 1 && missedLayer <= k_LayersAtTIDEnd)) {
545 hasMissingHits =
true;
548 else if (diffPreviousLayer == -2 && missedLayer > k_LayersAtTIDEnd && missedLayer <= k_LayersAtTECEnd) {
550 hasMissingHits =
true;
554 if ((
layer > k_LayersStart &&
layer <= k_LayersAtTIBEnd) && (previous_layer == 12)) {
556 hasMissingHits =
true;
560 if ((
layer > k_LayersAtTIBEnd &&
layer <= k_LayersAtTOBEnd) && (previous_layer == 15)) {
562 hasMissingHits =
true;
566 hasMissingHits =
true;
569 missedLayers.push_back(
layer);
570 previous_layer =
layer;
575 unsigned int prev_TKlayers = 0;
576 for (
auto itm = TMeas.cbegin(); itm != TMeas.cend(); ++itm) {
577 const auto theInHit = (*itm).recHit();
582 LogDebug(
"SiStripHitEfficiencyWorker") <<
"theInHit is valid = " << theInHit->isValid();
584 unsigned int iidd = theInHit->geographicalId().rawId();
586 unsigned int TKlayers = ::checkLayer(iidd, tTopo);
592 LogDebug(
"SiStripHitEfficiencyWorker") <<
"TKlayer from trajectory: " << TKlayers <<
" from module = " << iidd
593 <<
" matched/stereo/rphi = " << ((iidd & 0x3) == 0) <<
"/" 594 << ((iidd & 0x3) == 1) <<
"/" << ((iidd & 0x3) == 2);
598 bool isFirstMeas = (itm == (TMeas.end() - 1));
599 bool isLastMeas = (itm == (TMeas.begin()));
612 if (TKlayers == bounds::k_LayersAtTOBEnd || TKlayers == bounds::k_LayersAtTECEnd) {
613 LogDebug(
"SiStripHitEfficiencyWorker") <<
"skipping original TM for TOB 6 or TEC 9";
617 std::vector<TrajectoryAtInvalidHit> TMs;
625 if (::isDoubleSided(iidd, tTopo) && ((iidd & 0x3) == 0)) {
628 TMs.emplace_back(*itm, tTopo, tkgeom,
propagator, 1);
629 TMs.emplace_back(*itm, tTopo, tkgeom,
propagator, 2);
630 }
else if (::isDoubleSided(iidd, tTopo) && (!::check2DPartner(iidd, TMeas))) {
633 TMs.emplace_back(*itm, tTopo, tkgeom,
propagator, 1);
634 TMs.emplace_back(*itm, tTopo, tkgeom,
propagator, 2);
635 LogDebug(
"SiStripHitEfficiencyWorker") <<
" found a hit with a missing partner";
638 TMs.emplace_back(*itm, tTopo, tkgeom,
propagator);
641 bool missingHitAdded{
false};
642 std::vector<TrajectoryMeasurement> tmpTmeas;
643 unsigned int misLayer = TKlayers + 1;
646 if (
int(TKlayers) -
int(prev_TKlayers) == -2) {
647 const DetLayer* detlayer = itm->layer();
650 std::vector<DetLayer::DetWithState> compatDets = detlayer->
compatibleDets(tsos, prop, chi2Estimator);
652 if (misLayer > k_LayersAtTIDEnd && misLayer < k_LayersAtTECEnd) {
653 std::vector<ForwardDetLayer const*> negTECLayers = measTracker.geometricSearchTracker()->negTecLayers();
654 std::vector<ForwardDetLayer const*> posTECLayers = measTracker.geometricSearchTracker()->posTecLayers();
655 const DetLayer* tecLayerneg = negTECLayers[misLayer - k_LayersAtTIDEnd - 1];
656 const DetLayer* tecLayerpos = posTECLayers[misLayer - k_LayersAtTIDEnd - 1];
657 if (tTopo->tecSide(iidd) == 1) {
658 tmpTmeas = layerMeasurements.measurements(*tecLayerneg, tsos, prop, chi2Estimator);
659 }
else if (tTopo->tecSide(iidd) == 2) {
660 tmpTmeas = layerMeasurements.measurements(*tecLayerpos, tsos, prop, chi2Estimator);
664 else if (misLayer == (k_LayersAtTIDEnd - 1) ||
665 misLayer == k_LayersAtTIDEnd) {
666 std::vector<ForwardDetLayer const*> negTIDLayers = measTracker.geometricSearchTracker()->negTidLayers();
667 std::vector<ForwardDetLayer const*> posTIDLayers = measTracker.geometricSearchTracker()->posTidLayers();
668 const DetLayer* tidLayerneg = negTIDLayers[misLayer - k_LayersAtTOBEnd - 1];
669 const DetLayer* tidLayerpos = posTIDLayers[misLayer - k_LayersAtTOBEnd - 1];
671 if (tTopo->tidSide(iidd) == 1) {
672 tmpTmeas = layerMeasurements.measurements(*tidLayerneg, tsos, prop, chi2Estimator);
673 }
else if (tTopo->tidSide(iidd) == 2) {
674 tmpTmeas = layerMeasurements.measurements(*tidLayerpos, tsos, prop, chi2Estimator);
678 if (misLayer > k_LayersStart && misLayer < k_LayersAtTOBEnd) {
680 std::vector<BarrelDetLayer const*> barrelTIBLayers = measTracker.geometricSearchTracker()->tibLayers();
681 std::vector<BarrelDetLayer const*> barrelTOBLayers = measTracker.geometricSearchTracker()->tobLayers();
683 if (misLayer > k_LayersStart && misLayer <= k_LayersAtTIBEnd) {
684 const DetLayer* tibLayer = barrelTIBLayers[misLayer - k_LayersStart - 1];
685 tmpTmeas = layerMeasurements.measurements(*tibLayer, tsos, prop, chi2Estimator);
686 }
else if (misLayer > k_LayersAtTIBEnd && misLayer < k_LayersAtTOBEnd) {
687 const DetLayer* tobLayer = barrelTOBLayers[misLayer - k_LayersAtTIBEnd - 1];
688 tmpTmeas = layerMeasurements.measurements(*tobLayer, tsos, prop, chi2Estimator);
692 if ((
int(TKlayers) > k_LayersStart &&
int(TKlayers) <= k_LayersAtTIBEnd) &&
int(prev_TKlayers) == 12) {
693 const DetLayer* detlayer = itm->layer();
696 std::vector<DetLayer::DetWithState> compatDets = detlayer->
compatibleDets(tsos, prop, chi2Estimator);
697 std::vector<ForwardDetLayer const*> negTIDLayers = measTracker.geometricSearchTracker()->negTidLayers();
698 std::vector<ForwardDetLayer const*> posTIDLayers = measTracker.geometricSearchTracker()->posTidLayers();
700 const DetLayer* tidLayerneg = negTIDLayers[k_LayersStart];
701 const DetLayer* tidLayerpos = posTIDLayers[k_LayersStart];
702 if (tTopo->tidSide(iidd) == 1) {
703 tmpTmeas = layerMeasurements.measurements(*tidLayerneg, tsos, prop, chi2Estimator);
704 }
else if (tTopo->tidSide(iidd) == 2) {
705 tmpTmeas = layerMeasurements.measurements(*tidLayerpos, tsos, prop, chi2Estimator);
709 if ((
int(TKlayers) > k_LayersAtTIBEnd &&
int(TKlayers) <= k_LayersAtTOBEnd) &&
int(prev_TKlayers) == 15) {
710 const DetLayer* detlayer = itm->layer();
713 std::vector<DetLayer::DetWithState> compatDets = detlayer->
compatibleDets(tsos, prop, chi2Estimator);
715 std::vector<ForwardDetLayer const*> negTECLayers = measTracker.geometricSearchTracker()->negTecLayers();
716 std::vector<ForwardDetLayer const*> posTECLayers = measTracker.geometricSearchTracker()->posTecLayers();
718 const DetLayer* tecLayerneg = negTECLayers[k_LayersStart];
719 const DetLayer* tecLayerpos = posTECLayers[k_LayersStart];
720 if (tTopo->tecSide(iidd) == 1) {
721 tmpTmeas = layerMeasurements.measurements(*tecLayerneg, tsos, prop, chi2Estimator);
722 }
else if (tTopo->tecSide(iidd) == 2) {
723 tmpTmeas = layerMeasurements.measurements(*tecLayerpos, tsos, prop, chi2Estimator);
727 if (!tmpTmeas.empty()) {
729 unsigned int iidd_tmp = TM_tmp.
recHit()->geographicalId().rawId();
731 LogDebug(
"SiStripHitEfficiency:HitEff") <<
" hit actually being added to TM vector";
734 if (::isDoubleSided(iidd_tmp, tTopo)) {
739 missingHitAdded =
true;
745 prev_TKlayers = TKlayers;
750 bool hitsWithBias =
false;
751 for (
auto ilayer : missedLayers) {
752 if (ilayer < TKlayers)
765 const auto nextId = (itm + 1 != TMeas.end()) ? (itm + 1)->recHit()->geographicalId() :
DetId{};
767 if (TKlayers == 9 && theInHit->isValid() && !((!nextId.null()) && (::checkLayer(nextId.rawId(), tTopo) == 9))) {
770 const DetLayer* tob6 = measTracker.geometricSearchTracker()->tobLayers().back();
773 const auto tmp = theLayerMeasurements.measurements(*tob6, tsosTOB5, prop, chi2Estimator);
776 LogDebug(
"SiStripHitEfficiencyWorker") <<
"size of TM from propagation = " <<
tmp.size();
781 const auto& tob6TM =
tmp.back();
782 const auto& tob6Hit = tob6TM.recHit();
783 if (tob6Hit->geographicalId().rawId() != 0) {
784 LogDebug(
"SiStripHitEfficiencyWorker") <<
"tob6 hit actually being added to TM vector";
785 TMs.emplace_back(tob6TM, tTopo, tkgeom,
propagator);
791 if (TKlayers == 21 && theInHit->isValid() &&
792 !((!nextId.null()) && (::checkLayer(nextId.rawId(), tTopo) == 21))) {
793 const DetLayer* tec9pos = measTracker.geometricSearchTracker()->posTecLayers().back();
794 const DetLayer* tec9neg = measTracker.geometricSearchTracker()->negTecLayers().back();
801 LogDebug(
"SiStripHitEfficiencyWorker") <<
"there is a problem with TEC 9 extrapolation";
804 std::vector<TrajectoryMeasurement>
tmp;
805 if (tTopo->tecSide(iidd) == 1) {
806 tmp = theLayerMeasurements.measurements(*tec9neg, tsosTEC9, prop, chi2Estimator);
809 if (tTopo->tecSide(iidd) == 2) {
810 tmp = theLayerMeasurements.measurements(*tec9pos, tsosTEC9, prop, chi2Estimator);
818 const auto& tec9TM =
tmp.back();
819 const auto& tec9Hit = tec9TM.recHit();
821 const unsigned int tec9id = tec9Hit->geographicalId().rawId();
822 LogDebug(
"SiStripHitEfficiencyWorker")
823 <<
"tec9id = " << tec9id <<
" is Double sided = " << ::isDoubleSided(tec9id, tTopo)
824 <<
" and 0x3 = " << (tec9id & 0x3);
826 if (tec9Hit->geographicalId().rawId() != 0) {
827 LogDebug(
"SiStripHitEfficiencyWorker") <<
"tec9 hit actually being added to TM vector";
830 if (::isDoubleSided(tec9id, tTopo)) {
831 TMs.emplace_back(tec9TM, tTopo, tkgeom,
propagator, 1);
832 TMs.emplace_back(tec9TM, tTopo, tkgeom,
propagator, 2);
834 TMs.emplace_back(tec9TM, tTopo, tkgeom,
propagator);
839 for (
const auto& tm : TMs) {
853 LogDebug(
"SiStripHitEfficiencyWorker") <<
"After looping over TrajAtValidHit list";
855 LogDebug(
"SiStripHitEfficiencyWorker") <<
"end TMeasurement loop";
857 LogDebug(
"SiStripHitEfficiencyWorker") <<
"end of trajectories loop";
const edm::ESGetToken< MeasurementTracker, CkfComponentsRecord > measTrackerToken_
const edm::EDGetTokenT< DetIdVector > digisVec_token_
MonitorElement * h_nTracksVsPU
virtual std::vector< DetWithState > compatibleDets(const TrajectoryStateOnSurface &startingState, const Propagator &prop, const MeasurementEstimator &est) const
MonitorElement * h_instLumi
void fillForTraj(const TrajectoryAtInvalidHit &tm, const TrackerTopology *tTopo, const TrackerGeometry *tkgeom, const StripClusterParameterEstimator &stripCPE, const SiStripQuality &stripQuality, const DetIdVector &fedErrorIds, const edm::Handle< edm::DetSetVector< SiStripRawDigi >> &commonModeDigis, const edmNew::DetSetVector< SiStripCluster > &theClusters, int bunchCrossing, float instLumi, float PU, bool highPurity)
const edm::EDGetTokenT< std::vector< Trajectory > > trajectories_token_
std::vector< int > missHitPerLayer
float instLumi() const
Return the luminosity for the current nibble.
std::unordered_map< uint32_t, int > fedErrorCounts
std::unique_ptr< TkHistoMap > FEDErrorOccupancy
const edm::EDGetTokenT< MeasurementTrackerEvent > trackerEvent_token_
const edm::ESGetToken< StripClusterParameterEstimator, TkStripCPERecord > stripCPEToken_
const edm::ESGetToken< Propagator, TrackingComponentsRecord > propagatorToken_
const edm::EDGetTokenT< OnlineLuminosityRecord > metaDataToken_
const edm::ESGetToken< Chi2MeasurementEstimatorBase, TrackingComponentsRecord > chi2EstimatorToken_
const edm::ESGetToken< TrackerGeometry, TrackerDigiGeometryRecord > tkGeomToken_
unsigned int trackMultiplicityCut_
bool useAllHitsFromTracksWithMissingHits_
dqm::reco::MonitorElement * EventStats
const edm::ESGetToken< TrackerTopology, TrackerTopologyRcd > tTopoToken_
const edm::ESGetToken< MagneticField, IdealMagneticFieldRecord > magFieldToken_
SiStripHitEffData calibData_
const edm::EDGetTokenT< TrajTrackAssociationCollection > trajTrackAsso_token_
const edm::EDGetTokenT< DetIdCollection > digisCol_token_
MonitorElement * h_nTracks
const edm::EDGetTokenT< edmNew::DetSetVector< SiStripCluster > > clusters_token_
const edm::EDGetTokenT< reco::TrackCollection > combinatorialTracks_token_
const edm::ESGetToken< SiStripQuality, SiStripQualityRcd > stripQualityToken_
const edm::EDGetTokenT< LumiScalersCollection > scalerToken_
std::vector< unsigned int > hitRecoveryCounters
const edm::EDGetTokenT< edm::DetSetVector< SiStripRawDigi > > commonModeToken_
Log< level::Warning, false > LogWarning
float avgPileUp() const
Return the average pileup for th current nibble.
ConstRecHitPointer const & recHit() const
bool doMissingHitsRecovery_