CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
IsolatedTracksCone.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: IsolatedTracksCone
4 // Class: IsolatedTracksCone
5 //
6 
15 //
16 // Original Author: Jim Hirschauer (adaptation of Seema Sharma's
17 // IsolatedTracksNew)
18 // Created: Thu Nov 6 15:30:40 CST 2008
19 //
20 //
21 #define _DEBUG_QUIET
22 
25 
32 
36 
38 
39  //now do what ever initialization is needed
40  doMC = iConfig.getUntrackedParameter<bool> ("DoMC", false);
41  myverbose_ =
42  iConfig.getUntrackedParameter<int>( "Verbosity", 5 );
44  iConfig.getUntrackedParameter<bool>( "useJetTrigger", false);
46  iConfig.getUntrackedParameter<double>( "drLeadJetVeto", 1.2 );
47  ptMinLeadJet_ =
48  iConfig.getUntrackedParameter<double>( "ptMinLeadJet", 15.0 );
49 
50  debugTrks_ =
51  iConfig.getUntrackedParameter<int>("DebugTracks");
53  iConfig.getUntrackedParameter<bool>("PrintTrkHitPattern");
54 
55  minTrackP_ =
56  iConfig.getUntrackedParameter<double>( "minTrackP", 10.0);
57  maxTrackEta_ =
58  iConfig.getUntrackedParameter<double>( "maxTrackEta", 5.0);
60  iConfig.getUntrackedParameter<double>( "maxNearTrackP", 1.0);
61 
63  iConfig.getUntrackedParameter<int>("DebugEcalSimInfo");
64 
66  iConfig.getUntrackedParameter<bool>("ApplyEcalIsolation");
67 
69  consumes<l1extra::L1JetParticleCollection>(iConfig.getParameter<edm::InputTag>("L1extraTauJetSource"));
71  consumes<l1extra::L1JetParticleCollection>(iConfig.getParameter<edm::InputTag>("L1extraCenJetSource"));
73  consumes<l1extra::L1JetParticleCollection>(iConfig.getParameter<edm::InputTag>("L1extraFwdJetSource"));
74 
75  // hard coded collection access
76  tok_EB_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit","EcalRecHitsEB"));
77  tok_EE_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit","EcalRecHitsEE"));
78  tok_hbhe_ = consumes<HBHERecHitCollection>(edm::InputTag("hbhereco"));
79  tok_genTrack_ = consumes<reco::TrackCollection>(edm::InputTag("generalTracks"));
80  tok_simTk_ = consumes<edm::SimTrackContainer>(edm::InputTag("g4SimHits"));
81  tok_simVtx_ = consumes<edm::SimVertexContainer>(edm::InputTag("g4SimHits"));
82  tok_caloEB_ = consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", "EcalHitsEB"));
83  tok_caloEE_ = consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", "EcalHitsEE"));
84  tok_caloHH_ = consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", "HcalHits"));
85  tok_trigger_ = consumes<edm::TriggerResults>(edm::InputTag("TriggerResults","","HLT"));
86 
87 
89  iConfig.getParameter<edm::ParameterSet>("TrackAssociatorParameters");
91  parameters_.loadParameters( parameters, iC );
94 
95  if(myverbose_>=0) {
96  std::cout <<"Parameters read from config file \n"
97  << "myverbose_ " << myverbose_ << "\n"
98  << "useJetTrigger_ " << useJetTrigger_ << "\n"
99  << "drLeadJetVeto_ " << drLeadJetVeto_ << "\n"
100  << "minTrackP_ " << minTrackP_ << "\n"
101  << "maxTrackEta_ " << maxTrackEta_ << "\n"
102  << "maxNearTrackP_ " << maxNearTrackP_
103  << std::endl;
104  }
105 }
106 
108  delete trackAssociator_;
109 }
110 
112  const edm::EventSetup& iSetup) {
113 
114  unsigned int irun = (unsigned int)iEvent.id().run();
115  unsigned int ilum = (unsigned int)iEvent.getLuminosityBlock().luminosityBlock();
116  unsigned int ievt = (unsigned int)iEvent.id().event();
117 
118 
119 
121 
122  // check the L1 objects
123  bool L1Pass = false;
124  leadL1JetPT=-999, leadL1JetEta=-999, leadL1JetPhi=-999;
125  if( !useJetTrigger_ ) {
126  L1Pass = true;
127  } else {
129  iEvent.getByToken(tok_L1extTauJet_,l1TauHandle);
130  l1extra::L1JetParticleCollection::const_iterator itr;
131  for(itr = l1TauHandle->begin(); itr != l1TauHandle->end(); ++itr )
132  {
133  if( itr->pt()>leadL1JetPT ) {
134  leadL1JetPT = itr->pt();
135  leadL1JetEta = itr->eta();
136  leadL1JetPhi = itr->phi();
137  }
138  }
140  iEvent.getByToken(tok_L1extCenJet_,l1CenJetHandle);
141  for( itr = l1CenJetHandle->begin(); itr != l1CenJetHandle->end(); ++itr )
142  {
143  if( itr->pt()>leadL1JetPT ) {
144  leadL1JetPT = itr->pt();
145  leadL1JetEta = itr->eta();
146  leadL1JetPhi = itr->phi();
147  }
148  }
150  iEvent.getByToken(tok_L1extFwdJet_,l1FwdJetHandle);
151  for( itr = l1FwdJetHandle->begin(); itr != l1FwdJetHandle->end(); ++itr )
152  {
153  if( itr->pt()>leadL1JetPT ) {
154  leadL1JetPT = itr->pt();
155  leadL1JetEta = itr->eta();
156  leadL1JetPhi = itr->phi();
157  }
158  }
160  {
161  L1Pass = true;
162  }
163 
164  }
165 
166 
168  // Break now if L1Pass is false
170  if (!L1Pass) {
171  nEVT_failL1++;
172  // // std::cout << "L1Pass is false : " << L1Pass << std::endl;
173  // return;
174  }
175 
177  // Get the collection handles
179 
180 
182  iSetup.get<CaloGeometryRecord>().get(pG);
183  const CaloGeometry* geo = pG.product();
184  const CaloSubdetectorGeometry* gEB =
186  const CaloSubdetectorGeometry* gEE =
188  const CaloSubdetectorGeometry* gHB =
190  const CaloSubdetectorGeometry* gHE =
192 
193  edm::ESHandle<CaloTopology> theCaloTopology;
194  iSetup.get<CaloTopologyRecord>().get(theCaloTopology);
195  const CaloTopology *caloTopology = theCaloTopology.product();
196  // const CaloSubdetectorTopology* theEBTopology = theCaloTopology->getSubdetectorTopology(DetId::Ecal,EcalBarrel);
197  // const CaloSubdetectorTopology* theEETopology = theCaloTopology->getSubdetectorTopology(DetId::Ecal,EcalEndcap);
198 
200  iSetup.get<IdealGeometryRecord>().get(htopo);
201  const HcalTopology* theHBHETopology = htopo.product();
202 
203  edm::Handle<EcalRecHitCollection> barrelRecHitsHandle;
204  edm::Handle<EcalRecHitCollection> endcapRecHitsHandle;
205  iEvent.getByToken(tok_EB_,barrelRecHitsHandle);
206  iEvent.getByToken(tok_EE_,endcapRecHitsHandle);
207 
208  // Retrieve the good/bad ECAL channels from the DB
210  iSetup.get<EcalChannelStatusRcd>().get(ecalChStatus);
211  const EcalChannelStatus* theEcalChStatus = ecalChStatus.product();
212 
214  iEvent.getByToken(tok_hbhe_,hbhe);
215  const HBHERecHitCollection Hithbhe = *(hbhe.product());
216 
218  iEvent.getByToken(tok_genTrack_, trkCollection);
219  reco::TrackCollection::const_iterator trkItr;
220  if(debugTrks_>1){
221  std::cout << "Track Collection: " << std::endl;
222  std::cout << "Number of Tracks " << trkCollection->size() << std::endl;
223  }
224  std::string theTrackQuality = "highPurity";
225  reco::TrackBase::TrackQuality trackQuality_=
226  reco::TrackBase::qualityByName(theTrackQuality);
227 
228  //get Handles to SimTracks and SimHits
230  if (doMC) iEvent.getByToken(tok_simTk_,SimTk);
231  edm::SimTrackContainer::const_iterator simTrkItr;
232 
234  if (doMC) iEvent.getByToken(tok_simVtx_,SimVtx);
235 
236  //get Handles to PCaloHitContainers of eb/ee/hbhe
238  if (doMC) iEvent.getByToken(tok_caloEB_, pcaloeb);
239 
241  if (doMC) iEvent.getByToken(tok_caloEE_, pcaloee);
242 
244  if (doMC) iEvent.getByToken(tok_caloHH_, pcalohh);
245 
246 
247 
249  // Get HLT_IsoTrackHB/HE Information
251 
253  iEvent.getByToken( tok_trigger_, triggerResults);
254 
255 
256 
257  std::vector<int> v_hlTriggers;
258  int hltHB(-99);
259  int hltHE(-99);
260  int hltL1Jet15 (-99);
261  int hltJet30 (-99);
262  int hltJet50 (-99);
263  int hltJet80 (-99);
264  int hltJet110 (-99);
265  int hltJet140 (-99);
266  int hltJet180 (-99);
267  int hltL1SingleEG5 (-99);
268  int hltZeroBias (-99);
269  int hltMinBiasHcal (-99);
270  int hltMinBiasEcal (-99);
271  int hltMinBiasPixel (-99);
272  int hltSingleIsoTau30_Trk5 (-99);
273  int hltDoubleLooseIsoTau15_Trk5(-99);
274 
275  if (triggerResults.isValid()) {
276 
277  const edm::TriggerNames & triggerNames = iEvent.triggerNames(*triggerResults);
278  // TriggerNames class triggerNames.init(*triggerResults);
279 
280 
281  for (unsigned int i=0; i<triggerResults->size(); i++){
282  // std::cout << "triggerNames.triggerName(" << i << ") = " << triggerNames.triggerName(i) << std::endl;
283  if (triggerNames.triggerName(i) == "HLT_IsoTrackHE_1E31") hltHE = triggerResults->accept(i);
284  if (triggerNames.triggerName(i) == "HLT_IsoTrackHB_1E31") hltHB = triggerResults->accept(i);
285  if (triggerNames.triggerName(i) == "HLT_L1Jet15" ) hltL1Jet15 = triggerResults->accept(i);
286  if (triggerNames.triggerName(i) == "HLT_Jet30" ) hltJet30 = triggerResults->accept(i);
287  if (triggerNames.triggerName(i) == "HLT_Jet50" ) hltJet50 = triggerResults->accept(i);
288  if (triggerNames.triggerName(i) == "HLT_Jet80" ) hltJet80 = triggerResults->accept(i);
289  if (triggerNames.triggerName(i) == "HLT_Jet110" ) hltJet110 = triggerResults->accept(i);
290  if (triggerNames.triggerName(i) == "HLT_Jet140" ) hltJet140 = triggerResults->accept(i);
291  if (triggerNames.triggerName(i) == "HLT_Jet180" ) hltJet180 = triggerResults->accept(i);
292  if (triggerNames.triggerName(i) == "HLT_L1SingleEG5" ) hltL1SingleEG5 = triggerResults->accept(i);
293  if (triggerNames.triggerName(i) == "HLT_ZeroBias" ) hltZeroBias = triggerResults->accept(i);
294  if (triggerNames.triggerName(i) == "HLT_MinBiasHcal" ) hltMinBiasHcal = triggerResults->accept(i);
295  if (triggerNames.triggerName(i) == "HLT_MinBiasEcal" ) hltMinBiasEcal = triggerResults->accept(i);
296  if (triggerNames.triggerName(i) == "HLT_MinBiasPixel" ) hltMinBiasPixel = triggerResults->accept(i);
297  if (triggerNames.triggerName(i) == "HLT_SingleIsoTau30_Trk5" ) hltSingleIsoTau30_Trk5 = triggerResults->accept(i);
298  if (triggerNames.triggerName(i) == "HLT_DoubleLooseIsoTau15_Trk5" ) hltDoubleLooseIsoTau15_Trk5 = triggerResults->accept(i);
299  }
300  }
301 
302 
303 
304 
306  // Primary loop over tracks
308  TrackerHitAssociator* associate=0;
309  if (doMC) associate = new TrackerHitAssociator(iEvent);
310 
311 
312  nTRK = 0;
313  nRawTRK = 0;
314  nFailPt = 0;
315  nFailEta = 0;
317  nMissEcal = 0;
318  nMissHcal = 0;
319 
320  for( trkItr = trkCollection->begin();
321  trkItr != trkCollection->end(); ++trkItr)
322  {
323 
324  nRawTRK++;
325 
326  const reco::Track* pTrack = &(*trkItr);
327 
329  // Check for min Pt and max Eta P
331 
332  bool trkQual = pTrack->quality(trackQuality_);
333  bool goodPt = pTrack->p()>minTrackP_;
334  bool goodEta = std::abs(pTrack->momentum().eta())<maxTrackEta_;
335 
336  double eta1 = pTrack->momentum().eta();
337  double phi1 = pTrack->momentum().phi();
338  double pt1 = pTrack->pt();
339  double p1 = pTrack->p();
340 
341 
342  if (!goodEta){
343  nFailEta++;
344  }
345  if (!goodPt){
346  nFailPt++;
347  }
348  if (!trkQual){
350  }
351 
352  hRawPt ->Fill(pt1 );
353  hRawP ->Fill(p1 );
354  hRawEta->Fill(eta1);
355  hRawPhi->Fill(phi1);
356 
357  if( !goodEta || !goodPt || !trkQual ) continue; // Skip to next track
358 
360  // Find track trajectory
362 
363 
364  const FreeTrajectoryState fts1 =
365  trackAssociator_->getFreeTrajectoryState(iSetup, *pTrack);
366 
367 
368  TrackDetMatchInfo info1 =
369  trackAssociator_->associate(iEvent, iSetup, fts1, parameters_);
370 
371 
372 
374  // First confirm track makes it to Hcal
376 
377  if (info1.trkGlobPosAtHcal.x()==0 &&
378  info1.trkGlobPosAtHcal.y()==0 &&
379  info1.trkGlobPosAtHcal.z()==0)
380  {
381  nMissHcal++;
382  continue;
383  }
384 
385  const GlobalPoint hpoint1(info1.trkGlobPosAtHcal.x(),
386  info1.trkGlobPosAtHcal.y(),
387  info1.trkGlobPosAtHcal.z());
388 
389 
390 
392  // Get basic quantities
394 
395  const reco::HitPattern& hitp = pTrack->hitPattern();
396  int nLayersCrossed = hitp.trackerLayersWithMeasurement();
397  int nOuterHits = hitp.stripTOBLayersWithMeasurement()
398  +hitp.stripTECLayersWithMeasurement() ;
399 
400 
401  double simP = 0;
402  if (doMC) {
403  edm::SimTrackContainer::const_iterator matchedSimTrk =
404  spr::matchedSimTrack(iEvent, SimTk, SimVtx, pTrack, *associate, false);
405  simP = matchedSimTrk->momentum().P();
406  }
408  // Get Ecal Point
410 
411  const GlobalPoint point1(info1.trkGlobPosAtEcal.x(),
412  info1.trkGlobPosAtEcal.y(),
413  info1.trkGlobPosAtEcal.z());
414 
415 
416  // Sanity check that track hits Ecal
417 
418  if (info1.trkGlobPosAtEcal.x()==0 &&
419  info1.trkGlobPosAtEcal.y()==0 &&
420  info1.trkGlobPosAtEcal.z()==0)
421  {
422  std::cout << "Track doesn't reach Ecal." << std::endl;
423  nMissEcal++;
424  continue;
425  }
426 
427  // Get Track Momentum - make sure you have latest version of TrackDetMatchInfo
428 
429  GlobalVector trackMomAtEcal = info1.trkMomAtEcal;
430  GlobalVector trackMomAtHcal = info1.trkMomAtHcal;
431 
433  // If using Jet trigger, get distance from leading jet
435 
436  double drFromLeadJet = 999.0;
437  if( useJetTrigger_ ) {
438  double dphi = DeltaPhi(phi1, leadL1JetPhi);
439  double deta = eta1 - leadL1JetEta;
440  drFromLeadJet = sqrt(dphi*dphi + deta*deta);
441  }
442 
443 
445  // Define Arrays for sizes of Charge, Neut Iso Radii and
446  // Clustering Cone Radius.
448 
449  const int a_size = 7;
450  double a_coneR[a_size];
451  double a_charIsoR[a_size];
452  double a_neutIsoR[a_size];
453 
454  a_coneR[0] = 17.49; // = area of 2x2
455  a_coneR[1] = 26.23; // = area of 3x3
456  a_coneR[2] = 30.61;
457  a_coneR[3] = 34.98; // = area of 4x4
458  a_coneR[4] = 39.35;
459  a_coneR[5] = 43.72; // = area of 5x5
460  a_coneR[6] = 52.46; // = area of 6x6
461 
462  for (int i=0; i<a_size; i++){
463  a_charIsoR[i] = a_coneR[i]+28.9; // 28.9 gives 55.1 for 3x3 benchmark
464  a_neutIsoR[i] = a_charIsoR[i]*0.726; // Ecal radius = 0.726*Hcal radius
465  }
466 
468  // Do Neutral Iso in radius on Ecal surface.
470 
471  // NxN cluster
472  double e3x3=-999.0;
473  double trkEcalEne =-999.0;
475  iSetup.get<EcalSeverityLevelAlgoRcd>().get(sevlv);
476 
477  if(std::abs(point1.eta())<1.479) {
478  const DetId isoCell = gEB->getClosestCell(point1);
479  e3x3 = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology, sevlv.product(),1,1).first;
480  trkEcalEne = spr::eCaloSimInfo(iEvent, geo, pcaloeb, pcaloee, SimTk, SimVtx, pTrack, *associate);
481  } else {
482  const DetId isoCell = gEE->getClosestCell(point1);
483  e3x3 = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology, sevlv.product(),1,1).first;
484  trkEcalEne = spr::eCaloSimInfo(iEvent, geo, pcaloeb, pcaloee, SimTk, SimVtx, pTrack, *associate);
485  }
486 
487  // Cone cluster
488 
489  // Set up array of cone sizes for MIP cut
490  const int a_mip_size = 5;
491  double a_mipR[a_mip_size];
492  a_mipR[0] = 3.84; // = area of 3x3 ecal
493  a_mipR[1] = 14.0;
494  a_mipR[2] = 19.0;
495  a_mipR[3] = 24.0;
496  a_mipR[4] = 9.0; // = from standard analyzer
497 
498  std::vector<double> v_eDR;
499  for (int i = 0 ; i < a_size ; i++){
500  int nRH_eDR = 0;
501 
502  // Cone in ecal
503  double eDR = spr::eCone_ecal(geo,
504  barrelRecHitsHandle,
505  endcapRecHitsHandle,
506  hpoint1, point1,
507  a_neutIsoR[i],
508  trackMomAtEcal, nRH_eDR);
509  v_eDR.push_back(eDR);
510 
511  }
512 
513  std::vector<double> v_eMipDR;
514  for (int i = 0 ; i < a_mip_size ; i++){
515  int nRH_eMipDR = 0;
516  double eMipDR = spr::eCone_ecal(geo, barrelRecHitsHandle,
517  endcapRecHitsHandle,
518  hpoint1, point1,
519  a_mipR[i], trackMomAtEcal, nRH_eMipDR);
520 
521  v_eMipDR.push_back(eMipDR);
522  }
523 
524 
526  // Do charge isolation in radius at Hcal surface for 5 different
527  // radii defined above in a_charIso
529 
530  std::vector<double> v_hmaxNearP_goodTrk;
531  std::vector<double> v_hmaxNearP ;
532  std::vector<int> v_hnNearTRKs ;
533  std::vector<int> v_hnLayers_maxNearP;
534  std::vector<int> v_htrkQual_maxNearP;
535 
536  std::vector<double> v_cone_hmaxNearP_goodTrk;
537  std::vector<double> v_cone_hmaxNearP ;
538  std::vector<int> v_cone_hnNearTRKs ;
539  std::vector<int> v_cone_hnLayers_maxNearP;
540  std::vector<int> v_cone_htrkQual_maxNearP;
541 
542  for (int i = 0 ; i < a_size ; i++){
543 
544  double hmaxNearP = -999.0;
545  int hnNearTRKs = 0;
546  int hnLayers_maxNearP = 0;
547  int htrkQual_maxNearP = -1;
548  double hmaxNearP_goodTrk = -999.0;
549 
550  double conehmaxNearP = -999.0;
551  int conehnNearTRKs = 0;
552  int conehnLayers_maxNearP = 0;
553  int conehtrkQual_maxNearP = -1;
554  double conehmaxNearP_goodTrk = -999.0;
555 
556  conehmaxNearP = spr::coneChargeIsolation(iEvent, iSetup,
557  trkItr, trkCollection,
559  theTrackQuality, conehnNearTRKs,
560  conehnLayers_maxNearP,
561  conehtrkQual_maxNearP,
562  conehmaxNearP_goodTrk,
563  hpoint1, trackMomAtHcal,
564  a_charIsoR[i]);
565 
566  v_hmaxNearP_goodTrk.push_back(hmaxNearP_goodTrk);
567  v_hmaxNearP .push_back(hmaxNearP );
568  v_hnNearTRKs .push_back(hnNearTRKs );
569  v_hnLayers_maxNearP.push_back(hnLayers_maxNearP);
570  v_htrkQual_maxNearP.push_back(htrkQual_maxNearP);
571 
572  v_cone_hmaxNearP_goodTrk.push_back(conehmaxNearP_goodTrk);
573  v_cone_hmaxNearP .push_back(conehmaxNearP );
574  v_cone_hnNearTRKs .push_back(conehnNearTRKs );
575  v_cone_hnLayers_maxNearP.push_back(conehnLayers_maxNearP);
576  v_cone_htrkQual_maxNearP.push_back(conehtrkQual_maxNearP);
577 
578  }
579 
580  double h3x3=-999.0, h5x5=-999.0;
581  double hsim3x3=-999.0, hsim5x5=-999.0, trkHcalEne=-999.0;
582  std::map<std::string, double> hsimInfo3x3, hsimInfo5x5;
583  double distFromHotCell_h3x3 = -99.;
584  int ietaFromHotCell_h3x3 = -99;
585  int iphiFromHotCell_h3x3 = -99;
586  double distFromHotCell_h5x5 = -99.;
587  int ietaFromHotCell_h5x5 = -99;
588  int iphiFromHotCell_h5x5 = -99;
589 
590  GlobalPoint gPosHotCell_h3x3(0.,0.,0.);
591  GlobalPoint gPosHotCell_h5x5(0.,0.,0.);
592 
593  int nRH_h3x3(0), nRH_h5x5(0);
594 
595  // Hcal Energy Clustering
596 
597  // Get closetcell for ietaFromHotCell and iphiFromHotCell
598  DetId ClosestCell;
599  if( std::abs(pTrack->eta())<1.392 ) {
600  ClosestCell = gHB->getClosestCell(hpoint1);
601  } else {
602  ClosestCell = gHE->getClosestCell(hpoint1);
603  }
604  // Transform into HcalDetId so that I can get ieta, iphi later.
605  HcalDetId ClosestCell_HcalDetId(ClosestCell.rawId());
606 
607  // Using NxN Cluster
608  std::vector<int> v_RH_h3x3_ieta;
609  std::vector<int> v_RH_h3x3_iphi;
610  std::vector<double> v_RH_h3x3_ene;
611  std::vector<int> v_RH_h5x5_ieta;
612  std::vector<int> v_RH_h5x5_iphi;
613  std::vector<double> v_RH_h5x5_ene;
614 
615 
616  h3x3 = spr::eHCALmatrix(geo, theHBHETopology, ClosestCell, hbhe,1,1,
617  nRH_h3x3, v_RH_h3x3_ieta, v_RH_h3x3_iphi, v_RH_h3x3_ene,
618  gPosHotCell_h3x3);
619  distFromHotCell_h3x3 = spr::getDistInPlaneTrackDir(hpoint1, trackMomAtHcal, gPosHotCell_h3x3);
620 
621  h5x5 = spr::eHCALmatrix(geo, theHBHETopology, ClosestCell, hbhe,2,2,
622  nRH_h5x5, v_RH_h5x5_ieta, v_RH_h5x5_iphi, v_RH_h5x5_ene,
623  gPosHotCell_h5x5);
624  distFromHotCell_h5x5 = spr::getDistInPlaneTrackDir(hpoint1, trackMomAtHcal, gPosHotCell_h5x5);
625 
626 
627  // double heta = (double)hpoint1.eta();
628  // double hphi = (double)hpoint1.phi();
629  std::vector<int> multiplicity3x3;
630  std::vector<int> multiplicity5x5;
631  if (doMC) {
632  hsim3x3 = spr::eHCALmatrix(theHBHETopology, ClosestCell,
633  pcalohh,1,1);
634  hsim5x5 = spr::eHCALmatrix(theHBHETopology, ClosestCell,
635  pcalohh,2,2);
636 
637  hsimInfo3x3 = spr::eHCALSimInfo(iEvent, theHBHETopology, ClosestCell, pcalohh, SimTk, SimVtx, pTrack, *associate, 1,1, multiplicity3x3);
638  hsimInfo5x5 = spr::eHCALSimInfo(iEvent, theHBHETopology, ClosestCell, pcalohh, SimTk, SimVtx, pTrack, *associate, 2,2, multiplicity5x5);
639 
640  // Get energy from all simhits in hcal associated with iso track
641  trkHcalEne = spr::eCaloSimInfo(iEvent, geo, pcalohh, SimTk, SimVtx, pTrack, *associate);
642  }
643 
644  // Finally for cones of varying radii.
645  std::vector<double> v_hsimInfoConeMatched;
646  std::vector<double> v_hsimInfoConeRest ;
647  std::vector<double> v_hsimInfoConePhoton ;
648  std::vector<double> v_hsimInfoConeNeutHad;
649  std::vector<double> v_hsimInfoConeCharHad;
650  std::vector<double> v_hsimInfoConePdgMatched;
651  std::vector<double> v_hsimInfoConeTotal ;
652 
653  std::vector<int> v_hsimInfoConeNMatched;
654  std::vector<int> v_hsimInfoConeNTotal ;
655  std::vector<int> v_hsimInfoConeNNeutHad;
656  std::vector<int> v_hsimInfoConeNCharHad;
657  std::vector<int> v_hsimInfoConeNPhoton ;
658  std::vector<int> v_hsimInfoConeNRest ;
659 
660  std::vector<double> v_hsimCone ;
661  std::vector<double> v_hCone ;
662 
663  std::vector<int> v_nRecHitsCone ;
664  std::vector<int> v_nSimHitsCone ;
665 
666  std::vector<double> v_distFromHotCell;
667  std::vector<int> v_ietaFromHotCell;
668  std::vector<int> v_iphiFromHotCell;
669  GlobalPoint gposHotCell(0.,0.,0.);
670 
671 
672  std::vector<int> v_RH_r26_ieta;
673  std::vector<int> v_RH_r26_iphi;
674  std::vector<double> v_RH_r26_ene;
675  std::vector<int> v_RH_r44_ieta;
676  std::vector<int> v_RH_r44_iphi;
677  std::vector<double> v_RH_r44_ene;
678 
679 
680 
681  for (int i = 0 ; i < a_size ; i++){
682 
683 
684  std::map<std::string, double> hsimInfoCone;
685  double hsimCone = -999.0, hCone = -999.0;
686  double distFromHotCell = -99.0;
687  int ietaFromHotCell = -99;
688  int iphiFromHotCell = -99;
689  int ietaHotCell = -99;
690  int iphiHotCell = -99;
691  int nRecHitsCone = -999;
692  int nSimHitsCone = -999;
693 
694  std::vector<int> multiplicityCone;
695  std::vector<DetId> coneRecHitDetIds;
696  if (doMC)
697  hsimCone = spr::eCone_hcal(geo, pcalohh, hpoint1, point1,
698  a_coneR[i], trackMomAtHcal, nSimHitsCone);
699 
700  // If needed, get ieta and iphi of rechits for cones of 23.25
701  // and for hitmap for debugging
702  bool makeHitmaps = false;
703  if (a_coneR[i] == 26.23 && makeHitmaps)
704  {
705 
706  hCone = spr::eCone_hcal(geo, hbhe, hpoint1, point1,
707  a_coneR[i], trackMomAtHcal,nRecHitsCone,
708  v_RH_r26_ieta, v_RH_r26_iphi, v_RH_r26_ene,
709  coneRecHitDetIds, distFromHotCell,
710  ietaHotCell, iphiHotCell, gposHotCell);
711  }
712  else if (a_coneR[i] == 43.72 && makeHitmaps)
713  {
714 
715  hCone = spr::eCone_hcal(geo, hbhe, hpoint1, point1,
716  a_coneR[i], trackMomAtHcal,nRecHitsCone,
717  v_RH_r44_ieta, v_RH_r44_iphi, v_RH_r44_ene,
718  coneRecHitDetIds, distFromHotCell,
719  ietaHotCell, iphiHotCell, gposHotCell);
720  }
721  else
722  {
723 
724  hCone = spr::eCone_hcal(geo, hbhe, hpoint1, point1,
725  a_coneR[i], trackMomAtHcal, nRecHitsCone,
726  coneRecHitDetIds, distFromHotCell,
727  ietaHotCell, iphiHotCell, gposHotCell);
728  }
729 
730 
731 
732  if (ietaHotCell != 99){
733  ietaFromHotCell = ietaHotCell-ClosestCell_HcalDetId.ieta();
734  iphiFromHotCell = iphiHotCell-ClosestCell_HcalDetId.iphi();
735  }
736 
737  // SimHits NOT matched to RecHits
738  if (doMC) {
739  hsimInfoCone = spr::eHCALSimInfoCone(iEvent,pcalohh, SimTk, SimVtx, pTrack, *associate, geo, hpoint1, point1, a_coneR[i], trackMomAtHcal, multiplicityCone);
740 
741 
742 
743  // SimHits matched to RecHits
744  // hsimInfoCone = spr::eHCALSimInfoCone(iEvent,pcalohh, SimTk, SimVtx,
745  // pTrack, *associate,
746  // geo, hpoint1, point1,
747  // a_coneR[i], trackMomAtHcal,
748  // multiplicityCone,
749  // coneRecHitDetIds);
750 
751  v_hsimInfoConeMatched .push_back(hsimInfoCone["eMatched" ]);
752  v_hsimInfoConeRest .push_back(hsimInfoCone["eRest" ]);
753  v_hsimInfoConePhoton .push_back(hsimInfoCone["eGamma" ]);
754  v_hsimInfoConeNeutHad .push_back(hsimInfoCone["eNeutralHad"]);
755  v_hsimInfoConeCharHad .push_back(hsimInfoCone["eChargedHad"]);
756  v_hsimInfoConePdgMatched.push_back(hsimInfoCone["pdgMatched" ]);
757  v_hsimInfoConeTotal .push_back(hsimInfoCone["eTotal" ]);
758 
759  v_hsimInfoConeNMatched .push_back(multiplicityCone.at(0));
760 
761  v_hsimInfoConeNTotal .push_back(multiplicityCone.at(1));
762  v_hsimInfoConeNNeutHad .push_back(multiplicityCone.at(2));
763  v_hsimInfoConeNCharHad .push_back(multiplicityCone.at(3));
764  v_hsimInfoConeNPhoton .push_back(multiplicityCone.at(4));
765  v_hsimInfoConeNRest .push_back(multiplicityCone.at(5));
766 
767  v_hsimCone .push_back(hsimCone );
768  v_nSimHitsCone .push_back(nSimHitsCone );
769  }
770  v_hCone .push_back(hCone );
771  v_nRecHitsCone .push_back(nRecHitsCone );
772 
773  v_distFromHotCell .push_back(distFromHotCell );
774  v_ietaFromHotCell .push_back(ietaFromHotCell );
775  v_iphiFromHotCell .push_back(iphiFromHotCell );
776 
777 
778  }
779 
780 
782  // Fill Vectors that go into root file
784 
785  t_v_hnNearTRKs ->push_back(v_hnNearTRKs );
786  t_v_hnLayers_maxNearP ->push_back(v_hnLayers_maxNearP );
787  t_v_htrkQual_maxNearP ->push_back(v_htrkQual_maxNearP );
788  t_v_hmaxNearP_goodTrk ->push_back(v_hmaxNearP_goodTrk );
789  t_v_hmaxNearP ->push_back(v_hmaxNearP );
790 
791  t_v_cone_hnNearTRKs ->push_back(v_cone_hnNearTRKs );
792  t_v_cone_hnLayers_maxNearP ->push_back(v_cone_hnLayers_maxNearP);
793  t_v_cone_htrkQual_maxNearP ->push_back(v_cone_htrkQual_maxNearP);
794  t_v_cone_hmaxNearP_goodTrk ->push_back(v_cone_hmaxNearP_goodTrk);
795  t_v_cone_hmaxNearP ->push_back(v_cone_hmaxNearP );
796 
797  // t_hScale ->push_back(hScale );
798  t_trkNOuterHits ->push_back(nOuterHits );
799  t_trkNLayersCrossed ->push_back(nLayersCrossed );
800  t_dtFromLeadJet ->push_back(drFromLeadJet );
801  t_trkP ->push_back(p1 );
802  t_trkPt ->push_back(pt1 );
803  t_trkEta ->push_back(eta1 );
804  t_trkPhi ->push_back(phi1 );
805 
806  t_e3x3 ->push_back(e3x3 );
807  t_v_eDR ->push_back(v_eDR );
808  t_v_eMipDR ->push_back(v_eMipDR );
809 
810  t_h3x3 ->push_back(h3x3 );
811  t_h5x5 ->push_back(h5x5 );
812  t_nRH_h3x3 ->push_back(nRH_h3x3 );
813  t_nRH_h5x5 ->push_back(nRH_h5x5 );
814 
815  t_v_RH_h3x3_ieta ->push_back(v_RH_h3x3_ieta);
816  t_v_RH_h3x3_iphi ->push_back(v_RH_h3x3_iphi);
817  t_v_RH_h3x3_ene ->push_back(v_RH_h3x3_ene);
818  t_v_RH_h5x5_ieta ->push_back(v_RH_h5x5_ieta);
819  t_v_RH_h5x5_iphi ->push_back(v_RH_h5x5_iphi);
820  t_v_RH_h5x5_ene ->push_back(v_RH_h5x5_ene);
821 
822  if (doMC) {
823  t_simP ->push_back(simP );
824  t_hsim3x3 ->push_back(hsim3x3 );
825  t_hsim5x5 ->push_back(hsim5x5 );
826 
827  t_hsim3x3Matched ->push_back(hsimInfo3x3["eMatched"] );
828  t_hsim5x5Matched ->push_back(hsimInfo5x5["eMatched"] );
829  t_hsim3x3Rest ->push_back(hsimInfo3x3["eRest"] );
830  t_hsim5x5Rest ->push_back(hsimInfo5x5["eRest"] );
831  t_hsim3x3Photon ->push_back(hsimInfo3x3["eGamma"] );
832  t_hsim5x5Photon ->push_back(hsimInfo5x5["eGamma"] );
833  t_hsim3x3NeutHad ->push_back(hsimInfo3x3["eNeutralHad"] );
834  t_hsim5x5NeutHad ->push_back(hsimInfo5x5["eNeutralHad"] );
835  t_hsim3x3CharHad ->push_back(hsimInfo3x3["eChargedHad"] );
836  t_hsim5x5CharHad ->push_back(hsimInfo5x5["eChargedHad"] );
837  t_hsim3x3Total ->push_back(hsimInfo3x3["eTotal"] );
838  t_hsim5x5Total ->push_back(hsimInfo5x5["eTotal"] );
839  t_hsim3x3PdgMatched ->push_back(hsimInfo3x3["pdgMatched"] );
840  t_hsim5x5PdgMatched ->push_back(hsimInfo5x5["pdgMatched"] );
841 
842  t_hsim3x3NMatched ->push_back(multiplicity3x3.at(0));
843  t_hsim3x3NTotal ->push_back(multiplicity3x3.at(1));
844  t_hsim3x3NNeutHad ->push_back(multiplicity3x3.at(2));
845  t_hsim3x3NCharHad ->push_back(multiplicity3x3.at(3));
846  t_hsim3x3NPhoton ->push_back(multiplicity3x3.at(4));
847  t_hsim3x3NRest ->push_back(multiplicity3x3.at(5));
848 
849  t_hsim5x5NMatched ->push_back(multiplicity5x5.at(0));
850  t_hsim5x5NTotal ->push_back(multiplicity5x5.at(1));
851  t_hsim5x5NNeutHad ->push_back(multiplicity5x5.at(2));
852  t_hsim5x5NCharHad ->push_back(multiplicity5x5.at(3));
853  t_hsim5x5NPhoton ->push_back(multiplicity5x5.at(4));
854  t_hsim5x5NRest ->push_back(multiplicity5x5.at(5));
855  }
856 
857  t_distFromHotCell_h3x3->push_back(distFromHotCell_h3x3);
858  t_ietaFromHotCell_h3x3->push_back(ietaFromHotCell_h3x3);
859  t_iphiFromHotCell_h3x3->push_back(iphiFromHotCell_h3x3);
860  t_distFromHotCell_h5x5->push_back(distFromHotCell_h5x5);
861  t_ietaFromHotCell_h5x5->push_back(ietaFromHotCell_h5x5);
862  t_iphiFromHotCell_h5x5->push_back(iphiFromHotCell_h5x5);
863 
864  if (doMC) {
865  t_trkHcalEne ->push_back(trkHcalEne );
866  t_trkEcalEne ->push_back(trkEcalEne );
867 
868  t_v_hsimInfoConeMatched ->push_back(v_hsimInfoConeMatched );
869  t_v_hsimInfoConeRest ->push_back(v_hsimInfoConeRest );
870  t_v_hsimInfoConePhoton ->push_back(v_hsimInfoConePhoton );
871  t_v_hsimInfoConeNeutHad ->push_back(v_hsimInfoConeNeutHad );
872  t_v_hsimInfoConeCharHad ->push_back(v_hsimInfoConeCharHad );
873  t_v_hsimInfoConePdgMatched ->push_back(v_hsimInfoConePdgMatched);
874  t_v_hsimInfoConeTotal ->push_back(v_hsimInfoConeTotal );
875 
876  t_v_hsimInfoConeNMatched ->push_back(v_hsimInfoConeNMatched );
877  t_v_hsimInfoConeNTotal ->push_back(v_hsimInfoConeNTotal );
878  t_v_hsimInfoConeNNeutHad ->push_back(v_hsimInfoConeNNeutHad );
879  t_v_hsimInfoConeNCharHad ->push_back(v_hsimInfoConeNCharHad );
880  t_v_hsimInfoConeNPhoton ->push_back(v_hsimInfoConeNPhoton );
881  t_v_hsimInfoConeNRest ->push_back(v_hsimInfoConeNRest );
882 
883  t_v_hsimCone ->push_back(v_hsimCone );
884  t_v_hCone ->push_back(v_hCone );
885  t_v_nRecHitsCone->push_back(v_nRecHitsCone);
886  t_v_nSimHitsCone->push_back(v_nSimHitsCone);
887  }
888 
889 
890  t_v_distFromHotCell->push_back(v_distFromHotCell);
891  t_v_ietaFromHotCell->push_back(v_ietaFromHotCell);
892  t_v_iphiFromHotCell->push_back(v_iphiFromHotCell);
893 
894  t_v_RH_r26_ieta ->push_back(v_RH_r26_ieta);
895  t_v_RH_r26_iphi ->push_back(v_RH_r26_iphi);
896  t_v_RH_r26_ene ->push_back(v_RH_r26_ene);
897  t_v_RH_r44_ieta ->push_back(v_RH_r44_ieta);
898  t_v_RH_r44_iphi ->push_back(v_RH_r44_iphi);
899  t_v_RH_r44_ene ->push_back(v_RH_r44_ene);
900 
901 
902 
903  t_v_hlTriggers ->push_back(v_hlTriggers);
904  t_hltHB ->push_back(hltHB);
905  t_hltHE ->push_back(hltHE);
906  t_hltL1Jet15 ->push_back(hltL1Jet15 );
907  t_hltJet30 ->push_back(hltJet30 );
908  t_hltJet50 ->push_back(hltJet50 );
909  t_hltJet80 ->push_back(hltJet80 );
910  t_hltJet110 ->push_back(hltJet110 );
911  t_hltJet140 ->push_back(hltJet140 );
912  t_hltJet180 ->push_back(hltJet180 );
913  t_hltL1SingleEG5 ->push_back(hltL1SingleEG5 );
914  t_hltZeroBias ->push_back(hltZeroBias );
915  t_hltMinBiasHcal ->push_back(hltMinBiasHcal );
916  t_hltMinBiasEcal ->push_back(hltMinBiasEcal );
917  t_hltMinBiasPixel ->push_back(hltMinBiasPixel );
918  t_hltSingleIsoTau30_Trk5 ->push_back(hltSingleIsoTau30_Trk5 );
919  t_hltDoubleLooseIsoTau15_Trk5->push_back(hltDoubleLooseIsoTau15_Trk5);
920 
921  t_irun->push_back(irun);
922  t_ievt->push_back(ievt);
923  t_ilum->push_back(ilum);
924 
925  nTRK++;
926 
927 
928  } // Loop over track collection
929 
930  // std::cout << "nEVT= " << nEVT << std::endl;
931 
932  ntp->Fill();
933  nEVT++;
934 
935 
936  delete associate;
937 }
938 
939 
940 
942 
943  // hbScale = 120.0;
944  // heScale = 135.0;
945  nEVT=0;
946  nEVT_failL1=0;
947  nTRK=0;
948 
949  double tempgen_TH[22] = { 0.0, 1.0, 2.0, 3.0, 4.0,
950  5.0, 6.0, 7.0, 8.0, 9.0,
951  10.0, 12.0, 15.0, 20.0, 25.0,
952  30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 100};
953  for(int i=0; i<22; i++) genPartPBins[i] = tempgen_TH[i];
954 
955 
956  double tempgen_Eta[5] = {0.0, 0.5, 1.1, 1.7, 2.0};
957  for(int i=0; i<5; i++) genPartEtaBins[i] = tempgen_Eta[i];
958 
959  t_v_hnNearTRKs = new std::vector<std::vector<int> > ();
960  t_v_hnLayers_maxNearP = new std::vector<std::vector<int> > ();
961  t_v_htrkQual_maxNearP = new std::vector<std::vector<int> > ();
962  t_v_hmaxNearP_goodTrk = new std::vector<std::vector<double> >();
963  t_v_hmaxNearP = new std::vector<std::vector<double> >();
964 
965  t_v_cone_hnNearTRKs = new std::vector<std::vector<int> > ();
966  t_v_cone_hnLayers_maxNearP = new std::vector<std::vector<int> > ();
967  t_v_cone_htrkQual_maxNearP = new std::vector<std::vector<int> > ();
968  t_v_cone_hmaxNearP_goodTrk = new std::vector<std::vector<double> >();
969  t_v_cone_hmaxNearP = new std::vector<std::vector<double> >();
970 
971  // t_hScale = new std::vector<double>();
972  t_trkNOuterHits = new std::vector<double>();
973  t_trkNLayersCrossed = new std::vector<double>();
974  t_dtFromLeadJet = new std::vector<double>();
975  t_trkP = new std::vector<double>();
976  t_trkPt = new std::vector<double>();
977  t_trkEta = new std::vector<double>();
978  t_trkPhi = new std::vector<double>();
979 
980  t_e3x3 = new std::vector<double>();
981  t_v_eDR = new std::vector<std::vector<double> >();
982  t_v_eMipDR = new std::vector<std::vector<double> >();
983 
984  t_h3x3 = new std::vector<double>();
985  t_h5x5 = new std::vector<double>();
986 
987  t_nRH_h3x3 = new std::vector<double>();
988  t_nRH_h5x5 = new std::vector<double>();
989 
990  if (doMC) {
991  t_simP = new std::vector<double>();
992  t_hsim3x3 = new std::vector<double>();
993  t_hsim5x5 = new std::vector<double>();
994 
995  t_hsim3x3Matched = new std::vector<double>();
996  t_hsim5x5Matched = new std::vector<double>();
997  t_hsim3x3Rest = new std::vector<double>();
998  t_hsim5x5Rest = new std::vector<double>();
999  t_hsim3x3Photon = new std::vector<double>();
1000  t_hsim5x5Photon = new std::vector<double>();
1001  t_hsim3x3NeutHad = new std::vector<double>();
1002  t_hsim5x5NeutHad = new std::vector<double>();
1003  t_hsim3x3CharHad = new std::vector<double>();
1004  t_hsim5x5CharHad = new std::vector<double>();
1005  t_hsim3x3PdgMatched = new std::vector<double>();
1006  t_hsim5x5PdgMatched = new std::vector<double>();
1007  t_hsim3x3Total = new std::vector<double>();
1008  t_hsim5x5Total = new std::vector<double>();
1009 
1010  t_hsim3x3NMatched = new std::vector<int>();
1011  t_hsim3x3NTotal = new std::vector<int>();
1012  t_hsim3x3NNeutHad = new std::vector<int>();
1013  t_hsim3x3NCharHad = new std::vector<int>();
1014  t_hsim3x3NPhoton = new std::vector<int>();
1015  t_hsim3x3NRest = new std::vector<int>();
1016 
1017  t_hsim5x5NMatched = new std::vector<int>();
1018  t_hsim5x5NTotal = new std::vector<int>();
1019  t_hsim5x5NNeutHad = new std::vector<int>();
1020  t_hsim5x5NCharHad = new std::vector<int>();
1021  t_hsim5x5NPhoton = new std::vector<int>();
1022  t_hsim5x5NRest = new std::vector<int>();
1023 
1024  t_trkHcalEne = new std::vector<double>();
1025  t_trkEcalEne = new std::vector<double>();
1026  }
1027 
1028  t_distFromHotCell_h3x3 = new std::vector<double>();
1029  t_ietaFromHotCell_h3x3 = new std::vector<int>();
1030  t_iphiFromHotCell_h3x3 = new std::vector<int>();
1031  t_distFromHotCell_h5x5 = new std::vector<double>();
1032  t_ietaFromHotCell_h5x5 = new std::vector<int>();
1033  t_iphiFromHotCell_h5x5 = new std::vector<int>();
1034 
1035  if (doMC) {
1036  t_v_hsimInfoConeMatched = new std::vector<std::vector<double> >();
1037  t_v_hsimInfoConeRest = new std::vector<std::vector<double> >();
1038  t_v_hsimInfoConePhoton = new std::vector<std::vector<double> >();
1039  t_v_hsimInfoConeNeutHad = new std::vector<std::vector<double> >();
1040  t_v_hsimInfoConeCharHad = new std::vector<std::vector<double> >();
1041  t_v_hsimInfoConePdgMatched= new std::vector<std::vector<double> >();
1042  t_v_hsimInfoConeTotal = new std::vector<std::vector<double> >();
1043 
1044  t_v_hsimInfoConeNMatched = new std::vector<std::vector<int> >();
1045  t_v_hsimInfoConeNTotal = new std::vector<std::vector<int> >();
1046  t_v_hsimInfoConeNNeutHad = new std::vector<std::vector<int> >();
1047  t_v_hsimInfoConeNCharHad = new std::vector<std::vector<int> >();
1048  t_v_hsimInfoConeNPhoton = new std::vector<std::vector<int> >();
1049  t_v_hsimInfoConeNRest = new std::vector<std::vector<int> >();
1050 
1051  t_v_hsimCone = new std::vector<std::vector<double> >();
1052  }
1053 
1054  t_v_hCone = new std::vector<std::vector<double> >();
1055  t_v_nRecHitsCone= new std::vector<std::vector<int> >();
1056  t_v_nSimHitsCone= new std::vector<std::vector<int> >();
1057 
1058  t_v_distFromHotCell= new std::vector<std::vector<double> >();
1059  t_v_ietaFromHotCell= new std::vector<std::vector<int> >();
1060  t_v_iphiFromHotCell= new std::vector<std::vector<int> >();
1061 
1062  t_v_RH_h3x3_ieta = new std::vector<std::vector<int> >();
1063  t_v_RH_h3x3_iphi = new std::vector<std::vector<int> >();
1064  t_v_RH_h3x3_ene = new std::vector<std::vector<double> >();
1065  t_v_RH_h5x5_ieta = new std::vector<std::vector<int> >();
1066  t_v_RH_h5x5_iphi = new std::vector<std::vector<int> >();
1067  t_v_RH_h5x5_ene = new std::vector<std::vector<double> >();
1068  t_v_RH_r26_ieta = new std::vector<std::vector<int> >();
1069  t_v_RH_r26_iphi = new std::vector<std::vector<int> >();
1070  t_v_RH_r26_ene = new std::vector<std::vector<double> >();
1071  t_v_RH_r44_ieta = new std::vector<std::vector<int> >();
1072  t_v_RH_r44_iphi = new std::vector<std::vector<int> >();
1073  t_v_RH_r44_ene = new std::vector<std::vector<double> >();
1074 
1075 
1076  t_v_hlTriggers = new std::vector<std::vector<int> >();
1077 
1078  t_hltHE = new std::vector<int>();
1079  t_hltHB = new std::vector<int>();
1080  t_hltL1Jet15 = new std::vector<int>();
1081  t_hltJet30 = new std::vector<int>();
1082  t_hltJet50 = new std::vector<int>();
1083  t_hltJet80 = new std::vector<int>();
1084  t_hltJet110 = new std::vector<int>();
1085  t_hltJet140 = new std::vector<int>();
1086  t_hltJet180 = new std::vector<int>();
1087  t_hltL1SingleEG5 = new std::vector<int>();
1088  t_hltZeroBias = new std::vector<int>();
1089  t_hltMinBiasHcal = new std::vector<int>();
1090  t_hltMinBiasEcal = new std::vector<int>();
1091  t_hltMinBiasPixel = new std::vector<int>();
1092  t_hltSingleIsoTau30_Trk5 = new std::vector<int>();
1093  t_hltDoubleLooseIsoTau15_Trk5 = new std::vector<int>();
1094 
1095 
1096  t_irun = new std::vector<unsigned int>();
1097  t_ievt = new std::vector<unsigned int>();
1098  t_ilum = new std::vector<unsigned int>();
1099 
1100  BuildTree();
1101 }
1102 
1104 
1105  t_v_hnNearTRKs ->clear();
1106  t_v_hnLayers_maxNearP ->clear();
1107  t_v_htrkQual_maxNearP ->clear();
1108  t_v_hmaxNearP_goodTrk ->clear();
1109  t_v_hmaxNearP ->clear();
1110 
1111  t_v_cone_hnNearTRKs ->clear();
1112  t_v_cone_hnLayers_maxNearP ->clear();
1113  t_v_cone_htrkQual_maxNearP ->clear();
1114  t_v_cone_hmaxNearP_goodTrk ->clear();
1115  t_v_cone_hmaxNearP ->clear();
1116 
1117  // t_hScale ->clear();
1118  t_trkNOuterHits ->clear();
1119  t_trkNLayersCrossed ->clear();
1120  t_dtFromLeadJet ->clear();
1121  t_trkP ->clear();
1122  t_trkPt ->clear();
1123  t_trkEta ->clear();
1124  t_trkPhi ->clear();
1125  t_e3x3 ->clear();
1126  t_v_eDR ->clear();
1127  t_v_eMipDR ->clear();
1128  t_h3x3 ->clear();
1129  t_h5x5 ->clear();
1130  t_nRH_h3x3 ->clear();
1131  t_nRH_h5x5 ->clear();
1132 
1133  if (doMC) {
1134  t_simP ->clear();
1135  t_hsim3x3 ->clear();
1136  t_hsim5x5 ->clear();
1137  t_hsim3x3Matched ->clear();
1138  t_hsim5x5Matched ->clear();
1139  t_hsim3x3Rest ->clear();
1140  t_hsim5x5Rest ->clear();
1141  t_hsim3x3Photon ->clear();
1142  t_hsim5x5Photon ->clear();
1143  t_hsim3x3NeutHad ->clear();
1144  t_hsim5x5NeutHad ->clear();
1145  t_hsim3x3CharHad ->clear();
1146  t_hsim5x5CharHad ->clear();
1147  t_hsim3x3PdgMatched ->clear();
1148  t_hsim5x5PdgMatched ->clear();
1149  t_hsim3x3Total ->clear();
1150  t_hsim5x5Total ->clear();
1151 
1152  t_hsim3x3NMatched ->clear();
1153  t_hsim3x3NTotal ->clear();
1154  t_hsim3x3NNeutHad ->clear();
1155  t_hsim3x3NCharHad ->clear();
1156  t_hsim3x3NPhoton ->clear();
1157  t_hsim3x3NRest ->clear();
1158 
1159  t_hsim5x5NMatched ->clear();
1160  t_hsim5x5NTotal ->clear();
1161  t_hsim5x5NNeutHad ->clear();
1162  t_hsim5x5NCharHad ->clear();
1163  t_hsim5x5NPhoton ->clear();
1164  t_hsim5x5NRest ->clear();
1165 
1166  t_trkHcalEne ->clear();
1167  t_trkEcalEne ->clear();
1168  }
1169 
1170  t_distFromHotCell_h3x3 ->clear();
1171  t_ietaFromHotCell_h3x3 ->clear();
1172  t_iphiFromHotCell_h3x3 ->clear();
1173  t_distFromHotCell_h5x5 ->clear();
1174  t_ietaFromHotCell_h5x5 ->clear();
1175  t_iphiFromHotCell_h5x5 ->clear();
1176 
1177  if (doMC) {
1178  t_v_hsimInfoConeMatched ->clear();
1179  t_v_hsimInfoConeRest ->clear();
1180  t_v_hsimInfoConePhoton ->clear();
1181  t_v_hsimInfoConeNeutHad ->clear();
1182  t_v_hsimInfoConeCharHad ->clear();
1183  t_v_hsimInfoConePdgMatched->clear();
1184  t_v_hsimInfoConeTotal ->clear();
1185 
1186  t_v_hsimInfoConeNMatched ->clear();
1187  t_v_hsimInfoConeNRest ->clear();
1188  t_v_hsimInfoConeNPhoton ->clear();
1189  t_v_hsimInfoConeNNeutHad ->clear();
1190  t_v_hsimInfoConeNCharHad ->clear();
1191  t_v_hsimInfoConeNTotal ->clear();
1192 
1193  t_v_hsimCone ->clear();
1194  }
1195 
1196  t_v_hCone ->clear();
1197  t_v_nRecHitsCone->clear();
1198  t_v_nSimHitsCone->clear();
1199 
1200  t_v_distFromHotCell->clear();
1201  t_v_ietaFromHotCell->clear();
1202  t_v_iphiFromHotCell->clear();
1203 
1204  t_v_RH_h3x3_ieta ->clear();
1205  t_v_RH_h3x3_iphi ->clear();
1206  t_v_RH_h3x3_ene ->clear();
1207  t_v_RH_h5x5_ieta ->clear();
1208  t_v_RH_h5x5_iphi ->clear();
1209  t_v_RH_h5x5_ene ->clear();
1210  t_v_RH_r26_ieta ->clear();
1211  t_v_RH_r26_iphi ->clear();
1212  t_v_RH_r26_ene ->clear();
1213  t_v_RH_r44_ieta ->clear();
1214  t_v_RH_r44_iphi ->clear();
1215  t_v_RH_r44_ene ->clear();
1216 
1217  t_v_hlTriggers ->clear();
1218  t_hltHB->clear();
1219  t_hltHE->clear();
1220  t_hltL1Jet15 ->clear();
1221  t_hltJet30 ->clear();
1222  t_hltJet50 ->clear();
1223  t_hltJet80 ->clear();
1224  t_hltJet110 ->clear();
1225  t_hltJet140 ->clear();
1226  t_hltJet180 ->clear();
1227  t_hltL1SingleEG5 ->clear();
1228  t_hltZeroBias ->clear();
1229  t_hltMinBiasHcal ->clear();
1230  t_hltMinBiasEcal ->clear();
1231  t_hltMinBiasPixel ->clear();
1232  t_hltSingleIsoTau30_Trk5 ->clear();
1234 
1235  t_irun->clear();
1236  t_ievt->clear();
1237  t_ilum->clear();
1238 
1239 
1240 
1241 }
1242 
1243 // ----- method called once each job just after ending the event loop ----
1245 
1246  std::cout << "Number of Events Processed " << nEVT << " failed L1 "
1247  << nEVT_failL1 << std::endl;
1248 
1249 }
1250 
1251 
1253 
1254 
1255  hRawPt = fs->make<TH1F>("hRawPt ", "hRawPt ", 100, 0.0, 100.0);
1256  hRawP = fs->make<TH1F>("hRawP ", "hRawP ", 100, 0.0, 100.0);
1257  hRawEta = fs->make<TH1F>("hRawEta", "hRawEta", 15, 0.0, 3.0);
1258  hRawPhi = fs->make<TH1F>("hRawPhi", "hRawPhi", 100, -3.2, 3.2);
1259 
1260  ntp = fs->make<TTree>("ntp", "ntp");
1261 
1262 
1263  // Counters
1264  ntp->Branch("nEVT" , &nEVT , "nEVT/I" );
1265  ntp->Branch("leadL1JetPT" , &leadL1JetPT , "leadL1JetPT/D" );
1266  ntp->Branch("leadL1JetEta", &leadL1JetEta, "leadL1JetEta/D");
1267  ntp->Branch("leadL1JetPhi", &leadL1JetPhi, "leadL1JetPhi/D");
1268  ntp->Branch("nTRK", &nTRK, "nTRK/I");
1269  ntp->Branch("nRawTRK" , &nRawTRK ,"nRawTRK/I" );
1270  ntp->Branch("nFailHighPurityQaul", &nFailHighPurityQaul,"nFailHighPurityQaul/I");
1271  ntp->Branch("nFailPt" , &nFailPt ,"nFailPt/I" );
1272  ntp->Branch("nFailEta" , &nFailEta ,"nFailEta/I" );
1273  ntp->Branch("nMissEcal" , &nMissEcal ,"nMissEcal/I" );
1274  ntp->Branch("nMissHcal" , &nMissHcal ,"nMissHcal/I" );
1275 
1276  ntp->Branch("hnNearTRKs" ,"vector<vector<int> > ",&t_v_hnNearTRKs );
1277  ntp->Branch("hnLayers_maxNearP" ,"vector<vector<int> > ",&t_v_hnLayers_maxNearP );
1278  ntp->Branch("htrkQual_maxNearP" ,"vector<vector<int> > ",&t_v_htrkQual_maxNearP );
1279  ntp->Branch("hmaxNearP_goodTrk" ,"vector<vector<double> >",&t_v_hmaxNearP_goodTrk );
1280  ntp->Branch("hmaxNearP" ,"vector<vector<double> >",&t_v_hmaxNearP );
1281 
1282  ntp->Branch("cone_hnNearTRKs" ,"vector<vector<int> > ",&t_v_cone_hnNearTRKs );
1283  ntp->Branch("cone_hnLayers_maxNearP","vector<vector<int> > ",&t_v_cone_hnLayers_maxNearP);
1284  ntp->Branch("cone_htrkQual_maxNearP","vector<vector<int> > ",&t_v_cone_htrkQual_maxNearP);
1285  ntp->Branch("cone_hmaxNearP_goodTrk","vector<vector<double> >",&t_v_cone_hmaxNearP_goodTrk);
1286  ntp->Branch("cone_hmaxNearP" ,"vector<vector<double> >",&t_v_cone_hmaxNearP );
1287 
1288  // ntp->Branch("hScale" , "vector<double>", &t_hScale );
1289  ntp->Branch("trkNOuterHits" , "vector<double>", &t_trkNOuterHits );
1290  ntp->Branch("trkNLayersCrossed", "vector<double>", &t_trkNLayersCrossed);
1291  ntp->Branch("drFromLeadJet" , "vector<double>", &t_dtFromLeadJet );
1292  ntp->Branch("trkP" , "vector<double>", &t_trkP );
1293  ntp->Branch("trkPt" , "vector<double>", &t_trkPt );
1294  ntp->Branch("trkEta" , "vector<double>", &t_trkEta );
1295  ntp->Branch("trkPhi" , "vector<double>", &t_trkPhi );
1296  ntp->Branch("e3x3" , "vector<double>", &t_e3x3 );
1297 
1298  ntp->Branch("e3x3" , "vector<double>" , &t_e3x3 );
1299  ntp->Branch("v_eDR" , "vector<vector<double> >", &t_v_eDR);
1300  ntp->Branch("v_eMipDR" , "vector<vector<double> >", &t_v_eMipDR);
1301 
1302  ntp->Branch("h3x3" , "vector<double>", &t_h3x3 );
1303  ntp->Branch("h5x5" , "vector<double>", &t_h5x5 );
1304  ntp->Branch("nRH_h3x3" , "vector<double>", &t_nRH_h3x3 );
1305  ntp->Branch("nRH_h5x5" , "vector<double>", &t_nRH_h5x5 );
1306 
1307  if (doMC) {
1308  ntp->Branch("simP" , "vector<double>", &t_simP );
1309  ntp->Branch("hsim3x3" , "vector<double>", &t_hsim3x3 );
1310  ntp->Branch("hsim5x5" , "vector<double>", &t_hsim5x5 );
1311 
1312  ntp->Branch("hsim3x3Matched" , "vector<double>", &t_hsim3x3Matched );
1313  ntp->Branch("hsim5x5Matched" , "vector<double>", &t_hsim5x5Matched );
1314  ntp->Branch("hsim3x3Rest" , "vector<double>", &t_hsim3x3Rest );
1315  ntp->Branch("hsim5x5Rest" , "vector<double>", &t_hsim5x5Rest );
1316  ntp->Branch("hsim3x3Photon" , "vector<double>", &t_hsim3x3Photon );
1317  ntp->Branch("hsim5x5Photon" , "vector<double>", &t_hsim5x5Photon );
1318  ntp->Branch("hsim3x3NeutHad" , "vector<double>", &t_hsim3x3NeutHad );
1319  ntp->Branch("hsim5x5NeutHad" , "vector<double>", &t_hsim5x5NeutHad );
1320  ntp->Branch("hsim3x3CharHad" , "vector<double>", &t_hsim3x3CharHad );
1321  ntp->Branch("hsim5x5CharHad" , "vector<double>", &t_hsim5x5CharHad );
1322  ntp->Branch("hsim3x3PdgMatched", "vector<double>", &t_hsim3x3PdgMatched);
1323  ntp->Branch("hsim5x5PdgMatched", "vector<double>", &t_hsim5x5PdgMatched);
1324  ntp->Branch("hsim3x3Total" , "vector<double>", &t_hsim3x3Total );
1325  ntp->Branch("hsim5x5Total" , "vector<double>", &t_hsim5x5Total );
1326 
1327  ntp->Branch("hsim3x3NMatched" , "vector<int>", &t_hsim3x3NMatched );
1328  ntp->Branch("hsim3x3NRest" , "vector<int>", &t_hsim3x3NRest );
1329  ntp->Branch("hsim3x3NPhoton" , "vector<int>", &t_hsim3x3NPhoton );
1330  ntp->Branch("hsim3x3NNeutHad" , "vector<int>", &t_hsim3x3NNeutHad );
1331  ntp->Branch("hsim3x3NCharHad" , "vector<int>", &t_hsim3x3NCharHad );
1332  ntp->Branch("hsim3x3NTotal" , "vector<int>", &t_hsim3x3NTotal );
1333 
1334  ntp->Branch("hsim5x5NMatched" , "vector<int>", &t_hsim5x5NMatched );
1335  ntp->Branch("hsim5x5NRest" , "vector<int>", &t_hsim5x5NRest );
1336  ntp->Branch("hsim5x5NPhoton" , "vector<int>", &t_hsim5x5NPhoton );
1337  ntp->Branch("hsim5x5NNeutHad" , "vector<int>", &t_hsim5x5NNeutHad );
1338  ntp->Branch("hsim5x5NCharHad" , "vector<int>", &t_hsim5x5NCharHad );
1339  ntp->Branch("hsim5x5NTotal" , "vector<int>", &t_hsim5x5NTotal );
1340 
1341  ntp->Branch("trkHcalEne" , "vector<double>", &t_trkHcalEne );
1342  ntp->Branch("trkEcalEne" , "vector<double>", &t_trkEcalEne );
1343  }
1344 
1345  ntp->Branch("distFromHotCell_h3x3", "vector<double>", &t_distFromHotCell_h3x3);
1346  ntp->Branch("ietaFromHotCell_h3x3", "vector<int>", &t_ietaFromHotCell_h3x3);
1347  ntp->Branch("iphiFromHotCell_h3x3", "vector<int>", &t_iphiFromHotCell_h3x3);
1348  ntp->Branch("distFromHotCell_h5x5", "vector<double>", &t_distFromHotCell_h5x5);
1349  ntp->Branch("ietaFromHotCell_h5x5", "vector<int>", &t_ietaFromHotCell_h5x5);
1350  ntp->Branch("iphiFromHotCell_h5x5", "vector<int>", &t_iphiFromHotCell_h5x5);
1351 
1352  if (doMC) {
1353  ntp->Branch("v_hsimInfoConeMatched" ,"vector<vector<double> >",&t_v_hsimInfoConeMatched );
1354  ntp->Branch("v_hsimInfoConeRest" ,"vector<vector<double> >",&t_v_hsimInfoConeRest );
1355  ntp->Branch("v_hsimInfoConePhoton" ,"vector<vector<double> >",&t_v_hsimInfoConePhoton );
1356  ntp->Branch("v_hsimInfoConeNeutHad" ,"vector<vector<double> >",&t_v_hsimInfoConeNeutHad );
1357  ntp->Branch("v_hsimInfoConeCharHad" ,"vector<vector<double> >",&t_v_hsimInfoConeCharHad );
1358  ntp->Branch("v_hsimInfoConePdgMatched","vector<vector<double> >",&t_v_hsimInfoConePdgMatched);
1359  ntp->Branch("v_hsimInfoConeTotal" ,"vector<vector<double> >",&t_v_hsimInfoConeTotal );
1360 
1361  ntp->Branch("v_hsimInfoConeNMatched" ,"vector<vector<int> >" ,&t_v_hsimInfoConeNMatched );
1362  ntp->Branch("v_hsimInfoConeNRest" ,"vector<vector<int> >" ,&t_v_hsimInfoConeNRest );
1363  ntp->Branch("v_hsimInfoConeNPhoton" ,"vector<vector<int> >" ,&t_v_hsimInfoConeNPhoton );
1364  ntp->Branch("v_hsimInfoConeNNeutHad" ,"vector<vector<int> >" ,&t_v_hsimInfoConeNNeutHad );
1365  ntp->Branch("v_hsimInfoConeNCharHad" ,"vector<vector<int> >" ,&t_v_hsimInfoConeNCharHad );
1366  ntp->Branch("v_hsimInfoConeNTotal" ,"vector<vector<int> >" ,&t_v_hsimInfoConeNTotal );
1367 
1368  ntp->Branch("v_hsimCone" ,"vector<vector<double> >",&t_v_hsimCone );
1369  }
1370 
1371  ntp->Branch("v_hCone" ,"vector<vector<double> >",&t_v_hCone );
1372  ntp->Branch("v_nRecHitsCone" ,"vector<vector<int> >" ,&t_v_nRecHitsCone);
1373  ntp->Branch("v_nSimHitsCone" ,"vector<vector<int> >" ,&t_v_nSimHitsCone);
1374 
1375  ntp->Branch("v_distFromHotCell" ,"vector<vector<double> >",&t_v_distFromHotCell );
1376  ntp->Branch("v_ietaFromHotCell" ,"vector<vector<int> >",&t_v_ietaFromHotCell );
1377  ntp->Branch("v_iphiFromHotCell" ,"vector<vector<int> >",&t_v_iphiFromHotCell );
1378 
1379  ntp->Branch("v_RH_h3x3_ieta" ,"vector<vector<int> >" ,&t_v_RH_h3x3_ieta);
1380  ntp->Branch("v_RH_h3x3_iphi" ,"vector<vector<int> >" ,&t_v_RH_h3x3_iphi);
1381  ntp->Branch("v_RH_h3x3_ene" ,"vector<vector<double> >",&t_v_RH_h3x3_ene );
1382  ntp->Branch("v_RH_h5x5_ieta" ,"vector<vector<int> >" ,&t_v_RH_h5x5_ieta);
1383  ntp->Branch("v_RH_h5x5_iphi" ,"vector<vector<int> >" ,&t_v_RH_h5x5_iphi);
1384  ntp->Branch("v_RH_h5x5_ene" ,"vector<vector<double> >",&t_v_RH_h5x5_ene );
1385  ntp->Branch("v_RH_r26_ieta" ,"vector<vector<int> >" ,&t_v_RH_r26_ieta );
1386  ntp->Branch("v_RH_r26_iphi" ,"vector<vector<int> >" ,&t_v_RH_r26_iphi );
1387  ntp->Branch("v_RH_r26_ene" ,"vector<vector<double> >",&t_v_RH_r26_ene );
1388  ntp->Branch("v_RH_r44_ieta" ,"vector<vector<int> >" ,&t_v_RH_r44_ieta );
1389  ntp->Branch("v_RH_r44_iphi" ,"vector<vector<int> >" ,&t_v_RH_r44_iphi );
1390  ntp->Branch("v_RH_r44_ene" ,"vector<vector<double> >",&t_v_RH_r44_ene );
1391 
1392  ntp->Branch("v_hlTriggers" ,"vector<vector<int> >",&t_v_hlTriggers);
1393  ntp->Branch("v_hltHB" ,"vector<int>",&t_hltHB);
1394  ntp->Branch("v_hltHE" ,"vector<int>",&t_hltHE);
1395  ntp->Branch("v_hltL1Jet15" ,"vector<int>",&t_hltL1Jet15 );
1396  ntp->Branch("v_hltJet30" ,"vector<int>",&t_hltJet30 );
1397  ntp->Branch("v_hltJet50" ,"vector<int>",&t_hltJet50 );
1398  ntp->Branch("v_hltJet80" ,"vector<int>",&t_hltJet80 );
1399  ntp->Branch("v_hltJet110" ,"vector<int>",&t_hltJet110 );
1400  ntp->Branch("v_hltJet140" ,"vector<int>",&t_hltJet140 );
1401  ntp->Branch("v_hltJet180" ,"vector<int>",&t_hltJet180 );
1402  ntp->Branch("v_hltL1SingleEG5" ,"vector<int>",&t_hltL1SingleEG5 );
1403  ntp->Branch("v_hltZeroBias" ,"vector<int>",&t_hltZeroBias );
1404  ntp->Branch("v_hltMinBiasHcal" ,"vector<int>",&t_hltMinBiasHcal );
1405  ntp->Branch("v_hltMinBiasEcal" ,"vector<int>",&t_hltMinBiasEcal );
1406  ntp->Branch("v_hltMinBiasPixel" ,"vector<int>",&t_hltMinBiasPixel );
1407  ntp->Branch("v_hltSingleIsoTau30_Trk5" ,"vector<int>",&t_hltSingleIsoTau30_Trk5 );
1408  ntp->Branch("v_hltDoubleLooseIsoTau15_Trk5" ,"vector<int>",&t_hltDoubleLooseIsoTau15_Trk5);
1409 
1410  ntp->Branch("irun" ,"vector<unsigned int>", &t_irun);
1411  ntp->Branch("ievt" ,"vector<unsigned int>", &t_ievt);
1412  ntp->Branch("ilum" ,"vector<unsigned int>", &t_ilum);
1413 
1414 }
1415 
1416 
1418 
1419  std::string theTrackQuality = "highPurity";
1420  reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality);
1421 
1422  std::cout << " Reference Point " << pTrack->referencePoint() <<"\n"
1423  << " TrackMmentum " << pTrack->momentum()
1424  << " (pt,eta,phi)(" << pTrack->pt()<<","<<pTrack->eta()<<","<<pTrack->phi()<<")"
1425  << " p " << pTrack->p() << "\n"
1426  << " Normalized chi2 " << pTrack->normalizedChi2() <<" charge " << pTrack->charge()
1427  << " qoverp() " << pTrack->qoverp() <<"+-" << pTrack->qoverpError()
1428  << " d0 " << pTrack->d0() << "\n"
1429  << " NValidHits " << pTrack->numberOfValidHits() << " NLostHits " << pTrack->numberOfLostHits()
1430  << " TrackQuality " << pTrack->qualityName(trackQuality_) << " " << pTrack->quality(trackQuality_)
1431  << std::endl;
1432 
1433  if( printTrkHitPattern_ ) {
1434  const reco::HitPattern& p = pTrack->hitPattern();
1435 
1436  for (int i=0; i<p.numberOfHits(); i++) {
1438  }
1439  }
1440 
1441 }
1442 
1443 double IsolatedTracksCone::DeltaPhi(double v1, double v2) {
1444  // Computes the correctly normalized phi difference
1445  // v1, v2 = phi of object 1 and 2
1446 
1447  double pi = 3.141592654;
1448  double twopi = 6.283185307;
1449 
1450  double diff = std::abs(v2 - v1);
1451  double corr = twopi - diff;
1452  if (diff < pi){ return diff;} else { return corr;}
1453 }
1454 
1455 
1456 double IsolatedTracksCone::DeltaR(double eta1, double phi1,
1457  double eta2, double phi2) {
1458  double deta = eta1 - eta2;
1459  double dphi = DeltaPhi(phi1, phi2);
1460  return std::sqrt(deta*deta + dphi*dphi);
1461 }
1462 
1463 
1464 //define this as a plug-in
1466 
1467 
1468 
1469 
1470 
1471 
1472 
1473 
1474 
1475 
1476 
1477 
1478 
1479 
1480 
RunNumber_t run() const
Definition: EventID.h:42
std::vector< int > * t_hsim5x5NTotal
edm::EDGetTokenT< l1extra::L1JetParticleCollection > tok_L1extFwdJet_
std::vector< std::vector< double > > * t_v_distFromHotCell
double qoverp() const
q/p
Definition: TrackBase.h:113
std::vector< double > * t_hsim5x5CharHad
double p() const
momentum vector magnitude
Definition: TrackBase.h:127
T getParameter(std::string const &) const
const CaloSubdetectorGeometry * getSubdetectorGeometry(const DetId &id) const
access the subdetector geometry for the given subdetector directly
Definition: CaloGeometry.cc:43
virtual void beginJob()
Definition: EDAnalyzer.h:67
EventNumber_t event() const
Definition: EventID.h:44
std::vector< std::vector< double > > * t_v_hmaxNearP_goodTrk
T getUntrackedParameter(std::string const &, T const &) const
std::vector< std::vector< int > > * t_v_hlTriggers
int i
Definition: DBlmapReader.cc:9
dictionary parameters
Definition: Parameters.py:2
std::vector< double > * t_trkEcalEne
virtual edm::TriggerNames const & triggerNames(edm::TriggerResults const &triggerResults) const
Definition: Event.cc:204
std::vector< std::vector< double > > * t_v_hsimCone
TrackAssociatorParameters parameters_
static std::string qualityName(TrackQuality)
Definition: TrackBase.h:404
std::vector< int > * t_hltJet30
edm::EDGetTokenT< edm::SimTrackContainer > tok_simTk_
std::vector< double > * t_trkPt
std::vector< int > * t_ietaFromHotCell_h5x5
void eCaloSimInfo(std::vector< DetId > vdets, const CaloGeometry *geo, edm::Handle< T > &hitsEB, edm::Handle< T > &hitsEE, edm::Handle< edm::SimTrackContainer > &SimTk, edm::Handle< edm::SimVertexContainer > &SimVtx, edm::SimTrackContainer::const_iterator trkInfo, caloSimInfo &info, double timeCut=150, bool debug=false)
std::vector< double > * t_nRH_h3x3
std::vector< double > * t_distFromHotCell_h3x3
const Vector & momentum() const
track momentum vector
Definition: TrackBase.h:148
edm::EDGetTokenT< HBHERecHitCollection > tok_hbhe_
std::vector< std::vector< int > > * t_v_hsimInfoConeNCharHad
std::vector< int > * t_hsim5x5NNeutHad
edm::EDGetTokenT< EcalRecHitCollection > tok_EE_
std::vector< double > * t_hsim5x5Photon
double d0() const
dxy parameter in perigee convention (d0 = - dxy)
Definition: TrackBase.h:121
double getDistInPlaneTrackDir(const GlobalPoint &caloPoint, const GlobalVector &caloVector, const GlobalPoint &rechitPoint)
Definition: FindDistCone.cc:8
std::vector< int > * t_hltMinBiasHcal
std::vector< std::vector< double > > * t_v_cone_hmaxNearP_goodTrk
std::vector< double > * t_hsim3x3CharHad
std::vector< int > * t_hsim5x5NPhoton
std::vector< std::vector< double > > * t_v_RH_r44_ene
std::map< std::string, double > eHCALSimInfo(const edm::Event &, const HcalTopology *topology, const DetId &det, const CaloGeometry *geo, edm::Handle< T > &hits, edm::Handle< edm::SimTrackContainer > &SimTk, edm::Handle< edm::SimVertexContainer > &SimVtx, const reco::Track *pTrack, TrackerHitAssociator &associate, int ieta, int iphi, double timeCut=150, bool includeHO=false, bool debug=false)
std::vector< int > * t_iphiFromHotCell_h3x3
double normalizedChi2() const
chi-squared divided by n.d.o.f. (or chi-squared * 1e6 if n.d.o.f. is zero)
Definition: TrackBase.h:109
std::vector< std::vector< int > > * t_v_cone_hnLayers_maxNearP
std::vector< double > * t_hsim5x5PdgMatched
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:434
std::vector< double > * t_hsim3x3NeutHad
std::vector< double > * t_nRH_h5x5
TrackQuality
track quality
Definition: TrackBase.h:93
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
static FreeTrajectoryState getFreeTrajectoryState(const edm::EventSetup &, const reco::Track &)
get FreeTrajectoryState from different track representations
std::vector< std::vector< int > > * t_v_hnLayers_maxNearP
std::vector< int > * t_hltHB
edm::EDGetTokenT< edm::TriggerResults > tok_trigger_
std::vector< std::vector< double > > * t_v_hmaxNearP
std::vector< std::vector< double > > * t_v_RH_h3x3_ene
void useDefaultPropagator()
use the default propagator
edm::EDGetTokenT< reco::TrackCollection > tok_genTrack_
std::vector< int > * t_hsim5x5NCharHad
void loadParameters(const edm::ParameterSet &, edm::ConsumesCollector &)
std::vector< std::vector< int > > * t_v_hsimInfoConeNNeutHad
std::vector< std::vector< int > > * t_v_nRecHitsCone
double phi() const
azimuthal angle of momentum vector
Definition: TrackBase.h:137
std::vector< std::vector< int > > * t_v_hsimInfoConeNMatched
unsigned short numberOfLostHits() const
number of cases where track crossed a layer without getting a hit.
Definition: TrackBase.h:232
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
std::vector< double > * t_trkP
std::vector< int > * t_hltL1SingleEG5
std::vector< int > * t_hsim5x5NRest
std::vector< double > * t_hsim5x5NeutHad
std::vector< int > * t_hltJet50
std::vector< std::vector< double > > * t_v_cone_hmaxNearP
std::vector< double > * t_simP
const Point & referencePoint() const
Reference point on the track.
Definition: TrackBase.h:151
std::vector< int > * t_hsim3x3NCharHad
std::vector< std::vector< int > > * t_v_iphiFromHotCell
std::vector< std::vector< int > > * t_v_hnNearTRKs
std::vector< std::vector< int > > * t_v_cone_htrkQual_maxNearP
math::XYZPoint trkGlobPosAtHcal
std::vector< std::vector< double > > * t_v_RH_r26_ene
const Double_t pi
LuminosityBlockNumber_t luminosityBlock() const
uint32_t rawId() const
get the raw id
Definition: DetId.h:43
std::vector< int > * t_hltDoubleLooseIsoTau15_Trk5
double eHCALmatrix(const HcalTopology *topology, const DetId &det, edm::Handle< T > &hits, int ieta, int iphi, bool includeHO=false, bool algoNew=true, double hbThr=-100, double heThr=-100, double hfThr=-100, double hoThr=-100, double tMin=-500, double tMax=500, bool debug=false)
edm::EDGetTokenT< l1extra::L1JetParticleCollection > tok_L1extTauJet_
std::vector< int > * t_hltJet140
double eCone_hcal(const CaloGeometry *geo, edm::Handle< T > &hits, const GlobalPoint &hpoint1, const GlobalPoint &point1, double dR, const GlobalVector &trackMom, int &nRecHits, double hbThr=-100, double heThr=-100, double hfThr=-100, double hoThr=-100, double tMin=-500, double tMax=500, int detOnly=-1)
int iEvent
Definition: GenABIO.cc:230
double eta() const
pseudorapidity of momentum vector
Definition: TrackBase.h:139
std::vector< int > * t_hltMinBiasPixel
std::vector< double > * t_trkPhi
std::vector< double > * t_hsim3x3Total
edm::EDGetTokenT< edm::PCaloHitContainer > tok_caloEB_
std::vector< int > * t_hltHE
int trackerLayersWithMeasurement() const
Definition: HitPattern.h:740
ConsumesCollector consumesCollector()
Use a ConsumesCollector to gather consumes information from helper functions.
T sqrt(T t)
Definition: SSEVec.h:48
std::vector< std::vector< int > > * t_v_RH_h3x3_iphi
std::vector< int > * t_iphiFromHotCell_h5x5
std::vector< std::vector< double > > * t_v_hsimInfoConePdgMatched
double pt() const
track transverse momentum
Definition: TrackBase.h:129
std::vector< std::vector< int > > * t_v_htrkQual_maxNearP
std::vector< int > * t_hltMinBiasEcal
std::vector< std::vector< int > > * t_v_RH_r44_iphi
std::vector< double > * t_hsim5x5Matched
std::vector< std::vector< double > > * t_v_hsimInfoConeTotal
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
int numberOfHits() const
Definition: HitPattern.cc:211
edm::EDGetTokenT< EcalRecHitCollection > tok_EB_
std::vector< std::vector< double > > * t_v_RH_h5x5_ene
double DeltaPhi(double v1, double v2)
std::vector< std::vector< int > > * t_v_hsimInfoConeNRest
LuminosityBlock const & getLuminosityBlock() const
Definition: Event.h:80
edm::SimTrackContainer::const_iterator matchedSimTrack(const edm::Event &iEvent, edm::Handle< edm::SimTrackContainer > &SimTk, edm::Handle< edm::SimVertexContainer > &SimVtx, const reco::Track *pTrack, TrackerHitAssociator &associate, bool debug=false)
unsigned short numberOfValidHits() const
number of valid hits found
Definition: TrackBase.h:230
std::vector< std::vector< double > > * t_v_hsimInfoConeCharHad
static std::string const triggerResults
Definition: EdmProvDump.cc:41
const HitPattern & hitPattern() const
Access the hit pattern, indicating in which Tracker layers the track has hits.
Definition: TrackBase.h:221
bool first
Definition: L1TdeRCT.cc:75
std::vector< double > * t_hsim5x5Total
bool isValid() const
Definition: HandleBase.h:76
std::vector< double > * t_hsim3x3Photon
std::vector< double > * t_hsim3x3PdgMatched
std::vector< std::vector< int > > * t_v_ietaFromHotCell
std::vector< double > * t_hsim5x5Rest
std::vector< int > * t_hltJet110
JetCorrectorParameters corr
Definition: classes.h:5
std::vector< std::vector< double > > * t_v_eMipDR
std::vector< int > * t_hltL1Jet15
double qoverpError() const
error on signed transverse curvature
Definition: TrackBase.h:190
virtual DetId getClosestCell(const GlobalPoint &r) const
std::vector< int > * t_hsim3x3NRest
std::vector< std::vector< double > > * t_v_hsimInfoConeRest
std::vector< std::vector< int > > * t_v_hsimInfoConeNPhoton
std::vector< std::vector< int > > * t_v_RH_r26_ieta
std::vector< double > * t_dtFromLeadJet
std::vector< std::vector< int > > * t_v_RH_h5x5_ieta
std::vector< std::vector< int > > * t_v_RH_h3x3_ieta
Definition: DetId.h:18
std::vector< unsigned int > * t_ievt
std::vector< double > * t_trkNLayersCrossed
edm::EDGetTokenT< l1extra::L1JetParticleCollection > tok_L1extCenJet_
static TrackQuality qualityByName(const std::string &name)
Definition: TrackBase.cc:46
std::vector< std::vector< double > > * t_v_hsimInfoConeNeutHad
std::vector< std::vector< int > > * t_v_nSimHitsCone
std::vector< int > * t_hltSingleIsoTau30_Trk5
std::vector< double > * t_hsim5x5
GlobalVector trkMomAtEcal
std::string const & triggerName(unsigned int index) const
Definition: TriggerNames.cc:27
const T & get() const
Definition: EventSetup.h:55
std::vector< std::vector< int > > * t_v_RH_h5x5_iphi
T const * product() const
Definition: ESHandle.h:62
IsolatedTracksCone(const edm::ParameterSet &)
std::vector< std::vector< int > > * t_v_hsimInfoConeNTotal
std::vector< double > * t_hsim3x3Matched
T const * product() const
Definition: Handle.h:81
bool quality(const TrackQuality) const
Track quality.
Definition: TrackBase.h:375
std::vector< int > * t_hsim3x3NNeutHad
virtual void analyze(const edm::Event &, const edm::EventSetup &)
edm::Service< TFileService > fs
GlobalVector trkMomAtHcal
edm::EDGetTokenT< edm::SimVertexContainer > tok_simVtx_
std::vector< int > * t_hltZeroBias
edm::EventID id() const
Definition: EventBase.h:56
std::vector< double > * t_h5x5
std::vector< double > * t_trkNOuterHits
std::vector< std::vector< double > > * t_v_hCone
double p1[4]
Definition: TauolaWrapper.h:89
double eCone_ecal(const CaloGeometry *geo, edm::Handle< T > &barrelhits, edm::Handle< T > &endcaphits, const GlobalPoint &hpoint1, const GlobalPoint &point1, double dR, const GlobalVector &trackMom, int &nRecHits, double ebThr=-100, double eeThr=-100, double tMin=-500, double tMax=500)
std::vector< int > * t_hsim3x3NPhoton
std::vector< double > * t_hsim3x3Rest
std::vector< double > * t_hsim3x3
std::vector< int > * t_hsim3x3NMatched
std::vector< std::vector< int > > * t_v_RH_r44_ieta
std::vector< double > * t_trkHcalEne
std::vector< unsigned int > * t_irun
std::vector< int > * t_hltJet180
std::vector< int > * t_hsim3x3NTotal
TrackDetMatchInfo associate(const edm::Event &, const edm::EventSetup &, const FreeTrajectoryState &, const AssociatorParameters &)
void printTrack(const reco::Track *pTrack)
std::vector< int > * t_ietaFromHotCell_h3x3
math::XYZPoint trkGlobPosAtEcal
Track position at different parts of the calorimeter.
tuple cout
Definition: gather_cfg.py:121
int charge() const
track electric charge
Definition: TrackBase.h:111
std::vector< double > * t_distFromHotCell_h5x5
edm::EDGetTokenT< edm::PCaloHitContainer > tok_caloHH_
std::vector< double > * t_e3x3
void printHitPattern(int position, std::ostream &stream) const
Definition: HitPattern.cc:563
double DeltaR(double eta1, double phi1, double eta2, double phi2)
std::vector< double > * t_trkEta
std::vector< double > * t_h3x3
double coneChargeIsolation(const edm::Event &iEvent, const edm::EventSetup &iSetup, reco::TrackCollection::const_iterator trkItr, edm::Handle< reco::TrackCollection > trkCollection, TrackDetectorAssociator &associator, TrackAssociatorParameters &parameters_, std::string theTrackQuality, int &nNearTRKs, int &nLayers_maxNearP, int &trkQual_maxNearP, double &maxNearP_goodTrk, const GlobalPoint &hpoint1, const GlobalVector &trackMom, double dR)
std::vector< std::vector< double > > * t_v_hsimInfoConePhoton
std::vector< unsigned int > * t_ilum
edm::EDGetTokenT< edm::PCaloHitContainer > tok_caloEE_
TrackDetectorAssociator * trackAssociator_
std::vector< int > * t_hltJet80
std::vector< std::vector< double > > * t_v_eDR
std::vector< std::vector< double > > * t_v_hsimInfoConeMatched
std::vector< std::vector< int > > * t_v_RH_r26_iphi
std::vector< int > * t_hsim5x5NMatched
std::vector< std::vector< int > > * t_v_cone_hnNearTRKs
double eECALmatrix(const DetId &detId, edm::Handle< T > &hitsEB, edm::Handle< T > &hitsEE, const CaloGeometry *geo, const CaloTopology *caloTopology, int ieta, int iphi, double ebThr=-100, double eeThr=-100, double tMin=-500, double tMax=500, bool debug=false)