Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-06-16 22:18:32

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/SiStripHitEfficiencyHelpers.h"
0013 #include "CalibTracker/SiStripHitEfficiency/interface/TrajectoryAtInvalidHit.h"
0014 #include "DQM/SiStripCommon/interface/TkHistoMap.h"
0015 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0016 #include "DataFormats/Common/interface/DetSetVector.h"
0017 #include "DataFormats/Common/interface/DetSetVectorNew.h"
0018 #include "DataFormats/Common/interface/Handle.h"
0019 #include "DataFormats/DetId/interface/DetIdCollection.h"
0020 #include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementError.h"
0021 #include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementVector.h"
0022 #include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h"
0023 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0024 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0025 #include "DataFormats/GeometryVector/interface/LocalVector.h"
0026 #include "DataFormats/OnlineMetaData/interface/OnlineLuminosityRecord.h"
0027 #include "DataFormats/Scalers/interface/LumiScalers.h"
0028 #include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
0029 #include "DataFormats/SiStripCommon/interface/ConstantsForHardwareSystems.h" /* for STRIPS_PER_APV*/
0030 #include "DataFormats/SiStripDigi/interface/SiStripRawDigi.h"
0031 #include "DataFormats/TrackReco/interface/Track.h"
0032 #include "DataFormats/TrackReco/interface/TrackBase.h"
0033 #include "DataFormats/TrackReco/interface/TrackExtra.h"
0034 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0035 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0036 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0037 #include "FWCore/Framework/interface/Event.h"
0038 #include "FWCore/Framework/interface/EventSetup.h"
0039 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0040 #include "FWCore/ParameterSet/interface/ParameterDescription.h"
0041 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0042 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0043 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0044 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0045 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0046 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0047 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0048 #include "MagneticField/Engine/interface/MagneticField.h"
0049 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0050 #include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h"
0051 #include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h"
0052 #include "RecoTracker/Record/interface/CkfComponentsRecord.h"
0053 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0054 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
0055 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
0056 #include "TrackingTools/MaterialEffects/interface/PropagatorWithMaterial.h"
0057 #include "TrackingTools/MeasurementDet/interface/LayerMeasurements.h"
0058 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0059 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0060 
0061 class SiStripHitEfficiencyWorker : public DQMEDAnalyzer {
0062 public:
0063   explicit SiStripHitEfficiencyWorker(const edm::ParameterSet& conf);
0064   ~SiStripHitEfficiencyWorker() override = default;
0065   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0066 
0067 private:
0068   void beginJob();  // TODO remove
0069   void endJob();    // TODO remove
0070   void bookHistograms(DQMStore::IBooker& booker, const edm::Run& run, const edm::EventSetup& setup) override;
0071   void analyze(const edm::Event& e, const edm::EventSetup& c) override;
0072   void fillForTraj(const TrajectoryAtInvalidHit& tm,
0073                    const TrackerTopology* tTopo,
0074                    const TrackerGeometry* tkgeom,
0075                    const StripClusterParameterEstimator& stripCPE,
0076                    const SiStripQuality& stripQuality,
0077                    const DetIdCollection& fedErrorIds,
0078                    const edm::Handle<edm::DetSetVector<SiStripRawDigi>>& commonModeDigis,
0079                    const edmNew::DetSetVector<SiStripCluster>& theClusters,
0080                    int bunchCrossing,
0081                    float instLumi,
0082                    float PU,
0083                    bool highPurity);
0084 
0085   // ----------member data ---------------------------
0086 
0087   // event data tokens
0088   const edm::EDGetTokenT<LumiScalersCollection> scalerToken_;
0089   const edm::EDGetTokenT<OnlineLuminosityRecord> metaDataToken_;
0090   const edm::EDGetTokenT<edm::DetSetVector<SiStripRawDigi>> commonModeToken_;
0091   const edm::EDGetTokenT<reco::TrackCollection> combinatorialTracks_token_;
0092   const edm::EDGetTokenT<std::vector<Trajectory>> trajectories_token_;
0093   const edm::EDGetTokenT<TrajTrackAssociationCollection> trajTrackAsso_token_;
0094   const edm::EDGetTokenT<edmNew::DetSetVector<SiStripCluster>> clusters_token_;
0095   const edm::EDGetTokenT<DetIdCollection> digis_token_;
0096   const edm::EDGetTokenT<MeasurementTrackerEvent> trackerEvent_token_;
0097 
0098   // event setup tokens
0099   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> tTopoToken_;
0100   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> tkGeomToken_;
0101   const edm::ESGetToken<StripClusterParameterEstimator, TkStripCPERecord> stripCPEToken_;
0102   const edm::ESGetToken<SiStripQuality, SiStripQualityRcd> stripQualityToken_;
0103   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magFieldToken_;
0104   const edm::ESGetToken<MeasurementTracker, CkfComponentsRecord> measTrackerToken_;
0105   const edm::ESGetToken<Chi2MeasurementEstimatorBase, TrackingComponentsRecord> chi2EstimatorToken_;
0106   const edm::ESGetToken<Propagator, TrackingComponentsRecord> propagatorToken_;
0107   const edm::ESGetToken<TkDetMap, TrackerTopologyRcd> tkDetMapToken_;
0108 
0109   // configurable parameters
0110   std::string dqmDir_;
0111   unsigned int layers_;
0112   bool DEBUG_;
0113   bool addLumi_;
0114   bool addCommonMode_;
0115   bool cutOnTracks_;
0116   unsigned int trackMultiplicityCut_;
0117   bool useFirstMeas_;
0118   bool useLastMeas_;
0119   bool useAllHitsFromTracksWithMissingHits_;
0120   unsigned int clusterMatchingMethod_;
0121   float resXSig_;
0122   float clusterTracjDist_;
0123   float stripsApvEdge_;
0124   bool useOnlyHighPurityTracks_;
0125   int bunchX_;
0126   bool showRings_;
0127   bool showTOB6TEC9_;
0128   unsigned int nTEClayers_;
0129 
0130   // output file
0131   std::set<uint32_t> badModules_;
0132 
0133   // counters
0134   int events, EventTrackCKF;
0135 
0136   struct EffME1 {
0137     EffME1() : hTotal(nullptr), hFound(nullptr) {}
0138     EffME1(MonitorElement* total, MonitorElement* found) : hTotal(total), hFound(found) {}
0139 
0140     void fill(double x, bool found, float weight = 1.) {
0141       hTotal->Fill(x, weight);
0142       if (found) {
0143         hFound->Fill(x, weight);
0144       }
0145     }
0146 
0147     MonitorElement *hTotal, *hFound;
0148   };
0149   struct EffTkMap {
0150     EffTkMap() : hTotal(nullptr), hFound(nullptr) {}
0151     EffTkMap(std::unique_ptr<TkHistoMap>&& total, std::unique_ptr<TkHistoMap>&& found)
0152         : hTotal(std::move(total)), hFound(std::move(found)) {}
0153 
0154     void fill(uint32_t id, bool found, float weight = 1.) {
0155       hTotal->fill(id, weight);
0156       if (found) {
0157         hFound->fill(id, weight);
0158       }
0159     }
0160 
0161     std::unique_ptr<TkHistoMap> hTotal, hFound;
0162   };
0163 
0164   MonitorElement *h_bx, *h_instLumi, *h_PU;
0165   MonitorElement *h_nTracks, *h_nTracksVsPU;
0166   EffME1 h_goodLayer;
0167   EffME1 h_allLayer;
0168   EffME1 h_layer;
0169   std::vector<MonitorElement*> h_resolution;
0170   std::vector<EffME1> h_layer_vsLumi;
0171   std::vector<EffME1> h_layer_vsBx;
0172   std::vector<EffME1> h_layer_vsPU;
0173   std::vector<EffME1> h_layer_vsCM;
0174   std::vector<MonitorElement*> h_hotcold;
0175 
0176   EffTkMap h_module;
0177 };
0178 
0179 //
0180 // constructors and destructor
0181 //
0182 
0183 SiStripHitEfficiencyWorker::SiStripHitEfficiencyWorker(const edm::ParameterSet& conf)
0184     : scalerToken_(consumes<LumiScalersCollection>(conf.getParameter<edm::InputTag>("lumiScalers"))),
0185       metaDataToken_(consumes<OnlineLuminosityRecord>(conf.getParameter<edm::InputTag>("metadata"))),
0186       commonModeToken_(mayConsume<edm::DetSetVector<SiStripRawDigi>>(conf.getParameter<edm::InputTag>("commonMode"))),
0187       combinatorialTracks_token_(
0188           consumes<reco::TrackCollection>(conf.getParameter<edm::InputTag>("combinatorialTracks"))),
0189       trajectories_token_(consumes<std::vector<Trajectory>>(conf.getParameter<edm::InputTag>("trajectories"))),
0190       trajTrackAsso_token_(consumes<TrajTrackAssociationCollection>(conf.getParameter<edm::InputTag>("trajectories"))),
0191       clusters_token_(
0192           consumes<edmNew::DetSetVector<SiStripCluster>>(conf.getParameter<edm::InputTag>("siStripClusters"))),
0193       digis_token_(consumes<DetIdCollection>(conf.getParameter<edm::InputTag>("siStripDigis"))),
0194       trackerEvent_token_(consumes<MeasurementTrackerEvent>(conf.getParameter<edm::InputTag>("trackerEvent"))),
0195       tTopoToken_(esConsumes()),
0196       tkGeomToken_(esConsumes()),
0197       stripCPEToken_(esConsumes(edm::ESInputTag{"", "StripCPEfromTrackAngle"})),
0198       stripQualityToken_(esConsumes()),
0199       magFieldToken_(esConsumes()),
0200       measTrackerToken_(esConsumes()),
0201       chi2EstimatorToken_(esConsumes(edm::ESInputTag{"", "Chi2"})),
0202       propagatorToken_(esConsumes(edm::ESInputTag{"", "PropagatorWithMaterial"})),
0203       tkDetMapToken_(esConsumes<edm::Transition::BeginRun>()),
0204       dqmDir_(conf.getParameter<std::string>("dqmDir")),
0205       layers_(conf.getParameter<int>("Layer")),
0206       DEBUG_(conf.getUntrackedParameter<bool>("Debug", false)),
0207       addLumi_(conf.getUntrackedParameter<bool>("addLumi", false)),
0208       addCommonMode_(conf.getUntrackedParameter<bool>("addCommonMode", false)),
0209       cutOnTracks_(conf.getParameter<bool>("cutOnTracks")),
0210       trackMultiplicityCut_(conf.getParameter<unsigned int>("trackMultiplicity")),
0211       useFirstMeas_(conf.getParameter<bool>("useFirstMeas")),
0212       useLastMeas_(conf.getParameter<bool>("useLastMeas")),
0213       useAllHitsFromTracksWithMissingHits_(conf.getParameter<bool>("useAllHitsFromTracksWithMissingHits")),
0214       clusterMatchingMethod_(conf.getParameter<int>("ClusterMatchingMethod")),
0215       resXSig_(conf.getParameter<double>("ResXSig")),
0216       clusterTracjDist_(conf.getParameter<double>("ClusterTrajDist")),
0217       stripsApvEdge_(conf.getParameter<double>("StripsApvEdge")),
0218       useOnlyHighPurityTracks_(conf.getParameter<bool>("UseOnlyHighPurityTracks")),
0219       bunchX_(conf.getUntrackedParameter<int>("BunchCrossing", 0)),
0220       showRings_(conf.getUntrackedParameter<bool>("ShowRings", false)),
0221       showTOB6TEC9_(conf.getUntrackedParameter<bool>("ShowTOB6TEC9", false)) {
0222   nTEClayers_ = (showRings_ ? 7 : 9);  // number of rings or wheels
0223 
0224   const std::string badModulesFile = conf.getUntrackedParameter<std::string>("BadModulesFile", "");
0225   if (!badModulesFile.empty()) {
0226     std::ifstream badModules_file(badModulesFile);
0227     uint32_t badmodule_detid;
0228     int mods, fiber1, fiber2, fiber3;
0229     if (badModules_file.is_open()) {
0230       std::string line;
0231       while (getline(badModules_file, line)) {
0232         if (badModules_file.eof())
0233           continue;
0234         std::stringstream ss(line);
0235         ss >> badmodule_detid >> mods >> fiber1 >> fiber2 >> fiber3;
0236         if (badmodule_detid != 0 && mods == 1 && (fiber1 == 1 || fiber2 == 1 || fiber3 == 1))
0237           badModules_.insert(badmodule_detid);
0238       }
0239       badModules_file.close();
0240     }
0241   }
0242   if (!badModules_.empty())
0243     LogDebug("SiStripHitEfficiencyWorker") << "Remove additionnal bad modules from the analysis: ";
0244   for (const auto badMod : badModules_) {
0245     LogDebug("SiStripHitEfficiencyWorker") << " " << badMod;
0246   }
0247 }
0248 
0249 void SiStripHitEfficiencyWorker::beginJob() {
0250   // TODO convert to counters, or simply remove?
0251   events = 0;
0252   EventTrackCKF = 0;
0253 }
0254 
0255 void SiStripHitEfficiencyWorker::bookHistograms(DQMStore::IBooker& booker,
0256                                                 const edm::Run& run,
0257                                                 const edm::EventSetup& setup) {
0258   booker.setCurrentFolder(fmt::format("{}/EventInfo", dqmDir_));
0259   h_bx = booker.book1D("bx", "bx", 3600, 0, 3600);
0260   h_instLumi = booker.book1D("instLumi", "inst. lumi.", 250, 0, 25000);
0261   h_PU = booker.book1D("PU", "PU", 200, 0, 200);
0262   h_nTracks = booker.book1D("ntracks", "n.tracks;n. tracks;n.events", 500, -0.5, 499.5);
0263   h_nTracksVsPU = booker.bookProfile("nTracksVsPU", "n. tracks vs PU; PU; n.tracks ", 200, 0, 200, 500, -0.5, 499.5);
0264 
0265   booker.setCurrentFolder(dqmDir_);
0266   h_goodLayer = EffME1(booker.book1D("goodlayer_total", "goodlayer_total", 35, 0., 35.),
0267                        booker.book1D("goodlayer_found", "goodlayer_found", 35, 0., 35.));
0268   h_allLayer = EffME1(booker.book1D("alllayer_total", "alllayer_total", 35, 0., 35.),
0269                       booker.book1D("alllayer_found", "alllayer_found", 35, 0., 35.));
0270 
0271   h_layer = EffME1(
0272       booker.book1D(
0273           "layer_found", "layer_found", bounds::k_END_OF_LAYERS, 0., static_cast<float>(bounds::k_END_OF_LAYERS)),
0274       booker.book1D(
0275           "layer_total", "layer_total", bounds::k_END_OF_LAYERS, 0., static_cast<float>(bounds::k_END_OF_LAYERS)));
0276 
0277   for (int layer = 1; layer != bounds::k_END_OF_LAYERS; ++layer) {
0278     const auto lyrName = ::layerName(layer, showRings_, nTEClayers_);
0279 
0280     // book resolutions
0281     booker.setCurrentFolder(fmt::format("{}/Resolutions", dqmDir_));
0282     auto ihres = booker.book1D(Form("resol_layer_%i", layer), lyrName, 125, -125., 125.);
0283     ihres->setAxisTitle("trajX-clusX [strip unit]");
0284     h_resolution.push_back(ihres);
0285 
0286     // book plots vs Lumi
0287     booker.setCurrentFolder(fmt::format("{}/VsLumi", dqmDir_));
0288     h_layer_vsLumi.push_back(EffME1(booker.book1D(Form("layertotal_vsLumi_layer_%i", layer), lyrName, 100, 0, 25000),
0289                                     booker.book1D(Form("layerfound_vsLumi_layer_%i", layer), lyrName, 100, 0, 25000)));
0290 
0291     // book plots vs Lumi
0292     booker.setCurrentFolder(fmt::format("{}/VsPu", dqmDir_));
0293     h_layer_vsPU.push_back(EffME1(booker.book1D(Form("layertotal_vsPU_layer_%i", layer), lyrName, 45, 0, 90),
0294                                   booker.book1D(Form("layerfound_vsPU_layer_%i", layer), lyrName, 45, 0, 90)));
0295     if (addCommonMode_) {
0296       // book plots for common mode
0297       booker.setCurrentFolder(fmt::format("{}/CommonMode", dqmDir_));
0298       h_layer_vsCM.push_back(EffME1(booker.book1D(Form("layertotal_vsCM_layer_%i", layer), lyrName, 20, 0, 400),
0299                                     booker.book1D(Form("layerfound_vsCM_layer_%i", layer), lyrName, 20, 0, 400)));
0300     }
0301 
0302     // book plots vs Lumi
0303     booker.setCurrentFolder(fmt::format("{}/VsBx", dqmDir_));
0304     h_layer_vsBx.push_back(EffME1(
0305         booker.book1D(Form("totalVsBx_layer%i", layer), Form("layer %i (%s)", layer, lyrName.c_str()), 3565, 0, 3565),
0306         booker.book1D(Form("foundVsBx_layer%i", layer), Form("layer %i (%s)", layer, lyrName.c_str()), 3565, 0, 3565)));
0307 
0308     // book hot and cold
0309     booker.setCurrentFolder(fmt::format("{}/MissingHits", dqmDir_));
0310     if (layer <= bounds::k_LayersAtTOBEnd) {
0311       const bool isTIB = layer <= bounds::k_LayersAtTIBEnd;
0312       const auto partition = (isTIB ? "TIB" : "TOB");
0313       const auto yMax = (isTIB ? 100 : 120);
0314 
0315       const auto tit = Form("%s%i: Map of missing hits", partition, (isTIB ? layer : layer - bounds::k_LayersAtTIBEnd));
0316 
0317       auto ihhotcold = booker.book2D(tit, tit, 100, -1, 361, 100, -yMax, yMax);
0318       ihhotcold->setAxisTitle("#phi [deg]", 1);
0319       ihhotcold->setBinLabel(1, "360", 1);
0320       ihhotcold->setBinLabel(50, "180", 1);
0321       ihhotcold->setBinLabel(100, "0", 1);
0322       ihhotcold->setAxisTitle("Global Z [cm]", 2);
0323       ihhotcold->setOption("colz");
0324       h_hotcold.push_back(ihhotcold);
0325     } else {
0326       const bool isTID = layer <= bounds::k_LayersAtTIDEnd;
0327       const auto partitions =
0328           (isTID ? std::vector<std::string>{"TID-", "TID+"} : std::vector<std::string>{"TEC-", "TEC+"});
0329       const auto axMax = (isTID ? 100 : 120);
0330       for (const auto& part : partitions) {
0331         const auto tit = Form("%s%i: Map of missing hits",
0332                               part.c_str(),
0333                               (isTID ? layer - bounds::k_LayersAtTOBEnd : layer - bounds::k_LayersAtTIDEnd));
0334 
0335         auto ihhotcold = booker.book2D(tit, tit, 100, -axMax, axMax, 100, -axMax, axMax);
0336         ihhotcold->setAxisTitle("Global Y", 1);
0337         ihhotcold->setBinLabel(1, "+Y", 1);
0338         ihhotcold->setBinLabel(50, "0", 1);
0339         ihhotcold->setBinLabel(100, "-Y", 1);
0340         ihhotcold->setAxisTitle("Global X", 2);
0341         ihhotcold->setBinLabel(1, "-X", 2);
0342         ihhotcold->setBinLabel(50, "0", 2);
0343         ihhotcold->setBinLabel(100, "+X", 2);
0344         ihhotcold->setOption("colz");
0345         h_hotcold.push_back(ihhotcold);
0346       }
0347     }
0348   }
0349 
0350   // come back to the main folder
0351   booker.setCurrentFolder(dqmDir_);
0352   const auto tkDetMapFolder = fmt::format("{}/TkDetMaps", dqmDir_);
0353 
0354   const TkDetMap* tkDetMap = &setup.getData(tkDetMapToken_);
0355   h_module =
0356       EffTkMap(std::make_unique<TkHistoMap>(tkDetMap, booker, tkDetMapFolder, "perModule_total", 0, false, true),
0357                std::make_unique<TkHistoMap>(tkDetMap, booker, tkDetMapFolder, "perModule_found", 0, false, true));
0358 }
0359 
0360 void SiStripHitEfficiencyWorker::analyze(const edm::Event& e, const edm::EventSetup& es) {
0361   const auto tTopo = &es.getData(tTopoToken_);
0362 
0363   //  bool DEBUG_ = false;
0364 
0365   LogDebug("SiStripHitEfficiencyWorker") << "beginning analyze from HitEff";
0366 
0367   // Step A: Get Inputs
0368 
0369   // Luminosity informations
0370   edm::Handle<LumiScalersCollection> lumiScalers = e.getHandle(scalerToken_);
0371   edm::Handle<OnlineLuminosityRecord> metaData = e.getHandle(metaDataToken_);
0372 
0373   float instLumi = 0;
0374   float PU = 0;
0375   if (addLumi_) {
0376     if (lumiScalers.isValid() && !lumiScalers->empty()) {
0377       if (lumiScalers->begin() != lumiScalers->end()) {
0378         instLumi = lumiScalers->begin()->instantLumi();
0379         PU = lumiScalers->begin()->pileup();
0380       }
0381     } else if (metaData.isValid()) {
0382       instLumi = metaData->instLumi();
0383       PU = metaData->avgPileUp();
0384     } else {
0385       edm::LogWarning("SiStripHitEfficiencyWorker") << "could not find a source for the Luminosity and PU";
0386     }
0387   }
0388 
0389   h_bx->Fill(e.bunchCrossing());
0390   h_instLumi->Fill(instLumi);
0391   h_PU->Fill(PU);
0392 
0393   edm::Handle<edm::DetSetVector<SiStripRawDigi>> commonModeDigis;
0394   if (addCommonMode_)
0395     e.getByToken(commonModeToken_, commonModeDigis);
0396 
0397   edm::Handle<reco::TrackCollection> tracksCKF;
0398   e.getByToken(combinatorialTracks_token_, tracksCKF);
0399 
0400   edm::Handle<std::vector<Trajectory>> TrajectoryCollectionCKF;
0401   e.getByToken(trajectories_token_, TrajectoryCollectionCKF);
0402 
0403   edm::Handle<TrajTrackAssociationCollection> trajTrackAssociationHandle;
0404   e.getByToken(trajTrackAsso_token_, trajTrackAssociationHandle);
0405 
0406   edm::Handle<edmNew::DetSetVector<SiStripCluster>> theClusters;
0407   e.getByToken(clusters_token_, theClusters);
0408 
0409   edm::Handle<DetIdCollection> fedErrorIds;
0410   e.getByToken(digis_token_, fedErrorIds);
0411 
0412   edm::Handle<MeasurementTrackerEvent> measurementTrackerEvent;
0413   e.getByToken(trackerEvent_token_, measurementTrackerEvent);
0414 
0415   const auto tkgeom = &es.getData(tkGeomToken_);
0416   const auto& stripcpe = es.getData(stripCPEToken_);
0417   const auto& stripQuality = es.getData(stripQualityToken_);
0418   const auto& magField = es.getData(magFieldToken_);
0419   const auto& measTracker = es.getData(measTrackerToken_);
0420   const auto& chi2Estimator = es.getData(chi2EstimatorToken_);
0421   const auto& prop = es.getData(propagatorToken_);
0422 
0423   ++events;
0424 
0425   // Tracking
0426   LogDebug("SiStripHitEfficiencyWorker") << "number ckf tracks found = " << tracksCKF->size();
0427 
0428   h_nTracks->Fill(tracksCKF->size());
0429   h_nTracksVsPU->Fill(PU, tracksCKF->size());
0430 
0431   if (!tracksCKF->empty()) {
0432     if (cutOnTracks_ && (tracksCKF->size() >= trackMultiplicityCut_))
0433       return;
0434     if (cutOnTracks_)
0435       LogDebug("SiStripHitEfficiencyWorker")
0436           << "starting checking good event with < " << trackMultiplicityCut_ << " tracks";
0437 
0438     ++EventTrackCKF;
0439 
0440     // actually should do a loop over all the tracks in the event here
0441 
0442     // Looping over traj-track associations to be able to get traj & track informations
0443     for (const auto& trajTrack : *trajTrackAssociationHandle) {
0444       // for each track, fill some variables such as number of hits and momentum
0445 
0446       const bool highPurity = trajTrack.val->quality(reco::TrackBase::TrackQuality::highPurity);
0447       auto TMeas = trajTrack.key->measurements();
0448 
0449       const bool hasMissingHits = std::any_of(std::begin(TMeas), std::end(TMeas), [](const auto& tm) {
0450         return tm.recHit()->getType() == TrackingRecHit::Type::missing;
0451       });
0452 
0453       // Loop on each measurement and take it into consideration
0454       //--------------------------------------------------------
0455       for (auto itm = TMeas.cbegin(); itm != TMeas.cend(); ++itm) {
0456         const auto theInHit = (*itm).recHit();
0457 
0458         LogDebug("SiStripHitEfficiencyWorker") << "theInHit is valid = " << theInHit->isValid();
0459 
0460         unsigned int iidd = theInHit->geographicalId().rawId();
0461 
0462         unsigned int TKlayers = ::checkLayer(iidd, tTopo);
0463 
0464         // do not bother with pixel hits
0465         if (DetId(iidd).subdetId() < SiStripSubdetector::TIB)
0466           continue;
0467 
0468         LogDebug("SiStripHitEfficiencyWorker") << "TKlayer from trajectory: " << TKlayers << "  from module = " << iidd
0469                                                << "   matched/stereo/rphi = " << ((iidd & 0x3) == 0) << "/"
0470                                                << ((iidd & 0x3) == 1) << "/" << ((iidd & 0x3) == 2);
0471 
0472         // Test first and last points of the trajectory
0473         // the list of measurements starts from outer layers  !!! This could change -> should add a check
0474         if ((!useFirstMeas_ && (itm == (TMeas.end() - 1))) || (!useLastMeas_ && (itm == (TMeas.begin()))) ||
0475             // In case of missing hit in the track, check whether to use the other hits or not.
0476             (!useAllHitsFromTracksWithMissingHits_ && hasMissingHits &&
0477              theInHit->getType() != TrackingRecHit::Type::missing))
0478           continue;
0479         // If Trajectory measurement from TOB 6 or TEC 9, skip it because it's always valid they are filled later
0480         if (TKlayers == bounds::k_LayersAtTOBEnd || TKlayers == bounds::k_LayersAtTECEnd) {
0481           LogDebug("SiStripHitEfficiencyWorker") << "skipping original TM for TOB 6 or TEC 9";
0482           continue;
0483         }
0484 
0485         std::vector<TrajectoryAtInvalidHit> TMs;
0486 
0487         // Make AnalyticalPropagat // TODO where to save these?or to use in TAVH constructor
0488         AnalyticalPropagator propagator(&magField, anyDirection);
0489 
0490         // for double sided layers check both sensors--if no hit was found on either sensor surface,
0491         // the trajectory measurements only have one invalid hit entry on the matched surface
0492         // so get the TrajectoryAtInvalidHit for both surfaces and include them in the study
0493         if (::isDoubleSided(iidd, tTopo) && ((iidd & 0x3) == 0)) {
0494           // do hit eff check twice--once for each sensor
0495           //add a TM for each surface
0496           TMs.emplace_back(*itm, tTopo, tkgeom, propagator, 1);
0497           TMs.emplace_back(*itm, tTopo, tkgeom, propagator, 2);
0498         } else if (::isDoubleSided(iidd, tTopo) && (!::check2DPartner(iidd, TMeas))) {
0499           // if only one hit was found the trajectory measurement is on that sensor surface, and the other surface from
0500           // the matched layer should be added to the study as well
0501           TMs.emplace_back(*itm, tTopo, tkgeom, propagator, 1);
0502           TMs.emplace_back(*itm, tTopo, tkgeom, propagator, 2);
0503           LogDebug("SiStripHitEfficiencyWorker") << " found a hit with a missing partner";
0504         } else {
0505           //only add one TM for the single surface and the other will be added in the next iteration
0506           TMs.emplace_back(*itm, tTopo, tkgeom, propagator);
0507         }
0508 
0509         //////////////////////////////////////////////
0510         //Now check for tracks at TOB6 and TEC9
0511 
0512         // to make sure we only propagate on the last TOB5 hit check the next entry isn't also in TOB5
0513         // to avoid bias, make sure the TOB5 hit is valid (an invalid hit on TOB5 could only exist with a valid hit on TOB6)
0514         const auto nextId = (itm + 1 != TMeas.end()) ? (itm + 1)->recHit()->geographicalId() : DetId{};  // null if last
0515 
0516         if (TKlayers == 9 && theInHit->isValid() && !((!nextId.null()) && (::checkLayer(nextId.rawId(), tTopo) == 9))) {
0517           //      if ( TKlayers==9 && itm==TMeas.rbegin()) {
0518           //      if ( TKlayers==9 && (itm==TMeas.back()) ) {     // to check for only the last entry in the trajectory for propagation
0519           const DetLayer* tob6 = measTracker.geometricSearchTracker()->tobLayers().back();
0520           const LayerMeasurements theLayerMeasurements{measTracker, *measurementTrackerEvent};
0521           const TrajectoryStateOnSurface tsosTOB5 = itm->updatedState();
0522           const auto tmp = theLayerMeasurements.measurements(*tob6, tsosTOB5, prop, chi2Estimator);
0523 
0524           if (!tmp.empty()) {
0525             LogDebug("SiStripHitEfficiencyWorker") << "size of TM from propagation = " << tmp.size();
0526 
0527             // take the last of the TMs, which is always an invalid hit
0528             // if no detId is available, ie detId==0, then no compatible layer was crossed
0529             // otherwise, use that TM for the efficiency measurement
0530             const auto& tob6TM = tmp.back();
0531             const auto& tob6Hit = tob6TM.recHit();
0532             if (tob6Hit->geographicalId().rawId() != 0) {
0533               LogDebug("SiStripHitEfficiencyWorker") << "tob6 hit actually being added to TM vector";
0534               TMs.emplace_back(tob6TM, tTopo, tkgeom, propagator);
0535             }
0536           }
0537         }
0538 
0539         // same for TEC8
0540         if (TKlayers == 21 && theInHit->isValid() &&
0541             !((!nextId.null()) && (::checkLayer(nextId.rawId(), tTopo) == 21))) {
0542           const DetLayer* tec9pos = measTracker.geometricSearchTracker()->posTecLayers().back();
0543           const DetLayer* tec9neg = measTracker.geometricSearchTracker()->negTecLayers().back();
0544 
0545           const LayerMeasurements theLayerMeasurements{measTracker, *measurementTrackerEvent};
0546           const TrajectoryStateOnSurface tsosTEC9 = itm->updatedState();
0547 
0548           // check if track on positive or negative z
0549           if (!(iidd == SiStripSubdetector::TEC))
0550             LogDebug("SiStripHitEfficiencyWorker") << "there is a problem with TEC 9 extrapolation";
0551 
0552           //LogDebug("SiStripHitEfficiencyWorker") << " tec9 id = " << iidd << " and side = " << tTopo->tecSide(iidd) ;
0553           std::vector<TrajectoryMeasurement> tmp;
0554           if (tTopo->tecSide(iidd) == 1) {
0555             tmp = theLayerMeasurements.measurements(*tec9neg, tsosTEC9, prop, chi2Estimator);
0556             //LogDebug("SiStripHitEfficiencyWorker") << "on negative side" ;
0557           }
0558           if (tTopo->tecSide(iidd) == 2) {
0559             tmp = theLayerMeasurements.measurements(*tec9pos, tsosTEC9, prop, chi2Estimator);
0560             //LogDebug("SiStripHitEfficiencyWorker") << "on positive side" ;
0561           }
0562 
0563           if (!tmp.empty()) {
0564             // take the last of the TMs, which is always an invalid hit
0565             // if no detId is available, ie detId==0, then no compatible layer was crossed
0566             // otherwise, use that TM for the efficiency measurement
0567             const auto& tec9TM = tmp.back();
0568             const auto& tec9Hit = tec9TM.recHit();
0569 
0570             const unsigned int tec9id = tec9Hit->geographicalId().rawId();
0571             LogDebug("SiStripHitEfficiencyWorker")
0572                 << "tec9id = " << tec9id << " is Double sided = " << ::isDoubleSided(tec9id, tTopo)
0573                 << "  and 0x3 = " << (tec9id & 0x3);
0574 
0575             if (tec9Hit->geographicalId().rawId() != 0) {
0576               LogDebug("SiStripHitEfficiencyWorker") << "tec9 hit actually being added to TM vector";
0577               // in tec the hit can be single or doubled sided. whenever the invalid hit at the end of vector of TMs is
0578               // double sided it is always on the matched surface, so we need to split it into the true sensor surfaces
0579               if (::isDoubleSided(tec9id, tTopo)) {
0580                 TMs.emplace_back(tec9TM, tTopo, tkgeom, propagator, 1);
0581                 TMs.emplace_back(tec9TM, tTopo, tkgeom, propagator, 2);
0582               } else
0583                 TMs.emplace_back(tec9TM, tTopo, tkgeom, propagator);
0584             }
0585           }  //else LogDebug("SiStripHitEfficiencyWorker") << "tec9 tmp empty" ;
0586         }
0587 
0588         for (const auto& tm : TMs) {
0589           fillForTraj(tm,
0590                       tTopo,
0591                       tkgeom,
0592                       stripcpe,
0593                       stripQuality,
0594                       *fedErrorIds,
0595                       commonModeDigis,
0596                       *theClusters,
0597                       e.bunchCrossing(),
0598                       instLumi,
0599                       PU,
0600                       highPurity);
0601         }
0602         LogDebug("SiStripHitEfficiencyWorker") << "After looping over TrajAtValidHit list";
0603       }
0604       LogDebug("SiStripHitEfficiencyWorker") << "end TMeasurement loop";
0605     }
0606     LogDebug("SiStripHitEfficiencyWorker") << "end of trajectories loop";
0607   }
0608 }
0609 
0610 void SiStripHitEfficiencyWorker::fillForTraj(const TrajectoryAtInvalidHit& tm,
0611                                              const TrackerTopology* tTopo,
0612                                              const TrackerGeometry* tkgeom,
0613                                              const StripClusterParameterEstimator& stripCPE,
0614                                              const SiStripQuality& stripQuality,
0615                                              const DetIdCollection& fedErrorIds,
0616                                              const edm::Handle<edm::DetSetVector<SiStripRawDigi>>& commonModeDigis,
0617                                              const edmNew::DetSetVector<SiStripCluster>& theClusters,
0618                                              int bunchCrossing,
0619                                              float instLumi,
0620                                              float PU,
0621                                              bool highPurity) {
0622   // --> Get trajectory from combinatedStat& e
0623   const auto iidd = tm.monodet_id();
0624   LogDebug("SiStripHitEfficiencyWorker") << "setting iidd = " << iidd << " before checking efficiency and ";
0625 
0626   const auto xloc = tm.localX();
0627   const auto yloc = tm.localY();
0628 
0629   const auto xErr = tm.localErrorX();
0630   const auto yErr = tm.localErrorY();
0631 
0632   int TrajStrip = -1;
0633 
0634   // reget layer from iidd here, to account for TOB 6 and TEC 9 TKlayers being off
0635   const auto TKlayers = ::checkLayer(iidd, tTopo);
0636 
0637   const bool withinAcceptance =
0638       tm.withinAcceptance() && (!::isInBondingExclusionZone(iidd, TKlayers, yloc, yErr, tTopo));
0639 
0640   if (  // (TKlayers > 0) && // FIXME confirm this
0641       ((layers_ == TKlayers) ||
0642        (layers_ == bounds::k_LayersStart))) {  // Look at the layer not used to reconstruct the track
0643     LogDebug("SiStripHitEfficiencyWorker") << "Looking at layer under study";
0644     unsigned int ModIsBad = 2;
0645     unsigned int SiStripQualBad = 0;
0646     float commonMode = -100;
0647 
0648     // RPhi RecHit Efficiency
0649 
0650     if (!theClusters.empty()) {
0651       LogDebug("SiStripHitEfficiencyWorker") << "Checking clusters with size = " << theClusters.size();
0652       std::vector<::ClusterInfo> VCluster_info;  //fill with X residual, X residual pull, local X
0653       const auto idsv = theClusters.find(iidd);
0654       if (idsv != theClusters.end()) {
0655         //if (DEBUG_)      LogDebug("SiStripHitEfficiencyWorker") << "the ID from the dsv = " << dsv.id() ;
0656         LogDebug("SiStripHitEfficiencyWorker")
0657             << "found  (ClusterId == iidd) with ClusterId = " << idsv->id() << " and iidd = " << iidd;
0658         const auto stripdet = dynamic_cast<const StripGeomDetUnit*>(tkgeom->idToDetUnit(DetId(iidd)));
0659         const StripTopology& Topo = stripdet->specificTopology();
0660 
0661         float hbedge = 0.0;
0662         float htedge = 0.0;
0663         float hapoth = 0.0;
0664         float uylfac = 0.0;
0665         float uxlden = 0.0;
0666         if (TKlayers > bounds::k_LayersAtTOBEnd) {
0667           const BoundPlane& plane = stripdet->surface();
0668           const TrapezoidalPlaneBounds* trapezoidalBounds(
0669               dynamic_cast<const TrapezoidalPlaneBounds*>(&(plane.bounds())));
0670           std::array<const float, 4> const& parameterTrap = (*trapezoidalBounds).parameters();  // el bueno aqui
0671           hbedge = parameterTrap[0];
0672           htedge = parameterTrap[1];
0673           hapoth = parameterTrap[3];
0674           uylfac = (htedge - hbedge) / (htedge + hbedge) / hapoth;
0675           uxlden = 1 + yloc * uylfac;
0676         }
0677 
0678         // Need to know position of trajectory in strip number for selecting the right APV later
0679         if (TrajStrip == -1) {
0680           int nstrips = Topo.nstrips();
0681           float pitch = stripdet->surface().bounds().width() / nstrips;
0682           TrajStrip = xloc / pitch + nstrips / 2.0;
0683           // Need additionnal corrections for endcap
0684           if (TKlayers > bounds::k_LayersAtTOBEnd) {
0685             const float TrajLocXMid = xloc / (1 + (htedge - hbedge) * yloc / (htedge + hbedge) /
0686                                                       hapoth);  // radialy extrapolated x loc position at middle
0687             TrajStrip = TrajLocXMid / pitch + nstrips / 2.0;
0688           }
0689           //LogDebug("SiStripHitEfficiency")<<" Layer "<<TKlayers<<" TrajStrip: "<<nstrips<<" "<<pitch<<" "<<TrajStrip;;
0690         }
0691 
0692         for (const auto& clus : *idsv) {
0693           StripClusterParameterEstimator::LocalValues parameters = stripCPE.localParameters(clus, *stripdet);
0694           float res = (parameters.first.x() - xloc);
0695           float sigma = ::checkConsistency(parameters, xloc, xErr);
0696           // The consistency is probably more accurately measured with the Chi2MeasurementEstimator. To use it
0697           // you need a TransientTrackingRecHit instead of the cluster
0698           //theEstimator=       new Chi2MeasurementEstimator(30);
0699           //const Chi2MeasurementEstimator *theEstimator(100);
0700           //theEstimator->estimate(tm.tsos(), TransientTrackingRecHit);
0701 
0702           if (TKlayers > bounds::k_LayersAtTOBEnd) {
0703             res = parameters.first.x() - xloc / uxlden;  // radialy extrapolated x loc position at middle
0704             sigma = abs(res) / sqrt(parameters.second.xx() + xErr * xErr / uxlden / uxlden +
0705                                     yErr * yErr * xloc * xloc * uylfac * uylfac / uxlden / uxlden / uxlden / uxlden);
0706           }
0707 
0708           VCluster_info.emplace_back(res, sigma, parameters.first.x());
0709 
0710           LogDebug("SiStripHitEfficiencyWorker") << "Have ID match. residual = " << res << "  res sigma = " << sigma;
0711           //LogDebug("SiStripHitEfficiencyWorker")
0712           //    << "trajectory measurement compatability estimate = " << (*itm).estimate() ;
0713           LogDebug("SiStripHitEfficiencyWorker")
0714               << "hit position = " << parameters.first.x() << "  hit error = " << sqrt(parameters.second.xx())
0715               << "  trajectory position = " << xloc << "  traj error = " << xErr;
0716         }
0717       }
0718       ::ClusterInfo finalCluster{1000.0, 1000.0, 0.0};
0719       if (!VCluster_info.empty()) {
0720         LogDebug("SiStripHitEfficiencyWorker") << "found clusters > 0";
0721         if (VCluster_info.size() > 1) {
0722           //get the smallest one
0723           for (const auto& res : VCluster_info) {
0724             if (std::abs(res.xResidualPull) < std::abs(finalCluster.xResidualPull)) {
0725               finalCluster = res;
0726             }
0727             LogDebug("SiStripHitEfficiencyWorker")
0728                 << "iresidual = " << res.xResidual << "  isigma = " << res.xResidualPull
0729                 << "  and FinalRes = " << finalCluster.xResidual;
0730           }
0731         } else {
0732           finalCluster = VCluster_info[0];
0733         }
0734         VCluster_info.clear();
0735       }
0736 
0737       LogDebug("SiStripHitEfficiencyWorker") << "Final residual in X = " << finalCluster.xResidual << "+-"
0738                                              << (finalCluster.xResidual / finalCluster.xResidualPull);
0739       LogDebug("SiStripHitEfficiencyWorker")
0740           << "Checking location of trajectory: abs(yloc) = " << abs(yloc) << "  abs(xloc) = " << abs(xloc);
0741 
0742       //
0743       // fill ntuple varibles
0744 
0745       //if ( stripQuality->IsModuleBad(iidd) )
0746       if (stripQuality.getBadApvs(iidd) != 0) {
0747         SiStripQualBad = 1;
0748         LogDebug("SiStripHitEfficiencyWorker") << "strip is bad from SiStripQuality";
0749       } else {
0750         SiStripQualBad = 0;
0751         LogDebug("SiStripHitEfficiencyWorker") << "strip is good from SiStripQuality";
0752       }
0753 
0754       //check for FED-detected errors and include those in SiStripQualBad
0755       for (unsigned int ii = 0; ii < fedErrorIds.size(); ii++) {
0756         if (iidd == fedErrorIds[ii].rawId())
0757           SiStripQualBad = 1;
0758       }
0759 
0760       // CM of APV crossed by traj
0761       if (addCommonMode_)
0762         if (commonModeDigis.isValid() && TrajStrip >= 0 && TrajStrip <= 768) {
0763           const auto digiframe = commonModeDigis->find(iidd);
0764           if (digiframe != commonModeDigis->end())
0765             if ((unsigned)TrajStrip / sistrip::STRIPS_PER_APV < digiframe->data.size())
0766               commonMode = digiframe->data.at(TrajStrip / sistrip::STRIPS_PER_APV).adc();
0767         }
0768 
0769       LogDebug("SiStripHitEfficiencyWorker") << "before check good";
0770 
0771       if (finalCluster.xResidualPull < 999.0) {  //could make requirement on track/hit consistency, but for
0772         //now take anything with a hit on the module
0773         LogDebug("SiStripHitEfficiencyWorker")
0774             << "hit being counted as good " << finalCluster.xResidual << " FinalRecHit " << iidd << "   TKlayers  "
0775             << TKlayers << " xloc " << xloc << " yloc  " << yloc << " module " << iidd
0776             << "   matched/stereo/rphi = " << ((iidd & 0x3) == 0) << "/" << ((iidd & 0x3) == 1) << "/"
0777             << ((iidd & 0x3) == 2);
0778         ModIsBad = 0;
0779       } else {
0780         LogDebug("SiStripHitEfficiencyWorker")
0781             << "hit being counted as bad   ######### Invalid RPhi FinalResX " << finalCluster.xResidual
0782             << " FinalRecHit " << iidd << "   TKlayers  " << TKlayers << " xloc " << xloc << " yloc  " << yloc
0783             << " module " << iidd << "   matched/stereo/rphi = " << ((iidd & 0x3) == 0) << "/" << ((iidd & 0x3) == 1)
0784             << "/" << ((iidd & 0x3) == 2);
0785         ModIsBad = 1;
0786         LogDebug("SiStripHitEfficiencyWorker")
0787             << " RPhi Error " << sqrt(xErr * xErr + yErr * yErr) << " ErrorX " << xErr << " yErr " << yErr;
0788       }
0789 
0790       LogDebug("SiStripHitEfficiencyWorker")
0791           << "To avoid them staying unused: ModIsBad=" << ModIsBad << ", SiStripQualBad=" << SiStripQualBad
0792           << ", commonMode=" << commonMode << ", highPurity=" << highPurity
0793           << ", withinAcceptance=" << withinAcceptance;
0794 
0795       unsigned int layer = TKlayers;
0796       if (showRings_ && layer > bounds::k_LayersAtTOBEnd) {  // use rings instead of wheels
0797         if (layer <= bounds::k_LayersAtTIDEnd) {             // TID
0798           layer = bounds::k_LayersAtTOBEnd +
0799                   tTopo->tidRing(iidd);  // ((iidd >> 9) & 0x3);  // 3 disks and also 3 rings -> use the same container
0800         } else {                         // TEC
0801           layer = bounds::k_LayersAtTIDEnd + tTopo->tecRing(iidd);  // ((iidd >> 5) & 0x7);
0802         }
0803       }
0804       unsigned int layerWithSide = layer;
0805       if (layer > bounds::k_LayersAtTOBEnd && layer <= bounds::k_LayersAtTIDEnd) {
0806         const auto side = tTopo->tidSide(iidd);  //(iidd >> 13) & 0x3;  // TID
0807         if (side == 2)
0808           layerWithSide = layer + 3;
0809       } else if (layer > bounds::k_LayersAtTIDEnd) {
0810         const auto side = tTopo->tecSide(iidd);  // (iidd >> 18) & 0x3;  // TEC
0811         if (side == 1) {
0812           layerWithSide = layer + 3;
0813         } else if (side == 2) {
0814           layerWithSide = layer + 3 + (showRings_ ? 7 : 9);
0815         }
0816       }
0817 
0818       if ((bunchX_ > 0 && bunchX_ != bunchCrossing) || (!withinAcceptance) ||
0819           (useOnlyHighPurityTracks_ && !highPurity) ||
0820           (!showTOB6TEC9_ && (TKlayers == bounds::k_LayersAtTOBEnd || TKlayers == bounds::k_LayersAtTECEnd)) ||
0821           (badModules_.end() != badModules_.find(iidd)))
0822         return;
0823 
0824       const bool badquality = (SiStripQualBad == 1);
0825 
0826       //Now that we have a good event, we need to look at if we expected it or not, and the location
0827       //if we didn't
0828       //Fill the missing hit information first
0829       bool badflag = false;  // true for hits that are expected but not found
0830       if (resXSig_ < 0) {
0831         if (ModIsBad == 1)
0832           badflag = true;  // isBad set to false in the tree when resxsig<999.0
0833       } else {
0834         if (ModIsBad == 1 || finalCluster.xResidualPull > resXSig_)
0835           badflag = true;
0836       }
0837 
0838       // Conversion of positions in strip unit
0839       int nstrips = -9;
0840       float Pitch = -9.0;
0841       const StripGeomDetUnit* stripdet = nullptr;
0842       if (finalCluster.xResidualPull ==
0843           1000.0) {      // special treatment, no GeomDetUnit associated in some cases when no cluster found
0844         Pitch = 0.0205;  // maximum
0845         nstrips = 768;   // maximum
0846       } else {
0847         stripdet = dynamic_cast<const StripGeomDetUnit*>(tkgeom->idToDetUnit(iidd));
0848         const StripTopology& Topo = stripdet->specificTopology();
0849         nstrips = Topo.nstrips();
0850         Pitch = stripdet->surface().bounds().width() / Topo.nstrips();
0851       }
0852       double stripTrajMid = xloc / Pitch + nstrips / 2.0;
0853       double stripCluster = finalCluster.xLocal / Pitch + nstrips / 2.0;
0854       // For trapezoidal modules: extrapolation of x trajectory position to the y middle of the module
0855       //  for correct comparison with cluster position
0856       if (stripdet && layer > bounds::k_LayersAtTOBEnd) {
0857         const auto& trapezoidalBounds = dynamic_cast<const TrapezoidalPlaneBounds&>(stripdet->surface().bounds());
0858         std::array<const float, 4> const& parameters = trapezoidalBounds.parameters();
0859         const float hbedge = parameters[0];
0860         const float htedge = parameters[1];
0861         const float hapoth = parameters[3];
0862         const float TrajLocXMid = xloc / (1 + (htedge - hbedge) * yloc / (htedge + hbedge) /
0863                                                   hapoth);  // radialy extrapolated x loc position at middle
0864         stripTrajMid = TrajLocXMid / Pitch + nstrips / 2.0;
0865       }
0866 
0867       if ((!badquality) && (layer < h_resolution.size())) {
0868         LogDebug("SiStripHitEfficiencyWorker")
0869             << "layer " << layer << " vector index " << layer - 1 << " before filling h_resolution" << std::endl;
0870         h_resolution[layer - 1]->Fill(finalCluster.xResidualPull != 1000.0 ? stripTrajMid - stripCluster : 1000);
0871       }
0872 
0873       // New matching methods
0874       if (clusterMatchingMethod_ >= 1) {
0875         badflag = false;
0876         if (finalCluster.xResidualPull == 1000.0) {
0877           LogDebug("SiStripHitEfficiencyWorker") << "Marking bad for resxsig=1000";
0878           badflag = true;
0879         } else {
0880           if (clusterMatchingMethod_ == 2 || clusterMatchingMethod_ == 4) {
0881             // check the distance between cluster and trajectory position
0882             if (std::abs(stripCluster - stripTrajMid) > clusterTracjDist_) {
0883               LogDebug("SiStripHitEfficiencyWorker") << "Marking bad for cluster-to-traj distance";
0884               badflag = true;
0885             }
0886           }
0887           if (clusterMatchingMethod_ == 3 || clusterMatchingMethod_ == 4) {
0888             // cluster and traj have to be in the same APV (don't take edges into accounts)
0889             const int tapv = (int)stripTrajMid / sistrip::STRIPS_PER_APV;
0890             const int capv = (int)stripCluster / sistrip::STRIPS_PER_APV;
0891             float stripInAPV = stripTrajMid - tapv * sistrip::STRIPS_PER_APV;
0892             if (stripInAPV < stripsApvEdge_ || stripInAPV > sistrip::STRIPS_PER_APV - stripsApvEdge_) {
0893               LogDebug("SiStripHitEfficiencyWorker") << "Too close to the edge: " << stripInAPV;
0894               return;
0895             }
0896             if (tapv != capv) {
0897               LogDebug("SiStripHitEfficiencyWorker") << "Marking bad for tapv!=capv";
0898               badflag = true;
0899             }
0900           }
0901         }
0902       }
0903       if (!badquality) {
0904         LogDebug("SiStripHitEfficiencyWorker")
0905             << "Filling measurement for " << iidd << " in layer " << layer << " histograms with bx=" << bunchCrossing
0906             << ", lumi=" << instLumi << ", PU=" << PU << "; bad flag=" << badflag;
0907 
0908         // hot/cold maps of hits that are expected but not found
0909         if (badflag) {
0910           if (layer > bounds::k_LayersStart && layer <= bounds::k_LayersAtTIBEnd) {
0911             //We are in the TIB
0912             float phi = ::calcPhi(tm.globalX(), tm.globalY());
0913             h_hotcold[layer - 1]->Fill(360. - phi, tm.globalZ(), 1.);
0914           } else if (layer > bounds::k_LayersAtTIBEnd && layer <= bounds::k_LayersAtTOBEnd) {
0915             //We are in the TOB
0916             float phi = ::calcPhi(tm.globalX(), tm.globalY());
0917             h_hotcold[layer - 1]->Fill(360. - phi, tm.globalZ(), 1.);
0918           } else if (layer > bounds::k_LayersAtTOBEnd && layer <= bounds::k_LayersAtTIDEnd) {
0919             //We are in the TID
0920             //There are 2 different maps here
0921             int side = tTopo->tidSide(iidd);
0922             if (side == 1)
0923               h_hotcold[(layer - 1) + (layer - 11)]->Fill(-tm.globalY(), tm.globalX(), 1.);
0924             else if (side == 2)
0925               h_hotcold[(layer - 1) + (layer - 10)]->Fill(-tm.globalY(), tm.globalX(), 1.);
0926           } else if (layer > bounds::k_LayersAtTIDEnd) {
0927             //We are in the TEC
0928             //There are 2 different maps here
0929             int side = tTopo->tecSide(iidd);
0930             if (side == 1)
0931               h_hotcold[(layer + 2) + (layer - 14)]->Fill(-tm.globalY(), tm.globalX(), 1.);
0932             else if (side == 2)
0933               h_hotcold[(layer + 2) + (layer - 13)]->Fill(-tm.globalY(), tm.globalX(), 1.);
0934           }
0935         }
0936 
0937         LogDebug("SiStripHitEfficiencyWorker")
0938             << "layer " << layer << " vector index " << layer - 1 << " before filling h_layer_vsSmthg" << std::endl;
0939         h_layer_vsBx[layer - 1].fill(bunchCrossing, !badflag);
0940         if (addLumi_) {
0941           h_layer_vsLumi[layer - 1].fill(instLumi, !badflag);
0942           h_layer_vsPU[layer - 1].fill(PU, !badflag);
0943         }
0944         if (addCommonMode_) {
0945           h_layer_vsCM[layer - 1].fill(commonMode, !badflag);
0946         }
0947         h_goodLayer.fill(layerWithSide, !badflag);
0948 
0949         // efficiency with bad modules excluded
0950         if (TKlayers) {
0951           h_module.fill(iidd, !badflag);
0952         }
0953       }
0954       // efficiency without bad modules excluded
0955       h_allLayer.fill(layerWithSide, !badflag);
0956 
0957       /* Used in SiStripHitEffFromCalibTree:
0958        * run              -> "run"              -> run              // e.id().run()
0959        * event            -> "event"            -> evt              // e.id().event()
0960        * ModIsBad         -> "ModIsBad"         -> isBad
0961        * SiStripQualBad   -> "SiStripQualBad""  -> quality
0962        * Id               -> "Id"               -> id               // iidd
0963        * withinAcceptance -> "withinAcceptance" -> accept
0964        * whatlayer        -> "layer"            -> layer_wheel      // Tklayers
0965        * highPurity       -> "highPurity"       -> highPurity
0966        * TrajGlbX         -> "TrajGlbX"         -> x                // tm.globalX()
0967        * TrajGlbY         -> "TrajGlbY"         -> y                // tm.globalY()
0968        * TrajGlbZ         -> "TrajGlbZ"         -> z                // tm.globalZ()
0969        * ResXSig          -> "ResXSig"          -> resxsig          // finalCluster.xResidualPull;
0970        * TrajLocX         -> "TrajLocX"         -> TrajLocX         // xloc
0971        * TrajLocY         -> "TrajLocY"         -> TrajLocY         // yloc
0972        * ClusterLocX      -> "ClusterLocX"      -> ClusterLocX      // finalCluster.xLocal
0973        * bunchx           -> "bunchx"           -> bx               // e.bunchCrossing()
0974        * instLumi         -> "instLumi"         -> instLumi         ## if addLumi_
0975        * PU               -> "PU"               -> PU               ## if addLumi_
0976        * commonMode       -> "commonMode"       -> CM               ## if addCommonMode_ / _useCM
0977        */
0978       LogDebug("SiStripHitEfficiencyWorker") << "after good location check";
0979     }
0980     LogDebug("SiStripHitEfficiencyWorker") << "after list of clusters";
0981   }
0982   LogDebug("SiStripHitEfficiencyWorker") << "After layers=TKLayers if with TKlayers=" << TKlayers
0983                                          << ", layers=" << layers_;
0984 }
0985 
0986 void SiStripHitEfficiencyWorker::endJob() {
0987   LogDebug("SiStripHitEfficiencyWorker") << " Events Analysed             " << events;
0988   LogDebug("SiStripHitEfficiencyWorker") << " Number Of Tracked events    " << EventTrackCKF;
0989 }
0990 
0991 void SiStripHitEfficiencyWorker::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0992   edm::ParameterSetDescription desc;
0993   desc.add<std::string>("dqmDir", "AlCaReco/SiStripHitEfficiency");
0994   desc.add<bool>("UseOnlyHighPurityTracks", true);
0995   desc.add<bool>("cutOnTracks", false);
0996   desc.add<bool>("useAllHitsFromTracksWithMissingHits", false);
0997   desc.add<bool>("useFirstMeas", false);
0998   desc.add<bool>("useLastMeas", false);
0999   desc.add<double>("ClusterTrajDist", 64.0);
1000   desc.add<double>("ResXSig", -1);
1001   desc.add<double>("StripsApvEdge", 10.0);
1002   desc.add<edm::InputTag>("combinatorialTracks", edm::InputTag{"generalTracks"});
1003   desc.add<edm::InputTag>("commonMode", edm::InputTag{"siStripDigis", "CommonMode"});
1004   desc.add<edm::InputTag>("lumiScalers", edm::InputTag{"scalersRawToDigi"});
1005   desc.add<edm::InputTag>("metadata", edm::InputTag{"onlineMetaDataDigis"});
1006   desc.add<edm::InputTag>("siStripClusters", edm::InputTag{"siStripClusters"});
1007   desc.add<edm::InputTag>("siStripDigis", edm::InputTag{"siStripDigis"});
1008   desc.add<edm::InputTag>("trackerEvent", edm::InputTag{"MeasurementTrackerEvent"});
1009   desc.add<edm::InputTag>("trajectories", edm::InputTag{"generalTracks"});
1010   desc.add<int>("ClusterMatchingMethod", 0);
1011   desc.add<int>("Layer", 0);
1012   desc.add<unsigned int>("trackMultiplicity", 100);
1013   desc.addUntracked<bool>("Debug", false);
1014   desc.addUntracked<bool>("ShowRings", false);
1015   desc.addUntracked<bool>("ShowTOB6TEC9", false);
1016   desc.addUntracked<bool>("addCommonMode", false);
1017   desc.addUntracked<bool>("addLumi", true);
1018   desc.addUntracked<int>("BunchCrossing", 0);
1019   desc.addUntracked<std::string>("BadModulesFile", "");
1020   descriptions.addWithDefaultLabel(desc);
1021 }
1022 
1023 #include "FWCore/Framework/interface/MakerMacros.h"
1024 DEFINE_FWK_MODULE(SiStripHitEfficiencyWorker);