Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:35:00

0001 #include <cmath>
0002 #include <map>
0003 
0004 #include <TNtuple.h>
0005 
0006 #include "FWCore/Framework/interface/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::EDAnalyzer {
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   edm::Service<TFileService> fs;
0068   ntuple = fs->make<TNtuple>("jets", "", "mc:tag:et:eta");
0069 }
0070 
0071 PileupJetAnalyzer::~PileupJetAnalyzer() {}
0072 
0073 void PileupJetAnalyzer::analyze(const edm::Event &event, const edm::EventSetup &es) {
0074   edm::Handle<JetTracksAssociationCollection> jetTracksAssoc;
0075   event.getByToken(jetTracksAssocToken, jetTracksAssoc);
0076 
0077   edm::Handle<JetFloatAssociation::Container> jetTags;
0078   event.getByToken(jetTagToken, jetTags);
0079 
0080   edm::Handle<edm::SimTrackContainer> simTracks;
0081   event.getByToken(tokenSimTrack, simTracks);
0082 
0083   edm::Handle<edm::SimVertexContainer> simVertices;
0084   event.getByToken(tokenSimVertex, simVertices);
0085 
0086   // find vertices in all jets
0087   for (JetTracksAssociationCollection::const_iterator iter = jetTracksAssoc->begin(); iter != jetTracksAssoc->end();
0088        ++iter) {
0089     const edm::RefToBase<Jet> &jet = iter->first;
0090     if (jet->energy() < jetMinE || jet->et() < jetMinEt || std::abs(jet->eta()) > jetMaxEta)
0091       continue;
0092 
0093     double tag = (*jetTags)[jet];
0094 
0095     double sigSum = 0.;
0096     double bkgSum = 0.;
0097 
0098     // note: the following code is FastSim specific right now
0099 
0100     const TrackRefVector &tracks = iter->second;
0101     for (TrackRefVector::const_iterator track = tracks.begin(); track != tracks.end(); ++track) {
0102       TrackingRecHitRef hitRef = (*track)->recHit(0);
0103       bool signal = false;
0104       unsigned int trackId;
0105 
0106       const FastTrackerRecHit *hit = dynamic_cast<const FastTrackerRecHit *>(&*hitRef);
0107 
0108       if (!hit)
0109         continue;
0110       trackId = hit->simTrackId(0);
0111 
0112       for (;;) {
0113         const SimTrack &simTrack = (*simTracks)[trackId];
0114         int genPartIndex = simTrack.genpartIndex();
0115         if (genPartIndex >= 0) {
0116           signal = true;
0117           break;
0118         }
0119 
0120         int vertIndex = simTrack.vertIndex();
0121         if (vertIndex < 0)
0122           break;
0123 
0124         const SimVertex &simVertex = (*simVertices)[vertIndex];
0125         int parentIndex = simVertex.parentIndex();
0126         if (parentIndex < 0)
0127           break;
0128 
0129         trackId = (unsigned int)parentIndex;
0130       }
0131 
0132       *(signal ? &sigSum : &bkgSum) += std::min((*track)->pt(), trackPtLimit);
0133     }
0134 
0135     if (sigSum + bkgSum < 1.0e-9)
0136       continue;
0137     double signal = sigSum / (sigSum + bkgSum);
0138 
0139     ntuple->Fill(signal, tag, jet->et(), jet->eta());
0140   }
0141 }
0142 
0143 DEFINE_FWK_MODULE(PileupJetAnalyzer);