CMS 3D CMS Logo

PatBasicFWLiteJetAnalyzer_Selector.cc
Go to the documentation of this file.
1 /* A macro for making a histogram of Jet Pt with cuts
2 This is a basic way to cut out jets of a certain Pt and Eta using an if statement
3 This example creates a histogram of Jet Pt, using Jets with Pt above 30 and ETA above -2.1 and below 2.1
4 */
5 
6 #include "TFile.h"
7 #include "TH1.h"
8 #include "TH2.h"
9 #include "TCanvas.h"
10 #include "TLegend.h"
11 #include "TSystem.h"
12 #include "TLorentzVector.h"
13 
30 
31 
32 #include <iostream>
33 #include <memory>
34 #include <cmath> //necessary for absolute function fabs()
35 
36 using namespace std;
37 
41 
42 // -*- C++ -*-
43 
44 // CMS includes
50 
51 
52 // Root includes
53 #include "TROOT.h"
54 
55 using namespace std;
56 
57 
59 public:
60  JetIDStudiesSelector( edm::ParameterSet const & caloJetIdParams,
61  edm::ParameterSet const & pfJetIdParams,
62  edm::ParameterSet const & params ) :
63  jetSel_ (new JetIDSelectionFunctor(caloJetIdParams)),
64  pfJetSel_ (new PFJetIDSelectionFunctor(pfJetIdParams)),
65  jetSrc_ (params.getParameter<edm::InputTag>("jetSrc")),
66  pfJetSrc_ (params.getParameter<edm::InputTag>("pfJetSrc"))
67  {
68  bool useCalo = params.getParameter<bool>("useCalo");
69 
70  push_back("Calo Cuts");
71  push_back("Calo Kin Cuts");
72  push_back("Calo Delta Phi");
73  push_back("Calo Jet ID");
74  push_back("PF Cuts");
75  push_back("PF Kin Cuts");
76  push_back("PF Delta Phi");
77  push_back("PF Jet ID");
78 
79  set("Calo Cuts", useCalo);
80  set("Calo Kin Cuts", useCalo);
81  set("Calo Delta Phi", useCalo);
82  set("Calo Jet ID", useCalo);
83  set("PF Cuts", !useCalo);
84  set("PF Kin Cuts", !useCalo);
85  set("PF Delta Phi", !useCalo);
86  set("PF Jet ID", !useCalo);
87 
88  // Indices for fast caching
89  caloCuts_ = index_type(&bits_, std::string("Calo Cuts") );
90  caloKin_ = index_type(&bits_, std::string("Calo Kin Cuts"));
91  caloDeltaPhi_ = index_type(&bits_, std::string("Calo Delta Phi"));
92  caloJetID_ = index_type(&bits_, std::string("Calo Jet ID"));
93  pfCuts_ = index_type(&bits_, std::string("PF Cuts"));
94  pfKin_ = index_type(&bits_, std::string("PF Kin Cuts"));
95  pfDeltaPhi_ = index_type(&bits_, std::string("PF Delta Phi"));
96  pfJetID_ = index_type(&bits_, std::string("PF Jet ID"));
97 
98  }
99 
100  ~JetIDStudiesSelector() override {}
101 
102  bool operator()( edm::EventBase const & event, pat::strbitset & ret) override{
103 
104  pat::strbitset retCaloJet = jetSel_->getBitTemplate();
105  pat::strbitset retPFJet = pfJetSel_->getBitTemplate();
106 
107 
108  if ( considerCut(caloCuts_) ) {
109  passCut(ret, caloCuts_);
110  event.getByLabel( jetSrc_, h_jets_ );
111  // Calo Cuts
112  if ( h_jets_->size() >= 2 || ignoreCut(caloKin_) ) {
113  passCut(ret, caloKin_);
114  pat::Jet const & jet0 = h_jets_->at(0);
115  pat::Jet const & jet1 = h_jets_->at(1);
116  double dphi = fabs(deltaPhi<double>( jet0.phi(),
117  jet1.phi() ) );
118 
119  if ( fabs(dphi - TMath::Pi()) < 1.0 || ignoreCut(caloDeltaPhi_) ) {
120  passCut(ret, caloDeltaPhi_);
121 
122 
123  retCaloJet.set(false);
124  bool pass0 = (*jetSel_)( jet0, retCaloJet );
125  retCaloJet.set(false);
126  bool pass1 = (*jetSel_)( jet1, retCaloJet );
127  if ( (pass0 && pass1) || ignoreCut(caloJetID_) ) {
128  passCut(ret, caloJetID_);
129  caloJet0_ = edm::Ptr<pat::Jet>( h_jets_, 0);
130  caloJet1_ = edm::Ptr<pat::Jet>( h_jets_, 1);
131 
132  return (bool)ret;
133  }// end if found 2 "loose" jet ID jets
134  }// end if delta phi
135  }// end calo kin cuts
136  }// end if calo cuts
137 
138 
139  if ( considerCut(pfCuts_) ) {
140 
141  passCut(ret, pfCuts_);
142  event.getByLabel( pfJetSrc_, h_pfjets_ );
143  // PF Cuts
144  if ( h_pfjets_->size() >= 2 || ignoreCut(pfKin_) ) {
145  passCut( ret, pfKin_);
146  pat::Jet const & jet0 = h_pfjets_->at(0);
147  pat::Jet const & jet1 = h_pfjets_->at(1);
148  double dphi = fabs(deltaPhi<double>( jet0.phi(),
149  jet1.phi() ) );
150 
151  if ( fabs(dphi - TMath::Pi()) < 1.0 || ignoreCut(pfDeltaPhi_) ) {
152  passCut(ret, pfDeltaPhi_);
153 
154 
155  retPFJet.set(false);
156  bool pass0 = (*pfJetSel_)( jet0, retPFJet );
157  retPFJet.set(false);
158  bool pass1 = (*pfJetSel_)( jet1, retPFJet );
159  if ( (pass0 && pass1) || ignoreCut(pfJetID_) ) {
160  passCut(ret, pfJetID_);
161  pfJet0_ = edm::Ptr<pat::Jet>( h_pfjets_, 0);
162  pfJet1_ = edm::Ptr<pat::Jet>( h_pfjets_, 1);
163 
164  return (bool)ret;
165  }// end if found 2 "loose" jet ID jets
166  }// end if delta phi
167  }// end pf kin cuts
168  }// end if pf cuts
169 
170 
171  setIgnored(ret);
172 
173  return false;
174  }// end of method
175 
176  std::shared_ptr<JetIDSelectionFunctor> const & jetSel() const { return jetSel_;}
177  std::shared_ptr<PFJetIDSelectionFunctor> const & pfJetSel() const { return pfJetSel_;}
178 
179  vector<pat::Jet> const & allCaloJets () const { return *h_jets_; }
180  vector<pat::Jet> const & allPFJets () const { return *h_pfjets_; }
181 
182  pat::Jet const & caloJet0() const { return *caloJet0_; }
183  pat::Jet const & caloJet1() const { return *caloJet1_; }
184 
185  pat::Jet const & pfJet0() const { return *pfJet0_; }
186  pat::Jet const & pfJet1() const { return *pfJet1_; }
187 
188 
189  // Fast caching indices
190  index_type const & caloCuts() const {return caloCuts_;}
191  index_type const & caloKin() const {return caloKin_;}
192  index_type const & caloDeltaPhi() const {return caloDeltaPhi_;}
193  index_type const & caloJetID() const {return caloJetID_;}
194  index_type const & pfCuts() const {return pfCuts_;}
195  index_type const & pfKin() const {return pfKin_;}
196  index_type const & pfDeltaPhi() const {return pfDeltaPhi_;}
197  index_type const & pfJetID() const {return pfJetID_;}
198 
199 
200 protected:
201  std::shared_ptr<JetIDSelectionFunctor> jetSel_;
202  std::shared_ptr<PFJetIDSelectionFunctor> pfJetSel_;
205 
208 
211 
214 
215  // Fast caching indices
216  index_type caloCuts_;
217  index_type caloKin_;
218  index_type caloDeltaPhi_;
219  index_type caloJetID_;
220  index_type pfCuts_;
221  index_type pfKin_;
222  index_type pfDeltaPhi_;
223  index_type pfJetID_;
224 
225 };
226 
228 // ///////////////////// //
229 // // Main Subroutine // //
230 // ///////////////////// //
232 
233 int main (int argc, char* argv[])
234 {
235 
236 
237  if ( argc < 2 ) {
238  std::cout << "Usage : " << argv[0] << " [parameters.py]" << std::endl;
239  return 0;
240  }
241 
242  cout << "Hello from " << argv[0] << "!" << endl;
243 
244 
245  // load framework libraries
246  gSystem->Load( "libFWCoreFWLite" );
248 
249 
250  cout << "Getting parameters" << endl;
251  // Get the python configuration
252  PythonProcessDesc builder(argv[1]);
253  std::shared_ptr<edm::ProcessDesc> b = builder.processDesc();
254  std::shared_ptr<edm::ParameterSet> parameters = b->getProcessPSet();
255  parameters->registerIt();
256 
257  edm::ParameterSet const& jetStudiesParams = parameters->getParameter<edm::ParameterSet>("jetStudies");
258  edm::ParameterSet const& pfJetStudiesParams = parameters->getParameter<edm::ParameterSet>("pfJetStudies");
259  edm::ParameterSet const& caloJetIDParameters = parameters->getParameter<edm::ParameterSet>("jetIDSelector");
260  edm::ParameterSet const& pfJetIDParameters = parameters->getParameter<edm::ParameterSet>("pfJetIDSelector");
261  edm::ParameterSet const& plotParameters = parameters->getParameter<edm::ParameterSet>("plotParameters");
262  edm::ParameterSet const& inputs = parameters->getParameter<edm::ParameterSet>("inputs");
263  edm::ParameterSet const& outputs = parameters->getParameter<edm::ParameterSet>("outputs");
264 
265  cout << "Making RunLumiSelector" << endl;
266  RunLumiSelector runLumiSel( inputs );
267 
268  cout << "setting up TFileService" << endl;
269  // book a set of histograms
270  fwlite::TFileService fs = fwlite::TFileService( outputs.getParameter<std::string>("outputName") );
271  TFileDirectory theDir = fs.mkdir( "histos" );
272 
273  cout << "Setting up chain event" << endl;
274  // This object 'event' is used both to get all information from the
275  // event as well as to store histograms, etc.
276  fwlite::ChainEvent ev ( inputs.getParameter<std::vector<std::string> > ("fileNames") );
277 
278 
279  cout << "Booking histograms" << endl;
280  // Book histograms
281 
282  std::map<std::string, TH1*> hists;
283 
284  hists["hist_nJet" ] = theDir.make<TH1D>( "hist_nJet" ,"Number of Calo Jets", 10, 0, 10 ) ;
285  hists["hist_nPFJet" ] = theDir.make<TH1D>( "hist_nPFJet" ,"Number of PF Jets", 10, 0, 10 ) ;
286 
287  hists["hist_jetPt" ] = theDir.make<TH1D>( "hist_jetPt" , "Jet p_{T}", 400, 0, 400 ) ;
288  hists["hist_jetEtaVsPhi" ] = theDir.make<TH2D>( "hist_jetEtaVsPhi" , "Jet #phi versus #eta;#eta;#phi", 50, -5.0, 5.0, 50, -TMath::Pi(), TMath::Pi() ) ;
289  hists["hist_jetNTracks" ] = theDir.make<TH1D>( "hist_jetNTracks" , "Jet N_{TRACKS}", 20, 0, 20 ) ;
290  hists["hist_jetNTracksVsPt" ] = theDir.make<TH2D>( "hist_jetNTracksVsPt" , "Number of Tracks versus Jet p_{T};Jet p_{T}(GeV/c) ;N_{Tracks}",20, 0, 200, 20, 0, 20 ) ;
291  hists["hist_jetEMF" ] = theDir.make<TH1D>( "hist_jetEMF" , "Jet EMF", 200, 0, 1) ;
292  hists["hist_jetPdgID" ] = theDir.make<TH1D>( "hist_jetPdgID" , "PDG Id of Jet Constituents", 10000, 0, 10000 ) ;
293  hists["hist_jetGenEmE" ] = theDir.make<TH1D>( "hist_jetGenEmE" , "Gen Jet EM Energy", 200, 0, 200 ) ;
294  hists["hist_jetGenHadE" ] = theDir.make<TH1D>( "hist_jetGenHadE" , "Gen Jet HAD Energy", 200, 0, 200 ) ;
295  hists["hist_jetGenEMF" ] = theDir.make<TH1D>( "hist_jetGenEMF" , "Gen Jet EMF", 200, 0, 1) ;
296  hists["hist_jetEoverGenE" ] = theDir.make<TH1D>( "hist_jetEoverGenE" , "Energy of reco Jet / Energy of gen Jet", 200, 0, 2.0) ;
297  hists["hist_jetCorr" ] = theDir.make<TH1D>( "hist_jetCorr" , "Jet Correction Factor", 200, 0, 1.0 ) ;
298  hists["hist_n90Hits" ] = theDir.make<TH1D>( "hist_n90Hits" , "Jet n90Hits", 20, 0, 20) ;
299  hists["hist_fHPD" ] = theDir.make<TH1D>( "hist_fHPD" , "Jet fHPD", 200, 0, 1) ;
300  hists["hist_nConstituents" ] = theDir.make<TH1D>( "hist_nConstituents" , "Jet nConstituents", 20, 0, 20 ) ;
301  hists["hist_jetCHF" ] = theDir.make<TH1D>( "hist_jetCHF" , "Jet Charged Hadron Fraction", 200, 0, 1.0) ;
302 
303  hists["hist_good_jetPt" ] = theDir.make<TH1D>( "hist_good_jetPt" , "Jet p_{T}", 400, 0, 400 ) ;
304  hists["hist_good_jetEtaVsPhi" ] = theDir.make<TH2D>( "hist_good_jetEtaVsPhi" , "Jet #phi versus #eta;#eta;#phi", 50, -5.0, 5.0, 50, -TMath::Pi(), TMath::Pi() ) ;
305  hists["hist_good_jetNTracks" ] = theDir.make<TH1D>( "hist_good_jetNTracks" , "Jet N_{TRACKS}", 20, 0, 20 ) ;
306  hists["hist_good_jetNTracksVsPt" ] = theDir.make<TH2D>( "hist_good_jetNTracksVsPt" , "Number of Tracks versus Jet p_{T};Jet p_{T}(GeV/c) ;N_{Tracks}",20, 0, 200, 20, 0, 20 ) ;
307  hists["hist_good_jetEMF" ] = theDir.make<TH1D>( "hist_good_jetEMF" , "Jet EMF", 200, 0, 1) ;
308  hists["hist_good_jetPdgID" ] = theDir.make<TH1D>( "hist_good_jetPdgID" , "PDG Id of Jet Constituents", 10000, 0, 10000 ) ;
309  hists["hist_good_jetGenEmE" ] = theDir.make<TH1D>( "hist_good_jetGenEmE" , "Gen Jet EM Energy", 200, 0, 200 ) ;
310  hists["hist_good_jetGenHadE" ] = theDir.make<TH1D>( "hist_good_jetGenHadE" , "Gen Jet HAD Energy", 200, 0, 200 ) ;
311  hists["hist_good_jetGenEMF" ] = theDir.make<TH1D>( "hist_good_jetGenEMF" , "Gen Jet EMF", 200, 0, 1) ;
312  hists["hist_good_jetEoverGenE" ] = theDir.make<TH1D>( "hist_good_jetEoverGenE" , "Energy of reco Jet / Energy of gen Jet", 200, 0, 2.0) ;
313  hists["hist_good_jetCorr" ] = theDir.make<TH1D>( "hist_good_jetCorr" , "Jet Correction Factor", 200, 0, 1.0 ) ;
314  hists["hist_good_n90Hits" ] = theDir.make<TH1D>( "hist_good_n90Hits" , "Jet n90Hits", 20, 0, 20) ;
315  hists["hist_good_fHPD" ] = theDir.make<TH1D>( "hist_good_fHPD" , "Jet fHPD", 200, 0, 1) ;
316  hists["hist_good_nConstituents" ] = theDir.make<TH1D>( "hist_good_nConstituents" , "Jet nConstituents", 20, 0, 20 ) ;
317  hists["hist_good_jetCHF" ] = theDir.make<TH1D>( "hist_good_jetCHF" , "Jet Charged Hadron Fraction", 200, 0, 1.0) ;
318 
319 
320  hists["hist_pf_jetPt" ] = theDir.make<TH1D>( "hist_pf_jetPt" , "PFJet p_{T}", 400, 0, 400 ) ;
321  hists["hist_pf_jetEtaVsPhi" ] = theDir.make<TH2D>( "hist_pf_jetEtaVsPhi" , "PFJet #phi versus #eta;#eta;#phi", 50, -5.0, 5.0, 50, -TMath::Pi(), TMath::Pi() ) ;
322  hists["hist_pf_jetNTracks" ] = theDir.make<TH1D>( "hist_pf_jetNTracks" , "PFJet N_{TRACKS}", 20, 0, 20 ) ;
323  hists["hist_pf_jetNTracksVsPt" ] = theDir.make<TH2D>( "hist_pf_jetNTracksVsPt" , "Number of Tracks versus Jet p_{T};Jet p_{T}(GeV/c) ;N_{Tracks}",20, 0, 200, 20, 0, 20 ) ;
324  hists["hist_pf_jetEMF" ] = theDir.make<TH1D>( "hist_pf_jetEMF" , "PFJet EMF", 200, 0, 1) ;
325  hists["hist_pf_jetCHF" ] = theDir.make<TH1D>( "hist_pf_jetCHF" , "PFJet CHF", 200, 0, 1) ;
326  hists["hist_pf_jetNHF" ] = theDir.make<TH1D>( "hist_pf_jetNHF" , "PFJet NHF", 200, 0, 1) ;
327  hists["hist_pf_jetCEF" ] = theDir.make<TH1D>( "hist_pf_jetCEF" , "PFJet CEF", 200, 0, 1) ;
328  hists["hist_pf_jetNEF" ] = theDir.make<TH1D>( "hist_pf_jetNEF" , "PFJet NEF", 200, 0, 1) ;
329  hists["hist_pf_jetPdgID" ] = theDir.make<TH1D>( "hist_pf_jetPdgID" , "PDG Id of Jet Constituents", 10000, 0, 10000 ) ;
330  hists["hist_pf_jetGenEmE" ] = theDir.make<TH1D>( "hist_pf_jetGenEmE" , "Gen Jet EM Energy", 200, 0, 200 ) ;
331  hists["hist_pf_jetGenHadE" ] = theDir.make<TH1D>( "hist_pf_jetGenHadE" , "Gen Jet HAD Energy", 200, 0, 200 ) ;
332  hists["hist_pf_jetGenEMF" ] = theDir.make<TH1D>( "hist_pf_jetGenEMF" , "Gen Jet EMF", 200, 0, 1) ;
333  hists["hist_pf_jetEoverGenE" ] = theDir.make<TH1D>( "hist_pf_jetEoverGenE" , "Energy of reco Jet / Energy of gen Jet", 200, 0, 2.0) ;
334  hists["hist_pf_jetCorr" ] = theDir.make<TH1D>( "hist_pf_jetCorr" , "PFJet Correction Factor", 200, 0, 1.0 ) ;
335  hists["hist_pf_nConstituents" ] = theDir.make<TH1D>( "hist_pf_nConstituents" , "PFJet nConstituents", 20, 0, 20 ) ;
336 
337  hists["hist_pf_good_jetPt" ] = theDir.make<TH1D>( "hist_pf_good_jetPt" , "PFJet p_{T}", 400, 0, 400 ) ;
338  hists["hist_pf_good_jetEtaVsPhi" ] = theDir.make<TH2D>( "hist_pf_good_jetEtaVsPhi" , "PFJet #phi versus #eta;#eta;#phi", 50, -5.0, 5.0, 50, -TMath::Pi(), TMath::Pi() ) ;
339  hists["hist_pf_good_jetNTracks" ] = theDir.make<TH1D>( "hist_pf_good_jetNTracks" , "PFJet N_{TRACKS}", 20, 0, 20 ) ;
340  hists["hist_pf_good_jetNTracksVsPt" ] = theDir.make<TH2D>( "hist_pf_good_jetNTracksVsPt", "Number of Tracks versus Jet p_{T};Jet p_{T}(GeV/c) ;N_{Tracks}",20, 0, 200, 20, 0, 20 ) ;
341  hists["hist_pf_good_jetEMF" ] = theDir.make<TH1D>( "hist_pf_good_jetEMF" , "PFJet EMF", 200, 0, 1) ;
342  hists["hist_pf_good_jetCHF" ] = theDir.make<TH1D>( "hist_pf_good_jetCHF" , "PFJet CHF", 200, 0, 1) ;
343  hists["hist_pf_good_jetNHF" ] = theDir.make<TH1D>( "hist_pf_good_jetNHF" , "PFJet NHF", 200, 0, 1) ;
344  hists["hist_pf_good_jetCEF" ] = theDir.make<TH1D>( "hist_pf_good_jetCEF" , "PFJet CEF", 200, 0, 1) ;
345  hists["hist_pf_good_jetNEF" ] = theDir.make<TH1D>( "hist_pf_good_jetNEF" , "PFJet NEF", 200, 0, 1) ;
346  hists["hist_pf_good_jetPdgID" ] = theDir.make<TH1D>( "hist_pf_good_jetPdgID" , "PDG Id of Jet Constituents", 10000, 0, 10000 ) ;
347  hists["hist_pf_good_jetGenEmE" ] = theDir.make<TH1D>( "hist_pf_good_jetGenEmE" , "Gen Jet EM Energy", 200, 0, 200 ) ;
348  hists["hist_pf_good_jetGenHadE" ] = theDir.make<TH1D>( "hist_pf_good_jetGenHadE" , "Gen Jet HAD Energy", 200, 0, 200 ) ;
349  hists["hist_pf_good_jetGenEMF" ] = theDir.make<TH1D>( "hist_pf_good_jetGenEMF" , "Gen Jet EMF", 200, 0, 1) ;
350  hists["hist_pf_good_jetEoverGenE" ] = theDir.make<TH1D>( "hist_pf_good_jetEoverGenE" , "Energy of reco Jet / Energy of gen Jet", 200, 0, 2.0) ;
351  hists["hist_pf_good_jetCorr" ] = theDir.make<TH1D>( "hist_pf_good_jetCorr" , "PFJet Correction Factor", 200, 0, 1.0 ) ;
352  hists["hist_pf_good_nConstituents" ] = theDir.make<TH1D>( "hist_pf_good_nConstituents" , "PFJet nConstituents", 20, 0, 20 ) ;
353 
354  hists["hist_mjj" ] = theDir.make<TH1D>( "hist_mjj" , "Dijet mass", 300, 0, 300 ) ;
355  hists["hist_dR_jj" ] = theDir.make<TH1D>( "hist_dR_jj" , "#Delta R_{JJ}", 200, 0, TMath::TwoPi() ) ;
356  hists["hist_imbalance_jj" ] = theDir.make<TH1D>( "hist_imbalance_jj" , "Dijet imbalance", 200, -1.0, 1.0 ) ;
357 
358  hists["hist_pf_mjj" ] = theDir.make<TH1D>( "hist_pf_mjj" , "Dijet mass", 300, 0, 300 ) ;
359  hists["hist_pf_dR_jj" ] = theDir.make<TH1D>( "hist_pf_dR_jj" , "#Delta R_{JJ}", 200, 0, TMath::TwoPi() ) ;
360  hists["hist_pf_imbalance_jj" ] = theDir.make<TH1D>( "hist_pf_imbalance_jj" , "Dijet imbalance", 200, -1.0, 1.0 ) ;
361 
362 
363 
364  cout << "Making functors" << endl;
365  JetIDStudiesSelector caloSelector( caloJetIDParameters,
366  pfJetIDParameters,
367  jetStudiesParams );
368 
369  JetIDStudiesSelector pfSelector( caloJetIDParameters,
370  pfJetIDParameters,
371  pfJetStudiesParams );
372 
373  bool doTracks = plotParameters.getParameter<bool>("doTracks");
374  bool useMC = plotParameters.getParameter<bool>("useMC");
375 
376  cout << "About to begin looping" << endl;
377 
378  int nev = 0;
379  //loop through each event
380  for (ev.toBegin(); ! ev.atEnd(); ++ev, ++nev) {
381 
382 
383  edm::EventBase const & event = ev;
384 
385  if ( runLumiSel(ev) == false ) continue;
386 
387  if ( nev % 100 == 0 ) cout << "Processing run " << event.id().run() << ", lumi " << event.id().luminosityBlock() << ", event " << event.id().event() << endl;
388 
389 
390 
391 
392  pat::strbitset retCalo = caloSelector.getBitTemplate();
393  caloSelector( event, retCalo );
394 
395 
396  pat::strbitset retPF = pfSelector.getBitTemplate();
397  pfSelector( event, retPF );
398 
402  if ( retCalo.test( caloSelector.caloKin() ) ) {
403  if ( retCalo.test( caloSelector.caloDeltaPhi() ) ) {
404  vector<pat::Jet> const & allCaloJets = caloSelector.allCaloJets();
405 
406  for ( std::vector<pat::Jet>::const_iterator jetBegin = allCaloJets.begin(),
407  jetEnd = jetBegin + 2, ijet = jetBegin;
408  ijet != jetEnd; ++ijet ) {
409 
410  const pat::Jet & jet = *ijet;
411 
412  double pt = jet.pt();
413 
415 
416  hists["hist_jetPt"]->Fill( pt );
417  hists["hist_jetEtaVsPhi"]->Fill( jet.eta(), jet.phi() );
418  hists["hist_jetNTracks"]->Fill( jetTracks.size() );
419  hists["hist_jetNTracksVsPt"]->Fill( pt, jetTracks.size() );
420  hists["hist_jetEMF"]->Fill( jet.emEnergyFraction() );
421  hists["hist_jetCorr"]->Fill( jet.jecFactor("Uncorrected") );
422  hists["hist_n90Hits"]->Fill( static_cast<int>(jet.jetID().n90Hits) );
423  hists["hist_fHPD"]->Fill( jet.jetID().fHPD );
424  hists["hist_nConstituents"]->Fill( jet.nConstituents() );
425 
426  if ( useMC && jet.genJet() != nullptr ) {
427  hists["hist_jetGenEmE"]->Fill( jet.genJet()->emEnergy() );
428  hists["hist_jetGenHadE"]->Fill( jet.genJet()->hadEnergy() );
429  hists["hist_jetEoverGenE"]->Fill( jet.energy() / jet.genJet()->energy() );
430  hists["hist_jetGenEMF"]->Fill( jet.genJet()->emEnergy() / jet.genJet()->energy() );
431  }
432  if ( doTracks ) {
433  TLorentzVector p4_tracks(0,0,0,0);
434  for ( reco::TrackRefVector::const_iterator itrk = jetTracks.begin(),
435  itrkEnd = jetTracks.end();
436  itrk != itrkEnd; ++itrk ) {
437  TLorentzVector p4_trk;
438  double M_PION = 0.140;
439  p4_trk.SetPtEtaPhiM( (*itrk)->pt(), (*itrk)->eta(), (*itrk)->phi(), M_PION );
440  p4_tracks += p4_trk;
441  }
442  hists["hist_jetCHF"]->Fill( p4_tracks.Energy() / jet.energy() );
443  }
444 
445 
446  } // end loop over jets
447 
448 
449 
450  if ( retCalo.test( caloSelector.caloJetID() ) ) {
451  pat::Jet const & jet0 = caloSelector.caloJet0();
452  pat::Jet const & jet1 = caloSelector.caloJet1();
453 
454  TLorentzVector p4_j0( jet0.px(), jet0.py(), jet0.pz(), jet0.energy() );
455  TLorentzVector p4_j1( jet1.px(), jet1.py(), jet1.pz(), jet1.energy() );
456 
457  TLorentzVector p4_jj = p4_j0 + p4_j1;
458 
459  hists["hist_mjj"]->Fill( p4_jj.M() );
460  hists["hist_dR_jj"]->Fill( p4_j0.DeltaR( p4_j1 ) );
461  hists["hist_imbalance_jj"]->Fill( (p4_j0.Perp() - p4_j1.Perp() ) /
462  (p4_j0.Perp() + p4_j1.Perp() ) );
463 
464  hists["hist_good_jetPt"]->Fill( jet0.pt() );
465  hists["hist_good_jetEtaVsPhi"]->Fill( jet0.eta(), jet0.phi() );
466  hists["hist_good_jetNTracks"]->Fill( jet0.associatedTracks().size() );
467  hists["hist_good_jetNTracksVsPt"]->Fill( jet0.pt(), jet0.associatedTracks().size() );
468  hists["hist_good_jetEMF"]->Fill( jet0.emEnergyFraction() );
469  hists["hist_good_jetCorr"]->Fill( jet0.jecFactor("Uncorrected") );
470  hists["hist_good_n90Hits"]->Fill( static_cast<int>(jet0.jetID().n90Hits) );
471  hists["hist_good_fHPD"]->Fill( jet0.jetID().fHPD );
472  hists["hist_good_nConstituents"]->Fill( jet0.nConstituents() );
473 
474 
475  hists["hist_good_jetPt"]->Fill( jet1.pt() );
476  hists["hist_good_jetEtaVsPhi"]->Fill( jet1.eta(), jet1.phi() );
477  hists["hist_good_jetNTracks"]->Fill( jet1.associatedTracks().size() );
478  hists["hist_good_jetNTracksVsPt"]->Fill( jet1.pt(), jet1.associatedTracks().size() );
479  hists["hist_good_jetEMF"]->Fill( jet1.emEnergyFraction() );
480  hists["hist_good_jetCorr"]->Fill( jet1.jecFactor("Uncorrected") );
481  hists["hist_good_n90Hits"]->Fill( static_cast<int>(jet1.jetID().n90Hits) );
482  hists["hist_good_fHPD"]->Fill( jet1.jetID().fHPD );
483  hists["hist_good_nConstituents"]->Fill( jet1.nConstituents() );
484 
485  }// end if passed calo jet id
486  }// end if passed dphi cuts
487  }// end if passed kin cuts
488 
489 
493  if ( retPF.test( pfSelector.pfDeltaPhi() ) ) {
494 
495  vector<pat::Jet> const & allPFJets = pfSelector.allPFJets();
496 
497  for ( std::vector<pat::Jet>::const_iterator jetBegin = allPFJets.begin(),
498  jetEnd = jetBegin + 2, ijet = jetBegin;
499  ijet != jetEnd; ++ijet ) {
500 
501  const pat::Jet & jet = *ijet;
502 
503  double pt = jet.pt();
504 
505  hists["hist_pf_jetPt"]->Fill( pt );
506  hists["hist_pf_jetEtaVsPhi"]->Fill( jet.eta(), jet.phi() );
507  hists["hist_pf_nConstituents"]->Fill( jet.nConstituents() );
508  hists["hist_pf_jetCEF"]->Fill( jet.chargedEmEnergyFraction() );
509  hists["hist_pf_jetNEF"]->Fill( jet.neutralEmEnergyFraction() );
510  hists["hist_pf_jetCHF"]->Fill( jet.chargedHadronEnergyFraction() );
511  hists["hist_pf_jetNHF"]->Fill( jet.neutralHadronEnergyFraction() );
512 
513 
514  if ( useMC && jet.genJet() != nullptr ) {
515  hists["hist_pf_jetGenEmE"]->Fill( jet.genJet()->emEnergy() );
516  hists["hist_pf_jetGenHadE"]->Fill( jet.genJet()->hadEnergy() );
517  hists["hist_pf_jetEoverGenE"]->Fill( jet.energy() / jet.genJet()->energy() );
518 
519  hists["hist_pf_jetGenEMF"]->Fill( jet.genJet()->emEnergy() / jet.genJet()->energy() );
520  }
521 
522  } // end loop over jets
523 
524 
525 
526  if ( retPF.test( pfSelector.pfJetID() ) ) {
527  pat::Jet const & jet0 = pfSelector.pfJet0();
528  pat::Jet const & jet1 = pfSelector.pfJet1();
529 
530  TLorentzVector p4_j0( jet0.px(), jet0.py(), jet0.pz(), jet0.energy() );
531  TLorentzVector p4_j1( jet1.px(), jet1.py(), jet1.pz(), jet1.energy() );
532 
533  TLorentzVector p4_jj = p4_j0 + p4_j1;
534 
535  hists["hist_pf_mjj"]->Fill( p4_jj.M() );
536  hists["hist_pf_dR_jj"]->Fill( p4_j0.DeltaR( p4_j1 ) );
537  hists["hist_pf_imbalance_jj"]->Fill( (p4_j0.Perp() - p4_j1.Perp() ) /
538  (p4_j0.Perp() + p4_j1.Perp() ) );
539 
540  hists["hist_pf_good_jetPt"]->Fill( jet0.pt() );
541  hists["hist_pf_good_jetEtaVsPhi"]->Fill( jet0.eta(), jet0.phi() );
542  hists["hist_pf_good_nConstituents"]->Fill( jet0.nConstituents() );
543  hists["hist_pf_good_jetCEF"]->Fill( jet0.chargedEmEnergyFraction() );
544  hists["hist_pf_good_jetNEF"]->Fill( jet0.neutralEmEnergyFraction() );
545  hists["hist_pf_good_jetCHF"]->Fill( jet0.chargedHadronEnergyFraction() );
546  hists["hist_pf_good_jetNHF"]->Fill( jet0.neutralHadronEnergyFraction() );
547 
548 
549  hists["hist_pf_good_jetPt"]->Fill( jet1.pt() );
550  hists["hist_pf_good_jetEtaVsPhi"]->Fill( jet1.eta(), jet1.phi() );
551  hists["hist_pf_good_nConstituents"]->Fill( jet1.nConstituents() );
552  hists["hist_pf_good_jetCEF"]->Fill( jet1.chargedEmEnergyFraction() );
553  hists["hist_pf_good_jetNEF"]->Fill( jet1.neutralEmEnergyFraction() );
554  hists["hist_pf_good_jetCHF"]->Fill( jet1.chargedHadronEnergyFraction() );
555  hists["hist_pf_good_jetNHF"]->Fill( jet1.neutralHadronEnergyFraction() );
556 
557  } // end if 2 good PF jets
558 
559  }// end if delta phi pf cuts
560 
561  } // end loop over events
562 
563  cout << "Calo jet selection" << endl;
564  caloSelector.print(std::cout);
565  cout << "PF jet selection" << endl;
566  pfSelector.print(std::cout);
567 
568 
569 
570 
571  return 0;
572 }
573 
ChainEvent const & toBegin() override
Definition: ChainEvent.cc:198
const double TwoPi
const double Pi
T getParameter(std::string const &) const
float neutralHadronEnergyFraction() const
neutralHadronEnergyFraction (relative to uncorrected jet energy)
Definition: Jet.h:374
float hadEnergy() const
Definition: GenJet.h:59
double eta() const final
momentum pseudorapidity
float emEnergy() const
Definition: GenJet.h:57
Selector< edm::EventBase > EventSelector
Definition: EventSelector.h:20
float fHPD
Definition: JetID.h:45
float chargedHadronEnergyFraction() const
chargedHadronEnergyFraction (relative to uncorrected jet energy)
Definition: Jet.h:372
index_type const & pfDeltaPhi() const
const reco::TrackRefVector & associatedTracks() const
method to return a vector of refs to the tracks associated to this jet
double px() const final
x coordinate of momentum vector
float chargedEmEnergyFraction() const
chargedEmEnergyFraction (relative to uncorrected jet energy)
Definition: Jet.h:376
int main(int argc, char *argv[])
double pt() const final
transverse momentum
bool ev
#define nullptr
const_iterator end() const
Termination of iteration.
Definition: RefVector.h:253
const_iterator begin() const
Initialize an iterator over the RefVector.
Definition: RefVector.h:248
double pz() const final
z coordinate of momentum vector
static void enable()
enable automatic library loading
std::shared_ptr< JetIDSelectionFunctor > const & jetSel() const
double energy() const final
energy
index_type const & caloJetID() const
bool atEnd() const override
Definition: ChainEvent.cc:301
T * make(const Args &...args) const
make new ROOT object
PF Jet selector for pat::Jets.
float emEnergyFraction() const
returns the jet electromagnetic energy fraction
Definition: Jet.h:300
edm::Handle< vector< pat::Jet > > h_jets_
index_type const & caloDeltaPhi() const
bool operator()(edm::EventBase const &event, pat::strbitset &ret) override
Jet selector for pat::Jets and for CaloJets.
std::shared_ptr< JetIDSelectionFunctor > jetSel_
reco::JetID const & jetID() const
accessing Jet ID information
Definition: Jet.h:456
vector< pat::Jet > const & allPFJets() const
TFileDirectory mkdir(const std::string &dir, const std::string &descr="")
create a new subdirectory
strbitset & set(bool val=true)
set method of all bits
Definition: strbitset.h:144
std::shared_ptr< PFJetIDSelectionFunctor > const & pfJetSel() const
double b
Definition: hdecay.h:120
double py() const final
y coordinate of momentum vector
Analysis-level calorimeter jet class.
Definition: Jet.h:80
std::shared_ptr< PFJetIDSelectionFunctor > pfJetSel_
bool test(std::string s) const
test
Definition: strbitset.h:331
vector< pat::Jet > const & allCaloJets() const
index_type const & caloCuts() const
HLT enums.
virtual int nConstituents() const
of constituents
Definition: Jet.h:65
const reco::GenJet * genJet() const
return the matched generated jet
edm::Handle< vector< pat::Jet > > h_pfjets_
size_type size() const
Size of the RefVector.
Definition: RefVector.h:107
index_type const & pfJetID() const
float jecFactor(const std::string &level, const std::string &flavor="none", const std::string &set="") const
float neutralEmEnergyFraction() const
neutralEmEnergyFraction (relative to uncorrected jet energy)
Definition: Jet.h:378
JetIDStudiesSelector(edm::ParameterSet const &caloJetIdParams, edm::ParameterSet const &pfJetIdParams, edm::ParameterSet const &params)
double phi() const final
momentum azimuthal angle
std::shared_ptr< edm::ProcessDesc > processDesc() const
index_type const & caloKin() const
index_type const & pfCuts() const
short n90Hits
Definition: JetID.h:47
Definition: event.py:1