JetAssociationProducers/test/PileupJetAnalyzer.cc

0001 #include <cmath>
0002 #include <map>
0003
0004 #include <TNtuple.h>
0005
0006 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0007 #include "FWCore/Framework/interface/Event.h"
0008 #include "FWCore/Framework/interface/EventSetup.h"
0009 #include "FWCore/Framework/interface/MakerMacros.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011 #include "FWCore/Utilities/interface/InputTag.h"
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013 #include "FWCore/Utilities/interface/Exception.h"
0014 #include "FWCore/ServiceRegistry/interface/Service.h"
0015
0016 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0017 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0018
0019 #include "DataFormats/Common/interface/Handle.h"
0020 #include "DataFormats/Common/interface/RefToBase.h"
0021 #include "DataFormats/Common/interface/RefToBaseProd.h"
0022 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0023 #include "DataFormats/JetReco/interface/Jet.h"
0024 #include "DataFormats/JetReco/interface/JetTracksAssociation.h"
0025 #include "DataFormats/JetReco/interface/JetFloatAssociation.h"
0026 #include "DataFormats/TrackReco/interface/Track.h"
0027 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0028 #include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h"
0029 #include "DataFormats/TrackerRecHit2D/interface/FastTrackerRecHit.h"
0030
0031 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0032
0033 using namespace reco;
0034
0035 class PileupJetAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0036 public:
0037   PileupJetAnalyzer(const edm::ParameterSet &params);
0038   ~PileupJetAnalyzer();
0039
0040   virtual void analyze(const edm::Event &event, const edm::EventSetup &es);
0041
0042 private:
0043   TNtuple *ntuple;
0044
0045   edm::EDGetTokenT<JetTracksAssociationCollection> jetTracksAssocToken;
0046   edm::EDGetTokenT<JetFloatAssociation::Container> jetTagToken;
0047   edm::EDGetTokenT<edm::SimTrackContainer> tokenSimTrack;
0048   edm::EDGetTokenT<edm::SimVertexContainer> tokenSimVertex;
0049   const double signalFraction;
0050   const double jetMinE;
0051   const double jetMinEt;
0052   const double jetMaxEta;
0053   const double trackPtLimit;
0054 };
0055
0056 PileupJetAnalyzer::PileupJetAnalyzer(const edm::ParameterSet &params)
0057     : signalFraction(params.getParameter<double>("signalFraction")),
0058       jetMinE(params.getParameter<double>("jetMinE")),
0059       jetMinEt(params.getParameter<double>("jetMinEt")),
0060       jetMaxEta(params.getParameter<double>("jetMaxEta")),
0061       trackPtLimit(params.getParameter<double>("trackPtLimit")) {
0062   jetTracksAssocToken = consumes<JetTracksAssociationCollection>(params.getParameter<edm::InputTag>("jetTracksAssoc"));
0063   jetTagToken = consumes<JetFloatAssociation::Container>(params.getParameter<edm::InputTag>("jetTagLabel"));
0064   tokenSimTrack = consumes<edm::SimTrackContainer>(edm::InputTag("fastSimProducer"));
0065   tokenSimVertex = consumes<edm::SimVertexContainer>(edm::InputTag("fastSimProducer"));
0066
0067   usesResource(TFileService::kSharedResource);
0068   edm::Service<TFileService> fs;
0069   ntuple = fs->make<TNtuple>("jets", "", "mc:tag:et:eta");
0070 }
0071
0072 PileupJetAnalyzer::~PileupJetAnalyzer() {}
0073
0074 void PileupJetAnalyzer::analyze(const edm::Event &event, const edm::EventSetup &es) {
0075   edm::Handle<JetTracksAssociationCollection> jetTracksAssoc;
0076   event.getByToken(jetTracksAssocToken, jetTracksAssoc);
0077
0078   edm::Handle<JetFloatAssociation::Container> jetTags;
0079   event.getByToken(jetTagToken, jetTags);
0080
0081   edm::Handle<edm::SimTrackContainer> simTracks;
0082   event.getByToken(tokenSimTrack, simTracks);
0083
0084   edm::Handle<edm::SimVertexContainer> simVertices;
0085   event.getByToken(tokenSimVertex, simVertices);
0086
0087   // find vertices in all jets
0088   for (JetTracksAssociationCollection::const_iterator iter = jetTracksAssoc->begin(); iter != jetTracksAssoc->end();
0089        ++iter) {
0090     const edm::RefToBase<Jet> &jet = iter->first;
0091     if (jet->energy() < jetMinE || jet->et() < jetMinEt || std::abs(jet->eta()) > jetMaxEta)
0092       continue;
0093
0094     double tag = (*jetTags)[jet];
0095
0096     double sigSum = 0.;
0097     double bkgSum = 0.;
0098
0099     // note: the following code is FastSim specific right now
0100
0101     const TrackRefVector &tracks = iter->second;
0102     for (TrackRefVector::const_iterator track = tracks.begin(); track != tracks.end(); ++track) {
0103       TrackingRecHitRef hitRef = (*track)->recHit(0);
0104       bool signal = false;
0105       unsigned int trackId;
0106
0107       const FastTrackerRecHit *hit = dynamic_cast<const FastTrackerRecHit *>(&*hitRef);
0108
0109       if (!hit)
0110         continue;
0111       trackId = hit->simTrackId(0);
0112
0113       for (;;) {
0114         const SimTrack &simTrack = (*simTracks)[trackId];
0115         int genPartIndex = simTrack.genpartIndex();
0116         if (genPartIndex >= 0) {
0117           signal = true;
0118           break;
0119         }
0120
0121         int vertIndex = simTrack.vertIndex();
0122         if (vertIndex < 0)
0123           break;
0124
0125         const SimVertex &simVertex = (*simVertices)[vertIndex];
0126         int parentIndex = simVertex.parentIndex();
0127         if (parentIndex < 0)
0128           break;
0129
0130         trackId = (unsigned int)parentIndex;
0131       }
0132
0133       *(signal ? &sigSum : &bkgSum) += std::min((*track)->pt(), trackPtLimit);
0134     }
0135
0136     if (sigSum + bkgSum < 1.0e-9)
0137       continue;
0138     double signal = sigSum / (sigSum + bkgSum);
0139
0140     ntuple->Fill(signal, tag, jet->et(), jet->eta());
0141   }
0142 }
0143
0144 DEFINE_FWK_MODULE(PileupJetAnalyzer);