Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2025-03-31 22:27:14

0001 #ifndef MuonTrackValidatorBase_h
0002 #define MuonTrackValidatorBase_h
0003 
0004 /** \class MuonTrackValidatorBase
0005 * Base class for analyzers that produces histograms to validate Muon Track Reconstruction performances
0006 *
0007 */
0008 
0009 #include <memory>
0010 
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "FWCore/Framework/interface/ESHandle.h"
0013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0014 #include "FWCore/Framework/interface/ConsumesCollector.h"
0015 
0016 #include "MagneticField/Engine/interface/MagneticField.h"
0017 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0018 
0019 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0020 #include "DataFormats/TrackReco/interface/Track.h"
0021 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0022 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
0023 #include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"
0024 
0025 #include "DQMServices/Core/interface/DQMStore.h"
0026 #include "FWCore/ServiceRegistry/interface/Service.h"
0027 
0028 #include "SimTracker/Common/interface/TrackingParticleSelector.h"
0029 #include "CommonTools/RecoAlgos/interface/CosmicTrackingParticleSelector.h"
0030 
0031 #include <iostream>
0032 #include <sstream>
0033 #include <string>
0034 #include <TH1F.h>
0035 #include <TH2F.h>
0036 
0037 struct HistoParams {
0038   HistoParams(edm::ParameterSet muonHistoParameters) {
0039     mindR = muonHistoParameters.getParameter<double>("mindR");
0040     maxdR = muonHistoParameters.getParameter<double>("maxdR");
0041     nintdR = muonHistoParameters.getParameter<int>("nintdR");
0042     minEta = muonHistoParameters.getParameter<double>("minEta");
0043     maxEta = muonHistoParameters.getParameter<double>("maxEta");
0044     nintEta = muonHistoParameters.getParameter<int>("nintEta");
0045     minPt = muonHistoParameters.getParameter<double>("minPt");
0046     maxPt = muonHistoParameters.getParameter<double>("maxPt");
0047     nintPt = muonHistoParameters.getParameter<int>("nintPt");
0048     useLogPt = muonHistoParameters.getUntrackedParameter<bool>("useLogPt", false);
0049     minNHit = muonHistoParameters.getParameter<double>("minNHit");
0050     maxNHit = muonHistoParameters.getParameter<double>("maxNHit");
0051     nintNHit = muonHistoParameters.getParameter<int>("nintNHit");
0052     //
0053     minDTHit = muonHistoParameters.getParameter<double>("minDTHit");
0054     maxDTHit = muonHistoParameters.getParameter<double>("maxDTHit");
0055     nintDTHit = muonHistoParameters.getParameter<int>("nintDTHit");
0056     //
0057     minCSCHit = muonHistoParameters.getParameter<double>("minCSCHit");
0058     maxCSCHit = muonHistoParameters.getParameter<double>("maxCSCHit");
0059     nintCSCHit = muonHistoParameters.getParameter<int>("nintCSCHit");
0060     //
0061     minRPCHit = muonHistoParameters.getParameter<double>("minRPCHit");
0062     maxRPCHit = muonHistoParameters.getParameter<double>("maxRPCHit");
0063     nintRPCHit = muonHistoParameters.getParameter<int>("nintRPCHit");
0064     //
0065 
0066     minNTracks = muonHistoParameters.getParameter<int>("minNTracks");
0067     maxNTracks = muonHistoParameters.getParameter<int>("maxNTracks");
0068     nintNTracks = muonHistoParameters.getParameter<int>("nintNTracks");
0069     minFTracks = muonHistoParameters.getParameter<int>("minFTracks");
0070     maxFTracks = muonHistoParameters.getParameter<int>("maxFTracks");
0071     nintFTracks = muonHistoParameters.getParameter<int>("nintFTracks");
0072     minLayers = muonHistoParameters.getParameter<double>("minLayers");
0073     maxLayers = muonHistoParameters.getParameter<double>("maxLayers");
0074     nintLayers = muonHistoParameters.getParameter<int>("nintLayers");
0075     minPixels = muonHistoParameters.getParameter<double>("minPixels");
0076     maxPixels = muonHistoParameters.getParameter<double>("maxPixels");
0077     nintPixels = muonHistoParameters.getParameter<int>("nintPixels");
0078     minPhi = muonHistoParameters.getParameter<double>("minPhi");
0079     maxPhi = muonHistoParameters.getParameter<double>("maxPhi");
0080     nintPhi = muonHistoParameters.getParameter<int>("nintPhi");
0081     minDxy = muonHistoParameters.getParameter<double>("minDxy");
0082     maxDxy = muonHistoParameters.getParameter<double>("maxDxy");
0083     nintDxy = muonHistoParameters.getParameter<int>("nintDxy");
0084     minDz = muonHistoParameters.getParameter<double>("minDz");
0085     maxDz = muonHistoParameters.getParameter<double>("maxDz");
0086     nintDz = muonHistoParameters.getParameter<int>("nintDz");
0087     minRpos = muonHistoParameters.getParameter<double>("minRpos");
0088     maxRpos = muonHistoParameters.getParameter<double>("maxRpos");
0089     nintRpos = muonHistoParameters.getParameter<int>("nintRpos");
0090     minZpos = muonHistoParameters.getParameter<double>("minZpos");
0091     maxZpos = muonHistoParameters.getParameter<double>("maxZpos");
0092     nintZpos = muonHistoParameters.getParameter<int>("nintZpos");
0093     minPU = muonHistoParameters.getParameter<double>("minPU");
0094     maxPU = muonHistoParameters.getParameter<double>("maxPU");
0095     nintPU = muonHistoParameters.getParameter<int>("nintPU");
0096     //
0097     ptRes_rangeMin = muonHistoParameters.getParameter<double>("ptRes_rangeMin");
0098     ptRes_rangeMax = muonHistoParameters.getParameter<double>("ptRes_rangeMax");
0099     ptRes_nbin = muonHistoParameters.getParameter<int>("ptRes_nbin");
0100     etaRes_rangeMin = muonHistoParameters.getParameter<double>("etaRes_rangeMin");
0101     etaRes_rangeMax = muonHistoParameters.getParameter<double>("etaRes_rangeMax");
0102     etaRes_nbin = muonHistoParameters.getParameter<int>("etaRes_nbin");
0103     phiRes_rangeMin = muonHistoParameters.getParameter<double>("phiRes_rangeMin");
0104     phiRes_rangeMax = muonHistoParameters.getParameter<double>("phiRes_rangeMax");
0105     phiRes_nbin = muonHistoParameters.getParameter<int>("phiRes_nbin");
0106     cotThetaRes_rangeMin = muonHistoParameters.getParameter<double>("cotThetaRes_rangeMin");
0107     cotThetaRes_rangeMax = muonHistoParameters.getParameter<double>("cotThetaRes_rangeMax");
0108     cotThetaRes_nbin = muonHistoParameters.getParameter<int>("cotThetaRes_nbin");
0109     dxyRes_rangeMin = muonHistoParameters.getParameter<double>("dxyRes_rangeMin");
0110     dxyRes_rangeMax = muonHistoParameters.getParameter<double>("dxyRes_rangeMax");
0111     dxyRes_nbin = muonHistoParameters.getParameter<int>("dxyRes_nbin");
0112     dzRes_rangeMin = muonHistoParameters.getParameter<double>("dzRes_rangeMin");
0113     dzRes_rangeMax = muonHistoParameters.getParameter<double>("dzRes_rangeMax");
0114     dzRes_nbin = muonHistoParameters.getParameter<int>("dzRes_nbin");
0115     //
0116     usetracker = muonHistoParameters.getParameter<bool>("usetracker");
0117     usemuon = muonHistoParameters.getParameter<bool>("usemuon");
0118     do_TRKhitsPlots = muonHistoParameters.getParameter<bool>("do_TRKhitsPlots");
0119     do_MUOhitsPlots = muonHistoParameters.getParameter<bool>("do_MUOhitsPlots");
0120 
0121     if (useLogPt) {
0122       minPt = log10(std::max(0.01, minPt));
0123       maxPt = log10(maxPt);
0124     }
0125   }
0126 
0127   int minNTracks, maxNTracks, nintNTracks;
0128   int minFTracks, maxFTracks, nintFTracks;
0129   double mindR, maxdR;
0130   int nintdR;
0131   double minEta, maxEta;
0132   int nintEta;
0133   double minPt, maxPt;
0134   int nintPt;
0135   bool useLogPt;
0136   double minNHit, maxNHit;
0137   int nintNHit;
0138   double minDTHit, maxDTHit;
0139   int nintDTHit;
0140   double minCSCHit, maxCSCHit;
0141   int nintCSCHit;
0142   double minRPCHit, maxRPCHit;
0143   int nintRPCHit;
0144   double minLayers, maxLayers;
0145   int nintLayers;
0146   double minPixels, maxPixels;
0147   int nintPixels;
0148   double minPhi, maxPhi;
0149   int nintPhi;
0150   double minDxy, maxDxy;
0151   int nintDxy;
0152   double minDz, maxDz;
0153   int nintDz;
0154   double minRpos, maxRpos;
0155   int nintRpos;
0156   double minZpos, maxZpos;
0157   int nintZpos;
0158   double minPU, maxPU;
0159   int nintPU;
0160   //
0161   double ptRes_rangeMin, ptRes_rangeMax;
0162   int ptRes_nbin;
0163   double etaRes_rangeMin, etaRes_rangeMax;
0164   int etaRes_nbin;
0165   double phiRes_rangeMin, phiRes_rangeMax;
0166   int phiRes_nbin;
0167   double cotThetaRes_rangeMin, cotThetaRes_rangeMax;
0168   int cotThetaRes_nbin;
0169   double dxyRes_rangeMin, dxyRes_rangeMax;
0170   int dxyRes_nbin;
0171   double dzRes_rangeMin, dzRes_rangeMax;
0172   int dzRes_nbin;
0173   bool usetracker, usemuon;
0174   bool do_TRKhitsPlots, do_MUOhitsPlots;
0175 };
0176 
0177 class MuonTrackValidatorBase {
0178   typedef dqm::legacy::DQMStore DQMStore;
0179   typedef dqm::reco::MonitorElement MonitorElement;
0180 
0181 public:
0182   /// Constructor
0183   MuonTrackValidatorBase(const edm::ParameterSet& pset, edm::ConsumesCollector iC) : MuonTrackValidatorBase(pset) {
0184     bsSrc_Token = iC.consumes<reco::BeamSpot>(bsSrc);
0185 
0186     if (label_tp_refvector)
0187       tp_refvector_Token = iC.consumes<TrackingParticleRefVector>(label_tp);
0188     else
0189       tp_Token = iC.consumes<TrackingParticleCollection>(label_tp);
0190 
0191     pileupinfo_Token = iC.consumes<std::vector<PileupSummaryInfo>>(label_pileupinfo);
0192 
0193     histoParameters.reserve(label.size());
0194     for (unsigned int www = 0; www < label.size(); ++www) {
0195       track_Collection_Token[www] = iC.consumes<edm::View<reco::Track>>(label[www]);
0196       // Get histo parameters
0197       histoParameters.push_back(HistoParams(muonHistoParameters[www]));
0198     }
0199   }
0200 
0201   MuonTrackValidatorBase(const edm::ParameterSet& pset)
0202       : label(pset.getParameter<std::vector<edm::InputTag>>("label")),
0203         bsSrc(pset.getParameter<edm::InputTag>("beamSpot")),
0204         label_tp(pset.getParameter<edm::InputTag>("label_tp")),
0205         label_tp_refvector(pset.getParameter<bool>("label_tp_refvector")),
0206         label_pileupinfo(pset.getParameter<edm::InputTag>("label_pileupinfo")),
0207         associators(pset.getParameter<std::vector<std::string>>("associators")),
0208         out(pset.getParameter<std::string>("outputFile")),
0209         parametersDefiner(pset.getParameter<std::string>("parametersDefiner")),
0210         ignoremissingtkcollection_(pset.getUntrackedParameter<bool>("ignoremissingtrackcollection", false)),
0211         muonHistoParameters(pset.getParameter<std::vector<edm::ParameterSet>>("muonHistoParameters")) {
0212     histoParameters.reserve(label.size());
0213     for (size_t www = 0; www < label.size(); ++www) {
0214       histoParameters.emplace_back(muonHistoParameters[www]);
0215     }
0216   }
0217   /// Destructor
0218   virtual ~MuonTrackValidatorBase() noexcept(false) {}
0219 
0220   template <typename T>
0221   void fillPlotNoFlow(MonitorElement* h, T val) {
0222     h->Fill(std::min(std::max(val, ((T)h->getTH1()->GetXaxis()->GetXmin())), ((T)h->getTH1()->GetXaxis()->GetXmax())));
0223   }
0224 
0225   void doProfileX(TH2* th2, MonitorElement* me) {
0226     if (th2->GetNbinsX() == me->getNbinsX()) {
0227       TProfile* p1 = (TProfile*)th2->ProfileX();
0228       p1->Copy(*me->getTProfile());
0229       delete p1;
0230     } else {
0231       throw cms::Exception("MuonTrackValidator") << "Different number of bins!";
0232     }
0233   }
0234 
0235   void doProfileX(MonitorElement* th2m, MonitorElement* me) { doProfileX(th2m->getTH2F(), me); }
0236 
0237   //  virtual double getEta(double eta) {
0238   double getEta(double eta) { return eta; }
0239 
0240   //  virtual double getPt(double pt) {
0241   double getPt(double pt) { return pt; }
0242 
0243   void BinLogX(TH1* h) {
0244     TAxis* axis = h->GetXaxis();
0245     int bins = axis->GetNbins();
0246 
0247     float from = axis->GetXmin();
0248     float to = axis->GetXmax();
0249     float width = (to - from) / bins;
0250     float* new_bins = new float[bins + 1];
0251 
0252     for (int i = 0; i <= bins; i++) {
0253       new_bins[i] = TMath::Power(10, from + i * width);
0254     }
0255     axis->Set(bins, new_bins);
0256     delete[] new_bins;
0257   }
0258 
0259 protected:
0260   std::vector<edm::InputTag> label;
0261   edm::InputTag bsSrc;
0262   edm::InputTag label_tp;
0263   bool label_tp_refvector;
0264   edm::InputTag label_pileupinfo;
0265   std::vector<std::string> associators;
0266   std::string out;
0267   std::string parametersDefiner;
0268   std::vector<edm::EDGetTokenT<edm::View<reco::Track>>> track_Collection_Token;
0269   edm::EDGetTokenT<reco::BeamSpot> bsSrc_Token;
0270   edm::EDGetTokenT<TrackingParticleCollection> tp_Token;
0271   edm::EDGetTokenT<TrackingParticleRefVector> tp_refvector_Token;
0272   edm::EDGetTokenT<std::vector<PileupSummaryInfo>> pileupinfo_Token;
0273   edm::ESHandle<MagneticField> theMF;
0274 
0275   bool ignoremissingtkcollection_;
0276 
0277   std::vector<edm::ParameterSet> muonHistoParameters;
0278 
0279   std::vector<HistoParams> histoParameters;
0280 
0281   // collections
0282   std::vector<MonitorElement*> h_assoc_coll, h_simul_coll, h_reco_coll, h_assoc2_coll;
0283   //1D
0284   std::vector<MonitorElement*> h_tracks, h_fakes, h_nhits, h_charge;
0285   std::vector<MonitorElement*> h_recoeta, h_assoceta, h_assoc2eta, h_simuleta, h_misideta;
0286   std::vector<MonitorElement*> h_recopT, h_assocpT, h_assoc2pT, h_simulpT, h_misidpT;
0287   std::vector<MonitorElement*> h_recohit, h_assochit, h_assoc2hit, h_simulhit, h_misidhit;
0288   std::vector<MonitorElement*> h_recophi, h_assocphi, h_assoc2phi, h_simulphi, h_misidphi;
0289   std::vector<MonitorElement*> h_recodxy, h_assocdxy, h_assoc2dxy, h_simuldxy, h_misiddxy;
0290   std::vector<MonitorElement*> h_recodz, h_assocdz, h_assoc2dz, h_simuldz, h_misiddz;
0291   std::vector<MonitorElement*> h_recopu, h_assocpu, h_assoc2pu, h_simulpu, h_misidpu;
0292   std::vector<MonitorElement*> h_recodR, h_assocdR, h_assoc2dR, h_simuldR, h_misiddR;
0293 
0294   std::vector<MonitorElement*> h_assocpTB, h_simulpTB;
0295   std::vector<MonitorElement*> h_assocphiB, h_simulphiB;
0296 
0297   std::vector<MonitorElement*> h_assocpTO, h_simulpTO;
0298   std::vector<MonitorElement*> h_assocphiO, h_simulphiO;
0299 
0300   std::vector<MonitorElement*> h_assocpTE, h_simulpTE;
0301   std::vector<MonitorElement*> h_assocphiE, h_simulphiE;
0302 
0303   std::vector<MonitorElement*> h_assocRpos, h_simulRpos, h_assocZpos, h_simulZpos;
0304   std::vector<MonitorElement*> h_etaRes;
0305 
0306   std::vector<MonitorElement*> h_nchi2, h_nchi2_prob, h_losthits;
0307   std::vector<MonitorElement*> h_nmisslayers_inner, h_nmisslayers_outer, h_nlosthits;
0308   std::vector<MonitorElement*> h_assochi2, h_assochi2_prob;
0309   std::vector<MonitorElement*> h_assocFraction, h_assocSharedHit;
0310 
0311   //2D
0312   std::vector<MonitorElement*> nRecHits_vs_nSimHits;
0313   std::vector<MonitorElement*> h_PurityVsQuality;
0314   std::vector<MonitorElement*> chi2_vs_nhits, etares_vs_eta;
0315   std::vector<MonitorElement*> ptres_vs_phi, chi2_vs_phi, nhits_vs_phi, phires_vs_phi;
0316 
0317   std::vector<MonitorElement*> nhits_vs_eta, nDThits_vs_eta, nCSChits_vs_eta, nRPChits_vs_eta, nGEMhits_vs_eta,
0318       nME0hits_vs_eta;
0319   std::vector<MonitorElement*> chi2_vs_eta, nlosthits_vs_eta;
0320   std::vector<MonitorElement*> nTRK_LayersWithMeas_vs_eta, nPixel_LayersWithMeas_vs_eta;
0321 
0322   std::vector<MonitorElement*> dxyres_vs_eta, ptres_vs_eta, dzres_vs_eta, phires_vs_eta, thetaCotres_vs_eta;
0323   std::vector<MonitorElement*> dxyres_vs_pt, ptres_vs_pt, dzres_vs_pt, phires_vs_pt, thetaCotres_vs_pt;
0324 
0325   std::vector<MonitorElement*> dxypull_vs_eta, ptpull_vs_eta, dzpull_vs_eta, phipull_vs_eta, thetapull_vs_eta;
0326   std::vector<MonitorElement*> ptpull_vs_phi, phipull_vs_phi, thetapull_vs_phi;
0327   std::vector<MonitorElement*> h_dxypulleta, h_ptpulleta, h_dzpulleta, h_phipulleta, h_thetapulleta;
0328   std::vector<MonitorElement*> h_ptpullphi, h_phipullphi, h_thetapullphi;
0329   std::vector<MonitorElement*> h_ptpull, h_qoverppull, h_thetapull, h_phipull, h_dxypull, h_dzpull;
0330 };
0331 
0332 #endif