CMS 3D CMS Logo

IsoTrackCalib.cc
Go to the documentation of this file.
1 // system include files
2 #include <memory>
3 
4 // Root objects
5 #include "TROOT.h"
6 #include "TSystem.h"
7 #include "TFile.h"
8 #include "TProfile.h"
9 #include "TDirectory.h"
10 #include "TTree.h"
11 #include "TLorentzVector.h"
12 #include "TInterpreter.h"
13 
19 
20  //L1 trigger Menus etc
26 
30 
31 //Tracks
36 // Vertices
40 // Jets
44 
45 //Triggers
52 
62 
65 
74 
77 
81 
82 class IsoTrackCalib : public edm::one::EDAnalyzer<edm::one::WatchRuns,edm::one::SharedResources> {
83 
84 public:
85  explicit IsoTrackCalib(const edm::ParameterSet&);
86  ~IsoTrackCalib() override;
87 
88  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
89 
90 private:
91  void beginJob() override ;
92  void analyze(const edm::Event&, const edm::EventSetup&) override;
93  void endJob() override {}
94  void beginRun(edm::Run const&, edm::EventSetup const&) override;
95  void endRun(edm::Run const&, edm::EventSetup const&) override;
96 
100  double deltaR(double eta1,double eta2,double phi1,double phi2);
101 
106  const int verbosity_;
107  const std::vector<std::string> l1Names_;
110  const double a_coneR_, a_charIsoR_, a_mipR_;
111  std::vector<bool> *t_l1bits;
112 
125 
126  TTree *tree;
131 
133  std::vector<unsigned int> *t_DetIds;
134  std::vector<double> *t_HitEnergies, pbin;
136  TH1I *h_iEta, *h_tkEta0[5], *h_tkEta1[5], *h_tkEta2[5];
137  TH1I *h_tkEta3[5], *h_tkEta4[5], *h_tkEta5[5];
139  TH1F *h_jetpt[4];
140  TH1I *h_tketa0[6], *h_tketa1[6], *h_tketa2[6];
141  TH1I *h_tketa3[6], *h_tketa4[6],*h_tketa5[6];
142  std::map< std::pair<unsigned int,std::string>, int> l1AlgoMap_;
143 };
144 
145 static const bool useL1GtTriggerMenuLite(true);
148  verbosity_(iConfig.getUntrackedParameter<int>("Verbosity",0)),
149  l1Names_(iConfig.getUntrackedParameter<std::vector<std::string> >("L1Seed")),
150  theTrackQuality_(iConfig.getUntrackedParameter<std::string>("TrackQuality","highPurity")),
151  a_coneR_(iConfig.getUntrackedParameter<double>("ConeRadius",34.98)),
152  a_charIsoR_(a_coneR_+28.9),
153  a_mipR_(iConfig.getUntrackedParameter<double>("ConeRadiusMIP",14.0)),
154  tok_genTrack_(consumes<reco::TrackCollection>(edm::InputTag("generalTracks"))),
155  tok_recVtx_(consumes<reco::VertexCollection>(edm::InputTag("offlinePrimaryVertices"))),
156  tok_bs_(consumes<reco::BeamSpot>(edm::InputTag("offlineBeamSpot"))),
157  tok_ew_(consumes<GenEventInfoProduct>(edm::InputTag("generatorSmeared"))),
158  tok_jets_(consumes<reco::GenJetCollection>(iConfig.getParameter<edm::InputTag>("JetSource"))),
159  tok_pfjets_(consumes<reco::PFJetCollection>(edm::InputTag("ak5PFJets"))) {
160 
161  usesResource(TFileService::kSharedResource);
162 
163  //now do whatever initialization is needed
165  selectionParameters_.minPt = iConfig.getUntrackedParameter<double>("MinTrackPt", 10.0);
166  selectionParameters_.minQuality = trackQuality_;
167  selectionParameters_.maxDxyPV = iConfig.getUntrackedParameter<double>("MaxDxyPV", 0.2);
168  selectionParameters_.maxDzPV = iConfig.getUntrackedParameter<double>("MaxDzPV", 5.0);
169  selectionParameters_.maxChi2 = iConfig.getUntrackedParameter<double>("MaxChi2", 5.0);
170  selectionParameters_.maxDpOverP = iConfig.getUntrackedParameter<double>("MaxDpOverP", 0.1);
171  selectionParameters_.minOuterHit = iConfig.getUntrackedParameter<int>("MinOuterHit", 4);
172  selectionParameters_.minLayerCrossed= iConfig.getUntrackedParameter<int>("MinLayerCrossed", 8);
173  selectionParameters_.maxInMiss = iConfig.getUntrackedParameter<int>("MaxInMiss", 0);
174  selectionParameters_.maxOutMiss = iConfig.getUntrackedParameter<int>("MaxOutMiss", 0);
175  bool isItAOD = iConfig.getUntrackedParameter<bool>("IsItAOD", false);
176  edm::InputTag L1extraTauJetSource_ = iConfig.getParameter<edm::InputTag> ("L1extraTauJetSource" );
177  edm::InputTag L1extraCenJetSource_ = iConfig.getParameter<edm::InputTag> ("L1extraCenJetSource" );
178  edm::InputTag L1extraFwdJetSource_ = iConfig.getParameter<edm::InputTag> ("L1extraFwdJetSource" );
179 
180  // define tokens for access
181  if (isItAOD) {
182  tok_EB_ = consumes<EcalRecHitCollection>(edm::InputTag("reducedEcalRecHitsEB"));
183  tok_EE_ = consumes<EcalRecHitCollection>(edm::InputTag("reducedEcalRecHitsEE"));
184  tok_hbhe_ = consumes<HBHERecHitCollection>(edm::InputTag("reducedHcalRecHits", "hbhereco"));
185  } else {
186  tok_EB_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit","EcalRecHitsEB"));
187  tok_EE_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit","EcalRecHitsEE"));
188  tok_hbhe_ = consumes<HBHERecHitCollection>(edm::InputTag("hbhereco"));
189  }
190  tok_L1extTauJet_ = consumes<l1extra::L1JetParticleCollection>(L1extraTauJetSource_);
191  tok_L1extCenJet_ = consumes<l1extra::L1JetParticleCollection>(L1extraCenJetSource_);
192  tok_L1extFwdJet_ = consumes<l1extra::L1JetParticleCollection>(L1extraFwdJetSource_);
193  if (verbosity_>=0) {
194  edm::LogInfo("IsoTrack")
195  <<"Parameters read from config file \n"
196  <<"\t minPt " << selectionParameters_.minPt
197  <<"\t theTrackQuality " << theTrackQuality_
198  <<"\t minQuality " << selectionParameters_.minQuality
199  <<"\t maxDxyPV " << selectionParameters_.maxDxyPV
200  <<"\t maxDzPV " << selectionParameters_.maxDzPV
201  <<"\t maxChi2 " << selectionParameters_.maxChi2
202  <<"\t maxDpOverP " << selectionParameters_.maxDpOverP
203  <<"\t minOuterHit " << selectionParameters_.minOuterHit
204  <<"\t minLayerCrossed " << selectionParameters_.minLayerCrossed
205  <<"\t maxInMiss " << selectionParameters_.maxInMiss
206  <<"\t maxOutMiss " << selectionParameters_.maxOutMiss
207  <<"\t a_coneR " << a_coneR_
208  <<"\t a_charIsoR " << a_charIsoR_
209  <<"\t a_mipR " << a_mipR_
210  <<"\t isItAOD " << isItAOD;
211  edm::LogInfo("IsoTrack") << l1Names_.size() << " triggers to be studied";
212  for (unsigned int k=0; k<l1Names_.size(); ++k)
213  edm::LogInfo("IsoTrack") << "[" << k << "]: " << l1Names_[k];
214  }
215 
216 }
217 
219  // do anything here that needs to be done at desctruction time
220  // (e.g. close files, deallocate resources etc.)
221 
222 }
223 
224 // ------------ method fills 'descriptions' with the allowed parameters for the module ------------
226 
227  std::vector<std::string> seeds = {"L1_SingleJet36","L1_SingleJet52",
228  "L1_SingleJet68","L1_SingleJet92",
229  "L1_SingleJet128"};
231  desc.addUntracked<int>("Verbosity", 0);
232  desc.addUntracked<std::vector<std::string> >("L1Seed", seeds);
233  desc.addUntracked<std::string>("TrackQuality", "highPurity");
234  desc.addUntracked<double>("MinTrackPt", 10.0);
235  desc.addUntracked<double>("MaxDxyPV", 0.02);
236  desc.addUntracked<double>("MaxDzPV", 0.02);
237  desc.addUntracked<double>("MaxChi2", 5.0);
238  desc.addUntracked<double>("MaxDpOverP", 0.1);
239  desc.addUntracked<int>("MinOuterHit", 4);
240  desc.addUntracked<int>("MinLayerCrossed", 8);
241  desc.addUntracked<int>("MaxInMiss", 0);
242  desc.addUntracked<int>("MaxOutMiss", 0);
243  desc.addUntracked<double>("ConeRadius", 34.98);
244  desc.addUntracked<double>("ConeRadiusMIP",14.0);
245  desc.addUntracked<bool>("IsItAOD", false);
246  descriptions.add("isoTrackCalib",desc);
247 }
248 
250  const edm::EventSetup& iSetup) {
251 
252  t_Run = iEvent.id().run();
253  t_Event = iEvent.id().event();
254  if (verbosity_%10 > 0)
255  edm::LogInfo("IsoTrack")
256  << "Run " << t_Run << " Event " << t_Event << " Luminosity "
257  << iEvent.luminosityBlock() << " Bunch " << iEvent.bunchCrossing()
258  << " starts ==========";
259 
260  //Get magnetic field and ECAL channel status
262  iSetup.get<IdealMagneticFieldRecord>().get(bFieldH);
263  const MagneticField *bField = bFieldH.product();
264 
266  iSetup.get<EcalSeverityLevelAlgoRcd>().get(sevlv);
267 
268  // get handles to calogeometry and calotopology
270  iSetup.get<CaloGeometryRecord>().get(pG);
271  const CaloGeometry* geo = pG.product();
272 
273  //Get track collection
275  iEvent.getByToken(tok_genTrack_, trkCollection);
276 
277  //event weight for FLAT sample
278  t_EventWeight = 1.0;
280  iEvent.getByToken(tok_ew_, genEventInfo);
281  if (genEventInfo.isValid()) t_EventWeight = genEventInfo->weight();
282 
283  // genJet information
285  iEvent.getByToken(tok_jets_, genJets);
286  if (genJets.isValid()) {
287  for (unsigned iGenJet = 0; iGenJet < genJets->size(); ++iGenJet) {
288  const reco::GenJet& genJet = (*genJets) [iGenJet];
289  double genJetPt = genJet.pt();
290  double genJetEta = genJet.eta();
291  h_jetpt[0]->Fill(genJetPt);
292  h_jetpt[1]->Fill(genJetPt,t_EventWeight);
293  if (genJetEta>-2.5 && genJetEta<2.5) {
294  h_jetpt[2]->Fill(genJetPt);
295  h_jetpt[3]->Fill(genJetPt,t_EventWeight);
296  }
297  break;
298  }
299  }
300 
301  //pf jets
303  iEvent.getByToken(tok_pfjets_, pfJets);
304 
305  //Define the best vertex and the beamspot
307  iEvent.getByToken(tok_recVtx_, recVtxs);
308  edm::Handle<reco::BeamSpot> beamSpotH;
309  iEvent.getByToken(tok_bs_, beamSpotH);
310  math::XYZPoint leadPV(0,0,0);
311  if (!recVtxs->empty() && !((*recVtxs)[0].isFake())) {
312  leadPV = math::XYZPoint( (*recVtxs)[0].x(),(*recVtxs)[0].y(), (*recVtxs)[0].z() );
313  } else if (beamSpotH.isValid()) {
314  leadPV = beamSpotH->position();
315  }
316  if (verbosity_>10) {
317  if ((verbosity_%100)/10>2)
318  edm::LogInfo("IsoTrack") << "Primary Vertex " << leadPV;
319  if (beamSpotH.isValid()) edm::LogInfo("IsoTrack") << " Beam Spot "
320  << beamSpotH->position();
321  }
322 
323  // RecHits
324  edm::Handle<EcalRecHitCollection> barrelRecHitsHandle;
325  edm::Handle<EcalRecHitCollection> endcapRecHitsHandle;
326  iEvent.getByToken(tok_EB_, barrelRecHitsHandle);
327  iEvent.getByToken(tok_EE_, endcapRecHitsHandle);
329  iEvent.getByToken(tok_hbhe_, hbhe);
331 
332  for (rhitItr=hbhe->begin();rhitItr!=hbhe->end();rhitItr++) {
333  double rec_energy = rhitItr->energy();
334  int rec_ieta = rhitItr->id().ieta();
335  int rec_depth = rhitItr->id().depth();
336  int rec_zside = rhitItr->id().zside();
337  double num1_1 = rec_zside*(rec_ieta+0.2*(rec_depth-1));
338  if (verbosity_%10>0)
339  edm::LogInfo("IsoTrack") << "detid/rechit/ieta/zside/depth/num " << " = "
340  << rhitItr->id() << "/" << rec_energy << "/"
341  << rec_ieta << "/" << rec_zside << "/"
342  << rec_depth << "/" << num1_1;
343  h_iEta->Fill(rec_ieta);
344  h_Rechit_E->Fill(rec_energy);
345  h_RecHit_iEta->Fill(rec_ieta,rec_energy);
346  h_RecHit_num->Fill(num1_1,rec_energy);
347  }
348 
349  //Propagate tracks to calorimeter surface)
350  std::vector<spr::propagatedTrackDirection> trkCaloDirections;
351  spr::propagateCALO(trkCollection, geo, bField, theTrackQuality_,
352  trkCaloDirections, ((verbosity_/100)%10>2));
353  std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
354  for (trkDetItr = trkCaloDirections.begin();
355  trkDetItr != trkCaloDirections.end(); trkDetItr++) {
356  if (trkDetItr->okHCAL) {
357  HcalDetId detId = (HcalDetId)(trkDetItr->detIdHCAL);
358  int tk_ieta = detId.ieta();
359  const reco::Track* pTrack = &(*(trkDetItr->trkItr));
360  double tk_p = pTrack->p();
361  h_tketa0[0]->Fill(tk_ieta);
362  for (unsigned int k=1; k<pbin.size(); ++k) {
363  if (tk_p >= pbin[k-1] && tk_p < pbin[k]) {
364  h_tketa0[k]->Fill(tk_ieta);
365  break;
366  }
367  }
368  }
369  }
371  t_l1bits->clear();
372  for (unsigned int i=0; i<l1Names_.size(); ++i) t_l1bits->push_back(false);
373  bool useL1EventSetup = true;
374  bool useL1GtTriggerMenuLite = true;
375 
376  m_l1GtUtils.getL1GtRunCache(iEvent, iSetup, useL1EventSetup, useL1GtTriggerMenuLite);
377  int iErrorCode = -1;
378  int l1ConfCode = -1;
379  const bool l1Conf = m_l1GtUtils.availableL1Configuration(iErrorCode, l1ConfCode);
380  if( !l1Conf) {
381  edm::LogInfo("IsoTrack") << "\nL1 configuration code:" << l1ConfCode
382  << "\nNo valid L1 trigger configuration available."
383  << "\nSee text above for error code interpretation"
384  << "\nNo return here, in order to test each method"
385  << ", protected against configuration error.";
386  }
387 
388  const AlgorithmMap& algorithmMap = m_l1GtMenu->gtAlgorithmMap();
389  const std::string& menuName = m_l1GtMenu->gtTriggerMenuName();
390  if (verbosity_%10>0)
391  edm::LogInfo("IsoTrack") << "menuName " << menuName << std::endl;
392 
393  std::vector<int> algbits;
394  for (CItAlgo itAlgo = algorithmMap.begin(); itAlgo != algorithmMap.end();
395  itAlgo++) {
396  std::string algName = itAlgo->first;
397  int algBitNumber = ( itAlgo->second ).algoBitNumber();
398  bool decision = m_l1GtUtils.decision(iEvent, itAlgo->first, iErrorCode);
399 
400  bool l1ok(false);
401  if (verbosity_%10>0)
402  edm::LogInfo("IsoTrack") << algName <<" "<< algBitNumber << " "
403  << decision;
404  for (unsigned int i=0; i<l1Names_.size(); ++i) {
405  if (algName.find(l1Names_[i])!=std::string::npos){
406  if (verbosity_%10>0)
407  edm::LogInfo("IsoTrack") << "match found" << " " << algName << " "
408  << decision;
409  t_l1bits->at(i) = (decision>0);
410  if (decision > 0) l1ok = true;
411  }
412  }
413  if (verbosity_%10>0) edm::LogInfo("IsoTrack") << "l1 ok =" << l1ok;
414 
415  if(l1ok){
417  iEvent.getByToken(tok_L1extTauJet_,l1TauHandle);
418  l1extra::L1JetParticleCollection::const_iterator itr;
419  double ptTriggered = -10;
420  double etaTriggered = -100;
421  double phiTriggered = -100;
422 
423  for(itr = l1TauHandle->begin(); itr != l1TauHandle->end(); ++itr) {
424  if(itr->pt()>ptTriggered){
425  ptTriggered = itr->pt();
426  etaTriggered = itr->eta();
427  phiTriggered = itr->phi() ;
428  }
429  if (verbosity_%10>0)
430  edm::LogInfo("IsoTrack") << "tauJ pt " << itr->pt() << " eta/phi "
431  << itr->eta() << " " << itr->phi();
432  }
434  iEvent.getByToken(tok_L1extCenJet_,l1CenJetHandle);
435  for (itr = l1CenJetHandle->begin(); itr != l1CenJetHandle->end(); ++itr){
436  if (itr->pt()>ptTriggered) {
437  ptTriggered = itr->pt();
438  etaTriggered = itr->eta();
439  phiTriggered = itr->phi() ;
440  }
441  if (verbosity_%10>0)
442  edm::LogInfo("IsoTrack") << "cenJ pt " << itr->pt()
443  << " eta/phi " << itr->eta() << " "
444  << itr->phi();
445  }
447  iEvent.getByToken(tok_L1extFwdJet_,l1FwdJetHandle);
448  for (itr = l1FwdJetHandle->begin(); itr != l1FwdJetHandle->end(); ++itr) {
449  if (itr->pt()>ptTriggered) {
450  ptTriggered = itr->pt();
451  etaTriggered = itr->eta();
452  phiTriggered = itr->phi() ;
453  }
454  if (verbosity_%10>0)
455  edm::LogInfo("IsoTrack") << "forJ pt " << itr->pt() << " eta/phi "
456  << itr->eta() << " " << itr->phi();
457  }
458  if (verbosity_%10>0)
459  edm::LogInfo("IsoTrack") << "jets pt/eta/phi = " << ptTriggered << "/"
460  << etaTriggered << "/" <<phiTriggered;
462  unsigned int nTracks(0),nselTracks(0);
463  for (trkDetItr = trkCaloDirections.begin(),nTracks=0;
464  trkDetItr != trkCaloDirections.end(); trkDetItr++,nTracks++) {
465  const reco::Track* pTrack = &(*(trkDetItr->trkItr));
466  math::XYZTLorentzVector v4(pTrack->px(), pTrack->py(),
467  pTrack->pz(), pTrack->p());
468 
469  t_mindR1 = deltaR(etaTriggered,v4.eta(),phiTriggered,v4.phi());
470  t_mindR2 = -999;
471  if (verbosity_%10>0)
472  edm::LogInfo("IsoTrack") << "This track : " << nTracks << " (pt/eta/phi/p) :"
473  << pTrack->pt() << "/" << pTrack->eta() << "/"
474  << pTrack->phi() << "/" << pTrack->p();
475 
476  if (verbosity_%10>0) edm::LogInfo("IsoTrack") << "dr values are = " << t_mindR1;
477 
478  t_l1pt = ptTriggered;
479  t_l1eta = etaTriggered;
480  t_l1phi = phiTriggered;
481  t_l3pt = -999;
482  t_l3eta = -999;
483  t_l3phi = -999;
484 
485  //Selection of good track
486  t_selectTk = spr::goodTrack(pTrack,leadPV,selectionParameters_,((verbosity_/100)%10>2));
488  oneCutParameters.maxDxyPV = 10;
489  oneCutParameters.maxDzPV = 100;
490  oneCutParameters.maxInMiss = 2;
491  oneCutParameters.maxOutMiss= 2;
492  bool qltyFlag = spr::goodTrack(pTrack,leadPV,oneCutParameters,((verbosity_/100)%10>2));
493  oneCutParameters = selectionParameters_;
494  oneCutParameters.maxDxyPV = 10;
495  oneCutParameters.maxDzPV = 100;
496  t_qltyMissFlag = spr::goodTrack(pTrack,leadPV,oneCutParameters,((verbosity_/100)%10>2));
497  oneCutParameters = selectionParameters_;
498  oneCutParameters.maxInMiss = 2;
499  oneCutParameters.maxOutMiss= 2;
500  t_qltyPVFlag = spr::goodTrack(pTrack,leadPV,oneCutParameters,((verbosity_/100)%10>2));
501  t_ieta = 0;
502  if (trkDetItr->okHCAL) {
503  HcalDetId detId = (HcalDetId)(trkDetItr->detIdHCAL);
504  t_ieta = detId.ieta();
505  }
506  if (verbosity_%10 > 0)
507  edm::LogInfo("IsoTrack") << "qltyFlag|okECAL|okHCAL : " << qltyFlag << "|"
508  << trkDetItr->okECAL << "/" << trkDetItr->okHCAL;
509  t_qltyFlag = (qltyFlag && trkDetItr->okECAL && trkDetItr->okHCAL);
510  t_p = pTrack->p();
511  h_tketa1[0]->Fill(t_ieta);
512  for (unsigned int k=1; k<pbin.size(); ++k) {
513  if (t_p >= pbin[k-1] && t_p < pbin[k]) {
514  h_tketa1[k]->Fill(t_ieta);
515  break;
516  }
517  }
518  if (t_qltyFlag) {
519  nselTracks++;
520  h_tketa2[0]->Fill(t_ieta);
521  for (unsigned int k=1; k<pbin.size(); ++k) {
522  if (t_p >= pbin[k-1] && t_p < pbin[k]) {
523  h_tketa2[k]->Fill(t_ieta);
524  break;
525  }
526  }
527  int nRH_eMipDR(0), nNearTRKs(0), nRecHits(-999);
528  t_eMipDR = spr::eCone_ecal(geo, barrelRecHitsHandle,
529  endcapRecHitsHandle,
530  trkDetItr->pointHCAL,
531  trkDetItr->pointECAL,
532  a_mipR_, trkDetItr->directionECAL,
533  nRH_eMipDR);
534  t_DetIds->clear(); t_HitEnergies->clear();
535  std::vector<DetId> ids;
536  t_eHcal = spr::eCone_hcal(geo, hbhe, trkDetItr->pointHCAL,
537  trkDetItr->pointECAL, a_coneR_,
538  trkDetItr->directionHCAL,nRecHits,
539  ids, *t_HitEnergies);
540  for (unsigned int k=0; k<ids.size(); ++k) {
541  t_DetIds->push_back(ids[k].rawId());
542  }
543  t_hmaxNearP = spr::chargeIsolationCone(nTracks,trkCaloDirections,
544  a_charIsoR_, nNearTRKs,
545  ((verbosity_/100)%10>2));
546  if (t_hmaxNearP < 2) {
547  h_tketa3[0]->Fill(t_ieta);
548  for (unsigned int k=1; k<pbin.size(); ++k) {
549  if (t_p >= pbin[k-1] && t_p < pbin[k]) {
550  h_tketa3[k]->Fill(t_ieta);
551  break;
552  }
553  }
554  if (t_eMipDR < 1) {
555  h_tketa4[0]->Fill(t_ieta);
556  for (unsigned int k=1; k<pbin.size(); ++k) {
557  if (t_p >= pbin[k-1] && t_p < pbin[k]) {
558  h_tketa4[k]->Fill(t_ieta);
559  break;
560  }
561  }
562  if (t_mindR1 > 1) {
563  h_tketa5[0]->Fill(t_ieta);
564  for (unsigned int k=1; k<pbin.size(); ++k) {
565  if (t_p >= pbin[k-1] && t_p < pbin[k]) {
566  h_tketa5[k]->Fill(t_ieta);
567  break;
568  }
569  }
570  }
571  }
572  }
573  if (verbosity_%10 > 0) {
574  edm::LogInfo("IsoTrack") << "This track : " << nTracks << " (pt/eta/phi/p) :"
575  << pTrack->pt() << "/" << pTrack->eta() << "/"
576  << pTrack->phi() << "/" << t_p;
577  edm::LogInfo("IsoTrack") << "e_MIP " << t_eMipDR << " Chg Isolation "
578  << t_hmaxNearP << " eHcal" << t_eHcal << " ieta "
579  << t_ieta << " Quality " << t_qltyMissFlag
580  << ":" << t_qltyPVFlag << ":" << t_selectTk;
581  for (unsigned int lll=0;lll<t_DetIds->size();lll++) {
582  edm::LogInfo("IsoTrack") << "det id is = " << t_DetIds->at(lll) << " "
583  << " hit enery is = " << t_HitEnergies->at(lll);
584  }
585  }
586  if (t_p>20.0 && t_eMipDR<2.0 && t_hmaxNearP<10.0) {
587  tree->Fill();
588  }
589  }
590  }
591  }
592  }
593 }
594 
596  h_RecHit_iEta = fs_->make<TProfile>("rechit_ieta","Rec hit vs. ieta",60,-30,30,0,1000);
597  h_RecHit_num = fs_->make<TProfile>("rechit_num","Rec hit vs. num",100,0,20,0,1000);
598  h_iEta = fs_->make<TH1I>("iEta","iEta",60,-30,30);
599  h_Rechit_E = fs_->make<TH1F>("Rechit_E","Rechit_E",100,0,1000);
600 
601  double prange[5] = {20,30,40,60,100};
602  for (int k=0; k<5; ++k) pbin.push_back(prange[k]);
603  std::string type[6] = {"All", "Trigger OK", "Tree Selected",
604  "Charge Isolation", "MIP Cut", "L1 Cut"};
605  for (unsigned int k=0; k<pbin.size(); ++k) {
606  char name[20], namp[20], title[100];
607  if (k == 0) sprintf (namp, "all momentum");
608  else sprintf (namp, "p = %4.0f:%4.0f GeV", pbin[k-1], pbin[k]);
609  sprintf (name, "TrackEta0%d", k);
610  sprintf (title, "Track #eta for tracks with %s (%s)",namp,type[0].c_str());
611  h_tketa0[k] = fs_->make<TH1I>(name, title, 60, -30, 30);
612  sprintf (name, "TrackEta1%d", k);
613  sprintf (title, "Track #eta for tracks with %s (%s)",namp,type[1].c_str());
614  h_tketa1[k] = fs_->make<TH1I>(name, title, 60, -30, 30);
615  sprintf (name, "TrackEta2%d", k);
616  sprintf (title, "Track #eta for tracks with %s (%s)",namp,type[2].c_str());
617  h_tketa2[k] = fs_->make<TH1I>(name, title, 60, -30, 30);
618  sprintf (name, "TrackEta3%d", k);
619  sprintf (title, "Track #eta for tracks with %s (%s)",namp,type[3].c_str());
620  h_tketa3[k] = fs_->make<TH1I>(name, title, 60, -30, 30);
621  sprintf (name, "TrackEta4%d", k);
622  sprintf (title, "Track #eta for tracks with %s (%s)",namp,type[4].c_str());
623  h_tketa4[k] = fs_->make<TH1I>(name, title, 60, -30, 30);
624  sprintf (name, "TrackEta5%d", k);
625  sprintf (title, "Track #eta for tracks with %s (%s)",namp,type[5].c_str());
626  h_tketa5[k] = fs_->make<TH1I>(name, title, 60, -30, 30);
627  }
628  h_jetpt[0] = fs_->make<TH1F>("Jetpt0","Jet p_T (All)", 500,0.,2500.);
629  h_jetpt[1] = fs_->make<TH1F>("Jetpt1","Jet p_T (All Weighted)", 500,0.,2500.);
630  h_jetpt[2] = fs_->make<TH1F>("Jetpt2","Jet p_T (|#eta| < 2.5)", 500,0.,2500.);
631  h_jetpt[3] = fs_->make<TH1F>("Jetpt3","Jet p_T (|#eta| < 2.5 Weighted)", 500,0.,2500.);
632 
633  tree = fs_->make<TTree>("CalibTree", "CalibTree");
634 
635  tree->Branch("t_Run", &t_Run, "t_Run/I");
636  tree->Branch("t_Event", &t_Event, "t_Event/I");
637  tree->Branch("t_ieta", &t_ieta, "t_ieta/I");
638  tree->Branch("t_EventWeight", &t_EventWeight, "t_EventWeight/D");
639  tree->Branch("t_l1pt", &t_l1pt, "t_l1pt/D");
640  tree->Branch("t_l1eta", &t_l1eta, "t_l1eta/D");
641  tree->Branch("t_l1phi", &t_l1phi, "t_l1phi/D");
642  tree->Branch("t_l3pt", &t_l3pt, "t_l3pt/D");
643  tree->Branch("t_l3eta", &t_l3eta, "t_l3eta/D");
644  tree->Branch("t_l3phi", &t_l3phi, "t_l3phi/D");
645  tree->Branch("t_p", &t_p, "t_p/D");
646  tree->Branch("t_mindR1", &t_mindR1, "t_mindR1/D");
647  tree->Branch("t_mindR2", &t_mindR2, "t_mindR2/D");
648  tree->Branch("t_eMipDR", &t_eMipDR, "t_eMipDR/D");
649  tree->Branch("t_eHcal", &t_eHcal, "t_eHcal/D");
650  tree->Branch("t_hmaxNearP", &t_hmaxNearP, "t_hmaxNearP/D");
651  tree->Branch("t_selectTk", &t_selectTk, "t_selectTk/O");
652  tree->Branch("t_qltyFlag", &t_qltyFlag, "t_qltyFlag/O");
653  tree->Branch("t_qltyMissFlag",&t_qltyMissFlag,"t_qltyMissFlag/O");
654  tree->Branch("t_qltyPVFlag", &t_qltyPVFlag, "t_qltyPVFlag/O)");
655 
656  t_DetIds = new std::vector<unsigned int>();
657  t_HitEnergies = new std::vector<double>();
658  t_l1bits = new std::vector<bool>();
659  tree->Branch("t_DetIds", "std::vector<unsigned int>", &t_DetIds);
660  tree->Branch("t_HitEnergies", "std::vector<double>", &t_HitEnergies);
661  tree->Branch("t_l1bits","std::vector<bool>", &t_l1bits);
662 }
663 
664 // ------------ method called when starting to processes a run ------------
665 void IsoTrackCalib::beginRun(edm::Run const& iRun, edm::EventSetup const& iSetup) {
666  bool changed = false;
667  bool ok = hltConfig_.init(iRun,iSetup,"HLT",changed);
668  edm::LogInfo("IsoTrack") << "Run " << iRun.run() << " hltconfig.init " << ok;
669 
670  int iErrorCode = -1;
672  const AlgorithmMap& algorithmMap = m_l1GtMenu->gtAlgorithmMap();
673  const std::string& menuName = m_l1GtMenu->gtTriggerMenuName();
674 
675  if (verbosity_%10>0) edm::LogInfo("IsoTrack") << "menuName " << menuName;
676  for (CItAlgo itAlgo = algorithmMap.begin(); itAlgo != algorithmMap.end();
677  itAlgo++) {
678  std::string algName = itAlgo->first;
679  int algBitNumber = ( itAlgo->second ).algoBitNumber();
680  l1AlgoMap_.insert( std::pair<std::pair<unsigned int,std::string>,int>( std::pair<unsigned int,std::string>(algBitNumber, algName) , 0) ) ;
681  }
682  std::map< std::pair<unsigned int,std::string>, int>::iterator itr;
683  for (itr=l1AlgoMap_.begin(); itr!=l1AlgoMap_.end(); itr++) {
684  if (verbosity_%10>0)
685  edm::LogInfo("IsoTrack") << " ********** " << (itr->first).first << " "
686  << (itr->first).second <<" "<<itr->second;
687  }
688 }
689 
690 // ------------ method called when ending the processing of a run ------------
691 void IsoTrackCalib::endRun(edm::Run const& iRun, edm::EventSetup const&) {
692  edm::LogInfo("IsoTrack") << "endRun " << iRun.run() << std::endl;
693 }
694 
696  return (vec1.eta()-vec2.eta());
697 }
698 
700 
701  double phi1 = vec1.phi();
702  if (phi1 < 0) phi1 += 2.0*M_PI;
703  double phi2 = vec2.phi();
704  if (phi2 < 0) phi2 += 2.0*M_PI;
705  double dphi = phi1-phi2;
706  if (dphi > M_PI) dphi -= 2.*M_PI;
707  else if (dphi < -M_PI) dphi += 2.*M_PI;
708  return dphi;
709 }
710 
712  double deta = dEta(vec1,vec2);
713  double dphi = dPhi(vec1,vec2);
714  return std::sqrt(deta*deta + dphi*dphi);
715 }
716 
717 double IsoTrackCalib::deltaR(double eta1,double eta2,double phi1,double phi2){
718  double deta = eta1-eta2;
719  if (phi1 < 0) phi1 += 2.0*M_PI;
720  if (phi2 < 0) phi2 += 2.0*M_PI;
721  double dphi = phi1-phi2;
722  if (dphi > M_PI) dphi -= 2.*M_PI;
723  else if (dphi < -M_PI) dphi += 2.*M_PI;
724  return std::sqrt(deta*deta + dphi*dphi);
725 }
726 
727 
728 
729 //define this as a plug-in
731 
RunNumber_t run() const
Definition: EventID.h:39
static const std::string kSharedResource
Definition: TFileService.h:76
double p() const
momentum vector magnitude
Definition: TrackBase.h:648
void beginRun(edm::Run const &, edm::EventSetup const &) override
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 useRaw=0, bool debug=false)
type
Definition: HCALResponse.h:21
const std::string & gtTriggerMenuName() const
T getParameter(std::string const &) const
EventNumber_t event() const
Definition: EventID.h:41
spr::trackSelectionParameters selectionParameters_
std::vector< double > * t_HitEnergies
T getUntrackedParameter(std::string const &, T const &) const
const edm::EDGetTokenT< reco::TrackCollection > tok_genTrack_
void beginJob() override
std::vector< spr::propagatedTrackID > propagateCALO(edm::Handle< reco::TrackCollection > &trkCollection, const CaloGeometry *geo, const MagneticField *bField, const std::string &theTrackQuality, bool debug=false)
TH1I * h_tkEta0[5]
const unsigned int nTracks(const reco::Vertex &sv)
const L1GtTriggerMenu * m_l1GtMenu
const edm::EDGetTokenT< GenEventInfoProduct > tok_ew_
std::vector< unsigned int > * t_DetIds
double eta() const final
momentum pseudorapidity
const double a_mipR_
edm::EDGetTokenT< EcalRecHitCollection > tok_EE_
RunNumber_t run() const
Definition: RunBase.h:40
ParameterDescriptionBase * addUntracked(U const &iLabel, T const &value)
TProfile * h_RecHit_num
HLTConfigProvider hltConfig_
edm::EDGetTokenT< l1extra::L1JetParticleCollection > tok_L1extCenJet_
const int verbosity_
std::map< std::pair< unsigned int, std::string >, int > l1AlgoMap_
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:517
TrackQuality
track quality
Definition: TrackBase.h:151
void endJob() override
std::map< std::string, L1GtAlgorithm > AlgorithmMap
map containing the algorithms
const bool availableL1Configuration(int &errorCode, int &l1ConfCode) const
Definition: L1GtUtils.cc:1948
const edm::EDGetTokenT< reco::BeamSpot > tok_bs_
std::vector< GenJet > GenJetCollection
collection of GenJet objects
TH1I * h_tketa0[6]
std::vector< Track > TrackCollection
collection of Tracks
Definition: TrackFwd.h:15
std::vector< T >::const_iterator const_iterator
int bunchCrossing() const
Definition: EventBase.h:64
TH1I * h_tketa3[6]
edm::LuminosityBlockNumber_t luminosityBlock() const
Definition: EventBase.h:61
double phi() const
azimuthal angle of momentum vector
Definition: TrackBase.h:678
double pt() const final
transverse momentum
const bool decision(const edm::Event &iEvent, const std::string &nameAlgoTechTrig, int &errorCode) const
Definition: L1GtUtils.cc:1208
std::vector< Vertex > VertexCollection
collection of Vertex objects
Definition: VertexFwd.h:9
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
double chargeIsolationCone(unsigned int trkIndex, std::vector< spr::propagatedTrackDirection > &trkDirs, double dR, int &nNearTRKs, bool debug=false)
double px() const
x coordinate of momentum vector
Definition: TrackBase.h:660
double weight() const
edm::EDGetTokenT< l1extra::L1JetParticleCollection > tok_L1extTauJet_
edm::EDGetTokenT< HBHERecHitCollection > tok_hbhe_
const std::string theTrackQuality_
TH1I * h_tketa5[6]
TH1I * h_tketa1[6]
EDGetTokenT< ProductType > consumes(edm::InputTag const &tag)
AlgorithmMap::const_iterator CItAlgo
iterators through map containing the algorithms
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
TH1I * h_tketa4[6]
int iEvent
Definition: GenABIO.cc:224
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
double eta() const
pseudorapidity of momentum vector
Definition: TrackBase.h:684
std::vector< double > pbin
TH1I * h_tkEta2[5]
TH1I * h_tkEta4[5]
static const bool useL1GtTriggerMenuLite(true)
edm::Service< TFileService > fs_
double deltaR(double eta1, double eta2, double phi1, double phi2)
bool goodTrack(const reco::Track *pTrack, math::XYZPoint leadPV, trackSelectionParameters parameters, bool debug=false)
std::vector< bool > * t_l1bits
std::vector< double > vec1
Definition: HCALResponse.h:15
ConsumesCollector consumesCollector()
Use a ConsumesCollector to gather consumes information from helper functions.
T sqrt(T t)
Definition: SSEVec.h:18
double pt() const
track transverse momentum
Definition: TrackBase.h:654
TH1I * h_tketa2[6]
int ieta() const
get the cell ieta
Definition: HcalDetId.h:159
edm::EDGetTokenT< EcalRecHitCollection > tok_EB_
const double a_charIsoR_
Jets made from MC generator particles.
Definition: GenJet.h:25
const edm::EDGetTokenT< reco::VertexCollection > tok_recVtx_
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)
~IsoTrackCalib() override
edm::EDGetTokenT< l1extra::L1JetParticleCollection > tok_L1extFwdJet_
const edm::EDGetTokenT< reco::GenJetCollection > tok_jets_
bool isValid() const
Definition: HandleBase.h:74
#define M_PI
double pz() const
z coordinate of momentum vector
Definition: TrackBase.h:672
int k[5][pyjets_maxn]
const_iterator end() const
double t_EventWeight
double dEta(math::XYZTLorentzVector &, math::XYZTLorentzVector &)
static TrackQuality qualityByName(const std::string &name)
Definition: TrackBase.cc:134
L1GtUtils m_l1GtUtils
TH1I * h_tkEta3[5]
TProfile * h_RecHit_iEta
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
bool init(const edm::Run &iRun, const edm::EventSetup &iSetup, const std::string &processName, bool &changed)
d&#39;tor
void add(std::string const &label, ParameterSetDescription const &psetDescription)
TH1I * h_tkEta5[5]
TH1I * h_tkEta1[5]
edm::EventID id() const
Definition: EventBase.h:59
std::vector< PFJet > PFJetCollection
collection of PFJet objects
fixed size matrix
HLT enums.
void getL1GtRunCache(const edm::Run &, const edm::EventSetup &, const bool, const bool)
get all the run-constant quantities for L1 trigger and cache them
Definition: L1GtUtils.cc:319
T get() const
Definition: EventSetup.h:71
reco::TrackBase::TrackQuality minQuality
void endRun(edm::Run const &, edm::EventSetup const &) override
const AlgorithmMap & gtAlgorithmMap() const
get / set the algorithm map (by name)
void analyze(const edm::Event &, const edm::EventSetup &) override
const Point & position() const
position
Definition: BeamSpot.h:62
double dR(math::XYZTLorentzVector &, math::XYZTLorentzVector &)
double dPhi(math::XYZTLorentzVector &, math::XYZTLorentzVector &)
Definition: tree.py:1
const edm::EDGetTokenT< reco::PFJetCollection > tok_pfjets_
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
const L1GtTriggerMenu * ptrL1TriggerMenuEventSetup(int &errorCode)
return a pointer to the L1 trigger menu from event setup
Definition: L1GtUtils.cc:1883
std::vector< vec1 > vec2
Definition: HCALResponse.h:16
const double a_coneR_
TH1F * h_jetpt[4]
T const * product() const
Definition: ESHandle.h:86
const std::vector< std::string > l1Names_
const_iterator begin() const
double py() const
y coordinate of momentum vector
Definition: TrackBase.h:666
static const bool useL1EventSetup(true)
Definition: Run.h:45
IsoTrackCalib(const edm::ParameterSet &)