/data/refman/pasoursint/CMSSW_4_4_5_patch3/src/PhysicsTools/PatExamples/plugins/PatTrackAnalyzer.cc

Go to the documentation of this file.

#include <iostream>
#include <cmath>
#include <vector>
#include <string>

#include <TH1.h>
#include <TProfile.h>

#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/EDAnalyzer.h"
#include "FWCore/Utilities/interface/InputTag.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ServiceRegistry/interface/Service.h"

#include "CommonTools/UtilAlgos/interface/TFileService.h"

#include "DataFormats/Common/interface/View.h"
#include "DataFormats/BeamSpot/interface/BeamSpot.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/HitPattern.h"
#include "DataFormats/PatCandidates/interface/Muon.h"

class PatTrackAnalyzer : public edm::EDAnalyzer  {
    public: 
        PatTrackAnalyzer(const edm::ParameterSet &params);
        ~PatTrackAnalyzer();

        // virtual methods called from base class EDAnalyzer
        virtual void beginJob();
        virtual void analyze(const edm::Event &event, const edm::EventSetup &es);

    private:
        // configuration parameters
        edm::InputTag src_;
        edm::InputTag beamSpot_;

        // the list of track quality cuts to demand from the tracking
        std::vector<std::string> qualities_;

        // holder for the histograms, one set per quality flag
        struct Plots {
                TH1 *eta, *phi;
                TH1 *pt, *ptErr;
                TH1 *invPt, *invPtErr;
                TH1 *d0, *d0Err;
                TH1 *nHits;

                TProfile *pxbHitsEta, *pxeHitsEta;
                TProfile *tibHitsEta, *tobHitsEta;
                TProfile *tidHitsEta, *tecHitsEta;
        };

        std::vector<Plots> plots_;
};


PatTrackAnalyzer::PatTrackAnalyzer(const edm::ParameterSet &params) :
        src_(params.getParameter<edm::InputTag>("src")),
        beamSpot_(params.getParameter<edm::InputTag>("beamSpot")),
        qualities_(params.getParameter< std::vector<std::string> >("qualities"))
{
}

PatTrackAnalyzer::~PatTrackAnalyzer()
{
}

void PatTrackAnalyzer::beginJob()
{
        // retrieve handle to auxiliary service
        //  used for storing histograms into ROOT file
        edm::Service<TFileService> fs;

        // now book the histograms, for each category
        unsigned int nQualities = qualities_.size();

        plots_.resize(nQualities);

        for(unsigned int i = 0; i < nQualities; ++i) {
                // the name of the quality flag
                const char *quality = qualities_[i].c_str();

                // the set of plots
                Plots &plots = plots_[i];

                plots.eta = fs->make<TH1F>(Form("eta_%s", quality),
                                           Form("track \\eta (%s)", quality),
                                           100, -3, 3);
                plots.phi = fs->make<TH1F>(Form("phi_%s", quality),
                                           Form("track \\phi (%s)", quality),
                                           100, -M_PI, +M_PI);
                plots.pt = fs->make<TH1F>(Form("pt_%s", quality),
                                          Form("track p_{T} (%s)", quality),
                                          100, 0, 10);
                plots.ptErr = fs->make<TH1F>(Form("ptErr_%s", quality),
                                             Form("track p_{T} error (%s)", quality),
                                             100, 0, 1);
                plots.invPt = fs->make<TH1F>(Form("invPt_%s", quality),
                                             Form("track 1/p_{T} (%s)", quality),
                                             100, -5, 5);
                plots.invPtErr = fs->make<TH1F>(Form("invPtErr_%s", quality),
                                                Form("track 1/p_{T} error (%s)", quality),
                                                100, 0, 0.1);
                plots.d0 = fs->make<TH1F>(Form("d0_%s", quality),
                                          Form("track d0 (%s)", quality),
                                                100, 0, 0.1);
                plots.d0Err = fs->make<TH1F>(Form("d0Err_%s", quality),
                                             Form("track d0 error (%s)", quality),
                                                   100, 0, 0.1);
                plots.nHits = fs->make<TH1F>(Form("nHits_%s", quality),
                                             Form("track number of total hits (%s)", quality),
                                                   60, 0, 60);

                plots.pxbHitsEta = fs->make<TProfile>(Form("pxbHitsEta_%s", quality),
                                                      Form("#hits in Pixel Barrel (%s)", quality),
                                                           100, 0, 3);
                plots.pxeHitsEta = fs->make<TProfile>(Form("pxeHitsEta_%s", quality),
                                                      Form("#hits in Pixel Endcap (%s)", quality),
                                                           100, 0, 3);
                plots.tibHitsEta = fs->make<TProfile>(Form("tibHitsEta_%s", quality),
                                                      Form("#hits in Tracker Inner Barrel (%s)", quality),
                                                           100, 0, 3);
                plots.tobHitsEta = fs->make<TProfile>(Form("tobHitsEta_%s", quality),
                                                      Form("#hits in Tracker Outer Barrel (%s)", quality),
                                                           100, 0, 3);
                plots.tidHitsEta = fs->make<TProfile>(Form("tidHitsEta_%s", quality),
                                                      Form("#hits in Tracker Inner Disk (%s)", quality),
                                                           100, 0, 3);
                plots.tecHitsEta = fs->make<TProfile>(Form("tecHitsEta_%s", quality),
                                                      Form("#hits in Tracker Endcap (%s)", quality),
                                                           100, 0, 3);
        }
}

void PatTrackAnalyzer::analyze(const edm::Event &event, const edm::EventSetup &es)
{  
        // handles to kinds of data we might want to read
        edm::Handle<reco::BeamSpot> beamSpot;
        edm::Handle< edm::View<reco::Track> > tracksHandle;
        edm::Handle< pat::MuonCollection > muonsHandle;

        // read the beam spot
        event.getByLabel(beamSpot_, beamSpot);

        // our internal copy of track points
        // (we need this in order to able to simultaneously access tracks
        //  directly or embedded in PAT objects, like muons, normally you
        //  would iterate over the handle directly)
        std::vector<const reco::Track*> tracks;

        event.getByLabel(src_, tracksHandle);
        if (tracksHandle.isValid()) {
                // framework was able to read the collection as a view of
                // tracks, no copy them to our "tracks" variable
                for(edm::View<reco::Track>::const_iterator iter = tracksHandle->begin();
                    iter != tracksHandle->end(); ++iter)
                        tracks.push_back(&*iter);
        } else {
                // does not exist or is not a track collection
                // let's assume it is a collection of PAT muons
                event.getByLabel(src_, muonsHandle);

                // and copy them over
                // NOTE: We are using ->globalTrack() here
                //       This means we are using the global fit over both
                //       the inner tracker and the muon stations!
                //       other alternatives are: innerTrack(), outerTrack()
                for(pat::MuonCollection::const_iterator iter = muonsHandle->begin();
                    iter != muonsHandle->end(); ++iter) {
                        reco::TrackRef track = iter->globalTrack();
                        // the muon might not be a "global" muon
                        if (track.isNonnull())
                                tracks.push_back(&*track);
                }
        }

        // we are done filling the tracks into our "tracks" vector.
        // now analyze them, once for each track quality category

        unsigned int nQualities = qualities_.size();
        for(unsigned int i = 0; i < nQualities; ++i) {
                // we convert the quality flag from its name as a string
                // to the enumeration value used by the tracking code
                // (which is essentially an integer number)
                reco::Track::TrackQuality quality = reco::Track::qualityByName(qualities_[i]);

                // our set of plots
                Plots &plots = plots_[i];

                // now loop over the tracks
                for(std::vector<const reco::Track*>::const_iterator iter = tracks.begin();
                    iter != tracks.end(); ++iter) {
                        // this is our track
                        const reco::Track &track = **iter;

                        // ignore tracks that fail the quality cut
                        if (!track.quality(quality))
                                continue;

                        // and fill all the plots
                        plots.eta->Fill(track.eta());
                        plots.phi->Fill(track.phi());

                        plots.pt->Fill(track.pt());
                        plots.ptErr->Fill(track.ptError());

                        plots.invPt->Fill(track.qoverp());
                        plots.invPtErr->Fill(track.qoverpError());

                        // the transverse IP is taken with respect to
                        // the beam spot instead of (0, 0)
                        // because the beam spot in CMS is not at (0, 0)
                        plots.d0->Fill(track.dxy(beamSpot->position()));
                        plots.d0Err->Fill(track.dxyError());

                        plots.nHits->Fill(track.numberOfValidHits());

                        // the hit pattern contains information about
                        // which modules of the detector have been hit
                        const reco::HitPattern &hits = track.hitPattern();

                        double absEta = std::abs(track.eta());
                        // now fill the number of hits in a layer depending on eta
                        plots.pxbHitsEta->Fill(absEta, hits.numberOfValidPixelBarrelHits());
                        plots.pxeHitsEta->Fill(absEta, hits.numberOfValidPixelEndcapHits());
                        plots.tibHitsEta->Fill(absEta, hits.numberOfValidStripTIBHits());
                        plots.tobHitsEta->Fill(absEta, hits.numberOfValidStripTOBHits());
                        plots.tidHitsEta->Fill(absEta, hits.numberOfValidStripTIDHits());
                        plots.tecHitsEta->Fill(absEta, hits.numberOfValidStripTECHits());
                }
        }
}
        
#include "FWCore/Framework/interface/MakerMacros.h"

DEFINE_FWK_MODULE(PatTrackAnalyzer);