PatExamples/plugins/PatTrackAnalyzer.cc

0001 #include <iostream>
0002 #include <cmath>
0003 #include <vector>
0004 #include <string>
0005
0006 #include <TH1.h>
0007 #include <TProfile.h>
0008
0009 #include "FWCore/Framework/interface/Frameworkfwd.h"
0010 #include "FWCore/Framework/interface/Event.h"
0011 #include "FWCore/Framework/interface/EventSetup.h"
0012 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0013 #include "FWCore/Utilities/interface/InputTag.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "FWCore/ServiceRegistry/interface/Service.h"
0016
0017 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0018
0019 #include "DataFormats/Common/interface/View.h"
0020 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0021 #include "DataFormats/TrackReco/interface/Track.h"
0022 #include "DataFormats/TrackReco/interface/HitPattern.h"
0023 #include "DataFormats/PatCandidates/interface/Muon.h"
0024
0025 class PatTrackAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0026 public:
0027   /// constructor and destructor
0028   PatTrackAnalyzer(const edm::ParameterSet &params);
0029   ~PatTrackAnalyzer() override;
0030
0031   // virtual methods called from base class EDAnalyzer
0032   void beginJob() override;
0033   void analyze(const edm::Event &event, const edm::EventSetup &es) override;
0034
0035 private:
0036   // configuration parameters
0037   edm::EDGetTokenT<edm::View<reco::Track> > srcTracksToken_;
0038   edm::EDGetTokenT<pat::MuonCollection> srcMuonsToken_;
0039   edm::EDGetTokenT<reco::BeamSpot> beamSpotToken_;
0040
0041   // the list of track quality cuts to demand from the tracking
0042   std::vector<std::string> qualities_;
0043
0044   // holder for the histograms, one set per quality flag
0045   struct Plots {
0046     TH1 *eta, *phi;
0047     TH1 *pt, *ptErr;
0048     TH1 *invPt, *invPtErr;
0049     TH1 *d0, *d0Err;
0050     TH1 *nHits;
0051
0052     TProfile *pxbHitsEta, *pxeHitsEta;
0053     TProfile *tibHitsEta, *tobHitsEta;
0054     TProfile *tidHitsEta, *tecHitsEta;
0055   };
0056
0057   std::vector<Plots> plots_;
0058 };
0059
0060 PatTrackAnalyzer::PatTrackAnalyzer(const edm::ParameterSet &params)
0061     : srcTracksToken_(mayConsume<edm::View<reco::Track> >(params.getParameter<edm::InputTag>("src"))),
0062       srcMuonsToken_(mayConsume<pat::MuonCollection>(params.getParameter<edm::InputTag>("src"))),
0063       beamSpotToken_(consumes<reco::BeamSpot>(params.getParameter<edm::InputTag>("beamSpot"))),
0064       qualities_(params.getParameter<std::vector<std::string> >("qualities")) {
0065   usesResource(TFileService::kSharedResource);
0066 }
0067
0068 PatTrackAnalyzer::~PatTrackAnalyzer() {}
0069
0070 void PatTrackAnalyzer::beginJob() {
0071   // retrieve handle to auxiliary service
0072   //  used for storing histograms into ROOT file
0073   edm::Service<TFileService> fs;
0074
0075   // now book the histograms, for each category
0076   unsigned int nQualities = qualities_.size();
0077
0078   plots_.resize(nQualities);
0079
0080   for (unsigned int i = 0; i < nQualities; ++i) {
0081     // the name of the quality flag
0082     const char *quality = qualities_[i].c_str();
0083
0084     // the set of plots
0085     Plots &plots = plots_[i];
0086
0087     plots.eta = fs->make<TH1F>(Form("eta_%s", quality), Form("track \\eta (%s)", quality), 100, -3, 3);
0088     plots.phi = fs->make<TH1F>(Form("phi_%s", quality), Form("track \\phi (%s)", quality), 100, -M_PI, +M_PI);
0089     plots.pt = fs->make<TH1F>(Form("pt_%s", quality), Form("track p_{T} (%s)", quality), 100, 0, 10);
0090     plots.ptErr = fs->make<TH1F>(Form("ptErr_%s", quality), Form("track p_{T} error (%s)", quality), 100, 0, 1);
0091     plots.invPt = fs->make<TH1F>(Form("invPt_%s", quality), Form("track 1/p_{T} (%s)", quality), 100, -5, 5);
0092     plots.invPtErr =
0093         fs->make<TH1F>(Form("invPtErr_%s", quality), Form("track 1/p_{T} error (%s)", quality), 100, 0, 0.1);
0094     plots.d0 = fs->make<TH1F>(Form("d0_%s", quality), Form("track d0 (%s)", quality), 100, 0, 0.1);
0095     plots.d0Err = fs->make<TH1F>(Form("d0Err_%s", quality), Form("track d0 error (%s)", quality), 100, 0, 0.1);
0096     plots.nHits =
0097         fs->make<TH1F>(Form("nHits_%s", quality), Form("track number of total hits (%s)", quality), 60, 0, 60);
0098
0099     plots.pxbHitsEta =
0100         fs->make<TProfile>(Form("pxbHitsEta_%s", quality), Form("#hits in Pixel Barrel (%s)", quality), 100, 0, 3);
0101     plots.pxeHitsEta =
0102         fs->make<TProfile>(Form("pxeHitsEta_%s", quality), Form("#hits in Pixel Endcap (%s)", quality), 100, 0, 3);
0103     plots.tibHitsEta = fs->make<TProfile>(
0104         Form("tibHitsEta_%s", quality), Form("#hits in Tracker Inner Barrel (%s)", quality), 100, 0, 3);
0105     plots.tobHitsEta = fs->make<TProfile>(
0106         Form("tobHitsEta_%s", quality), Form("#hits in Tracker Outer Barrel (%s)", quality), 100, 0, 3);
0107     plots.tidHitsEta = fs->make<TProfile>(
0108         Form("tidHitsEta_%s", quality), Form("#hits in Tracker Inner Disk (%s)", quality), 100, 0, 3);
0109     plots.tecHitsEta =
0110         fs->make<TProfile>(Form("tecHitsEta_%s", quality), Form("#hits in Tracker Endcap (%s)", quality), 100, 0, 3);
0111   }
0112 }
0113
0114 void PatTrackAnalyzer::analyze(const edm::Event &event, const edm::EventSetup &es) {
0115   // handles to kinds of data we might want to read
0116   edm::Handle<reco::BeamSpot> beamSpot;
0117   edm::Handle<edm::View<reco::Track> > tracksHandle;
0118   edm::Handle<pat::MuonCollection> muonsHandle;
0119
0120   // read the beam spot
0121   event.getByToken(beamSpotToken_, beamSpot);
0122
0123   // our internal copy of track points
0124   // (we need this in order to able to simultaneously access tracks
0125   //  directly or embedded in PAT objects, like muons, normally you
0126   //  would iterate over the handle directly)
0127   std::vector<const reco::Track *> tracks;
0128
0129   event.getByToken(srcTracksToken_, tracksHandle);
0130   if (tracksHandle.isValid()) {
0131     // framework was able to read the collection as a view of
0132     // tracks, no copy them to our "tracks" variable
0133     for (edm::View<reco::Track>::const_iterator iter = tracksHandle->begin(); iter != tracksHandle->end(); ++iter)
0134       tracks.push_back(&*iter);
0135   } else {
0136     // does not exist or is not a track collection
0137     // let's assume it is a collection of PAT muons
0138     event.getByToken(srcMuonsToken_, muonsHandle);
0139
0140     // and copy them over
0141     // NOTE: We are using ->globalTrack() here
0142     //       This means we are using the global fit over both
0143     //       the inner tracker and the muon stations!
0144     //       other alternatives are: innerTrack(), outerTrack()
0145     for (pat::MuonCollection::const_iterator iter = muonsHandle->begin(); iter != muonsHandle->end(); ++iter) {
0146       reco::TrackRef track = iter->globalTrack();
0147       // the muon might not be a "global" muon
0148       if (track.isNonnull())
0149         tracks.push_back(&*track);
0150     }
0151   }
0152
0153   // we are done filling the tracks into our "tracks" vector.
0154   // now analyze them, once for each track quality category
0155
0156   unsigned int nQualities = qualities_.size();
0157   for (unsigned int i = 0; i < nQualities; ++i) {
0158     // we convert the quality flag from its name as a string
0159     // to the enumeration value used by the tracking code
0160     // (which is essentially an integer number)
0161     reco::Track::TrackQuality quality = reco::Track::qualityByName(qualities_[i]);
0162
0163     // our set of plots
0164     Plots &plots = plots_[i];
0165
0166     // now loop over the tracks
0167     for (std::vector<const reco::Track *>::const_iterator iter = tracks.begin(); iter != tracks.end(); ++iter) {
0168       // this is our track
0169       const reco::Track &track = **iter;
0170
0171       // ignore tracks that fail the quality cut
0172       if (!track.quality(quality))
0173         continue;
0174
0175       // and fill all the plots
0176       plots.eta->Fill(track.eta());
0177       plots.phi->Fill(track.phi());
0178
0179       plots.pt->Fill(track.pt());
0180       plots.ptErr->Fill(track.ptError());
0181
0182       plots.invPt->Fill(track.qoverp());
0183       plots.invPtErr->Fill(track.qoverpError());
0184
0185       // the transverse IP is taken with respect to
0186       // the beam spot instead of (0, 0)
0187       // because the beam spot in CMS is not at (0, 0)
0188       plots.d0->Fill(track.dxy(beamSpot->position()));
0189       plots.d0Err->Fill(track.dxyError());
0190
0191       plots.nHits->Fill(track.numberOfValidHits());
0192
0193       // the hit pattern contains information about
0194       // which modules of the detector have been hit
0195       const reco::HitPattern &hits = track.hitPattern();
0196
0197       double absEta = std::abs(track.eta());
0198       // now fill the number of hits in a layer depending on eta
0199       plots.pxbHitsEta->Fill(absEta, hits.numberOfValidPixelBarrelHits());
0200       plots.pxeHitsEta->Fill(absEta, hits.numberOfValidPixelEndcapHits());
0201       plots.tibHitsEta->Fill(absEta, hits.numberOfValidStripTIBHits());
0202       plots.tobHitsEta->Fill(absEta, hits.numberOfValidStripTOBHits());
0203       plots.tidHitsEta->Fill(absEta, hits.numberOfValidStripTIDHits());
0204       plots.tecHitsEta->Fill(absEta, hits.numberOfValidStripTECHits());
0205     }
0206   }
0207 }
0208
0209 #include "FWCore/Framework/interface/MakerMacros.h"
0210
0211 DEFINE_FWK_MODULE(PatTrackAnalyzer);