Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:10:56

0001 #include "DQM/Physics/src/FSQDQM.h"
0002 #include <memory>
0003 
0004 // DQM
0005 #include "DQMServices/Core/interface/DQMStore.h"
0006 
0007 // Framework
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include "FWCore/Framework/interface/EventSetup.h"
0010 #include "FWCore/Framework/interface/LuminosityBlock.h"
0011 #include "FWCore/ServiceRegistry/interface/Service.h"
0012 #include "FWCore/ParameterSet/interface/FileInPath.h"
0013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0014 #include "FWCore/Framework/interface/Run.h"
0015 #include "DataFormats/Provenance/interface/EventID.h"
0016 
0017 // Candidate handling
0018 #include "DataFormats/Candidate/interface/Candidate.h"
0019 #include "DataFormats/Candidate/interface/CandidateFwd.h"
0020 #include "DataFormats/Candidate/interface/OverlapChecker.h"
0021 #include "DataFormats/Candidate/interface/CompositeCandidate.h"
0022 #include "DataFormats/Candidate/interface/CompositeCandidateFwd.h"
0023 #include "DataFormats/Candidate/interface/CandMatchMap.h"
0024 
0025 // Vertex utilities
0026 #include "DataFormats/VertexReco/interface/Vertex.h"
0027 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0028 
0029 // Other
0030 #include "DataFormats/TrackReco/interface/Track.h"
0031 #include "DataFormats/DetId/interface/DetId.h"
0032 #include "DataFormats/Common/interface/RefToBase.h"
0033 
0034 // Math
0035 #include "DataFormats/Math/interface/Vector3D.h"
0036 #include "DataFormats/Math/interface/LorentzVector.h"
0037 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0038 #include "DataFormats/Common/interface/AssociationVector.h"
0039 #include "DataFormats/Math/interface/Point3D.h"
0040 #include "DataFormats/Math/interface/deltaR.h"
0041 #include "DataFormats/Math/interface/deltaPhi.h"
0042 
0043 // vertexing
0044 
0045 // Transient tracks
0046 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0047 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0048 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0049 #include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h"
0050 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0051 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0052 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0053 
0054 // JetCorrection
0055 #include "JetMETCorrections/Objects/interface/JetCorrector.h"
0056 
0057 // Substructure
0058 #include "RecoJets/JetAlgorithms/interface/CATopJetHelper.h"
0059 #include "DataFormats/BTauReco/interface/CATopJetTagInfo.h"
0060 
0061 // ROOT
0062 #include "TLorentzVector.h"
0063 
0064 // STDLIB
0065 #include <iostream>
0066 #include <iomanip>
0067 #include <cstdio>
0068 #include <string>
0069 #include <sstream>
0070 #include <cmath>
0071 
0072 using namespace edm;
0073 using namespace std;
0074 using namespace reco;
0075 using namespace trigger;
0076 
0077 struct SortByPt
0078 
0079 {
0080   bool operator()(const TLorentzVector& a, const TLorentzVector& b) const { return a.Pt() > b.Pt(); }
0081 };
0082 
0083 FSQDQM::FSQDQM(const edm::ParameterSet& iConfig) {
0084   edm::LogInfo("FSQDQM") << " Creating FSQDQM "
0085                          << "\n";
0086   pvs_ = consumes<edm::View<reco::Vertex> >(iConfig.getParameter<edm::InputTag>("pvs"));
0087   labelPFJet_ = iConfig.getParameter<std::string>("LabelPFJet");
0088   labelCastorJet_ = iConfig.getParameter<std::string>("LabelCastorJet");
0089   labelTrack_ = iConfig.getParameter<std::string>("LabelTrack");
0090   tok_pfjet_ = consumes<reco::PFJetCollection>(labelPFJet_);
0091   tok_track_ = consumes<reco::TrackCollection>(labelTrack_);
0092   tok_castorjet_ = consumes<reco::BasicJetCollection>(labelCastorJet_);
0093 }
0094 
0095 FSQDQM::~FSQDQM() {
0096   edm::LogInfo("FSQDQM") << " Deleting FSQDQM "
0097                          << "\n";
0098   // do anything here that needs to be done at desctruction time
0099   // (e.g. close files, deallocate resources etc.)
0100 }
0101 void FSQDQM::bookHistograms(DQMStore::IBooker& bei, edm::Run const&, edm::EventSetup const&) {
0102   bei.setCurrentFolder("Physics/FSQ");
0103   PFJetpt = bei.book1D("PFJetpt", ";p_{T}(PFJet)", 100, 0.0, 100);
0104   PFJeteta = bei.book1D("PFJeteta", ";#eta(PFJet)", 50, -2.5, 2.5);
0105   PFJetphi = bei.book1D("PFJetphi", ";#phi(PFJet)", 50, -3.14, 3.14);
0106   PFJetMulti = bei.book1D("PFJetMulti", ";No. of PFJets", 10, -0.5, 9.5);
0107   PFJetRapidity = bei.book1D("PFJetRapidity", ";PFJetRapidity", 50, -6.0, 6.0);
0108   CastorJetphi = bei.book1D("CastorJetphi", ";#phi(CastorJet)", 50, -3.14, 3.14);
0109   CastorJetMulti = bei.book1D("CastorJetMulti", ";No. of CastorJets", 10, -0.5, 9.5);
0110   Track_HP_Phi = bei.book1D("Track_HP_Phi", ";#phi(HPtrack)", 50, -3.14, 3.14);
0111   Track_HP_Eta = bei.book1D("Track_HP_Eta", ";#eta(HPtrack)", 50, -2.5, 2.5);
0112   Track_HP_Pt = bei.book1D("Track_HP_Pt", ";p_{T}(HPtrack)", 500, 0.0, 50);
0113   Track_HP_ptErr_over_pt = bei.book1D("Track_HP_ptErr_over_pt", ";{p_{T}Err}/p_{T}", 100, 0, 0.1);
0114   Track_HP_dzvtx_over_dzerr = bei.book1D("Track_HP_dzvtx_over_dzerr", ";dZerr/dZ", 100, -10, 10);
0115   Track_HP_dxyvtx_over_dxyerror = bei.book1D("Track_HP_dxyvtx_over_dxyerror", ";dxyErr/dxy", 100, -10, 10);
0116   NPV = bei.book1D("NPV", ";NPV", 10, -0.5, 9.5);
0117   PV_chi2 = bei.book1D("PV_chi2", ";PV_chi2", 100, 0.0, 2.0);
0118   PV_d0 = bei.book1D("PV_d0", ";PV_d0", 100, -10.0, 10.0);
0119   PV_numTrks = bei.book1D("PV_numTrks", ";PV_numTrks", 100, -0.5, 99.5);
0120   PV_sumTrks = bei.book1D("PV_sumTrks", ";PV_sumTrks", 100, 0, 100);
0121   h_ptsum_towards = bei.book1D("h_ptsum_towards", ";h_ptsum_towards", 100, 0, 100);
0122   h_ptsum_transverse = bei.book1D("h_ptsum_transverse", ";h_ptsum_transverse", 100, 0, 100);
0123 
0124   h_ntracks = bei.book1D("h_ntracks", ";h_ntracks", 50, -0.5, 49.5);
0125   h_trkptsum = bei.book1D("h_trkptsum", ";h_trkptsum", 100, 0, 100);
0126   h_ptsum_away = bei.book1D("h_ptsum_away", ";h_ptsum_away", 100, 0, 100);
0127   h_ntracks_towards = bei.book1D("h_ntracks_towards", ";h_ntracks_towards", 50, -0.5, 49.5);
0128   h_ntracks_transverse = bei.book1D("h_ntracks_transverse", ";h_ntracks_transverse", 50, -0.5, 49.5);
0129   h_ntracks_away = bei.book1D("h_ntracks_away", ";h_ntracks_away", 50, -0.5, 49.5);
0130 
0131   h_leadingtrkpt_ntrk_away =
0132       bei.bookProfile("h_leadingtrkpt_ntrk_away", "h_leadingtrkpt_ntrk_away", 50, 0.0, 50, 0, 30, " ");
0133   h_leadingtrkpt_ntrk_towards =
0134       bei.bookProfile("h_leadingtrkpt_ntrk_towards", "h_leadingtrkpt_ntrk_towards", 50, 0.0, 50, 0, 30, " ");
0135   h_leadingtrkpt_ntrk_transverse =
0136       bei.bookProfile("h_leadingtrkpt_ntrk_transverse", "h_leadingtrkpt_ntrk_transverse", 50, 0.0, 50, 0, 30, " ");
0137   h_leadingtrkpt_ptsum_away =
0138       bei.bookProfile("h_leadingtrkpt_ptsum_away", "h_leadingtrkpt_ptsum_away", 50, 0.0, 50, 0, 30, " ");
0139   h_leadingtrkpt_ptsum_towards =
0140       bei.bookProfile("h_leadingtrkpt_ptsum_towards", "h_leadingtrkpt_ptsum_towards", 50, 0.0, 50, 0, 30, " ");
0141   h_leadingtrkpt_ptsum_transverse =
0142       bei.bookProfile("h_leadingtrkpt_ptsum_transverse", "h_leadingtrkpt_ptsum_transverse", 50, 0.0, 50, 0, 30, " ");
0143 }
0144 
0145 // ------------ method called for each event  ------------
0146 void FSQDQM::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0147   using namespace edm;
0148   using namespace std;
0149   using namespace reco;
0150 
0151   runNumber_ = iEvent.id().run();
0152   eventNumber_ = iEvent.id().event();
0153   lumiNumber_ = iEvent.id().luminosityBlock();
0154   bxNumber_ = iEvent.bunchCrossing();
0155 
0156   edm::Handle<edm::View<reco::Vertex> > privtxs;
0157   //    if(!iEvent.getByToken(pvs_, privtxs)) return;
0158   iEvent.getByToken(pvs_, privtxs);
0159   double bestvz = -999.9, bestvx = -999.9, bestvy = -999.9;
0160   double bestvzError = -999.9, bestvxError = -999.9, bestvyError = -999.9;
0161   if (privtxs.isValid()) {
0162     const reco::Vertex& pvtx = privtxs->front();
0163     NPV->Fill(privtxs->size());
0164     if (privtxs->begin() != privtxs->end() && !(pvtx.isFake()) && pvtx.position().Rho() <= 2. &&
0165         fabs(pvtx.position().z()) <= 24) {
0166       bestvz = pvtx.z();
0167       bestvx = pvtx.x();
0168       bestvy = pvtx.y();
0169       bestvzError = pvtx.zError();
0170       bestvxError = pvtx.xError();
0171       bestvyError = pvtx.yError();
0172       PV_chi2->Fill(pvtx.normalizedChi2());
0173       PV_d0->Fill(sqrt(pvtx.x() * pvtx.x() + pvtx.y() * pvtx.y()));
0174       PV_numTrks->Fill(pvtx.tracksSize());
0175       double vertex_sumTrks = 0.0;
0176       for (reco::Vertex::trackRef_iterator iTrack = pvtx.tracks_begin(); iTrack != pvtx.tracks_end(); iTrack++) {
0177         vertex_sumTrks += (*iTrack)->pt();
0178       }
0179       PV_sumTrks->Fill(vertex_sumTrks);
0180     }
0181   }
0182 
0183   std::vector<Jet> recoPFJets;
0184   recoPFJets.clear();
0185   int nPFCHSJet = 0;
0186   edm::Handle<PFJetCollection> pfjetchscoll;
0187   //     if(!iEvent.getByToken(tok_pfjet_, pfjetchscoll)) return;
0188   iEvent.getByToken(tok_pfjet_, pfjetchscoll);
0189   if (pfjetchscoll.isValid()) {
0190     const reco::PFJetCollection* pfchsjets = pfjetchscoll.product();
0191     reco::PFJetCollection::const_iterator pfjetchsclus = pfchsjets->begin();
0192     for (pfjetchsclus = pfchsjets->begin(); pfjetchsclus != pfchsjets->end(); ++pfjetchsclus) {
0193       PFJetpt->Fill(pfjetchsclus->pt());
0194       PFJeteta->Fill(pfjetchsclus->eta());
0195       PFJetphi->Fill(pfjetchsclus->phi());
0196       PFJetRapidity->Fill(pfjetchsclus->rapidity());
0197       nPFCHSJet++;
0198     }
0199     PFJetMulti->Fill(nPFCHSJet);
0200   }
0201 
0202   std::vector<Jet> recoCastorJets;
0203   recoCastorJets.clear();
0204   edm::Handle<BasicJetCollection> castorJets;
0205   //     if(!iEvent.getByToken(tok_castorjet_, castorJets)) return;
0206   iEvent.getByToken(tok_castorjet_, castorJets);
0207   if (castorJets.isValid()) {
0208     for (unsigned ijet = 0; ijet < castorJets->size(); ijet++) {
0209       recoCastorJets.push_back((*castorJets)[ijet]);
0210     }
0211     for (unsigned ijet = 0; ijet < recoCastorJets.size(); ijet++) {
0212       //       cout<<recoCastorJets[ijet].pt()<<endl;
0213       CastorJetphi->Fill(recoCastorJets[ijet].phi());
0214 
0215       CastorJetMulti->Fill(recoCastorJets.size());
0216     }
0217   }
0218   edm::Handle<reco::TrackCollection> itracks;
0219   //     if(!iEvent.getByToken(tok_track_, itracks)) return;
0220   iEvent.getByToken(tok_track_, itracks);
0221   if (itracks.isValid()) {
0222     reco::TrackBase::TrackQuality hiPurity = reco::TrackBase::qualityByName("highPurity");
0223     std::vector<TLorentzVector> T_trackRec_P4;
0224 
0225     int ntracks = 0;
0226     int ntracks_towards = 0;
0227     int ntracks_transverse = 0;
0228     int ntracks_away = 0;
0229     float ptsum = 0;
0230     float dphi = 0;
0231     float ptsum_towards = 0;
0232     float ptsum_transverse = 0;
0233     float ptsum_away = 0;
0234 
0235     T_trackRec_P4.clear();
0236     for (reco::TrackCollection::const_iterator iT = itracks->begin(); iT != itracks->end(); ++iT) {
0237       if (iT->quality(hiPurity)) {
0238         math::XYZPoint bestvtx(bestvx, bestvy, bestvz);
0239         double dzvtx = iT->dz(bestvtx);
0240         double dxyvtx = iT->dxy(bestvtx);
0241         double dzerror = sqrt(iT->dzError() * iT->dzError() + bestvzError * bestvzError);
0242         double dxyerror = sqrt(iT->d0Error() * iT->d0Error() + bestvxError * bestvyError);
0243         if ((iT->ptError()) / iT->pt() < 0.05 && dzvtx < 3.0 && dxyvtx < 3.0) {
0244           TLorentzVector trk;
0245           trk.SetPtEtaPhiE(iT->pt(), iT->eta(), iT->phi(), iT->p());
0246           T_trackRec_P4.push_back(trk);
0247           Track_HP_Eta->Fill(iT->eta());
0248           Track_HP_Phi->Fill(iT->phi());
0249           Track_HP_Pt->Fill(iT->pt());
0250           Track_HP_ptErr_over_pt->Fill((iT->ptError()) / iT->pt());
0251           Track_HP_dzvtx_over_dzerr->Fill(dzvtx / dzerror);
0252           Track_HP_dxyvtx_over_dxyerror->Fill(dxyvtx / dxyerror);
0253         }
0254       }
0255     }
0256 
0257     float highest_pt_track = -999;
0258     int index = -999;
0259     for (unsigned int k = 0; k < T_trackRec_P4.size(); k++) {
0260       if (T_trackRec_P4.at(k).Pt() > highest_pt_track) {
0261         highest_pt_track = T_trackRec_P4.at(k).Pt();
0262         index = k;
0263       }
0264     }
0265     unsigned int finalid = abs(index);
0266     //   if(T_trackRec_P4.at(index).Pt()!=0.){
0267     if (finalid < T_trackRec_P4.size()) {
0268       //     std::sort(T_trackRec_P4.begin(), T_trackRec_P4.end(), SortByPt());
0269       for (unsigned int itrk = 0; itrk < T_trackRec_P4.size(); itrk++) {
0270         ++ntracks;
0271         ptsum = ptsum + T_trackRec_P4.at(itrk).Pt();
0272         dphi = deltaPhi(T_trackRec_P4.at(itrk).Phi(), T_trackRec_P4.at(index).Phi());
0273         if (fabs(dphi) < 1.05) {
0274           ++ntracks_towards;
0275           ptsum_towards = ptsum_towards + T_trackRec_P4.at(itrk).Pt();
0276         }
0277         if (fabs(dphi) > 1.05 && fabs(dphi) < 2.09) {
0278           ++ntracks_transverse;
0279           ptsum_transverse = ptsum_transverse + T_trackRec_P4.at(itrk).Pt();
0280         }
0281         if (fabs(dphi) > 2.09) {
0282           ++ntracks_away;
0283           ptsum_away = ptsum_away + T_trackRec_P4.at(itrk).Pt();
0284         }
0285       }
0286 
0287       h_ntracks->Fill(ntracks);
0288       h_trkptsum->Fill(ptsum);
0289       h_ptsum_towards->Fill(ptsum_towards);
0290       h_ptsum_transverse->Fill(ptsum_transverse);
0291       h_ptsum_away->Fill(ptsum_away);
0292       h_ntracks_towards->Fill(ntracks_towards);
0293       h_ntracks_transverse->Fill(ntracks_transverse);
0294       h_ntracks_away->Fill(ntracks_away);
0295 
0296       if (!T_trackRec_P4.empty()) {
0297         h_leadingtrkpt_ntrk_towards->Fill(T_trackRec_P4.at(index).Pt(), ntracks_towards / 8.37);
0298         h_leadingtrkpt_ntrk_transverse->Fill(T_trackRec_P4.at(index).Pt(), ntracks_transverse / 8.37);
0299         h_leadingtrkpt_ntrk_away->Fill(T_trackRec_P4.at(index).Pt(), ntracks_away / 8.37);
0300         h_leadingtrkpt_ptsum_towards->Fill(T_trackRec_P4.at(index).Pt(), ptsum_towards / 8.37);
0301         h_leadingtrkpt_ptsum_transverse->Fill(T_trackRec_P4.at(index).Pt(), ptsum_transverse / 8.37);
0302         h_leadingtrkpt_ptsum_away->Fill(T_trackRec_P4.at(index).Pt(), ptsum_away / 8.37);
0303       }
0304     }
0305   }
0306 }  //analyze
0307 
0308 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0309 
0310 //define this as a plug-in
0311 //DEFINE_FWK_MODULE(FSQDQM);
0312 
0313 //  LocalWords:  TH1F ptsum fs ntracks