OfflineValidation/plugins/ShortenedTrackValidation.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    Alignment/OfflineValidation
0004 // Class:      ShortenedTrackValidation
0005 //
0006 /*
0007  *\class ShortenedTrackValidation ShortenedTrackValidation.cc Alignment/OfflineValidation/plugins/ShortenedTrackValidation.cc
0008
0009  Description: This module is meant to monitor the track pT resolution using the amputated tracks method, by comparing the performance using different alignments.
0010
0011  Implementation: The implemenation takes advantage of the existing implementation in the DQM/TrackingMonitorSource.
0012
0013 */
0014 //
0015 // Original Author:  Marco Musich
0016 //         Created:  Fri, 05 Jan 2023 11:41:00 GMT
0017 //
0018 //
0019
0020 // ROOT includes files
0021 #include "TMath.h"
0022 #include "TFile.h"
0023 #include "TH1D.h"
0024 #include "TH1I.h"
0025 #include "TH2D.h"
0026 #include "TProfile.h"
0027 #include "TLorentzVector.h"
0028
0029 // standard includes
0030 #include <fmt/printf.h>
0031
0032 // user includes
0033 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0034 #include "DataFormats/TrackReco/interface/Track.h"
0035 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
0036 #include "DataFormats/TrackerRecHit2D/interface/ProjectedSiStripRecHit2D.h"
0037 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
0038 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit1D.h"
0039 #include "DataFormats/TrackerRecHit2D/interface/SiTrackerMultiRecHit.h"
0040 #include "DataFormats/VertexReco/interface/Vertex.h"
0041 #include "FWCore/Framework/interface/Event.h"
0042 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0043 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0044 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0045 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0046 #include "FWCore/ServiceRegistry/interface/Service.h"
0047 #include "FWCore/Utilities/interface/transform.h"  // for edm::vector_transform
0048
0049 #define CREATE_HIST_1D(varname, nbins, first, last, fs) fs.make<TH1D>(#varname, #varname, nbins, first, last)
0050
0051 #define CREATE_HIST_2D(varname, nbins, first, last, fs) \
0052   fs.make<TH2D>(#varname, #varname, nbins, first, last, nbins, first, last)
0053
0054 const int kBPIX = PixelSubdetector::PixelBarrel;
0055 const int kFPIX = PixelSubdetector::PixelEndcap;
0056
0057 class ShortenedTrackValidation : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0058   class trackingMon {
0059   public:
0060     trackingMon() {}
0061     ~trackingMon() = default;
0062
0063     void book(const TFileDirectory &fs) {
0064       h_chi2ndof = CREATE_HIST_1D(h_chi2ndof, 100, 0.0, 10.0, fs);
0065       h_trkQuality = CREATE_HIST_1D(h_trkQuality, 6, -1, 5, fs);
0066       h_trkAlgo = CREATE_HIST_1D(h_trkAlgo, reco::TrackBase::algoSize, 0.0, double(reco::TrackBase::algoSize), fs);
0067       h_trkOriAlgo =
0068           CREATE_HIST_1D(h_trkOriAlgo, reco::TrackBase::algoSize, 0.0, double(reco::TrackBase::algoSize), fs);
0069       h_P = CREATE_HIST_1D(h_P, 100, 0.0, 200.0, fs);
0070       h_Pt = CREATE_HIST_1D(h_Pt, 100, 0.0, 100.0, fs);
0071       h_nHit = CREATE_HIST_1D(h_nHit, 50, -0.5, 49.5, fs);
0072       h_nHit2D = CREATE_HIST_1D(h_nHit2D, 20, -0.5, 19.5, fs);
0073       h_Charge = CREATE_HIST_1D(h_Charge, 3, -1.5, 1.5, fs);
0074       h_QoverP = CREATE_HIST_1D(h_QoverP, 100, -1.0, 1.0, fs);
0075       h_QoverPZoom = CREATE_HIST_1D(h_QoverPZoom, 100, -0.1, 0.1, fs);
0076       h_Eta = CREATE_HIST_1D(h_Eta, 100, -3., 3., fs);
0077       h_Phi = CREATE_HIST_1D(h_Phi, 100, -M_PI, M_PI, fs);
0078       h_vx = CREATE_HIST_1D(h_vx, 100, -0.5, 0.5, fs);
0079       h_vy = CREATE_HIST_1D(h_vy, 100, -0.5, 0.5, fs);
0080       h_vz = CREATE_HIST_1D(h_vz, 100, -20.0, 20.0, fs);
0081       h_d0 = CREATE_HIST_1D(h_d0, 100, -0.5, 0.5, fs);
0082       h_dz = CREATE_HIST_1D(h_dz, 100, -20.0, 20.0, fs);
0083       h_dxy = CREATE_HIST_1D(h_dxy, 100, -0.5, 0.5, fs);
0084       h_nhpxb = CREATE_HIST_1D(h_nhpxb, 10, -0.5, 9.5, fs);
0085       h_nhpxe = CREATE_HIST_1D(h_nhpxe, 10, -0.5, 9.5, fs);
0086       h_nhTIB = CREATE_HIST_1D(h_nhTIB, 20, -0.5, 19.5, fs);
0087       h_nhTID = CREATE_HIST_1D(h_nhTID, 20, -0.5, 19.5, fs);
0088       h_nhTOB = CREATE_HIST_1D(h_nhTOB, 20, -0.5, 19.5, fs);
0089       h_nhTEC = CREATE_HIST_1D(h_nhTEC, 20, -0.5, 19.5, fs);
0090       h_dxyBS = CREATE_HIST_1D(h_dxyBS, 100, -0.05, 0.05, fs);
0091       h_d0BS = CREATE_HIST_1D(h_d0BS, 100, -0.05, 0.05, fs);
0092       h_dzBS = CREATE_HIST_1D(h_dzBS, 100, -20.0, 20., fs);
0093       h_dxyPV = CREATE_HIST_1D(h_dxyPV, 100, -0.05, 0.05, fs);
0094       h_d0PV = CREATE_HIST_1D(h_d0PV, 100, -0.05, 0.05, fs);
0095       h_dzPV = CREATE_HIST_1D(h_dzPV, 100, -0.05, 0.05, fs);
0096
0097       edm::LogInfo("trackingMonitoring") << "done booking";
0098     }
0099
0100     //____________________________________________________________
0101     int trackQual(const reco::Track &track) {
0102       int myquality = -99;
0103       if (track.quality(reco::TrackBase::undefQuality))
0104         myquality = -1;
0105       if (track.quality(reco::TrackBase::loose))
0106         myquality = 0;
0107       if (track.quality(reco::TrackBase::tight))
0108         myquality = 1;
0109       if (track.quality(reco::TrackBase::highPurity))
0110         myquality = 2;
0111       if (track.quality(reco::TrackBase::goodIterative))
0112         myquality = 3;
0113
0114       return myquality;
0115     }
0116
0117     //____________________________________________________________
0118     static bool isHit2D(const TrackingRecHit &hit) {
0119       if (hit.dimension() < 2) {
0120         return false;  // some (muon...) stuff really has RecHit1D
0121       } else {
0122         const DetId detId(hit.geographicalId());
0123         if (detId.det() == DetId::Tracker) {
0124           if (detId.subdetId() == kBPIX || detId.subdetId() == kFPIX) {
0125             return true;  // pixel is always 2D
0126           } else {        // should be SiStrip now
0127             if (dynamic_cast<const SiStripRecHit2D *>(&hit))
0128               return false;  // normal hit
0129             else if (dynamic_cast<const SiStripMatchedRecHit2D *>(&hit))
0130               return true;  // matched is 2D
0131             else if (dynamic_cast<const ProjectedSiStripRecHit2D *>(&hit))
0132               return false;  // crazy hit...
0133             else {
0134               edm::LogError("UnknownType") << "@SUB=CalibrationTrackSelector::isHit2D"
0135                                            << "Tracker hit not in pixel and neither SiStripRecHit2D nor "
0136                                            << "SiStripMatchedRecHit2D nor ProjectedSiStripRecHit2D.";
0137               return false;
0138             }
0139           }
0140         } else {  // not tracker??
0141           edm::LogWarning("DetectorMismatch") << "@SUB=CalibrationTrackSelector::isHit2D"
0142                                               << "Hit not in tracker with 'official' dimension >=2.";
0143           return true;  // dimension() >= 2 so accept that...
0144         }
0145       }
0146       // never reached...
0147     }
0148
0149     //____________________________________________________________
0150     unsigned int count2DHits(const reco::Track &track) {
0151       unsigned int nHit2D = 0;
0152       for (auto iHit = track.recHitsBegin(); iHit != track.recHitsEnd(); ++iHit) {
0153         if (isHit2D(**iHit)) {
0154           ++nHit2D;
0155         }
0156       }
0157       return nHit2D;
0158     }
0159
0160     //____________________________________________________________
0161     void fill(const reco::Track &track, const reco::BeamSpot &beamSpot, const reco::Vertex &pvtx) {
0162       h_chi2ndof->Fill(track.normalizedChi2());
0163       h_trkQuality->Fill(trackQual(track));
0164       h_trkAlgo->Fill(static_cast<float>(track.algo()));
0165       h_trkOriAlgo->Fill(static_cast<float>(track.originalAlgo()));
0166       h_P->Fill(track.p());
0167       h_Pt->Fill(track.pt());
0168       h_nHit->Fill(track.numberOfValidHits());
0169       h_nHit2D->Fill(count2DHits(track));
0170       h_Charge->Fill(track.charge());
0171       h_QoverP->Fill(track.qoverp());
0172       h_QoverPZoom->Fill(track.qoverp());
0173       h_Eta->Fill(track.eta());
0174       h_Phi->Fill(track.phi());
0175       h_vx->Fill(track.vx());
0176       h_vy->Fill(track.vy());
0177       h_vz->Fill(track.vz());
0178       h_d0->Fill(track.d0());
0179       h_dz->Fill(track.dz());
0180       h_dxy->Fill(track.dxy());
0181       h_nhpxb->Fill(track.hitPattern().numberOfValidPixelBarrelHits());
0182       h_nhpxe->Fill(track.hitPattern().numberOfValidPixelEndcapHits());
0183       h_nhTIB->Fill(track.hitPattern().numberOfValidStripTIBHits());
0184       h_nhTID->Fill(track.hitPattern().numberOfValidStripTIDHits());
0185       h_nhTOB->Fill(track.hitPattern().numberOfValidStripTOBHits());
0186       h_nhTEC->Fill(track.hitPattern().numberOfValidStripTECHits());
0187
0188       math::XYZPoint BS(beamSpot.x0(), beamSpot.y0(), beamSpot.z0());
0189       h_dxyBS->Fill(track.dxy(BS));
0190       h_d0BS->Fill(-track.dxy(BS));
0191       h_dzBS->Fill(track.dz(BS));
0192
0193       math::XYZPoint PV(pvtx.x(), pvtx.y(), pvtx.z());
0194       h_dxyPV->Fill(track.dxy(PV));
0195       h_d0PV->Fill(-track.dxy(PV));
0196       h_dzPV->Fill(track.dz(PV));
0197     }
0198
0199   private:
0200     TH1D *h_chi2ndof;
0201     TH1D *h_trkQuality;
0202     TH1D *h_trkAlgo;
0203     TH1D *h_trkOriAlgo;
0204     TH1D *h_P;
0205     TH1D *h_Pt;
0206     TH1D *h_nHit;
0207     TH1D *h_nHit2D;
0208     TH1D *h_Charge;
0209     TH1D *h_QoverP;
0210     TH1D *h_QoverPZoom;
0211     TH1D *h_Eta;
0212     TH1D *h_Phi;
0213     TH1D *h_vx;
0214     TH1D *h_vy;
0215     TH1D *h_vz;
0216     TH1D *h_d0;
0217     TH1D *h_dz;
0218     TH1D *h_dxy;
0219     TH1D *h_nhpxb;
0220     TH1D *h_nhpxe;
0221     TH1D *h_nhTIB;
0222     TH1D *h_nhTID;
0223     TH1D *h_nhTOB;
0224     TH1D *h_nhTEC;
0225     TH1D *h_dxyBS;
0226     TH1D *h_d0BS;
0227     TH1D *h_dzBS;
0228     TH1D *h_dxyPV;
0229     TH1D *h_d0PV;
0230     TH1D *h_dzPV;
0231   };
0232
0233   class trackComparator {
0234   public:
0235     trackComparator() {}
0236     ~trackComparator() = default;
0237
0238     //__________________________________________________
0239     void book(const TFileDirectory &fs) {
0240       h2_chi2ndof = CREATE_HIST_2D(h2_chi2ndof, 100, 0.0, 10.0, fs);
0241       h2_trkAlgo = CREATE_HIST_2D(h2_trkAlgo, reco::TrackBase::algoSize, 0.0, double(reco::TrackBase::algoSize), fs);
0242       h2_trkOriAlgo =
0243           CREATE_HIST_2D(h2_trkOriAlgo, reco::TrackBase::algoSize, 0.0, double(reco::TrackBase::algoSize), fs);
0244       h2_P = CREATE_HIST_2D(h2_P, 100, 0.0, 200.0, fs);
0245       h2_Pt = CREATE_HIST_2D(h2_Pt, 100, 0.0, 100.0, fs);
0246       h2_nHit = CREATE_HIST_2D(h2_nHit, 50, -0.5, 49.5, fs);
0247       h2_Charge = CREATE_HIST_2D(h2_Charge, 3, -1.5, 1.5, fs);
0248       h2_QoverPZoom = CREATE_HIST_2D(h2_QoverPZoom, 100, -0.1, 0.1, fs);
0249       h2_Eta = CREATE_HIST_2D(h2_Eta, 100, -3., 3., fs);
0250       h2_Phi = CREATE_HIST_2D(h2_Phi, 100, -M_PI, M_PI, fs);
0251       h2_vx = CREATE_HIST_2D(h2_vx, 100, -0.5, 0.5, fs);
0252       h2_vy = CREATE_HIST_2D(h2_vy, 100, -0.5, 0.5, fs);
0253       h2_vz = CREATE_HIST_2D(h2_vz, 100, -20.0, 20.0, fs);
0254       h2_d0 = CREATE_HIST_2D(h2_d0, 100, -0.5, 0.5, fs);
0255       h2_dz = CREATE_HIST_2D(h2_dz, 100, -20.0, 20.0, fs);
0256       h2_nhpxb = CREATE_HIST_2D(h2_nhpxb, 10, -0.5, 9.5, fs);
0257       h2_nhpxe = CREATE_HIST_2D(h2_nhpxe, 10, -0.5, 9.5, fs);
0258       h2_nhTIB = CREATE_HIST_2D(h2_nhTIB, 20, -0.5, 19.5, fs);
0259       h2_nhTID = CREATE_HIST_2D(h2_nhTID, 20, -0.5, 19.5, fs);
0260       h2_nhTOB = CREATE_HIST_2D(h2_nhTOB, 20, -0.5, 19.5, fs);
0261       h2_nhTEC = CREATE_HIST_2D(h2_nhTEC, 20, -0.5, 19.5, fs);
0262     }
0263
0264     //__________________________________________________
0265     void fill(const reco::Track &tk1, const reco::Track &tk2) {
0266       h2_chi2ndof->Fill(tk1.normalizedChi2(), tk2.normalizedChi2());
0267       h2_trkAlgo->Fill(static_cast<float>(tk1.algo()), static_cast<float>(tk2.algo()));
0268       h2_trkOriAlgo->Fill(static_cast<float>(tk1.originalAlgo()), static_cast<float>(tk2.originalAlgo()));
0269       h2_P->Fill(tk1.p(), tk2.p());
0270       h2_Pt->Fill(tk1.pt(), tk2.p());
0271       h2_nHit->Fill(tk1.numberOfValidHits(), tk2.numberOfValidHits());
0272       h2_Charge->Fill(tk1.charge(), tk2.charge());
0273       h2_QoverPZoom->Fill(tk1.qoverp(), tk2.qoverp());
0274       h2_Eta->Fill(tk1.eta(), tk2.eta());
0275       h2_Phi->Fill(tk1.phi(), tk2.phi());
0276       h2_vx->Fill(tk1.vx(), tk2.vx());
0277       h2_vy->Fill(tk1.vy(), tk2.vy());
0278       h2_vz->Fill(tk1.vz(), tk2.vz());
0279       h2_d0->Fill(tk1.d0(), tk2.d0());
0280       h2_dz->Fill(tk2.dz(), tk2.dz());
0281       h2_nhpxb->Fill(tk1.hitPattern().numberOfValidPixelBarrelHits(), tk2.hitPattern().numberOfValidPixelBarrelHits());
0282       h2_nhpxe->Fill(tk1.hitPattern().numberOfValidPixelEndcapHits(), tk2.hitPattern().numberOfValidPixelEndcapHits());
0283       h2_nhTIB->Fill(tk1.hitPattern().numberOfValidStripTIBHits(), tk2.hitPattern().numberOfValidStripTIBHits());
0284       h2_nhTID->Fill(tk1.hitPattern().numberOfValidStripTIDHits(), tk2.hitPattern().numberOfValidStripTIDHits());
0285       h2_nhTOB->Fill(tk1.hitPattern().numberOfValidStripTOBHits(), tk2.hitPattern().numberOfValidStripTOBHits());
0286       h2_nhTEC->Fill(tk1.hitPattern().numberOfValidStripTECHits(), tk2.hitPattern().numberOfValidStripTECHits());
0287     }
0288
0289   private:
0290     TH2D *h2_chi2ndof;
0291     TH2D *h2_trkAlgo;
0292     TH2D *h2_trkOriAlgo;
0293     TH2D *h2_P;
0294     TH2D *h2_Pt;
0295     TH2D *h2_nHit;
0296     TH2D *h2_Charge;
0297     TH2D *h2_QoverPZoom;
0298     TH2D *h2_Eta;
0299     TH2D *h2_Phi;
0300     TH2D *h2_vx;
0301     TH2D *h2_vy;
0302     TH2D *h2_vz;
0303     TH2D *h2_d0;
0304     TH2D *h2_dz;
0305     TH2D *h2_nhpxb;
0306     TH2D *h2_nhpxe;
0307     TH2D *h2_nhTIB;
0308     TH2D *h2_nhTID;
0309     TH2D *h2_nhTOB;
0310     TH2D *h2_nhTEC;
0311   };
0312
0313 public:
0314   explicit ShortenedTrackValidation(const edm::ParameterSet &);
0315   ~ShortenedTrackValidation() override = default;
0316   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0317
0318 private:
0319   template <typename T, typename... Args>
0320   T *book(const TFileDirectory &dir, const Args &...args) const;
0321   void beginJob() override;
0322   void analyze(edm::Event const &iEvent, edm::EventSetup const &iSetup) override;
0323
0324   // ----------member data ---------------------------
0325   edm::Service<TFileService> fs_;
0326   const std::string folderName_;
0327   const std::vector<std::string> hitsRemain_;
0328   const double minTracksEta_;
0329   const double maxTracksEta_;
0330   const double minTracksPt_;
0331   const double maxTracksPt_;
0332
0333   const double maxDr_;
0334   const edm::InputTag tracksTag_;
0335   const std::vector<edm::InputTag> tracksRerecoTag_;
0336   const edm::InputTag BeamSpotTag_;
0337   const edm::InputTag VerticesTag_;
0338   const edm::EDGetTokenT<std::vector<reco::Track>> tracksToken_;
0339   const std::vector<edm::EDGetTokenT<std::vector<reco::Track>>> tracksRerecoToken_;
0340   const edm::EDGetTokenT<reco::BeamSpot> beamspotToken_;
0341   const edm::EDGetTokenT<reco::VertexCollection> vertexToken_;
0342
0343   // monitoring histograms
0344   std::vector<TH1F *> histsPtRatioAll_;
0345   std::vector<TH1F *> histsPtDiffAll_;
0346   std::vector<TH1F *> histsEtaDiffAll_;
0347   std::vector<TH1F *> histsPhiDiffAll_;
0348   std::vector<TH2F *> histsPtRatioVsDeltaRAll_;
0349   std::vector<TH1F *> histsDeltaPtOverPtAll_;
0350   std::vector<TH1F *> histsPtAll_;
0351   std::vector<TH1F *> histsNhitsAll_;
0352   std::vector<TH1F *> histsDeltaRAll_;
0353
0354   trackingMon originalTrack;
0355   std::vector<trackComparator *> comparators_;
0356   static constexpr double muMass = 0.105658;
0357 };
0358
0359 // -----------------------------
0360 // constructors and destructor
0361 // -----------------------------
0362 ShortenedTrackValidation::ShortenedTrackValidation(const edm::ParameterSet &ps)
0363     : folderName_(ps.getUntrackedParameter<std::string>("folderName", "TrackRefitting")),
0364       hitsRemain_(ps.getUntrackedParameter<std::vector<std::string>>("hitsRemainInput")),
0365       minTracksEta_(ps.getUntrackedParameter<double>("minTracksEtaInput", 0.0)),
0366       maxTracksEta_(ps.getUntrackedParameter<double>("maxTracksEtaInput", 2.2)),
0367       minTracksPt_(ps.getUntrackedParameter<double>("minTracksPtInput", 15.0)),
0368       maxTracksPt_(ps.getUntrackedParameter<double>("maxTracksPtInput", 99999.9)),
0369       maxDr_(ps.getUntrackedParameter<double>("maxDrInput", 0.01)),
0370       tracksTag_(ps.getUntrackedParameter<edm::InputTag>("tracksInputTag", edm::InputTag("generalTracks", "", "DQM"))),
0371       tracksRerecoTag_(ps.getUntrackedParameter<std::vector<edm::InputTag>>("tracksRerecoInputTag")),
0372       BeamSpotTag_(ps.getUntrackedParameter<edm::InputTag>("BeamSpotTag", edm::InputTag("offlineBeamSpot"))),
0373       VerticesTag_(ps.getUntrackedParameter<edm::InputTag>("VerticesTag", edm::InputTag("offlinePrimaryVertices"))),
0374       tracksToken_(consumes<std::vector<reco::Track>>(tracksTag_)),
0375       tracksRerecoToken_(edm::vector_transform(
0376           tracksRerecoTag_, [this](edm::InputTag const &tag) { return consumes<std::vector<reco::Track>>(tag); })),
0377       beamspotToken_(consumes<reco::BeamSpot>(BeamSpotTag_)),
0378       vertexToken_(consumes<reco::VertexCollection>(VerticesTag_)) {
0379   usesResource(TFileService::kSharedResource);
0380   histsPtRatioAll_.clear();
0381   histsPtDiffAll_.clear();
0382   histsEtaDiffAll_.clear();
0383   histsPhiDiffAll_.clear();
0384   histsPtRatioVsDeltaRAll_.clear();
0385   histsDeltaPtOverPtAll_.clear();
0386   histsPtAll_.clear();
0387   histsNhitsAll_.clear();
0388   histsDeltaRAll_.clear();
0389   comparators_.clear();
0390
0391   comparators_.reserve(hitsRemain_.size());
0392   for (unsigned int i = 0; i < hitsRemain_.size(); ++i) {
0393     comparators_.push_back(new trackComparator());
0394   }
0395 }
0396
0397 //__________________________________________________________________________________
0398 template <typename T, typename... Args>
0399 T *ShortenedTrackValidation::book(const TFileDirectory &dir, const Args &...args) const {
0400   T *t = dir.make<T>(args...);
0401   return t;
0402 }
0403
0404 //__________________________________________________________________________________
0405 void ShortenedTrackValidation::beginJob() {
0406   std::string currentFolder = folderName_ + "/Resolutions";
0407   TFileDirectory ShortTrackResolution = fs_->mkdir(currentFolder);
0408   currentFolder = folderName_ + "/Tracks";
0409   TFileDirectory TrackQuals = fs_->mkdir(currentFolder);
0410
0411   for (unsigned int i = 0; i < hitsRemain_.size(); ++i) {
0412     histsPtRatioAll_.push_back(
0413         book<TH1F>(ShortTrackResolution,
0414                    fmt::sprintf("trackPtRatio_%s", hitsRemain_[i]).c_str(),
0415                    fmt::sprintf("Short Track p_{T} / Full Track p_{T} - %s layers;p_{T}^{short}/p_{T}^{full};n. tracks",
0416                                 hitsRemain_[i])
0417                        .c_str(),
0418                    100,
0419                    0.5,
0420                    1.5));
0421
0422     histsPtDiffAll_.push_back(book<TH1F>(
0423         ShortTrackResolution,
0424         fmt::sprintf("trackPtDiff_%s", hitsRemain_[i]).c_str(),
0425         fmt::sprintf("Short Track p_{T} - Full Track p_{T} - %s layers;p_{T}^{short} - p_{T}^{full} [GeV];n. tracks",
0426                      hitsRemain_[i])
0427             .c_str(),
0428         100,
0429         -10.,
0430         10.));
0431
0432     histsEtaDiffAll_.push_back(
0433         book<TH1F>(ShortTrackResolution,
0434                    fmt::sprintf("trackEtaDiff_%s", hitsRemain_[i]).c_str(),
0435                    fmt::sprintf("Short Track #eta - Full Track #eta - %s layers;#eta^{short} - #eta^{full};n. tracks",
0436                                 hitsRemain_[i])
0437                        .c_str(),
0438                    100,
0439                    -0.001,
0440                    0.001));
0441
0442     histsPhiDiffAll_.push_back(
0443         book<TH1F>(ShortTrackResolution,
0444                    fmt::sprintf("trackPhiDiff_%s", hitsRemain_[i]).c_str(),
0445                    fmt::sprintf("Short Track #phi - Full Track #phi - %s layers;#phi^{short} - #phi^{full};n. tracks",
0446                                 hitsRemain_[i])
0447                        .c_str(),
0448                    100,
0449                    -0.001,
0450                    0.001));
0451
0452     histsPtRatioVsDeltaRAll_.push_back(
0453         book<TH2F>(ShortTrackResolution,
0454                    fmt::sprintf("trackPtRatioVsDeltaR_%s", hitsRemain_[i]).c_str(),
0455                    fmt::sprintf("Short Track p_{T} / Full Track p_{T} - %s layers vs "
0456                                 "#DeltaR;#DeltaR(short,full);p_{T}^{short}/p_{T}^{full} [GeV];n. tracks",
0457                                 hitsRemain_[i])
0458                        .c_str(),
0459                    100,
0460                    0.,
0461                    0.01,
0462                    101,
0463                    -0.05,
0464                    2.05));
0465
0466     histsDeltaPtOverPtAll_.push_back(
0467         book<TH1F>(ShortTrackResolution,
0468                    fmt::sprintf("trackDeltaPtOverPt_%s", hitsRemain_[i]).c_str(),
0469                    fmt::sprintf("Short Track p_{T} - Full Track p_{T} / Full Track p_{T} - %s layers;p_{T}^{short} - "
0470                                 "p_{T}^{full} / p^{full}_{T};n. tracks",
0471                                 hitsRemain_[i])
0472                        .c_str(),
0473                    101,
0474                    -5.,
0475                    5.));
0476
0477     histsPtAll_.push_back(
0478         book<TH1F>(TrackQuals,
0479                    fmt::sprintf("trackPt_%s", hitsRemain_[i]).c_str(),
0480                    fmt::sprintf("Short Track p_{T} - %s layers;p_{T}^{short} [GeV];n. tracks", hitsRemain_[i]).c_str(),
0481                    100,
0482                    0.,
0483                    100.));
0484
0485     histsNhitsAll_.push_back(
0486         book<TH1F>(TrackQuals,
0487                    fmt::sprintf("trackNhits_%s", hitsRemain_[i]).c_str(),
0488                    fmt::sprintf("Short Track n. hits - %s layers; n. hits per track;n. tracks", hitsRemain_[i]).c_str(),
0489                    20,
0490                    -0.5,
0491                    19.5));
0492
0493     histsDeltaRAll_.push_back(book<TH1F>(
0494         TrackQuals,
0495         fmt::sprintf("trackDeltaR_%s", hitsRemain_[i]).c_str(),
0496         fmt::sprintf("Short Track / Long Track #DeltaR %s layers;#DeltaR(short,long);n. tracks", hitsRemain_[i]).c_str(),
0497         100,
0498         0.,
0499         0.005));
0500
0501     currentFolder = fmt::sprintf("%s/Compare_%sHit", folderName_, hitsRemain_[i]);
0502     comparators_[i]->book(fs_->mkdir(currentFolder));
0503   }
0504
0505   currentFolder = folderName_ + "/OriginalTrack";
0506   TFileDirectory original = fs_->mkdir(currentFolder);
0507   originalTrack.book(original);
0508 }
0509
0510 //__________________________________________________________________________________
0511 void ShortenedTrackValidation::analyze(edm::Event const &iEvent, edm::EventSetup const &iSetup) {
0512   const auto &tracks = iEvent.getHandle(tracksToken_);
0513
0514   if (!tracks.isValid()) {
0515     edm::LogError("ShortenedTrackValidation") << "Missing input track collection " << tracksTag_.encode() << std::endl;
0516     return;
0517   }
0518
0519   reco::BeamSpot beamSpot;
0520   edm::Handle<reco::BeamSpot> beamSpotHandle = iEvent.getHandle(beamspotToken_);
0521   if (beamSpotHandle.isValid()) {
0522     beamSpot = *beamSpotHandle;
0523   } else {
0524     beamSpot = reco::BeamSpot();
0525   }
0526
0527   reco::Vertex pvtx;
0528   edm::Handle<reco::VertexCollection> vertexHandle = iEvent.getHandle(vertexToken_);
0529   if (vertexHandle.isValid()) {
0530     pvtx = (*vertexHandle).at(0);
0531   } else {
0532     pvtx = reco::Vertex();
0533   }
0534
0535   // the original long track
0536   for (const auto &track : *tracks) {
0537     const reco::HitPattern &hp = track.hitPattern();
0538     if (int(int(hp.numberOfValidHits()) - int(hp.numberOfAllHits(reco::HitPattern::TRACK_HITS))) != 0) {
0539       break;
0540     }
0541
0542     // fill the original track properties monitoring
0543     originalTrack.fill(track, beamSpot, pvtx);
0544
0545     TLorentzVector tvec;
0546     tvec.SetPtEtaPhiM(track.pt(), track.eta(), track.phi(), muMass);
0547
0548     int i = 0;  // token index
0549     // loop on the re-recoed shortened track collections
0550     for (const auto &token : tracksRerecoToken_) {
0551       const auto &tracks_rereco = iEvent.getHandle(token);
0552
0553       for (const auto &track_rereco : *tracks_rereco) {
0554         TLorentzVector trerecovec;
0555         trerecovec.SetPtEtaPhiM(track_rereco.pt(), track_rereco.eta(), track_rereco.phi(), 0.0);
0556         double deltaR = tvec.DeltaR(trerecovec);
0557
0558         if (deltaR < maxDr_) {
0559           if (track_rereco.pt() >= minTracksPt_ && track_rereco.pt() <= maxTracksPt_ &&
0560               std::abs(track_rereco.eta()) >= minTracksEta_ && std::abs(track_rereco.eta()) <= maxTracksEta_) {
0561             // fill the 2D comparisons per track
0562             comparators_[i]->fill(track, track_rereco);
0563
0564             histsPtRatioAll_[i]->Fill(1.0 * track_rereco.pt() / track.pt());
0565             histsPtDiffAll_[i]->Fill(track_rereco.pt() - track.pt());
0566             histsDeltaPtOverPtAll_[i]->Fill((track_rereco.pt() - track.pt()) / track.pt());
0567             histsEtaDiffAll_[i]->Fill(track_rereco.eta() - track.eta());
0568             histsPhiDiffAll_[i]->Fill(track_rereco.phi() - track.phi());
0569             histsPtRatioVsDeltaRAll_[i]->Fill(deltaR, track_rereco.pt() / track.pt());
0570             histsPtAll_[i]->Fill(track_rereco.pt());
0571             histsNhitsAll_[i]->Fill(track_rereco.numberOfValidHits());
0572             histsDeltaRAll_[i]->Fill(deltaR);
0573           }
0574         }
0575       }
0576       ++i;
0577     }
0578   }
0579 }
0580
0581 //__________________________________________________________________________________
0582 void ShortenedTrackValidation::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0583   edm::ParameterSetDescription desc;
0584   desc.addUntracked<std::string>("folderName", "TrackRefitting");
0585   desc.addUntracked<std::vector<std::string>>("hitsRemainInput", {});
0586   desc.addUntracked<double>("minTracksEtaInput", 0.0);
0587   desc.addUntracked<double>("maxTracksEtaInput", 2.2);
0588   desc.addUntracked<double>("minTracksPtInput", 15.0);
0589   desc.addUntracked<double>("maxTracksPtInput", 99999.9);
0590   desc.addUntracked<double>("maxDrInput", 0.01);
0591   desc.addUntracked<edm::InputTag>("tracksInputTag", edm::InputTag("generalTracks", "", "DQM"));
0592   desc.addUntracked<std::vector<edm::InputTag>>("tracksRerecoInputTag", {});
0593   desc.addUntracked<edm::InputTag>("BeamSpotTag", edm::InputTag("offlineBeamSpot"));
0594   desc.addUntracked<edm::InputTag>("VerticesTag", edm::InputTag("offlinePrimaryVertices"));
0595   descriptions.addWithDefaultLabel(desc);
0596 }
0597
0598 // Define this as a plug-in
0599 #include "FWCore/Framework/interface/MakerMacros.h"
0600 DEFINE_FWK_MODULE(ShortenedTrackValidation);