Physics/src/FSQDQM.cc

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