/data/refman/pasoursint/CMSSW_5_3_3/src/PhysicsTools/PatExamples/plugins/PatBJetVertexAnalyzer.cc

Go to the documentation of this file.

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

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

#include <Math/VectorUtil.h>
#include <Math/GenVector/PxPyPzE4D.h>
#include <Math/GenVector/PxPyPzM4D.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/BTauReco/interface/SecondaryVertexTagInfo.h"
#include "DataFormats/Math/interface/LorentzVector.h"
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "DataFormats/VertexReco/interface/VertexFwd.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "DataFormats/PatCandidates/interface/Jet.h"
#include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"

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

        // 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 jets_;

        double jetPtCut_;               // minimum (uncorrected) jet energy
        double jetEtaCut_;              // maximum |eta| for jet
        double maxDeltaR_;              // angle between jet and tracks

        // jet flavour constants

        enum Flavour {
                ALL_JETS = 0,
                UDSG_JETS,
                C_JETS,
                B_JETS,
                NONID_JETS,
                N_JET_TYPES
        };

        TH1 * flavours_;

        // one group of plots per jet flavour;
        struct Plots {
                TH1 *nVertices;
                TH1 *deltaR, *mass, *dist, *distErr, *distSig;
                TH1 *nTracks, *chi2;
        } plots_[N_JET_TYPES];
};

PatBJetVertexAnalyzer::PatBJetVertexAnalyzer(const edm::ParameterSet &params) :
        jets_(params.getParameter<edm::InputTag>("jets")),
        jetPtCut_(params.getParameter<double>("jetPtCut")),
        jetEtaCut_(params.getParameter<double>("jetEtaCut"))
{
}

PatBJetVertexAnalyzer::~PatBJetVertexAnalyzer()
{
}

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

        flavours_ = fs->make<TH1F>("flavours", "jet flavours", 5, 0, 5);

        // book histograms for all jet flavours
        for(unsigned int i = 0; i < N_JET_TYPES; i++) {
                Plots &plots = plots_[i];
                const char *flavour, *name;

                switch((Flavour)i) {
                    case ALL_JETS:
                        flavour = "all jets";
                        name = "all";
                        break;
                    case UDSG_JETS:
                        flavour = "light flavour jets";
                        name = "udsg";
                        break;
                    case C_JETS:
                        flavour = "charm jets";
                        name = "c";
                        break;
                    case B_JETS:
                        flavour = "bottom jets";
                        name = "b";
                        break;
                    default:
                        flavour = "unidentified jets";
                        name = "ni";
                        break;
                }

                plots.nVertices = fs->make<TH1F>(Form("nVertices_%s", name),
                                                 Form("number of secondary vertices in %s", flavour),
                                                 5, 0, 5);
                plots.deltaR = fs->make<TH1F>(Form("deltaR_%s", name),
                                              Form("\\DeltaR between vertex direction and jet direction in %s", flavour),
                                              100, 0, 0.5);
                plots.mass = fs->make<TH1F>(Form("mass_%s", name),
                                            Form("vertex mass in %s", flavour),
                                            100, 0, 10);
                plots.dist = fs->make<TH1F>(Form("dist_%s", name),
                                            Form("Transverse distance between PV and SV in %s", flavour),
                                            100, 0, 2);
                plots.distErr = fs->make<TH1F>(Form("distErr_%s", name),
                                               Form("Transverse distance error between PV and SV in %s", flavour),
                                               100, 0, 0.5);
                plots.distSig = fs->make<TH1F>(Form("distSig_%s", name),
                                               Form("Transverse distance significance between PV and SV in %s", flavour),
                                               100, 0, 50);
                plots.nTracks = fs->make<TH1F>(Form("nTracks_%s", name),
                                               Form("number of tracks at secondary vertex in %s", flavour),
                                               20, 0, 20);
                plots.chi2 = fs->make<TH1F>(Form("chi2_%s", name),
                                            Form("secondary vertex fit \\chi^{2} in %s", flavour),
                                            100, 0, 50);
        }
}

// helper function to sort the tracks by impact parameter significance

void PatBJetVertexAnalyzer::analyze(const edm::Event &event, const edm::EventSetup &es)
{  
        // handle to the jets collection
        edm::Handle<pat::JetCollection> jetsHandle;
        event.getByLabel(jets_, jetsHandle);

        // now go through all jets
        for(pat::JetCollection::const_iterator jet = jetsHandle->begin();
            jet != jetsHandle->end(); ++jet) {

                // only look at jets that pass the pt and eta cut
                if (jet->pt() < jetPtCut_ ||
                    std::abs(jet->eta()) > jetEtaCut_)
                        continue;

                Flavour flavour;
                // find out the jet flavour (differs between quark and anti-quark)
                switch(std::abs(jet->partonFlavour())) {
                    case 1:
                    case 2:
                    case 3:
                    case 21:
                        flavour = UDSG_JETS;
                        break;
                    case 4:
                        flavour = C_JETS;
                        break;
                    case 5:
                        flavour = B_JETS;
                        break;
                    default:
                        flavour = NONID_JETS;
                }

                // simply count the number of accepted jets
                flavours_->Fill(ALL_JETS);
                flavours_->Fill(flavour);

                // retrieve the "secondary vertex tag infos"
                // this is the output of the b-tagging reconstruction code
                // and contains secondary vertices in the jets
                const reco::SecondaryVertexTagInfo &svTagInfo =
                                        *jet->tagInfoSecondaryVertex();

                // count the number of secondary vertices
                plots_[ALL_JETS].nVertices->Fill(svTagInfo.nVertices());
                plots_[flavour].nVertices->Fill(svTagInfo.nVertices());

                // ignore jets without SV from now on
                if (svTagInfo.nVertices() < 1)
                        continue;

                // pick the first secondary vertex (the "best" one)
                const reco::Vertex &sv = svTagInfo.secondaryVertex(0);

                // and plot number of tracks and chi^2
                plots_[ALL_JETS].nTracks->Fill(sv.tracksSize());
                plots_[flavour].nTracks->Fill(sv.tracksSize());

                plots_[ALL_JETS].chi2->Fill(sv.chi2());
                plots_[flavour].chi2->Fill(sv.chi2());

                // the precomputed transverse distance to the primary vertex
                Measurement1D distance = svTagInfo.flightDistance(0, true);

                plots_[ALL_JETS].dist->Fill(distance.value());
                plots_[flavour].dist->Fill(distance.value());

                plots_[ALL_JETS].distErr->Fill(distance.error());
                plots_[flavour].distErr->Fill(distance.error());

                plots_[ALL_JETS].distSig->Fill(distance.significance());
                plots_[flavour].distSig->Fill(distance.significance());


                // the precomputed direction with respect to the primary vertex
                GlobalVector dir = svTagInfo.flightDirection(0);

                // unfortunately CMSSW hsa all kinds of vectors,
                // and sometimes we need to convert them *sigh*
                math::XYZVector dir2(dir.x(), dir.y(), dir.z());

                // compute a few variables that we are plotting
                double deltaR = ROOT::Math::VectorUtil::DeltaR(
                                                jet->momentum(), dir2);

                plots_[ALL_JETS].deltaR->Fill(deltaR);
                plots_[flavour].deltaR->Fill(deltaR);

                // compute the invariant mass from a four-vector sum
                math::XYZTLorentzVector trackFourVectorSum;

                // loop over all tracks in the vertex
                for(reco::Vertex::trackRef_iterator track = sv.tracks_begin();
                    track != sv.tracks_end(); ++track) {
                        ROOT::Math::LorentzVector<ROOT::Math::PxPyPzM4D<double> > vec;
                        vec.SetPx((*track)->px());
                        vec.SetPy((*track)->py());
                        vec.SetPz((*track)->pz());
                        vec.SetM(0.13957);      // pion mass
                        trackFourVectorSum += vec;
                }

                // get the invariant mass: sqrt(E² - px² - py² - pz²)
                double vertexMass = trackFourVectorSum.M();

                plots_[ALL_JETS].mass->Fill(vertexMass);
                plots_[flavour].mass->Fill(vertexMass);
        }
}
        
#include "FWCore/Framework/interface/MakerMacros.h"

DEFINE_FWK_MODULE(PatBJetVertexAnalyzer);