CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
DiJetAnalyzer.cc
Go to the documentation of this file.
1 //
2 // DiJetAnalyzer.cc
3 //
4 // description: Create an ntuple necessary for dijet balance calibration for the HCAL
5 //
6 
8 
10 {
11  // set parameters
12  pfJetCollName_ = iConfig.getParameter<std::string>("pfJetCollName");
13  pfJetCorrName_ = iConfig.getParameter<std::string>("pfJetCorrName");
14  hbheRecHitName_ = iConfig.getParameter<std::string>("hbheRecHitName");
15  hfRecHitName_ = iConfig.getParameter<std::string>("hfRecHitName");
16  hoRecHitName_ = iConfig.getParameter<std::string>("hoRecHitName");
17  pvCollName_ = iConfig.getParameter<std::string>("pvCollName");
18  rootHistFilename_ = iConfig.getParameter<std::string>("rootHistFilename");
19  maxDeltaEta_ = iConfig.getParameter<double>("maxDeltaEta");
20  minTagJetEta_ = iConfig.getParameter<double>("minTagJetEta");
21  maxTagJetEta_ = iConfig.getParameter<double>("maxTagJetEta");
22  minSumJetEt_ = iConfig.getParameter<double>("minSumJetEt");
23  minJetEt_ = iConfig.getParameter<double>("minJetEt");
24  maxThirdJetEt_ = iConfig.getParameter<double>("maxThirdJetEt");
25  debug_ = iConfig.getUntrackedParameter<bool>("debug", false);
26 
27  tok_PFJet_ = consumes<reco::PFJetCollection>(pfJetCollName_);
28  tok_HBHE_ = consumes<edm::SortedCollection<HBHERecHit,edm::StrictWeakOrdering<HBHERecHit> > >(hbheRecHitName_);
29  tok_HF_ = consumes<edm::SortedCollection<HFRecHit,edm::StrictWeakOrdering<HFRecHit> > >(hfRecHitName_);
30  tok_HO_ = consumes<edm::SortedCollection<HORecHit,edm::StrictWeakOrdering<HORecHit> > >(hoRecHitName_);
31  tok_Vertex_ = consumes<reco::VertexCollection>(pvCollName_);
32 }
33 
35 {
36 }
37 
38 //
39 // member functions
40 //
41 
42 // ------------ method called to for each event ------------
43 void
45 {
46  pf_Run_ = iEvent.id().run();
47  pf_Lumi_ = iEvent.id().luminosityBlock();
48  pf_Event_ = iEvent.id().event();
49 
50  // Get PFJets
52  iEvent.getByToken(tok_PFJet_,pfjets);
53  if(!pfjets.isValid()) {
55  << " could not find PFJetCollection named " << pfJetCollName_ << ".\n";
56  return;
57  }
58 
59  // Get RecHits in HB and HE
61  iEvent.getByToken(tok_HBHE_,hbhereco);
62  if(!hbhereco.isValid()) {
64  << " could not find HBHERecHit named " << hbheRecHitName_ << ".\n";
65  return;
66  }
67 
68  // Get RecHits in HF
70  iEvent.getByToken(tok_HF_,hfreco);
71  if(!hfreco.isValid()) {
73  << " could not find HFRecHit named " << hfRecHitName_ << ".\n";
74  return;
75  }
76 
77  // Get RecHits in HO
79  iEvent.getByToken(tok_HO_,horeco);
80  if(!horeco.isValid()) {
82  << " could not find HORecHit named " << hoRecHitName_ << ".\n";
83  return;
84  }
85 
86  // Get geometry
88  evSetup.get<CaloGeometryRecord>().get(geoHandle);
89  const CaloSubdetectorGeometry *HBGeom = geoHandle->getSubdetectorGeometry(DetId::Hcal, 1);
90  const CaloSubdetectorGeometry *HEGeom = geoHandle->getSubdetectorGeometry(DetId::Hcal, 2);
91  const CaloSubdetectorGeometry *HOGeom = geoHandle->getSubdetectorGeometry(DetId::Hcal, 3);
92  const CaloSubdetectorGeometry *HFGeom = geoHandle->getSubdetectorGeometry(DetId::Hcal, 4);
93 
94  int HBHE_n = 0;
95  for(edm::SortedCollection<HBHERecHit,edm::StrictWeakOrdering<HBHERecHit>>::const_iterator ith=hbhereco->begin(); ith!=hbhereco->end(); ++ith){
96  HBHE_n++;
97  }
98  int HF_n = 0;
99  for(edm::SortedCollection<HFRecHit,edm::StrictWeakOrdering<HFRecHit>>::const_iterator ith=hfreco->begin(); ith!=hfreco->end(); ++ith){
100  HF_n++;
101  }
102  int HO_n = 0;
103  for(edm::SortedCollection<HORecHit,edm::StrictWeakOrdering<HORecHit>>::const_iterator ith=horeco->begin(); ith!=horeco->end(); ++ith){
104  HO_n++;
105  }
106 
107  // Get primary vertices
109  iEvent.getByToken(tok_Vertex_,pv);
110  if(!pv.isValid()) {
112  << " could not find Vertex named " << pvCollName_ << ".\n";
113  return;
114  }
115  pf_NPV_ = 0;
116  for(std::vector<reco::Vertex>::const_iterator it=pv->begin(); it!=pv->end(); ++it){
117  if(!it->isFake() && it->ndof() > 4) ++pf_NPV_;
118  }
119 
120  // Get jet corrections
121  const JetCorrector* correctorPF = JetCorrector::getJetCorrector(pfJetCorrName_,evSetup);
122 
124  // Event Selection
126 
127  // determine which cut results in failure
128  int passSelPF=0;
129 
130  // sort jets by corrected et
131  std::set<JetCorretPair, JetCorretPairComp> pfjetcorretpairset;
132  for(reco::PFJetCollection::const_iterator it=pfjets->begin(); it!=pfjets->end(); ++it) {
133  const reco::PFJet* jet=&(*it);
134  double jec = correctorPF->correction(*it, iEvent, evSetup);
135  pfjetcorretpairset.insert(JetCorretPair(jet, jec));
136  }
137 
138  JetCorretPair pf_tag, pf_probe;
142  int cntr=0;
143  for(std::set<JetCorretPair, JetCorretPairComp>::const_iterator it=pfjetcorretpairset.begin(); it!=pfjetcorretpairset.end(); ++it) {
144  JetCorretPair jet=(*it);
145  ++cntr;
146  if(cntr==1) pf_tag=jet;
147  else if(cntr==2) pf_probe=jet;
148  else {
149  pf_thirdjet_px_ += jet.scale()*jet.jet()->px();
150  pf_thirdjet_py_ += jet.scale()*jet.jet()->py();
151  if(cntr==3){
152  pf_realthirdjet_px_ = jet.jet()->px();
153  pf_realthirdjet_py_ = jet.jet()->py();
155  }
156  }
157  }
158 
159  if(pf_tag.jet() && pf_probe.jet()){
160  // require that the first two jets are above some minimum,
161  // and the rest are below some maximum
162  if((pf_tag.jet()->et()+pf_probe.jet()->et())<minSumJetEt_) passSelPF |= 0x1;
163  if(pf_tag.jet()->et()<minJetEt_ || pf_probe.jet()->et()<minJetEt_) passSelPF |= 0x2;
165 
166  // force the tag jet to have the smaller |eta|
167  if(std::fabs(pf_tag.jet()->eta())>std::fabs(pf_probe.jet()->eta())) {
168  JetCorretPair temp=pf_tag;
169  pf_tag=pf_probe;
170  pf_probe=temp;
171  }
172 
173  // eta cuts
174  double dAbsEta=std::fabs(std::fabs(pf_tag.jet()->eta())-std::fabs(pf_probe.jet()->eta()));
175  if(dAbsEta>maxDeltaEta_) passSelPF |= 0x8;
176  if(fabs(pf_tag.jet()->eta())<minTagJetEta_) passSelPF |= 0x10;
177  if(fabs(pf_tag.jet()->eta())>maxTagJetEta_) passSelPF |= 0x10;
178  }
179  else{
180  passSelPF = 0x40;
181  }
182 
183  h_PassSelPF_->Fill(passSelPF);
184  if(passSelPF) return;
185  // dump
186  if(debug_) {
187  std::cout << "Run: " << iEvent.id().run() << "; Event: " << iEvent.id().event() << std::endl;
188  for(reco::PFJetCollection::const_iterator it=pfjets->begin(); it!=pfjets->end(); ++it) {
189  const reco::PFJet *jet=&(*it);
190  std::cout << "istag=" << (jet==pf_tag.jet()) << "; isprobe=" << (jet==pf_probe.jet()) << "; et=" << jet->et() << "; eta=" << jet->eta() << std::endl;
191  }
192  }
193 
194  // Reset particle variables
200  tpfjet_had_n_ = 0;
201  tpfjet_cluster_n_ = 0;
207  ppfjet_had_n_ = 0;
208  ppfjet_cluster_n_ = 0;
209 
210  tpfjet_had_E_.clear();
211  tpfjet_had_px_.clear();
212  tpfjet_had_py_.clear();
213  tpfjet_had_pz_.clear();
214  tpfjet_had_EcalE_.clear();
215  tpfjet_had_rawHcalE_.clear();
216  tpfjet_had_emf_.clear();
217  tpfjet_had_id_.clear();
218  tpfjet_had_candtrackind_.clear();
219  tpfjet_had_ntwrs_.clear();
220  tpfjet_twr_ieta_.clear();
221  tpfjet_twr_iphi_.clear();
222  tpfjet_twr_depth_.clear();
223  tpfjet_twr_subdet_.clear();
224  tpfjet_twr_candtrackind_.clear();
225  tpfjet_twr_hadind_.clear();
226  tpfjet_twr_elmttype_.clear();
227  tpfjet_twr_hade_.clear();
228  tpfjet_twr_frac_.clear();
229  tpfjet_twr_dR_.clear();
230  tpfjet_twr_clusterind_.clear();
231  tpfjet_cluster_eta_.clear();
232  tpfjet_cluster_phi_.clear();
233  tpfjet_cluster_dR_.clear();
234  tpfjet_candtrack_px_.clear();
235  tpfjet_candtrack_py_.clear();
236  tpfjet_candtrack_pz_.clear();
237  tpfjet_candtrack_EcalE_.clear();
238  ppfjet_had_E_.clear();
239  ppfjet_had_px_.clear();
240  ppfjet_had_py_.clear();
241  ppfjet_had_pz_.clear();
242  ppfjet_had_EcalE_.clear();
243  ppfjet_had_rawHcalE_.clear();
244  ppfjet_had_emf_.clear();
245  ppfjet_had_id_.clear();
246  ppfjet_had_candtrackind_.clear();
247  ppfjet_had_ntwrs_.clear();
248  ppfjet_twr_ieta_.clear();
249  ppfjet_twr_iphi_.clear();
250  ppfjet_twr_depth_.clear();
251  ppfjet_twr_subdet_.clear();
252  ppfjet_twr_candtrackind_.clear();
253  ppfjet_twr_hadind_.clear();
254  ppfjet_twr_elmttype_.clear();
255  ppfjet_twr_hade_.clear();
256  ppfjet_twr_frac_.clear();
257  ppfjet_twr_dR_.clear();
258  ppfjet_twr_clusterind_.clear();
259  ppfjet_cluster_eta_.clear();
260  ppfjet_cluster_phi_.clear();
261  ppfjet_cluster_dR_.clear();
262  ppfjet_candtrack_px_.clear();
263  ppfjet_candtrack_py_.clear();
264  ppfjet_candtrack_pz_.clear();
265  ppfjet_candtrack_EcalE_.clear();
266 
267  std::map<int,std::pair<int,std::set<float>>> tpfjet_rechits;
268  std::map<int,std::pair<int,std::set<float>>> ppfjet_rechits;
269  std::map<float,int> tpfjet_clusters;
270  std::map<float,int> ppfjet_clusters;
271 
272  // fill tag jet variables
273  tpfjet_pt_ = pf_tag.jet()->pt();
274  tpfjet_p_ = pf_tag.jet()->p();
275  tpfjet_E_ = pf_tag.jet()->energy();
276  tpfjet_eta_ = pf_tag.jet()->eta();
277  tpfjet_phi_ = pf_tag.jet()->phi();
278  tpfjet_scale_ = pf_tag.scale();
279  tpfjet_area_ = pf_tag.jet()->jetArea();
280  tpfjet_ntwrs_=0;
282 
283  tpfjet_jetID_ = 0; // Not a loose, medium, or tight jet
284  if(fabs(pf_tag.jet()->eta()) < 2.4){
285  if(pf_tag.jet()->chargedHadronEnergyFraction() > 0 &&
286  pf_tag.jet()->chargedMultiplicity() > 0 &&
287  pf_tag.jet()->chargedEmEnergyFraction() < 0.99 &&
288  (pf_tag.jet()->chargedMultiplicity() + pf_tag.jet()->neutralMultiplicity()) > 1){
289  if(pf_tag.jet()->neutralHadronEnergyFraction() < 0.9 &&
290  pf_tag.jet()->neutralEmEnergyFraction() < 0.9){
291  tpfjet_jetID_ = 3; // Tight jet
292  }
293  else if(pf_tag.jet()->neutralHadronEnergyFraction() < 0.95 &&
294  pf_tag.jet()->neutralEmEnergyFraction() < 0.95){
295  tpfjet_jetID_ = 2; // Medium jet
296  }
297  else if(pf_tag.jet()->neutralHadronEnergyFraction() < 0.99 &&
298  pf_tag.jet()->neutralEmEnergyFraction() < 0.99){
299  tpfjet_jetID_ = 1; // Loose jet
300  }
301  }
302  }
303  else if((pf_tag.jet()->chargedMultiplicity() + pf_tag.jet()->neutralMultiplicity()) > 1){
304  if(pf_tag.jet()->neutralHadronEnergyFraction() < 0.9 &&
305  pf_tag.jet()->neutralEmEnergyFraction() < 0.9){
306  tpfjet_jetID_ = 3; // Tight jet
307  }
308  else if(pf_tag.jet()->neutralHadronEnergyFraction() < 0.95 &&
309  pf_tag.jet()->neutralEmEnergyFraction() < 0.95){
310  tpfjet_jetID_ = 2; // Medium jet
311  }
312  else if(pf_tag.jet()->neutralHadronEnergyFraction() < 0.99 &&
313  pf_tag.jet()->neutralEmEnergyFraction() < 0.99){
314  tpfjet_jetID_ = 1; // Loose jet
315  }
316  }
317 
319  // Get PF constituents and fill HCAL towers
321 
322  // Get tag PFCandidates
323  std::vector<reco::PFCandidatePtr> tagconst=pf_tag.jet()->getPFConstituents();
324  for(std::vector<reco::PFCandidatePtr>::const_iterator it=tagconst.begin(); it!=tagconst.end(); ++it){
325  bool hasTrack = false;
326  // Test PFCandidate type
327  reco::PFCandidate::ParticleType candidateType = (*it)->particleId();
328  switch(candidateType){
330  tpfjet_unkown_E_ += (*it)->energy();
331  tpfjet_unkown_px_ += (*it)->px();
332  tpfjet_unkown_py_ += (*it)->py();
333  tpfjet_unkown_pz_ += (*it)->pz();
334  tpfjet_unkown_EcalE_ += (*it)->ecalEnergy();
336  continue;
338  {
339  tpfjet_had_E_.push_back((*it)->energy());
340  tpfjet_had_px_.push_back((*it)->px());
341  tpfjet_had_py_.push_back((*it)->py());
342  tpfjet_had_pz_.push_back((*it)->pz());
343  tpfjet_had_EcalE_.push_back((*it)->ecalEnergy());
344  tpfjet_had_rawHcalE_.push_back((*it)->rawHcalEnergy());
345  tpfjet_had_id_.push_back(0);
346  tpfjet_had_ntwrs_.push_back(0);
347  tpfjet_had_n_++;
348 
349  reco::TrackRef trackRef = (*it)->trackRef();
350  if(trackRef.isNonnull()){
351  reco::Track track = *trackRef;
352  tpfjet_candtrack_px_.push_back(track.px());
353  tpfjet_candtrack_py_.push_back(track.py());
354  tpfjet_candtrack_pz_.push_back(track.pz());
355  tpfjet_candtrack_EcalE_.push_back((*it)->ecalEnergy());
357  hasTrack = true;
359  }
360  else{
361  tpfjet_had_candtrackind_.push_back(-2);
362  }
363  }
364  break;
366  tpfjet_electron_E_ += (*it)->energy();
367  tpfjet_electron_px_ += (*it)->px();
368  tpfjet_electron_py_ += (*it)->py();
369  tpfjet_electron_pz_ += (*it)->pz();
370  tpfjet_electron_EcalE_ += (*it)->ecalEnergy();
372  continue;
374  tpfjet_muon_E_ += (*it)->energy();
375  tpfjet_muon_px_ += (*it)->px();
376  tpfjet_muon_py_ += (*it)->py();
377  tpfjet_muon_pz_ += (*it)->pz();
378  tpfjet_muon_EcalE_ += (*it)->ecalEnergy();
379  tpfjet_muon_n_++;
380  continue;
382  tpfjet_photon_E_ += (*it)->energy();
383  tpfjet_photon_px_ += (*it)->px();
384  tpfjet_photon_py_ += (*it)->py();
385  tpfjet_photon_pz_ += (*it)->pz();
386  tpfjet_photon_EcalE_ += (*it)->ecalEnergy();
388  continue;
390  {
391  tpfjet_had_E_.push_back((*it)->energy());
392  tpfjet_had_px_.push_back((*it)->px());
393  tpfjet_had_py_.push_back((*it)->py());
394  tpfjet_had_pz_.push_back((*it)->pz());
395  tpfjet_had_EcalE_.push_back((*it)->ecalEnergy());
396  tpfjet_had_rawHcalE_.push_back((*it)->rawHcalEnergy());
397  tpfjet_had_id_.push_back(1);
398  tpfjet_had_candtrackind_.push_back(-1);
399  tpfjet_had_ntwrs_.push_back(0);
400  tpfjet_had_n_++;
401  break;
402  }
404  {
405  tpfjet_had_E_.push_back((*it)->energy());
406  tpfjet_had_px_.push_back((*it)->px());
407  tpfjet_had_py_.push_back((*it)->py());
408  tpfjet_had_pz_.push_back((*it)->pz());
409  tpfjet_had_EcalE_.push_back((*it)->ecalEnergy());
410  tpfjet_had_rawHcalE_.push_back((*it)->rawHcalEnergy());
411  tpfjet_had_id_.push_back(2);
412  tpfjet_had_candtrackind_.push_back(-1);
413  tpfjet_had_ntwrs_.push_back(0);
414  tpfjet_had_n_++;
415  break;
416  }
418  {
419  tpfjet_had_E_.push_back((*it)->energy());
420  tpfjet_had_px_.push_back((*it)->px());
421  tpfjet_had_py_.push_back((*it)->py());
422  tpfjet_had_pz_.push_back((*it)->pz());
423  tpfjet_had_EcalE_.push_back((*it)->ecalEnergy());
424  tpfjet_had_rawHcalE_.push_back((*it)->rawHcalEnergy());
425  tpfjet_had_id_.push_back(3);
426  tpfjet_had_candtrackind_.push_back(-1);
427  tpfjet_had_ntwrs_.push_back(0);
428  tpfjet_had_n_++;
429  break;
430  }
431  }
432 
433  std::map<int,int> twrietas;
434  float HFHAD_E = 0;
435  float HFEM_E = 0;
436  int maxElement=(*it)->elementsInBlocks().size();
437  for(int e=0; e<maxElement; ++e){
438  // Get elements from block
439  reco::PFBlockRef blockRef = (*it)->elementsInBlocks()[e].first;
440  const edm::OwnVector<reco::PFBlockElement>& elements = blockRef->elements();
441  for(unsigned iEle=0; iEle<elements.size(); iEle++) {
442  if(elements[iEle].index() == (*it)->elementsInBlocks()[e].second){
443  if(elements[iEle].type() == reco::PFBlockElement::HCAL){ // Element is HB or HE
444  // Get cluster and hits
445  reco::PFClusterRef clusterref = elements[iEle].clusterRef();
446  reco::PFCluster cluster = *clusterref;
447  double cluster_dR = deltaR(tpfjet_eta_,tpfjet_phi_,cluster.eta(),cluster.phi());
448  if(tpfjet_clusters.count(cluster_dR) == 0){
449  tpfjet_clusters[cluster_dR] = tpfjet_cluster_n_;
450  tpfjet_cluster_eta_.push_back(cluster.eta());
451  tpfjet_cluster_phi_.push_back(cluster.phi());
452  tpfjet_cluster_dR_.push_back(cluster_dR);
454  }
455  int cluster_ind = tpfjet_clusters[cluster_dR];
456 
457  std::vector<std::pair<DetId,float>> hitsAndFracs = cluster.hitsAndFractions();
458 
459  // Run over hits and match
460  int nHits = hitsAndFracs.size();
461  for(int iHit=0; iHit<nHits; iHit++){
462  int etaPhiPF = getEtaPhi(hitsAndFracs[iHit].first);
463 
464  for(edm::SortedCollection<HBHERecHit,edm::StrictWeakOrdering<HBHERecHit>>::const_iterator ith=hbhereco->begin(); ith!=hbhereco->end(); ++ith){
465  int etaPhiRecHit = getEtaPhi((*ith).id());
466  if(etaPhiPF == etaPhiRecHit){
468  if(tpfjet_rechits.count((*ith).id()) == 0){
469  tpfjet_twr_ieta_.push_back((*ith).id().ieta());
470  tpfjet_twr_iphi_.push_back((*ith).id().iphi());
471  tpfjet_twr_depth_.push_back((*ith).id().depth());
472  tpfjet_twr_subdet_.push_back((*ith).id().subdet());
473  if(hitsAndFracs[iHit].second > 0.05 && (*ith).energy() > 0.0) twrietas[(*ith).id().ieta()]++;
474  tpfjet_twr_hade_.push_back((*ith).energy());
475  tpfjet_twr_frac_.push_back(hitsAndFracs[iHit].second);
476  tpfjet_rechits[(*ith).id()].second.insert(hitsAndFracs[iHit].second);
477  tpfjet_twr_hadind_.push_back(tpfjet_had_n_ - 1);
478  tpfjet_twr_elmttype_.push_back(0);
479  tpfjet_twr_clusterind_.push_back(cluster_ind);
480  if(hasTrack){
482  }
483  else{
484  tpfjet_twr_candtrackind_.push_back(-1);
485  }
486  switch((*ith).id().subdet()){
488  {
489  const CaloCellGeometry *thisCell = HBGeom->getGeometry((*ith).id().rawId());
490  const CaloCellGeometry::CornersVec& cv = thisCell->getCorners();
491  float avgeta = (cv[0].eta() + cv[2].eta())/2.0;
492  float avgphi = (static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi()))/2.0;
493  if(cv[0].phi() < cv[2].phi()) avgphi = (2.0*3.141592653 + static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi()))/2.0;
494  tpfjet_twr_dR_.push_back(deltaR(tpfjet_eta_,tpfjet_phi_,avgeta,avgphi));
495  break;
496  }
498  {
499  const CaloCellGeometry *thisCell = HEGeom->getGeometry((*ith).id().rawId());
500  const CaloCellGeometry::CornersVec& cv = thisCell->getCorners();
501  float avgeta = (cv[0].eta() + cv[2].eta())/2.0;
502  float avgphi = (static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi()))/2.0;
503  if(cv[0].phi() < cv[2].phi()) avgphi = (2.0*3.141592653 + static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi()))/2.0;
504  tpfjet_twr_dR_.push_back(deltaR(tpfjet_eta_,tpfjet_phi_,avgeta,avgphi));
505  break;
506  }
507  default:
508  tpfjet_twr_dR_.push_back(-1);
509  break;
510  }
511  tpfjet_rechits[(*ith).id()].first = tpfjet_ntwrs_;
512  ++tpfjet_ntwrs_;
513  }
514  else if(tpfjet_rechits[(*ith).id()].second.count(hitsAndFracs[iHit].second) == 0){
515  tpfjet_twr_frac_.at(tpfjet_rechits[(*ith).id()].first) += hitsAndFracs[iHit].second;
516  if(cluster_dR < tpfjet_cluster_dR_.at(tpfjet_twr_clusterind_.at(tpfjet_rechits[(*ith).id()].first))){
517  tpfjet_twr_clusterind_.at(tpfjet_rechits[(*ith).id()].first) = cluster_ind;
518  }
519  tpfjet_rechits[(*ith).id()].second.insert(hitsAndFracs[iHit].second);
520  }
521  } // Test if ieta,iphi matches
522  } // Loop over rechits
523  } // Loop over hits
524  } // Test if element is from HCAL
525  else if(elements[iEle].type() == reco::PFBlockElement::HFHAD){ // Element is HF
526  for(edm::SortedCollection<HFRecHit,edm::StrictWeakOrdering<HFRecHit>>::const_iterator ith=hfreco->begin(); ith!=hfreco->end(); ++ith){
527  if((*ith).id().depth() == 1) continue; // Remove long fibers
528  const CaloCellGeometry *thisCell = HFGeom->getGeometry((*ith).id().rawId());
529  const CaloCellGeometry::CornersVec& cv = thisCell->getCorners();
530 
531  bool passMatch = false;
532  if((*it)->eta() < cv[0].eta() && (*it)->eta() > cv[2].eta()){
533  if((*it)->phi() < cv[0].phi() && (*it)->phi() > cv[2].phi()) passMatch = true;
534  else if(cv[0].phi() < cv[2].phi()){
535  if((*it)->phi() < cv[0].phi()) passMatch = true;
536  else if((*it)->phi() > cv[2].phi()) passMatch = true;
537  }
538  }
539 
540  if(passMatch){
542  tpfjet_twr_ieta_.push_back((*ith).id().ieta());
543  tpfjet_twr_iphi_.push_back((*ith).id().iphi());
544  tpfjet_twr_depth_.push_back((*ith).id().depth());
545  tpfjet_twr_subdet_.push_back((*ith).id().subdet());
546  tpfjet_twr_hade_.push_back((*ith).energy());
547  tpfjet_twr_frac_.push_back(1.0);
548  tpfjet_twr_hadind_.push_back(tpfjet_had_n_ - 1);
549  tpfjet_twr_elmttype_.push_back(1);
550  tpfjet_twr_clusterind_.push_back(-1);
551  tpfjet_twr_candtrackind_.push_back(-1);
552  float avgeta = (cv[0].eta() + cv[2].eta())/2.0;
553  float avgphi = (static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi()))/2.0;
554  if(cv[0].phi() < cv[2].phi()) avgphi = (2.0*3.141592653 + static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi()))/2.0;
555  tpfjet_twr_dR_.push_back(deltaR(tpfjet_eta_,tpfjet_phi_,avgeta,avgphi));
556  ++tpfjet_ntwrs_;
557  HFHAD_E += (*ith).energy();
558  }
559  }
560  }
561  else if(elements[iEle].type() == reco::PFBlockElement::HFEM){ // Element is HF
562  for(edm::SortedCollection<HFRecHit,edm::StrictWeakOrdering<HFRecHit>>::const_iterator ith=hfreco->begin(); ith!=hfreco->end(); ++ith){
563  if((*ith).id().depth() == 2) continue; // Remove short fibers
564  const CaloCellGeometry *thisCell = HFGeom->getGeometry((*ith).id().rawId());
565  const CaloCellGeometry::CornersVec& cv = thisCell->getCorners();
566 
567  bool passMatch = false;
568  if((*it)->eta() < cv[0].eta() && (*it)->eta() > cv[2].eta()){
569  if((*it)->phi() < cv[0].phi() && (*it)->phi() > cv[2].phi()) passMatch = true;
570  else if(cv[0].phi() < cv[2].phi()){
571  if((*it)->phi() < cv[0].phi()) passMatch = true;
572  else if((*it)->phi() > cv[2].phi()) passMatch = true;
573  }
574  }
575 
576  if(passMatch){
578  tpfjet_twr_ieta_.push_back((*ith).id().ieta());
579  tpfjet_twr_iphi_.push_back((*ith).id().iphi());
580  tpfjet_twr_depth_.push_back((*ith).id().depth());
581  tpfjet_twr_subdet_.push_back((*ith).id().subdet());
582  tpfjet_twr_hade_.push_back((*ith).energy());
583  tpfjet_twr_frac_.push_back(1.0);
584  tpfjet_twr_hadind_.push_back(tpfjet_had_n_ - 1);
585  tpfjet_twr_elmttype_.push_back(2);
586  tpfjet_twr_clusterind_.push_back(-1);
587  tpfjet_twr_candtrackind_.push_back(-1);
588  float avgeta = (cv[0].eta() + cv[2].eta())/2.0;
589  float avgphi = (static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi()))/2.0;
590  if(cv[0].phi() < cv[2].phi()) avgphi = (2.0*3.141592653 + static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi()))/2.0;
591  tpfjet_twr_dR_.push_back(deltaR(tpfjet_eta_,tpfjet_phi_,avgeta,avgphi));
592  ++tpfjet_ntwrs_;
593  HFEM_E += (*ith).energy();
594  }
595  }
596  }
597  else if(elements[iEle].type() == reco::PFBlockElement::HO){ // Element is HO
598  reco::PFClusterRef clusterref = elements[iEle].clusterRef();
599  reco::PFCluster cluster = *clusterref;
600  double cluster_dR = deltaR(tpfjet_eta_,tpfjet_phi_,cluster.eta(),cluster.phi());
601  if(tpfjet_clusters.count(cluster_dR) == 0){
602  tpfjet_clusters[cluster_dR] = tpfjet_cluster_n_;
603  tpfjet_cluster_eta_.push_back(cluster.eta());
604  tpfjet_cluster_phi_.push_back(cluster.phi());
605  tpfjet_cluster_dR_.push_back(cluster_dR);
607  }
608  int cluster_ind = tpfjet_clusters[cluster_dR];
609 
610  std::vector<std::pair<DetId,float>> hitsAndFracs = cluster.hitsAndFractions();
611  int nHits = hitsAndFracs.size();
612  for(int iHit=0; iHit<nHits; iHit++){
613  int etaPhiPF = getEtaPhi(hitsAndFracs[iHit].first);
614 
615  for(edm::SortedCollection<HORecHit,edm::StrictWeakOrdering<HORecHit>>::const_iterator ith=horeco->begin(); ith!=horeco->end(); ++ith){
616  int etaPhiRecHit = getEtaPhi((*ith).id());
617  if(etaPhiPF == etaPhiRecHit){
619  if(tpfjet_rechits.count((*ith).id()) == 0){
620  tpfjet_twr_ieta_.push_back((*ith).id().ieta());
621  tpfjet_twr_iphi_.push_back((*ith).id().iphi());
622  tpfjet_twr_depth_.push_back((*ith).id().depth());
623  tpfjet_twr_subdet_.push_back((*ith).id().subdet());
624  if(hitsAndFracs[iHit].second > 0.05 && (*ith).energy() > 0.0) twrietas[(*ith).id().ieta()]++;
625  tpfjet_twr_hade_.push_back((*ith).energy());
626  tpfjet_twr_frac_.push_back(hitsAndFracs[iHit].second);
627  tpfjet_rechits[(*ith).id()].second.insert(hitsAndFracs[iHit].second);
628  tpfjet_twr_hadind_.push_back(tpfjet_had_n_ - 1);
629  tpfjet_twr_elmttype_.push_back(3);
630  tpfjet_twr_clusterind_.push_back(cluster_ind);
631  if(hasTrack){
633  }
634  else{
635  tpfjet_twr_candtrackind_.push_back(-1);
636  }
637  const CaloCellGeometry *thisCell = HOGeom->getGeometry((*ith).id().rawId());
638  const CaloCellGeometry::CornersVec& cv = thisCell->getCorners();
639  float avgeta = (cv[0].eta() + cv[2].eta())/2.0;
640  float avgphi = (static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi()))/2.0;
641  if(cv[0].phi() < cv[2].phi()) avgphi = (2.0*3.141592653 + static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi()))/2.0;
642  tpfjet_twr_dR_.push_back(deltaR(tpfjet_eta_,tpfjet_phi_,avgeta,avgphi));
643  tpfjet_rechits[(*ith).id()].first = tpfjet_ntwrs_;
644  ++tpfjet_ntwrs_;
645  }
646  else if(tpfjet_rechits[(*ith).id()].second.count(hitsAndFracs[iHit].second) == 0){
647  tpfjet_twr_frac_.at(tpfjet_rechits[(*ith).id()].first) += hitsAndFracs[iHit].second;
648  if(cluster_dR < tpfjet_cluster_dR_.at(tpfjet_twr_clusterind_.at(tpfjet_rechits[(*ith).id()].first))){
649  tpfjet_twr_clusterind_.at(tpfjet_rechits[(*ith).id()].first) = cluster_ind;
650  }
651  tpfjet_rechits[(*ith).id()].second.insert(hitsAndFracs[iHit].second);
652  }
653  } // Test if ieta,iphi match
654  } // Loop over rechits
655  } // Loop over hits
656  } // Test if element is HO
657  } // Test for right element index
658  } // Loop over elements
659  } // Loop over elements in blocks
660 
661  switch(candidateType){
663  tpfjet_had_emf_.push_back(HFEM_E/(HFEM_E + HFHAD_E));
664  break;
666  tpfjet_had_emf_.push_back(-1);
667  break;
668  default:
669  tpfjet_had_emf_.push_back(-1);
670  break;
671  }
672  } // Loop over PF constitutents
673 
674  int tag_had_EcalE = 0;
675  int tag_had_rawHcalE = 0;
676  for(int i=0; i<tpfjet_had_n_; i++){
677  tag_had_EcalE += tpfjet_had_EcalE_[i];
678  tag_had_rawHcalE += tpfjet_had_rawHcalE_[i];
679  }
680  tpfjet_EMfrac_ = 1.0 - tag_had_rawHcalE/(tag_had_rawHcalE + tag_had_EcalE + tpfjet_unkown_E_ + tpfjet_electron_E_ + tpfjet_muon_E_ + tpfjet_photon_E_);
681  tpfjet_hadEcalEfrac_ = tag_had_EcalE/(tag_had_rawHcalE + tag_had_EcalE + tpfjet_unkown_E_ + tpfjet_electron_E_ + tpfjet_muon_E_ + tpfjet_photon_E_);
682 
683  // fill probe jet variables
684  ppfjet_pt_ = pf_probe.jet()->pt();
685  ppfjet_p_ = pf_probe.jet()->p();
686  ppfjet_E_ = pf_probe.jet()->energy();
687  ppfjet_eta_ = pf_probe.jet()->eta();
688  ppfjet_phi_ = pf_probe.jet()->phi();
689  ppfjet_scale_ = pf_probe.scale();
690  ppfjet_area_ = pf_probe.jet()->jetArea();
691  ppfjet_ntwrs_=0;
693 
694  ppfjet_jetID_ = 0; // Not a loose, medium, or tight jet
695  if(fabs(pf_probe.jet()->eta()) < 2.4){
696  if(pf_probe.jet()->chargedHadronEnergyFraction() > 0 &&
697  pf_probe.jet()->chargedMultiplicity() > 0 &&
698  pf_probe.jet()->chargedEmEnergyFraction() < 0.99 &&
699  (pf_probe.jet()->chargedMultiplicity() + pf_probe.jet()->neutralMultiplicity()) > 1){
700  if(pf_probe.jet()->neutralHadronEnergyFraction() < 0.9 &&
701  pf_probe.jet()->neutralEmEnergyFraction() < 0.9){
702  ppfjet_jetID_ = 3; // Tight jet
703  }
704  else if(pf_probe.jet()->neutralHadronEnergyFraction() < 0.95 &&
705  pf_probe.jet()->neutralEmEnergyFraction() < 0.95){
706  ppfjet_jetID_ = 2; // Medium jet
707  }
708  else if(pf_probe.jet()->neutralHadronEnergyFraction() < 0.99 &&
709  pf_probe.jet()->neutralEmEnergyFraction() < 0.99){
710  ppfjet_jetID_ = 1; // Loose jet
711  }
712  }
713  }
714  else if((pf_probe.jet()->chargedMultiplicity() + pf_probe.jet()->neutralMultiplicity()) > 1){
715  if(pf_probe.jet()->neutralHadronEnergyFraction() < 0.9 &&
716  pf_probe.jet()->neutralEmEnergyFraction() < 0.9){
717  ppfjet_jetID_ = 3; // Tight jet
718  }
719  else if(pf_probe.jet()->neutralHadronEnergyFraction() < 0.95 &&
720  pf_probe.jet()->neutralEmEnergyFraction() < 0.95){
721  ppfjet_jetID_ = 2; // Medium jet
722  }
723  else if(pf_probe.jet()->neutralHadronEnergyFraction() < 0.99 &&
724  pf_probe.jet()->neutralEmEnergyFraction() < 0.99){
725  ppfjet_jetID_ = 1; // Loose jet
726  }
727  }
728 
729  // Get PF constituents and fill HCAL towers
730  std::vector<reco::PFCandidatePtr> probeconst=pf_probe.jet()->getPFConstituents();
731  for(std::vector<reco::PFCandidatePtr>::const_iterator it=probeconst.begin(); it!=probeconst.end(); ++it){
732  bool hasTrack = false;
733  reco::PFCandidate::ParticleType candidateType = (*it)->particleId();
734  switch(candidateType){
736  ppfjet_unkown_E_ += (*it)->energy();
737  ppfjet_unkown_px_ += (*it)->px();
738  ppfjet_unkown_py_ += (*it)->py();
739  ppfjet_unkown_pz_ += (*it)->pz();
740  ppfjet_unkown_EcalE_ += (*it)->ecalEnergy();
742  continue;
744  {
745  ppfjet_had_E_.push_back((*it)->energy());
746  ppfjet_had_px_.push_back((*it)->px());
747  ppfjet_had_py_.push_back((*it)->py());
748  ppfjet_had_pz_.push_back((*it)->pz());
749  ppfjet_had_EcalE_.push_back((*it)->ecalEnergy());
750  ppfjet_had_rawHcalE_.push_back((*it)->rawHcalEnergy());
751  ppfjet_had_id_.push_back(0);
752  ppfjet_had_ntwrs_.push_back(0);
753  ppfjet_had_n_++;
754 
755  reco::TrackRef trackRef = (*it)->trackRef();
756  if(trackRef.isNonnull()){
757  reco::Track track = *trackRef;
758  ppfjet_candtrack_px_.push_back(track.px());
759  ppfjet_candtrack_py_.push_back(track.py());
760  ppfjet_candtrack_pz_.push_back(track.pz());
761  ppfjet_candtrack_EcalE_.push_back((*it)->ecalEnergy());
763  hasTrack = true;
765  }
766  else{
767  ppfjet_had_candtrackind_.push_back(-2);
768  }
769  }
770  break;
772  ppfjet_electron_E_ += (*it)->energy();
773  ppfjet_electron_px_ += (*it)->px();
774  ppfjet_electron_py_ += (*it)->py();
775  ppfjet_electron_pz_ += (*it)->pz();
776  ppfjet_electron_EcalE_ += (*it)->ecalEnergy();
778  continue;
780  ppfjet_muon_E_ += (*it)->energy();
781  ppfjet_muon_px_ += (*it)->px();
782  ppfjet_muon_py_ += (*it)->py();
783  ppfjet_muon_pz_ += (*it)->pz();
784  ppfjet_muon_EcalE_ += (*it)->ecalEnergy();
785  ppfjet_muon_n_++;
786  continue;
788  ppfjet_photon_E_ += (*it)->energy();
789  ppfjet_photon_px_ += (*it)->px();
790  ppfjet_photon_py_ += (*it)->py();
791  ppfjet_photon_pz_ += (*it)->pz();
792  ppfjet_photon_EcalE_ += (*it)->ecalEnergy();
794  continue;
796  {
797  ppfjet_had_E_.push_back((*it)->energy());
798  ppfjet_had_px_.push_back((*it)->px());
799  ppfjet_had_py_.push_back((*it)->py());
800  ppfjet_had_pz_.push_back((*it)->pz());
801  ppfjet_had_EcalE_.push_back((*it)->ecalEnergy());
802  ppfjet_had_rawHcalE_.push_back((*it)->rawHcalEnergy());
803  ppfjet_had_id_.push_back(1);
804  ppfjet_had_candtrackind_.push_back(-1);
805  ppfjet_had_ntwrs_.push_back(0);
806  ppfjet_had_n_++;
807  break;
808  }
810  {
811  ppfjet_had_E_.push_back((*it)->energy());
812  ppfjet_had_px_.push_back((*it)->px());
813  ppfjet_had_py_.push_back((*it)->py());
814  ppfjet_had_pz_.push_back((*it)->pz());
815  ppfjet_had_EcalE_.push_back((*it)->ecalEnergy());
816  ppfjet_had_rawHcalE_.push_back((*it)->rawHcalEnergy());
817  ppfjet_had_id_.push_back(2);
818  ppfjet_had_candtrackind_.push_back(-1);
819  ppfjet_had_ntwrs_.push_back(0);
820  ppfjet_had_n_++;
821  break;
822  }
824  {
825  ppfjet_had_E_.push_back((*it)->energy());
826  ppfjet_had_px_.push_back((*it)->px());
827  ppfjet_had_py_.push_back((*it)->py());
828  ppfjet_had_pz_.push_back((*it)->pz());
829  ppfjet_had_EcalE_.push_back((*it)->ecalEnergy());
830  ppfjet_had_rawHcalE_.push_back((*it)->rawHcalEnergy());
831  ppfjet_had_id_.push_back(3);
832  ppfjet_had_candtrackind_.push_back(-1);
833  ppfjet_had_ntwrs_.push_back(0);
834  ppfjet_had_n_++;
835  break;
836  }
837  }
838 
839  float HFHAD_E = 0;
840  float HFEM_E = 0;
841  int maxElement=(*it)->elementsInBlocks().size();
842  for(int e=0; e<maxElement; ++e){
843  // Get elements from block
844  reco::PFBlockRef blockRef = (*it)->elementsInBlocks()[e].first;
845  const edm::OwnVector<reco::PFBlockElement>& elements = blockRef->elements();
846  for(unsigned iEle=0; iEle<elements.size(); iEle++) {
847  if(elements[iEle].index() == (*it)->elementsInBlocks()[e].second){
848  if(elements[iEle].type() == reco::PFBlockElement::HCAL){ // Element is HB or HE
849  // Get cluster and hits
850  reco::PFClusterRef clusterref = elements[iEle].clusterRef();
851  reco::PFCluster cluster = *clusterref;
852  double cluster_dR = deltaR(ppfjet_eta_,ppfjet_phi_,cluster.eta(),cluster.phi());
853  if(ppfjet_clusters.count(cluster_dR) == 0){
854  ppfjet_clusters[cluster_dR] = ppfjet_cluster_n_;
855  ppfjet_cluster_eta_.push_back(cluster.eta());
856  ppfjet_cluster_phi_.push_back(cluster.phi());
857  ppfjet_cluster_dR_.push_back(cluster_dR);
859  }
860  int cluster_ind = ppfjet_clusters[cluster_dR];
861  std::vector<std::pair<DetId,float>> hitsAndFracs = cluster.hitsAndFractions();
862 
863  // Run over hits and match
864  int nHits = hitsAndFracs.size();
865  for(int iHit=0; iHit<nHits; iHit++){
866  int etaPhiPF = getEtaPhi(hitsAndFracs[iHit].first);
867 
868  for(edm::SortedCollection<HBHERecHit,edm::StrictWeakOrdering<HBHERecHit>>::const_iterator ith=hbhereco->begin(); ith!=hbhereco->end(); ++ith){
869  int etaPhiRecHit = getEtaPhi((*ith).id());
870  if(etaPhiPF == etaPhiRecHit){
872  if(ppfjet_rechits.count((*ith).id()) == 0){
873  ppfjet_twr_ieta_.push_back((*ith).id().ieta());
874  ppfjet_twr_iphi_.push_back((*ith).id().iphi());
875  ppfjet_twr_depth_.push_back((*ith).id().depth());
876  ppfjet_twr_subdet_.push_back((*ith).id().subdet());
877  ppfjet_twr_hade_.push_back((*ith).energy());
878  ppfjet_twr_frac_.push_back(hitsAndFracs[iHit].second);
879  ppfjet_rechits[(*ith).id()].second.insert(hitsAndFracs[iHit].second);
880  ppfjet_twr_hadind_.push_back(ppfjet_had_n_ - 1);
881  ppfjet_twr_elmttype_.push_back(0);
882  ppfjet_twr_clusterind_.push_back(cluster_ind);
883  if(hasTrack){
885  }
886  else{
887  ppfjet_twr_candtrackind_.push_back(-1);
888  }
889  switch((*ith).id().subdet()){
891  {
892  const CaloCellGeometry *thisCell = HBGeom->getGeometry((*ith).id().rawId());
893  const CaloCellGeometry::CornersVec& cv = thisCell->getCorners();
894  float avgeta = (cv[0].eta() + cv[2].eta())/2.0;
895  float avgphi = (static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi()))/2.0;
896  if(cv[0].phi() < cv[2].phi()) avgphi = (2.0*3.141592653 + static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi()))/2.0;
897  ppfjet_twr_dR_.push_back(deltaR(ppfjet_eta_,ppfjet_phi_,avgeta,avgphi));
898  break;
899  }
901  {
902  const CaloCellGeometry *thisCell = HEGeom->getGeometry((*ith).id().rawId());
903  const CaloCellGeometry::CornersVec& cv = thisCell->getCorners();
904  float avgeta = (cv[0].eta() + cv[2].eta())/2.0;
905  float avgphi = (static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi()))/2.0;
906  if(cv[0].phi() < cv[2].phi()) avgphi = (2.0*3.141592653 + static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi()))/2.0;
907  ppfjet_twr_dR_.push_back(deltaR(ppfjet_eta_,ppfjet_phi_,avgeta,avgphi));
908  break;
909  }
910  default:
911  ppfjet_twr_dR_.push_back(-1);
912  break;
913  }
914  ppfjet_rechits[(*ith).id()].first = ppfjet_ntwrs_;
915  ++ppfjet_ntwrs_;
916  }
917  else if(ppfjet_rechits[(*ith).id()].second.count(hitsAndFracs[iHit].second) == 0){
918  ppfjet_twr_frac_.at(ppfjet_rechits[(*ith).id()].first) += hitsAndFracs[iHit].second;
919  if(cluster_dR < ppfjet_cluster_dR_.at(ppfjet_twr_clusterind_.at(ppfjet_rechits[(*ith).id()].first))){
920  ppfjet_twr_clusterind_.at(ppfjet_rechits[(*ith).id()].first) = cluster_ind;
921  }
922  ppfjet_rechits[(*ith).id()].second.insert(hitsAndFracs[iHit].second);
923  }
924  } // Test if ieta,iphi matches
925  } // Loop over rechits
926  } // Loop over hits
927  } // Test if element is from HCAL
928  else if(elements[iEle].type() == reco::PFBlockElement::HFHAD){ // Element is HF
929  for(edm::SortedCollection<HFRecHit,edm::StrictWeakOrdering<HFRecHit>>::const_iterator ith=hfreco->begin(); ith!=hfreco->end(); ++ith){
930  if((*ith).id().depth() == 1) continue; // Remove long fibers
931  const CaloCellGeometry *thisCell = HFGeom->getGeometry((*ith).id().rawId());
932  const CaloCellGeometry::CornersVec& cv = thisCell->getCorners();
933 
934  bool passMatch = false;
935  if((*it)->eta() < cv[0].eta() && (*it)->eta() > cv[2].eta()){
936  if((*it)->phi() < cv[0].phi() && (*it)->phi() > cv[2].phi()) passMatch = true;
937  else if(cv[0].phi() < cv[2].phi()){
938  if((*it)->phi() < cv[0].phi()) passMatch = true;
939  else if((*it)->phi() > cv[2].phi()) passMatch = true;
940  }
941  }
942 
943  if(passMatch){
945  ppfjet_twr_ieta_.push_back((*ith).id().ieta());
946  ppfjet_twr_iphi_.push_back((*ith).id().iphi());
947  ppfjet_twr_depth_.push_back((*ith).id().depth());
948  ppfjet_twr_subdet_.push_back((*ith).id().subdet());
949  ppfjet_twr_hade_.push_back((*ith).energy());
950  ppfjet_twr_frac_.push_back(1.0);
951  ppfjet_twr_hadind_.push_back(ppfjet_had_n_ - 1);
952  ppfjet_twr_elmttype_.push_back(1);
953  ppfjet_twr_clusterind_.push_back(-1);
954  ppfjet_twr_candtrackind_.push_back(-1);
955  float avgeta = (cv[0].eta() + cv[2].eta())/2.0;
956  float avgphi = (static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi()))/2.0;
957  if(cv[0].phi() < cv[2].phi()) avgphi = (2.0*3.141592653 + static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi()))/2.0;
958  ppfjet_twr_dR_.push_back(deltaR(ppfjet_eta_,ppfjet_phi_,avgeta,avgphi));
959  ++ppfjet_ntwrs_;
960  HFHAD_E += (*ith).energy();
961  }
962  }
963  }
964  else if(elements[iEle].type() == reco::PFBlockElement::HFEM){ // Element is HF
965  for(edm::SortedCollection<HFRecHit,edm::StrictWeakOrdering<HFRecHit>>::const_iterator ith=hfreco->begin(); ith!=hfreco->end(); ++ith){
966  if((*ith).id().depth() == 2) continue; // Remove short fibers
967  const CaloCellGeometry *thisCell = HFGeom->getGeometry((*ith).id().rawId());
968  const CaloCellGeometry::CornersVec& cv = thisCell->getCorners();
969 
970  bool passMatch = false;
971  if((*it)->eta() < cv[0].eta() && (*it)->eta() > cv[2].eta()){
972  if((*it)->phi() < cv[0].phi() && (*it)->phi() > cv[2].phi()) passMatch = true;
973  else if(cv[0].phi() < cv[2].phi()){
974  if((*it)->phi() < cv[0].phi()) passMatch = true;
975  else if((*it)->phi() > cv[2].phi()) passMatch = true;
976  }
977  }
978 
979  if(passMatch){
981  ppfjet_twr_ieta_.push_back((*ith).id().ieta());
982  ppfjet_twr_iphi_.push_back((*ith).id().iphi());
983  ppfjet_twr_depth_.push_back((*ith).id().depth());
984  ppfjet_twr_subdet_.push_back((*ith).id().subdet());
985  ppfjet_twr_hade_.push_back((*ith).energy());
986  ppfjet_twr_frac_.push_back(1.0);
987  ppfjet_twr_hadind_.push_back(ppfjet_had_n_ - 1);
988  ppfjet_twr_elmttype_.push_back(2);
989  ppfjet_twr_clusterind_.push_back(-1);
990  ppfjet_twr_candtrackind_.push_back(-1);
991  float avgeta = (cv[0].eta() + cv[2].eta())/2.0;
992  float avgphi = (static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi()))/2.0;
993  if(cv[0].phi() < cv[2].phi()) avgphi = (2.0*3.141592653 + static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi()))/2.0;
994  ppfjet_twr_dR_.push_back(deltaR(ppfjet_eta_,ppfjet_phi_,avgeta,avgphi));
995  ++ppfjet_ntwrs_;
996  HFEM_E += (*ith).energy();
997  }
998  }
999  }
1000  else if(elements[iEle].type() == reco::PFBlockElement::HO){ // Element is HO
1001  reco::PFClusterRef clusterref = elements[iEle].clusterRef();
1002  reco::PFCluster cluster = *clusterref;
1003  double cluster_dR = deltaR(ppfjet_eta_,ppfjet_phi_,cluster.eta(),cluster.phi());
1004  if(ppfjet_clusters.count(cluster_dR) == 0){
1005  ppfjet_clusters[cluster_dR] = ppfjet_cluster_n_;
1006  ppfjet_cluster_eta_.push_back(cluster.eta());
1007  ppfjet_cluster_phi_.push_back(cluster.phi());
1008  ppfjet_cluster_dR_.push_back(cluster_dR);
1010  }
1011  int cluster_ind = ppfjet_clusters[cluster_dR];
1012 
1013  std::vector<std::pair<DetId,float>> hitsAndFracs = cluster.hitsAndFractions();
1014  int nHits = hitsAndFracs.size();
1015  for(int iHit=0; iHit<nHits; iHit++){
1016  int etaPhiPF = getEtaPhi(hitsAndFracs[iHit].first);
1017 
1018  for(edm::SortedCollection<HORecHit,edm::StrictWeakOrdering<HORecHit>>::const_iterator ith=horeco->begin(); ith!=horeco->end(); ++ith){
1019  int etaPhiRecHit = getEtaPhi((*ith).id());
1020  if(etaPhiPF == etaPhiRecHit){
1021  ppfjet_had_ntwrs_.at(ppfjet_had_n_ - 1)++;
1022  if(ppfjet_rechits.count((*ith).id()) == 0){
1023  ppfjet_twr_ieta_.push_back((*ith).id().ieta());
1024  ppfjet_twr_iphi_.push_back((*ith).id().iphi());
1025  ppfjet_twr_depth_.push_back((*ith).id().depth());
1026  ppfjet_twr_subdet_.push_back((*ith).id().subdet());
1027  ppfjet_twr_hade_.push_back((*ith).energy());
1028  ppfjet_twr_frac_.push_back(hitsAndFracs[iHit].second);
1029  ppfjet_rechits[(*ith).id()].second.insert(hitsAndFracs[iHit].second);
1030  ppfjet_twr_hadind_.push_back(ppfjet_had_n_ - 1);
1031  ppfjet_twr_elmttype_.push_back(3);
1032  ppfjet_twr_clusterind_.push_back(cluster_ind);
1033  if(hasTrack){
1035  }
1036  else{
1037  ppfjet_twr_candtrackind_.push_back(-1);
1038  }
1039  const CaloCellGeometry *thisCell = HOGeom->getGeometry((*ith).id().rawId());
1040  const CaloCellGeometry::CornersVec& cv = thisCell->getCorners();
1041  float avgeta = (cv[0].eta() + cv[2].eta())/2.0;
1042  float avgphi = (static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi()))/2.0;
1043  if(cv[0].phi() < cv[2].phi()) avgphi = (2.0*3.141592653 + static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi()))/2.0;
1044  ppfjet_twr_dR_.push_back(deltaR(ppfjet_eta_,ppfjet_phi_,avgeta,avgphi));
1045  ppfjet_rechits[(*ith).id()].first = ppfjet_ntwrs_;
1046  ++ppfjet_ntwrs_;
1047  }
1048  else if(ppfjet_rechits[(*ith).id()].second.count(hitsAndFracs[iHit].second) == 0){
1049  ppfjet_twr_frac_.at(ppfjet_rechits[(*ith).id()].first) += hitsAndFracs[iHit].second;
1050  if(cluster_dR < ppfjet_cluster_dR_.at(ppfjet_twr_clusterind_.at(ppfjet_rechits[(*ith).id()].first))){
1051  ppfjet_twr_clusterind_.at(ppfjet_rechits[(*ith).id()].first) = cluster_ind;
1052  }
1053  ppfjet_rechits[(*ith).id()].second.insert(hitsAndFracs[iHit].second);
1054  }
1055  } // Test if ieta,iphi match
1056  } // Loop over rechits
1057  } // Loop over hits
1058  } // Test if element is from HO
1059  } // Test for right element index
1060  } // Loop over elements
1061  } // Loop over elements in blocks
1062  switch(candidateType){
1064  ppfjet_had_emf_.push_back(HFEM_E/(HFEM_E + HFHAD_E));
1065  break;
1067  ppfjet_had_emf_.push_back(-1);
1068  break;
1069  default:
1070  ppfjet_had_emf_.push_back(-1);
1071  break;
1072  }
1073  } // Loop over PF constitutents
1074 
1075  int probe_had_EcalE = 0;
1076  int probe_had_rawHcalE = 0;
1077  for(int i=0; i<ppfjet_had_n_; i++){
1078  probe_had_EcalE += ppfjet_had_EcalE_[i];
1079  probe_had_rawHcalE += ppfjet_had_rawHcalE_[i];
1080  }
1081  ppfjet_EMfrac_ = 1.0 - probe_had_rawHcalE/(probe_had_rawHcalE + probe_had_EcalE + ppfjet_unkown_E_ + ppfjet_electron_E_ + ppfjet_muon_E_ + ppfjet_photon_E_);
1082  ppfjet_hadEcalEfrac_ = probe_had_EcalE/(probe_had_rawHcalE + probe_had_EcalE + ppfjet_unkown_E_ + ppfjet_electron_E_ + ppfjet_muon_E_ + ppfjet_photon_E_);
1083 
1084  // fill dijet variables
1085  pf_dijet_deta_=std::fabs(std::fabs(pf_tag.jet()->eta())-std::fabs(pf_probe.jet()->eta()));
1086  pf_dijet_dphi_=pf_tag.jet()->phi()-pf_probe.jet()->phi();
1087  if(pf_dijet_dphi_>3.1415) pf_dijet_dphi_ = 6.2832-pf_dijet_dphi_;
1089 
1090  tree_->Fill();
1091 
1092  return;
1093 }
1094 
1095 // ------------ method called once each job just before starting event loop ------------
1097 {
1098  // book histograms
1099  rootfile_ = new TFile(rootHistFilename_.c_str(), "RECREATE");
1100 
1101  h_PassSelPF_ = new TH1D("h_PassSelectionPF", "Selection Pass Failures PFJets",200,-0.5,199.5);
1102 
1103  tree_ = new TTree("dijettree", "tree for dijet balancing");
1104 
1105  tree_->Branch("tpfjet_pt",&tpfjet_pt_, "tpfjet_pt/F");
1106  tree_->Branch("tpfjet_p",&tpfjet_p_, "tpfjet_p/F");
1107  tree_->Branch("tpfjet_E",&tpfjet_E_, "tpfjet_E/F");
1108  tree_->Branch("tpfjet_eta",&tpfjet_eta_, "tpfjet_eta/F");
1109  tree_->Branch("tpfjet_phi",&tpfjet_phi_, "tpfjet_phi/F");
1110  tree_->Branch("tpfjet_EMfrac",&tpfjet_EMfrac_, "tpfjet_EMfrac/F");
1111  tree_->Branch("tpfjet_hadEcalEfrac",&tpfjet_hadEcalEfrac_, "tpfjet_hadEcalEfrac/F");
1112  tree_->Branch("tpfjet_scale",&tpfjet_scale_, "tpfjet_scale/F");
1113  tree_->Branch("tpfjet_area",&tpfjet_area_, "tpfjet_area/F");
1114  tree_->Branch("tpfjet_jetID",&tpfjet_jetID_, "tpfjet_jetID/I");
1115  tree_->Branch("tpfjet_unkown_E",&tpfjet_unkown_E_, "tpfjet_unkown_E/F");
1116  tree_->Branch("tpfjet_electron_E",&tpfjet_electron_E_, "tpfjet_electron_E/F");
1117  tree_->Branch("tpfjet_muon_E",&tpfjet_muon_E_, "tpfjet_muon_E/F");
1118  tree_->Branch("tpfjet_photon_E",&tpfjet_photon_E_, "tpfjet_photon_E/F");
1119  tree_->Branch("tpfjet_unkown_px",&tpfjet_unkown_px_, "tpfjet_unkown_px/F");
1120  tree_->Branch("tpfjet_electron_px",&tpfjet_electron_px_, "tpfjet_electron_px/F");
1121  tree_->Branch("tpfjet_muon_px",&tpfjet_muon_px_, "tpfjet_muon_px/F");
1122  tree_->Branch("tpfjet_photon_px",&tpfjet_photon_px_, "tpfjet_photon_px/F");
1123  tree_->Branch("tpfjet_unkown_py",&tpfjet_unkown_py_, "tpfjet_unkown_py/F");
1124  tree_->Branch("tpfjet_electron_py",&tpfjet_electron_py_, "tpfjet_electron_py/F");
1125  tree_->Branch("tpfjet_muon_py",&tpfjet_muon_py_, "tpfjet_muon_py/F");
1126  tree_->Branch("tpfjet_photon_py",&tpfjet_photon_py_, "tpfjet_photon_py/F");
1127  tree_->Branch("tpfjet_unkown_pz",&tpfjet_unkown_pz_, "tpfjet_unkown_pz/F");
1128  tree_->Branch("tpfjet_electron_pz",&tpfjet_electron_pz_, "tpfjet_electron_pz/F");
1129  tree_->Branch("tpfjet_muon_pz",&tpfjet_muon_pz_, "tpfjet_muon_pz/F");
1130  tree_->Branch("tpfjet_photon_pz",&tpfjet_photon_pz_, "tpfjet_photon_pz/F");
1131  tree_->Branch("tpfjet_unkown_EcalE",&tpfjet_unkown_EcalE_, "tpfjet_unkown_EcalE/F");
1132  tree_->Branch("tpfjet_electron_EcalE",&tpfjet_electron_EcalE_, "tpfjet_electron_EcalE/F");
1133  tree_->Branch("tpfjet_muon_EcalE",&tpfjet_muon_EcalE_, "tpfjet_muon_EcalE/F");
1134  tree_->Branch("tpfjet_photon_EcalE",&tpfjet_photon_EcalE_, "tpfjet_photon_EcalE/F");
1135  tree_->Branch("tpfjet_unkown_n",&tpfjet_unkown_n_, "tpfjet_unkown_n/I");
1136  tree_->Branch("tpfjet_electron_n",&tpfjet_electron_n_, "tpfjet_electron_n/I");
1137  tree_->Branch("tpfjet_muon_n",&tpfjet_muon_n_, "tpfjet_muon_n/I");
1138  tree_->Branch("tpfjet_photon_n",&tpfjet_photon_n_, "tpfjet_photon_n/I");
1139  tree_->Branch("tpfjet_had_n",&tpfjet_had_n_, "tpfjet_had_n/I");
1140  tree_->Branch("tpfjet_had_E",&tpfjet_had_E_);
1141  tree_->Branch("tpfjet_had_px",&tpfjet_had_px_);
1142  tree_->Branch("tpfjet_had_py",&tpfjet_had_py_);
1143  tree_->Branch("tpfjet_had_pz",&tpfjet_had_pz_);
1144  tree_->Branch("tpfjet_had_EcalE",&tpfjet_had_EcalE_);
1145  tree_->Branch("tpfjet_had_rawHcalE",&tpfjet_had_rawHcalE_);
1146  tree_->Branch("tpfjet_had_emf",&tpfjet_had_emf_);
1147  tree_->Branch("tpfjet_had_id",&tpfjet_had_id_);
1148  tree_->Branch("tpfjet_had_candtrackind",&tpfjet_had_candtrackind_);
1149  tree_->Branch("tpfjet_had_ntwrs",&tpfjet_had_ntwrs_);
1150  tree_->Branch("tpfjet_ntwrs",&tpfjet_ntwrs_, "tpfjet_ntwrs/I");
1151  tree_->Branch("tpfjet_twr_ieta",&tpfjet_twr_ieta_);
1152  tree_->Branch("tpfjet_twr_iphi",&tpfjet_twr_iphi_);
1153  tree_->Branch("tpfjet_twr_depth",&tpfjet_twr_depth_);
1154  tree_->Branch("tpfjet_twr_subdet",&tpfjet_twr_subdet_);
1155  tree_->Branch("tpfjet_twr_hade",&tpfjet_twr_hade_);
1156  tree_->Branch("tpfjet_twr_frac",&tpfjet_twr_frac_);
1157  tree_->Branch("tpfjet_twr_candtrackind",&tpfjet_twr_candtrackind_);
1158  tree_->Branch("tpfjet_twr_hadind",&tpfjet_twr_hadind_);
1159  tree_->Branch("tpfjet_twr_elmttype",&tpfjet_twr_elmttype_);
1160  tree_->Branch("tpfjet_twr_dR",&tpfjet_twr_dR_);
1161  tree_->Branch("tpfjet_twr_clusterind",&tpfjet_twr_clusterind_);
1162  tree_->Branch("tpfjet_cluster_n",&tpfjet_cluster_n_, "tpfjet_cluster_n/I");
1163  tree_->Branch("tpfjet_cluster_eta",&tpfjet_cluster_eta_);
1164  tree_->Branch("tpfjet_cluster_phi",&tpfjet_cluster_phi_);
1165  tree_->Branch("tpfjet_cluster_dR",&tpfjet_cluster_dR_);
1166  tree_->Branch("tpfjet_ncandtracks",&tpfjet_ncandtracks_, "tpfjet_ncandtracks/I");
1167  tree_->Branch("tpfjet_candtrack_px",&tpfjet_candtrack_px_);
1168  tree_->Branch("tpfjet_candtrack_py",&tpfjet_candtrack_py_);
1169  tree_->Branch("tpfjet_candtrack_pz",&tpfjet_candtrack_pz_);
1170  tree_->Branch("tpfjet_candtrack_EcalE",&tpfjet_candtrack_EcalE_);
1171  tree_->Branch("ppfjet_pt",&ppfjet_pt_, "ppfjet_pt/F");
1172  tree_->Branch("ppfjet_p",&ppfjet_p_, "ppfjet_p/F");
1173  tree_->Branch("ppfjet_E",&ppfjet_E_, "ppfjet_E/F");
1174  tree_->Branch("ppfjet_eta",&ppfjet_eta_, "ppfjet_eta/F");
1175  tree_->Branch("ppfjet_phi",&ppfjet_phi_, "ppfjet_phi/F");
1176  tree_->Branch("ppfjet_EMfrac",&ppfjet_EMfrac_, "ppfjet_EMfrac/F");
1177  tree_->Branch("ppfjet_hadEcalEfrac",&ppfjet_hadEcalEfrac_, "ppfjet_hadEcalEfrac/F");
1178  tree_->Branch("ppfjet_scale",&ppfjet_scale_, "ppfjet_scale/F");
1179  tree_->Branch("ppfjet_area",&ppfjet_area_, "ppfjet_area/F");
1180  tree_->Branch("ppfjet_jetID",&ppfjet_jetID_, "ppfjet_jetID/I");
1181  tree_->Branch("ppfjet_unkown_E",&ppfjet_unkown_E_, "ppfjet_unkown_E/F");
1182  tree_->Branch("ppfjet_electron_E",&ppfjet_electron_E_, "ppfjet_electron_E/F");
1183  tree_->Branch("ppfjet_muon_E",&ppfjet_muon_E_, "ppfjet_muon_E/F");
1184  tree_->Branch("ppfjet_photon_E",&ppfjet_photon_E_, "ppfjet_photon_E/F");
1185  tree_->Branch("ppfjet_unkown_px",&ppfjet_unkown_px_, "ppfjet_unkown_px/F");
1186  tree_->Branch("ppfjet_electron_px",&ppfjet_electron_px_, "ppfjet_electron_px/F");
1187  tree_->Branch("ppfjet_muon_px",&ppfjet_muon_px_, "ppfjet_muon_px/F");
1188  tree_->Branch("ppfjet_photon_px",&ppfjet_photon_px_, "ppfjet_photon_px/F");
1189  tree_->Branch("ppfjet_unkown_py",&ppfjet_unkown_py_, "ppfjet_unkown_py/F");
1190  tree_->Branch("ppfjet_electron_py",&ppfjet_electron_py_, "ppfjet_electron_py/F");
1191  tree_->Branch("ppfjet_muon_py",&ppfjet_muon_py_, "ppfjet_muon_py/F");
1192  tree_->Branch("ppfjet_photon_py",&ppfjet_photon_py_, "ppfjet_photon_py/F");
1193  tree_->Branch("ppfjet_unkown_pz",&ppfjet_unkown_pz_, "ppfjet_unkown_pz/F");
1194  tree_->Branch("ppfjet_electron_pz",&ppfjet_electron_pz_, "ppfjet_electron_pz/F");
1195  tree_->Branch("ppfjet_muon_pz",&ppfjet_muon_pz_, "ppfjet_muon_pz/F");
1196  tree_->Branch("ppfjet_photon_pz",&ppfjet_photon_pz_, "ppfjet_photon_pz/F");
1197  tree_->Branch("ppfjet_unkown_EcalE",&ppfjet_unkown_EcalE_, "ppfjet_unkown_EcalE/F");
1198  tree_->Branch("ppfjet_electron_EcalE",&ppfjet_electron_EcalE_, "ppfjet_electron_EcalE/F");
1199  tree_->Branch("ppfjet_muon_EcalE",&ppfjet_muon_EcalE_, "ppfjet_muon_EcalE/F");
1200  tree_->Branch("ppfjet_photon_EcalE",&ppfjet_photon_EcalE_, "ppfjet_photon_EcalE/F");
1201  tree_->Branch("ppfjet_unkown_n",&ppfjet_unkown_n_, "ppfjet_unkown_n/I");
1202  tree_->Branch("ppfjet_electron_n",&ppfjet_electron_n_, "ppfjet_electron_n/I");
1203  tree_->Branch("ppfjet_muon_n",&ppfjet_muon_n_, "ppfjet_muon_n/I");
1204  tree_->Branch("ppfjet_photon_n",&ppfjet_photon_n_, "ppfjet_photon_n/I");
1205  tree_->Branch("ppfjet_had_n",&ppfjet_had_n_, "ppfjet_had_n/I");
1206  tree_->Branch("ppfjet_had_E",&ppfjet_had_E_);
1207  tree_->Branch("ppfjet_had_px",&ppfjet_had_px_);
1208  tree_->Branch("ppfjet_had_py",&ppfjet_had_py_);
1209  tree_->Branch("ppfjet_had_pz",&ppfjet_had_pz_);
1210  tree_->Branch("ppfjet_had_EcalE",&ppfjet_had_EcalE_);
1211  tree_->Branch("ppfjet_had_rawHcalE",&ppfjet_had_rawHcalE_);
1212  tree_->Branch("ppfjet_had_emf",&ppfjet_had_emf_);
1213  tree_->Branch("ppfjet_had_id",&ppfjet_had_id_);
1214  tree_->Branch("ppfjet_had_candtrackind",&ppfjet_had_candtrackind_);
1215  tree_->Branch("ppfjet_had_ntwrs",&ppfjet_had_ntwrs_);
1216  tree_->Branch("ppfjet_ntwrs",&ppfjet_ntwrs_, "ppfjet_ntwrs/I");
1217  tree_->Branch("ppfjet_twr_ieta",&ppfjet_twr_ieta_);
1218  tree_->Branch("ppfjet_twr_iphi",&ppfjet_twr_iphi_);
1219  tree_->Branch("ppfjet_twr_depth",&ppfjet_twr_depth_);
1220  tree_->Branch("ppfjet_twr_subdet",&ppfjet_twr_subdet_);
1221  tree_->Branch("ppfjet_twr_hade",&ppfjet_twr_hade_);
1222  tree_->Branch("ppfjet_twr_frac",&ppfjet_twr_frac_);
1223  tree_->Branch("ppfjet_twr_candtrackind",&ppfjet_twr_candtrackind_);
1224  tree_->Branch("ppfjet_twr_hadind",&ppfjet_twr_hadind_);
1225  tree_->Branch("ppfjet_twr_elmttype",&ppfjet_twr_elmttype_);
1226  tree_->Branch("ppfjet_twr_dR",&ppfjet_twr_dR_);
1227  tree_->Branch("ppfjet_twr_clusterind",&ppfjet_twr_clusterind_);
1228  tree_->Branch("ppfjet_cluster_n",&ppfjet_cluster_n_, "ppfjet_cluster_n/I");
1229  tree_->Branch("ppfjet_cluster_eta",&ppfjet_cluster_eta_);
1230  tree_->Branch("ppfjet_cluster_phi",&ppfjet_cluster_phi_);
1231  tree_->Branch("ppfjet_cluster_dR",&ppfjet_cluster_dR_);
1232  tree_->Branch("ppfjet_ncandtracks",&ppfjet_ncandtracks_, "ppfjet_ncandtracks/I");
1233  tree_->Branch("ppfjet_candtrack_px",&ppfjet_candtrack_px_);
1234  tree_->Branch("ppfjet_candtrack_py",&ppfjet_candtrack_py_);
1235  tree_->Branch("ppfjet_candtrack_pz",&ppfjet_candtrack_pz_);
1236  tree_->Branch("ppfjet_candtrack_EcalE",&ppfjet_candtrack_EcalE_);
1237  tree_->Branch("pf_dijet_deta",&pf_dijet_deta_, "pf_dijet_deta/F");
1238  tree_->Branch("pf_dijet_dphi",&pf_dijet_dphi_, "pf_dijet_dphi/F");
1239  tree_->Branch("pf_dijet_balance",&pf_dijet_balance_, "pf_dijet_balance/F");
1240  tree_->Branch("pf_thirdjet_px",&pf_thirdjet_px_, "pf_thirdjet_px/F");
1241  tree_->Branch("pf_thirdjet_py",&pf_thirdjet_py_, "pf_thirdjet_py/F");
1242  tree_->Branch("pf_realthirdjet_px",&pf_realthirdjet_px_, "pf_realthirdjet_px/F");
1243  tree_->Branch("pf_realthirdjet_py",&pf_realthirdjet_py_, "pf_realthirdjet_py/F");
1244  tree_->Branch("pf_realthirdjet_scale",&pf_realthirdjet_scale_, "pf_realthirdjet_scale/F");
1245  tree_->Branch("pf_Run",&pf_Run_, "pf_Run/I");
1246  tree_->Branch("pf_Lumi",&pf_Lumi_, "pf_Lumi/I");
1247  tree_->Branch("pf_Event",&pf_Event_, "pf_Event/I");
1248  tree_->Branch("pf_NPV",&pf_NPV_, "pf_NPV/I");
1249 
1250  return;
1251 }
1252 
1253 // ------------ method called once each job just after ending the event loop ------------
1254 void
1256  // write histograms
1257  rootfile_->cd();
1258 
1259  h_PassSelPF_->Write();
1260  tree_->Write();
1261 
1262  rootfile_->Close();
1263 
1264 }
1265 
1266 // helper function
1267 
1268 double DiJetAnalyzer::deltaR(const reco::Jet* j1, const reco::Jet* j2)
1269 {
1270  double deta = j1->eta()-j2->eta();
1271  double dphi = std::fabs(j1->phi()-j2->phi());
1272  if(dphi>3.1415927) dphi = 2*3.1415927 - dphi;
1273  return std::sqrt(deta*deta + dphi*dphi);
1274 }
1275 
1276 double DiJetAnalyzer::deltaR(const double eta1, const double phi1, const double eta2, const double phi2)
1277 {
1278  double deta = eta1 - eta2;
1279  double dphi = std::fabs(phi1 - phi2);
1280  if(dphi>3.1415927) dphi = 2*3.1415927 - dphi;
1281  return std::sqrt(deta*deta + dphi*dphi);
1282 }
1283 
1285 {
1286  return id.rawId() & 0x3FFF; // Get 14 least-significant digits
1287 }
1288 
1290 {
1291  return id.rawId() & 0x3FFF; // Get 14 least-significant digits
1292 }
1293 
1294 
1295 //define this as a plug-in
TH1D * h_PassSelPF_
float tpfjet_unkown_pz_
RunNumber_t run() const
Definition: EventID.h:39
std::vector< int > tpfjet_twr_elmttype_
float tpfjet_photon_py_
type
Definition: HCALResponse.h:21
T getParameter(std::string const &) const
EventNumber_t event() const
Definition: EventID.h:41
std::vector< float > tpfjet_had_emf_
T getUntrackedParameter(std::string const &, T const &) const
std::vector< float > tpfjet_had_E_
int i
Definition: DBlmapReader.cc:9
std::vector< int > ppfjet_had_candtrackind_
const reco::PFJet * jet(void) const
Definition: DiJetAnalyzer.h:72
std::vector< float > tpfjet_candtrack_pz_
std::vector< float > tpfjet_cluster_eta_
std::vector< float > tpfjet_twr_hade_
ParticleType
particle types
Definition: PFCandidate.h:44
std::vector< float > tpfjet_twr_dR_
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:252
double minTagJetEta_
std::vector< int > ppfjet_twr_depth_
std::vector< float > tpfjet_cluster_dR_
std::vector< float > ppfjet_twr_frac_
float tpfjet_muon_E_
Particle flow cluster, see clustering algorithm in PFClusterAlgo.
Definition: PFCluster.h:47
float ppfjet_photon_py_
virtual double energy() const final
energy
float pf_thirdjet_py_
std::vector< float > ppfjet_candtrack_py_
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:462
float ppfjet_muon_py_
std::vector< int > tpfjet_had_candtrackind_
double maxDeltaEta_
virtual double correction(const LorentzVector &fJet) const =0
get correction using Jet information only
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
float ppfjet_photon_pz_
std::vector< float > tpfjet_had_py_
float chargedHadronEnergyFraction() const
chargedHadronEnergyFraction
Definition: PFJet.h:100
Base class for all types of Jets.
Definition: Jet.h:20
float tpfjet_photon_EcalE_
std::vector< float > ppfjet_had_emf_
std::vector< int > ppfjet_twr_subdet_
float tpfjet_muon_EcalE_
size_type size() const
Definition: OwnVector.h:264
float tpfjet_electron_px_
virtual void endJob()
float tpfjet_photon_px_
std::vector< int > tpfjet_twr_depth_
float pf_realthirdjet_px_
virtual double phi() const final
momentum azimuthal angle
std::string hfRecHitName_
Definition: DiJetAnalyzer.h:98
std::vector< int > ppfjet_twr_clusterind_
const std::vector< std::pair< DetId, float > > & hitsAndFractions() const
Definition: CaloCluster.h:192
std::vector< float > tpfjet_candtrack_EcalE_
TFile * rootfile_
double px() const
x coordinate of momentum vector
Definition: TrackBase.h:622
std::vector< int > tpfjet_twr_clusterind_
float pf_realthirdjet_scale_
std::vector< int > ppfjet_twr_candtrackind_
std::vector< float > ppfjet_twr_hade_
edm::EDGetTokenT< reco::VertexCollection > tok_Vertex_
LuminosityBlockNumber_t luminosityBlock() const
Definition: EventID.h:40
float tpfjet_unkown_py_
int chargedMultiplicity() const
chargedMultiplicity
Definition: PFJet.h:155
Jets made from PFObjects.
Definition: PFJet.h:21
float tpfjet_EMfrac_
double eta() const
pseudorapidity of cluster centroid
Definition: CaloCluster.h:163
float ppfjet_EMfrac_
virtual const CaloCellGeometry * getGeometry(const DetId &id) const
Get the cell geometry of a given detector id. Should return false if not found.
float ppfjet_unkown_px_
std::string hbheRecHitName_
Definition: DiJetAnalyzer.h:97
float ppfjet_electron_py_
dictionary elements
float ppfjet_muon_E_
U second(std::pair< T, U > const &p)
std::string pfJetCorrName_
Definition: DiJetAnalyzer.h:96
virtual void analyze(const edm::Event &, const edm::EventSetup &)
std::vector< float > ppfjet_cluster_eta_
float tpfjet_muon_py_
float neutralEmEnergyFraction() const
neutralEmEnergyFraction
Definition: PFJet.h:152
std::vector< float > tpfjet_had_px_
std::string pvCollName_
int iEvent
Definition: GenABIO.cc:230
float ppfjet_unkown_pz_
edm::EDGetTokenT< edm::SortedCollection< HFRecHit, edm::StrictWeakOrdering< HFRecHit > > > tok_HF_
float tpfjet_muon_px_
std::vector< float > ppfjet_candtrack_pz_
std::vector< float > tpfjet_had_EcalE_
edm::EDGetTokenT< reco::PFJetCollection > tok_PFJet_
virtual void beginJob()
float tpfjet_photon_E_
T sqrt(T t)
Definition: SSEVec.h:18
float ppfjet_electron_pz_
std::vector< float > tpfjet_had_pz_
std::vector< int > tpfjet_twr_ieta_
std::vector< float > tpfjet_cluster_phi_
float ppfjet_unkown_EcalE_
std::vector< float > ppfjet_cluster_dR_
float neutralHadronEnergyFraction() const
neutralHadronEnergyFraction
Definition: PFJet.h:104
virtual double py() const final
y coordinate of momentum vector
std::vector< float > ppfjet_had_rawHcalE_
bool isValid() const
Definition: HandleBase.h:75
std::vector< float > tpfjet_candtrack_px_
std::vector< float > ppfjet_cluster_phi_
float pf_realthirdjet_py_
int neutralMultiplicity() const
neutralMultiplicity
Definition: PFJet.h:157
dictionary cv
Definition: cuy.py:362
double pz() const
z coordinate of momentum vector
Definition: TrackBase.h:634
std::vector< int > tpfjet_twr_candtrackind_
std::vector< float > ppfjet_had_px_
std::vector< int > ppfjet_twr_hadind_
std::vector< float > ppfjet_had_pz_
float ppfjet_muon_px_
DiJetAnalyzer(const edm::ParameterSet &)
Definition: DiJetAnalyzer.cc:9
Definition: DetId.h:18
std::vector< int > tpfjet_twr_iphi_
float ppfjet_hadEcalEfrac_
float ppfjet_muon_pz_
std::vector< int > ppfjet_twr_ieta_
float chargedEmEnergyFraction() const
chargedEmEnergyFraction
Definition: PFJet.h:144
float tpfjet_muon_pz_
float tpfjet_electron_EcalE_
float tpfjet_unkown_E_
std::string hoRecHitName_
Definition: DiJetAnalyzer.h:99
edm::EDGetTokenT< edm::SortedCollection< HORecHit, edm::StrictWeakOrdering< HORecHit > > > tok_HO_
double maxThirdJetEt_
const T & get() const
Definition: EventSetup.h:56
static const JetCorrector * getJetCorrector(const std::string &fName, const edm::EventSetup &fSetup)
retrieve corrector from the event setup. troughs exception if something is missing ...
Definition: JetCorrector.cc:50
double maxTagJetEta_
float tpfjet_hadEcalEfrac_
std::vector< float > ppfjet_had_E_
double scale(void) const
Definition: DiJetAnalyzer.h:74
std::vector< float > ppfjet_twr_dR_
float pf_thirdjet_px_
virtual double p() const final
magnitude of momentum vector
float tpfjet_electron_pz_
std::vector< float > tpfjet_had_rawHcalE_
float ppfjet_electron_px_
std::vector< float > ppfjet_had_py_
float tpfjet_electron_py_
edm::EventID id() const
Definition: EventBase.h:59
float ppfjet_muon_EcalE_
virtual float jetArea() const
get jet area
Definition: Jet.h:105
float tpfjet_unkown_EcalE_
float pf_dijet_balance_
std::string rootHistFilename_
float tpfjet_photon_pz_
double deltaR(const reco::Jet *j1, const reco::Jet *j2)
float pf_dijet_deta_
float ppfjet_photon_E_
std::vector< int > ppfjet_had_id_
virtual std::vector< reco::PFCandidatePtr > getPFConstituents() const
get all constituents
Definition: PFJet.cc:52
virtual double px() const final
x coordinate of momentum vector
virtual double et() const final
transverse energy
int getEtaPhi(const DetId id)
std::vector< int > ppfjet_had_ntwrs_
tuple cout
Definition: gather_cfg.py:145
std::vector< float > tpfjet_twr_frac_
virtual double eta() const final
momentum pseudorapidity
float tpfjet_unkown_px_
float ppfjet_photon_px_
std::vector< float > tpfjet_candtrack_py_
float tpfjet_electron_E_
std::vector< int > tpfjet_had_id_
double phi() const
azimuthal angle of cluster centroid
Definition: CaloCluster.h:166
std::vector< float > ppfjet_had_EcalE_
std::vector< float > ppfjet_candtrack_px_
float ppfjet_electron_EcalE_
std::vector< int > tpfjet_had_ntwrs_
std::vector< int > ppfjet_twr_elmttype_
const CornersVec & getCorners() const
Returns the corner points of this cell&#39;s volume.
std::vector< int > tpfjet_twr_subdet_
std::string pfJetCollName_
Definition: DiJetAnalyzer.h:95
double minSumJetEt_
std::vector< float > ppfjet_candtrack_EcalE_
float ppfjet_photon_EcalE_
float ppfjet_unkown_E_
std::vector< int > tpfjet_twr_hadind_
float ppfjet_electron_E_
double py() const
y coordinate of momentum vector
Definition: TrackBase.h:628
float ppfjet_unkown_py_
virtual double pt() const final
transverse momentum
edm::EDGetTokenT< edm::SortedCollection< HBHERecHit, edm::StrictWeakOrdering< HBHERecHit > > > tok_HBHE_
std::vector< int > ppfjet_twr_iphi_
float pf_dijet_dphi_