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