File indexing completed on 2023-03-17 11:18:27
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 ¶ms);
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 ¶ms)
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
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
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);