OfflineValidation/plugins/DMRChecker.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    Alignment/OfflineValidation
0004 // Class:      DMRChecker
0005 //
0006 /**\class DMRChecker DMRChecker.cc Alignment/OfflineValidation/plugins/DMRChecker.cc
0007
0008 */
0009 //
0010 // Original Author:  Marco Musich
0011 //         Created:  Mon, 10 Aug 2020 15:45:00 GMT
0012 //
0013 //
0014
0015 // ROOT includes
0016
0017 #include "RtypesCore.h"
0018 #include "TAxis.h"
0019 #include "TCanvas.h"
0020 #include "TColor.h"
0021 #include "TH1.h"
0022 #include "TH2.h"
0023 #include "TLatex.h"
0024 #include "TMath.h"
0025 #include "TProfile.h"
0026 #include "TString.h"
0027 #include "TStyle.h"
0028 #include "TVirtualPad.h"
0029
0030 // STL includes
0031
0032 #include <algorithm>
0033 #include <array>
0034 #include <cmath>
0035 #include <cstdint>
0036 #include <cstdlib>
0037 #include <ctime>
0038 #include <ext/alloc_traits.h>
0039 #include <fmt/printf.h>
0040 #include <iomanip>
0041 #include <iostream>
0042 #include <map>
0043 #include <memory>
0044 #include <string>
0045 #include <type_traits>
0046 #include <utility>
0047 #include <vector>
0048 #include "boost/range/adaptor/indexed.hpp"
0049
0050 // user system includes
0051
0052 #include "CalibTracker/StandaloneTrackerTopology/interface/StandaloneTrackerTopology.h"
0053 #include "CommonTools/TrackerMap/interface/TrackerMap.h"
0054 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0055 #include "CommonTools/Utils/interface/TFileDirectory.h"
0056 #include "DQM/TrackerRemapper/interface/Phase1PixelMaps.h"
0057 #include "DQM/TrackerRemapper/interface/Phase1PixelSummaryMap.h"
0058 #include "CondCore/SiPixelPlugins/interface/PixelRegionContainers.h"
0059 #include "CondCore/SiPixelPlugins/interface/SiPixelPayloadInspectorHelper.h"
0060 #include "CondFormats/DataRecord/interface/RunSummaryRcd.h"
0061 #include "CondFormats/DataRecord/interface/SiStripCondDataRecords.h"
0062 #include "CondFormats/RunInfo/interface/RunInfo.h"
0063 #include "CondFormats/SiStripObjects/interface/SiStripLatency.h"
0064 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0065 #include "DataFormats/Common/interface/Handle.h"
0066 #include "DataFormats/Common/interface/RefTraits.h"
0067 #include "DataFormats/Common/interface/TriggerResults.h"
0068 #include "DataFormats/DetId/interface/DetId.h"
0069 #include "DataFormats/GeometrySurface/interface/Plane.h"
0070 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0071 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0072 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0073 #include "DataFormats/GeometryVector/interface/Phi.h"
0074 #include "DataFormats/GeometryVector/interface/Theta.h"
0075 #include "DataFormats/Math/interface/Point3D.h"
0076 #include "DataFormats/Math/interface/deltaPhi.h"
0077 #include "DataFormats/Provenance/interface/RunLumiEventNumber.h"
0078 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0079 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
0080 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0081 #include "DataFormats/TrackReco/interface/HitPattern.h"
0082 #include "DataFormats/TrackReco/interface/Track.h"
0083 #include "DataFormats/TrackReco/interface/TrackBase.h"
0084 #include "DataFormats/TrackReco/interface/TrackExtraFwd.h"
0085 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0086 #include "DataFormats/TrackReco/interface/TrackResiduals.h"
0087 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0088 #include "DataFormats/TrackerRecHit2D/interface/ProjectedSiStripRecHit2D.h"
0089 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
0090 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
0091 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit1D.h"
0092 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
0093 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0094 #include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h"
0095 #include "DataFormats/VertexReco/interface/Vertex.h"
0096 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0097 #include "FWCore/Common/interface/TriggerNames.h"
0098 #include "FWCore/Framework/interface/Event.h"
0099 #include "FWCore/Framework/interface/EventSetup.h"
0100 #include "FWCore/Framework/interface/MakerMacros.h"
0101 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0102 #include "FWCore/MessageLogger/interface/ErrorObj.h"
0103 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0104 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0105 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0106 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0107 #include "FWCore/ParameterSet/interface/ParameterSetDescriptionFiller.h"
0108 #include "FWCore/ServiceRegistry/interface/Service.h"
0109 #include "FWCore/Utilities/interface/EDGetToken.h"
0110 #include "FWCore/Utilities/interface/ESGetToken.h"
0111 #include "FWCore/Utilities/interface/InputTag.h"
0112 #include "Geometry/CommonTopologies/interface/GeomDet.h"
0113 #include "Geometry/CommonTopologies/interface/TrackerGeomDet.h"
0114 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0115 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0116 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0117 #include "MagneticField/Engine/interface/MagneticField.h"
0118 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0119
0120 #define DEBUG 0
0121
0122 using namespace std;
0123 using namespace edm;
0124
0125 const int kBPIX = PixelSubdetector::PixelBarrel;
0126 const int kFPIX = PixelSubdetector::PixelEndcap;
0127 constexpr float cmToUm = 10000.;
0128
0129 /**
0130  * Auxilliary POD to store the data for
0131  * the running mean algorithm.
0132  */
0133
0134 namespace running {
0135   struct Estimators {
0136     int rDirection;
0137     int zDirection;
0138     int rOrZDirection;
0139     int hitCount;
0140     float runningMeanOfRes_;
0141     float runningVarOfRes_;
0142     float runningNormMeanOfRes_;
0143     float runningNormVarOfRes_;
0144   };
0145
0146   using estimatorMap = std::map<uint32_t, running::Estimators>;
0147
0148 }  // namespace running
0149
0150 class DMRChecker : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0151 public:
0152   DMRChecker(const edm::ParameterSet &pset)
0153       : geomToken_(esConsumes<edm::Transition::BeginRun>()),
0154         runInfoToken_(esConsumes<edm::Transition::BeginRun>()),
0155         magFieldToken_(esConsumes<edm::Transition::BeginRun>()),
0156         topoToken_(esConsumes<edm::Transition::BeginRun>()),
0157         isCosmics_(pset.getParameter<bool>("isCosmics")),
0158         doLatencyAnalysis_(pset.getParameter<bool>("doLatencyAnalysis")) {
0159     usesResource(TFileService::kSharedResource);
0160
0161     if (doLatencyAnalysis_) {
0162       latencyToken_ = esConsumes<edm::Transition::BeginRun>();
0163     } else {
0164       latencyToken_ = edm::ESGetToken<SiStripLatency, SiStripLatencyRcd>();
0165     }
0166
0167     TkTag_ = pset.getParameter<edm::InputTag>("TkTag");
0168     theTrackCollectionToken_ = consumes<reco::TrackCollection>(TkTag_);
0169
0170     TriggerResultsTag_ = pset.getParameter<edm::InputTag>("TriggerResultsTag");
0171     hltresultsToken_ = consumes<edm::TriggerResults>(TriggerResultsTag_);
0172
0173     BeamSpotTag_ = pset.getParameter<edm::InputTag>("BeamSpotTag");
0174     beamspotToken_ = consumes<reco::BeamSpot>(BeamSpotTag_);
0175
0176     VerticesTag_ = pset.getParameter<edm::InputTag>("VerticesTag");
0177     vertexToken_ = consumes<reco::VertexCollection>(VerticesTag_);
0178
0179     // initialize conventional Tracker maps
0180
0181     pmap = std::make_unique<TrackerMap>("Pixel");
0182     pmap->onlyPixel(true);
0183     pmap->setTitle("Pixel Hit entries");
0184     pmap->setPalette(1);
0185
0186     tmap = std::make_unique<TrackerMap>("Strip");
0187     tmap->setTitle("Strip Hit entries");
0188     tmap->setPalette(1);
0189
0190     // initialize Phase1 Pixel Maps
0191
0192     pixelmap = std::make_unique<Phase1PixelMaps>("COLZ0 L");
0193     pixelmap->bookBarrelHistograms("DMRsX", "Median Residuals x-direction", "Median Residuals");
0194     pixelmap->bookForwardHistograms("DMRsX", "Median Residuals x-direction", "Median Residuals");
0195
0196     pixelmap->bookBarrelHistograms("DMRsY", "Median Residuals y-direction", "Median Residuals");
0197     pixelmap->bookForwardHistograms("DMRsY", "Median Residuals y-direction", "Median Residuals");
0198
0199     // set no rescale
0200     pixelmap->setNoRescale();
0201
0202     // initialize Full Pixel Map
0203     fullPixelmapXDMR = std::make_unique<Phase1PixelSummaryMap>("", "DMR-x", "median of residuals [#mum]");
0204     fullPixelmapYDMR = std::make_unique<Phase1PixelSummaryMap>("", "DMR-y", "median of residuals [#mum]");
0205   }
0206
0207   static void fillDescriptions(edm::ConfigurationDescriptions &);
0208
0209   ~DMRChecker() override = default;
0210
0211   /*_______________________________________________________
0212   //
0213   // auxilliary method to retrieve certain histogram types
0214   //_______________________________________________________
0215   */
0216   template <class OBJECT_TYPE>
0217   int index(const std::vector<OBJECT_TYPE *> &vec, const std::string &name) {
0218     for (const auto &iter : vec | boost::adaptors::indexed(0)) {
0219       if (iter.value() && iter.value()->GetName() == name) {
0220         return iter.index();
0221       }
0222     }
0223     edm::LogError("Alignment") << "@SUB=DMRChecker::index"
0224                                << " could not find " << name;
0225     return -1;
0226   }
0227
0228   template <typename T, typename... Args>
0229   T *book(const Args &...args) const {
0230     T *t = fs->make<T>(args...);
0231     return t;
0232   }
0233
0234 private:
0235   // tokens for the event setup
0236   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geomToken_;
0237   const edm::ESGetToken<RunInfo, RunInfoRcd> runInfoToken_;
0238   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magFieldToken_;
0239   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> topoToken_;
0240   // not const
0241   edm::ESGetToken<SiStripLatency, SiStripLatencyRcd> latencyToken_;
0242
0243   const MagneticField *magneticField_;
0244   const TrackerGeometry *trackerGeometry_;
0245   const TrackerTopology *trackerTopology_;
0246
0247   edm::Service<TFileService> fs;
0248
0249   std::unique_ptr<Phase1PixelMaps> pixelmap;
0250   std::unique_ptr<Phase1PixelSummaryMap> fullPixelmapXDMR;
0251   std::unique_ptr<Phase1PixelSummaryMap> fullPixelmapYDMR;
0252
0253   std::unique_ptr<PixelRegions::PixelRegionContainers> PixelDMRS_x_ByLayer;
0254   std::unique_ptr<PixelRegions::PixelRegionContainers> PixelDMRS_y_ByLayer;
0255
0256   std::unique_ptr<TrackerMap> tmap;
0257   std::unique_ptr<TrackerMap> pmap;
0258
0259   TH1D *hchi2ndof;
0260   TH1D *hNtrk;
0261   TH1D *hNtrkZoom;
0262   TH1I *htrkQuality;
0263   TH1I *htrkAlgo;
0264   TH1D *hNhighPurity;
0265   TH1D *hP;
0266   TH1D *hPPlus;
0267   TH1D *hPMinus;
0268   TH1D *hPt;
0269   TH1D *hMinPt;
0270   TH1D *hPtPlus;
0271   TH1D *hPtMinus;
0272   TH1D *hHit;
0273   TH1D *hHit2D;
0274
0275   TH1D *hBPixResXPrime;
0276   TH1D *hFPixResXPrime;
0277   TH1D *hFPixZPlusResXPrime;
0278   TH1D *hFPixZMinusResXPrime;
0279
0280   TH1D *hBPixResYPrime;
0281   TH1D *hFPixResYPrime;
0282   TH1D *hFPixZPlusResYPrime;
0283   TH1D *hFPixZMinusResYPrime;
0284
0285   TH1D *hBPixResXPull;
0286   TH1D *hFPixResXPull;
0287   TH1D *hFPixZPlusResXPull;
0288   TH1D *hFPixZMinusResXPull;
0289
0290   TH1D *hBPixResYPull;
0291   TH1D *hFPixResYPull;
0292   TH1D *hFPixZPlusResYPull;
0293   TH1D *hFPixZMinusResYPull;
0294
0295   TH1D *hTIBResXPrime;
0296   TH1D *hTOBResXPrime;
0297   TH1D *hTIDResXPrime;
0298   TH1D *hTECResXPrime;
0299
0300   TH1D *hTIBResXPull;
0301   TH1D *hTOBResXPull;
0302   TH1D *hTIDResXPull;
0303   TH1D *hTECResXPull;
0304
0305   TH1D *hHitCountVsXBPix;
0306   TH1D *hHitCountVsXFPix;
0307   TH1D *hHitCountVsYBPix;
0308   TH1D *hHitCountVsYFPix;
0309   TH1D *hHitCountVsZBPix;
0310   TH1D *hHitCountVsZFPix;
0311
0312   TH1D *hHitCountVsThetaBPix;
0313   TH1D *hHitCountVsPhiBPix;
0314
0315   TH1D *hHitCountVsThetaFPix;
0316   TH1D *hHitCountVsPhiFPix;
0317
0318   TH1D *hHitCountVsXFPixPlus;
0319   TH1D *hHitCountVsXFPixMinus;
0320   TH1D *hHitCountVsYFPixPlus;
0321   TH1D *hHitCountVsYFPixMinus;
0322   TH1D *hHitCountVsZFPixPlus;
0323   TH1D *hHitCountVsZFPixMinus;
0324
0325   TH1D *hHitCountVsThetaFPixPlus;
0326   TH1D *hHitCountVsPhiFPixPlus;
0327
0328   TH1D *hHitCountVsThetaFPixMinus;
0329   TH1D *hHitCountVsPhiFPixMinus;
0330
0331   TH1D *hHitPlus;
0332   TH1D *hHitMinus;
0333
0334   TH1D *hPhp;
0335   TH1D *hPthp;
0336   TH1D *hHithp;
0337   TH1D *hEtahp;
0338   TH1D *hPhihp;
0339   TH1D *hchi2ndofhp;
0340   TH1D *hchi2Probhp;
0341
0342   TH1D *hCharge;
0343   TH1D *hQoverP;
0344   TH1D *hQoverPZoom;
0345   TH1D *hEta;
0346   TH1D *hEtaPlus;
0347   TH1D *hEtaMinus;
0348   TH1D *hPhi;
0349   TH1D *hPhiBarrel;
0350   TH1D *hPhiOverlapPlus;
0351   TH1D *hPhiOverlapMinus;
0352   TH1D *hPhiEndcapPlus;
0353   TH1D *hPhiEndcapMinus;
0354   TH1D *hPhiPlus;
0355   TH1D *hPhiMinus;
0356
0357   TH1D *hDeltaPhi;
0358   TH1D *hDeltaEta;
0359   TH1D *hDeltaR;
0360
0361   TH1D *hvx;
0362   TH1D *hvy;
0363   TH1D *hvz;
0364   TH1D *hd0;
0365   TH1D *hdz;
0366   TH1D *hdxy;
0367
0368   TH2D *hd0PVvsphi;
0369   TH2D *hd0PVvseta;
0370   TH2D *hd0PVvspt;
0371
0372   TH2D *hd0vsphi;
0373   TH2D *hd0vseta;
0374   TH2D *hd0vspt;
0375
0376   TH1D *hnhpxb;
0377   TH1D *hnhpxe;
0378   TH1D *hnhTIB;
0379   TH1D *hnhTID;
0380   TH1D *hnhTOB;
0381   TH1D *hnhTEC;
0382
0383   TH1D *hHitComposition;
0384
0385   TProfile *pNBpixHitsVsVx;
0386   TProfile *pNBpixHitsVsVy;
0387   TProfile *pNBpixHitsVsVz;
0388
0389   TH1D *hMultCand;
0390
0391   TH1D *hdxyBS;
0392   TH1D *hd0BS;
0393   TH1D *hdzBS;
0394   TH1D *hdxyPV;
0395   TH1D *hd0PV;
0396   TH1D *hdzPV;
0397   TH1D *hrun;
0398   TH1D *hlumi;
0399
0400   std::vector<TH1 *> vTrackHistos_;
0401   std::vector<TH1 *> vTrackProfiles_;
0402   std::vector<TH1 *> vTrack2DHistos_;
0403
0404   TH1D *tksByTrigger_;
0405   TH1D *evtsByTrigger_;
0406
0407   TH1D *modeByRun_;
0408   TH1D *fieldByRun_;
0409
0410   TH1D *tracksByRun_;
0411   TH1D *hitsByRun_;
0412
0413   TH1D *trackRatesByRun_;
0414   TH1D *eventRatesByRun_;
0415
0416   TH1D *hitsinBPixByRun_;
0417   TH1D *hitsinFPixByRun_;
0418
0419   // Pixel
0420
0421   TH1D *DMRBPixX_;
0422   TH1D *DMRBPixY_;
0423
0424   TH1D *DMRFPixX_;
0425   TH1D *DMRFPixY_;
0426
0427   TH1D *DRnRBPixX_;
0428   TH1D *DRnRBPixY_;
0429
0430   TH1D *DRnRFPixX_;
0431   TH1D *DRnRFPixY_;
0432
0433   // Strips
0434
0435   TH1D *DMRTIB_;
0436   TH1D *DMRTOB_;
0437
0438   TH1D *DMRTID_;
0439   TH1D *DMRTEC_;
0440
0441   TH1D *DRnRTIB_;
0442   TH1D *DRnRTOB_;
0443
0444   TH1D *DRnRTID_;
0445   TH1D *DRnRTEC_;
0446
0447   // Split DMRs
0448
0449   std::array<TH1D *, 2> DMRBPixXSplit_;
0450   std::array<TH1D *, 2> DMRBPixYSplit_;
0451
0452   std::array<TH1D *, 2> DMRFPixXSplit_;
0453   std::array<TH1D *, 2> DMRFPixYSplit_;
0454
0455   std::array<TH1D *, 2> DMRTIBSplit_;
0456   std::array<TH1D *, 2> DMRTOBSplit_;
0457
0458   // residuals
0459
0460   std::map<unsigned int, TH1D *> barrelLayersResidualsX;
0461   std::map<unsigned int, TH1D *> barrelLayersPullsX;
0462   std::map<unsigned int, TH1D *> barrelLayersResidualsY;
0463   std::map<unsigned int, TH1D *> barrelLayersPullsY;
0464
0465   std::map<unsigned int, TH1D *> endcapDisksResidualsX;
0466   std::map<unsigned int, TH1D *> endcapDisksPullsX;
0467   std::map<unsigned int, TH1D *> endcapDisksResidualsY;
0468   std::map<unsigned int, TH1D *> endcapDisksPullsY;
0469
0470   int ievt;
0471   int itrks;
0472   int mode;
0473
0474   SiPixelPI::phase phase_;
0475   float etaMax_;
0476
0477   const bool isCosmics_;
0478   const bool doLatencyAnalysis_;
0479
0480   edm::InputTag TkTag_;
0481   edm::InputTag TriggerResultsTag_;
0482   edm::InputTag BeamSpotTag_;
0483   edm::InputTag VerticesTag_;
0484
0485   edm::EDGetTokenT<reco::TrackCollection> theTrackCollectionToken_;
0486   edm::EDGetTokenT<edm::TriggerResults> hltresultsToken_;
0487   edm::EDGetTokenT<reco::BeamSpot> beamspotToken_;
0488   edm::EDGetTokenT<reco::VertexCollection> vertexToken_;
0489
0490   std::map<std::string, std::pair<int, int> > triggerMap_;
0491   std::map<int, std::pair<int, float> > conditionsMap_;
0492   std::map<int, std::pair<int, int> > runInfoMap_;
0493   std::map<int, std::array<int, 6> > runHitsMap_;
0494
0495   std::map<int, float> timeMap_;
0496
0497   // Pixel
0498
0499   running::estimatorMap resDetailsBPixX_;
0500   running::estimatorMap resDetailsBPixY_;
0501   running::estimatorMap resDetailsFPixX_;
0502   running::estimatorMap resDetailsFPixY_;
0503
0504   // Strips
0505
0506   running::estimatorMap resDetailsTIB_;
0507   running::estimatorMap resDetailsTOB_;
0508   running::estimatorMap resDetailsTID_;
0509   running::estimatorMap resDetailsTEC_;
0510
0511   void analyze(const edm::Event &event, const edm::EventSetup &setup) override {
0512     ievt++;
0513
0514     edm::Handle<reco::TrackCollection> trackCollection = event.getHandle(theTrackCollectionToken_);
0515
0516     GlobalPoint zeroPoint(0, 0, 0);
0517     if (DEBUG)
0518       edm::LogVerbatim("DMRChecker") << "event #" << ievt << " Event ID = " << event.id()
0519                                      << " magnetic field: " << magneticField_->inTesla(zeroPoint) << std::endl;
0520
0521     const reco::TrackCollection tC = *(trackCollection.product());
0522     itrks += tC.size();
0523
0524     runInfoMap_[event.run()].first += 1;
0525     runInfoMap_[event.run()].second += tC.size();
0526
0527     if (DEBUG)
0528       edm::LogVerbatim("DMRChecker") << "Reconstructed " << tC.size() << " tracks" << std::endl;
0529
0530     edm::Handle<edm::TriggerResults> hltresults = event.getHandle(hltresultsToken_);
0531     if (hltresults.isValid()) {
0532       const edm::TriggerNames &triggerNames_ = event.triggerNames(*hltresults);
0533       int ntrigs = hltresults->size();
0534
0535       for (int itrig = 0; itrig != ntrigs; ++itrig) {
0536         const string &trigName = triggerNames_.triggerName(itrig);
0537         bool accept = hltresults->accept(itrig);
0538         if (accept == 1) {
0539           if (DEBUG)
0540             edm::LogVerbatim("DMRChecker") << trigName << " " << accept << " ,track size: " << tC.size() << endl;
0541           triggerMap_[trigName].first += 1;
0542           triggerMap_[trigName].second += tC.size();
0543         }
0544       }
0545     }
0546
0547     hrun->Fill(event.run());
0548     hlumi->Fill(event.luminosityBlock());
0549
0550     int nHighPurityTracks = 0;
0551
0552     for (const auto &track : tC) {
0553       auto const &residuals = track.extra()->residuals();
0554
0555       unsigned int nHit2D = 0;
0556       int h_index = 0;
0557       for (trackingRecHit_iterator iHit = track.recHitsBegin(); iHit != track.recHitsEnd(); ++iHit, ++h_index) {
0558         if (this->isHit2D(**iHit))
0559           ++nHit2D;
0560
0561         double resX = residuals.residualX(h_index);
0562         double resY = residuals.residualY(h_index);
0563         double pullX = residuals.pullX(h_index);
0564         double pullY = residuals.pullY(h_index);
0565
0566         const DetId &detId = (*iHit)->geographicalId();
0567
0568         unsigned int subid = detId.subdetId();
0569         uint32_t detid_db = detId.rawId();
0570
0571         const GeomDet *geomDet(trackerGeometry_->idToDet(detId));
0572
0573         float uOrientation(-999.F), vOrientation(-999.F);
0574         LocalPoint lPModule(0., 0., 0.), lUDirection(1., 0., 0.), lVDirection(0., 1., 0.), lWDirection(0., 0., 1.);
0575
0576         // do all the transformations here
0577         GlobalPoint gUDirection = geomDet->surface().toGlobal(lUDirection);
0578         GlobalPoint gVDirection = geomDet->surface().toGlobal(lVDirection);
0579         GlobalPoint gWDirection = geomDet->surface().toGlobal(lWDirection);
0580         GlobalPoint gPModule = geomDet->surface().toGlobal(lPModule);
0581
0582         if (!(*iHit)->detUnit())
0583           continue;  // is it a single physical module?
0584
0585         if ((*iHit)->isValid() && (subid > PixelSubdetector::PixelEndcap)) {
0586           if (phase_ != SiPixelPI::phase::two)
0587             tmap->fill(detid_db, 1);
0588
0589           //LocalPoint lp = (*iHit)->localPosition();
0590           //LocalError le = (*iHit)->localPositionError();
0591
0592           // fill DMRs and DrNRs
0593           if (subid == StripSubdetector::TIB) {
0594             uOrientation = deltaPhi(gUDirection.barePhi(), gPModule.barePhi()) >= 0. ? +1.F : -1.F;
0595             //vOrientation = gVDirection.z() - gPModule.z() >= 0 ? +1.F : -1.F; // not used for Strips
0596
0597             // if the detid has never occcurred yet, set the local orientations
0598             if (resDetailsTIB_.find(detid_db) == resDetailsTIB_.end()) {
0599               resDetailsTIB_[detid_db].rDirection = gWDirection.perp() - gPModule.perp() >= 0 ? +1 : -1;
0600               resDetailsTIB_[detid_db].zDirection = gVDirection.z() - gPModule.z() >= 0 ? +1 : -1;
0601               resDetailsTIB_[detid_db].rOrZDirection = resDetailsTIB_[detid_db].rDirection;  // barrel (split in r)
0602             }
0603
0604             hTIBResXPrime->Fill(uOrientation * resX * cmToUm);
0605             hTIBResXPull->Fill(pullX);
0606
0607             // update residuals
0608             this->updateOnlineMomenta(resDetailsTIB_, detid_db, uOrientation * resX * cmToUm, pullX);
0609
0610           } else if (subid == StripSubdetector::TOB) {
0611             uOrientation = deltaPhi(gUDirection.barePhi(), gPModule.barePhi()) >= 0. ? +1.F : -1.F;
0612             //vOrientation = gVDirection.z() - gPModule.z() >= 0 ? +1.F : -1.F; // not used for Strips
0613
0614             hTOBResXPrime->Fill(uOrientation * resX * cmToUm);
0615             hTOBResXPull->Fill(pullX);
0616
0617             // if the detid has never occcurred yet, set the local orientations
0618             if (resDetailsTOB_.find(detid_db) == resDetailsTOB_.end()) {
0619               resDetailsTOB_[detid_db].rDirection = gWDirection.perp() - gPModule.perp() >= 0 ? +1 : -1;
0620               resDetailsTOB_[detid_db].zDirection = gVDirection.z() - gPModule.z() >= 0 ? +1 : -1;
0621               resDetailsTOB_[detid_db].rOrZDirection = resDetailsTOB_[detid_db].rDirection;  // barrel (split in r)
0622             }
0623
0624             // update residuals
0625             this->updateOnlineMomenta(resDetailsTOB_, detid_db, uOrientation * resX * cmToUm, pullX);
0626
0627           } else if (subid == StripSubdetector::TID) {
0628             uOrientation = deltaPhi(gUDirection.barePhi(), gPModule.barePhi()) >= 0. ? +1.F : -1.F;
0629             //vOrientation = gVDirection.perp() - gPModule.perp() >= 0. ? +1.F : -1.F; // not used for Strips
0630
0631             hTIDResXPrime->Fill(uOrientation * resX * cmToUm);
0632             hTIDResXPull->Fill(pullX);
0633
0634             // update residuals
0635             this->updateOnlineMomenta(resDetailsTID_, detid_db, uOrientation * resX * cmToUm, pullX);
0636
0637           } else if (subid == StripSubdetector::TEC) {
0638             uOrientation = deltaPhi(gUDirection.barePhi(), gPModule.barePhi()) >= 0. ? +1.F : -1.F;
0639             //vOrientation = gVDirection.perp() - gPModule.perp() >= 0. ? +1.F : -1.F; // not used for Strips
0640
0641             hTECResXPrime->Fill(uOrientation * resX * cmToUm);
0642             hTECResXPull->Fill(pullX);
0643
0644             // update residuals
0645             this->updateOnlineMomenta(resDetailsTEC_, detid_db, uOrientation * resX * cmToUm, pullX);
0646           }
0647         }
0648
0649         const SiPixelRecHit *pixhit = dynamic_cast<const SiPixelRecHit *>(*iHit);
0650
0651         if (pixhit) {
0652           if (pixhit->isValid()) {
0653             if (phase_ == SiPixelPI::phase::zero) {
0654               pmap->fill(detid_db, 1);
0655             }
0656
0657             LocalPoint lp = (*iHit)->localPosition();
0658             //LocalError le = (*iHit)->localPositionError();
0659             GlobalPoint GP = geomDet->surface().toGlobal(lp);
0660
0661             if ((subid == PixelSubdetector::PixelBarrel) || (subid == PixelSubdetector::PixelEndcap)) {
0662               // 1 = PXB, 2 = PXF
0663               if (subid == PixelSubdetector::PixelBarrel) {
0664                 int layer_num = trackerTopology_->pxbLayer(detid_db);
0665
0666                 uOrientation = deltaPhi(gUDirection.barePhi(), gPModule.barePhi()) >= 0. ? +1.F : -1.F;
0667                 vOrientation = gVDirection.z() - gPModule.z() >= 0 ? +1.F : -1.F;
0668
0669                 // if the detid has never occcurred yet, set the local orientations
0670                 if (resDetailsBPixX_.find(detid_db) == resDetailsBPixX_.end()) {
0671                   resDetailsBPixX_[detid_db].rDirection = gWDirection.perp() - gPModule.perp() >= 0 ? +1 : -1;
0672                   resDetailsBPixX_[detid_db].zDirection = gVDirection.z() - gPModule.z() >= 0 ? +1 : -1;
0673                   resDetailsBPixX_[detid_db].rOrZDirection =
0674                       resDetailsBPixX_[detid_db].rDirection;  // barrel (split in r)
0675                 }
0676
0677                 // if the detid has never occcurred yet, set the local orientations
0678                 if (resDetailsBPixY_.find(detid_db) == resDetailsBPixY_.end()) {
0679                   resDetailsBPixY_[detid_db].rDirection = gWDirection.perp() - gPModule.perp() >= 0 ? +1 : -1;
0680                   resDetailsBPixY_[detid_db].zDirection = gVDirection.z() - gPModule.z() >= 0 ? +1 : -1;
0681                   resDetailsBPixY_[detid_db].rOrZDirection =
0682                       resDetailsBPixY_[detid_db].rDirection;  // barrel (split in r)
0683                 }
0684
0685                 hHitCountVsThetaBPix->Fill(GP.theta());
0686                 hHitCountVsPhiBPix->Fill(GP.phi());
0687
0688                 hHitCountVsZBPix->Fill(GP.z());
0689                 hHitCountVsXBPix->Fill(GP.x());
0690                 hHitCountVsYBPix->Fill(GP.y());
0691
0692                 hBPixResXPrime->Fill(uOrientation * resX * cmToUm);
0693                 hBPixResYPrime->Fill(vOrientation * resY * cmToUm);
0694                 hBPixResXPull->Fill(pullX);
0695                 hBPixResYPull->Fill(pullY);
0696
0697                 if (DEBUG)
0698                   edm::LogVerbatim("DMRChecker") << "layer: " << layer_num << std::endl;
0699
0700                 // update residuals X
0701                 this->updateOnlineMomenta(resDetailsBPixX_, detid_db, uOrientation * resX * cmToUm, pullX);
0702
0703                 // update residuals Y
0704                 this->updateOnlineMomenta(resDetailsBPixY_, detid_db, vOrientation * resY * cmToUm, pullY);
0705
0706                 fillByIndex(barrelLayersResidualsX, layer_num, uOrientation * resX * cmToUm);
0707                 fillByIndex(barrelLayersPullsX, layer_num, pullX);
0708                 fillByIndex(barrelLayersResidualsY, layer_num, vOrientation * resY * cmToUm);
0709                 fillByIndex(barrelLayersPullsY, layer_num, pullY);
0710
0711               } else if (subid == PixelSubdetector::PixelEndcap) {
0712                 uOrientation = gUDirection.perp() - gPModule.perp() >= 0 ? +1.F : -1.F;
0713                 vOrientation = deltaPhi(gVDirection.barePhi(), gPModule.barePhi()) >= 0. ? +1.F : -1.F;
0714
0715                 int side_num = trackerTopology_->pxfSide(detid_db);
0716                 int disk_num = trackerTopology_->pxfDisk(detid_db);
0717
0718                 int packedTopo = disk_num + 3 * (side_num - 1);
0719
0720                 if (DEBUG)
0721                   edm::LogVerbatim("DMRChecker") << "side: " << side_num << " disk: " << disk_num
0722                                                  << " packedTopo: " << packedTopo << " GP.z(): " << GP.z() << std::endl;
0723
0724                 hHitCountVsThetaFPix->Fill(GP.theta());
0725                 hHitCountVsPhiFPix->Fill(GP.phi());
0726
0727                 hHitCountVsZFPix->Fill(GP.z());
0728                 hHitCountVsXFPix->Fill(GP.x());
0729                 hHitCountVsYFPix->Fill(GP.y());
0730
0731                 hFPixResXPrime->Fill(uOrientation * resX * cmToUm);
0732                 hFPixResYPrime->Fill(vOrientation * resY * cmToUm);
0733                 hFPixResXPull->Fill(pullX);
0734                 hFPixResYPull->Fill(pullY);
0735
0736                 fillByIndex(endcapDisksResidualsX, packedTopo, uOrientation * resX * cmToUm);
0737                 fillByIndex(endcapDisksPullsX, packedTopo, pullX);
0738                 fillByIndex(endcapDisksResidualsY, packedTopo, vOrientation * resY * cmToUm);
0739                 fillByIndex(endcapDisksPullsY, packedTopo, pullY);
0740
0741                 // if the detid has never occcurred yet, set the local orientations
0742                 if (resDetailsFPixX_.find(detid_db) == resDetailsFPixX_.end()) {
0743                   resDetailsFPixX_[detid_db].rDirection = gUDirection.perp() - gPModule.perp() >= 0 ? +1 : -1;
0744                   resDetailsFPixX_[detid_db].zDirection = gWDirection.z() - gPModule.z() >= 0 ? +1 : -1;
0745                   resDetailsFPixX_[detid_db].rOrZDirection =
0746                       resDetailsFPixX_[detid_db].zDirection;  // endcaps (split in z)
0747                 }
0748
0749                 // if the detid has never occcurred yet, set the local orientations
0750                 if (resDetailsFPixY_.find(detid_db) == resDetailsFPixY_.end()) {
0751                   resDetailsFPixY_[detid_db].rDirection = gUDirection.perp() - gPModule.perp() >= 0 ? +1 : -1;
0752                   resDetailsFPixY_[detid_db].zDirection = gWDirection.z() - gPModule.z() >= 0 ? +1 : -1;
0753                   resDetailsFPixY_[detid_db].rOrZDirection =
0754                       resDetailsFPixY_[detid_db].zDirection;  // endcaps (split in z)
0755                 }
0756
0757                 // update residuals X
0758                 this->updateOnlineMomenta(resDetailsFPixX_, detid_db, uOrientation * resX * cmToUm, pullX);
0759
0760                 // update residuals Y
0761                 this->updateOnlineMomenta(resDetailsFPixY_, detid_db, vOrientation * resY * cmToUm, pullY);
0762
0763                 if (side_num == 1) {
0764                   hHitCountVsXFPixMinus->Fill(GP.x());
0765                   hHitCountVsYFPixMinus->Fill(GP.y());
0766                   hHitCountVsZFPixMinus->Fill(GP.z());
0767                   hHitCountVsThetaFPixMinus->Fill(GP.theta());
0768                   hHitCountVsPhiFPixMinus->Fill(GP.phi());
0769
0770                   hFPixZMinusResXPrime->Fill(uOrientation * resX * cmToUm);
0771                   hFPixZMinusResYPrime->Fill(vOrientation * resY * cmToUm);
0772                   hFPixZMinusResXPull->Fill(pullX);
0773                   hFPixZMinusResYPull->Fill(pullY);
0774
0775                 } else {
0776                   hHitCountVsXFPixPlus->Fill(GP.x());
0777                   hHitCountVsYFPixPlus->Fill(GP.y());
0778                   hHitCountVsZFPixPlus->Fill(GP.z());
0779                   hHitCountVsThetaFPixPlus->Fill(GP.theta());
0780                   hHitCountVsPhiFPixPlus->Fill(GP.phi());
0781
0782                   hFPixZPlusResXPrime->Fill(uOrientation * resX * cmToUm);
0783                   hFPixZPlusResYPrime->Fill(vOrientation * resY * cmToUm);
0784                   hFPixZPlusResXPull->Fill(pullX);
0785                   hFPixZPlusResYPull->Fill(pullY);
0786                 }
0787               }
0788             }
0789           }
0790         }
0791       }
0792
0793       hHit2D->Fill(nHit2D);
0794       hHit->Fill(track.numberOfValidHits());
0795       hnhpxb->Fill(track.hitPattern().numberOfValidPixelBarrelHits());
0796       hnhpxe->Fill(track.hitPattern().numberOfValidPixelEndcapHits());
0797       hnhTIB->Fill(track.hitPattern().numberOfValidStripTIBHits());
0798       hnhTID->Fill(track.hitPattern().numberOfValidStripTIDHits());
0799       hnhTOB->Fill(track.hitPattern().numberOfValidStripTOBHits());
0800       hnhTEC->Fill(track.hitPattern().numberOfValidStripTECHits());
0801
0802       runHitsMap_[event.run()][0] += track.hitPattern().numberOfValidPixelBarrelHits();
0803       runHitsMap_[event.run()][1] += track.hitPattern().numberOfValidPixelEndcapHits();
0804       runHitsMap_[event.run()][2] += track.hitPattern().numberOfValidStripTIBHits();
0805       runHitsMap_[event.run()][3] += track.hitPattern().numberOfValidStripTIDHits();
0806       runHitsMap_[event.run()][4] += track.hitPattern().numberOfValidStripTOBHits();
0807       runHitsMap_[event.run()][5] += track.hitPattern().numberOfValidStripTECHits();
0808
0809       // fill hit composition histogram
0810       if (track.hitPattern().numberOfValidPixelBarrelHits() != 0) {
0811         hHitComposition->Fill(0., track.hitPattern().numberOfValidPixelBarrelHits());
0812
0813         pNBpixHitsVsVx->Fill(track.vx(), track.hitPattern().numberOfValidPixelBarrelHits());
0814         pNBpixHitsVsVy->Fill(track.vy(), track.hitPattern().numberOfValidPixelBarrelHits());
0815         pNBpixHitsVsVz->Fill(track.vz(), track.hitPattern().numberOfValidPixelBarrelHits());
0816       }
0817       if (track.hitPattern().numberOfValidPixelEndcapHits() != 0) {
0818         hHitComposition->Fill(1., track.hitPattern().numberOfValidPixelEndcapHits());
0819       }
0820       if (track.hitPattern().numberOfValidStripTIBHits() != 0) {
0821         hHitComposition->Fill(2., track.hitPattern().numberOfValidStripTIBHits());
0822       }
0823       if (track.hitPattern().numberOfValidStripTIDHits() != 0) {
0824         hHitComposition->Fill(3., track.hitPattern().numberOfValidStripTIDHits());
0825       }
0826       if (track.hitPattern().numberOfValidStripTOBHits() != 0) {
0827         hHitComposition->Fill(4., track.hitPattern().numberOfValidStripTOBHits());
0828       }
0829       if (track.hitPattern().numberOfValidStripTECHits() != 0) {
0830         hHitComposition->Fill(5., track.hitPattern().numberOfValidStripTECHits());
0831       }
0832
0833       hCharge->Fill(track.charge());
0834       hQoverP->Fill(track.qoverp());
0835       hQoverPZoom->Fill(track.qoverp());
0836       hPt->Fill(track.pt());
0837       hP->Fill(track.p());
0838       hchi2ndof->Fill(track.normalizedChi2());
0839       hEta->Fill(track.eta());
0840       hPhi->Fill(track.phi());
0841
0842       if (fabs(track.eta()) < 0.8)
0843         hPhiBarrel->Fill(track.phi());
0844       if (track.eta() > 0.8 && track.eta() < 1.4)
0845         hPhiOverlapPlus->Fill(track.phi());
0846       if (track.eta() < -0.8 && track.eta() > -1.4)
0847         hPhiOverlapMinus->Fill(track.phi());
0848       if (track.eta() > 1.4)
0849         hPhiEndcapPlus->Fill(track.phi());
0850       if (track.eta() < -1.4)
0851         hPhiEndcapMinus->Fill(track.phi());
0852
0853       hd0->Fill(track.d0());
0854       hdz->Fill(track.dz());
0855       hdxy->Fill(track.dxy());
0856       hvx->Fill(track.vx());
0857       hvy->Fill(track.vy());
0858       hvz->Fill(track.vz());
0859
0860       htrkAlgo->Fill(track.algo());
0861
0862       int myquality = -99;
0863       if (track.quality(reco::TrackBase::undefQuality)) {
0864         myquality = -1;
0865         htrkQuality->Fill(myquality);
0866       }
0867       if (track.quality(reco::TrackBase::loose)) {
0868         myquality = 0;
0869         htrkQuality->Fill(myquality);
0870       }
0871       if (track.quality(reco::TrackBase::tight)) {
0872         myquality = 1;
0873         htrkQuality->Fill(myquality);
0874       }
0875       if (track.quality(reco::TrackBase::highPurity) && (!isCosmics_)) {
0876         myquality = 2;
0877         htrkQuality->Fill(myquality);
0878         hPhp->Fill(track.p());
0879         hPthp->Fill(track.pt());
0880         hHithp->Fill(track.numberOfValidHits());
0881         hEtahp->Fill(track.eta());
0882         hPhihp->Fill(track.phi());
0883         hchi2ndofhp->Fill(track.normalizedChi2());
0884         hchi2Probhp->Fill(TMath::Prob(track.chi2(), track.ndof()));
0885         nHighPurityTracks++;
0886       }
0887       if (track.quality(reco::TrackBase::confirmed)) {
0888         myquality = 3;
0889         htrkQuality->Fill(myquality);
0890       }
0891       if (track.quality(reco::TrackBase::goodIterative)) {
0892         myquality = 4;
0893         htrkQuality->Fill(myquality);
0894       }
0895
0896       // Fill 1D track histos
0897       static const int etaindex = this->index(vTrackHistos_, "h_tracketa");
0898       vTrackHistos_[etaindex]->Fill(track.eta());
0899       static const int phiindex = this->index(vTrackHistos_, "h_trackphi");
0900       vTrackHistos_[phiindex]->Fill(track.phi());
0901       static const int numOfValidHitsindex = this->index(vTrackHistos_, "h_trackNumberOfValidHits");
0902       vTrackHistos_[numOfValidHitsindex]->Fill(track.numberOfValidHits());
0903       static const int numOfLostHitsindex = this->index(vTrackHistos_, "h_trackNumberOfLostHits");
0904       vTrackHistos_[numOfLostHitsindex]->Fill(track.numberOfLostHits());
0905
0906       GlobalPoint gPoint(track.vx(), track.vy(), track.vz());
0907       double theLocalMagFieldInInverseGeV = magneticField_->inInverseGeV(gPoint).z();
0908       double kappa = -track.charge() * theLocalMagFieldInInverseGeV / track.pt();
0909
0910       static const int kappaindex = this->index(vTrackHistos_, "h_curvature");
0911       vTrackHistos_[kappaindex]->Fill(kappa);
0912       static const int kappaposindex = this->index(vTrackHistos_, "h_curvature_pos");
0913       if (track.charge() > 0)
0914         vTrackHistos_[kappaposindex]->Fill(fabs(kappa));
0915       static const int kappanegindex = this->index(vTrackHistos_, "h_curvature_neg");
0916       if (track.charge() < 0)
0917         vTrackHistos_[kappanegindex]->Fill(fabs(kappa));
0918
0919       double chi2Prob = TMath::Prob(track.chi2(), track.ndof());
0920       double normchi2 = track.normalizedChi2();
0921
0922       static const int normchi2index = this->index(vTrackHistos_, "h_normchi2");
0923       vTrackHistos_[normchi2index]->Fill(normchi2);
0924       static const int chi2index = this->index(vTrackHistos_, "h_chi2");
0925       vTrackHistos_[chi2index]->Fill(track.chi2());
0926       static const int chi2Probindex = this->index(vTrackHistos_, "h_chi2Prob");
0927       vTrackHistos_[chi2Probindex]->Fill(chi2Prob);
0928       static const int ptindex = this->index(vTrackHistos_, "h_pt");
0929       static const int pt2index = this->index(vTrackHistos_, "h_ptrebin");
0930       vTrackHistos_[ptindex]->Fill(track.pt());
0931       vTrackHistos_[pt2index]->Fill(track.pt());
0932       if (track.ptError() != 0.) {
0933         static const int ptResolutionindex = this->index(vTrackHistos_, "h_ptResolution");
0934         vTrackHistos_[ptResolutionindex]->Fill(track.ptError() / track.pt());
0935       }
0936       // Fill track profiles
0937       static const int d0phiindex = this->index(vTrackProfiles_, "p_d0_vs_phi");
0938       vTrackProfiles_[d0phiindex]->Fill(track.phi(), track.d0());
0939       static const int dzphiindex = this->index(vTrackProfiles_, "p_dz_vs_phi");
0940       vTrackProfiles_[dzphiindex]->Fill(track.phi(), track.dz());
0941       static const int d0etaindex = this->index(vTrackProfiles_, "p_d0_vs_eta");
0942       vTrackProfiles_[d0etaindex]->Fill(track.eta(), track.d0());
0943       static const int dzetaindex = this->index(vTrackProfiles_, "p_dz_vs_eta");
0944       vTrackProfiles_[dzetaindex]->Fill(track.eta(), track.dz());
0945       static const int chiProbphiindex = this->index(vTrackProfiles_, "p_chi2Prob_vs_phi");
0946       vTrackProfiles_[chiProbphiindex]->Fill(track.phi(), chi2Prob);
0947       static const int chiProbabsd0index = this->index(vTrackProfiles_, "p_chi2Prob_vs_d0");
0948       vTrackProfiles_[chiProbabsd0index]->Fill(fabs(track.d0()), chi2Prob);
0949       static const int chiProbabsdzindex = this->index(vTrackProfiles_, "p_chi2Prob_vs_dz");
0950       vTrackProfiles_[chiProbabsdzindex]->Fill(track.dz(), chi2Prob);
0951       static const int chiphiindex = this->index(vTrackProfiles_, "p_chi2_vs_phi");
0952       vTrackProfiles_[chiphiindex]->Fill(track.phi(), track.chi2());
0953       static const int normchiphiindex = this->index(vTrackProfiles_, "p_normchi2_vs_phi");
0954       vTrackProfiles_[normchiphiindex]->Fill(track.phi(), normchi2);
0955       static const int chietaindex = this->index(vTrackProfiles_, "p_chi2_vs_eta");
0956       vTrackProfiles_[chietaindex]->Fill(track.eta(), track.chi2());
0957       static const int normchiptindex = this->index(vTrackProfiles_, "p_normchi2_vs_pt");
0958       vTrackProfiles_[normchiptindex]->Fill(track.pt(), normchi2);
0959       static const int normchipindex = this->index(vTrackProfiles_, "p_normchi2_vs_p");
0960       vTrackProfiles_[normchipindex]->Fill(track.p(), normchi2);
0961       static const int chiProbetaindex = this->index(vTrackProfiles_, "p_chi2Prob_vs_eta");
0962       vTrackProfiles_[chiProbetaindex]->Fill(track.eta(), chi2Prob);
0963       static const int normchietaindex = this->index(vTrackProfiles_, "p_normchi2_vs_eta");
0964       vTrackProfiles_[normchietaindex]->Fill(track.eta(), normchi2);
0965       static const int kappaphiindex = this->index(vTrackProfiles_, "p_kappa_vs_phi");
0966       vTrackProfiles_[kappaphiindex]->Fill(track.phi(), kappa);
0967       static const int kappaetaindex = this->index(vTrackProfiles_, "p_kappa_vs_eta");
0968       vTrackProfiles_[kappaetaindex]->Fill(track.eta(), kappa);
0969       static const int ptResphiindex = this->index(vTrackProfiles_, "p_ptResolution_vs_phi");
0970       vTrackProfiles_[ptResphiindex]->Fill(track.phi(), track.ptError() / track.pt());
0971       static const int ptResetaindex = this->index(vTrackProfiles_, "p_ptResolution_vs_eta");
0972       vTrackProfiles_[ptResetaindex]->Fill(track.eta(), track.ptError() / track.pt());
0973
0974       // Fill 2D track histos
0975       static const int etaphiindex_2d = this->index(vTrack2DHistos_, "h2_phi_vs_eta");
0976       vTrack2DHistos_[etaphiindex_2d]->Fill(track.eta(), track.phi());
0977       static const int d0phiindex_2d = this->index(vTrack2DHistos_, "h2_d0_vs_phi");
0978       vTrack2DHistos_[d0phiindex_2d]->Fill(track.phi(), track.d0());
0979       static const int dzphiindex_2d = this->index(vTrack2DHistos_, "h2_dz_vs_phi");
0980       vTrack2DHistos_[dzphiindex_2d]->Fill(track.phi(), track.dz());
0981       static const int d0etaindex_2d = this->index(vTrack2DHistos_, "h2_d0_vs_eta");
0982       vTrack2DHistos_[d0etaindex_2d]->Fill(track.eta(), track.d0());
0983       static const int dzetaindex_2d = this->index(vTrack2DHistos_, "h2_dz_vs_eta");
0984       vTrack2DHistos_[dzetaindex_2d]->Fill(track.eta(), track.dz());
0985       static const int chiphiindex_2d = this->index(vTrack2DHistos_, "h2_chi2_vs_phi");
0986       vTrack2DHistos_[chiphiindex_2d]->Fill(track.phi(), track.chi2());
0987       static const int chiProbphiindex_2d = this->index(vTrack2DHistos_, "h2_chi2Prob_vs_phi");
0988       vTrack2DHistos_[chiProbphiindex_2d]->Fill(track.phi(), chi2Prob);
0989       static const int chiProbabsd0index_2d = this->index(vTrack2DHistos_, "h2_chi2Prob_vs_d0");
0990       vTrack2DHistos_[chiProbabsd0index_2d]->Fill(fabs(track.d0()), chi2Prob);
0991       static const int normchiphiindex_2d = this->index(vTrack2DHistos_, "h2_normchi2_vs_phi");
0992       vTrack2DHistos_[normchiphiindex_2d]->Fill(track.phi(), normchi2);
0993       static const int chietaindex_2d = this->index(vTrack2DHistos_, "h2_chi2_vs_eta");
0994       vTrack2DHistos_[chietaindex_2d]->Fill(track.eta(), track.chi2());
0995       static const int chiProbetaindex_2d = this->index(vTrack2DHistos_, "h2_chi2Prob_vs_eta");
0996       vTrack2DHistos_[chiProbetaindex_2d]->Fill(track.eta(), chi2Prob);
0997       static const int normchietaindex_2d = this->index(vTrack2DHistos_, "h2_normchi2_vs_eta");
0998       vTrack2DHistos_[normchietaindex_2d]->Fill(track.eta(), normchi2);
0999       static const int kappaphiindex_2d = this->index(vTrack2DHistos_, "h2_kappa_vs_phi");
1000       vTrack2DHistos_[kappaphiindex_2d]->Fill(track.phi(), kappa);
1001       static const int kappaetaindex_2d = this->index(vTrack2DHistos_, "h2_kappa_vs_eta");
1002       vTrack2DHistos_[kappaetaindex_2d]->Fill(track.eta(), kappa);
1003       static const int normchi2kappa_2d = this->index(vTrack2DHistos_, "h2_normchi2_vs_kappa");
1004       vTrack2DHistos_[normchi2kappa_2d]->Fill(normchi2, kappa);
1005
1006       if (DEBUG)
1007         edm::LogVerbatim("DMRChecker") << "filling histos" << std::endl;
1008
1009       //dxy with respect to the beamspot
1010       reco::BeamSpot beamSpot;
1011       edm::Handle<reco::BeamSpot> beamSpotHandle = event.getHandle(beamspotToken_);
1012       if (beamSpotHandle.isValid()) {
1013         beamSpot = *beamSpotHandle;
1014         math::XYZPoint point(beamSpot.x0(), beamSpot.y0(), beamSpot.z0());
1015         double dxy = track.dxy(point);
1016         double dz = track.dz(point);
1017         hdxyBS->Fill(dxy);
1018         hd0BS->Fill(-dxy);
1019         hdzBS->Fill(dz);
1020       }
1021
1022       //dxy with respect to the primary vertex
1023       reco::Vertex pvtx;
1024       edm::Handle<reco::VertexCollection> vertexHandle = event.getHandle(vertexToken_);
1025       double mindxy = 100.;
1026       double dz = 100;
1027       if (vertexHandle.isValid() && !isCosmics_) {
1028         for (reco::VertexCollection::const_iterator pvtx = vertexHandle->begin(); pvtx != vertexHandle->end(); ++pvtx) {
1029           math::XYZPoint mypoint(pvtx->x(), pvtx->y(), pvtx->z());
1030           if (abs(mindxy) > abs(track.dxy(mypoint))) {
1031             mindxy = track.dxy(mypoint);
1032             dz = track.dz(mypoint);
1033             if (DEBUG)
1034               edm::LogVerbatim("DMRChecker") << "dxy: " << mindxy << " dz: " << dz << std::endl;
1035           }
1036         }
1037
1038         hdxyPV->Fill(mindxy);
1039         hd0PV->Fill(-mindxy);
1040         hdzPV->Fill(dz);
1041
1042         hd0PVvsphi->Fill(track.phi(), -mindxy);
1043         hd0PVvseta->Fill(track.eta(), -mindxy);
1044         hd0PVvspt->Fill(track.pt(), -mindxy);
1045
1046       } else {
1047         hdxyPV->Fill(100);
1048         hd0PV->Fill(100);
1049         hdzPV->Fill(100);
1050       }
1051
1052       if (DEBUG)
1053         edm::LogVerbatim("DMRChecker") << "end of track loop" << std::endl;
1054     }
1055
1056     if (DEBUG)
1057       edm::LogVerbatim("DMRChecker") << "end of analysis" << std::endl;
1058
1059     hNtrk->Fill(tC.size());
1060     hNtrkZoom->Fill(tC.size());
1061     hNhighPurity->Fill(nHighPurityTracks);
1062
1063     if (DEBUG)
1064       edm::LogVerbatim("DMRChecker") << "end of analysis" << std::endl;
1065   }
1066
1067   //*************************************************************
1068   void beginRun(edm::Run const &run, edm::EventSetup const &setup) override
1069   //*************************************************************
1070   {
1071     // initialize runInfoMap_
1072     runInfoMap_[run.run()].first = 0;
1073     runInfoMap_[run.run()].second = 0;
1074
1075     // initialize runHitsMap
1076     for (int n : {0, 1, 2, 3, 4, 5}) {
1077       runHitsMap_[run.run()][n] = 0;  // 6 subdets
1078     }
1079
1080     // Magnetic Field setup
1081     magneticField_ = &setup.getData(magFieldToken_);
1082     float B_ = magneticField_->inTesla(GlobalPoint(0, 0, 0)).mag();
1083
1084     edm::LogInfo("DMRChecker") << "run number:" << run.run() << " magnetic field: " << B_ << " [T]" << std::endl;
1085
1086     const RunInfo *summary = &setup.getData(runInfoToken_);
1087     time_t start_time = summary->m_start_time_ll;
1088     ctime(&start_time);
1089     time_t end_time = summary->m_stop_time_ll;
1090     ctime(&end_time);
1091
1092     float average_current = summary->m_avg_current;
1093     float uptimeInSeconds = summary->m_run_intervall_micros;
1094     edm::LogVerbatim("DMRChecker") << " start_time: " << start_time << " ( " << summary->m_start_time_str << " )"
1095                                    << " | end_time: " << end_time << " ( " << summary->m_stop_time_str << " )"
1096                                    << " | average current: " << average_current
1097                                    << " | uptime in seconds: " << uptimeInSeconds << std::endl;
1098
1099     double seconds = difftime(end_time, start_time) / 1.0e+6;  // convert from micros-seconds
1100     edm::LogVerbatim("DMRChecker") << "time difference: " << seconds << " s" << std::endl;
1101     timeMap_[run.run()] = seconds;
1102
1103     // set geometry and topology
1104     trackerGeometry_ = &setup.getData(geomToken_);
1105     if (trackerGeometry_->isThere(GeomDetEnumerators::P2PXB) || trackerGeometry_->isThere(GeomDetEnumerators::P2PXEC)) {
1106       phase_ = SiPixelPI::phase::two;
1107     } else if (trackerGeometry_->isThere(GeomDetEnumerators::P1PXB) ||
1108                trackerGeometry_->isThere(GeomDetEnumerators::P1PXEC)) {
1109       phase_ = SiPixelPI::phase::one;
1110     } else {
1111       phase_ = SiPixelPI::phase::zero;
1112     }
1113
1114     trackerTopology_ = &setup.getData(topoToken_);
1115
1116     // if it's a phase-2 geometry there are no phase-1 conditions
1117     if (phase_ == SiPixelPI::phase::two) {
1118       mode = 0;
1119     } else {
1120       if (doLatencyAnalysis_) {
1121         //SiStrip Latency
1122         const SiStripLatency *apvlat = &setup.getData(latencyToken_);
1123         if (apvlat->singleReadOutMode() == 1) {
1124           mode = 1;  // peak mode
1125         } else if (apvlat->singleReadOutMode() == 0) {
1126           mode = -1;  // deco mode
1127         }
1128       } else {
1129         mode = 0.;
1130       }
1131     }
1132
1133     conditionsMap_[run.run()].first = mode;
1134     conditionsMap_[run.run()].second = B_;
1135   }
1136
1137   //*************************************************************
1138   void endRun(edm::Run const &, edm::EventSetup const &) override {}
1139   //*************************************************************
1140
1141   void beginJob() override {
1142     if (DEBUG)
1143       edm::LogVerbatim("DMRChecker") << __LINE__ << endl;
1144
1145     TH1D::SetDefaultSumw2(kTRUE);
1146
1147     etaMax_ = 3.;  // assign max value to eta
1148
1149     // intialize counters
1150     ievt = 0;
1151     itrks = 0;
1152
1153     hrun = book<TH1D>("h_run", "run", 100000, 230000, 240000);
1154     hlumi = book<TH1D>("h_lumi", "lumi", 1000, 0, 1000);
1155
1156     // clang-format off
1157
1158     hchi2ndof = book<TH1D>("h_chi2ndof", "chi2/ndf;#chi^{2}/ndf;tracks", 100, 0, 5.);
1159     hCharge = book<TH1D>("h_charge", "charge;Charge of the track;tracks", 5, -2.5, 2.5);
1160     hNtrk = book<TH1D>("h_Ntrk", "ntracks;Number of Tracks;events", 200, 0., 200.);
1161     hNtrkZoom = book<TH1D>("h_NtrkZoom", "Number of tracks; number of tracks;events", 10, 0., 10.);
1162     hNhighPurity = book<TH1D>("h_NhighPurity", "n. high purity tracks;Number of high purity tracks;events", 200, 0., 200.);
1163
1164     int nAlgos = reco::TrackBase::algoSize;
1165     htrkAlgo = book<TH1I>("h_trkAlgo", "tracking step;iterative tracking step;tracks", nAlgos, -0.5, nAlgos - 0.5);
1166     for (int nbin = 1; nbin <= htrkAlgo->GetNbinsX(); nbin++) {
1167       htrkAlgo->GetXaxis()->SetBinLabel(nbin, reco::TrackBase::algoNames[nbin - 1].c_str());
1168     }
1169
1170     htrkQuality = book<TH1I>("h_trkQuality", "track quality;track quality;tracks", 6, -1, 5);
1171     std::string qualities[7] = {"undef", "loose", "tight", "highPurity", "confirmed", "goodIterative"};
1172     for (int nbin = 1; nbin <= htrkQuality->GetNbinsX(); nbin++) {
1173       htrkQuality->GetXaxis()->SetBinLabel(nbin, qualities[nbin - 1].c_str());
1174     }
1175
1176     hP = book<TH1D>("h_P", "Momentum;track momentum [GeV];tracks", 100, 0., 100.);
1177     hQoverP = book<TH1D>("h_qoverp", "Track q/p; track q/p [GeV^{-1}];tracks", 100, -1., 1.);
1178     hQoverPZoom = book<TH1D>("h_qoverpZoom", "Track q/p; track q/p [GeV^{-1}];tracks", 100, -0.1, 0.1);
1179     hPt = book<TH1D>("h_Pt", "Transverse Momentum;track p_{T} [GeV];tracks", 100, 0., 100.);
1180     hHit = book<TH1D>("h_nHits", "Number of hits;track n. hits;tracks", 50, -0.5, 49.5);
1181     hHit2D = book<TH1D>("h_nHit2D", "Number of 2D hits; number of 2D hits;tracks", 20, 0, 20);
1182
1183     // Pixel
1184
1185     hBPixResXPrime = book<TH1D>("h_BPixResXPrime", "BPix track X-residuals;res_{X'} [#mum];hits", 100, -1000., 1000.);
1186     hFPixResXPrime = book<TH1D>("h_FPixResXPrime", "FPix track X-residuals;res_{X'} [#mum];hits", 100, -1000., 1000.);
1187     hFPixZPlusResXPrime = book<TH1D>("h_FPixZPlusResXPrime", "FPix (Z+) track X-residuals;res_{X'} [#mum];hits", 100, -1000., 1000.);
1188     hFPixZMinusResXPrime = book<TH1D>("h_FPixZMinusResXPrime", "FPix (Z-) track X-residuals;res_{X'} [#mum];hits", 100, -1000., 1000.);
1189
1190     hBPixResYPrime = book<TH1D>("h_BPixResYPrime", "BPix track Y-residuals;res_{Y'} [#mum];hits", 100, -1000., 1000.);
1191     hFPixResYPrime = book<TH1D>("h_FPixResYPrime", "FPix track Y-residuals;res_{Y'} [#mum];hits", 100, -1000., 1000.);
1192     hFPixZPlusResYPrime = book<TH1D>("h_FPixZPlusResYPrime", "FPix (Z+) track Y-residuals;res_{Y'} [#mum];hits", 100, -1000., 1000.);
1193     hFPixZMinusResYPrime = book<TH1D>("h_FPixZMinusResYPrime", "FPix (Z-) track Y-residuals;res_{Y'} [#mum];hits", 100, -1000., 1000.);
1194
1195     hBPixResXPull = book<TH1D>("h_BPixResXPull", "BPix track X-pulls;res_{X'}/#sigma_{res_{X'}};hits", 100, -5., 5.);
1196     hFPixResXPull = book<TH1D>("h_FPixResXPull", "FPix track X-pulls;res_{X'}/#sigma_{res_{X'}};hits", 100, -5., 5.);
1197     hFPixZPlusResXPull = book<TH1D>("h_FPixZPlusResXPull", "FPix (Z+) track X-pulls;res_{X'}/#sigma_{res_{X'}};hits", 100, -5., 5.);
1198     hFPixZMinusResXPull = book<TH1D>("h_FPixZMinusResXPull", "FPix (Z-) track X-pulls;res_{X'}/#sigma_{res_{X'}};hits", 100, -5., 5.);
1199
1200     hBPixResYPull = book<TH1D>("h_BPixResYPull", "BPix track Y-pulls;res_{Y'}/#sigma_{res_{Y'}};hits", 100, -5., 5.);
1201     hFPixResYPull = book<TH1D>("h_FPixResYPull", "FPix track Y-pulls;res_{Y'}/#sigma_{res_{Y'}};hits", 100, -5., 5.);
1202     hFPixZPlusResYPull = book<TH1D>("h_FPixZPlusResYPull", "FPix (Z+) track Y-pulls;res_{Y'}/#sigma_{res_{Y'}};hits", 100, -5., 5.);
1203     hFPixZMinusResYPull = book<TH1D>("h_FPixZMinusResYPull", "FPix (Z-) track Y-pulls;res_{Y'}/#sigma_{res_{Y'}};hits", 100, -5., 5.);
1204
1205     // Strips
1206
1207     hTIBResXPrime = book<TH1D>("h_TIBResXPrime", "TIB track X-residuals;res_{X'} [#mum];hits", 100, -1000., 1000.);
1208     hTOBResXPrime = book<TH1D>("h_TOBResXPrime", "TOB track X-residuals;res_{X'} [#mum];hits", 100, -1000., 1000.);
1209     hTIDResXPrime = book<TH1D>("h_TIDResXPrime", "TID track X-residuals;res_{X'} [#mum];hits", 100, -1000., 1000.);
1210     hTECResXPrime = book<TH1D>("h_TECResXPrime", "TEC track X-residuals;res_{X'} [#mum];hits", 100, -1000., 1000.);
1211
1212     hTIBResXPull = book<TH1D>("h_TIBResXPull", "TIB track X-pulls;res_{X'}/#sigma_{res_{X'}};hits", 100, -5., 5.);
1213     hTOBResXPull = book<TH1D>("h_TOBResXPull", "TOB track X-pulls;res_{X'}/#sigma_{res_{X'}};hits", 100, -5., 5.);
1214     hTIDResXPull = book<TH1D>("h_TIDResXPull", "TID track X-pulls;res_{X'}/#sigma_{res_{X'}};hits", 100, -5., 5.);
1215     hTECResXPull = book<TH1D>("h_TECResXPull", "TEC track X-pulls;res_{X'}/#sigma_{res_{X'}};hits", 100, -5., 5.);
1216
1217     // hit counts
1218
1219     hHitCountVsZBPix = book<TH1D>("h_HitCountVsZBpix", "Number of BPix hits vs z;hit global z;hits", 60, -30, 30);
1220     hHitCountVsZFPix = book<TH1D>("h_HitCountVsZFpix", "Number of FPix hits vs z;hit global z;hits", 100, -100, 100);
1221
1222     hHitCountVsXBPix = book<TH1D>("h_HitCountVsXBpix", "Number of BPix hits vs x;hit global x;hits", 20, -20, 20);
1223     hHitCountVsXFPix = book<TH1D>("h_HitCountVsXFpix", "Number of FPix hits vs x;hit global x;hits", 20, -20, 20);
1224
1225     hHitCountVsYBPix = book<TH1D>("h_HitCountVsYBpix", "Number of BPix hits vs y;hit global y;hits", 20, -20, 20);
1226     hHitCountVsYFPix = book<TH1D>("h_HitCountVsYFpix", "Number of FPix hits vs y;hit global y;hits", 20, -20, 20);
1227
1228     hHitCountVsThetaBPix = book<TH1D>("h_HitCountVsThetaBpix", "Number of BPix hits vs #theta;hit global #theta;hits", 20, 0., M_PI);
1229     hHitCountVsPhiBPix = book<TH1D>("h_HitCountVsPhiBpix", "Number of BPix hits vs #phi;hit global #phi;hits", 20, -M_PI, M_PI);
1230
1231     hHitCountVsThetaFPix = book<TH1D>("h_HitCountVsThetaFpix", "Number of FPix hits vs #theta;hit global #theta;hits", 40, 0., M_PI);
1232     hHitCountVsPhiFPix = book<TH1D>("h_HitCountVsPhiFpix", "Number of FPix hits vs #phi;hit global #phi;hits", 20, -M_PI, M_PI);
1233
1234     // two sides of FPix
1235
1236     hHitCountVsZFPixPlus = book<TH1D>("h_HitCountVsZFPixPlus", "Number of FPix(Z+) hits vs z;hit global z;hits", 60, 15., 60);
1237     hHitCountVsZFPixMinus = book<TH1D>("h_HitCountVsZFPixMinus", "Number of FPix(Z-) hits vs z;hit global z;hits", 100, -60., -15.);
1238
1239     hHitCountVsXFPixPlus = book<TH1D>("h_HitCountVsXFPixPlus", "Number of FPix(Z+) hits vs x;hit global x;hits", 20, -20, 20);
1240     hHitCountVsXFPixMinus = book<TH1D>("h_HitCountVsXFPixMinus", "Number of FPix(Z-) hits vs x;hit global x;hits", 20, -20, 20);
1241
1242     hHitCountVsYFPixPlus = book<TH1D>("h_HitCountVsYFPixPlus", "Number of FPix(Z+) hits vs y;hit global y;hits", 20, -20, 20);
1243     hHitCountVsYFPixMinus = book<TH1D>("h_HitCountVsYFPixMinus", "Number of FPix(Z-) hits vs y;hit global y;hits", 20, -20, 20);
1244
1245     hHitCountVsThetaFPixPlus = book<TH1D>("h_HitCountVsThetaFPixPlus", "Number of FPix(Z+) hits vs #theta;hit global #theta;hits", 20, 0., M_PI);
1246     hHitCountVsPhiFPixPlus = book<TH1D>("h_HitCountVsPhiFPixPlus","Number of FPix(Z+) hits vs #phi;hit global #phi;hits",20,-M_PI,M_PI);
1247
1248     hHitCountVsThetaFPixMinus = book<TH1D>("h_HitCountVsThetaFPixMinus", "Number of FPix(Z+) hits vs #theta;hit global #theta;hits", 40, 0., M_PI);
1249     hHitCountVsPhiFPixMinus = book<TH1D>("h_HitCountVsPhiFPixMinus","Number of FPix(Z+) hits vs #phi;hit global #phi;hits",20,-M_PI,M_PI);
1250
1251     TFileDirectory ByLayerResiduals = fs->mkdir("ByLayerResiduals");
1252     barrelLayersResidualsX = bookResidualsHistogram(ByLayerResiduals, 4, "X", "Res", "BPix");
1253     endcapDisksResidualsX = bookResidualsHistogram(ByLayerResiduals, 6, "X", "Res", "FPix");
1254     barrelLayersResidualsY = bookResidualsHistogram(ByLayerResiduals, 4, "Y", "Res", "BPix");
1255     endcapDisksResidualsY = bookResidualsHistogram(ByLayerResiduals, 6, "Y", "Res", "FPix");
1256
1257     TFileDirectory ByLayerPulls = fs->mkdir("ByLayerPulls");
1258     barrelLayersPullsX = bookResidualsHistogram(ByLayerPulls, 4, "X", "Pull", "BPix");
1259     endcapDisksPullsX = bookResidualsHistogram(ByLayerPulls, 6, "X", "Pull", "FPix");
1260     barrelLayersPullsY = bookResidualsHistogram(ByLayerPulls, 4, "Y", "Pull", "BPix");
1261     endcapDisksPullsY = bookResidualsHistogram(ByLayerPulls, 6, "Y", "Pull", "FPix");
1262
1263     hEta = book<TH1D>("h_Eta", "Track pseudorapidity; track #eta;tracks", 100, -etaMax_, etaMax_);
1264     hPhi = book<TH1D>("h_Phi", "Track azimuth; track #phi;tracks", 100, -M_PI, M_PI);
1265
1266     hPhiBarrel = book<TH1D>("h_PhiBarrel", "hPhiBarrel (0<|#eta|<0.8);track #Phi;tracks", 100, -M_PI, M_PI);
1267     hPhiOverlapPlus = book<TH1D>("h_PhiOverlapPlus", "hPhiOverlapPlus (0.8<#eta<1.4);track #phi;tracks", 100, -M_PI, M_PI);
1268     hPhiOverlapMinus = book<TH1D>("h_PhiOverlapMinus", "hPhiOverlapMinus (-1.4<#eta<-0.8);track #phi;tracks", 100, -M_PI, M_PI);
1269     hPhiEndcapPlus = book<TH1D>("h_PhiEndcapPlus", "hPhiEndcapPlus (#eta>1.4);track #phi;track", 100, -M_PI, M_PI);
1270     hPhiEndcapMinus = book<TH1D>("h_PhiEndcapMinus", "hPhiEndcapMinus (#eta<1.4);track #phi;tracks", 100, -M_PI, M_PI);
1271
1272     if (!isCosmics_) {
1273       hPhp = book<TH1D>("h_P_hp", "Momentum (high purity);track momentum [GeV];tracks", 100, 0., 100.);
1274       hPthp = book<TH1D>("h_Pt_hp", "Transverse Momentum (high purity);track p_{T} [GeV];tracks", 100, 0., 100.);
1275       hHithp = book<TH1D>("h_nHit_hp", "Number of hits (high purity);track n. hits;tracks", 30, 0, 30);
1276       hEtahp = book<TH1D>("h_Eta_hp", "Track pseudorapidity (high purity); track #eta;tracks", 100, -etaMax_, etaMax_);
1277       hPhihp = book<TH1D>("h_Phi_hp", "Track azimuth (high purity); track #phi;tracks", 100, -M_PI, M_PI);
1278       hchi2ndofhp = book<TH1D>("h_chi2ndof_hp", "chi2/ndf (high purity);#chi^{2}/ndf;tracks", 100, 0, 5.);
1279       hchi2Probhp = book<TH1D>("hchi2_Prob_hp", "#chi^{2} probability (high purity);#chi^{2}prob_{Track};Number of Tracks", 100, 0.0, 1.);
1280
1281       hvx = book<TH1D>("h_vx", "Track v_{x} ; track v_{x} [cm];tracks", 100, -1.5, 1.5);
1282       hvy = book<TH1D>("h_vy", "Track v_{y} ; track v_{y} [cm];tracks", 100, -1.5, 1.5);
1283       hvz = book<TH1D>("h_vz", "Track v_{z} ; track v_{z} [cm];tracks", 100, -20., 20.);
1284       hd0 = book<TH1D>("h_d0", "Track d_{0} ; track d_{0} [cm];tracks", 100, -1., 1.);
1285       hdxy = book<TH1D>("h_dxy", "Track d_{xy}; track d_{xy} [cm]; tracks", 100, -0.5, 0.5);
1286       hdz = book<TH1D>("h_dz", "Track d_{z} ; track d_{z} [cm]; tracks", 100, -20, 20);
1287
1288       hd0PVvsphi = book<TH2D>("h2_d0PVvsphi", "hd0PVvsphi;track #phi;track d_{0}(PV) [cm]", 160, -M_PI, M_PI, 100, -1., 1.);
1289       hd0PVvseta = book<TH2D>("h2_d0PVvseta", "hdPV0vseta;track #eta;track d_{0}(PV) [cm]", 160, -etaMax_, etaMax_, 100, -1., 1.);
1290       hd0PVvspt = book<TH2D>("h2_d0PVvspt", "hdPV0vspt;track p_{T};d_{0}(PV) [cm]", 50, 0., 100., 100, -1, 1.);
1291
1292       hdxyBS = book<TH1D>("h_dxyBS", "hdxyBS; track d_{xy}(BS) [cm];tracks", 100, -0.1, 0.1);
1293       hd0BS = book<TH1D>("h_d0BS", "hd0BS ; track d_{0}(BS) [cm];tracks", 100, -0.1, 0.1);
1294       hdzBS = book<TH1D>("h_dzBS", "hdzBS ; track d_{z}(BS) [cm];tracks", 100, -12, 12);
1295       hdxyPV = book<TH1D>("h_dxyPV", "hdxyPV; track d_{xy}(PV) [cm];tracks", 100, -0.1, 0.1);
1296       hd0PV = book<TH1D>("h_d0PV", "hd0PV ; track d_{0}(PV) [cm];tracks", 100, -0.15, 0.15);
1297       hdzPV = book<TH1D>("h_dzPV", "hdzPV ; track d_{z}(PV) [cm];tracks", 100, -0.1, 0.1);
1298
1299       hnhTIB = book<TH1D>("h_nHitTIB", "nhTIB;# hits in TIB; tracks", 20, 0., 20.);
1300       hnhTID = book<TH1D>("h_nHitTID", "nhTID;# hits in TID; tracks", 20, 0., 20.);
1301       hnhTOB = book<TH1D>("h_nHitTOB", "nhTOB;# hits in TOB; tracks", 20, 0., 20.);
1302       hnhTEC = book<TH1D>("h_nHitTEC", "nhTEC;# hits in TEC; tracks", 20, 0., 20.);
1303
1304     } else {
1305       hvx = book<TH1D>("h_vx", "Track v_{x};track v_{x} [cm];tracks", 100, -100., 100.);
1306       hvy = book<TH1D>("h_vy", "Track v_{y};track v_{y} [cm];tracks", 100, -100., 100.);
1307       hvz = book<TH1D>("h_vz", "Track v_{z};track v_{z} [cm];track", 100, -100., 100.);
1308       hd0 = book<TH1D>("h_d0", "Track d_{0};track d_{0} [cm];track", 100, -100., 100.);
1309       hdxy = book<TH1D>("h_dxy", "Track d_{xy};track d_{xy} [cm];tracks", 100, -100, 100);
1310       hdz = book<TH1D>("h_dz", "Track d_{z};track d_{z} [cm];tracks", 100, -200, 200);
1311
1312       hd0vsphi = book<TH2D>("h2_d0vsphi", "Track d_{0} vs #phi; track #phi;track d_{0} [cm]", 160, -M_PI, M_PI, 100, -100., 100.);
1313       hd0vseta = book<TH2D>("h2_d0vseta", "Track d_{0} vs #eta; track #eta;track d_{0} [cm]", 160, -etaMax_, etaMax_, 100, -100., 100.);
1314       hd0vspt = book<TH2D>("h2_d0vspt", "Track d_{0} vs p_{T};track p_{T};track d_{0} [cm]", 50, 0., 100., 100, -100, 100);
1315
1316       hdxyBS = book<TH1D>("h_dxyBS", "Track d_{xy}(BS);d_{xy}(BS) [cm];tracks", 100, -100., 100.);
1317       hd0BS = book<TH1D>("h_d0BS", "Track d_{0}(BS);d_{0}(BS) [cm];tracks", 100, -100., 100.);
1318       hdzBS = book<TH1D>("h_dzBS", "Track d_{z}(BS);d_{z}(BS) [cm];tracks", 100, -100., 100.);
1319       hdxyPV = book<TH1D>("h_dxyPV", "Track d_{xy}(PV); d_{xy}(PV) [cm];tracks", 100, -100., 100.);
1320       hd0PV = book<TH1D>("h_d0PV", "Track d_{0}(PV); d_{0}(PV) [cm];tracks", 100, -100., 100.);
1321       hdzPV = book<TH1D>("h_dzPV", "Track d_{z}(PV); d_{z}(PV) [cm];tracks", 100, -100., 100.);
1322
1323       hnhTIB = book<TH1D>("h_nHitTIB", "nhTIB;# hits in TIB; tracks", 30, 0., 30.);
1324       hnhTID = book<TH1D>("h_nHitTID", "nhTID;# hits in TID; tracks", 30, 0., 30.);
1325       hnhTOB = book<TH1D>("h_nHitTOB", "nhTOB;# hits in TOB; tracks", 30, 0., 30.);
1326       hnhTEC = book<TH1D>("h_nHitTEC", "nhTEC;# hits in TEC; tracks", 30, 0., 30.);
1327     }
1328
1329     hnhpxb = book<TH1D>("h_nHitPXB", "nhpxb;# hits in Pixel Barrel; tracks", 10, 0., 10.);
1330     hnhpxe = book<TH1D>("h_nHitPXE", "nhpxe;# hits in Pixel Endcap; tracks", 10, 0., 10.);
1331
1332     hHitComposition = book<TH1D>("h_hitcomposition", "track hit composition;;# hits", 6, -0.5, 5.5);
1333     pNBpixHitsVsVx = book<TProfile>("p_NpixHits_vs_Vx", "n. Barrel Pixel hits vs. v_{x};v_{x} (cm);n. BPix hits", 20, -20, 20);
1334     pNBpixHitsVsVy = book<TProfile>("p_NpixHits_vs_Vy", "n. Barrel Pixel hits vs. v_{y};v_{y} (cm);n. BPix hits", 20, -20, 20);
1335     pNBpixHitsVsVz = book<TProfile>("p_NpixHits_vs_Vz", "n. Barrel Pixel hits vs. v_{z};v_{z} (cm);n. BPix hits", 20, -100, 100);
1336
1337     std::string dets[6] = {"PXB", "PXF", "TIB", "TID", "TOB", "TEC"};
1338     for (int i = 1; i <= hHitComposition->GetNbinsX(); i++) {
1339       hHitComposition->GetXaxis()->SetBinLabel(i, dets[i - 1].c_str());
1340     }
1341
1342     vTrackHistos_.push_back(book<TH1F>("h_tracketa", "Track #eta;#eta_{Track};Number of Tracks", 90, -etaMax_, etaMax_));
1343     vTrackHistos_.push_back(book<TH1F>("h_trackphi", "Track #phi;#phi_{Track};Number of Tracks", 90, -M_PI, M_PI));
1344     vTrackHistos_.push_back(book<TH1F>("h_trackNumberOfValidHits", "Track # of valid hits;# of valid hits _{Track};Number of Tracks", 40, 0., 40.));
1345     vTrackHistos_.push_back(book<TH1F>("h_trackNumberOfLostHits", "Track # of lost hits;# of lost hits _{Track};Number of Tracks", 10, 0., 10.));
1346     vTrackHistos_.push_back(book<TH1F>("h_curvature", "Curvature #kappa;#kappa_{Track};Number of Tracks", 100, -.05, .05));
1347     vTrackHistos_.push_back(book<TH1F>("h_curvature_pos", "Curvature |#kappa| Positive Tracks;|#kappa_{pos Track}|;Number of Tracks", 100, .0, .05));
1348     vTrackHistos_.push_back(book<TH1F>("h_curvature_neg", "Curvature |#kappa| Negative Tracks;|#kappa_{neg Track}|;Number of Tracks", 100, .0, .05));
1349     vTrackHistos_.push_back(book<TH1F>("h_diff_curvature", "Curvature |#kappa| Tracks Difference;|#kappa_{Track}|;# Pos Tracks - # Neg Tracks", 100,.0,.05));
1350
1351     vTrackHistos_.push_back(book<TH1F>("h_chi2", "Track #chi^{2};#chi^{2}_{Track};Number of Tracks", 500, -0.01, 500.));
1352     vTrackHistos_.push_back(book<TH1F>("h_chi2Prob", "#chi^{2} probability;Track Prob(#chi^{2},ndof);Number of Tracks", 100, 0.0, 1.));
1353     vTrackHistos_.push_back(book<TH1F>("h_normchi2", "#chi^{2}/ndof;#chi^{2}/ndof;Number of Tracks", 100, -0.01, 10.));
1354
1355     //variable binning for chi2/ndof vs. pT
1356     double xBins[19] = {0., 0.15, 0.5, 1., 1.5, 2., 2.5, 3., 3.5, 4., 4.5, 5., 7., 10., 15., 25., 40., 100., 200.};
1357     vTrackHistos_.push_back(book<TH1F>("h_pt", "Track p_{T};p_{T}^{track} [GeV];Number of Tracks", 250, 0., 250));
1358     vTrackHistos_.push_back(book<TH1F>("h_ptrebin", "Track p_{T};p_{T}^{track} [GeV];Number of Tracks", 18, xBins));
1359
1360     vTrackHistos_.push_back(book<TH1F>("h_ptResolution", "#delta_{p_{T}}/p_{T}^{track};#delta_{p_{T}}/p_{T}^{track};Number of Tracks", 100, 0., 0.5));
1361     vTrackProfiles_.push_back(book<TProfile>("p_d0_vs_phi", "Transverse Impact Parameter vs. #phi;#phi_{Track};#LT d_{0} #GT [cm]", 100, -M_PI, M_PI));
1362     vTrackProfiles_.push_back(book<TProfile>("p_dz_vs_phi", "Longitudinal Impact Parameter vs. #phi;#phi_{Track};#LT d_{z} #GT [cm]", 100, -M_PI, M_PI));
1363     vTrackProfiles_.push_back(book<TProfile>("p_d0_vs_eta", "Transverse Impact Parameter vs. #eta;#eta_{Track};#LT d_{0} #GT [cm]", 100, -etaMax_, etaMax_));
1364     vTrackProfiles_.push_back(book<TProfile>("p_dz_vs_eta", "Longitudinal Impact Parameter vs. #eta;#eta_{Track};#LT d_{z} #GT [cm]", 100, -etaMax_, etaMax_));
1365     vTrackProfiles_.push_back(book<TProfile>("p_chi2_vs_phi", "#chi^{2} vs. #phi;#phi_{Track};#LT #chi^{2} #GT", 100, -M_PI, M_PI));
1366     vTrackProfiles_.push_back(book<TProfile>("p_chi2Prob_vs_phi","#chi^{2} probablility vs. #phi;#phi_{Track};#LT #chi^{2} probability#GT",100,-M_PI,M_PI));
1367     vTrackProfiles_.push_back(book<TProfile>("p_chi2Prob_vs_d0", "#chi^{2} probablility vs. |d_{0}|;|d_{0}|[cm];#LT #chi^{2} probability#GT", 100, 0, 80));
1368     vTrackProfiles_.push_back(book<TProfile>("p_chi2Prob_vs_dz", "#chi^{2} probablility vs. dz;d_{z} [cm];#LT #chi^{2} probability#GT", 100, -30, 30));
1369     vTrackProfiles_.push_back(book<TProfile>("p_normchi2_vs_phi", "#chi^{2}/ndof vs. #phi;#phi_{Track};#LT #chi^{2}/ndof #GT", 100, -M_PI, M_PI));
1370     vTrackProfiles_.push_back(book<TProfile>("p_chi2_vs_eta", "#chi^{2} vs. #eta;#eta_{Track};#LT #chi^{2} #GT", 100, -etaMax_, etaMax_));
1371     vTrackProfiles_.push_back(book<TProfile>("p_normchi2_vs_pt", "norm #chi^{2} vs. p_{T}_{Track}; p_{T}_{Track};#LT #chi^{2}/ndof #GT", 18, xBins));
1372     vTrackProfiles_.push_back(book<TProfile>("p_normchi2_vs_p", "#chi^{2}/ndof vs. p_{Track};p_{Track};#LT #chi^{2}/ndof #GT", 18, xBins));
1373     vTrackProfiles_.push_back(book<TProfile>("p_chi2Prob_vs_eta","#chi^{2} probability vs. #eta;#eta_{Track};#LT #chi^{2} probability #GT",100,-etaMax_, etaMax_));
1374     vTrackProfiles_.push_back(book<TProfile>("p_normchi2_vs_eta", "#chi^{2}/ndof vs. #eta;#eta_{Track};#LT #chi^{2}/ndof #GT", 100, -etaMax_, etaMax_));
1375     vTrackProfiles_.push_back(book<TProfile>("p_kappa_vs_phi", "#kappa vs. #phi;#phi_{Track};#kappa", 100, -M_PI, M_PI));
1376     vTrackProfiles_.push_back(book<TProfile>("p_kappa_vs_eta", "#kappa vs. #eta;#eta_{Track};#kappa", 100, -etaMax_, etaMax_));
1377     vTrackProfiles_.push_back(book<TProfile>("p_ptResolution_vs_phi","#delta_{p_{T}}/p_{T}^{track};#phi^{track};#delta_{p_{T}}/p_{T}^{track}", 100,-M_PI,M_PI));
1378     vTrackProfiles_.push_back(book<TProfile>("p_ptResolution_vs_eta","#delta_{p_{T}}/p_{T}^{track};#eta^{track};#delta_{p_{T}}/p_{T}^{track}", 100, -etaMax_, etaMax_));
1379     vTrack2DHistos_.push_back(book<TH2F>("h2_d0_vs_phi","Transverse Impact Parameter vs. #phi;#phi_{Track};d_{0} [cm]", 100, -M_PI, M_PI, 100, -1., 1.));
1380     vTrack2DHistos_.push_back(book<TH2F>("h2_phi_vs_eta", "Track #phi vs. #eta;#eta_{Track};#phi_{Track}",50, -etaMax_, etaMax_, 50, -M_PI, M_PI));
1381     vTrack2DHistos_.push_back(book<TH2F>("h2_dz_vs_phi","Longitudinal Impact Parameter vs. #phi;#phi_{Track};d_{z} [cm]",100,-M_PI,M_PI,100,-100.,100.));
1382     vTrack2DHistos_.push_back(book<TH2F>("h2_d0_vs_eta","Transverse Impact Parameter vs. #eta;#eta_{Track};d_{0} [cm]", 100, -etaMax_, etaMax_, 100, -1., 1.));
1383     vTrack2DHistos_.push_back(book<TH2F>("h2_dz_vs_eta","Longitudinal Impact Parameter vs. #eta;#eta_{Track};d_{z} [cm]",100, -etaMax_, etaMax_, 100,-100.,100.));
1384     vTrack2DHistos_.push_back(book<TH2F>("h2_chi2_vs_phi", "#chi^{2} vs. #phi;#phi_{Track};#chi^{2}", 100, -M_PI, M_PI, 500, 0., 500.));
1385     vTrack2DHistos_.push_back(book<TH2F>("h2_chi2Prob_vs_phi","#chi^{2} probability vs. #phi;#phi_{Track};#chi^{2} probability",100,-M_PI,M_PI,100,0.,1.));
1386     vTrack2DHistos_.push_back(book<TH2F>("h2_chi2Prob_vs_d0","#chi^{2} probability vs. |d_{0}|;|d_{0}| [cm];#chi^{2} probability",100,0,80,100,0.,1.));
1387     vTrack2DHistos_.push_back(book<TH2F>("h2_normchi2_vs_phi", "#chi^{2}/ndof vs. #phi;#phi_{Track};#chi^{2}/ndof", 100, -M_PI, M_PI, 100, 0., 10.));
1388     vTrack2DHistos_.push_back(book<TH2F>("h2_chi2_vs_eta", "#chi^{2} vs. #eta;#eta_{Track};#chi^{2}", 100, -etaMax_, etaMax_, 500, 0., 500.));
1389     vTrack2DHistos_.push_back(book<TH2F>("h2_chi2Prob_vs_eta","#chi^{2} probaility vs. #eta;#eta_{Track};#chi^{2} probability",100,-M_PI,M_PI,100,0.,1.));
1390     vTrack2DHistos_.push_back(book<TH2F>("h2_normchi2_vs_eta", "#chi^{2}/ndof vs. #eta;#eta_{Track};#chi^{2}/ndof", 100, -etaMax_, etaMax_, 100, 0., 10.));
1391     vTrack2DHistos_.push_back(book<TH2F>("h2_kappa_vs_phi", "#kappa vs. #phi;#phi_{Track};#kappa", 100, -M_PI, M_PI, 100, .0, .05));
1392     vTrack2DHistos_.push_back(book<TH2F>("h2_kappa_vs_eta", "#kappa vs. #eta;#eta_{Track};#kappa", 100, -etaMax_, etaMax_, 100, .0, .05));
1393     vTrack2DHistos_.push_back(book<TH2F>("h2_normchi2_vs_kappa", "#kappa vs. #chi^{2}/ndof;#chi^{2}/ndof;#kappa", 100, 0., 10, 100, -.03, .03));
1394     // clang-format on
1395
1396     // create the full maps
1397     fullPixelmapXDMR->createTrackerBaseMap();
1398     fullPixelmapYDMR->createTrackerBaseMap();
1399
1400   }  //beginJob
1401
1402   void endJob() override {
1403     edm::LogPrint("DMRChecker") << "*******************************" << std::endl;
1404     edm::LogPrint("DMRChecker") << "Events run in total: " << ievt << std::endl;
1405     edm::LogPrint("DMRChecker") << "n. tracks: " << itrks << std::endl;
1406     edm::LogPrint("DMRChecker") << "*******************************" << std::endl;
1407
1408     int nFiringTriggers = !triggerMap_.empty() ? triggerMap_.size() : 1;
1409     edm::LogPrint("DMRChecker") << "firing triggers: " << triggerMap_.size() << std::endl;
1410     edm::LogPrint("DMRChecker") << "*******************************" << std::endl;
1411
1412     tksByTrigger_ =
1413         book<TH1D>("tksByTrigger", "tracks by HLT path;;% of # traks", nFiringTriggers, -0.5, nFiringTriggers - 0.5);
1414     evtsByTrigger_ =
1415         book<TH1D>("evtsByTrigger", "events by HLT path;;% of # events", nFiringTriggers, -0.5, nFiringTriggers - 0.5);
1416
1417     if (DEBUG)
1418       edm::LogPrint("DMRChecker") << __FILE__ << "@" << __FUNCTION__ << " L-" << __LINE__ << std::endl;
1419
1420     int i = 0;
1421     for (const auto &it : triggerMap_) {
1422       i++;
1423
1424       double trkpercent = ((it.second).second) * 100. / double(itrks);
1425       double evtpercent = ((it.second).first) * 100. / double(ievt);
1426
1427       std::cout.precision(4);
1428
1429       edm::LogPrint("DMRChecker") << "HLT path: " << std::setw(60) << left << it.first << " | events firing: " << right
1430                                   << std::setw(8) << (it.second).first << " (" << setw(8) << fixed << evtpercent << "%)"
1431                                   << " | tracks collected: " << std::setw(10) << (it.second).second << " (" << setw(8)
1432                                   << fixed << trkpercent << "%)";
1433
1434       tksByTrigger_->SetBinContent(i, trkpercent);
1435       tksByTrigger_->GetXaxis()->SetBinLabel(i, (it.first).c_str());
1436
1437       evtsByTrigger_->SetBinContent(i, evtpercent);
1438       evtsByTrigger_->GetXaxis()->SetBinLabel(i, (it.first).c_str());
1439     }
1440
1441     if (DEBUG)
1442       edm::LogPrint("DMRChecker") << __FILE__ << "@" << __FUNCTION__ << " L-" << __LINE__ << std::endl;
1443
1444     int nRuns = conditionsMap_.size();
1445     if (nRuns < 1)
1446       return;
1447
1448     vector<int> theRuns_;
1449     theRuns_.reserve(conditionsMap_.size());
1450     for (const auto &it : conditionsMap_) {
1451       theRuns_.push_back(it.first);
1452     }
1453
1454     std::sort(theRuns_.begin(), theRuns_.end());
1455     int runRange = theRuns_.back() - theRuns_.front() + 1;
1456
1457     edm::LogPrint("DMRChecker") << "*******************************" << std::endl;
1458     edm::LogPrint("DMRChecker") << "first run: " << theRuns_.front() << std::endl;
1459     edm::LogPrint("DMRChecker") << "last run:  " << theRuns_.back() << std::endl;
1460     edm::LogPrint("DMRChecker") << "considered runs: " << nRuns << std::endl;
1461     edm::LogPrint("DMRChecker") << "*******************************" << std::endl;
1462
1463     modeByRun_ = book<TH1D>("modeByRun",
1464                             "Strip APV mode by run number;;APV mode (-1=deco,+1=peak)",
1465                             runRange,
1466                             theRuns_.front() - 0.5,
1467                             theRuns_.back() + 0.5);
1468     fieldByRun_ = book<TH1D>("fieldByRun",
1469                              "CMS B-field intensity by run number;;B-field intensity [T]",
1470                              runRange,
1471                              theRuns_.front() - 0.5,
1472                              theRuns_.back() + 0.5);
1473
1474     tracksByRun_ = book<TH1D>("tracksByRun",
1475                               "n. AlCaReco Tracks by run number;;n. of tracks",
1476                               runRange,
1477                               theRuns_.front() - 0.5,
1478                               theRuns_.back() + 0.5);
1479     hitsByRun_ = book<TH1D>(
1480         "histByRun", "n. of hits by run number;;n. of hits", runRange, theRuns_.front() - 0.5, theRuns_.back() + 0.5);
1481
1482     trackRatesByRun_ = book<TH1D>("trackRatesByRun",
1483                                   "rate of AlCaReco Tracks by run number;;n. of tracks/s",
1484                                   runRange,
1485                                   theRuns_.front() - 0.5,
1486                                   theRuns_.back() + 0.5);
1487     eventRatesByRun_ = book<TH1D>("eventRatesByRun",
1488                                   "rate of AlCaReco Events by run number;;n. of events/s",
1489                                   runRange,
1490                                   theRuns_.front() - 0.5,
1491                                   theRuns_.back() + 0.5);
1492
1493     hitsinBPixByRun_ = book<TH1D>("histinBPixByRun",
1494                                   "n. of hits in BPix by run number;;n. of BPix hits",
1495                                   runRange,
1496                                   theRuns_.front() - 0.5,
1497                                   theRuns_.back() + 0.5);
1498     hitsinFPixByRun_ = book<TH1D>("histinFPixByRun",
1499                                   "n. of hits in FPix by run number;;n. of FPix hits",
1500                                   runRange,
1501                                   theRuns_.front() - 0.5,
1502                                   theRuns_.back() + 0.5);
1503
1504     for (const auto &the_r : theRuns_) {
1505       if (conditionsMap_.find(the_r)->second.first != 0) {
1506         auto indexing = (the_r - theRuns_.front()) + 1;
1507         double runTime = timeMap_.find(the_r)->second;
1508
1509         edm::LogPrint("DMRChecker") << "run:" << the_r << " | isPeak: " << std::setw(4)
1510                                     << conditionsMap_.find(the_r)->second.first
1511                                     << "| B-field: " << conditionsMap_.find(the_r)->second.second << " [T]"
1512                                     << "| events: " << setw(10) << runInfoMap_.find(the_r)->second.first
1513                                     << "(rate: " << setw(10) << (runInfoMap_.find(the_r)->second.first) / runTime
1514                                     << " ev/s)"
1515                                     << ", tracks " << setw(10) << runInfoMap_.find(the_r)->second.second
1516                                     << "(rate: " << setw(10) << (runInfoMap_.find(the_r)->second.second) / runTime
1517                                     << " trk/s)" << std::endl;
1518
1519         // int the_bin = modeByRun_->GetXaxis()->FindBin(the_r);
1520         modeByRun_->SetBinContent(indexing, conditionsMap_.find(the_r)->second.first);
1521         modeByRun_->GetXaxis()->SetBinLabel(indexing, Form("%d", the_r));
1522         fieldByRun_->SetBinContent(indexing, conditionsMap_.find(the_r)->second.second);
1523         fieldByRun_->GetXaxis()->SetBinLabel(indexing, Form("%d", the_r));
1524
1525         tracksByRun_->SetBinContent(indexing, runInfoMap_.find(the_r)->second.first);
1526         tracksByRun_->GetXaxis()->SetBinLabel(indexing, Form("%d", the_r));
1527         hitsByRun_->SetBinContent(indexing, runInfoMap_.find(the_r)->second.second);
1528         hitsByRun_->GetXaxis()->SetBinLabel(indexing, Form("%d", the_r));
1529
1530         hitsinBPixByRun_->SetBinContent(indexing, (runHitsMap_.find(the_r)->second)[0]);
1531         hitsinBPixByRun_->GetXaxis()->SetBinLabel(indexing, Form("%d", the_r));
1532         hitsinFPixByRun_->SetBinContent(indexing, (runHitsMap_.find(the_r)->second)[1]);
1533         hitsinFPixByRun_->GetXaxis()->SetBinLabel(indexing, Form("%d", the_r));
1534
1535         trackRatesByRun_->SetBinContent(indexing, (runInfoMap_.find(the_r)->second.second) / runTime);
1536         trackRatesByRun_->GetXaxis()->SetBinLabel(indexing, Form("%d", the_r));
1537         eventRatesByRun_->SetBinContent(indexing, (runInfoMap_.find(the_r)->second.first) / runTime);
1538         eventRatesByRun_->GetXaxis()->SetBinLabel(indexing, Form("%d", the_r));
1539
1540         constexpr const char *subdets[]{"BPix", "FPix", "TIB", "TID", "TOB", "TEC"};
1541
1542         edm::LogPrint("DMRChecker") << "*******************************" << std::endl;
1543         edm::LogPrint("DMRChecker") << "Hits by SubDetector" << std::endl;
1544         int si = 0;
1545         for (const auto &entry : runHitsMap_.find(the_r)->second) {
1546           edm::LogPrint("DMRChecker") << subdets[si] << " " << entry << std::endl;
1547           si++;
1548         }
1549         edm::LogPrint("DMRChecker") << "*******************************" << std::endl;
1550       }
1551
1552       // modeByRun_->GetXaxis()->SetBinLabel(the_r-theRuns_[0]+1,(const char*)the_r);
1553     }
1554
1555     if (DEBUG)
1556       edm::LogPrint("DMRChecker") << __FILE__ << "@" << __FUNCTION__ << " L-" << __LINE__ << std::endl;
1557
1558     // DMRs
1559
1560     TFileDirectory DMeanR = fs->mkdir("DMRs");
1561
1562     DMRBPixX_ = DMeanR.make<TH1D>("DMRBPix-X", "DMR of BPix-X;mean of X-residuals;modules", 100., -200, 200);
1563     DMRBPixY_ = DMeanR.make<TH1D>("DMRBPix-Y", "DMR of BPix-Y;mean of Y-residuals;modules", 100., -200, 200);
1564
1565     DMRFPixX_ = DMeanR.make<TH1D>("DMRFPix-X", "DMR of FPix-X;mean of X-residuals;modules", 100., -200, 200);
1566     DMRFPixY_ = DMeanR.make<TH1D>("DMRFPix-Y", "DMR of FPix-Y;mean of Y-residuals;modules", 100., -200, 200);
1567
1568     DMRTIB_ = DMeanR.make<TH1D>("DMRTIB", "DMR of TIB;mean of X-residuals;modules", 100., -200, 200);
1569     DMRTOB_ = DMeanR.make<TH1D>("DMRTOB", "DMR of TOB;mean of X-residuals;modules", 100., -200, 200);
1570
1571     DMRTID_ = DMeanR.make<TH1D>("DMRTID", "DMR of TID;mean of X-residuals;modules", 100., -200, 200);
1572     DMRTEC_ = DMeanR.make<TH1D>("DMRTEC", "DMR of TEC;mean of X-residuals;modules", 100., -200, 200);
1573
1574     TFileDirectory DMeanRSplit = fs->mkdir("SplitDMRs");
1575
1576     DMRBPixXSplit_ = bookSplitDMRHistograms(DMeanRSplit, "BPix", "X", true);
1577     DMRBPixYSplit_ = bookSplitDMRHistograms(DMeanRSplit, "BPix", "Y", true);
1578
1579     DMRFPixXSplit_ = bookSplitDMRHistograms(DMeanRSplit, "FPix", "X", false);
1580     DMRFPixYSplit_ = bookSplitDMRHistograms(DMeanRSplit, "FPix", "Y", false);
1581
1582     DMRTIBSplit_ = bookSplitDMRHistograms(DMeanRSplit, "TIB", "X", true);
1583     DMRTOBSplit_ = bookSplitDMRHistograms(DMeanRSplit, "TOB", "X", true);
1584
1585     // DRnRs
1586     TFileDirectory DRnRs = fs->mkdir("DRnRs");
1587
1588     DRnRBPixX_ = DRnRs.make<TH1D>("DRnRBPix-X", "DRnR of BPix-X;rms of normalized X-residuals;modules", 100., 0., 3.);
1589     DRnRBPixY_ = DRnRs.make<TH1D>("DRnRBPix-Y", "DRnR of BPix-Y;rms of normalized Y-residuals;modules", 100., 0., 3.);
1590
1591     DRnRFPixX_ = DRnRs.make<TH1D>("DRnRFPix-X", "DRnR of FPix-X;rms of normalized X-residuals;modules", 100., 0., 3.);
1592     DRnRFPixY_ = DRnRs.make<TH1D>("DRnRFPix-Y", "DRnR of FPix-Y;rms of normalized Y-residuals;modules", 100., 0., 3.);
1593
1594     DRnRTIB_ = DRnRs.make<TH1D>("DRnRTIB", "DRnR of TIB;rms of normalized X-residuals;modules", 100., 0., 3.);
1595     DRnRTOB_ = DRnRs.make<TH1D>("DRnRTOB", "DRnR of TOB;rms of normalized Y-residuals;modules", 100., 0., 3.);
1596
1597     DRnRTID_ = DRnRs.make<TH1D>("DRnRTID", "DRnR of TID;rms of normalized X-residuals;modules", 100., 0., 3.);
1598     DRnRTEC_ = DRnRs.make<TH1D>("DRnRTEC", "DRnR of TEC;rms of normalized Y-residuals;modules", 100., 0., 3.);
1599
1600     // initialize the topology first
1601
1602     SiPixelPI::PhaseInfo ph_info(phase_);
1603     const char *path_toTopologyXML = ph_info.pathToTopoXML();
1604     const TrackerTopology standaloneTopo =
1605         StandaloneTrackerTopology::fromTrackerParametersXMLFile(edm::FileInPath(path_toTopologyXML).fullPath());
1606
1607     // initialize PixelRegionContainers
1608     PixelDMRS_x_ByLayer = std::make_unique<PixelRegions::PixelRegionContainers>(&standaloneTopo, ph_info.phase());
1609     PixelDMRS_y_ByLayer = std::make_unique<PixelRegions::PixelRegionContainers>(&standaloneTopo, ph_info.phase());
1610
1611     // book PixelRegionContainers
1612     PixelDMRS_x_ByLayer->bookAll("Barrel Pixel DMRs", "median(x'_{pred}-x'_{hit}) [#mum]", "# modules", 100, -50, 50);
1613     PixelDMRS_y_ByLayer->bookAll("Barrel Pixel DMRs", "median(y'_{pred}-y'_{hit}) [#mum]", "# modules", 100, -50, 50);
1614
1615     if (DEBUG) {
1616       auto dets = PixelRegions::attachedDets(PixelRegions::PixelId::L1, &standaloneTopo, phase_);
1617       for (const auto &det : dets) {
1618         auto myLocalTopo = PixelDMRS_x_ByLayer->getTheTopo();
1619         edm::LogVerbatim("DMRChecker") << myLocalTopo->print(det) << std::endl;
1620       }
1621     }
1622
1623     // pixel
1624
1625     for (auto &bpixid : resDetailsBPixX_) {
1626       DMRBPixX_->Fill(bpixid.second.runningMeanOfRes_);
1627       if (phase_ == SiPixelPI::phase::one) {
1628         pixelmap->fillBarrelBin("DMRsX", bpixid.first, bpixid.second.runningMeanOfRes_);
1629         fullPixelmapXDMR->fillTrackerMap(bpixid.first, bpixid.second.runningMeanOfRes_);
1630       }
1631
1632       if (DEBUG) {
1633         auto myLocalTopo = PixelDMRS_x_ByLayer->getTheTopo();
1634         edm::LogPrint("DMRChecker") << myLocalTopo->print(bpixid.first) << std::endl;
1635       }
1636
1637       PixelDMRS_x_ByLayer->fill(bpixid.first, bpixid.second.runningMeanOfRes_);
1638
1639       // split DMR
1640       if (bpixid.second.rDirection > 0) {
1641         DMRBPixXSplit_[0]->Fill(bpixid.second.runningMeanOfRes_);
1642       } else {
1643         DMRBPixXSplit_[1]->Fill(bpixid.second.runningMeanOfRes_);
1644       }
1645
1646       if (bpixid.second.hitCount < 2)
1647         DRnRBPixX_->Fill(-1);
1648       else
1649         DRnRBPixX_->Fill(sqrt(bpixid.second.runningNormVarOfRes_ / (bpixid.second.hitCount - 1)));
1650     }
1651
1652     for (auto &bpixid : resDetailsBPixY_) {
1653       DMRBPixY_->Fill(bpixid.second.runningMeanOfRes_);
1654       if (phase_ == SiPixelPI::phase::one) {
1655         pixelmap->fillBarrelBin("DMRsY", bpixid.first, bpixid.second.runningMeanOfRes_);
1656         fullPixelmapYDMR->fillTrackerMap(bpixid.first, bpixid.second.runningMeanOfRes_);
1657       }
1658
1659       PixelDMRS_y_ByLayer->fill(bpixid.first, bpixid.second.runningMeanOfRes_);
1660
1661       // split DMR
1662       if (bpixid.second.rDirection > 0) {
1663         DMRBPixYSplit_[0]->Fill(bpixid.second.runningMeanOfRes_);
1664       } else {
1665         DMRBPixYSplit_[1]->Fill(bpixid.second.runningMeanOfRes_);
1666       }
1667
1668       if (bpixid.second.hitCount < 2)
1669         DRnRBPixY_->Fill(-1);
1670       else
1671         DRnRBPixY_->Fill(sqrt(bpixid.second.runningNormVarOfRes_ / (bpixid.second.hitCount - 1)));
1672     }
1673
1674     for (auto &fpixid : resDetailsFPixX_) {
1675       DMRFPixX_->Fill(fpixid.second.runningMeanOfRes_);
1676       if (phase_ == SiPixelPI::phase::one) {
1677         pixelmap->fillForwardBin("DMRsX", fpixid.first, fpixid.second.runningMeanOfRes_);
1678         fullPixelmapXDMR->fillTrackerMap(fpixid.first, fpixid.second.runningMeanOfRes_);
1679       }
1680       PixelDMRS_x_ByLayer->fill(fpixid.first, fpixid.second.runningMeanOfRes_);
1681
1682       // split DMR
1683       if (fpixid.second.zDirection > 0) {
1684         DMRFPixXSplit_[0]->Fill(fpixid.second.runningMeanOfRes_);
1685       } else {
1686         DMRFPixXSplit_[1]->Fill(fpixid.second.runningMeanOfRes_);
1687       }
1688
1689       if (fpixid.second.hitCount < 2)
1690         DRnRFPixX_->Fill(-1);
1691       else
1692         DRnRFPixX_->Fill(sqrt(fpixid.second.runningNormVarOfRes_ / (fpixid.second.hitCount - 1)));
1693     }
1694
1695     for (auto &fpixid : resDetailsFPixY_) {
1696       DMRFPixY_->Fill(fpixid.second.runningMeanOfRes_);
1697       if (phase_ == SiPixelPI::phase::one) {
1698         pixelmap->fillForwardBin("DMRsY", fpixid.first, fpixid.second.runningMeanOfRes_);
1699         fullPixelmapXDMR->fillTrackerMap(fpixid.first, fpixid.second.runningMeanOfRes_);
1700       }
1701       PixelDMRS_y_ByLayer->fill(fpixid.first, fpixid.second.runningMeanOfRes_);
1702
1703       // split DMR
1704       if (fpixid.second.zDirection > 0) {
1705         DMRFPixYSplit_[0]->Fill(fpixid.second.runningMeanOfRes_);
1706       } else {
1707         DMRFPixYSplit_[1]->Fill(fpixid.second.runningMeanOfRes_);
1708       }
1709
1710       if (fpixid.second.hitCount < 2)
1711         DRnRFPixY_->Fill(-1);
1712       else
1713         DRnRFPixY_->Fill(sqrt(fpixid.second.runningNormVarOfRes_ / (fpixid.second.hitCount - 1)));
1714     }
1715
1716     // strips
1717
1718     for (auto &tibid : resDetailsTIB_) {
1719       DMRTIB_->Fill(tibid.second.runningMeanOfRes_);
1720
1721       // split DMR
1722       if (tibid.second.rDirection > 0) {
1723         DMRTIBSplit_[0]->Fill(tibid.second.runningMeanOfRes_);
1724       } else {
1725         DMRTIBSplit_[1]->Fill(tibid.second.runningMeanOfRes_);
1726       }
1727
1728       if (tibid.second.hitCount < 2)
1729         DRnRTIB_->Fill(-1);
1730       else
1731         DRnRTIB_->Fill(sqrt(tibid.second.runningNormVarOfRes_ / (tibid.second.hitCount - 1)));
1732     }
1733
1734     for (auto &tobid : resDetailsTOB_) {
1735       DMRTOB_->Fill(tobid.second.runningMeanOfRes_);
1736
1737       // split DMR
1738       if (tobid.second.rDirection > 0) {
1739         DMRTOBSplit_[0]->Fill(tobid.second.runningMeanOfRes_);
1740       } else {
1741         DMRTOBSplit_[1]->Fill(tobid.second.runningMeanOfRes_);
1742       }
1743
1744       if (tobid.second.hitCount < 2)
1745         DRnRTOB_->Fill(-1);
1746       else
1747         DRnRTOB_->Fill(sqrt(tobid.second.runningNormVarOfRes_ / (tobid.second.hitCount - 1)));
1748     }
1749
1750     for (auto &tidid : resDetailsTID_) {
1751       DMRTID_->Fill(tidid.second.runningMeanOfRes_);
1752
1753       if (tidid.second.hitCount < 2)
1754         DRnRTID_->Fill(-1);
1755       else
1756         DRnRTID_->Fill(sqrt(tidid.second.runningNormVarOfRes_ / (tidid.second.hitCount - 1)));
1757     }
1758
1759     for (auto &tecid : resDetailsTEC_) {
1760       DMRTEC_->Fill(tecid.second.runningMeanOfRes_);
1761
1762       if (tecid.second.hitCount < 2)
1763         DRnRTEC_->Fill(-1);
1764       else
1765         DRnRTEC_->Fill(sqrt(tecid.second.runningNormVarOfRes_ / (tecid.second.hitCount - 1)));
1766     }
1767
1768     edm::LogPrint("DMRChecker") << "n. of bpix modules " << resDetailsBPixX_.size() << std::endl;
1769     edm::LogPrint("DMRChecker") << "n. of fpix modules " << resDetailsFPixX_.size() << std::endl;
1770
1771     if (phase_ == SiPixelPI::phase::zero) {
1772       pmap->save(true, 0, 0, "PixelHitMap.pdf", 600, 800);
1773       pmap->save(true, 0, 0, "PixelHitMap.png", 500, 750);
1774     }
1775
1776     tmap->save(true, 0, 0, "StripHitMap.pdf");
1777     tmap->save(true, 0, 0, "StripHitMap.png");
1778
1779     if (phase_ == SiPixelPI::phase::one) {
1780       gStyle->SetPalette(kRainBow);
1781       pixelmap->beautifyAllHistograms();
1782
1783       TCanvas cBX("CanvXBarrel", "CanvXBarrel", 1200, 1000);
1784       pixelmap->drawBarrelMaps("DMRsX", cBX);
1785       cBX.SaveAs("pixelBarrelDMR_x.png");
1786
1787       TCanvas cFX("CanvXForward", "CanvXForward", 1600, 1000);
1788       pixelmap->drawForwardMaps("DMRsX", cFX);
1789       cFX.SaveAs("pixelForwardDMR_x.png");
1790
1791       TCanvas cBY("CanvYBarrel", "CanvYBarrel", 1200, 1000);
1792       pixelmap->drawBarrelMaps("DMRsY", cBY);
1793       cBY.SaveAs("pixelBarrelDMR_y.png");
1794
1795       TCanvas cFY("CanvXForward", "CanvXForward", 1600, 1000);
1796       pixelmap->drawForwardMaps("DMRsY", cFY);
1797       cFY.SaveAs("pixelForwardDMR_y.png");
1798
1799       TCanvas cFullPixelxDMR("CanvFullPixelX", "CanvFullPixelX", 3000, 2000);
1800       fullPixelmapXDMR->printTrackerMap(cFullPixelxDMR);
1801       cFullPixelxDMR.SaveAs("fullPixelDMR_x.png");
1802
1803       TCanvas cFullPixelyDMR("CanvFullPixelX", "CanvFullPixelY", 3000, 2000);
1804       fullPixelmapXDMR->printTrackerMap(cFullPixelyDMR);
1805       cFullPixelyDMR.SaveAs("fullPixelDMR_y.png");
1806     }
1807
1808     // take care now of the 1D histograms
1809     gStyle->SetOptStat("emr");
1810     PixelDMRS_x_ByLayer->beautify(2, 0);
1811     PixelDMRS_y_ByLayer->beautify(2, 0);
1812
1813     TCanvas DMRxBarrel("DMRxBarrelCanv", "x-coordinate", 1400, 1200);
1814     DMRxBarrel.Divide(2, 2);
1815     PixelDMRS_x_ByLayer->draw(DMRxBarrel, true, "HISTS");
1816     adjustCanvases(DMRxBarrel, true);
1817     for (unsigned int c = 1; c <= 4; c++) {
1818       DMRxBarrel.cd(c)->Update();
1819     }
1820     PixelDMRS_x_ByLayer->stats();
1821
1822     TCanvas DMRxForward("DMRxForwardCanv", "x-coordinate", 1400, 1200);
1823     DMRxForward.Divide(4, 3);
1824     PixelDMRS_x_ByLayer->draw(DMRxForward, false, "HISTS");
1825     adjustCanvases(DMRxForward, false);
1826     for (unsigned int c = 1; c <= 12; c++) {
1827       DMRxForward.cd(c)->Update();
1828     }
1829     PixelDMRS_x_ByLayer->stats();
1830
1831     DMRxBarrel.SaveAs("DMR_x_Barrel_ByLayer.png");
1832     DMRxForward.SaveAs("DMR_x_Forward_ByRing.png");
1833
1834     TCanvas DMRyBarrel("DMRyBarrelCanv", "y-coordinate", 1400, 1200);
1835     DMRyBarrel.Divide(2, 2);
1836     PixelDMRS_y_ByLayer->draw(DMRyBarrel, true, "HISTS");
1837     adjustCanvases(DMRyBarrel, true);
1838     for (unsigned int c = 1; c <= 4; c++) {
1839       DMRyBarrel.cd(c)->Update();
1840     }
1841     PixelDMRS_y_ByLayer->stats();
1842
1843     TCanvas DMRyForward("DMRyForwardCanv", "y-coordinate", 1400, 1200);
1844     DMRyForward.Divide(4, 3);
1845     PixelDMRS_y_ByLayer->draw(DMRyForward, false, "HISTS");
1846     adjustCanvases(DMRyForward, false);
1847     for (unsigned int c = 1; c <= 12; c++) {
1848       DMRyForward.cd(c)->Update();
1849     }
1850     PixelDMRS_y_ByLayer->stats();
1851
1852     DMRyBarrel.SaveAs("DMR_y_Barrel_ByLayer.png");
1853     DMRyForward.SaveAs("DMR_y_Forward_ByRing.png");
1854   }
1855
1856   //*************************************************************
1857   // Adjust canvas for DMRs
1858   //*************************************************************
1859   void adjustCanvases(TCanvas &canvas, bool isBarrel) {
1860     unsigned int maxPads = isBarrel ? 4 : 12;
1861     for (unsigned int c = 1; c <= maxPads; c++) {
1862       canvas.cd(c);
1863       SiPixelPI::adjustCanvasMargins(canvas.cd(c), 0.06, 0.12, 0.12, 0.05);
1864     }
1865
1866     auto ltx = TLatex();
1867     ltx.SetTextFont(62);
1868     ltx.SetTextSize(0.05);
1869     ltx.SetTextAlign(11);
1870
1871     std::string toAppend = canvas.GetTitle();
1872
1873     for (unsigned int c = 1; c <= maxPads; c++) {
1874       auto index = isBarrel ? c - 1 : c + 3;
1875       canvas.cd(c);
1876       ltx.DrawLatexNDC(gPad->GetLeftMargin(),
1877                        1 - gPad->GetTopMargin() + 0.01,
1878                        (PixelRegions::IDlabels.at(index) + " " + toAppend).c_str());
1879     }
1880   }
1881
1882   //*************************************************************
1883   // check if the hit is 2D
1884   //*************************************************************
1885   bool isHit2D(const TrackingRecHit &hit) {
1886     bool countStereoHitAs2D_ = true;
1887     // we count SiStrip stereo modules as 2D if selected via countStereoHitAs2D_
1888     // (since they provide theta information)
1889     if (!hit.isValid() ||
1890         (hit.dimension() < 2 && !countStereoHitAs2D_ && !dynamic_cast<const SiStripRecHit1D *>(&hit))) {
1891       return false;  // real RecHit1D - but SiStripRecHit1D depends on countStereoHitAs2D_
1892     } else {
1893       const DetId detId(hit.geographicalId());
1894       if (detId.det() == DetId::Tracker) {
1895         if (detId.subdetId() == kBPIX || detId.subdetId() == kFPIX) {
1896           return true;  // pixel is always 2D
1897         } else {        // should be SiStrip now
1898           const SiStripDetId stripId(detId);
1899           if (stripId.stereo())
1900             return countStereoHitAs2D_;  // stereo modules
1901           else if (dynamic_cast<const SiStripRecHit1D *>(&hit) || dynamic_cast<const SiStripRecHit2D *>(&hit))
1902             return false;  // rphi modules hit
1903           //the following two are not used any more since ages...
1904           else if (dynamic_cast<const SiStripMatchedRecHit2D *>(&hit))
1905             return true;  // matched is 2D
1906           else if (dynamic_cast<const ProjectedSiStripRecHit2D *>(&hit)) {
1907             const ProjectedSiStripRecHit2D *pH = static_cast<const ProjectedSiStripRecHit2D *>(&hit);
1908             return (countStereoHitAs2D_ && this->isHit2D(pH->originalHit()));  // depends on original...
1909           } else {
1910             edm::LogError("UnknownType") << "@SUB=DMRChecker::isHit2D"
1911                                          << "Tracker hit not in pixel, neither SiStripRecHit[12]D nor "
1912                                          << "SiStripMatchedRecHit2D nor ProjectedSiStripRecHit2D.";
1913             return false;
1914           }
1915         }
1916       } else {  // not tracker??
1917         edm::LogWarning("DetectorMismatch") << "@SUB=DMRChecker::isHit2D"
1918                                             << "Hit not in tracker with 'official' dimension >=2.";
1919         return true;  // dimension() >= 2 so accept that...
1920       }
1921     }
1922     // never reached...
1923   }
1924
1925   //*************************************************************
1926   // Generic booker of split DMRs
1927   //*************************************************************
1928   std::array<TH1D *, 2> bookSplitDMRHistograms(TFileDirectory dir,
1929                                                std::string subdet,
1930                                                std::string vartype,
1931                                                bool isBarrel) {
1932     TH1F::SetDefaultSumw2(kTRUE);
1933
1934     std::array<TH1D *, 2> out;
1935     std::array<std::string, 2> sign_name = {{"plus", "minus"}};
1936     std::array<std::string, 2> sign = {{">0", "<0"}};
1937     for (unsigned int i = 0; i < 2; i++) {
1938       const char *name_;
1939       const char *title_;
1940       const char *axisTitle_;
1941
1942       if (isBarrel) {
1943         name_ = Form("DMR%s_%s_rDir%s", subdet.c_str(), vartype.c_str(), sign_name[i].c_str());
1944         title_ = Form("Split DMR of %s-%s (rDir%s)", subdet.c_str(), vartype.c_str(), sign[i].c_str());
1945         axisTitle_ = Form("mean of %s-residuals (rDir%s);modules", vartype.c_str(), sign[i].c_str());
1946       } else {
1947         name_ = Form("DMR%s_%s_zDir%s", subdet.c_str(), vartype.c_str(), sign_name[i].c_str());
1948         title_ = Form("Split DMR of %s-%s (zDir%s)", subdet.c_str(), vartype.c_str(), sign[i].c_str());
1949         axisTitle_ = Form("mean of %s-residuals (zDir%s);modules", vartype.c_str(), sign[i].c_str());
1950       }
1951
1952       out[i] = dir.make<TH1D>(name_, fmt::sprintf("%s;%s", title_, axisTitle_).c_str(), 100., -200, 200);
1953     }
1954     return out;
1955   }
1956
1957   //*************************************************************
1958   // Generic booker function
1959   //*************************************************************
1960   std::map<unsigned int, TH1D *> bookResidualsHistogram(
1961       TFileDirectory dir, unsigned int theNLayers, std::string resType, std::string varType, std::string detType) {
1962     TH1F::SetDefaultSumw2(kTRUE);
1963
1964     std::pair<double, double> limits;
1965
1966     if (varType.find("Res") != std::string::npos) {
1967       limits = std::make_pair(-1000., 1000);
1968     } else {
1969       limits = std::make_pair(-3., 3.);
1970     }
1971
1972     std::map<unsigned int, TH1D *> h;
1973
1974     for (unsigned int i = 1; i <= theNLayers; i++) {
1975       const char *name_;
1976       const char *title_;
1977       std::string xAxisTitle_;
1978
1979       if (varType.find("Res") != std::string::npos) {
1980         xAxisTitle_ = fmt::sprintf("res_{%s'} [#mum]", resType);
1981       } else {
1982         xAxisTitle_ = fmt::sprintf("res_{%s'}/#sigma_{res_{%s`}}", resType, resType);
1983       }
1984
1985       unsigned int side = -1;
1986       if (detType.find("FPix") != std::string::npos) {
1987         side = (i - 1) / 3 + 1;
1988         unsigned int plane = (i - 1) % 3 + 1;
1989
1990         std::string theSide = "";
1991         if (side == 1) {
1992           theSide = "Z-";
1993         } else {
1994           theSide = "Z+";
1995         }
1996
1997         name_ = Form("h_%s%s%s_side%i_disk%i", detType.c_str(), varType.c_str(), resType.c_str(), side, plane);
1998         title_ = Form("%s (%s, disk %i) track %s-%s;%s;hits",
1999                       detType.c_str(),
2000                       theSide.c_str(),
2001                       plane,
2002                       resType.c_str(),
2003                       varType.c_str(),
2004                       xAxisTitle_.c_str());
2005
2006       } else {
2007         name_ = Form("h_%s%s%s_layer%i", detType.c_str(), varType.c_str(), resType.c_str(), i);
2008         title_ = Form("%s (layer %i) track %s-%s;%s;hits",
2009                       detType.c_str(),
2010                       i,
2011                       resType.c_str(),
2012                       varType.c_str(),
2013                       xAxisTitle_.c_str());
2014
2015         //edm::LogPrint("DMRChecker")<<"bookResidualsHistogram(): "<<i<<" layer:"<<i<<std::endl;
2016       }
2017
2018       h[i] = dir.make<TH1D>(name_, title_, 100, limits.first, limits.second);
2019     }
2020
2021     return h;
2022   }
2023
2024   //*************************************************************
2025   // Generic filler function
2026   //*************************************************************
2027   void fillByIndex(std::map<unsigned int, TH1D *> &h, unsigned int index, double x) {
2028     if (h.count(index) != 0) {
2029       //if(TString(h[index]->GetName()).Contains("BPix"))
2030       //edm::LogPrint("DMRChecker")<<"fillByIndex() index: "<< index << " filling histogram: "<< h[index]->GetName() << std::endl;
2031
2032       double min = h[index]->GetXaxis()->GetXmin();
2033       double max = h[index]->GetXaxis()->GetXmax();
2034       if (x < min)
2035         h[index]->Fill(min);
2036       else if (x >= max)
2037         h[index]->Fill(0.99999 * max);
2038       else
2039         h[index]->Fill(x);
2040     }
2041   }
2042
2043   //*************************************************************
2044   // Implementation of the online variance algorithm
2045   // as in https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Online_algorithm
2046   //*************************************************************
2047   void updateOnlineMomenta(running::estimatorMap &myDetails, uint32_t theID, float the_data, float the_pull) {
2048     myDetails[theID].hitCount += 1;
2049
2050     float delta = 0;
2051     float n_delta = 0;
2052
2053     if (myDetails[theID].hitCount != 1) {
2054       delta = the_data - myDetails[theID].runningMeanOfRes_;
2055       n_delta = the_pull - myDetails[theID].runningNormMeanOfRes_;
2056       myDetails[theID].runningMeanOfRes_ += (delta / myDetails[theID].hitCount);
2057       myDetails[theID].runningNormMeanOfRes_ += (n_delta / myDetails[theID].hitCount);
2058     } else {
2059       myDetails[theID].runningMeanOfRes_ = the_data;
2060       myDetails[theID].runningNormMeanOfRes_ = the_pull;
2061     }
2062
2063     float delta2 = the_data - myDetails[theID].runningMeanOfRes_;
2064     float n_delta2 = the_pull - myDetails[theID].runningNormMeanOfRes_;
2065
2066     myDetails[theID].runningVarOfRes_ += delta * delta2;
2067     myDetails[theID].runningNormVarOfRes_ += n_delta * n_delta2;
2068   }
2069
2070   //*************************************************************
2071   // Fill the histograms using the running::estimatorMap
2072   //**************************************************************
2073   void fillDMRs(const running::estimatorMap &myDetails,
2074                 TH1D *DMR,
2075                 TH1D *DRnR,
2076                 std::array<TH1D *, 2> DMRSplit,
2077                 std::unique_ptr<PixelRegions::PixelRegionContainers> regionalDMR) {
2078     for (const auto &element : myDetails) {
2079       // DMR
2080       DMR->Fill(element.second.runningMeanOfRes_);
2081
2082       // DMR by layer
2083       if (regionalDMR.get()) {
2084         regionalDMR->fill(element.first, element.second.runningMeanOfRes_);
2085       }
2086
2087       // split DMR
2088       if (element.second.rOrZDirection > 0) {
2089         DMRSplit[0]->Fill(element.second.runningMeanOfRes_);
2090       } else {
2091         DMRSplit[1]->Fill(element.second.runningMeanOfRes_);
2092       }
2093
2094       // DRnR
2095       if (element.second.hitCount < 2) {
2096         DRnR->Fill(-1);
2097       } else {
2098         DRnR->Fill(sqrt(element.second.runningNormVarOfRes_ / (element.second.hitCount - 1)));
2099       }
2100     }
2101   }
2102 };
2103
2104 //*************************************************************
2105 void DMRChecker::fillDescriptions(edm::ConfigurationDescriptions &descriptions)
2106 //*************************************************************
2107 {
2108   edm::ParameterSetDescription desc;
2109   desc.setComment("Generic track analyzer to check ALCARECO sample quantities / compute fast DMRs");
2110   desc.add<edm::InputTag>("TkTag", edm::InputTag("generalTracks"));
2111   desc.add<edm::InputTag>("TriggerResultsTag", edm::InputTag("TriggerResults", "", "HLT"));
2112   desc.add<edm::InputTag>("BeamSpotTag", edm::InputTag("offlineBeamSpot"));
2113   desc.add<edm::InputTag>("VerticesTag", edm::InputTag("offlinePrimaryVertices"));
2114   desc.add<bool>("isCosmics", false);
2115   desc.add<bool>("doLatencyAnalysis", true);
2116   descriptions.addWithDefaultLabel(desc);
2117 }
2118
2119 DEFINE_FWK_MODULE(DMRChecker);