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