Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-03-17 11:28:33

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