PatExamples/plugins/PatBJetVertexAnalyzer.cc

0001 #include <algorithm>
0002 #include <iostream>
0003 #include <cmath>
0004 #include <vector>
0005 #include <string>
0006
0007 #include <TH1.h>
0008 #include <TProfile.h>
0009
0010 #include <Math/VectorUtil.h>
0011 #include <Math/GenVector/PxPyPzE4D.h>
0012 #include <Math/GenVector/PxPyPzM4D.h>
0013
0014 #include "FWCore/Framework/interface/Frameworkfwd.h"
0015 #include "FWCore/Framework/interface/Event.h"
0016 #include "FWCore/Framework/interface/EventSetup.h"
0017 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0018 #include "FWCore/Utilities/interface/InputTag.h"
0019 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0020 #include "FWCore/ServiceRegistry/interface/Service.h"
0021
0022 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0023
0024 #include "DataFormats/BTauReco/interface/SecondaryVertexTagInfo.h"
0025 #include "DataFormats/Math/interface/LorentzVector.h"
0026 #include "DataFormats/VertexReco/interface/Vertex.h"
0027 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0028 #include "DataFormats/TrackReco/interface/Track.h"
0029 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0030 #include "DataFormats/PatCandidates/interface/Jet.h"
0031 #include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
0032
0033 class PatBJetVertexAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0034 public:
0035   /// constructor and destructor
0036   PatBJetVertexAnalyzer(const edm::ParameterSet &params);
0037   ~PatBJetVertexAnalyzer() override;
0038
0039   // virtual methods called from base class EDAnalyzer
0040   void beginJob() override;
0041   void analyze(const edm::Event &event, const edm::EventSetup &es) override;
0042
0043 private:
0044   // configuration parameters
0045   edm::EDGetTokenT<pat::JetCollection> jetsToken_;
0046
0047   double jetPtCut_;   // minimum (uncorrected) jet energy
0048   double jetEtaCut_;  // maximum |eta| for jet
0049   double maxDeltaR_;  // angle between jet and tracks
0050
0051   // jet flavour constants
0052
0053   enum Flavour { ALL_JETS = 0, UDSG_JETS, C_JETS, B_JETS, NONID_JETS, N_JET_TYPES };
0054
0055   TH1 *flavours_;
0056
0057   // one group of plots per jet flavour;
0058   struct Plots {
0059     TH1 *nVertices;
0060     TH1 *deltaR, *mass, *dist, *distErr, *distSig;
0061     TH1 *nTracks, *chi2;
0062   } plots_[N_JET_TYPES];
0063 };
0064
0065 PatBJetVertexAnalyzer::PatBJetVertexAnalyzer(const edm::ParameterSet &params)
0066     : jetsToken_(consumes<pat::JetCollection>(params.getParameter<edm::InputTag>("jets"))),
0067       jetPtCut_(params.getParameter<double>("jetPtCut")),
0068       jetEtaCut_(params.getParameter<double>("jetEtaCut")) {
0069   usesResource(TFileService::kSharedResource);
0070 }
0071
0072 PatBJetVertexAnalyzer::~PatBJetVertexAnalyzer() {}
0073
0074 void PatBJetVertexAnalyzer::beginJob() {
0075   // retrieve handle to auxiliary service
0076   //  used for storing histograms into ROOT file
0077   edm::Service<TFileService> fs;
0078
0079   flavours_ = fs->make<TH1F>("flavours", "jet flavours", 5, 0, 5);
0080
0081   // book histograms for all jet flavours
0082   for (unsigned int i = 0; i < N_JET_TYPES; i++) {
0083     Plots &plots = plots_[i];
0084     const char *flavour, *name;
0085
0086     switch ((Flavour)i) {
0087       case ALL_JETS:
0088         flavour = "all jets";
0089         name = "all";
0090         break;
0091       case UDSG_JETS:
0092         flavour = "light flavour jets";
0093         name = "udsg";
0094         break;
0095       case C_JETS:
0096         flavour = "charm jets";
0097         name = "c";
0098         break;
0099       case B_JETS:
0100         flavour = "bottom jets";
0101         name = "b";
0102         break;
0103       default:
0104         flavour = "unidentified jets";
0105         name = "ni";
0106         break;
0107     }
0108
0109     plots.nVertices =
0110         fs->make<TH1F>(Form("nVertices_%s", name), Form("number of secondary vertices in %s", flavour), 5, 0, 5);
0111     plots.deltaR = fs->make<TH1F>(Form("deltaR_%s", name),
0112                                   Form("\\DeltaR between vertex direction and jet direction in %s", flavour),
0113                                   100,
0114                                   0,
0115                                   0.5);
0116     plots.mass = fs->make<TH1F>(Form("mass_%s", name), Form("vertex mass in %s", flavour), 100, 0, 10);
0117     plots.dist =
0118         fs->make<TH1F>(Form("dist_%s", name), Form("Transverse distance between PV and SV in %s", flavour), 100, 0, 2);
0119     plots.distErr = fs->make<TH1F>(
0120         Form("distErr_%s", name), Form("Transverse distance error between PV and SV in %s", flavour), 100, 0, 0.5);
0121     plots.distSig = fs->make<TH1F>(
0122         Form("distSig_%s", name), Form("Transverse distance significance between PV and SV in %s", flavour), 100, 0, 50);
0123     plots.nTracks = fs->make<TH1F>(
0124         Form("nTracks_%s", name), Form("number of tracks at secondary vertex in %s", flavour), 20, 0, 20);
0125     plots.chi2 =
0126         fs->make<TH1F>(Form("chi2_%s", name), Form("secondary vertex fit \\chi^{2} in %s", flavour), 100, 0, 50);
0127   }
0128 }
0129
0130 // helper function to sort the tracks by impact parameter significance
0131
0132 void PatBJetVertexAnalyzer::analyze(const edm::Event &event, const edm::EventSetup &es) {
0133   // handle to the jets collection
0134   edm::Handle<pat::JetCollection> jetsHandle;
0135   event.getByToken(jetsToken_, jetsHandle);
0136
0137   // now go through all jets
0138   for (pat::JetCollection::const_iterator jet = jetsHandle->begin(); jet != jetsHandle->end(); ++jet) {
0139     // only look at jets that pass the pt and eta cut
0140     if (jet->pt() < jetPtCut_ || std::abs(jet->eta()) > jetEtaCut_)
0141       continue;
0142
0143     Flavour flavour;
0144     // find out the jet flavour (differs between quark and anti-quark)
0145     switch (std::abs(jet->partonFlavour())) {
0146       case 1:
0147       case 2:
0148       case 3:
0149       case 21:
0150         flavour = UDSG_JETS;
0151         break;
0152       case 4:
0153         flavour = C_JETS;
0154         break;
0155       case 5:
0156         flavour = B_JETS;
0157         break;
0158       default:
0159         flavour = NONID_JETS;
0160     }
0161
0162     // simply count the number of accepted jets
0163     flavours_->Fill(ALL_JETS);
0164     flavours_->Fill(flavour);
0165
0166     // retrieve the "secondary vertex tag infos"
0167     // this is the output of the b-tagging reconstruction code
0168     // and contains secondary vertices in the jets
0169     const reco::SecondaryVertexTagInfo &svTagInfo = *jet->tagInfoSecondaryVertex();
0170
0171     // count the number of secondary vertices
0172     plots_[ALL_JETS].nVertices->Fill(svTagInfo.nVertices());
0173     plots_[flavour].nVertices->Fill(svTagInfo.nVertices());
0174
0175     // ignore jets without SV from now on
0176     if (svTagInfo.nVertices() < 1)
0177       continue;
0178
0179     // pick the first secondary vertex (the "best" one)
0180     const reco::Vertex &sv = svTagInfo.secondaryVertex(0);
0181
0182     // and plot number of tracks and chi^2
0183     plots_[ALL_JETS].nTracks->Fill(sv.tracksSize());
0184     plots_[flavour].nTracks->Fill(sv.tracksSize());
0185
0186     plots_[ALL_JETS].chi2->Fill(sv.chi2());
0187     plots_[flavour].chi2->Fill(sv.chi2());
0188
0189     // the precomputed transverse distance to the primary vertex
0190     Measurement1D distance = svTagInfo.flightDistance(0, true);
0191
0192     plots_[ALL_JETS].dist->Fill(distance.value());
0193     plots_[flavour].dist->Fill(distance.value());
0194
0195     plots_[ALL_JETS].distErr->Fill(distance.error());
0196     plots_[flavour].distErr->Fill(distance.error());
0197
0198     plots_[ALL_JETS].distSig->Fill(distance.significance());
0199     plots_[flavour].distSig->Fill(distance.significance());
0200
0201     // the precomputed direction with respect to the primary vertex
0202     const GlobalVector &dir = svTagInfo.flightDirection(0);
0203
0204     // unfortunately CMSSW hsa all kinds of vectors,
0205     // and sometimes we need to convert them *sigh*
0206     math::XYZVector dir2(dir.x(), dir.y(), dir.z());
0207
0208     // compute a few variables that we are plotting
0209     double deltaR = ROOT::Math::VectorUtil::DeltaR(jet->momentum(), dir2);
0210
0211     plots_[ALL_JETS].deltaR->Fill(deltaR);
0212     plots_[flavour].deltaR->Fill(deltaR);
0213
0214     // compute the invariant mass from a four-vector sum
0215     math::XYZTLorentzVector trackFourVectorSum;
0216
0217     // loop over all tracks in the vertex
0218     for (reco::Vertex::trackRef_iterator track = sv.tracks_begin(); track != sv.tracks_end(); ++track) {
0219       ROOT::Math::LorentzVector<ROOT::Math::PxPyPzM4D<double> > vec;
0220       vec.SetPx((*track)->px());
0221       vec.SetPy((*track)->py());
0222       vec.SetPz((*track)->pz());
0223       vec.SetM(0.13957);  // pion mass
0224       trackFourVectorSum += vec;
0225     }
0226
0227     // get the invariant mass: sqrt(E² - px² - py² - pz²)
0228     double vertexMass = trackFourVectorSum.M();
0229
0230     plots_[ALL_JETS].mass->Fill(vertexMass);
0231     plots_[flavour].mass->Fill(vertexMass);
0232   }
0233 }
0234
0235 #include "FWCore/Framework/interface/MakerMacros.h"
0236
0237 DEFINE_FWK_MODULE(PatBJetVertexAnalyzer);