CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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 
101 
102  virtual 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() != 0 ) {
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() != 0 ) {
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 
const double TwoPi
const double Pi
T getParameter(std::string const &) const
dictionary parameters
Definition: Parameters.py:2
float neutralHadronEnergyFraction() const
neutralHadronEnergyFraction (relative to uncorrected jet energy)
Definition: Jet.h:368
float hadEnergy() const
Definition: GenJet.h:59
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:366
index_type const & pfDeltaPhi() const
const reco::TrackRefVector & associatedTracks() const
method to return a vector of refs to the tracks associated to this jet
float chargedEmEnergyFraction() const
chargedEmEnergyFraction (relative to uncorrected jet energy)
Definition: Jet.h:370
bool ev
const_iterator end() const
Termination of iteration.
Definition: RefVector.h:255
std::shared_ptr< edm::ProcessDesc > processDesc()
virtual double eta() const
momentum pseudorapidity
virtual double pt() const
transverse momentum
const_iterator begin() const
Initialize an iterator over the RefVector.
Definition: RefVector.h:250
virtual double energy() const
energy
static void enable()
enable automatic library loading
std::shared_ptr< JetIDSelectionFunctor > const & jetSel() const
index_type const & caloJetID() const
T * make(const Args &...args) const
make new ROOT object
How EventSelector::AcceptEvent() decides whether to accept an event for output otherwise it is excluding the probing of A single or multiple positive and the trigger will pass if any such matching triggers are PASS or EXCEPTION[A criterion thatmatches no triggers at all is detected and causes a throw.] A single negative with an expectation of appropriate bit checking in the decision and the trigger will pass if any such matching triggers are FAIL or EXCEPTION A wildcarded negative criterion that matches more than one trigger in the trigger but the state exists so we define the behavior If all triggers are the negative crieriion will lead to accepting the event(this again matches the behavior of"!*"before the partial wildcard feature was incorporated).The per-event"cost"of each negative criterion with multiple relevant triggers is about the same as!*was in the past
PF Jet selector for pat::Jets.
float emEnergyFraction() const
returns the jet electromagnetic energy fraction
Definition: Jet.h:294
virtual bool atEnd() const
Definition: ChainEvent.cc:303
edm::Handle< vector< pat::Jet > > h_jets_
index_type const & caloDeltaPhi() const
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:447
vector< pat::Jet > const & allPFJets() const
virtual double px() const
x coordinate of momentum vector
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
tuple argc
Definition: dir2webdir.py:38
double b
Definition: hdecay.h:120
Analysis-level calorimeter jet class.
Definition: Jet.h:77
std::shared_ptr< PFJetIDSelectionFunctor > pfJetSel_
virtual double pz() const
z coordinate of momentum vector
bool test(std::string s) const
test
Definition: strbitset.h:331
ChainEvent const & toBegin()
Definition: ChainEvent.cc:200
vector< pat::Jet > const & allCaloJets() const
index_type const & caloCuts() const
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:99
tuple cout
Definition: gather_cfg.py:121
index_type const & pfJetID() const
virtual bool operator()(edm::EventBase const &event, pat::strbitset &ret) override
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:372
volatile std::atomic< bool > shutdown_flag false
JetIDStudiesSelector(edm::ParameterSet const &caloJetIdParams, edm::ParameterSet const &pfJetIdParams, edm::ParameterSet const &params)
virtual double phi() const
momentum azimuthal angle
index_type const & caloKin() const
index_type const & pfCuts() const
virtual double py() const
y coordinate of momentum vector
short n90Hits
Definition: JetID.h:47