Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-05-23 03:12:46 

0001 ////////////////////////////////////////////////////////////////////////////////
0002 // Package:          CalibTracker/SiStripHitEfficiency
0003 // Class:            SiStripHitEfficiencyWorker
0004 // Original Author:  Pieter David
0005 //
0006 // Adapted from      HitEff (Keith Ulmer -- University of Colorado, keith.ulmer@colorado.edu
0007 //                   SiStripHitEffFromCalibTree (Christopher Edelmaier)
0008 ///////////////////////////////////////////////////////////////////////////////
0009 
0010 #include "CalibFormats/SiStripObjects/interface/SiStripQuality.h"
0011 #include "CalibTracker/Records/interface/SiStripQualityRcd.h"
0012 #include "CalibTracker/SiStripHitEfficiency/interface/SiStripHitEffData.h"
0013 #include "CalibTracker/SiStripHitEfficiency/interface/SiStripHitEfficiencyHelpers.h"
0014 #include "CalibTracker/SiStripHitEfficiency/interface/TrajectoryAtInvalidHit.h"
0015 #include "DQM/SiStripCommon/interface/TkHistoMap.h"
0016 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0017 #include "DataFormats/Common/interface/DetSetVector.h"
0018 #include "DataFormats/Common/interface/DetSetVectorNew.h"
0019 #include "DataFormats/Common/interface/Handle.h"
0020 #include "DataFormats/DetId/interface/DetIdCollection.h"
0021 #include "DataFormats/DetId/interface/DetIdVector.h"
0022 #include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementError.h"
0023 #include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementVector.h"
0024 #include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h"
0025 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0026 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0027 #include "DataFormats/GeometryVector/interface/LocalVector.h"
0028 #include "DataFormats/OnlineMetaData/interface/OnlineLuminosityRecord.h"
0029 #include "DataFormats/Scalers/interface/LumiScalers.h"
0030 #include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
0031 #include "DataFormats/SiStripCommon/interface/ConstantsForHardwareSystems.h" /* for STRIPS_PER_APV*/
0032 #include "DataFormats/SiStripDigi/interface/SiStripRawDigi.h"
0033 #include "DataFormats/TrackReco/interface/Track.h"
0034 #include "DataFormats/TrackReco/interface/TrackBase.h"
0035 #include "DataFormats/TrackReco/interface/TrackExtra.h"
0036 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0037 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0038 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0039 #include "FWCore/Framework/interface/Event.h"
0040 #include "FWCore/Framework/interface/EventSetup.h"
0041 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0042 #include "FWCore/ParameterSet/interface/ParameterDescription.h"
0043 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0044 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0045 #include "FWCore/Utilities/interface/Exception.h"
0046 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0047 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0048 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0049 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0050 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0051 #include "MagneticField/Engine/interface/MagneticField.h"
0052 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0053 #include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h"
0054 #include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h"
0055 #include "RecoTracker/Record/interface/CkfComponentsRecord.h"
0056 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0057 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
0058 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
0059 #include "TrackingTools/MaterialEffects/interface/PropagatorWithMaterial.h"
0060 #include "TrackingTools/MeasurementDet/interface/LayerMeasurements.h"
0061 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0062 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0063 
0064 class SiStripHitEfficiencyWorker : public DQMEDAnalyzer {
0065 public:
0066   explicit SiStripHitEfficiencyWorker(const edm::ParameterSet& conf);
0067   ~SiStripHitEfficiencyWorker() override = default;
0068   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0069 
0070 private:
0071   void bookHistograms(DQMStore::IBooker& booker, const edm::Run& run, const edm::EventSetup& setup) override;
0072   void analyze(const edm::Event& e, const edm::EventSetup& c) override;
0073   void fillForTraj(const TrajectoryAtInvalidHit& tm,
0074                    const TrackerTopology* tTopo,
0075                    const TrackerGeometry* tkgeom,
0076                    const StripClusterParameterEstimator& stripCPE,
0077                    const SiStripQuality& stripQuality,
0078                    const DetIdVector& fedErrorIds,
0079                    const edm::Handle<edm::DetSetVector<SiStripRawDigi>>& commonModeDigis,
0080                    const edmNew::DetSetVector<SiStripCluster>& theClusters,
0081                    int bunchCrossing,
0082                    float instLumi,
0083                    float PU,
0084                    bool highPurity);
0085 
0086   // ----------member data ---------------------------
0087   SiStripHitEffData calibData_;
0088 
0089   // event data tokens
0090   const edm::EDGetTokenT<LumiScalersCollection> scalerToken_;
0091   const edm::EDGetTokenT<OnlineLuminosityRecord> metaDataToken_;
0092   const edm::EDGetTokenT<edm::DetSetVector<SiStripRawDigi>> commonModeToken_;
0093   const edm::EDGetTokenT<reco::TrackCollection> combinatorialTracks_token_;
0094   const edm::EDGetTokenT<std::vector<Trajectory>> trajectories_token_;
0095   const edm::EDGetTokenT<TrajTrackAssociationCollection> trajTrackAsso_token_;
0096   const edm::EDGetTokenT<edmNew::DetSetVector<SiStripCluster>> clusters_token_;
0097   const edm::EDGetTokenT<DetIdCollection> digisCol_token_;
0098   const edm::EDGetTokenT<DetIdVector> digisVec_token_;
0099   const edm::EDGetTokenT<MeasurementTrackerEvent> trackerEvent_token_;
0100 
0101   // event setup tokens
0102   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> tTopoToken_;
0103   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> tkGeomToken_;
0104   const edm::ESGetToken<StripClusterParameterEstimator, TkStripCPERecord> stripCPEToken_;
0105   const edm::ESGetToken<SiStripQuality, SiStripQualityRcd> stripQualityToken_;
0106   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magFieldToken_;
0107   const edm::ESGetToken<MeasurementTracker, CkfComponentsRecord> measTrackerToken_;
0108   const edm::ESGetToken<Chi2MeasurementEstimatorBase, TrackingComponentsRecord> chi2EstimatorToken_;
0109   const edm::ESGetToken<Propagator, TrackingComponentsRecord> propagatorToken_;
0110   const edm::ESGetToken<TkDetMap, TrackerTopologyRcd> tkDetMapToken_;
0111 
0112   // configurable parameters
0113   std::string dqmDir_;
0114   unsigned int layers_;
0115   bool DEBUG_;
0116   bool addLumi_;
0117   bool addCommonMode_;
0118   bool cutOnTracks_;
0119   unsigned int trackMultiplicityCut_;
0120   bool useFirstMeas_;
0121   bool useLastMeas_;
0122   bool useAllHitsFromTracksWithMissingHits_;
0123   bool doMissingHitsRecovery_;
0124   unsigned int clusterMatchingMethod_;
0125   float resXSig_;
0126   float clusterTracjDist_;
0127   float stripsApvEdge_;
0128   bool useOnlyHighPurityTracks_;
0129   int bunchX_;
0130   bool showRings_;
0131   bool showTOB6TEC9_;
0132   unsigned int nTEClayers_;
0133 
0134   // output file
0135   std::set<uint32_t> badModules_;
0136 
0137   // for the missing hits recovery
0138   std::vector<unsigned int> hitRecoveryCounters;
0139   std::vector<unsigned int> hitTotalCounters;
0140   int totalNbHits;
0141   std::vector<int> missHitPerLayer;
0142 
0143   struct EffME1 {
0144     EffME1() : hTotal(nullptr), hFound(nullptr) {}
0145     EffME1(MonitorElement* total, MonitorElement* found) : hTotal(total), hFound(found) {}
0146 
0147     void fill(double x, bool found, float weight = 1.) {
0148       hTotal->Fill(x, weight);
0149       if (found) {
0150         hFound->Fill(x, weight);
0151       }
0152     }
0153 
0154     MonitorElement *hTotal, *hFound;
0155   };
0156   struct EffTkMap {
0157     EffTkMap() : hTotal(nullptr), hFound(nullptr) {}
0158     EffTkMap(std::unique_ptr<TkHistoMap>&& total, std::unique_ptr<TkHistoMap>&& found)
0159         : hTotal(std::move(total)), hFound(std::move(found)) {}
0160 
0161     void fill(uint32_t id, bool found, float weight = 1.) {
0162       hTotal->fill(id, weight);
0163       if (found) {
0164         hFound->fill(id, weight);
0165       }
0166     }
0167 
0168     bool check(uint32_t id) {
0169       if (hTotal->getValue(id) < hFound->getValue(id)) {
0170         return false;
0171       } else {
0172         return true;
0173       }
0174     }
0175 
0176     std::unique_ptr<TkHistoMap> hTotal, hFound;
0177   };
0178 
0179   MonitorElement *h_bx, *h_instLumi, *h_PU;
0180   MonitorElement *h_nTracks, *h_nTracksVsPU;
0181   EffME1 h_goodLayer;
0182   EffME1 h_allLayer;
0183   EffME1 h_layer;
0184   std::vector<MonitorElement*> h_resolution;
0185   std::vector<EffME1> h_layer_vsLumi;
0186   std::vector<EffME1> h_layer_vsBx;
0187   std::vector<EffME1> h_layer_vsPU;
0188   std::vector<EffME1> h_layer_vsCM;
0189   std::vector<MonitorElement*> h_hotcold;
0190 
0191   EffTkMap h_module;
0192 };
0193 
0194 //
0195 // constructors and destructor
0196 //
0197 
0198 SiStripHitEfficiencyWorker::SiStripHitEfficiencyWorker(const edm::ParameterSet& conf)
0199     : scalerToken_(consumes<LumiScalersCollection>(conf.getParameter<edm::InputTag>("lumiScalers"))),
0200       metaDataToken_(consumes<OnlineLuminosityRecord>(conf.getParameter<edm::InputTag>("metadata"))),
0201       commonModeToken_(mayConsume<edm::DetSetVector<SiStripRawDigi>>(conf.getParameter<edm::InputTag>("commonMode"))),
0202       combinatorialTracks_token_(
0203           consumes<reco::TrackCollection>(conf.getParameter<edm::InputTag>("combinatorialTracks"))),
0204       trajectories_token_(consumes<std::vector<Trajectory>>(conf.getParameter<edm::InputTag>("trajectories"))),
0205       trajTrackAsso_token_(consumes<TrajTrackAssociationCollection>(conf.getParameter<edm::InputTag>("trajectories"))),
0206       clusters_token_(
0207           consumes<edmNew::DetSetVector<SiStripCluster>>(conf.getParameter<edm::InputTag>("siStripClusters"))),
0208       digisCol_token_(consumes(conf.getParameter<edm::InputTag>("siStripDigis"))),
0209       digisVec_token_(consumes(conf.getParameter<edm::InputTag>("siStripDigis"))),
0210       trackerEvent_token_(consumes<MeasurementTrackerEvent>(conf.getParameter<edm::InputTag>("trackerEvent"))),
0211       tTopoToken_(esConsumes()),
0212       tkGeomToken_(esConsumes()),
0213       stripCPEToken_(esConsumes(edm::ESInputTag{"", "StripCPEfromTrackAngle"})),
0214       stripQualityToken_(esConsumes()),
0215       magFieldToken_(esConsumes()),
0216       measTrackerToken_(esConsumes()),
0217       chi2EstimatorToken_(esConsumes(edm::ESInputTag{"", "Chi2"})),
0218       propagatorToken_(esConsumes(edm::ESInputTag{"", "PropagatorWithMaterial"})),
0219       tkDetMapToken_(esConsumes<edm::Transition::BeginRun>()),
0220       dqmDir_(conf.getParameter<std::string>("dqmDir")),
0221       layers_(conf.getParameter<int>("Layer")),
0222       DEBUG_(conf.getUntrackedParameter<bool>("Debug", false)),
0223       addLumi_(conf.getUntrackedParameter<bool>("addLumi", false)),
0224       addCommonMode_(conf.getUntrackedParameter<bool>("addCommonMode", false)),
0225       cutOnTracks_(conf.getParameter<bool>("cutOnTracks")),
0226       trackMultiplicityCut_(conf.getParameter<unsigned int>("trackMultiplicity")),
0227       useFirstMeas_(conf.getParameter<bool>("useFirstMeas")),
0228       useLastMeas_(conf.getParameter<bool>("useLastMeas")),
0229       useAllHitsFromTracksWithMissingHits_(conf.getParameter<bool>("useAllHitsFromTracksWithMissingHits")),
0230       doMissingHitsRecovery_(conf.getParameter<bool>("doMissingHitsRecovery")),
0231       clusterMatchingMethod_(conf.getParameter<int>("ClusterMatchingMethod")),
0232       resXSig_(conf.getParameter<double>("ResXSig")),
0233       clusterTracjDist_(conf.getParameter<double>("ClusterTrajDist")),
0234       stripsApvEdge_(conf.getParameter<double>("StripsApvEdge")),
0235       useOnlyHighPurityTracks_(conf.getParameter<bool>("UseOnlyHighPurityTracks")),
0236       bunchX_(conf.getUntrackedParameter<int>("BunchCrossing", 0)),
0237       showRings_(conf.getUntrackedParameter<bool>("ShowRings", false)),
0238       showTOB6TEC9_(conf.getUntrackedParameter<bool>("ShowTOB6TEC9", false)) {
0239   hitRecoveryCounters.resize(k_END_OF_LAYERS, 0);
0240   hitTotalCounters.resize(k_END_OF_LAYERS, 0);
0241   missHitPerLayer.resize(k_END_OF_LAYERS, 0);
0242   totalNbHits = 0;
0243 
0244   nTEClayers_ = (showRings_ ? 7 : 9);  // number of rings or wheels
0245 
0246   const std::string badModulesFile = conf.getUntrackedParameter<std::string>("BadModulesFile", "");
0247   if (!badModulesFile.empty()) {
0248     std::ifstream badModules_file(badModulesFile);
0249     uint32_t badmodule_detid;
0250     int mods, fiber1, fiber2, fiber3;
0251     if (badModules_file.is_open()) {
0252       std::string line;
0253       while (getline(badModules_file, line)) {
0254         if (badModules_file.eof())
0255           continue;
0256         std::stringstream ss(line);
0257         ss >> badmodule_detid >> mods >> fiber1 >> fiber2 >> fiber3;
0258         if (badmodule_detid != 0 && mods == 1 && (fiber1 == 1 || fiber2 == 1 || fiber3 == 1))
0259           badModules_.insert(badmodule_detid);
0260       }
0261       badModules_file.close();
0262     }
0263   }
0264   if (!badModules_.empty())
0265     LogDebug("SiStripHitEfficiencyWorker") << "Remove additionnal bad modules from the analysis: ";
0266   for (const auto badMod : badModules_) {
0267     LogDebug("SiStripHitEfficiencyWorker") << " " << badMod;
0268   }
0269 }
0270 
0271 void SiStripHitEfficiencyWorker::bookHistograms(DQMStore::IBooker& booker,
0272                                                 const edm::Run& run,
0273                                                 const edm::EventSetup& setup) {
0274   booker.setCurrentFolder(fmt::format("{}/EventInfo", dqmDir_));
0275   h_bx = booker.book1D("bx", "bx", 3600, 0, 3600);
0276   h_instLumi = booker.book1D("instLumi", "inst. lumi.", 250, 0, 25000);
0277   h_PU = booker.book1D("PU", "PU", 200, 0, 200);
0278   h_nTracks = booker.book1D("ntracks", "n.tracks;n. tracks;n.events", 500, -0.5, 499.5);
0279   h_nTracksVsPU = booker.bookProfile("nTracksVsPU", "n. tracks vs PU; PU; n.tracks ", 200, 0, 200, 500, -0.5, 499.5);
0280 
0281   calibData_.EventStats = booker.book2I("EventStats", "Statistics", 3, -0.5, 2.5, 1, 0, 1);
0282   calibData_.EventStats->setBinLabel(1, "events count", 1);
0283   calibData_.EventStats->setBinLabel(2, "tracks count", 1);
0284   calibData_.EventStats->setBinLabel(3, "measurements count", 1);
0285 
0286   booker.setCurrentFolder(dqmDir_);
0287   h_goodLayer = EffME1(booker.book1D("goodlayer_total", "goodlayer_total", 35, 0., 35.),
0288                        booker.book1D("goodlayer_found", "goodlayer_found", 35, 0., 35.));
0289   h_allLayer = EffME1(booker.book1D("alllayer_total", "alllayer_total", 35, 0., 35.),
0290                       booker.book1D("alllayer_found", "alllayer_found", 35, 0., 35.));
0291 
0292   h_layer = EffME1(
0293       booker.book1D(
0294           "layer_found", "layer_found", bounds::k_END_OF_LAYERS, 0., static_cast<float>(bounds::k_END_OF_LAYERS)),
0295       booker.book1D(
0296           "layer_total", "layer_total", bounds::k_END_OF_LAYERS, 0., static_cast<float>(bounds::k_END_OF_LAYERS)));
0297 
0298   for (int layer = 1; layer != bounds::k_END_OF_LAYERS; ++layer) {
0299     const auto lyrName = ::layerName(layer, showRings_, nTEClayers_);
0300 
0301     // book resolutions
0302     booker.setCurrentFolder(fmt::format("{}/Resolutions", dqmDir_));
0303     auto ihres = booker.book1D(Form("resol_layer_%i", layer), lyrName, 125, -125., 125.);
0304     ihres->setAxisTitle("trajX-clusX [strip unit]");
0305     h_resolution.push_back(ihres);
0306 
0307     // book plots vs Lumi
0308     booker.setCurrentFolder(fmt::format("{}/VsLumi", dqmDir_));
0309     h_layer_vsLumi.push_back(EffME1(booker.book1D(Form("layertotal_vsLumi_layer_%i", layer), lyrName, 100, 0, 25000),
0310                                     booker.book1D(Form("layerfound_vsLumi_layer_%i", layer), lyrName, 100, 0, 25000)));
0311 
0312     // book plots vs Lumi
0313     booker.setCurrentFolder(fmt::format("{}/VsPu", dqmDir_));
0314     h_layer_vsPU.push_back(EffME1(booker.book1D(Form("layertotal_vsPU_layer_%i", layer), lyrName, 45, 0, 90),
0315                                   booker.book1D(Form("layerfound_vsPU_layer_%i", layer), lyrName, 45, 0, 90)));
0316     if (addCommonMode_) {
0317       // book plots for common mode
0318       booker.setCurrentFolder(fmt::format("{}/CommonMode", dqmDir_));
0319       h_layer_vsCM.push_back(EffME1(booker.book1D(Form("layertotal_vsCM_layer_%i", layer), lyrName, 20, 0, 400),
0320                                     booker.book1D(Form("layerfound_vsCM_layer_%i", layer), lyrName, 20, 0, 400)));
0321     }
0322 
0323     // book plots vs Lumi
0324     booker.setCurrentFolder(fmt::format("{}/VsBx", dqmDir_));
0325     h_layer_vsBx.push_back(EffME1(
0326         booker.book1D(Form("totalVsBx_layer%i", layer), Form("layer %i (%s)", layer, lyrName.c_str()), 3565, 0, 3565),
0327         booker.book1D(Form("foundVsBx_layer%i", layer), Form("layer %i (%s)", layer, lyrName.c_str()), 3565, 0, 3565)));
0328 
0329     // book hot and cold
0330     booker.setCurrentFolder(fmt::format("{}/MissingHits", dqmDir_));
0331     if (layer <= bounds::k_LayersAtTOBEnd) {
0332       const bool isTIB = layer <= bounds::k_LayersAtTIBEnd;
0333       const auto partition = (isTIB ? "TIB" : "TOB");
0334       const auto yMax = (isTIB ? 100 : 120);
0335 
0336       const auto& tit =
0337           Form("%s%i: Map of missing hits", partition, (isTIB ? layer : layer - bounds::k_LayersAtTIBEnd));
0338 
0339       // histogram name must not contain ":" otherwise it fails upload to the GUI
0340       // see https://github.com/cms-DQM/dqmgui_prod/blob/af0a388e8f57c60e51111585d298aeeea943367f/src/cpp/DQM/DQMStore.cc#L56
0341       std::string name{tit};
0342       ::replaceInString(name, ":", "");
0343 
0344       auto ihhotcold = booker.book2D(name, tit, 100, -1, 361, 100, -yMax, yMax);
0345       ihhotcold->setAxisTitle("#phi [deg]", 1);
0346       ihhotcold->setBinLabel(1, "360", 1);
0347       ihhotcold->setBinLabel(50, "180", 1);
0348       ihhotcold->setBinLabel(100, "0", 1);
0349       ihhotcold->setAxisTitle("Global Z [cm]", 2);
0350       ihhotcold->setOption("colz");
0351       h_hotcold.push_back(ihhotcold);
0352     } else {
0353       const bool isTID = layer <= bounds::k_LayersAtTIDEnd;
0354       const auto partitions =
0355           (isTID ? std::vector<std::string>{"TIDplus", "TIDminus"} : std::vector<std::string>{"TECplus", "TECminus"});
0356       const auto axMax = (isTID ? 100 : 120);
0357       for (const auto& part : partitions) {
0358         // create the title by replacing the minus/plus symbols
0359         std::string forTitle{part};
0360         ::replaceInString(forTitle, "minus", "-");
0361         ::replaceInString(forTitle, "plus", "+");
0362 
0363         // histogram name must not contain ":" otherwise it fails upload to the GUI
0364         // see https://github.com/cms-DQM/dqmgui_prod/blob/af0a388e8f57c60e51111585d298aeeea943367f/src/cpp/DQM/DQMStore.cc#L56
0365         const auto& name = Form("%s %i Map of Missing Hits",
0366                                 part.c_str(),
0367                                 (isTID ? layer - bounds::k_LayersAtTOBEnd : layer - bounds::k_LayersAtTIDEnd));
0368         const auto& tit = Form("%s%i: Map of Missing Hits",
0369                                forTitle.c_str(),
0370                                (isTID ? layer - bounds::k_LayersAtTOBEnd : layer - bounds::k_LayersAtTIDEnd));
0371 
0372         auto ihhotcold = booker.book2D(name, tit, 100, -axMax, axMax, 100, -axMax, axMax);
0373         ihhotcold->setAxisTitle("Global Y", 1);
0374         ihhotcold->setBinLabel(1, "+Y", 1);
0375         ihhotcold->setBinLabel(50, "0", 1);
0376         ihhotcold->setBinLabel(100, "-Y", 1);
0377         ihhotcold->setAxisTitle("Global X", 2);
0378         ihhotcold->setBinLabel(1, "-X", 2);
0379         ihhotcold->setBinLabel(50, "0", 2);
0380         ihhotcold->setBinLabel(100, "+X", 2);
0381         ihhotcold->setOption("colz");
0382         h_hotcold.push_back(ihhotcold);
0383       }
0384     }
0385   }
0386 
0387   // come back to the main folder
0388   booker.setCurrentFolder(dqmDir_);
0389   const auto tkDetMapFolder = fmt::format("{}/TkDetMaps", dqmDir_);
0390 
0391   const TkDetMap* tkDetMap = &setup.getData(tkDetMapToken_);
0392   h_module =
0393       EffTkMap(std::make_unique<TkHistoMap>(tkDetMap, booker, tkDetMapFolder, "perModule_total", 0, false, true),
0394                std::make_unique<TkHistoMap>(tkDetMap, booker, tkDetMapFolder, "perModule_found", 0, false, true));
0395 
0396   // fill the FED Errors
0397   booker.setCurrentFolder(dqmDir_);
0398   const auto FEDErrorMapFolder = fmt::format("{}/FEDErrorTkDetMaps", dqmDir_);
0399   calibData_.FEDErrorOccupancy =
0400       std::make_unique<TkHistoMap>(tkDetMap, booker, FEDErrorMapFolder, "perModule_FEDErrors", 0, false, true);
0401 }
0402 
0403 void SiStripHitEfficiencyWorker::analyze(const edm::Event& e, const edm::EventSetup& es) {
0404   const auto tTopo = &es.getData(tTopoToken_);
0405 
0406   //  bool DEBUG_ = false;
0407 
0408   LogDebug("SiStripHitEfficiencyWorker") << "beginning analyze from HitEff";
0409 
0410   // Step A: Get Inputs
0411 
0412   // Luminosity informations
0413   edm::Handle<LumiScalersCollection> lumiScalers = e.getHandle(scalerToken_);
0414   edm::Handle<OnlineLuminosityRecord> metaData = e.getHandle(metaDataToken_);
0415 
0416   float instLumi = 0;
0417   float PU = 0;
0418   if (addLumi_) {
0419     if (lumiScalers.isValid() && !lumiScalers->empty()) {
0420       if (lumiScalers->begin() != lumiScalers->end()) {
0421         instLumi = lumiScalers->begin()->instantLumi();
0422         PU = lumiScalers->begin()->pileup();
0423       }
0424     } else if (metaData.isValid()) {
0425       instLumi = metaData->instLumi();
0426       PU = metaData->avgPileUp();
0427     } else {
0428       edm::LogWarning("SiStripHitEfficiencyWorker") << "could not find a source for the Luminosity and PU";
0429     }
0430   }
0431 
0432   h_bx->Fill(e.bunchCrossing());
0433   h_instLumi->Fill(instLumi);
0434   h_PU->Fill(PU);
0435 
0436   edm::Handle<edm::DetSetVector<SiStripRawDigi>> commonModeDigis;
0437   if (addCommonMode_)
0438     e.getByToken(commonModeToken_, commonModeDigis);
0439 
0440   edm::Handle<reco::TrackCollection> tracksCKF;
0441   e.getByToken(combinatorialTracks_token_, tracksCKF);
0442 
0443   edm::Handle<std::vector<Trajectory>> TrajectoryCollectionCKF;
0444   e.getByToken(trajectories_token_, TrajectoryCollectionCKF);
0445 
0446   edm::Handle<TrajTrackAssociationCollection> trajTrackAssociationHandle;
0447   e.getByToken(trajTrackAsso_token_, trajTrackAssociationHandle);
0448 
0449   edm::Handle<edmNew::DetSetVector<SiStripCluster>> theClusters;
0450   e.getByToken(clusters_token_, theClusters);
0451 
0452   // get the list of module IDs with FED-detected errors
0453   //  - In Aug-2023, the data format was changed from DetIdCollection to DetIdVector.
0454   //  - To provide some level of backward-compatibility,
0455   //    the plugin checks for both types giving preference to the new format.
0456   //  - If only the old format is available, the collection is
0457   //    converted to the new format, then used downstream.
0458   auto const& fedErrorIdsCol_h = e.getHandle(digisCol_token_);
0459   auto const& fedErrorIdsVec_h = e.getHandle(digisVec_token_);
0460   if (not fedErrorIdsCol_h.isValid() and not fedErrorIdsVec_h.isValid()) {
0461     throw cms::Exception("InvalidProductSiStripDetIdsWithFEDErrors")
0462         << "no valid product for SiStrip DetIds with FED errors (see parameter \"siStripDigis\"), "
0463            "neither for new format (DetIdVector) nor old format (DetIdCollection)";
0464   }
0465   auto const& fedErrorIds = fedErrorIdsVec_h.isValid() ? *fedErrorIdsVec_h : fedErrorIdsCol_h->as_vector();
0466 
0467   // fill the calibData with the FEDErrors
0468   for (const auto& fedErr : fedErrorIds) {
0469     // fill the TkHistoMap occupancy map
0470     calibData_.FEDErrorOccupancy->fill(fedErr.rawId(), 1.);
0471 
0472     // fill the unordered map
0473     if (calibData_.fedErrorCounts.find(fedErr.rawId()) != calibData_.fedErrorCounts.end()) {
0474       calibData_.fedErrorCounts[fedErr.rawId()] += 1;
0475     } else {
0476       calibData_.fedErrorCounts.insert(std::make_pair(fedErr.rawId(), 1));
0477     }
0478   }
0479 
0480   edm::Handle<MeasurementTrackerEvent> measurementTrackerEvent;
0481   e.getByToken(trackerEvent_token_, measurementTrackerEvent);
0482 
0483   const auto tkgeom = &es.getData(tkGeomToken_);
0484   const auto& stripcpe = es.getData(stripCPEToken_);
0485   const auto& stripQuality = es.getData(stripQualityToken_);
0486   const auto& magField = es.getData(magFieldToken_);
0487   const auto& measTracker = es.getData(measTrackerToken_);
0488   const auto& chi2Estimator = es.getData(chi2EstimatorToken_);
0489   const auto& prop = es.getData(propagatorToken_);
0490 
0491   // Tracking
0492   LogDebug("SiStripHitEfficiencyWorker") << "number ckf tracks found = " << tracksCKF->size();
0493 
0494   h_nTracks->Fill(tracksCKF->size());
0495   h_nTracksVsPU->Fill(PU, tracksCKF->size());
0496 
0497   // bin 0: one entry for each event
0498   calibData_.EventStats->Fill(0., 0., 1);
0499   // bin 1: one entry for each track
0500   calibData_.EventStats->Fill(1., 0., tracksCKF->size());
0501 
0502   if (!tracksCKF->empty()) {
0503     if (cutOnTracks_ && (tracksCKF->size() >= trackMultiplicityCut_))
0504       return;
0505     if (cutOnTracks_)
0506       LogDebug("SiStripHitEfficiencyWorker")
0507           << "starting checking good event with < " << trackMultiplicityCut_ << " tracks";
0508 
0509     // actually should do a loop over all the tracks in the event here
0510 
0511     // Looping over traj-track associations to be able to get traj & track informations
0512     for (const auto& trajTrack : *trajTrackAssociationHandle) {
0513       // for each track, fill some variables such as number of hits and momentum
0514 
0515       const bool highPurity = trajTrack.val->quality(reco::TrackBase::TrackQuality::highPurity);
0516       auto TMeas = trajTrack.key->measurements();
0517       totalNbHits += int(TMeas.size());
0518 
0519       /*
0520       const bool hasMissingHits = std::any_of(std::begin(TMeas), std::end(TMeas), [](const auto& tm) {
0521         return tm.recHit()->getType() == TrackingRecHit::Type::missing;
0522       });
0523       */
0524 
0525       // Check whether the trajectory has some missing hits
0526       bool hasMissingHits{false};
0527       int previous_layer{999};
0528       std::vector<unsigned int> missedLayers;
0529 
0530       for (const auto& itm : TMeas) {
0531         auto theHit = itm.recHit();
0532         unsigned int iidd = theHit->geographicalId().rawId();
0533         int layer = ::checkLayer(iidd, tTopo);
0534         int missedLayer = layer + 1;
0535         int previousMissedLayer = (layer + 2);
0536         int diffPreviousLayer = (layer - previous_layer);
0537         if (doMissingHitsRecovery_) {
0538           //Layers from TIB + TOB
0539           if (diffPreviousLayer == -2 && missedLayer > k_LayersStart && missedLayer < k_LayersAtTOBEnd) {
0540             missHitPerLayer[missedLayer] += 1;
0541             hasMissingHits = true;
0542           }
0543           //Layers from TID
0544           else if (diffPreviousLayer == -2 && (missedLayer > k_LayersAtTOBEnd + 1 && missedLayer <= k_LayersAtTIDEnd)) {
0545             missHitPerLayer[missedLayer] += 1;
0546             hasMissingHits = true;
0547           }
0548           //Layers from TEC
0549           else if (diffPreviousLayer == -2 && missedLayer > k_LayersAtTIDEnd && missedLayer <= k_LayersAtTECEnd) {
0550             missHitPerLayer[missedLayer] += 1;
0551             hasMissingHits = true;
0552           }
0553 
0554           //##### TID Layer 11 in transition TID -> TIB : layer is in TIB, previous layer  = 12
0555           if ((layer > k_LayersStart && layer <= k_LayersAtTIBEnd) && (previous_layer == 12)) {
0556             missHitPerLayer[11] += 1;
0557             hasMissingHits = true;
0558           }
0559 
0560           //##### TEC Layer 14 in transition TEC -> TOB : layer is in TOB, previous layer =  15
0561           if ((layer > k_LayersAtTIBEnd && layer <= k_LayersAtTOBEnd) && (previous_layer == 15)) {
0562             missHitPerLayer[14] += 1;
0563             hasMissingHits = true;
0564           }
0565 
0566           //####### Consecutive missing hits case #######
0567 
0568           //##### Layers from TIB + TOB
0569           if (diffPreviousLayer == -3 && missedLayer > k_LayersStart && missedLayer < k_LayersAtTOBEnd &&
0570               previousMissedLayer > k_LayersStart && previousMissedLayer < k_LayersAtTOBEnd) {
0571             missHitPerLayer[missedLayer] += 1;
0572             missHitPerLayer[previousMissedLayer] += 1;
0573             hasMissingHits = true;
0574           }
0575 
0576           //##### Layers from TEC
0577           else if (diffPreviousLayer == -3 && missedLayer > k_LayersAtTIDEnd && missedLayer <= k_LayersAtTECEnd &&
0578                    previousMissedLayer > k_LayersAtTIDEnd && previousMissedLayer <= k_LayersAtTECEnd) {
0579             missHitPerLayer[missedLayer] += 1;
0580             missHitPerLayer[previousMissedLayer] += 1;
0581             hasMissingHits = true;
0582           }
0583         }
0584         if (theHit->getType() == TrackingRecHit::Type::missing)
0585           hasMissingHits = true;
0586 
0587         if (hasMissingHits)
0588           missedLayers.push_back(layer);
0589         previous_layer = layer;
0590       }
0591 
0592       // Loop on each measurement and take into consideration
0593       //--------------------------------------------------------
0594       unsigned int prev_TKlayers = 0;
0595       for (auto itm = TMeas.cbegin(); itm != TMeas.cend(); ++itm) {
0596         const auto theInHit = (*itm).recHit();
0597 
0598         //bin 2: one entry for each measurement
0599         calibData_.EventStats->Fill(2., 0., 1.);
0600 
0601         LogDebug("SiStripHitEfficiencyWorker") << "theInHit is valid = " << theInHit->isValid();
0602 
0603         unsigned int iidd = theInHit->geographicalId().rawId();
0604 
0605         unsigned int TKlayers = ::checkLayer(iidd, tTopo);
0606 
0607         bool foundConsMissingHits{false};
0608 
0609         // do not bother with pixel hits
0610         if (DetId(iidd).subdetId() < SiStripSubdetector::TIB)
0611           continue;
0612 
0613         LogDebug("SiStripHitEfficiencyWorker") << "TKlayer from trajectory: " << TKlayers << "  from module = " << iidd
0614                                                << "   matched/stereo/rphi = " << ((iidd & 0x3) == 0) << "/"
0615                                                << ((iidd & 0x3) == 1) << "/" << ((iidd & 0x3) == 2);
0616 
0617         // Test first and last points of the trajectory
0618         // the list of measurements starts from outer layers  !!! This could change -> should add a check
0619         bool isFirstMeas = (itm == (TMeas.end() - 1));
0620         bool isLastMeas = (itm == (TMeas.begin()));
0621 
0622         if (!useFirstMeas_ && isFirstMeas)
0623           continue;
0624         if (!useLastMeas_ && isLastMeas)
0625           continue;
0626 
0627         // In case of missing hit in the track, check whether to use the other hits or not.
0628         if (hasMissingHits && theInHit->getType() != TrackingRecHit::Type::missing &&
0629             !useAllHitsFromTracksWithMissingHits_)
0630           continue;
0631 
0632         // If Trajectory measurement from TOB 6 or TEC 9, skip it because it's always valid they are filled later
0633         if (TKlayers == bounds::k_LayersAtTOBEnd || TKlayers == bounds::k_LayersAtTECEnd) {
0634           LogDebug("SiStripHitEfficiencyWorker") << "skipping original TM for TOB 6 or TEC 9";
0635           continue;
0636         }
0637 
0638         std::vector<TrajectoryAtInvalidHit> TMs;
0639 
0640         // Make AnalyticalPropagat // TODO where to save these?or to use in TAVH constructor
0641         AnalyticalPropagator propagator(&magField, anyDirection);
0642 
0643         // for double sided layers check both sensors--if no hit was found on either sensor surface,
0644         // the trajectory measurements only have one invalid hit entry on the matched surface
0645         // so get the TrajectoryAtInvalidHit for both surfaces and include them in the study
0646         if (::isDoubleSided(iidd, tTopo) && ((iidd & 0x3) == 0)) {
0647           // do hit eff check twice--once for each sensor
0648           //add a TM for each surface
0649           TMs.emplace_back(*itm, tTopo, tkgeom, propagator, 1);
0650           TMs.emplace_back(*itm, tTopo, tkgeom, propagator, 2);
0651         } else if (::isDoubleSided(iidd, tTopo) && (!::check2DPartner(iidd, TMeas))) {
0652           // if only one hit was found the trajectory measurement is on that sensor surface, and the other surface from
0653           // the matched layer should be added to the study as well
0654           TMs.emplace_back(*itm, tTopo, tkgeom, propagator, 1);
0655           TMs.emplace_back(*itm, tTopo, tkgeom, propagator, 2);
0656           LogDebug("SiStripHitEfficiencyWorker") << " found a hit with a missing partner";
0657         } else {
0658           //only add one TM for the single surface and the other will be added in the next iteration
0659           TMs.emplace_back(*itm, tTopo, tkgeom, propagator);
0660         }
0661 
0662         bool missingHitAdded{false};
0663         std::vector<TrajectoryMeasurement> tmpTmeas, prev_tmpTmeas;
0664         unsigned int misLayer = TKlayers + 1;
0665         unsigned int previousMisLayer = TKlayers + 2;
0666         //Use bool doMissingHitsRecovery to add possible missing hits based on actual/previous hit
0667         if (doMissingHitsRecovery_) {
0668           if (int(TKlayers) - int(prev_TKlayers) == -2) {
0669             const DetLayer* detlayer = itm->layer();
0670             const LayerMeasurements layerMeasurements{measTracker, *measurementTrackerEvent};
0671             const TrajectoryStateOnSurface tsos = itm->updatedState();
0672             std::vector<DetLayer::DetWithState> compatDets = detlayer->compatibleDets(tsos, prop, chi2Estimator);
0673 
0674             if (misLayer > k_LayersAtTIDEnd && misLayer < k_LayersAtTECEnd) {  //TEC
0675               std::vector<ForwardDetLayer const*> negTECLayers = measTracker.geometricSearchTracker()->negTecLayers();
0676               std::vector<ForwardDetLayer const*> posTECLayers = measTracker.geometricSearchTracker()->posTecLayers();
0677               const DetLayer* tecLayerneg = negTECLayers[misLayer - k_LayersAtTIDEnd - 1];
0678               const DetLayer* tecLayerpos = posTECLayers[misLayer - k_LayersAtTIDEnd - 1];
0679               if (tTopo->tecSide(iidd) == 1) {
0680                 tmpTmeas = layerMeasurements.measurements(*tecLayerneg, tsos, prop, chi2Estimator);
0681               } else if (tTopo->tecSide(iidd) == 2) {
0682                 tmpTmeas = layerMeasurements.measurements(*tecLayerpos, tsos, prop, chi2Estimator);
0683               }
0684             }
0685 
0686             else if (misLayer == (k_LayersAtTIDEnd - 1) ||
0687                      misLayer == k_LayersAtTIDEnd) {  // This is for TID layers 12 and 13
0688               std::vector<ForwardDetLayer const*> negTIDLayers = measTracker.geometricSearchTracker()->negTidLayers();
0689               std::vector<ForwardDetLayer const*> posTIDLayers = measTracker.geometricSearchTracker()->posTidLayers();
0690               const DetLayer* tidLayerneg = negTIDLayers[misLayer - k_LayersAtTOBEnd - 1];
0691               const DetLayer* tidLayerpos = posTIDLayers[misLayer - k_LayersAtTOBEnd - 1];
0692 
0693               if (tTopo->tidSide(iidd) == 1) {
0694                 tmpTmeas = layerMeasurements.measurements(*tidLayerneg, tsos, prop, chi2Estimator);
0695               } else if (tTopo->tidSide(iidd) == 2) {
0696                 tmpTmeas = layerMeasurements.measurements(*tidLayerpos, tsos, prop, chi2Estimator);
0697               }
0698             }
0699 
0700             if (misLayer > k_LayersStart && misLayer < k_LayersAtTOBEnd) {  // Barrel
0701 
0702               std::vector<BarrelDetLayer const*> barrelTIBLayers = measTracker.geometricSearchTracker()->tibLayers();
0703               std::vector<BarrelDetLayer const*> barrelTOBLayers = measTracker.geometricSearchTracker()->tobLayers();
0704 
0705               if (misLayer > k_LayersStart && misLayer <= k_LayersAtTIBEnd) {
0706                 const DetLayer* tibLayer = barrelTIBLayers[misLayer - k_LayersStart - 1];
0707                 tmpTmeas = layerMeasurements.measurements(*tibLayer, tsos, prop, chi2Estimator);
0708               } else if (misLayer > k_LayersAtTIBEnd && misLayer < k_LayersAtTOBEnd) {
0709                 const DetLayer* tobLayer = barrelTOBLayers[misLayer - k_LayersAtTIBEnd - 1];
0710                 tmpTmeas = layerMeasurements.measurements(*tobLayer, tsos, prop, chi2Estimator);
0711               }
0712             }
0713           }
0714           if ((int(TKlayers) > k_LayersStart && int(TKlayers) <= k_LayersAtTIBEnd) && int(prev_TKlayers) == 12) {
0715             const DetLayer* detlayer = itm->layer();
0716             const LayerMeasurements layerMeasurements{measTracker, *measurementTrackerEvent};
0717             const TrajectoryStateOnSurface tsos = itm->updatedState();
0718             std::vector<DetLayer::DetWithState> compatDets = detlayer->compatibleDets(tsos, prop, chi2Estimator);
0719             std::vector<ForwardDetLayer const*> negTIDLayers = measTracker.geometricSearchTracker()->negTidLayers();
0720             std::vector<ForwardDetLayer const*> posTIDLayers = measTracker.geometricSearchTracker()->posTidLayers();
0721 
0722             const DetLayer* tidLayerneg = negTIDLayers[k_LayersStart];
0723             const DetLayer* tidLayerpos = posTIDLayers[k_LayersStart];
0724             if (tTopo->tidSide(iidd) == 1) {
0725               tmpTmeas = layerMeasurements.measurements(*tidLayerneg, tsos, prop, chi2Estimator);
0726             } else if (tTopo->tidSide(iidd) == 2) {
0727               tmpTmeas = layerMeasurements.measurements(*tidLayerpos, tsos, prop, chi2Estimator);
0728             }
0729           }
0730 
0731           if ((int(TKlayers) > k_LayersAtTIBEnd && int(TKlayers) <= k_LayersAtTOBEnd) && int(prev_TKlayers) == 15) {
0732             const DetLayer* detlayer = itm->layer();
0733             const LayerMeasurements layerMeasurements{measTracker, *measurementTrackerEvent};
0734             const TrajectoryStateOnSurface tsos = itm->updatedState();
0735             std::vector<DetLayer::DetWithState> compatDets = detlayer->compatibleDets(tsos, prop, chi2Estimator);
0736 
0737             std::vector<ForwardDetLayer const*> negTECLayers = measTracker.geometricSearchTracker()->negTecLayers();
0738             std::vector<ForwardDetLayer const*> posTECLayers = measTracker.geometricSearchTracker()->posTecLayers();
0739 
0740             const DetLayer* tecLayerneg = negTECLayers[k_LayersStart];
0741             const DetLayer* tecLayerpos = posTECLayers[k_LayersStart];
0742             if (tTopo->tecSide(iidd) == 1) {
0743               tmpTmeas = layerMeasurements.measurements(*tecLayerneg, tsos, prop, chi2Estimator);
0744             } else if (tTopo->tecSide(iidd) == 2) {
0745               tmpTmeas = layerMeasurements.measurements(*tecLayerpos, tsos, prop, chi2Estimator);
0746             }
0747           }
0748 
0749           //Test for two consecutive missing hits
0750           if (int(TKlayers) - int(prev_TKlayers) == -3) {
0751             foundConsMissingHits = true;
0752             const DetLayer* detlayer = itm->layer();
0753             const LayerMeasurements layerMeasurements{measTracker, *measurementTrackerEvent};
0754             const TrajectoryStateOnSurface tsos = itm->updatedState();
0755             std::vector<DetLayer::DetWithState> compatDets = detlayer->compatibleDets(tsos, prop, chi2Estimator);
0756 
0757             if (misLayer > k_LayersStart && misLayer <= k_LayersAtTOBEnd && previousMisLayer > k_LayersStart &&
0758                 previousMisLayer <= k_LayersAtTOBEnd) {  //Barrel case
0759               std::vector<BarrelDetLayer const*> barrelTIBLayers = measTracker.geometricSearchTracker()->tibLayers();
0760               std::vector<BarrelDetLayer const*> barrelTOBLayers = measTracker.geometricSearchTracker()->tobLayers();
0761               if (misLayer > k_LayersStart && misLayer < k_LayersAtTIBEnd) {
0762                 const DetLayer* tibLayer = barrelTIBLayers[misLayer - k_LayersStart - 1];
0763                 const DetLayer* prevTibLayer = barrelTIBLayers[previousMisLayer - k_LayersStart - 1];
0764 
0765                 tmpTmeas = layerMeasurements.measurements(*tibLayer, tsos, prop, chi2Estimator);
0766                 prev_tmpTmeas = layerMeasurements.measurements(*prevTibLayer, tsos, prop, chi2Estimator);
0767               } else if (misLayer > k_LayersAtTIBEnd && misLayer < k_LayersAtTOBEnd) {
0768                 const DetLayer* tobLayer = barrelTOBLayers[misLayer - k_LayersAtTIBEnd - 1];
0769                 const DetLayer* prevTobLayer = barrelTOBLayers[previousMisLayer - k_LayersAtTIBEnd - 1];
0770                 tmpTmeas = layerMeasurements.measurements(*tobLayer, tsos, prop, chi2Estimator);
0771                 prev_tmpTmeas = layerMeasurements.measurements(*prevTobLayer, tsos, prop, chi2Estimator);
0772               }
0773             } else if (misLayer > k_LayersAtTIDEnd && misLayer < k_LayersAtTECEnd &&
0774                        previousMisLayer > k_LayersAtTIDEnd && previousMisLayer < k_LayersAtTECEnd) {  //TEC
0775               std::vector<ForwardDetLayer const*> negTECLayers = measTracker.geometricSearchTracker()->negTecLayers();
0776               std::vector<ForwardDetLayer const*> posTECLayers = measTracker.geometricSearchTracker()->posTecLayers();
0777 
0778               const DetLayer* tecLayerneg = negTECLayers[misLayer - k_LayersAtTIDEnd - 1];
0779               const DetLayer* prevTecLayerneg = negTECLayers[previousMisLayer - k_LayersAtTIDEnd - 1];
0780 
0781               const DetLayer* tecLayerpos = posTECLayers[misLayer - k_LayersAtTIDEnd - 1];
0782               const DetLayer* prevTecLayerpos = posTECLayers[previousMisLayer - k_LayersAtTIDEnd - 1];
0783 
0784               if (tTopo->tecSide(iidd) == 1) {
0785                 tmpTmeas = layerMeasurements.measurements(*tecLayerneg, tsos, prop, chi2Estimator);
0786                 prev_tmpTmeas = layerMeasurements.measurements(*prevTecLayerneg, tsos, prop, chi2Estimator);
0787               } else if (tTopo->tecSide(iidd) == 2) {
0788                 tmpTmeas = layerMeasurements.measurements(*tecLayerpos, tsos, prop, chi2Estimator);
0789                 prev_tmpTmeas = layerMeasurements.measurements(*prevTecLayerpos, tsos, prop, chi2Estimator);
0790               }
0791             }
0792           }
0793           if (!tmpTmeas.empty() && !foundConsMissingHits) {
0794             TrajectoryMeasurement TM_tmp(tmpTmeas.back());
0795             unsigned int iidd_tmp = TM_tmp.recHit()->geographicalId().rawId();
0796             if (iidd_tmp != 0) {
0797               LogDebug("SiStripHitEfficiency:HitEff") << " hit actually being added to TM vector";
0798               if ((!useAllHitsFromTracksWithMissingHits_ || (!useFirstMeas_ && isFirstMeas)))
0799                 TMs.clear();
0800               if (::isDoubleSided(iidd_tmp, tTopo)) {
0801                 TMs.push_back(TrajectoryAtInvalidHit(TM_tmp, tTopo, tkgeom, propagator, 1));
0802                 TMs.push_back(TrajectoryAtInvalidHit(TM_tmp, tTopo, tkgeom, propagator, 2));
0803               } else
0804                 TMs.push_back(TrajectoryAtInvalidHit(TM_tmp, tTopo, tkgeom, propagator));
0805               missingHitAdded = true;
0806               hitRecoveryCounters[misLayer] += 1;
0807             }
0808           }
0809 
0810           if (!tmpTmeas.empty() && !prev_tmpTmeas.empty() &&
0811               foundConsMissingHits) {  //Found two consecutive missing hits
0812             TrajectoryMeasurement TM_tmp1(tmpTmeas.back());
0813             TrajectoryMeasurement TM_tmp2(prev_tmpTmeas.back());
0814             //Inner and outer hits module IDs
0815             unsigned int modIdInner = TM_tmp1.recHit()->geographicalId().rawId();
0816             unsigned int modIdOuter = TM_tmp2.recHit()->geographicalId().rawId();
0817             bool innerModInactive = false, outerModInactive = false;
0818             for (const auto& tm : tmpTmeas) {  //Check if inner module is inactive
0819               unsigned int tmModId = tm.recHit()->geographicalId().rawId();
0820               if (tmModId == modIdInner && tm.recHit()->getType() == 2) {
0821                 innerModInactive = true;
0822                 break;
0823               }
0824             }
0825             for (const auto& tm : prev_tmpTmeas) {  //Check if outer module is inactive
0826               unsigned int tmModId = tm.recHit()->geographicalId().rawId();
0827               if (tmModId == modIdOuter && tm.recHit()->getType() == 2) {
0828                 outerModInactive = true;
0829                 break;  //Found the inactive module
0830               }
0831             }
0832 
0833             if (outerModInactive) {  //If outer missing hit is in inactive module, recover the inner one
0834               if (modIdInner != 0) {
0835                 LogDebug("SiStripHitEfficiency:HitEff") << " hit actually being added to TM vector";
0836                 if ((!useAllHitsFromTracksWithMissingHits_ || (!useFirstMeas_ && isFirstMeas)))
0837                   TMs.clear();
0838                 if (::isDoubleSided(modIdInner, tTopo)) {
0839                   TMs.push_back(TrajectoryAtInvalidHit(TM_tmp1, tTopo, tkgeom, propagator, 1));
0840                   TMs.push_back(TrajectoryAtInvalidHit(TM_tmp1, tTopo, tkgeom, propagator, 2));
0841                 } else
0842                   TMs.push_back(TrajectoryAtInvalidHit(TM_tmp1, tTopo, tkgeom, propagator));
0843                 missingHitAdded = true;
0844                 hitRecoveryCounters[misLayer] += 1;
0845               }
0846             }
0847             if (innerModInactive) {  //If inner missing hit is in inactive module, recover the outer one
0848               if (modIdOuter != 0) {
0849                 LogDebug("SiStripHitEfficiency:HitEff") << " hit actually being added to TM vector";
0850                 if ((!useAllHitsFromTracksWithMissingHits_ || (!useFirstMeas_ && isFirstMeas)))
0851                   TMs.clear();
0852                 if (::isDoubleSided(modIdOuter, tTopo)) {
0853                   TMs.push_back(TrajectoryAtInvalidHit(TM_tmp2, tTopo, tkgeom, propagator, 1));
0854                   TMs.push_back(TrajectoryAtInvalidHit(TM_tmp2, tTopo, tkgeom, propagator, 2));
0855                 } else
0856                   TMs.push_back(TrajectoryAtInvalidHit(TM_tmp2, tTopo, tkgeom, propagator));
0857                 missingHitAdded = true;
0858                 hitRecoveryCounters[previousMisLayer] += 1;
0859               }
0860             }
0861           }
0862         }
0863 
0864         prev_TKlayers = TKlayers;
0865         if (!useFirstMeas_ && isFirstMeas && !missingHitAdded)
0866           continue;
0867         if (!useLastMeas_ && isLastMeas)
0868           continue;
0869         bool hitsWithBias = false;
0870         for (auto ilayer : missedLayers) {
0871           if (ilayer < TKlayers)
0872             hitsWithBias = true;
0873         }
0874         if (hasMissingHits && theInHit->getType() != TrackingRecHit::Type::missing && !missingHitAdded &&
0875             hitsWithBias && !useAllHitsFromTracksWithMissingHits_) {
0876           continue;
0877         }
0878 
0879         //////////////////////////////////////////////
0880         //Now check for tracks at TOB6 and TEC9
0881 
0882         // to make sure we only propagate on the last TOB5 hit check the next entry isn't also in TOB5
0883         // to avoid bias, make sure the TOB5 hit is valid (an invalid hit on TOB5 could only exist with a valid hit on TOB6)
0884         const auto nextId = (itm + 1 != TMeas.end()) ? (itm + 1)->recHit()->geographicalId() : DetId{};  // null if last
0885 
0886         if (TKlayers == 9 && theInHit->isValid() && !((!nextId.null()) && (::checkLayer(nextId.rawId(), tTopo) == 9))) {
0887           //      if ( TKlayers==9 && itm==TMeas.rbegin()) {
0888           //      if ( TKlayers==9 && (itm==TMeas.back()) ) {     // to check for only the last entry in the trajectory for propagation
0889           const DetLayer* tob6 = measTracker.geometricSearchTracker()->tobLayers().back();
0890           const LayerMeasurements theLayerMeasurements{measTracker, *measurementTrackerEvent};
0891           const TrajectoryStateOnSurface tsosTOB5 = itm->updatedState();
0892           const auto tmp = theLayerMeasurements.measurements(*tob6, tsosTOB5, prop, chi2Estimator);
0893 
0894           if (!tmp.empty()) {
0895             LogDebug("SiStripHitEfficiencyWorker") << "size of TM from propagation = " << tmp.size();
0896 
0897             // take the last of the TMs, which is always an invalid hit
0898             // if no detId is available, ie detId==0, then no compatible layer was crossed
0899             // otherwise, use that TM for the efficiency measurement
0900             const auto& tob6TM = tmp.back();
0901             const auto& tob6Hit = tob6TM.recHit();
0902             if (tob6Hit->geographicalId().rawId() != 0) {
0903               LogDebug("SiStripHitEfficiencyWorker") << "tob6 hit actually being added to TM vector";
0904               TMs.emplace_back(tob6TM, tTopo, tkgeom, propagator);
0905             }
0906           }
0907         }
0908 
0909         // same for TEC8
0910         if (TKlayers == 21 && theInHit->isValid() &&
0911             !((!nextId.null()) && (::checkLayer(nextId.rawId(), tTopo) == 21))) {
0912           const DetLayer* tec9pos = measTracker.geometricSearchTracker()->posTecLayers().back();
0913           const DetLayer* tec9neg = measTracker.geometricSearchTracker()->negTecLayers().back();
0914 
0915           const LayerMeasurements theLayerMeasurements{measTracker, *measurementTrackerEvent};
0916           const TrajectoryStateOnSurface tsosTEC9 = itm->updatedState();
0917 
0918           // check if track on positive or negative z
0919           if (!(iidd == SiStripSubdetector::TEC))
0920             LogDebug("SiStripHitEfficiencyWorker") << "there is a problem with TEC 9 extrapolation";
0921 
0922           //LogDebug("SiStripHitEfficiencyWorker") << " tec9 id = " << iidd << " and side = " << tTopo->tecSide(iidd) ;
0923           std::vector<TrajectoryMeasurement> tmp;
0924           if (tTopo->tecSide(iidd) == 1) {
0925             tmp = theLayerMeasurements.measurements(*tec9neg, tsosTEC9, prop, chi2Estimator);
0926             //LogDebug("SiStripHitEfficiencyWorker") << "on negative side" ;
0927           }
0928           if (tTopo->tecSide(iidd) == 2) {
0929             tmp = theLayerMeasurements.measurements(*tec9pos, tsosTEC9, prop, chi2Estimator);
0930             //LogDebug("SiStripHitEfficiencyWorker") << "on positive side" ;
0931           }
0932 
0933           if (!tmp.empty()) {
0934             // take the last of the TMs, which is always an invalid hit
0935             // if no detId is available, ie detId==0, then no compatible layer was crossed
0936             // otherwise, use that TM for the efficiency measurement
0937             const auto& tec9TM = tmp.back();
0938             const auto& tec9Hit = tec9TM.recHit();
0939 
0940             const unsigned int tec9id = tec9Hit->geographicalId().rawId();
0941             LogDebug("SiStripHitEfficiencyWorker")
0942                 << "tec9id = " << tec9id << " is Double sided = " << ::isDoubleSided(tec9id, tTopo)
0943                 << "  and 0x3 = " << (tec9id & 0x3);
0944 
0945             if (tec9Hit->geographicalId().rawId() != 0) {
0946               LogDebug("SiStripHitEfficiencyWorker") << "tec9 hit actually being added to TM vector";
0947               // in tec the hit can be single or doubled sided. whenever the invalid hit at the end of vector of TMs is
0948               // double sided it is always on the matched surface, so we need to split it into the true sensor surfaces
0949               if (::isDoubleSided(tec9id, tTopo)) {
0950                 TMs.emplace_back(tec9TM, tTopo, tkgeom, propagator, 1);
0951                 TMs.emplace_back(tec9TM, tTopo, tkgeom, propagator, 2);
0952               } else
0953                 TMs.emplace_back(tec9TM, tTopo, tkgeom, propagator);
0954             }
0955           }  //else LogDebug("SiStripHitEfficiencyWorker") << "tec9 tmp empty" ;
0956         }
0957 
0958         for (const auto& tm : TMs) {
0959           fillForTraj(tm,
0960                       tTopo,
0961                       tkgeom,
0962                       stripcpe,
0963                       stripQuality,
0964                       fedErrorIds,
0965                       commonModeDigis,
0966                       *theClusters,
0967                       e.bunchCrossing(),
0968                       instLumi,
0969                       PU,
0970                       highPurity);
0971         }
0972         LogDebug("SiStripHitEfficiencyWorker") << "After looping over TrajAtValidHit list";
0973       }
0974       LogDebug("SiStripHitEfficiencyWorker") << "end TMeasurement loop";
0975     }
0976     LogDebug("SiStripHitEfficiencyWorker") << "end of trajectories loop";
0977   }
0978 }
0979 
0980 void SiStripHitEfficiencyWorker::fillForTraj(const TrajectoryAtInvalidHit& tm,
0981                                              const TrackerTopology* tTopo,
0982                                              const TrackerGeometry* tkgeom,
0983                                              const StripClusterParameterEstimator& stripCPE,
0984                                              const SiStripQuality& stripQuality,
0985                                              const DetIdVector& fedErrorIds,
0986                                              const edm::Handle<edm::DetSetVector<SiStripRawDigi>>& commonModeDigis,
0987                                              const edmNew::DetSetVector<SiStripCluster>& theClusters,
0988                                              int bunchCrossing,
0989                                              float instLumi,
0990                                              float PU,
0991                                              bool highPurity) {
0992   // --> Get trajectory from combinatedStat& e
0993   const auto iidd = tm.monodet_id();
0994   LogDebug("SiStripHitEfficiencyWorker") << "setting iidd = " << iidd << " before checking efficiency and ";
0995 
0996   const auto xloc = tm.localX();
0997   const auto yloc = tm.localY();
0998 
0999   const auto xErr = tm.localErrorX();
1000   const auto yErr = tm.localErrorY();
1001 
1002   int TrajStrip = -1;
1003 
1004   // reget layer from iidd here, to account for TOB 6 and TEC 9 TKlayers being off
1005   const auto TKlayers = ::checkLayer(iidd, tTopo);
1006 
1007   const bool withinAcceptance =
1008       tm.withinAcceptance() && (!::isInBondingExclusionZone(iidd, TKlayers, yloc, yErr, tTopo));
1009 
1010   if (  // (TKlayers > 0) && // FIXME confirm this
1011       ((layers_ == TKlayers) ||
1012        (layers_ == bounds::k_LayersStart))) {  // Look at the layer not used to reconstruct the track
1013     LogDebug("SiStripHitEfficiencyWorker") << "Looking at layer under study";
1014     unsigned int ModIsBad = 2;
1015     unsigned int SiStripQualBad = 0;
1016     float commonMode = -100;
1017 
1018     // RPhi RecHit Efficiency
1019 
1020     if (!theClusters.empty()) {
1021       LogDebug("SiStripHitEfficiencyWorker") << "Checking clusters with size = " << theClusters.size();
1022       std::vector<::ClusterInfo> VCluster_info;  //fill with X residual, X residual pull, local X
1023       const auto idsv = theClusters.find(iidd);
1024       if (idsv != theClusters.end()) {
1025         //if (DEBUG_)      LogDebug("SiStripHitEfficiencyWorker") << "the ID from the dsv = " << dsv.id() ;
1026         LogDebug("SiStripHitEfficiencyWorker")
1027             << "found  (ClusterId == iidd) with ClusterId = " << idsv->id() << " and iidd = " << iidd;
1028         const auto stripdet = dynamic_cast<const StripGeomDetUnit*>(tkgeom->idToDetUnit(DetId(iidd)));
1029         const StripTopology& Topo = stripdet->specificTopology();
1030 
1031         float hbedge = 0.0;
1032         float htedge = 0.0;
1033         float hapoth = 0.0;
1034         float uylfac = 0.0;
1035         float uxlden = 0.0;
1036         if (TKlayers > bounds::k_LayersAtTOBEnd) {
1037           const BoundPlane& plane = stripdet->surface();
1038           const TrapezoidalPlaneBounds* trapezoidalBounds(
1039               dynamic_cast<const TrapezoidalPlaneBounds*>(&(plane.bounds())));
1040           std::array<const float, 4> const& parameterTrap = (*trapezoidalBounds).parameters();  // el bueno aqui
1041           hbedge = parameterTrap[0];
1042           htedge = parameterTrap[1];
1043           hapoth = parameterTrap[3];
1044           uylfac = (htedge - hbedge) / (htedge + hbedge) / hapoth;
1045           uxlden = 1 + yloc * uylfac;
1046         }
1047 
1048         // Need to know position of trajectory in strip number for selecting the right APV later
1049         if (TrajStrip == -1) {
1050           int nstrips = Topo.nstrips();
1051           float pitch = stripdet->surface().bounds().width() / nstrips;
1052           TrajStrip = xloc / pitch + nstrips / 2.0;
1053           // Need additionnal corrections for endcap
1054           if (TKlayers > bounds::k_LayersAtTOBEnd) {
1055             const float TrajLocXMid = xloc / (1 + (htedge - hbedge) * yloc / (htedge + hbedge) /
1056                                                       hapoth);  // radialy extrapolated x loc position at middle
1057             TrajStrip = TrajLocXMid / pitch + nstrips / 2.0;
1058           }
1059           //LogDebug("SiStripHitEfficiency")<<" Layer "<<TKlayers<<" TrajStrip: "<<nstrips<<" "<<pitch<<" "<<TrajStrip;;
1060         }
1061 
1062         for (const auto& clus : *idsv) {
1063           StripClusterParameterEstimator::LocalValues parameters = stripCPE.localParameters(clus, *stripdet);
1064           float res = (parameters.first.x() - xloc);
1065           float sigma = ::checkConsistency(parameters, xloc, xErr);
1066           // The consistency is probably more accurately measured with the Chi2MeasurementEstimator. To use it
1067           // you need a TransientTrackingRecHit instead of the cluster
1068           //theEstimator=       new Chi2MeasurementEstimator(30);
1069           //const Chi2MeasurementEstimator *theEstimator(100);
1070           //theEstimator->estimate(tm.tsos(), TransientTrackingRecHit);
1071 
1072           if (TKlayers > bounds::k_LayersAtTOBEnd) {
1073             res = parameters.first.x() - xloc / uxlden;  // radialy extrapolated x loc position at middle
1074             sigma = abs(res) / sqrt(parameters.second.xx() + xErr * xErr / uxlden / uxlden +
1075                                     yErr * yErr * xloc * xloc * uylfac * uylfac / uxlden / uxlden / uxlden / uxlden);
1076           }
1077 
1078           VCluster_info.emplace_back(res, sigma, parameters.first.x());
1079 
1080           LogDebug("SiStripHitEfficiencyWorker") << "Have ID match. residual = " << res << "  res sigma = " << sigma;
1081           //LogDebug("SiStripHitEfficiencyWorker")
1082           //    << "trajectory measurement compatability estimate = " << (*itm).estimate() ;
1083           LogDebug("SiStripHitEfficiencyWorker")
1084               << "hit position = " << parameters.first.x() << "  hit error = " << sqrt(parameters.second.xx())
1085               << "  trajectory position = " << xloc << "  traj error = " << xErr;
1086         }
1087       }
1088       ::ClusterInfo finalCluster{1000.0, 1000.0, 0.0};
1089       if (!VCluster_info.empty()) {
1090         LogDebug("SiStripHitEfficiencyWorker") << "found clusters > 0";
1091         if (VCluster_info.size() > 1) {
1092           //get the smallest one
1093           for (const auto& res : VCluster_info) {
1094             if (std::abs(res.xResidualPull) < std::abs(finalCluster.xResidualPull)) {
1095               finalCluster = res;
1096             }
1097             LogDebug("SiStripHitEfficiencyWorker")
1098                 << "iresidual = " << res.xResidual << "  isigma = " << res.xResidualPull
1099                 << "  and FinalRes = " << finalCluster.xResidual;
1100           }
1101         } else {
1102           finalCluster = VCluster_info[0];
1103         }
1104         VCluster_info.clear();
1105       }
1106 
1107       LogDebug("SiStripHitEfficiencyWorker") << "Final residual in X = " << finalCluster.xResidual << "+-"
1108                                              << (finalCluster.xResidual / finalCluster.xResidualPull);
1109       LogDebug("SiStripHitEfficiencyWorker")
1110           << "Checking location of trajectory: abs(yloc) = " << abs(yloc) << "  abs(xloc) = " << abs(xloc);
1111 
1112       //
1113       // fill ntuple varibles
1114 
1115       //if ( stripQuality->IsModuleBad(iidd) )
1116       if (stripQuality.getBadApvs(iidd) != 0) {
1117         SiStripQualBad = 1;
1118         LogDebug("SiStripHitEfficiencyWorker") << "strip is bad from SiStripQuality";
1119       } else {
1120         SiStripQualBad = 0;
1121         LogDebug("SiStripHitEfficiencyWorker") << "strip is good from SiStripQuality";
1122       }
1123 
1124       //check for FED-detected errors and include those in SiStripQualBad
1125       for (unsigned int ii = 0; ii < fedErrorIds.size(); ii++) {
1126         if (iidd == fedErrorIds[ii].rawId())
1127           SiStripQualBad = 1;
1128       }
1129 
1130       // CM of APV crossed by traj
1131       if (addCommonMode_)
1132         if (commonModeDigis.isValid() && TrajStrip >= 0 && TrajStrip <= 768) {
1133           const auto digiframe = commonModeDigis->find(iidd);
1134           if (digiframe != commonModeDigis->end())
1135             if ((unsigned)TrajStrip / sistrip::STRIPS_PER_APV < digiframe->data.size())
1136               commonMode = digiframe->data.at(TrajStrip / sistrip::STRIPS_PER_APV).adc();
1137         }
1138 
1139       LogDebug("SiStripHitEfficiencyWorker") << "before check good";
1140 
1141       if (finalCluster.xResidualPull < 999.0) {  //could make requirement on track/hit consistency, but for
1142         //now take anything with a hit on the module
1143         LogDebug("SiStripHitEfficiencyWorker")
1144             << "hit being counted as good " << finalCluster.xResidual << " FinalRecHit " << iidd << "   TKlayers  "
1145             << TKlayers << " xloc " << xloc << " yloc  " << yloc << " module " << iidd
1146             << "   matched/stereo/rphi = " << ((iidd & 0x3) == 0) << "/" << ((iidd & 0x3) == 1) << "/"
1147             << ((iidd & 0x3) == 2);
1148         ModIsBad = 0;
1149       } else {
1150         LogDebug("SiStripHitEfficiencyWorker")
1151             << "hit being counted as bad   ######### Invalid RPhi FinalResX " << finalCluster.xResidual
1152             << " FinalRecHit " << iidd << "   TKlayers  " << TKlayers << " xloc " << xloc << " yloc  " << yloc
1153             << " module " << iidd << "   matched/stereo/rphi = " << ((iidd & 0x3) == 0) << "/" << ((iidd & 0x3) == 1)
1154             << "/" << ((iidd & 0x3) == 2);
1155         ModIsBad = 1;
1156         LogDebug("SiStripHitEfficiencyWorker")
1157             << " RPhi Error " << sqrt(xErr * xErr + yErr * yErr) << " ErrorX " << xErr << " yErr " << yErr;
1158       }
1159 
1160       LogDebug("SiStripHitEfficiencyWorker")
1161           << "To avoid them staying unused: ModIsBad=" << ModIsBad << ", SiStripQualBad=" << SiStripQualBad
1162           << ", commonMode=" << commonMode << ", highPurity=" << highPurity
1163           << ", withinAcceptance=" << withinAcceptance;
1164 
1165       unsigned int layer = TKlayers;
1166       if (showRings_ && layer > bounds::k_LayersAtTOBEnd) {  // use rings instead of wheels
1167         if (layer <= bounds::k_LayersAtTIDEnd) {             // TID
1168           layer = bounds::k_LayersAtTOBEnd +
1169                   tTopo->tidRing(iidd);  // ((iidd >> 9) & 0x3);  // 3 disks and also 3 rings -> use the same container
1170         } else {                         // TEC
1171           layer = bounds::k_LayersAtTIDEnd + tTopo->tecRing(iidd);  // ((iidd >> 5) & 0x7);
1172         }
1173       }
1174       unsigned int layerWithSide = layer;
1175       if (layer > bounds::k_LayersAtTOBEnd && layer <= bounds::k_LayersAtTIDEnd) {
1176         const auto side = tTopo->tidSide(iidd);  //(iidd >> 13) & 0x3;  // TID
1177         if (side == 2)
1178           layerWithSide = layer + 3;
1179       } else if (layer > bounds::k_LayersAtTIDEnd) {
1180         const auto side = tTopo->tecSide(iidd);  // (iidd >> 18) & 0x3;  // TEC
1181         if (side == 1) {
1182           layerWithSide = layer + 3;
1183         } else if (side == 2) {
1184           layerWithSide = layer + 3 + (showRings_ ? 7 : 9);
1185         }
1186       }
1187 
1188       if ((bunchX_ > 0 && bunchX_ != bunchCrossing) || (!withinAcceptance) ||
1189           (useOnlyHighPurityTracks_ && !highPurity) ||
1190           (!showTOB6TEC9_ && (TKlayers == bounds::k_LayersAtTOBEnd || TKlayers == bounds::k_LayersAtTECEnd)) ||
1191           (badModules_.end() != badModules_.find(iidd)))
1192         return;
1193 
1194       const bool badquality = (SiStripQualBad == 1);
1195 
1196       //Now that we have a good event, we need to look at if we expected it or not, and the location
1197       //if we didn't
1198       //Fill the missing hit information first
1199       bool badflag = false;  // true for hits that are expected but not found
1200       if (resXSig_ < 0) {
1201         if (ModIsBad == 1)
1202           badflag = true;  // isBad set to false in the tree when resxsig<999.0
1203       } else {
1204         if (ModIsBad == 1 || finalCluster.xResidualPull > resXSig_)
1205           badflag = true;
1206       }
1207 
1208       // Conversion of positions in strip unit
1209       int nstrips = -9;
1210       float Pitch = -9.0;
1211       const StripGeomDetUnit* stripdet = nullptr;
1212       if (finalCluster.xResidualPull ==
1213           1000.0) {      // special treatment, no GeomDetUnit associated in some cases when no cluster found
1214         Pitch = 0.0205;  // maximum
1215         nstrips = 768;   // maximum
1216       } else {
1217         stripdet = dynamic_cast<const StripGeomDetUnit*>(tkgeom->idToDetUnit(iidd));
1218         const StripTopology& Topo = stripdet->specificTopology();
1219         nstrips = Topo.nstrips();
1220         Pitch = stripdet->surface().bounds().width() / Topo.nstrips();
1221       }
1222       double stripTrajMid = xloc / Pitch + nstrips / 2.0;
1223       double stripCluster = finalCluster.xLocal / Pitch + nstrips / 2.0;
1224       // For trapezoidal modules: extrapolation of x trajectory position to the y middle of the module
1225       //  for correct comparison with cluster position
1226       if (stripdet && layer > bounds::k_LayersAtTOBEnd) {
1227         const auto& trapezoidalBounds = dynamic_cast<const TrapezoidalPlaneBounds&>(stripdet->surface().bounds());
1228         std::array<const float, 4> const& parameters = trapezoidalBounds.parameters();
1229         const float hbedge = parameters[0];
1230         const float htedge = parameters[1];
1231         const float hapoth = parameters[3];
1232         const float TrajLocXMid = xloc / (1 + (htedge - hbedge) * yloc / (htedge + hbedge) /
1233                                                   hapoth);  // radialy extrapolated x loc position at middle
1234         stripTrajMid = TrajLocXMid / Pitch + nstrips / 2.0;
1235       }
1236 
1237       if ((!badquality) && (layer < h_resolution.size())) {
1238         LogDebug("SiStripHitEfficiencyWorker")
1239             << "layer " << layer << " vector index " << layer - 1 << " before filling h_resolution" << std::endl;
1240         h_resolution[layer - 1]->Fill(finalCluster.xResidualPull != 1000.0 ? stripTrajMid - stripCluster : 1000);
1241       }
1242 
1243       // New matching methods
1244       if (clusterMatchingMethod_ >= 1) {
1245         badflag = false;
1246         if (finalCluster.xResidualPull == 1000.0) {
1247           LogDebug("SiStripHitEfficiencyWorker") << "Marking bad for resxsig=1000";
1248           badflag = true;
1249         } else {
1250           if (clusterMatchingMethod_ == 2 || clusterMatchingMethod_ == 4) {
1251             // check the distance between cluster and trajectory position
1252             if (std::abs(stripCluster - stripTrajMid) > clusterTracjDist_) {
1253               LogDebug("SiStripHitEfficiencyWorker") << "Marking bad for cluster-to-traj distance";
1254               badflag = true;
1255             }
1256           }
1257           if (clusterMatchingMethod_ == 3 || clusterMatchingMethod_ == 4) {
1258             // cluster and traj have to be in the same APV (don't take edges into accounts)
1259             const int tapv = (int)stripTrajMid / sistrip::STRIPS_PER_APV;
1260             const int capv = (int)stripCluster / sistrip::STRIPS_PER_APV;
1261             float stripInAPV = stripTrajMid - tapv * sistrip::STRIPS_PER_APV;
1262             if (stripInAPV < stripsApvEdge_ || stripInAPV > sistrip::STRIPS_PER_APV - stripsApvEdge_) {
1263               LogDebug("SiStripHitEfficiencyWorker") << "Too close to the edge: " << stripInAPV;
1264               return;
1265             }
1266             if (tapv != capv) {
1267               LogDebug("SiStripHitEfficiencyWorker") << "Marking bad for tapv!=capv";
1268               badflag = true;
1269             }
1270           }
1271         }
1272       }
1273       if (!badquality) {
1274         LogDebug("SiStripHitEfficiencyWorker")
1275             << "Filling measurement for " << iidd << " in layer " << layer << " histograms with bx=" << bunchCrossing
1276             << ", lumi=" << instLumi << ", PU=" << PU << "; bad flag=" << badflag;
1277 
1278         // hot/cold maps of hits that are expected but not found
1279         if (badflag) {
1280           if (layer > bounds::k_LayersStart && layer <= bounds::k_LayersAtTIBEnd) {
1281             //We are in the TIB
1282             float phi = ::calcPhi(tm.globalX(), tm.globalY());
1283             h_hotcold[layer - 1]->Fill(360. - phi, tm.globalZ(), 1.);
1284           } else if (layer > bounds::k_LayersAtTIBEnd && layer <= bounds::k_LayersAtTOBEnd) {
1285             //We are in the TOB
1286             float phi = ::calcPhi(tm.globalX(), tm.globalY());
1287             h_hotcold[layer - 1]->Fill(360. - phi, tm.globalZ(), 1.);
1288           } else if (layer > bounds::k_LayersAtTOBEnd && layer <= bounds::k_LayersAtTIDEnd) {
1289             //We are in the TID
1290             //There are 2 different maps here
1291             int side = tTopo->tidSide(iidd);
1292             if (side == 1)
1293               h_hotcold[(layer - 1) + (layer - 11)]->Fill(-tm.globalY(), tm.globalX(), 1.);
1294             else if (side == 2)
1295               h_hotcold[(layer - 1) + (layer - 10)]->Fill(-tm.globalY(), tm.globalX(), 1.);
1296           } else if (layer > bounds::k_LayersAtTIDEnd) {
1297             //We are in the TEC
1298             //There are 2 different maps here
1299             int side = tTopo->tecSide(iidd);
1300             if (side == 1)
1301               h_hotcold[(layer + 2) + (layer - 14)]->Fill(-tm.globalY(), tm.globalX(), 1.);
1302             else if (side == 2)
1303               h_hotcold[(layer + 2) + (layer - 13)]->Fill(-tm.globalY(), tm.globalX(), 1.);
1304           }
1305         }
1306 
1307         LogDebug("SiStripHitEfficiencyWorker")
1308             << "layer " << layer << " vector index " << layer - 1 << " before filling h_layer_vsSmthg" << std::endl;
1309         h_layer_vsBx[layer - 1].fill(bunchCrossing, !badflag);
1310         if (addLumi_) {
1311           h_layer_vsLumi[layer - 1].fill(instLumi, !badflag);
1312           h_layer_vsPU[layer - 1].fill(PU, !badflag);
1313         }
1314         if (addCommonMode_) {
1315           h_layer_vsCM[layer - 1].fill(commonMode, !badflag);
1316         }
1317         h_goodLayer.fill(layerWithSide, !badflag);
1318       }
1319       // efficiency without bad modules excluded
1320       h_allLayer.fill(layerWithSide, !badflag);
1321 
1322       // efficiency without bad modules excluded
1323       if (TKlayers) {
1324         h_module.fill(iidd, !badflag);
1325         assert(h_module.check(iidd));
1326       }
1327 
1328       /* Used in SiStripHitEffFromCalibTree:
1329        * run              -> "run"              -> run              // e.id().run()
1330        * event            -> "event"            -> evt              // e.id().event()
1331        * ModIsBad         -> "ModIsBad"         -> isBad
1332        * SiStripQualBad   -> "SiStripQualBad""  -> quality
1333        * Id               -> "Id"               -> id               // iidd
1334        * withinAcceptance -> "withinAcceptance" -> accept
1335        * whatlayer        -> "layer"            -> layer_wheel      // Tklayers
1336        * highPurity       -> "highPurity"       -> highPurity
1337        * TrajGlbX         -> "TrajGlbX"         -> x                // tm.globalX()
1338        * TrajGlbY         -> "TrajGlbY"         -> y                // tm.globalY()
1339        * TrajGlbZ         -> "TrajGlbZ"         -> z                // tm.globalZ()
1340        * ResXSig          -> "ResXSig"          -> resxsig          // finalCluster.xResidualPull;
1341        * TrajLocX         -> "TrajLocX"         -> TrajLocX         // xloc
1342        * TrajLocY         -> "TrajLocY"         -> TrajLocY         // yloc
1343        * ClusterLocX      -> "ClusterLocX"      -> ClusterLocX      // finalCluster.xLocal
1344        * bunchx           -> "bunchx"           -> bx               // e.bunchCrossing()
1345        * instLumi         -> "instLumi"         -> instLumi         ## if addLumi_
1346        * PU               -> "PU"               -> PU               ## if addLumi_
1347        * commonMode       -> "commonMode"       -> CM               ## if addCommonMode_ / _useCM
1348        */
1349       LogDebug("SiStripHitEfficiencyWorker") << "after good location check";
1350     }
1351     LogDebug("SiStripHitEfficiencyWorker") << "after list of clusters";
1352   }
1353   LogDebug("SiStripHitEfficiencyWorker") << "After layers=TKLayers if with TKlayers=" << TKlayers
1354                                          << ", layers=" << layers_;
1355 }
1356 
1357 void SiStripHitEfficiencyWorker::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
1358   edm::ParameterSetDescription desc;
1359   desc.add<std::string>("dqmDir", "AlCaReco/SiStripHitEfficiency");
1360   desc.add<bool>("UseOnlyHighPurityTracks", true);
1361   desc.add<bool>("cutOnTracks", false);
1362   desc.add<bool>("doMissingHitsRecovery", false);
1363   desc.add<bool>("useAllHitsFromTracksWithMissingHits", false);
1364   desc.add<bool>("useFirstMeas", false);
1365   desc.add<bool>("useLastMeas", false);
1366   desc.add<double>("ClusterTrajDist", 64.0);
1367   desc.add<double>("ResXSig", -1);
1368   desc.add<double>("StripsApvEdge", 10.0);
1369   desc.add<edm::InputTag>("combinatorialTracks", edm::InputTag{"generalTracks"});
1370   desc.add<edm::InputTag>("commonMode", edm::InputTag{"siStripDigis", "CommonMode"});
1371   desc.add<edm::InputTag>("lumiScalers", edm::InputTag{"scalersRawToDigi"});
1372   desc.add<edm::InputTag>("metadata", edm::InputTag{"onlineMetaDataDigis"});
1373   desc.add<edm::InputTag>("siStripClusters", edm::InputTag{"siStripClusters"});
1374   desc.add<edm::InputTag>("siStripDigis", edm::InputTag{"siStripDigis"});
1375   desc.add<edm::InputTag>("trackerEvent", edm::InputTag{"MeasurementTrackerEvent"});
1376   desc.add<edm::InputTag>("trajectories", edm::InputTag{"generalTracks"});
1377   desc.add<int>("ClusterMatchingMethod", 0);
1378   desc.add<int>("Layer", 0);
1379   desc.add<unsigned int>("trackMultiplicity", 100);
1380   desc.addUntracked<bool>("Debug", false);
1381   desc.addUntracked<bool>("ShowRings", false);
1382   desc.addUntracked<bool>("ShowTOB6TEC9", false);
1383   desc.addUntracked<bool>("addCommonMode", false);
1384   desc.addUntracked<bool>("addLumi", true);
1385   desc.addUntracked<int>("BunchCrossing", 0);
1386   desc.addUntracked<std::string>("BadModulesFile", "");
1387   descriptions.addWithDefaultLabel(desc);
1388 }
1389 
1390 #include "FWCore/Framework/interface/MakerMacros.h"
1391 DEFINE_FWK_MODULE(SiStripHitEfficiencyWorker);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!