PatExamples/plugins/PatBJetTrackAnalyzer.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
0012 #include "FWCore/Framework/interface/Frameworkfwd.h"
0013 #include "FWCore/Framework/interface/Event.h"
0014 #include "FWCore/Framework/interface/EventSetup.h"
0015 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0016 #include "FWCore/Utilities/interface/InputTag.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0018 #include "FWCore/ServiceRegistry/interface/Service.h"
0019
0020 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0021
0022 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0023 #include "DataFormats/VertexReco/interface/Vertex.h"
0024 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0025 #include "DataFormats/TrackReco/interface/Track.h"
0026 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0027 #include "DataFormats/PatCandidates/interface/Jet.h"
0028 #include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
0029
0030 class PatBJetTrackAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0031 public:
0032   /// constructor and destructor
0033   PatBJetTrackAnalyzer(const edm::ParameterSet &params);
0034   ~PatBJetTrackAnalyzer() override;
0035
0036   // virtual methods called from base class EDAnalyzer
0037   void beginJob() override;
0038   void analyze(const edm::Event &event, const edm::EventSetup &es) override;
0039
0040 private:
0041   // configuration parameters
0042   edm::EDGetTokenT<pat::JetCollection> jetsToken_;
0043   edm::EDGetTokenT<reco::TrackCollection> tracksToken_;
0044   edm::EDGetTokenT<reco::BeamSpot> beamSpotToken_;
0045   edm::EDGetTokenT<reco::VertexCollection> primaryVerticesToken_;
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   double minPt_;               // track pt quality cut
0052   unsigned int minPixelHits_;  // minimum number of pixel hits
0053   unsigned int minTotalHits_;  // minimum number of total hits
0054
0055   unsigned int nThTrack_;  // n-th hightest track to choose
0056
0057   // jet flavour constants
0058
0059   enum Flavour { ALL_JETS = 0, UDSG_JETS, C_JETS, B_JETS, NONID_JETS, N_JET_TYPES };
0060
0061   TH1 *flavours_;
0062
0063   // one group of plots per jet flavour;
0064   struct Plots {
0065     TH1 *allIP, *allIPErr, *allIPSig;
0066     TH1 *trackIP, *trackIPErr, *trackIPSig;
0067     TH1 *negativeIP, *negativeIPErr, *negativeIPSig;
0068     TH1 *nTracks, *allDeltaR;
0069   } plots_[N_JET_TYPES];
0070 };
0071
0072 PatBJetTrackAnalyzer::PatBJetTrackAnalyzer(const edm::ParameterSet &params)
0073     : jetsToken_(consumes<pat::JetCollection>(params.getParameter<edm::InputTag>("jets"))),
0074       tracksToken_(consumes<reco::TrackCollection>(params.getParameter<edm::InputTag>("tracks"))),
0075       beamSpotToken_(consumes<reco::BeamSpot>(params.getParameter<edm::InputTag>("beamSpot"))),
0076       primaryVerticesToken_(consumes<reco::VertexCollection>(params.getParameter<edm::InputTag>("primaryVertices"))),
0077       jetPtCut_(params.getParameter<double>("jetPtCut")),
0078       jetEtaCut_(params.getParameter<double>("jetEtaCut")),
0079       maxDeltaR_(params.getParameter<double>("maxDeltaR")),
0080       minPt_(params.getParameter<double>("minPt")),
0081       minPixelHits_(params.getParameter<unsigned int>("minPixelHits")),
0082       minTotalHits_(params.getParameter<unsigned int>("minTotalHits")),
0083       nThTrack_(params.getParameter<unsigned int>("nThTrack")) {
0084   usesResource(TFileService::kSharedResource);
0085 }
0086
0087 PatBJetTrackAnalyzer::~PatBJetTrackAnalyzer() {}
0088
0089 void PatBJetTrackAnalyzer::beginJob() {
0090   // retrieve handle to auxiliary service
0091   //  used for storing histograms into ROOT file
0092   edm::Service<TFileService> fs;
0093
0094   flavours_ = fs->make<TH1F>("flavours", "jet flavours", 5, 0, 5);
0095
0096   // book histograms for all jet flavours
0097   for (unsigned int i = 0; i < N_JET_TYPES; i++) {
0098     Plots &plots = plots_[i];
0099     const char *flavour, *name;
0100
0101     switch ((Flavour)i) {
0102       case ALL_JETS:
0103         flavour = "all jets";
0104         name = "all";
0105         break;
0106       case UDSG_JETS:
0107         flavour = "light flavour jets";
0108         name = "udsg";
0109         break;
0110       case C_JETS:
0111         flavour = "charm jets";
0112         name = "c";
0113         break;
0114       case B_JETS:
0115         flavour = "bottom jets";
0116         name = "b";
0117         break;
0118       default:
0119         flavour = "unidentified jets";
0120         name = "ni";
0121         break;
0122     }
0123
0124     plots.allIP =
0125         fs->make<TH1F>(Form("allIP_%s", name), Form("signed IP for all tracks in %s", flavour), 100, -0.1, 0.2);
0126     plots.allIPErr = fs->make<TH1F>(
0127         Form("allIPErr_%s", name), Form("error of signed IP for all tracks in %s", flavour), 100, 0, 0.05);
0128     plots.allIPSig = fs->make<TH1F>(
0129         Form("allIPSig_%s", name), Form("signed IP significance for all tracks in %s", flavour), 100, -10, 20);
0130
0131     plots.trackIP = fs->make<TH1F>(
0132         Form("trackIP_%s", name), Form("signed IP for selected positive track in %s", flavour), 100, -0.1, 0.2);
0133     plots.trackIPErr = fs->make<TH1F>(Form("trackIPErr_%s", name),
0134                                       Form("error of signed IP for selected positive track in %s", flavour),
0135                                       100,
0136                                       0,
0137                                       0.05);
0138     plots.trackIPSig = fs->make<TH1F>(Form("trackIPSig_%s", name),
0139                                       Form("signed IP significance for selected positive track in %s", flavour),
0140                                       100,
0141                                       -10,
0142                                       20);
0143
0144     plots.negativeIP = fs->make<TH1F>(
0145         Form("negativeIP_%s", name), Form("signed IP for selected negative track in %s", flavour), 100, -0.2, 0.1);
0146     plots.negativeIPErr = fs->make<TH1F>(Form("negativeIPErr_%s", name),
0147                                          Form("error of signed IP for selected negative track in %s", flavour),
0148                                          100,
0149                                          0,
0150                                          0.05);
0151     plots.negativeIPSig = fs->make<TH1F>(Form("negativeIPSig_%s", name),
0152                                          Form("signed IP significance for selected negative track in %s", flavour),
0153                                          100,
0154                                          -20,
0155                                          10);
0156
0157     plots.nTracks = fs->make<TH1F>(Form("nTracks_%s", name), Form("number of usable tracks in %s", flavour), 30, 0, 30);
0158     plots.allDeltaR =
0159         fs->make<TH1F>(Form("allDeltaR_%s", name), Form("\\DeltaR between track and %s", flavour), 100, 0, 1);
0160   }
0161 }
0162
0163 // helper function to sort the tracks by impact parameter significance
0164
0165 static bool significanceHigher(const Measurement1D &meas1, const Measurement1D &meas2) {
0166   return meas1.significance() > meas2.significance();
0167 }
0168
0169 void PatBJetTrackAnalyzer::analyze(const edm::Event &event, const edm::EventSetup &es) {
0170   // handle to the primary vertex collection
0171   edm::Handle<reco::VertexCollection> pvHandle;
0172   event.getByToken(primaryVerticesToken_, pvHandle);
0173
0174   // handle to the tracks collection
0175   edm::Handle<reco::TrackCollection> tracksHandle;
0176   event.getByToken(tracksToken_, tracksHandle);
0177
0178   // handle to the jets collection
0179   edm::Handle<pat::JetCollection> jetsHandle;
0180   event.getByToken(jetsToken_, jetsHandle);
0181
0182   // handle to the beam spot
0183   edm::Handle<reco::BeamSpot> beamSpot;
0184   event.getByToken(beamSpotToken_, beamSpot);
0185
0186   // rare case of no reconstructed primary vertex
0187   if (pvHandle->empty())
0188     return;
0189
0190   // extract the position of the (most probable) reconstructed vertex
0191   math::XYZPoint pv = (*pvHandle)[0].position();
0192
0193   // now go through all jets
0194   for (pat::JetCollection::const_iterator jet = jetsHandle->begin(); jet != jetsHandle->end(); ++jet) {
0195     // only look at jets that pass the pt and eta cut
0196     if (jet->pt() < jetPtCut_ || std::abs(jet->eta()) > jetEtaCut_)
0197       continue;
0198
0199     Flavour flavour;
0200     // find out the jet flavour (differs between quark and anti-quark)
0201     switch (std::abs(jet->partonFlavour())) {
0202       case 1:
0203       case 2:
0204       case 3:
0205       case 21:
0206         flavour = UDSG_JETS;
0207         break;
0208       case 4:
0209         flavour = C_JETS;
0210         break;
0211       case 5:
0212         flavour = B_JETS;
0213         break;
0214       default:
0215         flavour = NONID_JETS;
0216     }
0217
0218     // simply count the number of accepted jets
0219     flavours_->Fill(ALL_JETS);
0220     flavours_->Fill(flavour);
0221
0222     // this vector will contain IP value / error pairs
0223     std::vector<Measurement1D> ipValErr;
0224
0225     // Note: PAT is also able to store associated tracks
0226     //       within the jet object, so we don't have to do the
0227     //       matching ourselves
0228     // (see ->associatedTracks() method)
0229     // However, using this we can't play with the DeltaR cone
0230     // withour rerunning the PAT producer
0231
0232     // now loop through all tracks
0233     for (reco::TrackCollection::const_iterator track = tracksHandle->begin(); track != tracksHandle->end(); ++track) {
0234       // check the quality criteria
0235       if (track->pt() < minPt_ || track->hitPattern().numberOfValidHits() < (int)minTotalHits_ ||
0236           track->hitPattern().numberOfValidPixelHits() < (int)minPixelHits_)
0237         continue;
0238
0239       // check the Delta R between jet axis and track
0240       // (Delta_R^2 = Delta_Eta^2 + Delta_Phi^2)
0241       double deltaR = ROOT::Math::VectorUtil::DeltaR(jet->momentum(), track->momentum());
0242
0243       plots_[ALL_JETS].allDeltaR->Fill(deltaR);
0244       plots_[flavour].allDeltaR->Fill(deltaR);
0245
0246       // only look at tracks in jet cone
0247       if (deltaR > maxDeltaR_)
0248         continue;
0249
0250       // What follows here is an approximation!
0251       //
0252       // The dxy() method of the tracks does a linear
0253       // extrapolation from the track reference position
0254       // given as the closest point to the beam spot
0255       // with respect to the given vertex.
0256       // Since we are using primary vertices, this
0257       // approximation works well
0258       //
0259       // In order to get better results, the
0260       // "TransientTrack" and specialised methods have
0261       // to be used.
0262       // Or look at the "impactParameterTagInfos",
0263       // which contains the precomputed information
0264       // from the official b-tagging algorithms
0265       //
0266       // see ->tagInfoTrackIP() method
0267
0268       double ipError = track->dxyError();
0269       double ipValue = std::abs(track->dxy(pv));
0270
0271       // in order to compute the sign, we check if
0272       // the point of closest approach to the vertex
0273       // is in front or behind the vertex.
0274       // Again, we a linear approximation
0275       //
0276       // dot product between reference point and jet axis
0277
0278       math::XYZVector closestPoint = track->referencePoint() - beamSpot->position();
0279       // only interested in transverse component, z -> 0
0280       closestPoint.SetZ(0.);
0281       double sign = closestPoint.Dot(jet->momentum());
0282
0283       if (sign < 0)
0284         ipValue = -ipValue;
0285
0286       ipValErr.push_back(Measurement1D(ipValue, ipError));
0287     }
0288
0289     // now order all tracks by significance (highest first)
0290     std::sort(ipValErr.begin(), ipValErr.end(), significanceHigher);
0291
0292     plots_[ALL_JETS].nTracks->Fill(ipValErr.size());
0293     plots_[flavour].nTracks->Fill(ipValErr.size());
0294
0295     // plot all tracks
0296
0297     for (std::vector<Measurement1D>::const_iterator iter = ipValErr.begin(); iter != ipValErr.end(); ++iter) {
0298       plots_[ALL_JETS].allIP->Fill(iter->value());
0299       plots_[flavour].allIP->Fill(iter->value());
0300
0301       plots_[ALL_JETS].allIPErr->Fill(iter->error());
0302       plots_[flavour].allIPErr->Fill(iter->error());
0303
0304       // significance (is really just value / error)
0305       plots_[ALL_JETS].allIPSig->Fill(iter->significance());
0306       plots_[flavour].allIPSig->Fill(iter->significance());
0307     }
0308
0309     // check if we have enough tracks to fulfill the
0310     // n-th track requirement
0311     if (ipValErr.size() < nThTrack_)
0312       continue;
0313
0314     // pick the n-th highest track
0315     const Measurement1D *trk = &ipValErr[nThTrack_ - 1];
0316
0317     plots_[ALL_JETS].trackIP->Fill(trk->value());
0318     plots_[flavour].trackIP->Fill(trk->value());
0319
0320     plots_[ALL_JETS].trackIPErr->Fill(trk->error());
0321     plots_[flavour].trackIPErr->Fill(trk->error());
0322
0323     plots_[ALL_JETS].trackIPSig->Fill(trk->significance());
0324     plots_[flavour].trackIPSig->Fill(trk->significance());
0325
0326     // here we define a "negative tagger", i.e. we take
0327     // the track with the n-lowest signed IP
0328     // (i.e. preferrably select tracks that appear to become
0329     //  from "behind" the primary vertex, supposedly mismeasured
0330     //  tracks really coming from the primary vertex, and
0331     //  the contamination of displaced tracks should be small)
0332     trk = &ipValErr[ipValErr.size() - nThTrack_];
0333
0334     plots_[ALL_JETS].negativeIP->Fill(trk->value());
0335     plots_[flavour].negativeIP->Fill(trk->value());
0336
0337     plots_[ALL_JETS].negativeIPErr->Fill(trk->error());
0338     plots_[flavour].negativeIPErr->Fill(trk->error());
0339
0340     plots_[ALL_JETS].negativeIPSig->Fill(trk->significance());
0341     plots_[flavour].negativeIPSig->Fill(trk->significance());
0342   }
0343 }
0344
0345 #include "FWCore/Framework/interface/MakerMacros.h"
0346
0347 DEFINE_FWK_MODULE(PatBJetTrackAnalyzer);