CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HcalIsoTrkAnalyzer.cc
Go to the documentation of this file.
1 // system include files
2 #include <atomic>
3 #include <memory>
4 #include <string>
5 #include <vector>
6 
7 // Root objects
8 #include "TROOT.h"
9 #include "TSystem.h"
10 #include "TFile.h"
11 #include "TProfile.h"
12 #include "TDirectory.h"
13 #include "TTree.h"
14 #include "TLorentzVector.h"
15 #include "TInterpreter.h"
16 
17 //Tracks
22 // Vertices
26 
27 //Triggers
35 
38 
46 
52 
65 
66 //#define DebugLog
67 
68 class HcalIsoTrkAnalyzer : public edm::one::EDAnalyzer<edm::one::WatchRuns,edm::one::SharedResources> {
69 
70 public:
71  explicit HcalIsoTrkAnalyzer(edm::ParameterSet const&);
73 
74  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
75 
76 private:
77  virtual void analyze(edm::Event const&, edm::EventSetup const&) override;
78  virtual void beginJob() ;
79  virtual void beginRun(edm::Run const&, edm::EventSetup const&) override;
80  virtual void endRun(edm::Run const&, edm::EventSetup const&) override;
81  virtual void beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&);
82  virtual void endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&);
83 
84  int fillTree(std::vector< math::XYZTLorentzVector>& vecL1,
85  std::vector< math::XYZTLorentzVector>& vecL3,
86  math::XYZPoint& leadPV,
87  std::vector<spr::propagatedTrackDirection>& trkCaloDirections,
88  const CaloGeometry* geo,
89  edm::Handle<EcalRecHitCollection>& barrelRecHitsHandle,
90  edm::Handle<EcalRecHitCollection>& endcapRecHitsHandle,
93 
94  unsigned int nRun_;
97  std::vector<std::string> trigNames_;
99  int dataType_;
116 
117  TTree *tree, *tree2;
123  std::vector<unsigned int> *t_DetIds;
124  std::vector<double> *t_HitEnergies, pbin;
125  std::vector<bool> *t_trgbits;
127  std::vector<int> *t_ietaAll, *t_ietaGood;
128 };
129 
131  nRun_(0) {
132 
133  usesResource("TFileService");
134 
135  //now do whatever initialization is needed
136  const double isolationRadius(28.9);
137  trigNames_ = iConfig.getParameter<std::vector<std::string> >("Triggers");
138  theTrackQuality_ = iConfig.getParameter<std::string>("TrackQuality");
139  processName_ = iConfig.getParameter<std::string>("ProcessName");
140  l1Filter_ = iConfig.getParameter<std::string>("L1Filter");
141  l2Filter_ = iConfig.getParameter<std::string>("L2Filter");
142  l3Filter_ = iConfig.getParameter<std::string>("L3Filter");
143  reco::TrackBase::TrackQuality trackQuality_=reco::TrackBase::qualityByName(theTrackQuality_);
144  maxRestrictionPt_ = iConfig.getParameter<double>("MinTrackPt");
145  slopeRestrictionPt_ = iConfig.getParameter<double>("SlopeTrackPt");
146  selectionParameter_.minPt = (slopeRestrictionPt_ > 0) ? (maxRestrictionPt_ - 2.5/slopeRestrictionPt_) : maxRestrictionPt_;
147  selectionParameter_.minQuality = trackQuality_;
148  selectionParameter_.maxDxyPV = iConfig.getParameter<double>("MaxDxyPV");
149  selectionParameter_.maxDzPV = iConfig.getParameter<double>("MaxDzPV");
150  selectionParameter_.maxChi2 = iConfig.getParameter<double>("MaxChi2");
151  selectionParameter_.maxDpOverP = iConfig.getParameter<double>("MaxDpOverP");
152  selectionParameter_.minOuterHit = iConfig.getParameter<int>("MinOuterHit");
153  selectionParameter_.minLayerCrossed = iConfig.getParameter<int>("MinLayerCrossed");
154  selectionParameter_.maxInMiss = iConfig.getParameter<int>("MaxInMiss");
155  selectionParameter_.maxOutMiss = iConfig.getParameter<int>("MaxOutMiss");
156  a_coneR_ = iConfig.getParameter<double>("ConeRadius");
157  a_charIsoR_ = a_coneR_ + isolationRadius;
158  a_mipR_ = iConfig.getParameter<double>("ConeRadiusMIP");
159  pTrackMin_ = iConfig.getParameter<double>("MinimumTrackP");
160  eEcalMax_ = iConfig.getParameter<double>("MaximumEcalEnergy");
161  eIsolation_ = iConfig.getParameter<double>("IsolationEnergy");
162  triggerEvent_ = iConfig.getParameter<edm::InputTag>("TriggerEventLabel");
163  theTriggerResultsLabel_ = iConfig.getParameter<edm::InputTag>("TriggerResultLabel");
164  labelGenTrack_ = iConfig.getParameter<std::string>("TrackLabel");
165  labelRecVtx_ = iConfig.getParameter<std::string>("VertexLabel");
166  labelEB_ = iConfig.getParameter<std::string>("EBRecHitLabel");
167  labelEE_ = iConfig.getParameter<std::string>("EERecHitLabel");
168  labelHBHE_ = iConfig.getParameter<std::string>("HBHERecHitLabel");
169  std::string labelBS = iConfig.getParameter<std::string>("BeamSpotLabel");
170  std::string modnam = iConfig.getUntrackedParameter<std::string>("ModuleName","");
171  std::string prdnam = iConfig.getUntrackedParameter<std::string>("ProducerName","");
172  ignoreTrigger_ = iConfig.getUntrackedParameter<bool>("IgnoreTriggers", false);
173  useRaw_ = iConfig.getUntrackedParameter<bool>("UseRaw", false);
174  hcalScale_ = iConfig.getUntrackedParameter<double>("HcalScale", 1.0);
175  dataType_ = iConfig.getUntrackedParameter<int>("DataType", 0);
176 
177  // define tokens for access
178  tok_trigEvt_ = consumes<trigger::TriggerEvent>(triggerEvent_);
179  tok_trigRes_ = consumes<edm::TriggerResults>(theTriggerResultsLabel_);
180  tok_bs_ = consumes<reco::BeamSpot>(labelBS);
181  tok_genTrack_ = consumes<reco::TrackCollection>(labelGenTrack_);
182  if (modnam == "") {
183  tok_recVtx_ = consumes<reco::VertexCollection>(labelRecVtx_);
184  tok_EB_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit",labelEB_));
185  tok_EE_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit",labelEE_));
186  tok_hbhe_ = consumes<HBHERecHitCollection>(labelHBHE_);
187  edm::LogInfo("HcalIsoTrack") << "Labels used " << triggerEvent_ << " "
188  << theTriggerResultsLabel_ << " "
189  << labelBS << " " << labelRecVtx_ << " "
190  << labelGenTrack_ << " "
191  << edm::InputTag("ecalRecHit",labelEB_) << " "
192  << edm::InputTag("ecalRecHit",labelEE_) << " "
193  << labelHBHE_;
194  } else {
195  tok_recVtx_ = consumes<reco::VertexCollection>(edm::InputTag(modnam,labelRecVtx_,prdnam));
196  tok_EB_ = consumes<EcalRecHitCollection>(edm::InputTag(modnam,labelEB_,prdnam));
197  tok_EE_ = consumes<EcalRecHitCollection>(edm::InputTag(modnam,labelEE_,prdnam));
198  tok_hbhe_ = consumes<HBHERecHitCollection>(edm::InputTag(modnam,labelHBHE_,prdnam));
199  edm::LogInfo("HcalIsoTrack") << "Labels used " << triggerEvent_
200  << "\n " << theTriggerResultsLabel_
201  << "\n " << labelBS
202  << "\n " << edm::InputTag(modnam,labelRecVtx_,prdnam)
203  << "\n " << labelGenTrack_
204  << "\n " << edm::InputTag(modnam,labelEB_,prdnam)
205  << "\n " << edm::InputTag(modnam,labelEE_,prdnam)
206  << "\n " << edm::InputTag(modnam,labelHBHE_,prdnam);
207  }
208 
209  edm::LogInfo("HcalIsoTrack") <<"Parameters read from config file \n"
210  <<"\t minPt " << maxRestrictionPt_ <<":" << slopeRestrictionPt_
211  <<"\t theTrackQuality " << theTrackQuality_
212  <<"\t minQuality " << selectionParameter_.minQuality
213  <<"\t maxDxyPV " << selectionParameter_.maxDxyPV
214  <<"\t maxDzPV " << selectionParameter_.maxDzPV
215  <<"\t maxChi2 " << selectionParameter_.maxChi2
216  <<"\t maxDpOverP " << selectionParameter_.maxDpOverP
217  <<"\t minOuterHit " << selectionParameter_.minOuterHit
218  <<"\t minLayerCrossed " << selectionParameter_.minLayerCrossed
219  <<"\t maxInMiss " << selectionParameter_.maxInMiss
220  <<"\t maxOutMiss " << selectionParameter_.maxOutMiss
221  <<"\t a_coneR " << a_coneR_
222  <<"\t a_charIsoR " << a_charIsoR_
223  <<"\t a_mipR " << a_mipR_
224  <<"\t hcalScale_ " << hcalScale_
225  <<"\t useRaw_ " << useRaw_
226  <<"\t ignoreTrigger_ " << ignoreTrigger_
227  <<"\t dataType_ " << dataType_;
228  edm::LogInfo("HcalIsoTrack") << "Process " << processName_ << " L1Filter:"
229  << l1Filter_ << " L2Filter:" << l2Filter_
230  << " L3Filter:" << l3Filter_;
231  for (unsigned int k=0; k<trigNames_.size(); ++k) {
232  edm::LogInfo("HcalIsoTrack") << "Trigger[" << k << "] " << trigNames_[k];
233  }
234 }
235 
237 
239 
240  t_Run = iEvent.id().run();
241  t_Event = iEvent.id().event();
243 #ifdef DebugLog
244  edm::LogInfo("HcalIsoTrack") << "Run " << t_Run << " Event " << t_Event
245  << " type " << t_DataType << " Luminosity "
246  << iEvent.luminosityBlock() << " Bunch "
247  << iEvent.bunchCrossing();
248 #endif
249  //Get magnetic field and ECAL channel status
251  iSetup.get<IdealMagneticFieldRecord>().get(bFieldH);
252  const MagneticField *bField = bFieldH.product();
253 
254  // get handles to calogeometry and calotopology
256  iSetup.get<CaloGeometryRecord>().get(pG);
257  const CaloGeometry* geo = pG.product();
258 
259  bool okC(true);
260  //Get track collection
262  iEvent.getByToken(tok_genTrack_, trkCollection);
263  reco::TrackCollection::const_iterator trkItr;
264  if (!trkCollection.isValid()) {
265  edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelGenTrack_;
266  okC = false;
267  }
268 
269  //event weight for FLAT sample
270  t_EventWeight = 1.0;
271 
272  //Define the best vertex and the beamspot
274  iEvent.getByToken(tok_recVtx_, recVtxs);
275  edm::Handle<reco::BeamSpot> beamSpotH;
276  iEvent.getByToken(tok_bs_, beamSpotH);
277  math::XYZPoint leadPV(0,0,0);
278  t_goodPV = 0;
279  if (recVtxs.isValid() && recVtxs->size()>0) {
280  for (unsigned int k=0; k<recVtxs->size(); ++k) {
281  if (!((*recVtxs)[k].isFake()) && ((*recVtxs)[k].ndof() > 4)) {
282  if (t_goodPV == 0) leadPV = math::XYZPoint((*recVtxs)[k].x(),(*recVtxs)[k].y(),(*recVtxs)[k].z());
283  t_goodPV++;
284  }
285  }
286  }
287  if (t_goodPV == 0 && beamSpotH.isValid()) {
288  leadPV = beamSpotH->position();
289  }
290 #ifdef DebugLog
291  edm::LogInfo("HcalIsoTrack") << "Primary Vertex " << leadPV << " out of "
292  << t_goodPV << " vertex";
293  if (beamSpotH.isValid()) {
294  edm::LogInfo("HcalIsoTrack") << " Beam Spot " << beamSpotH->position();
295  }
296 #endif
297  // RecHits
298  edm::Handle<EcalRecHitCollection> barrelRecHitsHandle;
299  iEvent.getByToken(tok_EB_, barrelRecHitsHandle);
300  if (!barrelRecHitsHandle.isValid()) {
301  edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelEB_;
302  okC = false;
303  }
304  edm::Handle<EcalRecHitCollection> endcapRecHitsHandle;
305  iEvent.getByToken(tok_EE_, endcapRecHitsHandle);
306  if (!endcapRecHitsHandle.isValid()) {
307  edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelEE_;
308  okC = false;
309  }
311  iEvent.getByToken(tok_hbhe_, hbhe);
312  if (!hbhe.isValid()) {
313  edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelHBHE_;
314  okC = false;
315  }
316 
317  //Propagate tracks to calorimeter surface)
318  std::vector<spr::propagatedTrackDirection> trkCaloDirections;
319  spr::propagateCALO(trkCollection, geo, bField, theTrackQuality_,
320  trkCaloDirections, false);
321  std::vector<math::XYZTLorentzVector> vecL1, vecL3;
322  t_Tracks = trkCollection->size();
323  t_TracksProp = trkCaloDirections.size();
324  t_ietaAll->clear(); t_ietaGood->clear();
325 #ifdef DebugLog
326  edm::LogInfo("HcalIsoTrack") << "# of propagated tracks " << t_TracksProp
327  << " out of " << t_Tracks << " with Trigger "
328  << ignoreTrigger_;
329 #endif
330  //Trigger
331  t_trgbits->clear();
332  t_TracksSaved = 0;
333  for (unsigned int i=0; i<trigNames_.size(); ++i) t_trgbits->push_back(false);
334  if (ignoreTrigger_) {
335  t_l1pt = t_l1eta = t_l1phi = 0;
336  t_l3pt = t_l3eta = t_l3phi = 0;
337  t_TracksSaved = fillTree(vecL1, vecL3, leadPV, trkCaloDirections, geo,
338  barrelRecHitsHandle, endcapRecHitsHandle, hbhe);
339  } else {
340  trigger::TriggerEvent triggerEvent;
341  edm::Handle<trigger::TriggerEvent> triggerEventHandle;
342  iEvent.getByToken(tok_trigEvt_, triggerEventHandle);
343  if (!triggerEventHandle.isValid()) {
344  edm::LogWarning("HcalIsoTrack") << "Error! Can't get the product "
345  << triggerEvent_.label() ;
346  } else if (okC) {
347  triggerEvent = *(triggerEventHandle.product());
348  const trigger::TriggerObjectCollection& TOC(triggerEvent.getObjects());
351  iEvent.getByToken(tok_trigRes_, triggerResults);
352  bool done(false);
353  if (triggerResults.isValid()) {
354  std::vector<std::string> modules;
355  const edm::TriggerNames & triggerNames = iEvent.triggerNames(*triggerResults);
356  const std::vector<std::string> & triggerNames_ = triggerNames.triggerNames();
357  int ntksave(0);
358  for (unsigned int iHLT=0; iHLT<triggerResults->size(); iHLT++) {
359  bool ok(false);
360  int hlt = triggerResults->accept(iHLT);
361  if (trigNames_.size() > 0) {
362  for (unsigned int i=0; i<trigNames_.size(); ++i) {
363  if (triggerNames_[iHLT].find(trigNames_[i].c_str())!=std::string::npos) {
364  t_trgbits->at(i) = (hlt>0);
365  if (hlt > 0) ok = true;
366 #ifdef DebugLog
367  edm::LogInfo("HcalIsoTrack") << "This trigger "
368  << triggerNames_[iHLT] << " Flag "
369  << hlt << ":" << ok;
370 #endif
371  }
372  }
373  } else {
374  ok = (hlt>0);
375  }
376 
377  if (ok) {
378  unsigned int triggerindx = hltConfig_.triggerIndex(triggerNames_[iHLT]);
379  const std::vector<std::string>& moduleLabels(hltConfig_.moduleLabels(triggerindx));
380  std::vector<math::XYZTLorentzVector> vecL2;
381  vecL1.clear(); vecL3.clear();
382  //loop over all trigger filters in event (i.e. filters passed)
383  for (unsigned int ifilter=0; ifilter<triggerEvent.sizeFilters();
384  ++ifilter) {
385  std::vector<int> Keys;
386  std::string label = triggerEvent.filterTag(ifilter).label();
387  //loop over keys to objects passing this filter
388  for (unsigned int imodule=0; imodule<moduleLabels.size(); imodule++) {
389  if (label.find(moduleLabels[imodule]) != std::string::npos) {
390 #ifdef DebugLog
391  edm::LogInfo("HcalIsoTrack") << "FilterName " << label;
392 #endif
393  for (unsigned int ifiltrKey=0; ifiltrKey<triggerEvent.filterKeys(ifilter).size(); ++ifiltrKey) {
394  Keys.push_back(triggerEvent.filterKeys(ifilter)[ifiltrKey]);
395  const trigger::TriggerObject& TO(TOC[Keys[ifiltrKey]]);
396  math::XYZTLorentzVector v4(TO.px(), TO.py(), TO.pz(), TO.energy());
397  if (label.find(l2Filter_) != std::string::npos) {
398  vecL2.push_back(v4);
399  } else if (label.find(l3Filter_) != std::string::npos) {
400  vecL3.push_back(v4);
401  } else if (label.find(l1Filter_) != std::string::npos || l1Filter_ == "") {
402  vecL1.push_back(v4);
403  }
404 #ifdef DebugLog
405  edm::LogInfo("HcalIsoTrack") << "key " <<ifiltrKey<<" : pt "
406  << TO.pt() << " eta "<<TO.eta()
407  << " phi " <<TO.phi()<<" mass "
408  << TO.mass() <<" Id "<<TO.id();
409 #endif
410  }
411 #ifdef DebugLog
412  edm::LogInfo("HcalIsoTrack") << "sizes " << vecL1.size() <<":"
413  << vecL2.size() <<":" <<vecL3.size();
414 #endif
415  }
416  }
417  }
419  math::XYZTLorentzVector mindRvec1;
420  double mindR1(999);
421  for (unsigned int i=0; i<vecL2.size(); i++) {
422  double dr = dR(vecL1[0],vecL2[i]);
423 #ifdef DebugLog
424  edm::LogInfo("HcalIsoTrack") << "lvl2[" << i << "] dR " << dr;
425 #endif
426  if (dr<mindR1) {
427  mindR1 = dr;
428  mindRvec1 = vecL2[i];
429  }
430  }
431 #ifdef DebugLog
432  edm::LogInfo("HcalIsoTrack") << "L2 object closest to L1 "
433  << mindRvec1 << " at Dr " << mindR1;
434 #endif
435 
436  if (vecL1.size()>0) {
437  t_l1pt = vecL1[0].pt();
438  t_l1eta = vecL1[0].eta();
439  t_l1phi = vecL1[0].phi();
440  } else {
441  t_l1pt = t_l1eta = t_l1phi = 0;
442  }
443  if (vecL3.size()>0) {
444  t_l3pt = vecL3[0].pt();
445  t_l3eta = vecL3[0].eta();
446  t_l3phi = vecL3[0].phi();
447  } else {
448  t_l3pt = t_l3eta = t_l3phi = 0;
449  }
450  // Now fill in the tree for each selected track
451  if (!done) {
452  ntksave = fillTree(vecL1, vecL3, leadPV, trkCaloDirections, geo,
453  barrelRecHitsHandle,endcapRecHitsHandle, hbhe);
454  done = true;
455  }
456  t_TracksSaved += ntksave;
457  }
458  }
459  }
460  }
461  }
462  tree2->Fill();
463 }
464 
466 
467  tree = fs->make<TTree>("CalibTree", "CalibTree");
468 
469  tree->Branch("t_Run", &t_Run, "t_Run/I");
470  tree->Branch("t_Event", &t_Event, "t_Event/I");
471  tree->Branch("t_DataType", &t_DataType, "t_DataType/I");
472  tree->Branch("t_ieta", &t_ieta, "t_ieta/I");
473  tree->Branch("t_EventWeight", &t_EventWeight, "t_EventWeight/D");
474  tree->Branch("t_goodPV", &t_goodPV, "t_goodPV/I");
475  tree->Branch("t_l1pt", &t_l1pt, "t_l1pt/D");
476  tree->Branch("t_l1eta", &t_l1eta, "t_l1eta/D");
477  tree->Branch("t_l1phi", &t_l1phi, "t_l1phi/D");
478  tree->Branch("t_l3pt", &t_l3pt, "t_l3pt/D");
479  tree->Branch("t_l3eta", &t_l3eta, "t_l3eta/D");
480  tree->Branch("t_l3phi", &t_l3phi, "t_l3phi/D");
481  tree->Branch("t_p", &t_p, "t_p/D");
482  tree->Branch("t_mindR1", &t_mindR1, "t_mindR1/D");
483  tree->Branch("t_mindR2", &t_mindR2, "t_mindR2/D");
484  tree->Branch("t_eMipDR", &t_eMipDR, "t_eMipDR/D");
485  tree->Branch("t_eHcal", &t_eHcal, "t_eHcal/D");
486  tree->Branch("t_hmaxNearP", &t_hmaxNearP, "t_hmaxNearP/D");
487  tree->Branch("t_selectTk", &t_selectTk, "t_selectTk/O");
488  tree->Branch("t_qltyFlag", &t_qltyFlag, "t_qltyFlag/O");
489  tree->Branch("t_qltyMissFlag",&t_qltyMissFlag,"t_qltyMissFlag/O");
490  tree->Branch("t_qltyPVFlag", &t_qltyPVFlag, "t_qltyPVFlag/O)");
491 
492  t_DetIds = new std::vector<unsigned int>();
493  t_HitEnergies = new std::vector<double>();
494  t_trgbits = new std::vector<bool>();
495  tree->Branch("t_DetIds", "std::vector<unsigned int>", &t_DetIds);
496  tree->Branch("t_HitEnergies", "std::vector<double>", &t_HitEnergies);
497  tree->Branch("t_trgbits", "std::vector<bool>", &t_trgbits);
498 
499  tree2 = fs->make<TTree>("EventInfo", "Event Information");
500 
501  tree2->Branch("t_Tracks", &t_Tracks, "t_Tracks/I");
502  tree2->Branch("t_TracksProp", &t_TracksProp, "t_TracksProp/I");
503  tree2->Branch("t_TracksSaved", &t_TracksSaved, "t_TracksSaved/I");
504  t_ietaAll = new std::vector<int>();
505  t_ietaGood = new std::vector<int>();
506  tree2->Branch("t_ietaAll", "std::vector<int>", &t_ietaAll);
507  tree2->Branch("t_ietaGood", "std::vector<int>", &t_ietaGood);
508 }
509 
510 
511 // ------------ method called when starting to processes a run ------------
512 void HcalIsoTrkAnalyzer::beginRun(edm::Run const& iRun, edm::EventSetup const& iSetup) {
513  bool changed_(true);
514  bool flag = hltConfig_.init(iRun,iSetup,processName_,changed_);
515  edm::LogInfo("HcalIsoTrack") << "Run[" << nRun_ << "] " << iRun.run()
516  << " process " << processName_ << " init flag "
517  << flag << " change flag " << changed_;
518 
519  // check if trigger names in (new) config
520  if (changed_) {
521  changed_ = false;
522 #ifdef DebugLog
523  edm::LogInfo("HcalIsoTrack") << "New trigger menu found !!!";
524 #endif
525  const unsigned int n(hltConfig_.size());
526  for (unsigned itrig=0; itrig<trigNames_.size(); itrig++) {
527  unsigned int triggerindx = hltConfig_.triggerIndex(trigNames_[itrig]);
528  if (triggerindx >= n) {
529  edm::LogWarning("HcalIsoTrack") << trigNames_[itrig] << " "
530  << triggerindx << " does not exist in "
531  << "the current menu";
532 #ifdef DebugLog
533  } else {
534  edm::LogInfo("HcalIsoTrack") << trigNames_[itrig] << " "
535  << triggerindx << " exists";
536 #endif
537  }
538  }
539  }
540 }
541 
542 // ------------ method called when ending the processing of a run ------------
544  nRun_++;
545  edm::LogInfo("HcalIsoTrack") << "endRun[" << nRun_ << "] " << iRun.run();
546 }
547 
548 // ------------ method called when starting to processes a luminosity block ------------
550 // ------------ method called when ending the processing of a luminosity block ------------
552 // ------------ method fills 'descriptions' with the allowed parameters for the module ------------
554  //The following says we do not know what parameters are allowed so do no validation
555  // Please change this to state exactly what you do use, even if it is no parameters
557  desc.setUnknown();
558  descriptions.addDefault(desc);
559 }
560 
561 int HcalIsoTrkAnalyzer::fillTree(std::vector< math::XYZTLorentzVector>& vecL1,
562  std::vector< math::XYZTLorentzVector>& vecL3,
563  math::XYZPoint& leadPV,
564  std::vector<spr::propagatedTrackDirection>& trkCaloDirections,
565  const CaloGeometry* geo,
566  edm::Handle<EcalRecHitCollection>& barrelRecHitsHandle,
567  edm::Handle<EcalRecHitCollection>& endcapRecHitsHandle,
569 
570  int nsave = 0;
571  //Loop over tracks
572  std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
573  unsigned int nTracks(0), nselTracks(0);
574  for (trkDetItr = trkCaloDirections.begin(),nTracks=0;
575  trkDetItr != trkCaloDirections.end(); trkDetItr++,nTracks++) {
576  const reco::Track* pTrack = &(*(trkDetItr->trkItr));
577  math::XYZTLorentzVector v4(pTrack->px(), pTrack->py(),
578  pTrack->pz(), pTrack->p());
579 #ifdef DebugLog
580  edm::LogInfo("HcalIsoTrack") << "This track : " << nTracks
581  << " (pt|eta|phi|p) :" << pTrack->pt()
582  << "|" << pTrack->eta() << "|"
583  << pTrack->phi() << "|" <<pTrack->p();
584 #endif
585  t_mindR2 = 999;
586  for (unsigned int k=0; k<vecL3.size(); ++k) {
587  double dr = dR(vecL3[k],v4);
588  if (dr<t_mindR2) {
589  t_mindR2 = dr;
590  }
591  }
592  t_mindR1 = (vecL1.size() > 0) ? dR(vecL1[0],v4) : 999;
593 #ifdef DebugLog
594  edm::LogInfo("HcalIsoTrack") << "Closest L3 object at dr :"
595  << t_mindR2 << " and from L1 " << t_mindR1;
596 #endif
597  //Selection of good track
599  t_selectTk = spr::goodTrack(pTrack,leadPV,selectionParameter_,false);
601  oneCutParameters.maxDxyPV = 10;
602  oneCutParameters.maxDzPV = 100;
603  oneCutParameters.maxInMiss = 2;
604  oneCutParameters.maxOutMiss= 2;
605  bool qltyFlag = spr::goodTrack(pTrack,leadPV,oneCutParameters,false);
606  oneCutParameters = selectionParameter_;
607  oneCutParameters.maxDxyPV = 10;
608  oneCutParameters.maxDzPV = 100;
609  t_qltyMissFlag = spr::goodTrack(pTrack,leadPV,oneCutParameters,false);
610  oneCutParameters = selectionParameter_;
611  oneCutParameters.maxInMiss = 2;
612  oneCutParameters.maxOutMiss= 2;
613  t_qltyPVFlag = spr::goodTrack(pTrack,leadPV,oneCutParameters,false);
614  t_ieta = 0;
615  t_p = pTrack->p();
616  if (trkDetItr->okHCAL) {
617  HcalDetId detId = (HcalDetId)(trkDetItr->detIdHCAL);
618  t_ieta = detId.ieta();
619  if (t_p > 40.0 && t_p <= 60.0) t_ietaAll->push_back(t_ieta);
620  }
621 #ifdef DebugLog
622  edm::LogInfo("HcalIsoTrack") << "qltyFlag|okECAL|okHCAL : "
623  << qltyFlag << "|" << trkDetItr->okECAL
624  << "|" << trkDetItr->okHCAL;
625 #endif
626  t_qltyFlag = (qltyFlag && trkDetItr->okECAL && trkDetItr->okHCAL);
627  if (t_qltyFlag) {
628  nselTracks++;
629  int nRH_eMipDR(0), nNearTRKs(0), nRecHits(-999);
630  t_eMipDR = spr::eCone_ecal(geo, barrelRecHitsHandle,
631  endcapRecHitsHandle, trkDetItr->pointHCAL,
632  trkDetItr->pointECAL, a_mipR_,
633  trkDetItr->directionECAL, nRH_eMipDR);
634  t_DetIds->clear(); t_HitEnergies->clear();
635  std::vector<DetId> ids;
636  t_eHcal = spr::eCone_hcal(geo, hbhe, trkDetItr->pointHCAL,
637  trkDetItr->pointECAL, a_coneR_,
638  trkDetItr->directionHCAL,nRecHits,
639  ids, *t_HitEnergies, useRaw_);
640  t_eHcal *= hcalScale_;
641  for (unsigned int k=0; k<ids.size(); ++k) {
642  t_DetIds->push_back(ids[k].rawId());
643  }
644  t_hmaxNearP = spr::chargeIsolationCone(nTracks, trkCaloDirections,
645  a_charIsoR_, nNearTRKs, false);
646 #ifdef DebugLog
647  edm::LogInfo("HcalIsoTrack") << "This track : " << nTracks
648  << " (pt|eta|phi|p) :" << pTrack->pt()
649  << "|" << pTrack->eta() << "|"
650  << pTrack->phi() << "|" << t_p;
651  edm::LogInfo("HcalIsoTrack") << "e_MIP " << t_eMipDR
652  << " Chg Isolation " << t_hmaxNearP
653  << " eHcal" << t_eHcal << " ieta "
654  << t_ieta << " Quality "
655  << t_qltyMissFlag << ":"
656  << t_qltyPVFlag << ":" << t_selectTk;
657  for (unsigned int lll=0;lll<t_DetIds->size();lll++) {
658  edm::LogInfo("HcalIsoTrack") << "det id is = " <<t_DetIds->at(lll)
659  << " hit enery is = "
660  << t_HitEnergies->at(lll) ;
661  }
662 #endif
663  if (t_p>pTrackMin_ && t_eMipDR<eEcalMax_ &&
665  tree->Fill();
666  nsave++;
667  if (t_p > 40.0 && t_p <= 60.0 && t_selectTk) t_ietaGood->push_back(t_ieta);
668 #ifdef DebugLog
669  for (unsigned int k=0; k<t_trgbits->size(); k++) {
670  edm::LogInfo("HcalIsoTrack") << "trigger bit is = "
671  << t_trgbits->at(k);
672  }
673 #endif
674  }
675  }
676  }
677  return nsave;
678 }
679 
681  return reco::deltaR(vec1.eta(),vec1.phi(),vec2.eta(),vec2.phi());
682 }
683 
684 //define this as a plug-in
RunNumber_t run() const
Definition: EventID.h:39
unsigned int size() const
number of trigger paths in trigger table
double p() const
momentum vector magnitude
Definition: TrackBase.h:610
spr::trackSelectionParameters selectionParameter_
T getParameter(std::string const &) const
EventNumber_t event() const
Definition: EventID.h:41
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
virtual edm::TriggerNames const & triggerNames(edm::TriggerResults const &triggerResults) const
Definition: Event.cc:215
const unsigned int nTracks(const reco::Vertex &sv)
std::vector< spr::propagatedTrackID > propagateCALO(edm::Handle< reco::TrackCollection > &trkCollection, const CaloGeometry *geo, const MagneticField *bField, std::string &theTrackQuality, bool debug=false)
edm::EDGetTokenT< EcalRecHitCollection > tok_EE_
int id() const
getters
Definition: TriggerObject.h:55
RunNumber_t run() const
Definition: RunBase.h:40
The single EDProduct to be saved for each event (AOD case)
Definition: TriggerEvent.h:25
trigger::size_type sizeFilters() const
Definition: TriggerEvent.h:135
std::vector< std::string > trigNames_
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:462
float phi() const
Definition: TriggerObject.h:58
TrackQuality
track quality
Definition: TrackBase.h:151
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
std::vector< int > * t_ietaAll
const Keys & filterKeys(trigger::size_type index) const
Definition: TriggerEvent.h:111
edm::InputTag theTriggerResultsLabel_
int bunchCrossing() const
Definition: EventBase.h:65
edm::LuminosityBlockNumber_t luminosityBlock() const
Definition: EventBase.h:62
double phi() const
azimuthal angle of momentum vector
Definition: TrackBase.h:640
virtual void beginLuminosityBlock(edm::LuminosityBlock const &, edm::EventSetup const &)
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
edm::EDGetTokenT< EcalRecHitCollection > tok_EB_
float energy() const
Definition: TriggerObject.h:65
std::vector< double > * t_HitEnergies
double chargeIsolationCone(unsigned int trkIndex, std::vector< spr::propagatedTrackDirection > &trkDirs, double dR, int &nNearTRKs, bool debug=false)
virtual void beginRun(edm::Run const &, edm::EventSetup const &) override
edm::EDGetTokenT< reco::TrackCollection > tok_genTrack_
float eta() const
Definition: TriggerObject.h:57
double px() const
x coordinate of momentum vector
Definition: TrackBase.h:622
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:7
Strings const & triggerNames() const
Definition: TriggerNames.cc:24
std::vector< unsigned int > * t_DetIds
edm::EDGetTokenT< HBHERecHitCollection > tok_hbhe_
std::vector< double > pbin
virtual void analyze(edm::Event const &, edm::EventSetup const &) override
unsigned int triggerIndex(const std::string &triggerName) const
slot position of trigger path in trigger table (0 to size-1)
Single trigger physics object (e.g., an isolated muon)
Definition: TriggerObject.h:22
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
int iEvent
Definition: GenABIO.cc:230
edm::EDGetTokenT< trigger::TriggerEvent > tok_trigEvt_
double eta() const
pseudorapidity of momentum vector
Definition: TrackBase.h:646
void addDefault(ParameterSetDescription const &psetDescription)
edm::EDGetTokenT< reco::BeamSpot > tok_bs_
edm::Service< TFileService > fs
bool goodTrack(const reco::Track *pTrack, math::XYZPoint leadPV, trackSelectionParameters parameters, bool debug=false)
const TriggerObjectCollection & getObjects() const
Definition: TriggerEvent.h:98
std::vector< double > vec1
Definition: HCALResponse.h:15
double pt() const
track transverse momentum
Definition: TrackBase.h:616
int ieta() const
get the cell ieta
Definition: HcalDetId.h:56
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
auto deltaR(const T1 &t1, const T2 &t2) -> decltype(t1.eta())
Definition: deltaR.h:28
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, bool debug=false)
static std::string const triggerResults
Definition: EdmProvDump.cc:41
virtual void endLuminosityBlock(edm::LuminosityBlock const &, edm::EventSetup const &)
HLTConfigProvider hltConfig_
bool isValid() const
Definition: HandleBase.h:75
double pz() const
z coordinate of momentum vector
Definition: TrackBase.h:634
const std::vector< std::string > & moduleLabels(unsigned int trigger) const
label(s) of module(s) on a trigger path
std::vector< bool > * t_trgbits
std::vector< TriggerObject > TriggerObjectCollection
collection of trigger physics objects (e.g., all isolated muons)
Definition: TriggerObject.h:81
const edm::InputTag filterTag(trigger::size_type index) const
Definition: TriggerEvent.h:103
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
static TrackQuality qualityByName(const std::string &name)
Definition: TrackBase.cc:125
T const * product() const
Definition: Handle.h:81
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
std::vector< size_type > Keys
bool init(const edm::Run &iRun, const edm::EventSetup &iSetup, const std::string &processName, bool &changed)
d&#39;tor
const T & get() const
Definition: EventSetup.h:56
T const * product() const
Definition: ESHandle.h:86
double dR(math::XYZTLorentzVector &, math::XYZTLorentzVector &)
std::string const & label() const
Definition: InputTag.h:36
edm::EventID id() const
Definition: EventBase.h:59
virtual void endRun(edm::Run const &, edm::EventSetup const &) override
HcalIsoTrkAnalyzer(edm::ParameterSet const &)
int fillTree(std::vector< math::XYZTLorentzVector > &vecL1, std::vector< math::XYZTLorentzVector > &vecL3, math::XYZPoint &leadPV, std::vector< spr::propagatedTrackDirection > &trkCaloDirections, const CaloGeometry *geo, edm::Handle< EcalRecHitCollection > &barrelRecHitsHandle, edm::Handle< EcalRecHitCollection > &endcapRecHitsHandle, edm::Handle< HBHERecHitCollection > &hbhe)
reco::TrackBase::TrackQuality minQuality
edm::EDGetTokenT< edm::TriggerResults > tok_trigRes_
float mass() const
Definition: TriggerObject.h:59
std::vector< int > * t_ietaGood
edm::EDGetTokenT< reco::VertexCollection > tok_recVtx_
edm::InputTag triggerEvent_
std::vector< vec1 > vec2
Definition: HCALResponse.h:16
double py() const
y coordinate of momentum vector
Definition: TrackBase.h:628
Definition: Run.h:43
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, bool useRaw=false, bool debug=false)