StandaloneTrackMonitor.cc

Go to the documentation of this file.

#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/EDAnalyzer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Utilities/interface/InputTag.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "DQMServices/Core/interface/DQMStore.h"
#include "DQMServices/Core/interface/DQMEDAnalyzer.h"
#include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "DataFormats/VertexReco/interface/VertexFwd.h"
#include "DataFormats/BeamSpot/interface/BeamSpot.h"
#include "DataFormats/Common/interface/DetSetVectorNew.h"
#include "RecoLocalTracker/SiStripClusterizer/interface/SiStripClusterInfo.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/BeamSpot/interface/BeamSpot.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "DataFormats/VertexReco/interface/VertexFwd.h"
#include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"
#include "Geometry/CommonDetUnit/interface/GluedGeomDet.h"
#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
#include "DataFormats/TrackerRecHit2D/interface/BaseTrackerRecHit.h"
#include "DataFormats/TrackerRecHit2D/interface/OmniClusterRef.h"
#include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
#include "DataFormats/DetId/interface/DetId.h"
#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"

#include "TFile.h"
#include "TH1.h"
#include "TMath.h"
#include "TPRegexp.h"

#include <string>
#include <vector>
#include <map>
#include <set>

class StandaloneTrackMonitor : public DQMEDAnalyzer {
public:
  StandaloneTrackMonitor(const edm::ParameterSet&);

protected:
  void analyze(edm::Event const& iEvent, edm::EventSetup const& iSetup) override;
  void processHit(const TrackingRecHit& recHit, edm::EventSetup const& iSetup, double wfac = 1);
  void processClusters(edm::Event const& iEvent, edm::EventSetup const& iSetup, double wfac = 1);
  void addClusterToMap(uint32_t detid, const SiStripCluster* cluster);
  void bookHistograms(DQMStore::IBooker&, edm::Run const&, edm::EventSetup const&) override;
  void dqmBeginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) override;

private:
  edm::ParameterSet parameters_;

  std::string moduleName_;
  std::string folderName_;
  const edm::InputTag trackTag_;
  const edm::InputTag bsTag_;
  const edm::InputTag vertexTag_;
  const edm::InputTag puSummaryTag_;
  const edm::InputTag clusterTag_;
  const edm::EDGetTokenT<reco::TrackCollection> trackToken_;
  const edm::EDGetTokenT<reco::BeamSpot> bsToken_;
  const edm::EDGetTokenT<reco::VertexCollection> vertexToken_;
  const edm::EDGetTokenT<std::vector<PileupSummaryInfo> > puSummaryToken_;
  const edm::EDGetTokenT<edmNew::DetSetVector<SiStripCluster> > clusterToken_;

  const std::string trackQuality_;
  SiStripClusterInfo siStripClusterInfo_;
  const bool doPUCorrection_;
  const bool isMC_;
  const bool haveAllHistograms_;
  const std::string puScaleFactorFile_;
  const bool verbose_;

  MonitorElement* trackEtaH_;
  MonitorElement* trackEtaerrH_;
  MonitorElement* trackCosThetaH_;
  MonitorElement* trackThetaerrH_;
  MonitorElement* trackPhiH_;
  MonitorElement* trackPhierrH_;
  MonitorElement* trackPH_;
  MonitorElement* trackPtH_;
  MonitorElement* trackPtUpto2GeVH_;
  MonitorElement* trackPtOver10GeVH_;
  MonitorElement* trackPterrH_;
  MonitorElement* trackqOverpH_;
  MonitorElement* trackqOverperrH_;
  MonitorElement* trackChargeH_;
  MonitorElement* trackChi2H_;
  MonitorElement* tracknDOFH_;
  MonitorElement* trackd0H_;
  MonitorElement* trackChi2bynDOFH_;

  MonitorElement* nlostHitsH_;
  MonitorElement* nvalidTrackerHitsH_;
  MonitorElement* nvalidPixelHitsH_;
  MonitorElement* nvalidStripHitsH_;
  MonitorElement* trkLayerwithMeasurementH_;
  MonitorElement* pixelLayerwithMeasurementH_;
  MonitorElement* stripLayerwithMeasurementH_;

  MonitorElement* beamSpotXYposH_;
  MonitorElement* beamSpotXYposerrH_;
  MonitorElement* beamSpotZposH_;
  MonitorElement* beamSpotZposerrH_;

  MonitorElement* vertexXposH_;
  MonitorElement* vertexYposH_;
  MonitorElement* vertexZposH_;
  MonitorElement* nVertexH_;

  MonitorElement* nPixBarrelH_;
  MonitorElement* nPixEndcapH_;
  MonitorElement* nStripTIBH_;
  MonitorElement* nStripTOBH_;
  MonitorElement* nStripTECH_;
  MonitorElement* nStripTIDH_;
  MonitorElement* nTracksH_;

  // MC only
  MonitorElement* bunchCrossingH_;
  MonitorElement* nPUH_;
  MonitorElement* trueNIntH_;

  // Exclusive Quantities
  MonitorElement* nHitsVspTH_;
  MonitorElement* nHitsVsnVtxH_;
  MonitorElement* nHitsVsEtaH_;
  MonitorElement* nHitsVsCosThetaH_;
  MonitorElement* nHitsVsPhiH_;
  MonitorElement* nLostHitsVspTH_;
  MonitorElement* nLostHitsVsEtaH_;
  MonitorElement* nLostHitsVsCosThetaH_;
  MonitorElement* nLostHitsVsPhiH_;

  MonitorElement* hOnTrkClusChargeThinH_;
  MonitorElement* hOnTrkClusWidthThinH_;
  MonitorElement* hOnTrkClusChargeThickH_;
  MonitorElement* hOnTrkClusWidthThickH_;

  MonitorElement* hOffTrkClusChargeThinH_;
  MonitorElement* hOffTrkClusWidthThinH_;
  MonitorElement* hOffTrkClusChargeThickH_;
  MonitorElement* hOffTrkClusWidthThickH_;

  unsigned long long m_cacheID_;

  std::vector<float> vpu_;
  std::map<uint32_t, std::set<const SiStripCluster*> > clusterMap_;
  const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geomToken_;
  const TrackerGeometry* tkGeom_ = nullptr;
};

// -----------------------------
// constructors and destructor
// -----------------------------
StandaloneTrackMonitor::StandaloneTrackMonitor(const edm::ParameterSet& ps)
    : parameters_(ps),
      moduleName_(parameters_.getUntrackedParameter<std::string>("moduleName", "StandaloneTrackMonitor")),
      folderName_(parameters_.getUntrackedParameter<std::string>("folderName", "highPurityTracks")),
      trackTag_(parameters_.getUntrackedParameter<edm::InputTag>("trackInputTag", edm::InputTag("generalTracks"))),
      bsTag_(parameters_.getUntrackedParameter<edm::InputTag>("offlineBeamSpot", edm::InputTag("offlineBeamSpot"))),
      vertexTag_(
          parameters_.getUntrackedParameter<edm::InputTag>("vertexTag", edm::InputTag("offlinePrimaryVertices"))),
      puSummaryTag_(parameters_.getUntrackedParameter<edm::InputTag>("puTag", edm::InputTag("addPileupInfo"))),
      clusterTag_(parameters_.getUntrackedParameter<edm::InputTag>("clusterTag", edm::InputTag("siStripClusters"))),
      trackToken_(consumes<reco::TrackCollection>(trackTag_)),
      bsToken_(consumes<reco::BeamSpot>(bsTag_)),
      vertexToken_(consumes<reco::VertexCollection>(vertexTag_)),
      puSummaryToken_(consumes<std::vector<PileupSummaryInfo> >(puSummaryTag_)),
      clusterToken_(consumes<edmNew::DetSetVector<SiStripCluster> >(clusterTag_)),
      trackQuality_(parameters_.getUntrackedParameter<std::string>("trackQuality", "highPurity")),
      siStripClusterInfo_(consumesCollector()),
      doPUCorrection_(parameters_.getUntrackedParameter<bool>("doPUCorrection", false)),
      isMC_(parameters_.getUntrackedParameter<bool>("isMC", false)),
      haveAllHistograms_(parameters_.getUntrackedParameter<bool>("haveAllHistograms", false)),
      puScaleFactorFile_(
          parameters_.getUntrackedParameter<std::string>("puScaleFactorFile", "PileupScaleFactor_run203002.root")),
      verbose_(parameters_.getUntrackedParameter<bool>("verbose", false)),
      geomToken_(esConsumes<TrackerGeometry, TrackerDigiGeometryRecord, edm::Transition::BeginRun>()) {
  trackEtaH_ = nullptr;
  trackEtaerrH_ = nullptr;
  trackCosThetaH_ = nullptr;
  trackThetaerrH_ = nullptr;
  trackPhiH_ = nullptr;
  trackPhierrH_ = nullptr;
  trackPH_ = nullptr;
  trackPtH_ = nullptr;
  trackPtUpto2GeVH_ = nullptr;
  trackPtOver10GeVH_ = nullptr;
  trackPterrH_ = nullptr;
  trackqOverpH_ = nullptr;
  trackqOverperrH_ = nullptr;
  trackChargeH_ = nullptr;
  nlostHitsH_ = nullptr;
  nvalidTrackerHitsH_ = nullptr;
  nvalidPixelHitsH_ = nullptr;
  nvalidStripHitsH_ = nullptr;
  trkLayerwithMeasurementH_ = nullptr;
  pixelLayerwithMeasurementH_ = nullptr;
  stripLayerwithMeasurementH_ = nullptr;
  beamSpotXYposH_ = nullptr;
  beamSpotXYposerrH_ = nullptr;
  beamSpotZposH_ = nullptr;
  beamSpotZposerrH_ = nullptr;
  trackChi2H_ = nullptr;
  tracknDOFH_ = nullptr;
  trackd0H_ = nullptr;
  trackChi2bynDOFH_ = nullptr;
  vertexXposH_ = nullptr;
  vertexYposH_ = nullptr;
  vertexZposH_ = nullptr;

  nPixBarrelH_ = nullptr;
  nPixEndcapH_ = nullptr;
  nStripTIBH_ = nullptr;
  nStripTOBH_ = nullptr;
  nStripTECH_ = nullptr;
  nStripTIDH_ = nullptr;
  nTracksH_ = nullptr;

  // for MC only
  nVertexH_ = nullptr;
  bunchCrossingH_ = nullptr;
  nPUH_ = nullptr;
  trueNIntH_ = nullptr;

  nHitsVspTH_ = nullptr;
  nHitsVsEtaH_ = nullptr;
  nHitsVsCosThetaH_ = nullptr;
  nHitsVsPhiH_ = nullptr;
  nHitsVsnVtxH_ = nullptr;
  nLostHitsVspTH_ = nullptr;
  nLostHitsVsEtaH_ = nullptr;
  nLostHitsVsCosThetaH_ = nullptr;
  nLostHitsVsPhiH_ = nullptr;

  hOnTrkClusChargeThinH_ = nullptr;
  hOnTrkClusWidthThinH_ = nullptr;
  hOnTrkClusChargeThickH_ = nullptr;
  hOnTrkClusWidthThickH_ = nullptr;

  hOffTrkClusChargeThinH_ = nullptr;
  hOffTrkClusWidthThinH_ = nullptr;
  hOffTrkClusChargeThickH_ = nullptr;
  hOffTrkClusWidthThickH_ = nullptr;

  // Read pileup weight factors
  if (isMC_ && doPUCorrection_) {
    vpu_.clear();
    TFile* f1 = TFile::Open(puScaleFactorFile_.c_str());
    TH1F* h1 = dynamic_cast<TH1F*>(f1->Get("pileupweight"));
    for (int i = 1; i <= h1->GetNbinsX(); ++i)
      vpu_.push_back(h1->GetBinContent(i));
    f1->Close();
  }
}

void StandaloneTrackMonitor::dqmBeginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
  tkGeom_ = &(iSetup.getData(geomToken_));
}

void StandaloneTrackMonitor::bookHistograms(DQMStore::IBooker& iBook,
                                            edm::Run const& iRun,
                                            edm::EventSetup const& iSetup) {
  edm::ParameterSet TrackEtaHistoPar = parameters_.getParameter<edm::ParameterSet>("trackEtaH");
  edm::ParameterSet TrackPtHistoPar = parameters_.getParameter<edm::ParameterSet>("trackPtH");

  std::string currentFolder = moduleName_ + "/" + folderName_;
  iBook.setCurrentFolder(currentFolder);

  // The following are common with the official tool
  if (haveAllHistograms_) {
    if (!trackEtaH_)
      trackEtaH_ = iBook.book1D("trackEta",
                                "Track Eta",
                                TrackEtaHistoPar.getParameter<int32_t>("Xbins"),
                                TrackEtaHistoPar.getParameter<double>("Xmin"),
                                TrackEtaHistoPar.getParameter<double>("Xmax"));
    if (!trackEtaerrH_)
      trackEtaerrH_ = iBook.book1D("trackEtaerr", "Track Eta Error", 50, 0.0, 1.0);
    if (!trackCosThetaH_)
      trackCosThetaH_ = iBook.book1D("trackCosTheta", "Track Cos(Theta)", 50, -1.0, 1.0);
    if (!trackThetaerrH_)
      trackThetaerrH_ = iBook.book1D("trackThetaerr", "Track Theta Error", 50, 0.0, 1.0);
    if (!trackPhiH_)
      trackPhiH_ = iBook.book1D("trackPhi", "Track Phi", 70, -3.5, 3.5);
    if (!trackPhierrH_)
      trackPhierrH_ = iBook.book1D("trackPhierr", "Track Phi Error", 50, 0.0, 1.0);

    if (!trackPH_)
      trackPH_ = iBook.book1D("trackP", "Track 4-momentum", 50, 0.0, 10.0);
    if (!trackPtH_)
      trackPtH_ = iBook.book1D("trackPt",
                               "Track Pt",
                               TrackPtHistoPar.getParameter<int32_t>("Xbins"),
                               TrackPtHistoPar.getParameter<double>("Xmin"),
                               TrackPtHistoPar.getParameter<double>("Xmax"));
    if (!trackPtUpto2GeVH_)
      trackPtUpto2GeVH_ = iBook.book1D("trackPtUpto2GeV", "Track Pt upto 2GeV", 100, 0, 2.0);
    if (!trackPtOver10GeVH_)
      trackPtOver10GeVH_ = iBook.book1D("trackPtOver10GeV", "Track Pt greater than 10 GeV", 100, 0, 100.0);
    if (!trackPterrH_)
      trackPterrH_ = iBook.book1D("trackPterr", "Track Pt Error", 100, 0.0, 100.0);
    if (!trackqOverpH_)
      trackqOverpH_ = iBook.book1D("trackqOverp", "q Over p", 40, -10.0, 10.0);
    if (!trackqOverperrH_)
      trackqOverperrH_ = iBook.book1D("trackqOverperr", "q Over p Error", 50, 0.0, 25.0);
    if (!trackChargeH_)
      trackChargeH_ = iBook.book1D("trackCharge", "Track Charge", 50, -5, 5);
    if (!trackChi2H_)
      trackChi2H_ = iBook.book1D("trackChi2", "Chi2", 100, 0.0, 100.0);
    if (!tracknDOFH_)
      tracknDOFH_ = iBook.book1D("tracknDOF", "nDOF", 100, 0.0, 100.0);
    if (!trackd0H_)
      trackd0H_ = iBook.book1D("trackd0", "Track d0", 100, -1, 1);
    if (!trackChi2bynDOFH_)
      trackChi2bynDOFH_ = iBook.book1D("trackChi2bynDOF", "Chi2 Over nDOF", 100, 0.0, 10.0);

    if (!nlostHitsH_)
      nlostHitsH_ = iBook.book1D("nlostHits", "No. of Lost Hits", 10, 0.0, 10.0);
    if (!nvalidTrackerHitsH_)
      nvalidTrackerHitsH_ = iBook.book1D("nvalidTrackerhits", "No. of Valid Tracker Hits", 35, 0.0, 35.0);
    if (!nvalidPixelHitsH_)
      nvalidPixelHitsH_ = iBook.book1D("nvalidPixelHits", "No. of Valid Hits in Pixel", 10, 0.0, 10.0);
    if (!nvalidStripHitsH_)
      nvalidStripHitsH_ = iBook.book1D("nvalidStripHits", "No.of Valid Hits in Strip", 25, 0.0, 25.0);

    if (!trkLayerwithMeasurementH_)
      trkLayerwithMeasurementH_ = iBook.book1D("trkLayerwithMeasurement", "No. of Layers per Track", 25, 0.0, 25.0);
    if (!pixelLayerwithMeasurementH_)
      pixelLayerwithMeasurementH_ =
          iBook.book1D("pixelLayerwithMeasurement", "No. of Pixel Layers per Track", 10, 0.0, 10.0);
    if (!stripLayerwithMeasurementH_)
      stripLayerwithMeasurementH_ =
          iBook.book1D("stripLayerwithMeasurement", "No. of Strip Layers per Track", 20, 0.0, 20.0);

    if (!beamSpotXYposH_)
      beamSpotXYposH_ = iBook.book1D("beamSpotXYpos", "XY position of beam spot", 40, -4.0, 4.0);
    if (!beamSpotXYposerrH_)
      beamSpotXYposerrH_ = iBook.book1D("beamSpotXYposerr", "Error in XY position of beam spot", 20, 0.0, 4.0);
    if (!beamSpotZposH_)
      beamSpotZposH_ = iBook.book1D("beamSpotZpos", "Z position of beam spot", 100, -20.0, 20.0);
    if (!beamSpotZposerrH_)
      beamSpotZposerrH_ = iBook.book1D("beamSpotZposerr", "Error in Z position of beam spot", 50, 0.0, 5.0);

    if (!vertexXposH_)
      vertexXposH_ = iBook.book1D("vertexXpos", "Vertex X position", 50, -1.0, 1.0);
    if (!vertexYposH_)
      vertexYposH_ = iBook.book1D("vertexYpos", "Vertex Y position", 50, -1.0, 1.0);
    if (!vertexZposH_)
      vertexZposH_ = iBook.book1D("vertexZpos", "Vertex Z position", 100, -20.0, 20.0);
    if (!nVertexH_)
      nVertexH_ = iBook.book1D("nVertex", "# of vertices", 60, -0.5, 59.5);

    if (!nPixBarrelH_)
      nPixBarrelH_ = iBook.book1D("nHitPixelBarrel", "No. of hits in Pixel Barrel per Track", 20, 0, 20.0);
    if (!nPixEndcapH_)
      nPixEndcapH_ = iBook.book1D("nHitPixelEndcap", "No. of hits in Pixel Endcap per Track", 20, 0, 20.0);
    if (!nStripTIBH_)
      nStripTIBH_ = iBook.book1D("nHitStripTIB", "No. of hits in Strip TIB per Track", 30, 0, 30.0);
    if (!nStripTOBH_)
      nStripTOBH_ = iBook.book1D("nHitStripTOB", "No. of hits in Strip TOB per Track", 30, 0, 30.0);
    if (!nStripTECH_)
      nStripTECH_ = iBook.book1D("nHitStripTEC", "No. of hits in Strip TEC per Track", 30, 0, 30.0);
    if (!nStripTIDH_)
      nStripTIDH_ = iBook.book1D("nHitStripTID", "No. of hits in Strip TID per Tracks", 30, 0, 30.0);

    if (!nTracksH_)
      nTracksH_ = iBook.book1D("nTracks", "No. of Tracks", 100, -0.5, 999.5);
  }
  if (isMC_) {
    if (!bunchCrossingH_)
      bunchCrossingH_ = iBook.book1D("bunchCrossing", "Bunch Crosssing", 60, 0, 60.0);
    if (!nPUH_)
      nPUH_ = iBook.book1D("nPU", "No of Pileup", 60, 0, 60.0);
    if (!trueNIntH_)
      trueNIntH_ = iBook.book1D("trueNInt", "True no of Interactions", 60, 0, 60.0);
  }
  // Exclusive histograms
  if (!nHitsVspTH_)
    nHitsVspTH_ = iBook.bookProfile("nHitsVspT",
                                    "Number of Hits Vs pT",
                                    TrackPtHistoPar.getParameter<int32_t>("Xbins"),
                                    TrackPtHistoPar.getParameter<double>("Xmin"),
                                    TrackPtHistoPar.getParameter<double>("Xmax"),
                                    0.0,
                                    0.0,
                                    "g");
  if (!nHitsVsnVtxH_)
    nHitsVsnVtxH_ = iBook.bookProfile("nHitsVsnVtx", "Number of Hits Vs Number of Vertex", 100, 0.0, 50, 0.0, 0.0, "g");
  if (!nHitsVsEtaH_)
    nHitsVsEtaH_ = iBook.bookProfile("nHitsVsEta",
                                     "Number of Hits Vs Eta",
                                     TrackEtaHistoPar.getParameter<int32_t>("Xbins"),
                                     TrackEtaHistoPar.getParameter<double>("Xmin"),
                                     TrackEtaHistoPar.getParameter<double>("Xmax"),
                                     0.0,
                                     0.0,
                                     "g");
  if (!nHitsVsCosThetaH_)
    nHitsVsCosThetaH_ =
        iBook.bookProfile("nHitsVsCosTheta", "Number of Hits Vs Cos(Theta)", 50, -1.0, 1.0, 0.0, 0.0, "g");
  if (!nHitsVsPhiH_)
    nHitsVsPhiH_ = iBook.bookProfile("nHitsVsPhi", "Number of Hits Vs Phi", 100, -3.5, 3.5, 0.0, 0.0, "g");

  if (!nLostHitsVspTH_)
    nLostHitsVspTH_ = iBook.bookProfile("nLostHitsVspT",
                                        "Number of Lost Hits Vs pT",
                                        TrackPtHistoPar.getParameter<int32_t>("Xbins"),
                                        TrackPtHistoPar.getParameter<double>("Xmin"),
                                        TrackPtHistoPar.getParameter<double>("Xmax"),
                                        0.0,
                                        0.0,
                                        "g");
  if (!nLostHitsVsEtaH_)
    nLostHitsVsEtaH_ = iBook.bookProfile("nLostHitsVsEta",
                                         "Number of Lost Hits Vs Eta",
                                         TrackEtaHistoPar.getParameter<int32_t>("Xbins"),
                                         TrackEtaHistoPar.getParameter<double>("Xmin"),
                                         TrackEtaHistoPar.getParameter<double>("Xmax"),
                                         0.0,
                                         0.0,
                                         "g");
  if (!nLostHitsVsCosThetaH_)
    nLostHitsVsCosThetaH_ =
        iBook.bookProfile("nLostHitsVsCosTheta", "Number of Lost Hits Vs Cos(Theta)", 50, -1.0, 1.0, 0.0, 0.0, "g");
  if (!nLostHitsVsPhiH_)
    nLostHitsVsPhiH_ = iBook.bookProfile("nLostHitsVsPhi", "Number of Lost Hits Vs Phi", 100, -3.5, 3.5, 0.0, 0.0, "g");

  // On and off-track cluster properties
  if (!hOnTrkClusChargeThinH_)
    hOnTrkClusChargeThinH_ =
        iBook.book1D("hOnTrkClusChargeThin", "On-track Cluster Charge (Thin Sensor)", 100, 0, 1000);
  if (!hOnTrkClusWidthThinH_)
    hOnTrkClusWidthThinH_ = iBook.book1D("hOnTrkClusWidthThin", "On-track Cluster Width (Thin Sensor)", 20, -0.5, 19.5);
  if (!hOnTrkClusChargeThickH_)
    hOnTrkClusChargeThickH_ =
        iBook.book1D("hOnTrkClusChargeThick", "On-track Cluster Charge (Thick Sensor)", 100, 0, 1000);
  if (!hOnTrkClusWidthThickH_)
    hOnTrkClusWidthThickH_ =
        iBook.book1D("hOnTrkClusWidthThick", "On-track Cluster Width (Thick Sensor)", 20, -0.5, 19.5);

  if (!hOffTrkClusChargeThinH_)
    hOffTrkClusChargeThinH_ =
        iBook.book1D("hOffTrkClusChargeThin", "Off-track Cluster Charge (Thin Sensor)", 100, 0, 1000);
  if (!hOffTrkClusWidthThinH_)
    hOffTrkClusWidthThinH_ =
        iBook.book1D("hOffTrkClusWidthThin", "Off-track Cluster Width (Thin Sensor)", 20, -0.5, 19.5);
  if (!hOffTrkClusChargeThickH_)
    hOffTrkClusChargeThickH_ =
        iBook.book1D("hOffTrkClusChargeThick", "Off-track Cluster Charge (Thick Sensor)", 100, 0, 1000);
  if (!hOffTrkClusWidthThickH_)
    hOffTrkClusWidthThickH_ =
        iBook.book1D("hOffTrkClusWidthThick", "Off-track Cluster Width (Thick Sensor)", 20, -0.5, 19.5);
}
void StandaloneTrackMonitor::analyze(edm::Event const& iEvent, edm::EventSetup const& iSetup) {
  // Get event setup (to get global transformation)

  siStripClusterInfo_.initEvent(iSetup);

  // Primary vertex collection
  edm::Handle<reco::VertexCollection> vertexColl;
  iEvent.getByToken(vertexToken_, vertexColl);

  // Beam spot
  edm::Handle<reco::BeamSpot> beamSpot;
  iEvent.getByToken(bsToken_, beamSpot);

  // Track collection
  edm::Handle<reco::TrackCollection> tracks;
  iEvent.getByToken(trackToken_, tracks);

  // Access PU information
  double wfac = 1.0;  // for data
  if (!iEvent.isRealData()) {
    edm::Handle<std::vector<PileupSummaryInfo> > PupInfo;
    iEvent.getByToken(puSummaryToken_, PupInfo);

    if (verbose_)
      edm::LogInfo("StandaloneTrackMonitor") << "nPUColl = " << PupInfo->size();
    for (auto const& v : *PupInfo) {
      int bx = v.getBunchCrossing();
      if (bunchCrossingH_)
        bunchCrossingH_->Fill(bx);
      if (bx == 0) {
        if (nPUH_)
          nPUH_->Fill(v.getPU_NumInteractions());
        int ntrueInt = v.getTrueNumInteractions();
        if (trueNIntH_)
          trueNIntH_->Fill(ntrueInt);
        if (doPUCorrection_)
          if (ntrueInt > -1 && ntrueInt < int(vpu_.size()))
            wfac = vpu_.at(ntrueInt);
      }
    }
  }
  if (verbose_)
    edm::LogInfo("StandaloneTrackMonitor") << "PU reweight factor = " << wfac;

  if (!vertexColl.isValid())
    edm::LogError("DqmTrackStudy") << "Error! Failed to get reco::Vertex Collection, " << vertexTag_;
  if (haveAllHistograms_) {
    int nvtx = (vertexColl.isValid() ? vertexColl->size() : 0);
    nVertexH_->Fill(nvtx);
  }

  int ntracks = 0;
  if (tracks.isValid()) {
    edm::LogInfo("StandaloneTrackMonitor") << "Total # of Tracks: " << tracks->size();
    if (verbose_)
      edm::LogInfo("StandaloneTrackMonitor") << "Total # of Tracks: " << tracks->size();
    reco::Track::TrackQuality quality = reco::Track::qualityByName(trackQuality_);
    for (auto const& track : *tracks) {
      if (!track.quality(quality))
        continue;

      ++ntracks;

      double eta = track.eta();
      double theta = track.theta();
      double phi = track.phi();
      double pt = track.pt();

      const reco::HitPattern& hitp = track.hitPattern();
      double nValidTrackerHits = hitp.numberOfValidTrackerHits();
      nHitsVsEtaH_->Fill(eta, nValidTrackerHits);
      nHitsVsCosThetaH_->Fill(std::cos(theta), nValidTrackerHits);
      nHitsVsPhiH_->Fill(phi, nValidTrackerHits);
      nHitsVspTH_->Fill(pt, nValidTrackerHits);
      nHitsVsnVtxH_->Fill(vertexColl->size(), nValidTrackerHits);

      int nLostHits = track.numberOfLostHits();
      nLostHitsVspTH_->Fill(pt, nLostHits);
      nLostHitsVsEtaH_->Fill(eta, nLostHits);
      nLostHitsVsCosThetaH_->Fill(std::cos(theta), nLostHits);
      nLostHitsVsPhiH_->Fill(phi, nLostHits);

      double nValidPixelHits = hitp.numberOfValidPixelHits();
      double nValidStripHits = hitp.numberOfValidStripHits();
      double pixelLayersWithMeasurement = hitp.pixelLayersWithMeasurement();
      double stripLayersWithMeasurement = hitp.stripLayersWithMeasurement();

      if (haveAllHistograms_) {
        double etaError = track.etaError();
        double thetaError = track.thetaError();
        double phiError = track.phiError();
        double p = track.p();
        double ptError = track.ptError();
        double qoverp = track.qoverp();
        double qoverpError = track.qoverpError();
        double charge = track.charge();

        double trackerLayersWithMeasurement = hitp.trackerLayersWithMeasurement();

        double dxy = track.dxy(beamSpot->position());
        double dxyError = track.dxyError();
        double dz = track.dz(beamSpot->position());
        double dzError = track.dzError();

        double trkd0 = track.d0();
        double chi2 = track.chi2();
        double ndof = track.ndof();
        double vx = track.vx();
        double vy = track.vy();
        double vz = track.vz();

        // Fill the histograms
        trackEtaH_->Fill(eta, wfac);
        trackEtaerrH_->Fill(etaError, wfac);
        trackCosThetaH_->Fill(std::cos(theta), wfac);
        trackThetaerrH_->Fill(thetaError, wfac);
        trackPhiH_->Fill(phi, wfac);
        trackPhierrH_->Fill(phiError, wfac);
        trackPH_->Fill(p, wfac);
        trackPtH_->Fill(pt, wfac);
        if (pt <= 2)
          trackPtUpto2GeVH_->Fill(pt, wfac);
        if (pt >= 10)
          trackPtOver10GeVH_->Fill(pt, wfac);
        trackPterrH_->Fill(ptError, wfac);
        trackqOverpH_->Fill(qoverp, wfac);
        trackqOverperrH_->Fill(qoverpError, wfac);
        trackChargeH_->Fill(charge, wfac);
        trackChi2H_->Fill(chi2, wfac);
        trackd0H_->Fill(trkd0, wfac);
        tracknDOFH_->Fill(ndof, wfac);
        trackChi2bynDOFH_->Fill(chi2 / ndof, wfac);

        nlostHitsH_->Fill(nLostHits, wfac);
        nvalidTrackerHitsH_->Fill(nValidTrackerHits, wfac);
        nvalidPixelHitsH_->Fill(nValidPixelHits, wfac);
        nvalidStripHitsH_->Fill(nValidStripHits, wfac);

        trkLayerwithMeasurementH_->Fill(trackerLayersWithMeasurement, wfac);
        pixelLayerwithMeasurementH_->Fill(pixelLayersWithMeasurement, wfac);
        stripLayerwithMeasurementH_->Fill(stripLayersWithMeasurement, wfac);

        beamSpotXYposH_->Fill(dxy, wfac);
        beamSpotXYposerrH_->Fill(dxyError, wfac);
        beamSpotZposH_->Fill(dz, wfac);
        beamSpotZposerrH_->Fill(dzError, wfac);

        vertexXposH_->Fill(vx, wfac);
        vertexYposH_->Fill(vy, wfac);
        vertexZposH_->Fill(vz, wfac);
      }
      int nPixBarrel = 0, nPixEndcap = 0, nStripTIB = 0, nStripTOB = 0, nStripTEC = 0, nStripTID = 0;
      for (auto it = track.recHitsBegin(); it != track.recHitsEnd(); ++it) {
        const TrackingRecHit& hit = (**it);
        if (hit.isValid()) {
          if (hit.geographicalId().det() == DetId::Tracker) {
            int subdetId = hit.geographicalId().subdetId();
            if (subdetId == PixelSubdetector::PixelBarrel)
              ++nPixBarrel;
            else if (subdetId == PixelSubdetector::PixelEndcap)
              ++nPixEndcap;
            else if (subdetId == StripSubdetector::TIB)
              ++nStripTIB;
            else if (subdetId == StripSubdetector::TOB)
              ++nStripTOB;
            else if (subdetId == StripSubdetector::TEC)
              ++nStripTEC;
            else if (subdetId == StripSubdetector::TID)
              ++nStripTID;

            // Find on-track clusters
            processHit(hit, iSetup, wfac);
          }
        }
      }
      if (verbose_)
        edm::LogInfo("StandaloneTrackMonitor")
            << " >>> HITs: nPixBarrel: " << nPixBarrel << " nPixEndcap: " << nPixEndcap << " nStripTIB: " << nStripTIB
            << " nStripTOB: " << nStripTOB << " nStripTEC: " << nStripTEC << " nStripTID: " << nStripTID;
      if (haveAllHistograms_) {
        nPixBarrelH_->Fill(nPixBarrel, wfac);
        nPixEndcapH_->Fill(nPixEndcap, wfac);
        nStripTIBH_->Fill(nStripTIB, wfac);
        nStripTOBH_->Fill(nStripTOB, wfac);
        nStripTECH_->Fill(nStripTEC, wfac);
        nStripTIDH_->Fill(nStripTID, wfac);
      }
    }
  } else {
    edm::LogError("DqmTrackStudy") << "Error! Failed to get reco::Track collection, " << trackTag_;
  }
  if (haveAllHistograms_)
    nTracksH_->Fill(ntracks);

  // off track cluster properties
  processClusters(iEvent, iSetup, wfac);
}
void StandaloneTrackMonitor::processClusters(edm::Event const& iEvent, edm::EventSetup const& iSetup, double wfac) {
  // SiStripClusters
  edm::Handle<edmNew::DetSetVector<SiStripCluster> > clusterHandle;
  iEvent.getByToken(clusterToken_, clusterHandle);

  if (clusterHandle.isValid()) {
    // Loop on Dets
    for (edmNew::DetSetVector<SiStripCluster>::const_iterator dsvit = clusterHandle->begin();
         dsvit != clusterHandle->end();
         ++dsvit) {
      uint32_t detId = dsvit->id();
      std::map<uint32_t, std::set<const SiStripCluster*> >::iterator jt = clusterMap_.find(detId);
      bool detid_found = (jt != clusterMap_.end()) ? true : false;

      // Loop on Clusters
      for (edmNew::DetSet<SiStripCluster>::const_iterator clusit = dsvit->begin(); clusit != dsvit->end(); ++clusit) {
        if (detid_found) {
          std::set<const SiStripCluster*>& s = jt->second;
          if (s.find(&*clusit) != s.end())
            continue;
        }

        siStripClusterInfo_.setCluster(*clusit, detId);
        float charge = siStripClusterInfo_.charge();
        float width = siStripClusterInfo_.width();

        const GeomDetUnit* detUnit = tkGeom_->idToDetUnit(detId);
        float thickness = detUnit->surface().bounds().thickness();  // unit cm
        if (thickness > 0.035) {
          hOffTrkClusChargeThickH_->Fill(charge, wfac);
          hOffTrkClusWidthThickH_->Fill(width, wfac);
        } else {
          hOffTrkClusChargeThinH_->Fill(charge, wfac);
          hOffTrkClusWidthThinH_->Fill(width, wfac);
        }
      }
    }
  } else {
    edm::LogError("StandaloneTrackMonitor") << "ClusterCollection " << clusterTag_ << " not valid!!" << std::endl;
  }
}
void StandaloneTrackMonitor::processHit(const TrackingRecHit& recHit, edm::EventSetup const& iSetup, double wfac) {
  uint32_t detid = recHit.geographicalId();
  const GeomDetUnit* detUnit = tkGeom_->idToDetUnit(detid);
  float thickness = detUnit->surface().bounds().thickness();  // unit cm

  auto const& thit = static_cast<BaseTrackerRecHit const&>(recHit);
  if (!thit.isValid())
    return;

  auto const& clus = thit.firstClusterRef();
  if (!clus.isValid())
    return;
  if (!clus.isStrip())
    return;

  if (thit.isMatched()) {
    const SiStripMatchedRecHit2D& matchedHit = dynamic_cast<const SiStripMatchedRecHit2D&>(recHit);

    auto& clusterM = matchedHit.monoCluster();
    siStripClusterInfo_.setCluster(clusterM, detid);
    if (thickness > 0.035) {
      hOnTrkClusChargeThickH_->Fill(siStripClusterInfo_.charge(), wfac);
      hOnTrkClusWidthThickH_->Fill(siStripClusterInfo_.width(), wfac);
    } else {
      hOnTrkClusChargeThinH_->Fill(siStripClusterInfo_.charge(), wfac);
      hOnTrkClusWidthThinH_->Fill(siStripClusterInfo_.width(), wfac);
    }
    addClusterToMap(detid, &clusterM);

    auto& clusterS = matchedHit.stereoCluster();
    siStripClusterInfo_.setCluster(clusterS, detid);
    if (thickness > 0.035) {
      hOnTrkClusChargeThickH_->Fill(siStripClusterInfo_.charge(), wfac);
      hOnTrkClusWidthThickH_->Fill(siStripClusterInfo_.width(), wfac);
    } else {
      hOnTrkClusChargeThinH_->Fill(siStripClusterInfo_.charge(), wfac);
      hOnTrkClusWidthThinH_->Fill(siStripClusterInfo_.width(), wfac);
    }
    addClusterToMap(detid, &clusterS);
  } else {
    auto& cluster = clus.stripCluster();
    siStripClusterInfo_.setCluster(cluster, detid);
    if (thickness > 0.035) {
      hOnTrkClusChargeThickH_->Fill(siStripClusterInfo_.charge(), wfac);
      hOnTrkClusWidthThickH_->Fill(siStripClusterInfo_.width(), wfac);
    } else {
      hOnTrkClusChargeThinH_->Fill(siStripClusterInfo_.charge(), wfac);
      hOnTrkClusWidthThinH_->Fill(siStripClusterInfo_.width(), wfac);
    }
    addClusterToMap(detid, &cluster);
  }
}
void StandaloneTrackMonitor::addClusterToMap(uint32_t detid, const SiStripCluster* cluster) {
  std::map<uint32_t, std::set<const SiStripCluster*> >::iterator it = clusterMap_.find(detid);
  if (it == clusterMap_.end()) {
    std::set<const SiStripCluster*> s;
    s.insert(cluster);
    clusterMap_.insert(std::pair<uint32_t, std::set<const SiStripCluster*> >(detid, s));
  } else {
    std::set<const SiStripCluster*>& s = it->second;
    s.insert(cluster);
  }
}
// Define this as a plug-in
#include "FWCore/Framework/interface/MakerMacros.h"
DEFINE_FWK_MODULE(StandaloneTrackMonitor);