CMS 3D CMS Logo
StudyHLT.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: IsolatedParticles
4 // Class: StudyHLT
5 //
13 //
14 // Original Author: Sunanda Banerjee
15 // Created: Thu Mar 4 18:52:02 CST 2011
16 //
17 //
18 
19 // system include files
20 #include <memory>
21 #include <string>
22 
23 // Root objects
24 #include "TH1.h"
25 #include "TH2.h"
26 #include "TTree.h"
27 
28 // user include files
29 #include "Calibration/IsolatedParticles/interface/FindCaloHit.h"
30 #include "Calibration/IsolatedParticles/interface/CaloPropagateTrack.h"
31 #include "Calibration/IsolatedParticles/interface/ChargeIsolation.h"
32 #include "Calibration/IsolatedParticles/interface/eECALMatrix.h"
33 #include "Calibration/IsolatedParticles/interface/eHCALMatrix.h"
34 #include "Calibration/IsolatedParticles/interface/TrackSelection.h"
35 
36 #include "CommonTools/UtilAlgos/interface/TFileService.h"
37 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
38 
39 #include "DataFormats/Common/interface/TriggerResults.h"
40 #include "DataFormats/DetId/interface/DetId.h"
41 #include "DataFormats/EcalDetId/interface/EBDetId.h"
42 #include "DataFormats/EcalDetId/interface/EEDetId.h"
43 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
44 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
45 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
46 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
47 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
48 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
49 #include "DataFormats/HLTReco/interface/TriggerObject.h"
50 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
51 #include "DataFormats/Luminosity/interface/LumiDetails.h"
52 #include "DataFormats/TrackReco/interface/Track.h"
53 #include "DataFormats/TrackReco/interface/TrackFwd.h"
54 #include "DataFormats/TrackReco/interface/HitPattern.h"
55 #include "DataFormats/VertexReco/interface/VertexFwd.h"
56 #include "DataFormats/VertexReco/interface/Vertex.h"
57 
58 #include "FWCore/Framework/interface/Frameworkfwd.h"
59 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
60 #include "FWCore/Framework/interface/Event.h"
61 #include "FWCore/Framework/interface/MakerMacros.h"
62 #include "FWCore/ServiceRegistry/interface/Service.h"
63 #include "FWCore/ParameterSet/interface/ParameterSet.h"
64 #include "FWCore/Common/interface/TriggerNames.h"
65 #include "FWCore/MessageLogger/interface/MessageLogger.h"
66 
67 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
68 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
69 #include "Geometry/CaloEventSetup/interface/CaloTopologyRecord.h"
70 #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
71 #include "Geometry/CaloTopology/interface/EcalTrigTowerConstituentsMap.h"
72 #include "Geometry/CaloTopology/interface/HcalTopology.h"
73 #include "Geometry/CaloTopology/interface/CaloTopology.h"
74 
75 #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
76 #include "MagneticField/Engine/interface/MagneticField.h"
77 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
78 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
79 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
80 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
81 
82 class StudyHLT : public edm::one::EDAnalyzer<edm::one::WatchRuns,edm::one::SharedResources> {
83 
84 public:
85  explicit StudyHLT(const edm::ParameterSet&);
86  ~StudyHLT() override { }
87 
88  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
89 
90 private:
91  void analyze(edm::Event const&, edm::EventSetup const&) override;
92  void beginJob() override;
93  void beginRun(edm::Run const&, edm::EventSetup const&) override;
94  void endRun(edm::Run const&, edm::EventSetup const&) override;
95  virtual void beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) { }
96  virtual void endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) { }
97 
98  void clear();
99  void fillTrack(int, double, double, double, double);
100  void fillIsolation(int, double, double, double);
101  void fillEnergy(int, int, double, double, double, double, double);
102  std::string truncate_str(const std::string&);
103  int trackPID(const reco::Track* ,
104  const edm::Handle<reco::GenParticleCollection>&);
105 
106  // ----------member data ---------------------------
107  static const int nPBin_=10, nEtaBin_=4, nPVBin_=4;
108  static const int nGen_=(nPVBin_+12);
109  HLTConfigProvider hltConfig_;
110  edm::Service<TFileService> fs_;
111  const int verbosity_;
112  const std::vector<std::string> trigNames_, newNames_;
113  const std::string theTrackQuality_;
114  const double minTrackP_, maxTrackEta_;
115  const double tMinE_, tMaxE_, tMinH_, tMaxH_;
116  const bool isItAOD_, doTree_;
117  const std::vector<double> puWeights_;
118  const edm::InputTag triggerEvent_, theTriggerResultsLabel_;
119  spr::trackSelectionParameters selectionParameters_;
120  std::vector<std::string> HLTNames_;
121  bool changed_, firstEvent_;
122 
123  edm::EDGetTokenT<LumiDetails> tok_lumi;
124  edm::EDGetTokenT<trigger::TriggerEvent> tok_trigEvt;
125  edm::EDGetTokenT<edm::TriggerResults> tok_trigRes;
126  edm::EDGetTokenT<reco::GenParticleCollection> tok_parts_;
127  edm::EDGetTokenT<reco::TrackCollection> tok_genTrack_;
128  edm::EDGetTokenT<reco::VertexCollection> tok_recVtx_;
129  edm::EDGetTokenT<EcalRecHitCollection> tok_EB_;
130  edm::EDGetTokenT<EcalRecHitCollection> tok_EE_;
131  edm::EDGetTokenT<HBHERecHitCollection> tok_hbhe_;
132  edm::EDGetTokenT<GenEventInfoProduct> tok_ew_;
133 
134  TH1I *h_nHLT, *h_HLTAccept, *h_HLTCorr, *h_numberPV;
135  TH1I *h_goodPV, *h_goodRun;
136  TH2I *h_nHLTvsRN;
137  std::vector<TH1I*> h_HLTAccepts;
138  TH1D *h_p[nGen_+2], *h_pt[nGen_+2];
139  TH1D *h_eta[nGen_+2], *h_phi[nGen_+2];
140  TH1I *h_ntrk[2];
141  TH1D *h_maxNearP[2], *h_ene1[2], *h_ene2[2], *h_ediff[2];
142  TH1D *h_energy[nPVBin_+8][nPBin_][nEtaBin_][6];
143  TTree *tree_;
144  int nRun_, etaBin_[nEtaBin_+1], pvBin_[nPVBin_+1];
145  double pBin_[nPBin_+1];
146  int tr_goodPV, tr_goodRun;
147  double tr_eventWeight;
148  std::vector<std::string> tr_TrigName;
149  std::vector<double> tr_TrkPt, tr_TrkP, tr_TrkEta, tr_TrkPhi;
150  std::vector<double> tr_MaxNearP31X31, tr_MaxNearHcalP7x7;
151  std::vector<double> tr_H3x3, tr_H5x5, tr_H7x7;
152  std::vector<double> tr_FE7x7P, tr_FE11x11P, tr_FE15x15P;
153  std::vector<bool> tr_SE7x7P, tr_SE11x11P, tr_SE15x15P;
154  std::vector<int> tr_iEta, tr_TrkID;
155 };
156 
157 StudyHLT::StudyHLT(const edm::ParameterSet& iConfig) :
158  verbosity_(iConfig.getUntrackedParameter<int>("verbosity",0)),
159  trigNames_(iConfig.getUntrackedParameter<std::vector<std::string> >("triggers")),
160  newNames_(iConfig.getUntrackedParameter<std::vector<std::string> >("newNames")),
161  theTrackQuality_(iConfig.getUntrackedParameter<std::string>("trackQuality","highPurity")),
162  minTrackP_(iConfig.getUntrackedParameter<double>("minTrackP",1.0)),
163  maxTrackEta_(iConfig.getUntrackedParameter<double>("maxTrackEta",2.5)),
164  tMinE_(iConfig.getUntrackedParameter<double>("timeMinCutECAL",-500.)),
165  tMaxE_(iConfig.getUntrackedParameter<double>("timeMaxCutECAL",500.)),
166  tMinH_(iConfig.getUntrackedParameter<double>("timeMinCutHCAL",-500.)),
167  tMaxH_(iConfig.getUntrackedParameter<double>("timeMaxCutHCAL",500.)),
168  isItAOD_(iConfig.getUntrackedParameter<bool>("isItAOD",false)),
169  doTree_(iConfig.getUntrackedParameter<bool>("doTree",false)),
170  puWeights_(iConfig.getUntrackedParameter<std::vector<double> >("puWeights")),
171  triggerEvent_(edm::InputTag("hltTriggerSummaryAOD","","HLT")),
172  theTriggerResultsLabel_(edm::InputTag("TriggerResults","","HLT")), nRun_(0) {
173 
174  usesResource(TFileService::kSharedResource);
175 
176  reco::TrackBase::TrackQuality trackQuality=reco::TrackBase::qualityByName(theTrackQuality_);
177  selectionParameters_.minPt = iConfig.getUntrackedParameter<double>("minTrackPt", 10.0);
178  selectionParameters_.minQuality = trackQuality;
179  selectionParameters_.maxDxyPV = iConfig.getUntrackedParameter<double>("maxDxyPV", 0.2);
180  selectionParameters_.maxDzPV = iConfig.getUntrackedParameter<double>("maxDzPV", 5.0);
181  selectionParameters_.maxChi2 = iConfig.getUntrackedParameter<double>("maxChi2", 5.0);
182  selectionParameters_.maxDpOverP = iConfig.getUntrackedParameter<double>("maxDpOverP", 0.1);
183  selectionParameters_.minOuterHit = iConfig.getUntrackedParameter<int>("minOuterHit", 4);
184  selectionParameters_.minLayerCrossed= iConfig.getUntrackedParameter<int>("minLayerCrossed", 8);
185  selectionParameters_.maxInMiss = iConfig.getUntrackedParameter<int>("maxInMiss", 0);
186  selectionParameters_.maxOutMiss = iConfig.getUntrackedParameter<int>("maxOutMiss", 0);
187 
188  // define tokens for access
189  tok_lumi = consumes<LumiDetails, edm::InLumi>(edm::InputTag("lumiProducer"));
190  tok_trigEvt = consumes<trigger::TriggerEvent>(triggerEvent_);
191  tok_trigRes = consumes<edm::TriggerResults>(theTriggerResultsLabel_);
192  tok_genTrack_ = consumes<reco::TrackCollection>(edm::InputTag("generalTracks"));
193  tok_recVtx_ = consumes<reco::VertexCollection>(edm::InputTag("offlinePrimaryVertices"));
194  tok_parts_ = consumes<reco::GenParticleCollection>(iConfig.getParameter<edm::InputTag>("particleSource"));
195 
196  if (isItAOD_) {
197  tok_EB_ = consumes<EcalRecHitCollection>(edm::InputTag("reducedEcalRecHitsEB"));
198  tok_EE_ = consumes<EcalRecHitCollection>(edm::InputTag("reducedEcalRecHitsEE"));
199  tok_hbhe_ = consumes<HBHERecHitCollection>(edm::InputTag("reducedHcalRecHits", "hbhereco"));
200  } else {
201  tok_EB_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit","EcalRecHitsEB"));
202  tok_EE_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit","EcalRecHitsEE"));
203  tok_hbhe_ = consumes<HBHERecHitCollection>(edm::InputTag("hbhereco"));
204  }
205  tok_ew_ = consumes<GenEventInfoProduct>(edm::InputTag("generator"));
206 
207  edm::LogInfo("IsoTrack") << "Verbosity " << verbosity_ << " with "
208  << trigNames_.size() << " triggers:";
209  for (unsigned int k=0; k<trigNames_.size(); ++k)
210  edm::LogInfo("IsoTrack") << " [" << k << "] " << trigNames_[k];
211  edm::LogInfo("IsoTrack") << "TrackQuality " << theTrackQuality_ << " Minpt "
212  << selectionParameters_.minPt << " maxDxy "
213  << selectionParameters_.maxDxyPV << " maxDz "
214  << selectionParameters_.maxDzPV << " maxChi2 "
215  << selectionParameters_.maxChi2 << " maxDp/p "
216  << selectionParameters_.maxDpOverP << " minOuterHit "
217  << selectionParameters_.minOuterHit << " minLayerCrossed "
218  << selectionParameters_.minLayerCrossed << " maxInMiss "
219  << selectionParameters_.maxInMiss << " maxOutMiss "
220  << selectionParameters_.maxOutMiss << " minTrackP "
221  << minTrackP_ << " maxTrackEta " << maxTrackEta_
222  << " tMinE_ " << tMinE_ << " tMaxE " << tMaxE_
223  << " tMinH_ " << tMinH_ << " tMaxH_ " << tMaxH_
224  << " isItAOD " << isItAOD_ << " doTree " << doTree_;
225 
226  double pBins[nPBin_+1] = {1.0,2.0,3.0,4.0,5.0,6.0,7.0,9.0,11.0,15.0,20.0};
227  int etaBins[nEtaBin_+1] = {1, 7, 13, 17, 23};
228  int pvBins[nPVBin_+1] = {1, 2, 3, 5, 100};
229  for (int i=0; i<=nPBin_; ++i) pBin_[i] = pBins[i];
230  for (int i=0; i<=nEtaBin_; ++i) etaBin_[i] = etaBins[i];
231  for (int i=0; i<=nPVBin_; ++i) pvBin_[i] = pvBins[i];
232 
233  firstEvent_ = true;
234  changed_ = false;
235 }
236 
237 void StudyHLT::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
238  std::vector<std::string> trig;
239  std::vector<double> weights;
240  std::vector<std::string> newNames = {"HLT", "PixelTracks_Multiplicity",
241  "HLT_Physics_","HLT_JetE",
242  "HLT_ZeroBias"};
243  edm::ParameterSetDescription desc;
244  desc.add<edm::InputTag>("particleSource",edm::InputTag("genParticles"));
245  desc.addUntracked<int>("verbosity", 0);
246  desc.addUntracked<std::vector<std::string> >("triggers",trig);
247  desc.addUntracked<std::vector<std::string> >("newNames",newNames);
248  desc.addUntracked<std::string>("trackQuality","highPurity");
249  desc.addUntracked<double>("minTrackPt", 1.0);
250  desc.addUntracked<double>("maxDxyPV", 0.02);
251  desc.addUntracked<double>("maxDzPV", 0.02);
252  desc.addUntracked<double>("maxChi2", 5.0);
253  desc.addUntracked<double>("maxDpOverP", 0.1);
254  desc.addUntracked<int>("minOuterHit", 4);
255  desc.addUntracked<int>("minLayerCrossed", 8);
256  desc.addUntracked<int>("maxInMiss", 0);
257  desc.addUntracked<int>("maxOutMiss", 0);
258  desc.addUntracked<double>("minTrackP", 1.0);
259  desc.addUntracked<double>("maxTrackEta", 2.6);
260  desc.addUntracked<double>("timeMinCutECAL",-500.0);
261  desc.addUntracked<double>("timeMaxCutECAL", 500.0);
262  desc.addUntracked<double>("timeMinCutHCAL",-500.0);
263  desc.addUntracked<double>("timeMaxCutHCAL", 500.0);
264  desc.addUntracked<bool>("isItAOD", false);
265  desc.addUntracked<bool>("doTree", false);
266  desc.addUntracked<std::vector<double> >("puWeights", weights);
267  descriptions.add("studyHLT",desc);
268 }
269 
270 void StudyHLT::analyze(edm::Event const& iEvent, edm::EventSetup const& iSetup) {
271  clear();
272  if (verbosity_ > 0)
273  edm::LogInfo("IsoTrack") << "Event starts====================================";
274  int RunNo = iEvent.id().run();
275  int EvtNo = iEvent.id().event();
276  int Lumi = iEvent.luminosityBlock();
277  int Bunch = iEvent.bunchCrossing();
278 
279  std::vector<int> newAccept(newNames_.size()+1,0);
280  float mybxlumi=-1;
281  /*
282  edm::Handle<LumiDetails> Lumid;
283  iEvent.getLuminosityBlock().getByToken(tok_lumi,Lumid);
284  if (Lumid.isValid()) mybxlumi=Lumid->lumiValue(LumiDetails::kOCC1,iEvent.bunchCrossing())*6.37;
285  */
286  if (verbosity_ > 0)
287  edm::LogInfo("IsoTrack") << "RunNo " << RunNo << " EvtNo " << EvtNo
288  << " Lumi " << Lumi << " Bunch " << Bunch
289  << " mybxlumi " << mybxlumi;
290 
291  trigger::TriggerEvent triggerEvent;
292  edm::Handle<trigger::TriggerEvent> triggerEventHandle;
293  iEvent.getByToken(tok_trigEvt,triggerEventHandle);
294 
295  bool ok(false);
296  if (!triggerEventHandle.isValid()) {
297  edm::LogWarning("IsoTrack") << "Error! Can't get the product "
298  << triggerEvent_.label();
299  } else {
300  triggerEvent = *(triggerEventHandle.product());
301 
303  edm::Handle<edm::TriggerResults> triggerResults;
304  iEvent.getByToken(tok_trigRes, triggerResults);
305 
306  if (triggerResults.isValid()) {
307  h_nHLT->Fill(triggerResults->size());
308  h_nHLTvsRN->Fill(RunNo, triggerResults->size());
309 
310  const edm::TriggerNames & triggerNames = iEvent.triggerNames(*triggerResults);
311  const std::vector<std::string> & triggerNames_ = triggerNames.triggerNames();
312  for (unsigned int iHLT=0; iHLT<triggerResults->size(); iHLT++) {
313  // unsigned int triggerindx = hltConfig_.triggerIndex(triggerNames_[iHLT]);
314  // const std::vector<std::string>& moduleLabels(hltConfig_.moduleLabels(triggerindx));
315  int ipos=-1;
316  std::string newtriggerName = truncate_str(triggerNames_[iHLT]);
317  for (unsigned int i=0; i<HLTNames_.size(); ++i) {
318  if (newtriggerName == HLTNames_[i]) {
319  ipos = i+1;
320  break;
321  }
322  }
323  if (ipos < 0) {
324  HLTNames_.push_back(newtriggerName);
325  ipos = (int)(HLTNames_.size());
326  if (ipos <= h_HLTAccept->GetNbinsX())
327  h_HLTAccept->GetXaxis()->SetBinLabel(ipos,newtriggerName.c_str());
328  }
329  if ((int)(iHLT+1) > h_HLTAccepts[nRun_]->GetNbinsX()) {
330  edm::LogInfo("IsoTrack") << "Wrong trigger " << RunNo << " Event "
331  << EvtNo << " Hlt " << iHLT;
332  } else {
333  if (firstEvent_) h_HLTAccepts[nRun_]->GetXaxis()->SetBinLabel(iHLT+1, newtriggerName.c_str());
334  }
335  int hlt = triggerResults->accept(iHLT);
336  if (hlt) {
337  h_HLTAccepts[nRun_]->Fill(iHLT+1);
338  h_HLTAccept->Fill(ipos);
339  }
340  if (trigNames_.empty()) {
341  ok = true;
342  } else {
343  for (unsigned int i=0; i<trigNames_.size(); ++i) {
344  if (newtriggerName.find(trigNames_[i]) != std::string::npos) {
345  if (verbosity_%10 > 0)
346  edm::LogInfo("IsoTrack") << newtriggerName;
347  if (hlt > 0) {
348  ok = true;
349  tr_TrigName.push_back(newtriggerName);
350  }
351  }
352  }
353  for (unsigned int i=0; i<newNames_.size(); ++i) {
354  if (newtriggerName.find(newNames_[i]) != std::string::npos) {
355  if (verbosity_%10 > 0)
356  edm::LogInfo("IsoTrack") << "[" << i << "] " << newNames_[i]
357  << " : " << newtriggerName;
358  if (hlt > 0) newAccept[i] = 1;
359  }
360  }
361  }
362  }
363  int iflg(0), indx(1);
364  for (unsigned int i=0; i<newNames_.size(); ++i) {
365  iflg += (indx*newAccept[i]); indx *= 2;
366  }
367  h_HLTCorr->Fill(iflg);
368  }
369  }
370 
371  //Look at the tracks
372  if (ok) {
373  h_goodRun->Fill(RunNo);
374  tr_goodRun = RunNo;
375  // get handles to calogeometry and calotopology
376  edm::ESHandle<CaloGeometry> pG;
377  iSetup.get<CaloGeometryRecord>().get(pG);
378  const CaloGeometry* geo = pG.product();
379 
380  edm::ESHandle<CaloTopology> theCaloTopology;
381  iSetup.get<CaloTopologyRecord>().get(theCaloTopology);
382  const CaloTopology *caloTopology = theCaloTopology.product();
383 
384  edm::ESHandle<HcalTopology> htopo;
385  iSetup.get<HcalRecNumberingRecord>().get(htopo);
386  const HcalTopology* theHBHETopology = htopo.product();
387 
388  edm::ESHandle<MagneticField> bFieldH;
389  iSetup.get<IdealMagneticFieldRecord>().get(bFieldH);
390  const MagneticField *bField = bFieldH.product();
391 
392  edm::ESHandle<EcalChannelStatus> ecalChStatus;
393  iSetup.get<EcalChannelStatusRcd>().get(ecalChStatus);
394  const EcalChannelStatus* theEcalChStatus = ecalChStatus.product();
395 
396  edm::Handle<reco::VertexCollection> recVtxs;
397  iEvent.getByToken(tok_recVtx_,recVtxs);
398  int ntrk(0), ngoodPV(0), nPV(-1);
399  int nvtxs = (int)(recVtxs->size());
400  for (int ind=0; ind<nvtxs; ind++) {
401  if (!((*recVtxs)[ind].isFake()) && (*recVtxs)[ind].ndof() > 4) ngoodPV++;
402  }
403  for (int i=0; i<nPVBin_; ++i) {
404  if (ngoodPV >= pvBin_[i] && ngoodPV < pvBin_[i+1]) {
405  nPV = i; break;
406  }
407  }
408 
409  tr_eventWeight = 1.0;
410  edm::Handle<GenEventInfoProduct> genEventInfo;
411  iEvent.getByToken(tok_ew_, genEventInfo);
412  if (genEventInfo.isValid()) tr_eventWeight = genEventInfo->weight();
413 
414  if ((verbosity_/10)%10 > 0)
415  edm::LogInfo("IsoTrack") << "Number of vertices: " << nvtxs
416  << " Good " << ngoodPV << " Bin " << nPV
417  << " Event weight " << tr_eventWeight;
418  h_numberPV->Fill(nvtxs,tr_eventWeight);
419  h_goodPV->Fill(ngoodPV,tr_eventWeight);
420  tr_goodPV = ngoodPV;
421 
422  if (!puWeights_.empty()) {
423  int npbin = h_goodPV->FindBin(ngoodPV);
424  if (npbin > 0 && npbin <= (int)(puWeights_.size()))
425  tr_eventWeight *= puWeights_[npbin-1];
426  else
427  tr_eventWeight = 0;
428  }
429 
430  edm::Handle<reco::TrackCollection> trkCollection;
431  iEvent.getByToken(tok_genTrack_, trkCollection);
432 
433  //=== genParticle information
434  edm::Handle<reco::GenParticleCollection> genParticles;
435  iEvent.getByToken(tok_parts_, genParticles);
436  if (genParticles.isValid()) {
437  for (const auto & p : (reco::GenParticleCollection)(*genParticles)) {
438  double pt1 = p.momentum().Rho();
439  double p1 = p.momentum().R();
440  double eta1 = p.momentum().Eta();
441  double phi1 = p.momentum().Phi();
442  fillTrack(nGen_, pt1,p1,eta1,phi1);
443  bool match(false);
444  double phi2(phi1);
445  if (phi2 < 0) phi2 += 2.0*M_PI;
446  for (const auto& trk : (reco::TrackCollection)(*trkCollection)) {
447  bool quality = trk.quality(selectionParameters_.minQuality);
448  if (quality) {
449  double dEta = trk.eta() - eta1;
450  double phi0 = trk.phi();
451  if (phi0 < 0) phi0 += 2.0*M_PI;
452  double dPhi = phi0-phi2;
453  if (dPhi > M_PI) dPhi -= 2.*M_PI;
454  else if (dPhi < -M_PI) dPhi += 2.*M_PI;
455  double dR = sqrt(dEta*dEta+dPhi*dPhi);
456  if (dR < 0.01) { match = true; break;}
457  }
458  }
459  if (match) fillTrack(nGen_+1, pt1,p1,eta1,phi1);
460  }
461  }
462 
463  reco::TrackCollection::const_iterator trkItr;
464  for (trkItr=trkCollection->begin(); trkItr != trkCollection->end(); ++trkItr,++ntrk) {
465  const reco::Track* pTrack = &(*trkItr);
466  double pt1 = pTrack->pt();
467  double p1 = pTrack->p();
468  double eta1 = pTrack->momentum().eta();
469  double phi1 = pTrack->momentum().phi();
470  bool quality = pTrack->quality(selectionParameters_.minQuality);
471  fillTrack(0, pt1,p1,eta1,phi1);
472  if (quality) fillTrack(1, pt1,p1,eta1,phi1);
473  }
474  h_ntrk[0]->Fill(ntrk,tr_eventWeight);
475 
476  std::vector<spr::propagatedTrackID> trkCaloDets;
477  spr::propagateCALO(trkCollection, geo, bField, theTrackQuality_, trkCaloDets, ((verbosity_/100)%10 > 0));
478  std::vector<spr::propagatedTrackID>::const_iterator trkDetItr;
479  for (trkDetItr = trkCaloDets.begin(),ntrk=0; trkDetItr != trkCaloDets.end(); trkDetItr++,ntrk++) {
480  const reco::Track* pTrack = &(*(trkDetItr->trkItr));
481  double pt1 = pTrack->pt();
482  double p1 = pTrack->p();
483  double eta1 = pTrack->momentum().eta();
484  double phi1 = pTrack->momentum().phi();
485  if ((verbosity_/10)%10 > 0)
486  edm::LogInfo("IsoTrack") << "track: p " << p1 << " pt " << pt1
487  << " eta " << eta1 << " phi " << phi1
488  << " okEcal " << trkDetItr->okECAL;
489  fillTrack(2, pt1,p1,eta1,phi1);
490  if (pt1>minTrackP_ && std::abs(eta1)<maxTrackEta_ && trkDetItr->okECAL) {
491  fillTrack(3, pt1,p1,eta1,phi1);
492  double maxNearP31x31 = spr::chargeIsolationEcal(ntrk, trkCaloDets, geo, caloTopology, 15, 15, ((verbosity_/1000)%10 > 0));
493 
494  edm::ESHandle<EcalSeverityLevelAlgo> sevlv;
495  iSetup.get<EcalSeverityLevelAlgoRcd>().get(sevlv);
496 
497  edm::Handle<EcalRecHitCollection> barrelRecHitsHandle;
498  edm::Handle<EcalRecHitCollection> endcapRecHitsHandle;
499  iEvent.getByToken(tok_EB_,barrelRecHitsHandle);
500  iEvent.getByToken(tok_EE_,endcapRecHitsHandle);
501  // get ECal Tranverse Profile
502  std::pair<double, bool> e7x7P, e11x11P, e15x15P;
503  const DetId isoCell = trkDetItr->detIdECAL;
504  e7x7P = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),3,3, 0.030, 0.150, tMinE_,tMaxE_, ((verbosity_/10000)%10 > 0));
505  e11x11P = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),5,5, 0.030, 0.150, tMinE_,tMaxE_, ((verbosity_/10000)%10 > 0));
506  e15x15P = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),7,7, 0.030, 0.150, tMinE_,tMaxE_, ((verbosity_/10000)%10 > 0));
507 
508  double maxNearHcalP7x7 = spr::chargeIsolationHcal(ntrk, trkCaloDets, theHBHETopology, 3,3, ((verbosity_/1000)%10 > 0));
509  int ieta(0);
510  double h3x3(0), h5x5(0), h7x7(0);
511  fillIsolation(0, maxNearP31x31,e11x11P.first,e15x15P.first);
512  if ((verbosity_/10)%10 > 0)
513  edm::LogInfo("IsoTrack") << "Accepted Tracks reaching Ecal maxNearP31x31 "
514  << maxNearP31x31 << " e11x11P "
515  << e11x11P.first << " e15x15P "
516  << e15x15P.first << " okHCAL "
517  << trkDetItr->okHCAL;
518 
519  int trackID = trackPID(pTrack,genParticles);
520  if (trkDetItr->okHCAL) {
521  edm::Handle<HBHERecHitCollection> hbhe;
522  iEvent.getByToken(tok_hbhe_, hbhe);
523  const DetId ClosestCell(trkDetItr->detIdHCAL);
524  ieta = ((HcalDetId)(ClosestCell)).ietaAbs();
525  h3x3 = spr::eHCALmatrix(theHBHETopology, ClosestCell, hbhe,1,1, false, true, 0.7, 0.8, -100.0, -100.0, tMinH_,tMaxH_, ((verbosity_/10000)%10 > 0));
526  h5x5 = spr::eHCALmatrix(theHBHETopology, ClosestCell, hbhe,2,2, false, true, 0.7, 0.8, -100.0, -100.0, tMinH_,tMaxH_, ((verbosity_/10000)%10 > 0) );
527  h7x7 = spr::eHCALmatrix(theHBHETopology, ClosestCell, hbhe,3,3, false, true, 0.7, 0.8, -100.0, -100.0, tMinH_,tMaxH_, ((verbosity_/10000)%10 > 0) );
528  fillIsolation(1, maxNearHcalP7x7,h5x5,h7x7);
529  if ((verbosity_/10)%10 > 0)
530  edm::LogInfo("IsoTrack") << "Tracks Reaching Hcal maxNearHcalP7x7/h5x5/h7x7 "
531  << maxNearHcalP7x7 << "/" << h5x5 << "/" << h7x7;
532  tr_TrkPt.push_back(pt1);
533  tr_TrkP.push_back(p1);
534  tr_TrkEta.push_back(eta1);
535  tr_TrkPhi.push_back(phi1);
536  tr_TrkID.push_back(trackID);
537  tr_MaxNearP31X31.push_back(maxNearP31x31);
538  tr_MaxNearHcalP7x7.push_back(maxNearHcalP7x7);
539  tr_FE7x7P.push_back(e7x7P.first);
540  tr_FE11x11P.push_back(e11x11P.first);
541  tr_FE15x15P.push_back(e15x15P.first);
542  tr_SE7x7P.push_back(e7x7P.second);
543  tr_SE11x11P.push_back(e11x11P.second);
544  tr_SE15x15P.push_back(e15x15P.second);
545  tr_iEta.push_back(ieta);
546  tr_H3x3.push_back(h3x3);
547  tr_H5x5.push_back(h5x5);
548  tr_H7x7.push_back(h7x7);
549 
550  if (maxNearP31x31 < 0) {
551  fillTrack(4, pt1,p1,eta1,phi1);
552  fillEnergy(0,ieta,p1,e7x7P.first,h3x3,e11x11P.first,h5x5);
553  if (maxNearHcalP7x7 < 0) {
554  fillTrack(5, pt1,p1,eta1,phi1);
555  fillEnergy(1,ieta,p1,e7x7P.first,h3x3,e11x11P.first,h5x5);
556  if ((e11x11P.second) && (e15x15P.second) &&
557  (e15x15P.first-e11x11P.first)<2.0) {
558  fillTrack(6, pt1,p1,eta1,phi1);
559  fillEnergy(2,ieta,p1,e7x7P.first,h3x3,e11x11P.first,h5x5);
560  if (h7x7-h5x5 < 2.0) {
561  fillTrack(7, pt1,p1,eta1,phi1);
562  fillEnergy(3,ieta,p1,e7x7P.first,h3x3,e11x11P.first,h5x5);
563  if (nPV >= 0) {
564  fillTrack(nPV+8, pt1,p1,eta1,phi1);
565  fillEnergy(nPV+4,ieta,p1,e7x7P.first,h3x3,e11x11P.first,h5x5);
566  }
567  if (trackID > 0) {
568  fillTrack(nPVBin_+trackID+7, pt1,p1,eta1,phi1);
569  fillEnergy(nPVBin_+trackID+3,ieta,p1,e7x7P.first,h3x3,e11x11P.first,h5x5);
570  }
571  }
572  }
573  }
574  }
575  }
576  }
577  }
578  h_ntrk[1]->Fill(ntrk,tr_eventWeight);
579  if ((!tr_TrkPt.empty()) && doTree_) tree_->Fill();
580  }
581  firstEvent_ = false;
582 }
583 
584 void StudyHLT::beginJob() {
585  // Book histograms
586  h_nHLT = fs_->make<TH1I>("h_nHLT" , "size of trigger Names", 1000, 0, 1000);
587  h_HLTAccept = fs_->make<TH1I>("h_HLTAccept", "HLT Accepts for all runs", 500, 0, 500);
588  for (int i=1; i<=500; ++i) h_HLTAccept->GetXaxis()->SetBinLabel(i," ");
589  h_nHLTvsRN = fs_->make<TH2I>("h_nHLTvsRN" , "size of trigger Names vs RunNo", 2168, 190949, 193116, 100, 400, 500);
590  h_HLTCorr = fs_->make<TH1I>("h_HLTCorr", "Correlation among different paths", 100, 0, 100);
591  h_numberPV = fs_->make<TH1I>("h_numberPV", "Number of Primary Vertex", 100, 0, 100);
592  h_goodPV = fs_->make<TH1I>("h_goodPV", "Number of good Primary Vertex", 100, 0, 100);
593  h_goodRun = fs_->make<TH1I>("h_goodRun","Number of accepted events for Run", 4000, 190000, 1940000);
594  char hname[60], htit[200];
595  std::string CollectionNames[2] = {"Reco", "Propagated"};
596  for (unsigned int i=0; i<2; i++) {
597  sprintf(hname, "h_nTrk_%s", CollectionNames[i].c_str());
598  sprintf(htit, "Number of %s tracks", CollectionNames[i].c_str());
599  h_ntrk[i] = fs_->make<TH1I>(hname, htit, 500, 0, 500);
600  }
601  std::string TrkNames[8] = {"All", "Quality", "NoIso", "okEcal", "EcalCharIso", "HcalCharIso", "EcalNeutIso", "HcalNeutIso"};
602  std::string particle[4] = {"Electron", "Pion", "Kaon", "Proton"};
603  for (unsigned int i=0; i<=nGen_+1; i++) {
604  if (i < 8) {
605  sprintf(hname, "h_pt_%s", TrkNames[i].c_str());
606  sprintf(htit, "p_{T} of %s tracks", TrkNames[i].c_str());
607  } else if (i < 8+nPVBin_) {
608  sprintf(hname, "h_pt_%s_%d", TrkNames[7].c_str(), i-8);
609  sprintf(htit, "p_{T} of %s tracks (PV=%d:%d)", TrkNames[7].c_str(), pvBin_[i-8], pvBin_[i-7]-1);
610  } else if (i >= nGen_) {
611  sprintf(hname, "h_pt_%s_%d", TrkNames[0].c_str(), i-nGen_);
612  sprintf(htit, "p_{T} of %s Generator tracks", TrkNames[0].c_str());
613  } else {
614  sprintf(hname, "h_pt_%s_%s", TrkNames[7].c_str(), particle[i-8-nPVBin_].c_str());
615  sprintf(htit, "p_{T} of %s tracks (%s)", TrkNames[7].c_str(), particle[i-8-nPVBin_].c_str());
616  }
617  h_pt[i] = fs_->make<TH1D>(hname, htit, 400, 0, 200.0);
618  h_pt[i]->Sumw2();
619 
620  if (i < 8) {
621  sprintf(hname, "h_p_%s", TrkNames[i].c_str());
622  sprintf(htit, "Momentum of %s tracks", TrkNames[i].c_str());
623  } else if (i < 8+nPVBin_) {
624  sprintf(hname, "h_p_%s_%d", TrkNames[7].c_str(), i-8);
625  sprintf(htit, "Momentum of %s tracks (PV=%d:%d)", TrkNames[7].c_str(), pvBin_[i-8], pvBin_[i-7]-1);
626  } else if (i >= nGen_) {
627  sprintf(hname, "h_p_%s_%d", TrkNames[0].c_str(), i-nGen_);
628  sprintf(htit, "Momentum of %s Generator tracks", TrkNames[0].c_str());
629  } else {
630  sprintf(hname, "h_p_%s_%s", TrkNames[7].c_str(), particle[i-8-nPVBin_].c_str());
631  sprintf(htit, "Momentum of %s tracks (%s)", TrkNames[7].c_str(), particle[i-8-nPVBin_].c_str());
632  }
633  h_p[i] = fs_->make<TH1D>(hname, htit, 400, 0, 200.0);
634  h_p[i]->Sumw2();
635 
636  if (i < 8) {
637  sprintf(hname, "h_eta_%s", TrkNames[i].c_str());
638  sprintf(htit, "Eta of %s tracks", TrkNames[i].c_str());
639  } else if (i < 8+nPVBin_) {
640  sprintf(hname, "h_eta_%s_%d", TrkNames[7].c_str(), i-8);
641  sprintf(htit, "Eta of %s tracks (PV=%d:%d)", TrkNames[7].c_str(), pvBin_[i-8], pvBin_[i-7]-1);
642  } else if (i >= nGen_) {
643  sprintf(hname, "h_eta_%s_%d", TrkNames[0].c_str(), i-nGen_);
644  sprintf(htit, "Eta of %s Generator tracks", TrkNames[0].c_str());
645  } else {
646  sprintf(hname, "h_eta_%s_%s", TrkNames[7].c_str(), particle[i-8-nPVBin_].c_str());
647  sprintf(htit, "Eta of %s tracks (%s)", TrkNames[7].c_str(), particle[i-8-nPVBin_].c_str());
648  }
649  h_eta[i] = fs_->make<TH1D>(hname, htit, 60, -3.0, 3.0);
650  h_eta[i]->Sumw2();
651 
652  if (i < 8) {
653  sprintf(hname, "h_phi_%s", TrkNames[i].c_str());
654  sprintf(htit, "Phi of %s tracks", TrkNames[i].c_str());
655  } else if (i < 8+nPVBin_) {
656  sprintf(hname, "h_phi_%s_%d", TrkNames[7].c_str(), i-8);
657  sprintf(htit, "Phi of %s tracks (PV=%d:%d)", TrkNames[7].c_str(), pvBin_[i-8], pvBin_[i-7]-1);
658  } else if (i >= nGen_) {
659  sprintf(hname, "h_phi_%s_%d", TrkNames[0].c_str(), i-nGen_);
660  sprintf(htit, "Phi of %s Generator tracks", TrkNames[0].c_str());
661  } else {
662  sprintf(hname, "h_phi_%s_%s", TrkNames[7].c_str(), particle[i-8-nPVBin_].c_str());
663  sprintf(htit, "Phi of %s tracks (%s)", TrkNames[7].c_str(), particle[i-8-nPVBin_].c_str());
664  }
665  h_phi[i] = fs_->make<TH1D>(hname, htit, 100, -3.15, 3.15);
666  h_phi[i]->Sumw2();
667  }
668  std::string IsolationNames[2] = {"Ecal", "Hcal"};
669  for (unsigned int i=0; i<2; i++) {
670  sprintf(hname, "h_maxNearP_%s", IsolationNames[i].c_str());
671  sprintf(htit, "Energy in ChargeIso region for %s", IsolationNames[i].c_str());
672  h_maxNearP[i] = fs_->make<TH1D>(hname, htit, 120, -1.5, 10.5);
673  h_maxNearP[i]->Sumw2();
674 
675  sprintf(hname, "h_ene1_%s", IsolationNames[i].c_str());
676  sprintf(htit, "Energy in smaller cone for %s", IsolationNames[i].c_str());
677  h_ene1[i] = fs_->make<TH1D>(hname, htit, 400, 0.0, 200.0);
678  h_ene1[i]->Sumw2();
679 
680  sprintf(hname, "h_ene2_%s", IsolationNames[i].c_str());
681  sprintf(htit, "Energy in bigger cone for %s", IsolationNames[i].c_str());
682  h_ene2[i] = fs_->make<TH1D>(hname, htit, 400, 0.0, 200.0);
683  h_ene2[i]->Sumw2();
684 
685  sprintf(hname, "h_ediff_%s", IsolationNames[i].c_str());
686  sprintf(htit, "Energy in NeutralIso region for %s", IsolationNames[i].c_str());
687  h_ediff[i] = fs_->make<TH1D>(hname, htit, 100, -0.5, 19.5);
688  h_ediff[i]->Sumw2();
689  }
690  std::string energyNames[6]={"E_{7x7}", "H_{3x3}", "(E_{7x7}+H_{3x3})",
691  "E_{11x11}", "H_{5x5}", "{E_{11x11}+H_{5x5})"};
692  for (int i=0; i<4+nPVBin_+4; ++i) {
693  for (int ip=0; ip<nPBin_; ++ip) {
694  for (int ie=0; ie<nEtaBin_; ++ie) {
695  for (int j=0; j<6; ++j) {
696  sprintf(hname, "h_energy_%d_%d_%d_%d", i, ip, ie, j);
697  if (i < 4) {
698  sprintf(htit,"%s/p (p=%4.1f:%4.1f; i#eta=%d:%d) for tracks with %s",
699  energyNames[j].c_str(),pBin_[ip],pBin_[ip+1],etaBin_[ie],
700  (etaBin_[ie+1]-1), TrkNames[i+4].c_str());
701  } else if (i < 4+nPVBin_) {
702  sprintf(htit,"%s/p (p=%4.1f:%4.1f, i#eta=%d:%d, PV=%d:%d) for tracks with %s",
703  energyNames[j].c_str(),pBin_[ip],pBin_[ip+1],etaBin_[ie],
704  (etaBin_[ie+1]-1), pvBin_[i-4], pvBin_[i-3],
705  TrkNames[7].c_str());
706  } else {
707  sprintf(htit,"%s/p (p=%4.1f:%4.1f, i#eta=%d:%d %s) for tracks with %s",
708  energyNames[j].c_str(),pBin_[ip],pBin_[ip+1],etaBin_[ie],
709  (etaBin_[ie+1]-1), particle[i-4-nPVBin_].c_str(),
710  TrkNames[7].c_str());
711  }
712  h_energy[i][ip][ie][j] = fs_->make<TH1D>(hname, htit, 500, -0.1, 4.9);
713  h_energy[i][ip][ie][j]->Sumw2();
714  }
715  }
716  }
717  }
718 
719  // Now the tree
720  if (doTree_) {
721  tree_ = fs_->make<TTree>("testTree", "new HLT Tree");
722  tree_->Branch("tr_goodRun", &tr_goodRun, "tr_goodRun/I");
723  tree_->Branch("tr_goodPV", &tr_goodPV, "tr_goodPV/I");
724  tree_->Branch("tr_eventWeight", &tr_eventWeight, "tr_eventWeight/D");
725  tree_->Branch("tr_tr_TrigName", &tr_TrigName);
726  tree_->Branch("tr_TrkPt", &tr_TrkPt);
727  tree_->Branch("tr_TrkP", &tr_TrkP);
728  tree_->Branch("tr_TrkEta", &tr_TrkEta);
729  tree_->Branch("tr_TrkPhi", &tr_TrkPhi);
730  tree_->Branch("tr_TrkID", &tr_TrkID);
731  tree_->Branch("tr_MaxNearP31X31", &tr_MaxNearP31X31);
732  tree_->Branch("tr_MaxNearHcalP7x7", &tr_MaxNearHcalP7x7);
733  tree_->Branch("tr_FE7x7P", &tr_FE7x7P);
734  tree_->Branch("tr_FE11x11P", &tr_FE11x11P);
735  tree_->Branch("tr_FE15x15P", &tr_FE15x15P);
736  tree_->Branch("tr_SE7x7P", &tr_SE7x7P);
737  tree_->Branch("tr_SE11x11P", &tr_SE11x11P);
738  tree_->Branch("tr_SE15x15P", &tr_SE15x15P);
739  tree_->Branch("tr_H3x3", &tr_H3x3);
740  tree_->Branch("tr_H5x5", &tr_H5x5);
741  tree_->Branch("tr_H7x7", &tr_H7x7);
742  tree_->Branch("tr_iEta", &tr_iEta);
743  }
744 }
745 
746 // ------------ method called when starting to processes a run ------------
747 void StudyHLT::beginRun(edm::Run const& iRun, edm::EventSetup const& iSetup) {
748  char hname[100], htit[400];
749  edm::LogInfo("IsoTrack") << "Run[" << nRun_ << "] " << iRun.run() << " hltconfig.init "
750  << hltConfig_.init(iRun,iSetup,"HLT",changed_);
751  sprintf(hname, "h_HLTAccepts_%i", iRun.run());
752  sprintf(htit, "HLT Accepts for Run No %i", iRun.run());
753  TH1I *hnew = fs_->make<TH1I>(hname, htit, 500, 0, 500);
754  for (int i=1; i<=500; ++i) hnew->GetXaxis()->SetBinLabel(i," ");
755  h_HLTAccepts.push_back(hnew);
756  edm::LogInfo("IsoTrack") << "beginrun " << iRun.run();
757  firstEvent_ = true;
758  changed_ = false;
759 }
760 
761 // ------------ method called when ending the processing of a run ------------
762 void StudyHLT::endRun(edm::Run const& iRun, edm::EventSetup const&) {
763  ++nRun_;
764  edm::LogInfo("IsoTrack") << "endrun[" << nRun_ << "] " << iRun.run();
765 }
766 
767 void StudyHLT::clear() {
768  tr_TrigName.clear();
769  tr_TrkPt.clear(); tr_TrkP.clear(); tr_TrkEta.clear(); tr_TrkPhi.clear();
770  tr_TrkID.clear(); tr_MaxNearP31X31.clear(); tr_MaxNearHcalP7x7.clear();
771  tr_FE7x7P.clear(); tr_FE11x11P.clear(); tr_FE15x15P.clear();
772  tr_SE7x7P.clear(); tr_SE11x11P.clear(); tr_SE15x15P.clear();
773  tr_H3x3.clear(); tr_H5x5.clear(); tr_H7x7.clear();
774  tr_iEta.clear();
775 }
776 
777 void StudyHLT::fillTrack(int i, double pt, double p, double eta, double phi){
778  h_pt[i]->Fill(pt,tr_eventWeight);
779  h_p[i]->Fill(p,tr_eventWeight);
780  h_eta[i]->Fill(eta,tr_eventWeight);
781  h_phi[i]->Fill(phi,tr_eventWeight);
782 }
783 
784 void StudyHLT::fillIsolation(int i, double emaxnearP, double eneutIso1, double eneutIso2){
785  h_maxNearP[i]->Fill(emaxnearP,tr_eventWeight);
786  h_ene1[i]->Fill(eneutIso1,tr_eventWeight);
787  h_ene2[i]->Fill(eneutIso2,tr_eventWeight);
788  h_ediff[i]->Fill(eneutIso2-eneutIso1,tr_eventWeight);
789 }
790 
791 void StudyHLT::fillEnergy(int flag, int ieta, double p, double enEcal1,
792  double enHcal1, double enEcal2, double enHcal2) {
793  int ip(-1), ie(-1);
794  for (int i=0; i<nPBin_; ++i) {
795  if (p >= pBin_[i] && p < pBin_[i+1]) { ip = i; break; }
796  }
797  for (int i=0; i<nEtaBin_; ++i) {
798  if (ieta >= etaBin_[i] && ieta < etaBin_[i+1]) { ie = i; break; }
799  }
800  if (ip >= 0 && ie >= 0 && enEcal1 > 0.02 && enHcal1 > 0.1) {
801  h_energy[flag][ip][ie][0]->Fill(enEcal1/p,tr_eventWeight);
802  h_energy[flag][ip][ie][1]->Fill(enHcal1/p,tr_eventWeight);
803  h_energy[flag][ip][ie][2]->Fill((enEcal1+enHcal1)/p,tr_eventWeight);
804  h_energy[flag][ip][ie][3]->Fill(enEcal2/p,tr_eventWeight);
805  h_energy[flag][ip][ie][4]->Fill(enHcal2/p,tr_eventWeight);
806  h_energy[flag][ip][ie][5]->Fill((enEcal2+enHcal2)/p,tr_eventWeight);
807  }
808 }
809 
810 std::string StudyHLT::truncate_str(const std::string& str) {
811  std::string truncated_str(str);
812  int length = str.length();
813  for (int i=0; i<length-2; i++){
814  if (str[i]=='_' && str[i+1]=='v' && isdigit(str.at(i+2))) {
815  int z = i+1;
816  truncated_str = str.substr(0,z);
817  }
818  }
819  return(truncated_str);
820 }
821 
822 int StudyHLT::trackPID(const reco::Track* pTrack, const edm::Handle<reco::GenParticleCollection> & genParticles) {
823  int id(0);
824  if (genParticles.isValid()) {
825  unsigned int indx;
826  reco::GenParticleCollection::const_iterator p;
827  double mindR(999.9);
828  for (p=genParticles->begin(),indx=0; p!=genParticles->end(); ++p,++indx) {
829  int pdgId = std::abs(p->pdgId());
830  int idx = (pdgId == 11) ? 1 :
831  ((pdgId == 211) ? 2 : ((pdgId == 321) ? 3 : ((pdgId == 2212) ? 4 : 0)));
832  if (idx > 0) {
833  double dEta = pTrack->eta() - p->momentum().Eta();
834  double phi1 = pTrack->phi();
835  double phi2 = p->momentum().Phi();
836  if (phi1 < 0) phi1 += 2.0*M_PI;
837  if (phi2 < 0) phi2 += 2.0*M_PI;
838  double dPhi = phi1-phi2;
839  if (dPhi > M_PI) dPhi -= 2.*M_PI;
840  else if (dPhi < -M_PI) dPhi += 2.*M_PI;
841  double dR = sqrt(dEta*dEta+dPhi*dPhi);
842  if (dR < mindR) {
843  mindR = dR; id = idx;
844  }
845  }
846  }
847  }
848  return id;
849 }
850 
851 DEFINE_FWK_MODULE(StudyHLT);
852 
edm::EventID::run
RunNumber_t run() const
Definition: EventID.h:39
StudyHLT::theTriggerResultsLabel_
const edm::InputTag theTriggerResultsLabel_
Definition: StudyHLT.cc:118
TFileService::kSharedResource
static const std::string kSharedResource
Definition: TFileService.h:76
StudyHLT::fillEnergy
void fillEnergy(int, int, double, double, double, double, double)
Definition: StudyHLT.cc:791
StudyHLT::tMinH_
const double tMinH_
Definition: StudyHLT.cc:115
reco::TrackBase::p
double p() const
momentum vector magnitude
Definition: TrackBase.h:615
particlelevel_cff.genEventInfo
genEventInfo
Definition: particlelevel_cff.py:15
reco::GenParticleCollection
std::vector< GenParticle > GenParticleCollection
collection of GenParticles
Definition: GenParticleFwd.h:10
jets_cff.quality
quality
Definition: jets_cff.py:21
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
edm::EventID::event
EventNumber_t event() const
Definition: EventID.h:41
StudyHLT::tok_EE_
edm::EDGetTokenT< EcalRecHitCollection > tok_EE_
Definition: StudyHLT.cc:130
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
StudyHLT::tr_TrkPt
std::vector< double > tr_TrkPt
Definition: StudyHLT.cc:149
StudyHLT::tr_H3x3
std::vector< double > tr_H3x3
Definition: StudyHLT.cc:151
spr::propagateCALO
std::vector< spr::propagatedTrackID > propagateCALO(edm::Handle< reco::TrackCollection > &trkCollection, const CaloGeometry *geo, const MagneticField *bField, const std::string &theTrackQuality, bool debug=false)
Definition: CaloPropagateTrack.cc:137
StudyHLT::theTrackQuality_
const std::string theTrackQuality_
Definition: StudyHLT.cc:113
StudyHLT::h_HLTCorr
TH1I * h_HLTCorr
Definition: StudyHLT.cc:134
StudyHLT::firstEvent_
bool firstEvent_
Definition: StudyHLT.cc:121
mps_fire.i
i
Definition: mps_fire.py:269
StudyHLT::h_numberPV
TH1I * h_numberPV
Definition: StudyHLT.cc:134
StudyHLT::h_nHLT
TH1I * h_nHLT
Definition: StudyHLT.cc:134
StudyHLT::triggerEvent_
const edm::InputTag triggerEvent_
Definition: StudyHLT.cc:118
StudyHLT::analyze
void analyze(edm::Event const &, edm::EventSetup const &) override
Definition: StudyHLT.cc:270
StudyHLT::tMinE_
const double tMinE_
Definition: StudyHLT.cc:115
edm::RunBase::run
RunNumber_t run() const
Definition: RunBase.h:40
edm::ParameterSetDescription::addUntracked
ParameterDescriptionBase * addUntracked(U const &iLabel, T const &value)
Definition: ParameterSetDescription.h:94
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
trigger::TriggerEvent
The single EDProduct to be saved for each event (AOD case)
Definition: TriggerEvent.h:25
StudyHLT::tr_MaxNearP31X31
std::vector< double > tr_MaxNearP31X31
Definition: StudyHLT.cc:150
ecalTB2006H4_GenSimDigiReco_cfg.bField
bField
Definition: ecalTB2006H4_GenSimDigiReco_cfg.py:235
StudyHLT::h_ntrk
TH1I * h_ntrk[2]
Definition: StudyHLT.cc:140
edm::Event::getByToken
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:519
reco::TrackBase::TrackQuality
TrackQuality
track quality
Definition: TrackBase.h:151
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
StudyHLT::h_HLTAccepts
std::vector< TH1I * > h_HLTAccepts
Definition: StudyHLT.cc:137
StudyHLT::tr_TrigName
std::vector< std::string > tr_TrigName
Definition: StudyHLT.cc:148
StudyHLT::nGen_
static const int nGen_
Definition: StudyHLT.cc:108
StudyHLT::tr_SE7x7P
std::vector< bool > tr_SE7x7P
Definition: StudyHLT.cc:153
edm::HLTGlobalStatus::accept
bool accept() const
Has at least one path accepted the event?
Definition: HLTGlobalStatus.h:52
StudyHLT::tok_ew_
edm::EDGetTokenT< GenEventInfoProduct > tok_ew_
Definition: StudyHLT.cc:132
reco::TrackCollection
std::vector< Track > TrackCollection
collection of Tracks
Definition: TrackFwd.h:14
StudyHLT::tok_parts_
edm::EDGetTokenT< reco::GenParticleCollection > tok_parts_
Definition: StudyHLT.cc:126
edm::EventBase::bunchCrossing
int bunchCrossing() const
Definition: EventBase.h:66
edm::EventBase::luminosityBlock
edm::LuminosityBlockNumber_t luminosityBlock() const
Definition: EventBase.h:63
StudyHLT::nPBin_
static const int nPBin_
Definition: StudyHLT.cc:107
reco::TrackBase::phi
double phi() const
azimuthal angle of momentum vector
Definition: TrackBase.h:645
StudyHLT::h_phi
TH1D * h_phi[nGen_+2]
Definition: StudyHLT.cc:139
TFileService::make
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
StudyHLT::h_HLTAccept
TH1I * h_HLTAccept
Definition: StudyHLT.cc:134
mix_E10TeV_L13E31_BX432_cfi.Lumi
Lumi
Definition: mix_E10TeV_L13E31_BX432_cfi.py:24
GenEventInfoProduct::weight
double weight() const
Definition: GenEventInfoProduct.h:31
reco::TrackBase::momentum
const Vector & momentum() const
track momentum vector
Definition: TrackBase.h:675
StudyHLT::HLTNames_
std::vector< std::string > HLTNames_
Definition: StudyHLT.cc:120
StudyHLT::beginLuminosityBlock
virtual void beginLuminosityBlock(edm::LuminosityBlock const &, edm::EventSetup const &)
Definition: StudyHLT.cc:95
StudyHLT::tMaxH_
const double tMaxH_
Definition: StudyHLT.cc:115
edm::TriggerNames::triggerNames
Strings const & triggerNames() const
Definition: TriggerNames.cc:24
StudyHLT::nPVBin_
static const int nPVBin_
Definition: StudyHLT.cc:107
StudyHLT::tr_TrkEta
std::vector< double > tr_TrkEta
Definition: StudyHLT.cc:149
StudyHLT::verbosity_
const int verbosity_
Definition: StudyHLT.cc:111
StudyHLT::tok_genTrack_
edm::EDGetTokenT< reco::TrackCollection > tok_genTrack_
Definition: StudyHLT.cc:127
StudyHLT::h_goodRun
TH1I * h_goodRun
Definition: StudyHLT.cc:135
StudyHLT::beginJob
void beginJob() override
Definition: StudyHLT.cc:584
StudyHLT::h_ene2
TH1D * h_ene2[2]
Definition: StudyHLT.cc:141
photonIsolationHIProducer_cfi.hbhe
hbhe
Definition: photonIsolationHIProducer_cfi.py:8
StudyHLT::fillIsolation
void fillIsolation(int, double, double, double)
Definition: StudyHLT.cc:784
StudyHLT::h_pt
TH1D * h_pt[nGen_+2]
Definition: StudyHLT.cc:138
StudyHLT::changed_
bool changed_
Definition: StudyHLT.cc:121
StudyHLT::tMaxE_
const double tMaxE_
Definition: StudyHLT.cc:115
StudyHLT::tok_hbhe_
edm::EDGetTokenT< HBHERecHitCollection > tok_hbhe_
Definition: StudyHLT.cc:131
StudyHLT::h_energy
TH1D * h_energy[nPVBin_+8][nPBin_][nEtaBin_][6]
Definition: StudyHLT.cc:142
spr::chargeIsolationEcal
double chargeIsolationEcal(unsigned int trkIndex, std::vector< spr::propagatedTrackID > &vdetIds, const CaloGeometry *geo, const CaloTopology *caloTopology, int ieta, int iphi, bool debug=false)
Definition: ChargeIsolation.cc:15
iEvent
int iEvent
Definition: GenABIO.cc:230
reco::TrackBase::eta
double eta() const
pseudorapidity of momentum vector
Definition: TrackBase.h:651
StudyHLT::~StudyHLT
~StudyHLT() override
Definition: StudyHLT.cc:86
StudyHLT::etaBin_
int etaBin_[nEtaBin_+1]
Definition: StudyHLT.cc:144
StudyHLT::tr_iEta
std::vector< int > tr_iEta
Definition: StudyHLT.cc:154
StudyHLT::selectionParameters_
spr::trackSelectionParameters selectionParameters_
Definition: StudyHLT.cc:119
StudyHLT::tok_recVtx_
edm::EDGetTokenT< reco::VertexCollection > tok_recVtx_
Definition: StudyHLT.cc:128
StudyHLT::fillTrack
void fillTrack(int, double, double, double, double)
Definition: StudyHLT.cc:777
StudyHLT::tr_MaxNearHcalP7x7
std::vector< double > tr_MaxNearHcalP7x7
Definition: StudyHLT.cc:150
StudyHLT::truncate_str
std::string truncate_str(const std::string &)
Definition: StudyHLT.cc:810
triggerObjects_cff.id
id
Definition: triggerObjects_cff.py:28
StudyHLT::puWeights_
const std::vector< double > puWeights_
Definition: StudyHLT.cc:117
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
edm::HLTGlobalStatus::size
unsigned int size() const
Get number of paths stored.
Definition: HLTGlobalStatus.h:39
StudyHLT::h_ene1
TH1D * h_ene1[2]
Definition: StudyHLT.cc:141
reco::TrackBase::pt
double pt() const
track transverse momentum
Definition: TrackBase.h:621
StudyHLT::nEtaBin_
static const int nEtaBin_
Definition: StudyHLT.cc:107
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
StudyHLT::tr_goodRun
int tr_goodRun
Definition: StudyHLT.cc:146
StudyHLT::hltConfig_
HLTConfigProvider hltConfig_
Definition: StudyHLT.cc:109
StudyHLT::tr_SE11x11P
std::vector< bool > tr_SE11x11P
Definition: StudyHLT.cc:153
StudyHLT::pvBin_
int pvBin_[nPVBin_+1]
Definition: StudyHLT.cc:144
triggerResults
static std::string const triggerResults
Definition: EdmProvDump.cc:41
create_public_lumi_plots.weights
weights
Definition: create_public_lumi_plots.py:1101
GenHFHadronMatcher_cfi.genParticles
genParticles
Definition: GenHFHadronMatcher_cfi.py:4
edm::ParameterSetDescription::add
ParameterDescriptionBase * add(U const &iLabel, T const &value)
Definition: ParameterSetDescription.h:89
edm::HandleBase::isValid
bool isValid() const
Definition: HandleBase.h:74
StudyHLT::h_ediff
TH1D * h_ediff[2]
Definition: StudyHLT.cc:141
StudyHLT::tr_goodPV
int tr_goodPV
Definition: StudyHLT.cc:146
StudyHLT::tr_FE11x11P
std::vector< double > tr_FE11x11P
Definition: StudyHLT.cc:152
M_PI
#define M_PI
Definition: BXVectorInputProducer.cc:51
gen::k
int k[5][pyjets_maxn]
Definition: Cascade2Hadronizer.cc:79
hiDetachedQuadStep_cff.pt1
pt1
Definition: hiDetachedQuadStep_cff.py:69
electrons_cff.bool
bool
Definition: electrons_cff.py:231
PFRecoTauDiscriminationAgainstElectronDeadECAL_cfi.dR
dR
Definition: PFRecoTauDiscriminationAgainstElectronDeadECAL_cfi.py:19
DetId
Definition: DetId.h:18
reco::TrackBase::qualityByName
static TrackQuality qualityByName(const std::string &name)
Definition: TrackBase.cc:125
StudyHLT::fs_
edm::Service< TFileService > fs_
Definition: StudyHLT.cc:110
StudyHLT::h_maxNearP
TH1D * h_maxNearP[2]
Definition: StudyHLT.cc:141
edm::Handle::product
T const * product() const
Definition: Handle.h:81
StudyHLT::clear
void clear()
Definition: StudyHLT.cc:767
StudyHLT::h_nHLTvsRN
TH2I * h_nHLTvsRN
Definition: StudyHLT.cc:136
common_cff.pdgId
pdgId
Definition: common_cff.py:54
StudyHLT::tok_EB_
edm::EDGetTokenT< EcalRecHitCollection > tok_EB_
Definition: StudyHLT.cc:129
StudyHLT::tr_TrkID
std::vector< int > tr_TrkID
Definition: StudyHLT.cc:154
StudyHLT::trackPID
int trackPID(const reco::Track *, const edm::Handle< reco::GenParticleCollection > &)
Definition: StudyHLT.cc:822
StudyHLT::maxTrackEta_
const double maxTrackEta_
Definition: StudyHLT.cc:114
HLTConfigProvider::init
bool init(const edm::Run &iRun, const edm::EventSetup &iSetup, const std::string &processName, bool &changed)
d&#39;tor
Definition: HLTConfigProvider.cc:42
edm::EventSetup::get
const T & get() const
Definition: EventSetup.h:59
training_settings.idx
idx
Definition: training_settings.py:16
StudyHLT::tree_
TTree * tree_
Definition: StudyHLT.cc:143
edm::ConfigurationDescriptions::add
void add(std::string const &label, ParameterSetDescription const &psetDescription)
Definition: ConfigurationDescriptions.cc:67
StudyHLT::h_goodPV
TH1I * h_goodPV
Definition: StudyHLT.cc:135
StudyHLT::tok_lumi
edm::EDGetTokenT< LumiDetails > tok_lumi
Definition: StudyHLT.cc:123
reco::TrackBase::quality
bool quality(const TrackQuality) const
Track quality.
Definition: TrackBase.h:510
StudyHLT::tr_H7x7
std::vector< double > tr_H7x7
Definition: StudyHLT.cc:151
edm::InputTag::label
std::string const & label() const
Definition: InputTag.h:36
StudyHLT::beginRun
void beginRun(edm::Run const &, edm::EventSetup const &) override
Definition: StudyHLT.cc:747
StudyHLT::tr_FE15x15P
std::vector< double > tr_FE15x15P
Definition: StudyHLT.cc:152
edm::EventBase::id
edm::EventID id() const
Definition: EventBase.h:60
edm
HLT enums.
Definition: AlignableModifier.h:17
p1
double p1[4]
Definition: TauolaWrapper.h:89
StudyHLT::tok_trigRes
edm::EDGetTokenT< edm::TriggerResults > tok_trigRes
Definition: StudyHLT.cc:125
StudyHLT::pBin_
double pBin_[nPBin_+1]
Definition: StudyHLT.cc:145
StudyHLT::tr_eventWeight
double tr_eventWeight
Definition: StudyHLT.cc:147
spr::trackSelectionParameters::minQuality
reco::TrackBase::TrackQuality minQuality
Definition: TrackSelection.h:22
StudyHLT::doTree_
const bool doTree_
Definition: StudyHLT.cc:116
StudyHLT::fillDescriptions
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
Definition: StudyHLT.cc:237
StudyHLT::h_eta
TH1D * h_eta[nGen_+2]
Definition: StudyHLT.cc:139
match
std::pair< typename Association::data_type::first_type, double > match(Reference key, Association association, bool bestMatchByMaxValue)
Generic matching function.
Definition: Utils.h:10
StudyHLT::minTrackP_
const double minTrackP_
Definition: StudyHLT.cc:114
StudyHLT::StudyHLT
StudyHLT(const edm::ParameterSet &)
Definition: StudyHLT.cc:157
StudyHLT::tr_FE7x7P
std::vector< double > tr_FE7x7P
Definition: StudyHLT.cc:152
StudyHLT::endLuminosityBlock
virtual void endLuminosityBlock(edm::LuminosityBlock const &, edm::EventSetup const &)
Definition: StudyHLT.cc:96
StudyHLT::tr_TrkPhi
std::vector< double > tr_TrkPhi
Definition: StudyHLT.cc:149
StudyHLT::isItAOD_
const bool isItAOD_
Definition: StudyHLT.cc:116
spr::chargeIsolationHcal
double chargeIsolationHcal(unsigned int trkIndex, std::vector< spr::propagatedTrackID > &vdetIds, const HcalTopology *topology, int ieta, int iphi, bool debug=false)
Definition: ChargeIsolation.cc:97
StudyHLT::tr_TrkP
std::vector< double > tr_TrkP
Definition: StudyHLT.cc:149
edm::ESHandle::product
T const * product() const
Definition: ESHandle.h:86
StudyHLT::tr_SE15x15P
std::vector< bool > tr_SE15x15P
Definition: StudyHLT.cc:153
StudyHLT::trigNames_
const std::vector< std::string > trigNames_
Definition: StudyHLT.cc:112
edm::Event::triggerNames
edm::TriggerNames const & triggerNames(edm::TriggerResults const &triggerResults) const override
Definition: Event.cc:301
StudyHLT::h_p
TH1D * h_p[nGen_+2]
Definition: StudyHLT.cc:138
photonIsolationHIProducer_cfi.trackQuality
trackQuality
Definition: photonIsolationHIProducer_cfi.py:14
StudyHLT::tr_H5x5
std::vector< double > tr_H5x5
Definition: StudyHLT.cc:151
StudyHLT::endRun
void endRun(edm::Run const &, edm::EventSetup const &) override
Definition: StudyHLT.cc:762
edm::Run
Definition: Run.h:43
StudyHLT::nRun_
int nRun_
Definition: StudyHLT.cc:144
StudyHLT::tok_trigEvt
edm::EDGetTokenT< trigger::TriggerEvent > tok_trigEvt
Definition: StudyHLT.cc:124
spr::eECALmatrix
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)
StudyHLT::newNames_
const std::vector< std::string > newNames_
Definition: StudyHLT.cc:112
spr::eHCALmatrix
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, int useRaw=0, bool debug=false)
Generated for CMSSW Reference Manual by doxygen 1.8.11