Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-09-10 01:02:40 

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