CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
GeneralPurposeTrackAnalyzer.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: Alignment/OfflineValidation
4 // Class: GeneralPurposeTrackAnalyzer
5 //
9 //
10 // Original Author: Marco Musich
11 // Created: Mon, 13 Jun 2016 15:07:11 GMT
12 //
13 //
14 
15 // ROOT includes files
16 
17 #include "TMath.h"
18 #include "TFile.h"
19 #include "TH1D.h"
20 #include "TH1I.h"
21 #include "TH2D.h"
22 #include "TProfile.h"
23 
24 // system includes files
25 
26 #include <cmath>
27 #include <fstream>
28 #include <iostream>
29 #include <map>
30 #include <string>
31 #include <vector>
32 #include <boost/range/adaptor/indexed.hpp>
33 
34 // user include files
81 
82 // toggle to enable debugging
83 #define DEBUG 0
84 
87 
88 class GeneralPurposeTrackAnalyzer : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
89 public:
92  magFieldToken_(esConsumes<edm::Transition::BeginRun>()),
93  latencyToken_(esConsumes<edm::Transition::BeginRun>()),
94  geomTokenBR_(esConsumes<edm::Transition::BeginRun>()),
97  usesResource(TFileService::kSharedResource);
98 
99  TkTag_ = pset.getParameter<edm::InputTag>("TkTag");
100  theTrackCollectionToken_ = consumes<reco::TrackCollection>(TkTag_);
101 
102  TriggerResultsTag_ = pset.getParameter<edm::InputTag>("TriggerResultsTag");
103  hltresultsToken_ = consumes<edm::TriggerResults>(TriggerResultsTag_);
104 
105  BeamSpotTag_ = pset.getParameter<edm::InputTag>("BeamSpotTag");
106  beamspotToken_ = consumes<reco::BeamSpot>(BeamSpotTag_);
107 
108  VerticesTag_ = pset.getParameter<edm::InputTag>("VerticesTag");
109  vertexToken_ = consumes<reco::VertexCollection>(VerticesTag_);
110 
111  isCosmics_ = pset.getParameter<bool>("isCosmics");
112 
113  pmap = std::make_unique<TrackerMap>("Pixel");
114  pmap->onlyPixel(true);
115  pmap->setTitle("Pixel Hit entries");
116  pmap->setPalette(1);
117 
118  tmap = std::make_unique<TrackerMap>("Strip");
119  tmap->setTitle("Strip Hit entries");
120  tmap->setPalette(1);
121 
122  pixelmap = std::make_unique<Phase1PixelMaps>("COLZ0 L");
123  pixelmap->bookBarrelHistograms("entriesBarrel", "# hits", "# pixel hits");
124  pixelmap->bookForwardHistograms("entriesForward", "# hits", "# pixel hits");
125 
126  pixelrocsmap_ = std::make_unique<Phase1PixelROCMaps>("");
127  }
128 
129  ~GeneralPurposeTrackAnalyzer() override = default;
130 
132 
133  template <class OBJECT_TYPE>
134  int GetIndex(const std::vector<OBJECT_TYPE *> &vec, const TString &name) {
135  for (const auto &iter : vec | boost::adaptors::indexed(0)) {
136  if (iter.value() && iter.value()->GetName() == name) {
137  return iter.index();
138  }
139  }
140  edm::LogError("GeneralPurposeTrackAnalyzer") << "@SUB=GeneralPurposeTrackAnalyzer::GetIndex"
141  << " could not find " << name;
142  return -1;
143  }
144 
145 private:
146  // tokens for the event setup
150 
154 
156 
158 
160 
161  std::unique_ptr<TrackerMap> tmap;
162  std::unique_ptr<TrackerMap> pmap;
163 
164  std::unique_ptr<Phase1PixelMaps> pixelmap;
165  std::unique_ptr<Phase1PixelROCMaps> pixelrocsmap_;
166 
167  TH1D *hchi2ndof;
168  TH1D *hNtrk;
169  TH1D *hNtrkZoom;
170  TH1I *htrkQuality;
171  TH1I *htrkAlgo;
172  TH1I *htrkOriAlgo;
174  TH1D *hP;
175  TH1D *hPPlus;
176  TH1D *hPMinus;
177  TH1D *hPt;
178  TH1D *hMinPt;
179  TH1D *hPtPlus;
180  TH1D *hPtMinus;
181  TH1D *hHit;
182  TH1D *hHit2D;
183 
190 
193 
196 
197  TH1D *hHitPlus;
198  TH1D *hHitMinus;
199 
200  TH1D *hPhp;
201  TH1D *hPthp;
202  TH1D *hHithp;
203  TH1D *hEtahp;
204  TH1D *hPhihp;
205  TH1D *hchi2ndofhp;
206  TH1D *hchi2Probhp;
207 
208  TH1D *hCharge;
209  TH1D *hQoverP;
210  TH1D *hQoverPZoom;
211  TH1D *hEta;
212  TH1D *hEtaPlus;
213  TH1D *hEtaMinus;
214  TH1D *hPhi;
215  TH1D *hPhiBarrel;
220  TH1D *hPhiPlus;
221  TH1D *hPhiMinus;
222 
223  TH1D *hDeltaPhi;
224  TH1D *hDeltaEta;
225  TH1D *hDeltaR;
226 
227  TH1D *hvx;
228  TH1D *hvy;
229  TH1D *hvz;
230  TH1D *hd0;
231  TH1D *hdz;
232  TH1D *hdxy;
233 
234  TH2D *hd0PVvsphi;
235  TH2D *hd0PVvseta;
236  TH2D *hd0PVvspt;
237 
238  TH2D *hd0vsphi;
239  TH2D *hd0vseta;
240  TH2D *hd0vspt;
241 
242  TH1D *hnhpxb;
243  TH1D *hnhpxe;
244  TH1D *hnhTIB;
245  TH1D *hnhTID;
246  TH1D *hnhTOB;
247  TH1D *hnhTEC;
248 
250 
251  TProfile *pNBpixHitsVsVx;
252  TProfile *pNBpixHitsVsVy;
253  TProfile *pNBpixHitsVsVz;
254 
255  TH1D *hMultCand;
256 
257  TH1D *hdxyBS;
258  TH1D *hd0BS;
259  TH1D *hdzBS;
260  TH1D *hdxyPV;
261  TH1D *hd0PV;
262  TH1D *hdzPV;
263  TH1D *hrun;
264  TH1D *hlumi;
265 
266  std::vector<TH1 *> vTrackHistos_;
267  std::vector<TH1 *> vTrackProfiles_;
268  std::vector<TH1 *> vTrack2DHistos_;
269 
272 
273  TH1D *modeByRun_;
274  TH1D *fieldByRun_;
275 
276  int ievt;
277  int itrks;
278  int mode;
280  bool isPhase1_;
281  float etaMax_;
282 
284 
289 
291 
296 
297  std::map<std::string, std::pair<int, int> > triggerMap_;
298  std::map<int, std::pair<int, float> > conditionsMap_;
299  std::map<int, std::pair<int, int> > runInfoMap_;
300 
301  //*************************************************************
302  void analyze(const edm::Event &event, const edm::EventSetup &setup) override
303  //*************************************************************
304  {
305  ievt++;
306 
308 
309  // geometry setup
310  const TrackerGeometry *theGeometry = &setup.getData(geomToken_);
311 
312  // switch on the phase1
313  if ((theGeometry->isThere(GeomDetEnumerators::P1PXB)) || (theGeometry->isThere(GeomDetEnumerators::P1PXEC))) {
314  isPhase1_ = true;
315  } else {
316  isPhase1_ = false;
317  }
318 
319  const reco::TrackCollection tC = *(trackCollection.product());
320  itrks += tC.size();
321 
322  runInfoMap_[event.run()].first += 1;
323  runInfoMap_[event.run()].second += tC.size();
324 
325  if (DEBUG) {
326  edm::LogInfo("GeneralPurposeTrackAnalyzer") << "Reconstructed " << tC.size() << " tracks" << std::endl;
327  }
328  //int iCounter=0;
329 
330  edm::Handle<edm::TriggerResults> hltresults = event.getHandle(hltresultsToken_);
331  if (hltresults.isValid()) {
332  const edm::TriggerNames &triggerNames_ = event.triggerNames(*hltresults);
333  int ntrigs = hltresults->size();
334  //const vector<std::string> &triggernames = triggerNames_.triggerNames();
335 
336  for (int itrig = 0; itrig != ntrigs; ++itrig) {
337  const std::string &trigName = triggerNames_.triggerName(itrig);
338  bool accept = hltresults->accept(itrig);
339  if (accept == 1) {
340  if (DEBUG) {
341  edm::LogInfo("GeneralPurposeTrackAnalyzer")
342  << trigName << " " << accept << " ,track size: " << tC.size() << std::endl;
343  }
344  triggerMap_[trigName].first += 1;
345  triggerMap_[trigName].second += tC.size();
346  // triggerInfo.push_back(pair <std::string, int> (trigName, accept));
347  }
348  }
349  }
350 
351  hrun->Fill(event.run());
352  hlumi->Fill(event.luminosityBlock());
353 
354  int nHighPurityTracks = 0;
355 
356  for (auto track = tC.cbegin(); track != tC.cend(); track++) {
357  unsigned int nHit2D = 0;
358  std::bitset<16> rocsToMask;
359  for (auto iHit = track->recHitsBegin(); iHit != track->recHitsEnd(); ++iHit) {
360  if (this->isHit2D(**iHit)) {
361  ++nHit2D;
362  }
363  // rest the ROCs for the map
364  rocsToMask.reset();
365  const DetId &detId = (*iHit)->geographicalId();
366  const GeomDet *geomDet(theGeometry->idToDet(detId));
367 
368  const SiPixelRecHit *pixhit = dynamic_cast<const SiPixelRecHit *>(*iHit);
369 
370  if (pixhit) {
371  if (pixhit->isValid()) {
372  unsigned int subid = detId.subdetId();
373  int detid_db = detId.rawId();
374 
375  // get the cluster
376  auto clustp = pixhit->cluster();
377 
378  if (clustp.isNull())
379  continue;
380  auto const &cluster = *clustp;
381  int row = cluster.x() - 0.5, col = cluster.y() - 0.5;
382  int rocId = coord_.roc(detId, std::make_pair(row, col));
383 
384  rocsToMask.set(rocId);
385  pixelrocsmap_->fillSelectedRocs(detid_db, rocsToMask, 1);
386 
387  if (!isPhase1_) {
388  pmap->fill(detid_db, 1);
389  } else {
390  if (subid == PixelSubdetector::PixelBarrel) {
391  pixelmap->fillBarrelBin("entriesBarrel", detid_db, 1);
392  } else {
393  pixelmap->fillForwardBin("entriesForward", detid_db, 1);
394  }
395  }
396 
397  LocalPoint lp = (*iHit)->localPosition();
398  //LocalError le = (*iHit)->localPositionError();
399 
400  GlobalPoint GP = geomDet->surface().toGlobal(lp);
401 
402  if ((subid == PixelSubdetector::PixelBarrel) || (subid == PixelSubdetector::PixelEndcap)) {
403  // 1 = PXB, 2 = PXF
404  if (subid == PixelSubdetector::PixelBarrel) {
405  hHitCountVsThetaBPix->Fill(GP.theta());
406  hHitCountVsPhiBPix->Fill(GP.phi());
407 
408  hHitCountVsZBPix->Fill(GP.z());
409  hHitCountVsXBPix->Fill(GP.x());
410  hHitCountVsYBPix->Fill(GP.y());
411 
412  } else if (subid == PixelSubdetector::PixelEndcap) {
413  hHitCountVsThetaFPix->Fill(GP.theta());
414  hHitCountVsPhiFPix->Fill(GP.phi());
415 
416  hHitCountVsZFPix->Fill(GP.z());
417  hHitCountVsXFPix->Fill(GP.x());
418  hHitCountVsYFPix->Fill(GP.y());
419  }
420  }
421  }
422  } else {
423  if ((*iHit)->isValid()) {
424  tmap->fill(detId.rawId(), 1);
425  }
426  }
427  }
428 
429  hHit2D->Fill(nHit2D);
430  hHit->Fill(track->numberOfValidHits());
431  hnhpxb->Fill(track->hitPattern().numberOfValidPixelBarrelHits());
432  hnhpxe->Fill(track->hitPattern().numberOfValidPixelEndcapHits());
433  hnhTIB->Fill(track->hitPattern().numberOfValidStripTIBHits());
434  hnhTID->Fill(track->hitPattern().numberOfValidStripTIDHits());
435  hnhTOB->Fill(track->hitPattern().numberOfValidStripTOBHits());
436  hnhTEC->Fill(track->hitPattern().numberOfValidStripTECHits());
437 
438  // fill hit composition histogram
439  if (track->hitPattern().numberOfValidPixelBarrelHits() != 0) {
440  hHitComposition->Fill(0., track->hitPattern().numberOfValidPixelBarrelHits());
441 
442  pNBpixHitsVsVx->Fill(track->vx(), track->hitPattern().numberOfValidPixelBarrelHits());
443  pNBpixHitsVsVy->Fill(track->vy(), track->hitPattern().numberOfValidPixelBarrelHits());
444  pNBpixHitsVsVz->Fill(track->vz(), track->hitPattern().numberOfValidPixelBarrelHits());
445  }
446  if (track->hitPattern().numberOfValidPixelEndcapHits() != 0) {
447  hHitComposition->Fill(1., track->hitPattern().numberOfValidPixelEndcapHits());
448  }
449  if (track->hitPattern().numberOfValidStripTIBHits() != 0) {
450  hHitComposition->Fill(2., track->hitPattern().numberOfValidStripTIBHits());
451  }
452  if (track->hitPattern().numberOfValidStripTIDHits() != 0) {
453  hHitComposition->Fill(3., track->hitPattern().numberOfValidStripTIDHits());
454  }
455  if (track->hitPattern().numberOfValidStripTOBHits() != 0) {
456  hHitComposition->Fill(4., track->hitPattern().numberOfValidStripTOBHits());
457  }
458  if (track->hitPattern().numberOfValidStripTECHits() != 0) {
459  hHitComposition->Fill(5., track->hitPattern().numberOfValidStripTECHits());
460  }
461 
462  hCharge->Fill(track->charge());
463  hQoverP->Fill(track->qoverp());
464  hQoverPZoom->Fill(track->qoverp());
465  hPt->Fill(track->pt());
466  hP->Fill(track->p());
467  hchi2ndof->Fill(track->normalizedChi2());
468  hEta->Fill(track->eta());
469  hPhi->Fill(track->phi());
470 
471  if (fabs(track->eta()) < 0.8) {
472  hPhiBarrel->Fill(track->phi());
473  }
474  if (track->eta() > 0.8 && track->eta() < 1.4) {
475  hPhiOverlapPlus->Fill(track->phi());
476  }
477  if (track->eta() < -0.8 && track->eta() > -1.4) {
478  hPhiOverlapMinus->Fill(track->phi());
479  }
480  if (track->eta() > 1.4) {
481  hPhiEndcapPlus->Fill(track->phi());
482  }
483  if (track->eta() < -1.4) {
484  hPhiEndcapMinus->Fill(track->phi());
485  }
486 
487  hd0->Fill(track->d0());
488  hdz->Fill(track->dz());
489  hdxy->Fill(track->dxy());
490  hvx->Fill(track->vx());
491  hvy->Fill(track->vy());
492  hvz->Fill(track->vz());
493 
494  htrkAlgo->Fill(static_cast<int>(track->algo()));
495  htrkOriAlgo->Fill(static_cast<int>(track->originalAlgo()));
496 
497  int myquality = -99;
498  if (track->quality(reco::TrackBase::undefQuality)) {
499  myquality = -1;
500  htrkQuality->Fill(myquality);
501  }
502  if (track->quality(reco::TrackBase::loose)) {
503  myquality = 0;
504  htrkQuality->Fill(myquality);
505  }
506  if (track->quality(reco::TrackBase::tight)) {
507  myquality = 1;
508  htrkQuality->Fill(myquality);
509  }
510  if (track->quality(reco::TrackBase::highPurity) && (!isCosmics_)) {
511  myquality = 2;
512  htrkQuality->Fill(myquality);
513  hPhp->Fill(track->p());
514  hPthp->Fill(track->pt());
515  hHithp->Fill(track->numberOfValidHits());
516  hEtahp->Fill(track->eta());
517  hPhihp->Fill(track->phi());
518  hchi2ndofhp->Fill(track->normalizedChi2());
519  hchi2Probhp->Fill(TMath::Prob(track->chi2(), track->ndof()));
520  nHighPurityTracks++;
521  }
522  if (track->quality(reco::TrackBase::confirmed)) {
523  myquality = 3;
524  htrkQuality->Fill(myquality);
525  }
526  if (track->quality(reco::TrackBase::goodIterative)) {
527  myquality = 4;
528  htrkQuality->Fill(myquality);
529  }
530 
531  // Fill 1D track histos
532  static const int etaindex = this->GetIndex(vTrackHistos_, "h_tracketa");
533  vTrackHistos_[etaindex]->Fill(track->eta());
534  static const int phiindex = this->GetIndex(vTrackHistos_, "h_trackphi");
535  vTrackHistos_[phiindex]->Fill(track->phi());
536  static const int numOfValidHitsindex = this->GetIndex(vTrackHistos_, "h_trackNumberOfValidHits");
537  vTrackHistos_[numOfValidHitsindex]->Fill(track->numberOfValidHits());
538  static const int numOfLostHitsindex = this->GetIndex(vTrackHistos_, "h_trackNumberOfLostHits");
539  vTrackHistos_[numOfLostHitsindex]->Fill(track->numberOfLostHits());
540 
541  GlobalPoint gPoint(track->vx(), track->vy(), track->vz());
542  double theLocalMagFieldInInverseGeV = magneticField_->inInverseGeV(gPoint).z();
543  double kappa = -track->charge() * theLocalMagFieldInInverseGeV / track->pt();
544 
545  static const int kappaindex = this->GetIndex(vTrackHistos_, "h_curvature");
546  vTrackHistos_[kappaindex]->Fill(kappa);
547  static const int kappaposindex = this->GetIndex(vTrackHistos_, "h_curvature_pos");
548  if (track->charge() > 0)
549  vTrackHistos_[kappaposindex]->Fill(fabs(kappa));
550  static const int kappanegindex = this->GetIndex(vTrackHistos_, "h_curvature_neg");
551  if (track->charge() < 0)
552  vTrackHistos_[kappanegindex]->Fill(fabs(kappa));
553 
554  double chi2Prob = TMath::Prob(track->chi2(), track->ndof());
555  double normchi2 = track->normalizedChi2();
556 
557  static const int normchi2index = this->GetIndex(vTrackHistos_, "h_normchi2");
558  vTrackHistos_[normchi2index]->Fill(normchi2);
559  static const int chi2index = this->GetIndex(vTrackHistos_, "h_chi2");
560  vTrackHistos_[chi2index]->Fill(track->chi2());
561  static const int chi2Probindex = this->GetIndex(vTrackHistos_, "h_chi2Prob");
562  vTrackHistos_[chi2Probindex]->Fill(chi2Prob);
563  static const int ptindex = this->GetIndex(vTrackHistos_, "h_pt");
564  static const int pt2index = this->GetIndex(vTrackHistos_, "h_ptrebin");
565  vTrackHistos_[ptindex]->Fill(track->pt());
566  vTrackHistos_[pt2index]->Fill(track->pt());
567  if (track->ptError() != 0.) {
568  static const int ptResolutionindex = this->GetIndex(vTrackHistos_, "h_ptResolution");
569  vTrackHistos_[ptResolutionindex]->Fill(track->ptError() / track->pt());
570  }
571  // Fill track profiles
572  static const int d0phiindex = this->GetIndex(vTrackProfiles_, "p_d0_vs_phi");
573  vTrackProfiles_[d0phiindex]->Fill(track->phi(), track->d0());
574  static const int dzphiindex = this->GetIndex(vTrackProfiles_, "p_dz_vs_phi");
575  vTrackProfiles_[dzphiindex]->Fill(track->phi(), track->dz());
576  static const int d0etaindex = this->GetIndex(vTrackProfiles_, "p_d0_vs_eta");
577  vTrackProfiles_[d0etaindex]->Fill(track->eta(), track->d0());
578  static const int dzetaindex = this->GetIndex(vTrackProfiles_, "p_dz_vs_eta");
579  vTrackProfiles_[dzetaindex]->Fill(track->eta(), track->dz());
580  static const int chiProbphiindex = this->GetIndex(vTrackProfiles_, "p_chi2Prob_vs_phi");
581  vTrackProfiles_[chiProbphiindex]->Fill(track->phi(), chi2Prob);
582  static const int chiProbabsd0index = this->GetIndex(vTrackProfiles_, "p_chi2Prob_vs_d0");
583  vTrackProfiles_[chiProbabsd0index]->Fill(fabs(track->d0()), chi2Prob);
584  static const int chiProbabsdzindex = this->GetIndex(vTrackProfiles_, "p_chi2Prob_vs_dz");
585  vTrackProfiles_[chiProbabsdzindex]->Fill(track->dz(), chi2Prob);
586  static const int chiphiindex = this->GetIndex(vTrackProfiles_, "p_chi2_vs_phi");
587  vTrackProfiles_[chiphiindex]->Fill(track->phi(), track->chi2());
588  static const int normchiphiindex = this->GetIndex(vTrackProfiles_, "p_normchi2_vs_phi");
589  vTrackProfiles_[normchiphiindex]->Fill(track->phi(), normchi2);
590  static const int chietaindex = this->GetIndex(vTrackProfiles_, "p_chi2_vs_eta");
591  vTrackProfiles_[chietaindex]->Fill(track->eta(), track->chi2());
592  static const int normchiptindex = this->GetIndex(vTrackProfiles_, "p_normchi2_vs_pt");
593  vTrackProfiles_[normchiptindex]->Fill(track->pt(), normchi2);
594  static const int normchipindex = this->GetIndex(vTrackProfiles_, "p_normchi2_vs_p");
595  vTrackProfiles_[normchipindex]->Fill(track->p(), normchi2);
596  static const int chiProbetaindex = this->GetIndex(vTrackProfiles_, "p_chi2Prob_vs_eta");
597  vTrackProfiles_[chiProbetaindex]->Fill(track->eta(), chi2Prob);
598  static const int normchietaindex = this->GetIndex(vTrackProfiles_, "p_normchi2_vs_eta");
599  vTrackProfiles_[normchietaindex]->Fill(track->eta(), normchi2);
600  static const int kappaphiindex = this->GetIndex(vTrackProfiles_, "p_kappa_vs_phi");
601  vTrackProfiles_[kappaphiindex]->Fill(track->phi(), kappa);
602  static const int kappaetaindex = this->GetIndex(vTrackProfiles_, "p_kappa_vs_eta");
603  vTrackProfiles_[kappaetaindex]->Fill(track->eta(), kappa);
604  static const int ptResphiindex = this->GetIndex(vTrackProfiles_, "p_ptResolution_vs_phi");
605  vTrackProfiles_[ptResphiindex]->Fill(track->phi(), track->ptError() / track->pt());
606  static const int ptResetaindex = this->GetIndex(vTrackProfiles_, "p_ptResolution_vs_eta");
607  vTrackProfiles_[ptResetaindex]->Fill(track->eta(), track->ptError() / track->pt());
608 
609  // Fill 2D track histos
610  static const int d0phiindex_2d = this->GetIndex(vTrack2DHistos_, "h2_d0_vs_phi");
611  vTrack2DHistos_[d0phiindex_2d]->Fill(track->phi(), track->d0());
612  static const int dzphiindex_2d = this->GetIndex(vTrack2DHistos_, "h2_dz_vs_phi");
613  vTrack2DHistos_[dzphiindex_2d]->Fill(track->phi(), track->dz());
614  static const int d0etaindex_2d = this->GetIndex(vTrack2DHistos_, "h2_d0_vs_eta");
615  vTrack2DHistos_[d0etaindex_2d]->Fill(track->eta(), track->d0());
616  static const int dzetaindex_2d = this->GetIndex(vTrack2DHistos_, "h2_dz_vs_eta");
617  vTrack2DHistos_[dzetaindex_2d]->Fill(track->eta(), track->dz());
618  static const int chiphiindex_2d = this->GetIndex(vTrack2DHistos_, "h2_chi2_vs_phi");
619  vTrack2DHistos_[chiphiindex_2d]->Fill(track->phi(), track->chi2());
620  static const int chiProbphiindex_2d = this->GetIndex(vTrack2DHistos_, "h2_chi2Prob_vs_phi");
621  vTrack2DHistos_[chiProbphiindex_2d]->Fill(track->phi(), chi2Prob);
622  static const int chiProbabsd0index_2d = this->GetIndex(vTrack2DHistos_, "h2_chi2Prob_vs_d0");
623  vTrack2DHistos_[chiProbabsd0index_2d]->Fill(fabs(track->d0()), chi2Prob);
624  static const int normchiphiindex_2d = this->GetIndex(vTrack2DHistos_, "h2_normchi2_vs_phi");
625  vTrack2DHistos_[normchiphiindex_2d]->Fill(track->phi(), normchi2);
626  static const int chietaindex_2d = this->GetIndex(vTrack2DHistos_, "h2_chi2_vs_eta");
627  vTrack2DHistos_[chietaindex_2d]->Fill(track->eta(), track->chi2());
628  static const int chiProbetaindex_2d = this->GetIndex(vTrack2DHistos_, "h2_chi2Prob_vs_eta");
629  vTrack2DHistos_[chiProbetaindex_2d]->Fill(track->eta(), chi2Prob);
630  static const int normchietaindex_2d = this->GetIndex(vTrack2DHistos_, "h2_normchi2_vs_eta");
631  vTrack2DHistos_[normchietaindex_2d]->Fill(track->eta(), normchi2);
632  static const int kappaphiindex_2d = this->GetIndex(vTrack2DHistos_, "h2_kappa_vs_phi");
633  vTrack2DHistos_[kappaphiindex_2d]->Fill(track->phi(), kappa);
634  static const int kappaetaindex_2d = this->GetIndex(vTrack2DHistos_, "h2_kappa_vs_eta");
635  vTrack2DHistos_[kappaetaindex_2d]->Fill(track->eta(), kappa);
636  static const int normchi2kappa_2d = this->GetIndex(vTrack2DHistos_, "h2_normchi2_vs_kappa");
637  vTrack2DHistos_[normchi2kappa_2d]->Fill(normchi2, kappa);
638 
639  //dxy with respect to the beamspot
641  edm::Handle<reco::BeamSpot> beamSpotHandle = event.getHandle(beamspotToken_);
642  if (beamSpotHandle.isValid()) {
643  beamSpot = *beamSpotHandle;
644  math::XYZPoint point(beamSpot.x0(), beamSpot.y0(), beamSpot.z0());
645  double dxy = track->dxy(point);
646  double dz = track->dz(point);
647  hdxyBS->Fill(dxy);
648  hd0BS->Fill(-dxy);
649  hdzBS->Fill(dz);
650  }
651 
652  //dxy with respect to the primary vertex
653  reco::Vertex pvtx;
654  edm::Handle<reco::VertexCollection> vertexHandle = event.getHandle(vertexToken_);
655  double mindxy = 100.;
656  double dz = 100;
657  if (vertexHandle.isValid() && !isCosmics_) {
658  for (auto pvtx = vertexHandle->cbegin(); pvtx != vertexHandle->cend(); ++pvtx) {
659  math::XYZPoint mypoint(pvtx->x(), pvtx->y(), pvtx->z());
660  if (abs(mindxy) > abs(track->dxy(mypoint))) {
661  mindxy = track->dxy(mypoint);
662  dz = track->dz(mypoint);
663  }
664  }
665 
666  hdxyPV->Fill(mindxy);
667  hd0PV->Fill(-mindxy);
668  hdzPV->Fill(dz);
669 
670  hd0PVvsphi->Fill(track->phi(), -mindxy);
671  hd0PVvseta->Fill(track->eta(), -mindxy);
672  hd0PVvspt->Fill(track->pt(), -mindxy);
673 
674  } else {
675  hdxyPV->Fill(100);
676  hd0PV->Fill(100);
677  hdzPV->Fill(100);
678 
679  hd0vsphi->Fill(track->phi(), -track->dxy());
680  hd0vseta->Fill(track->eta(), -track->dxy());
681  hd0vspt->Fill(track->pt(), -track->dxy());
682  }
683 
684  if (DEBUG) {
685  edm::LogInfo("GeneralPurposeTrackAnalyzer") << "end of track loop" << std::endl;
686  }
687  }
688 
689  hNtrk->Fill(tC.size());
690  hNtrkZoom->Fill(tC.size());
691  hNhighPurity->Fill(nHighPurityTracks);
692 
693  if (DEBUG) {
694  edm::LogInfo("GeneralPurposeTrackAnalyzer") << "end of analysis" << std::endl;
695  }
696  }
697 
698  //*************************************************************
699  void beginRun(edm::Run const &run, edm::EventSetup const &setup) override
700  //*************************************************************
701  {
702  // Magnetic Field setup
703  magneticField_ = setup.getHandle(magFieldToken_);
704  float B_ = magneticField_.product()->inTesla(GlobalPoint(0, 0, 0)).mag();
705 
706  if (DEBUG) {
707  edm::LogInfo("GeneralPurposeTrackAnalyzer")
708  << "run number:" << run.run() << " magnetic field: " << B_ << " [T]" << std::endl;
709  }
710 
711  //SiStrip Latency
713  if (apvlat->singleReadOutMode() == 1) {
714  mode = 1; // peak mode
715  } else if (apvlat->singleReadOutMode() == 0) {
716  mode = -1; // deco mode
717  }
718 
719  conditionsMap_[run.run()].first = mode;
720  conditionsMap_[run.run()].second = B_;
721 
722  // init the sipixel coordinates
723  const TrackerGeometry *trackerGeometry = &setup.getData(geomTokenBR_);
724  const TrackerTopology *trackerTopology = &setup.getData(trackerTopologyTokenBR_);
725  const SiPixelFedCablingMap *siPixelFedCablingMap = &setup.getData(siPixelFedCablingMapTokenBR_);
726 
727  // Pixel Phase-1 helper class
728  coord_.init(trackerTopology, trackerGeometry, siPixelFedCablingMap);
729  }
730 
731  //*************************************************************
732  void endRun(edm::Run const &, edm::EventSetup const &) override {}
733  //*************************************************************
734 
735  //*************************************************************
736  void beginJob() override
737  //*************************************************************
738  {
739  if (DEBUG) {
740  edm::LogInfo("GeneralPurposeTrackAnalyzer") << __LINE__ << std::endl;
741  }
742 
743  TH1D::SetDefaultSumw2(kTRUE);
744 
745  isPhase1_ = true;
746  etaMax_ = 3.0;
747 
748  ievt = 0;
749  itrks = 0;
750 
751  hrun = fs->make<TH1D>("h_run", "run", 100000, 230000, 240000);
752  hlumi = fs->make<TH1D>("h_lumi", "lumi", 1000, 0, 1000);
753 
754  hchi2ndof = fs->make<TH1D>("h_chi2ndof", "chi2/ndf;#chi^{2}/ndf;tracks", 100, 0, 5.);
755  hCharge = fs->make<TH1D>("h_charge", "charge;Charge of the track;tracks", 5, -2.5, 2.5);
756  hNtrk = fs->make<TH1D>("h_Ntrk", "ntracks;Number of Tracks;events", 200, 0., 200.);
757  hNtrkZoom = fs->make<TH1D>("h_NtrkZoom", "Number of tracks; number of tracks;events", 10, 0., 10.);
758  hNhighPurity =
759  fs->make<TH1D>("h_NhighPurity", "n. high purity tracks;Number of high purity tracks;events", 200, 0., 200.);
760 
761  htrkAlgo = fs->make<TH1I>("h_trkAlgo",
762  "tracking step;iterative tracking step;tracks",
764  0.,
765  double(reco::TrackBase::algoSize));
766 
767  htrkOriAlgo = fs->make<TH1I>("h_trkOriAlgo",
768  "original tracking step;original iterative tracking step;tracks",
770  0.,
771  double(reco::TrackBase::algoSize));
772 
773  for (size_t ibin = 0; ibin < reco::TrackBase::algoSize - 1; ibin++) {
774  htrkAlgo->GetXaxis()->SetBinLabel(ibin + 1, (reco::TrackBase::algoNames[ibin]).c_str());
775  htrkOriAlgo->GetXaxis()->SetBinLabel(ibin + 1, (reco::TrackBase::algoNames[ibin]).c_str());
776  }
777 
778  htrkQuality = fs->make<TH1I>("h_trkQuality", "track quality;track quality;tracks", 6, -1, 5);
779  std::string qualities[7] = {"undef", "loose", "tight", "highPurity", "confirmed", "goodIterative"};
780  for (int nbin = 1; nbin <= htrkQuality->GetNbinsX(); nbin++) {
781  htrkQuality->GetXaxis()->SetBinLabel(nbin, (qualities[nbin - 1]).c_str());
782  }
783 
784  hP = fs->make<TH1D>("h_P", "Momentum;track momentum [GeV];tracks", 100, 0., 100.);
785  hQoverP = fs->make<TH1D>("h_qoverp", "Track q/p; track q/p [GeV^{-1}];tracks", 100, -1., 1.);
786  hQoverPZoom = fs->make<TH1D>("h_qoverpZoom", "Track q/p; track q/p [GeV^{-1}];tracks", 100, -0.1, 0.1);
787  hPt = fs->make<TH1D>("h_Pt", "Transverse Momentum;track p_{T} [GeV];tracks", 100, 0., 100.);
788  hHit = fs->make<TH1D>("h_nHits", "Number of hits;track n. hits;tracks", 50, -0.5, 49.5);
789  hHit2D = fs->make<TH1D>("h_nHit2D", "Number of 2D hits; number of 2D hits;tracks", 20, 0, 20);
790 
791  hHitCountVsZBPix = fs->make<TH1D>("h_HitCountVsZBpix", "Number of BPix hits vs z;hit global z;hits", 60, -30, 30);
793  fs->make<TH1D>("h_HitCountVsZFpix", "Number of FPix hits vs z;hit global z;hits", 100, -100, 100);
794 
795  hHitCountVsXBPix = fs->make<TH1D>("h_HitCountVsXBpix", "Number of BPix hits vs x;hit global x;hits", 20, -20, 20);
796  hHitCountVsXFPix = fs->make<TH1D>("h_HitCountVsXFpix", "Number of FPix hits vs x;hit global x;hits", 20, -20, 20);
797 
798  hHitCountVsYBPix = fs->make<TH1D>("h_HitCountVsYBpix", "Number of BPix hits vs y;hit global y;hits", 20, -20, 20);
799  hHitCountVsYFPix = fs->make<TH1D>("h_HitCountVsYFpix", "Number of FPix hits vs y;hit global y;hits", 20, -20, 20);
800 
802  fs->make<TH1D>("h_HitCountVsThetaBpix", "Number of BPix hits vs #theta;hit global #theta;hits", 20, 0., M_PI);
804  fs->make<TH1D>("h_HitCountVsPhiBpix", "Number of BPix hits vs #phi;hit global #phi;hits", 20, -M_PI, M_PI);
805 
807  fs->make<TH1D>("h_HitCountVsThetaFpix", "Number of FPix hits vs #theta;hit global #theta;hits", 40, 0., M_PI);
809  fs->make<TH1D>("h_HitCountVsPhiFpix", "Number of FPix hits vs #phi;hit global #phi;hits", 20, -M_PI, M_PI);
810 
811  hEta = fs->make<TH1D>("h_Eta", "Track pseudorapidity; track #eta;tracks", 100, -etaMax_, etaMax_);
812  hPhi = fs->make<TH1D>("h_Phi", "Track azimuth; track #phi;tracks", 100, -M_PI, M_PI);
813 
814  hPhiBarrel = fs->make<TH1D>("h_PhiBarrel", "hPhiBarrel (0<|#eta|<0.8);track #Phi;tracks", 100, -M_PI, M_PI);
816  fs->make<TH1D>("h_PhiOverlapPlus", "hPhiOverlapPlus (0.8<#eta<1.4);track #phi;tracks", 100, -M_PI, M_PI);
818  fs->make<TH1D>("h_PhiOverlapMinus", "hPhiOverlapMinus (-1.4<#eta<-0.8);track #phi;tracks", 100, -M_PI, M_PI);
819  hPhiEndcapPlus = fs->make<TH1D>("h_PhiEndcapPlus", "hPhiEndcapPlus (#eta>1.4);track #phi;track", 100, -M_PI, M_PI);
821  fs->make<TH1D>("h_PhiEndcapMinus", "hPhiEndcapMinus (#eta<1.4);track #phi;tracks", 100, -M_PI, M_PI);
822 
823  if (!isCosmics_) {
824  hPhp = fs->make<TH1D>("h_P_hp", "Momentum (high purity);track momentum [GeV];tracks", 100, 0., 100.);
825  hPthp = fs->make<TH1D>("h_Pt_hp", "Transverse Momentum (high purity);track p_{T} [GeV];tracks", 100, 0., 100.);
826  hHithp = fs->make<TH1D>("h_nHit_hp", "Number of hits (high purity);track n. hits;tracks", 30, 0, 30);
827  hEtahp =
828  fs->make<TH1D>("h_Eta_hp", "Track pseudorapidity (high purity); track #eta;tracks", 100, -etaMax_, etaMax_);
829  hPhihp = fs->make<TH1D>("h_Phi_hp", "Track azimuth (high purity); track #phi;tracks", 100, -M_PI, M_PI);
830  hchi2ndofhp = fs->make<TH1D>("h_chi2ndof_hp", "chi2/ndf (high purity);#chi^{2}/ndf;tracks", 100, 0, 5.);
831  hchi2Probhp = fs->make<TH1D>(
832  "h_chi2_Prob_hp", "#chi^{2} probability (high purity);#chi^{2}prob_{Track};Number of Tracks", 100, 0.0, 1.);
833 
834  hvx = fs->make<TH1D>("h_vx", "Track v_{x} ; track v_{x} [cm];tracks", 100, -1.5, 1.5);
835  hvy = fs->make<TH1D>("h_vy", "Track v_{y} ; track v_{y} [cm];tracks", 100, -1.5, 1.5);
836  hvz = fs->make<TH1D>("h_vz", "Track v_{z} ; track v_{z} [cm];tracks", 100, -20., 20.);
837  hd0 = fs->make<TH1D>("h_d0", "Track d_{0} ; track d_{0} [cm];tracks", 100, -1., 1.);
838  hdxy = fs->make<TH1D>("h_dxy", "Track d_{xy}; track d_{xy} [cm]; tracks", 100, -0.5, 0.5);
839  hdz = fs->make<TH1D>("h_dz", "Track d_{z} ; track d_{z} [cm]; tracks", 100, -20, 20);
840 
841  hd0PVvsphi =
842  fs->make<TH2D>("h2_d0PVvsphi", "hd0PVvsphi;track #phi;track d_{0}(PV) [cm]", 160, -M_PI, M_PI, 100, -1., 1.);
843  hd0PVvseta =
844  fs->make<TH2D>("h2_d0PVvseta", "hdPV0vseta;track #eta;track d_{0}(PV) [cm]", 160, -2.5, 2.5, 100, -1., 1.);
845  hd0PVvspt = fs->make<TH2D>("h2_d0PVvspt", "hdPV0vspt;track p_{T};d_{0}(PV) [cm]", 50, 0., 100., 100, -1, 1.);
846 
847  hdxyBS = fs->make<TH1D>("h_dxyBS", "hdxyBS; track d_{xy}(BS) [cm];tracks", 100, -0.1, 0.1);
848  hd0BS = fs->make<TH1D>("h_d0BS", "hd0BS ; track d_{0}(BS) [cm];tracks", 100, -0.1, 0.1);
849  hdzBS = fs->make<TH1D>("h_dzBS", "hdzBS ; track d_{z}(BS) [cm];tracks", 100, -12, 12);
850  hdxyPV = fs->make<TH1D>("h_dxyPV", "hdxyPV; track d_{xy}(PV) [cm];tracks", 100, -0.1, 0.1);
851  hd0PV = fs->make<TH1D>("h_d0PV", "hd0PV ; track d_{0}(PV) [cm];tracks", 100, -0.15, 0.15);
852  hdzPV = fs->make<TH1D>("h_dzPV", "hdzPV ; track d_{z}(PV) [cm];tracks", 100, -0.1, 0.1);
853 
854  hnhTIB = fs->make<TH1D>("h_nHitTIB", "nhTIB;# hits in TIB; tracks", 20, 0., 20.);
855  hnhTID = fs->make<TH1D>("h_nHitTID", "nhTID;# hits in TID; tracks", 20, 0., 20.);
856  hnhTOB = fs->make<TH1D>("h_nHitTOB", "nhTOB;# hits in TOB; tracks", 20, 0., 20.);
857  hnhTEC = fs->make<TH1D>("h_nHitTEC", "nhTEC;# hits in TEC; tracks", 20, 0., 20.);
858 
859  } else {
860  hvx = fs->make<TH1D>("h_vx", "Track v_{x};track v_{x} [cm];tracks", 100, -100., 100.);
861  hvy = fs->make<TH1D>("h_vy", "Track v_{y};track v_{y} [cm];tracks", 100, -100., 100.);
862  hvz = fs->make<TH1D>("h_vz", "Track v_{z};track v_{z} [cm];track", 100, -100., 100.);
863  hd0 = fs->make<TH1D>("h_d0", "Track d_{0};track d_{0} [cm];track", 100, -100., 100.);
864  hdxy = fs->make<TH1D>("h_dxy", "Track d_{xy};track d_{xy} [cm];tracks", 100, -100, 100);
865  hdz = fs->make<TH1D>("h_dz", "Track d_{z};track d_{z} [cm];tracks", 100, -200, 200);
866 
867  hd0vsphi = fs->make<TH2D>(
868  "h2_d0vsphi", "Track d_{0} vs #phi; track #phi;track d_{0} [cm]", 160, -3.20, 3.20, 100, -100., 100.);
869  hd0vseta = fs->make<TH2D>(
870  "h2_d0vseta", "Track d_{0} vs #eta; track #eta;track d_{0} [cm]", 160, -3.20, 3.20, 100, -100., 100.);
871  hd0vspt = fs->make<TH2D>(
872  "h2_d0vspt", "Track d_{0} vs p_{T};track p_{T};track d_{0} [cm]", 50, 0., 100., 100, -100, 100);
873 
874  hdxyBS = fs->make<TH1D>("h_dxyBS", "Track d_{xy}(BS);d_{xy}(BS) [cm];tracks", 100, -100., 100.);
875  hd0BS = fs->make<TH1D>("h_d0BS", "Track d_{0}(BS);d_{0}(BS) [cm];tracks", 100, -100., 100.);
876  hdzBS = fs->make<TH1D>("h_dzBS", "Track d_{z}(BS);d_{z}(BS) [cm];tracks", 100, -100., 100.);
877  hdxyPV = fs->make<TH1D>("h_dxyPV", "Track d_{xy}(PV); d_{xy}(PV) [cm];tracks", 100, -100., 100.);
878  hd0PV = fs->make<TH1D>("h_d0PV", "Track d_{0}(PV); d_{0}(PV) [cm];tracks", 100, -100., 100.);
879  hdzPV = fs->make<TH1D>("h_dzPV", "Track d_{z}(PV); d_{z}(PV) [cm];tracks", 100, -100., 100.);
880 
881  hnhTIB = fs->make<TH1D>("h_nHitTIB", "nhTIB;# hits in TIB; tracks", 30, 0., 30.);
882  hnhTID = fs->make<TH1D>("h_nHitTID", "nhTID;# hits in TID; tracks", 30, 0., 30.);
883  hnhTOB = fs->make<TH1D>("h_nHitTOB", "nhTOB;# hits in TOB; tracks", 30, 0., 30.);
884  hnhTEC = fs->make<TH1D>("h_nHitTEC", "nhTEC;# hits in TEC; tracks", 30, 0., 30.);
885  }
886 
887  hnhpxb = fs->make<TH1D>("h_nHitPXB", "nhpxb;# hits in Pixel Barrel; tracks", 10, 0., 10.);
888  hnhpxe = fs->make<TH1D>("h_nHitPXE", "nhpxe;# hits in Pixel Endcap; tracks", 10, 0., 10.);
889 
890  hHitComposition = fs->make<TH1D>("h_hitcomposition", "track hit composition;;# hits", 6, -0.5, 5.5);
891 
893  fs->make<TProfile>("p_NpixHits_vs_Vx", "n. Barrel Pixel hits vs. v_{x};v_{x} (cm);n. BPix hits", 20, -20, 20);
894 
896  fs->make<TProfile>("p_NpixHits_vs_Vy", "n. Barrel Pixel hits vs. v_{y};v_{y} (cm);n. BPix hits", 20, -20, 20);
897 
899  fs->make<TProfile>("p_NpixHits_vs_Vz", "n. Barrel Pixel hits vs. v_{z};v_{z} (cm);n. BPix hits", 20, -100, 100);
900 
901  std::string dets[6] = {"PXB", "PXF", "TIB", "TID", "TOB", "TEC"};
902 
903  for (int i = 1; i <= hHitComposition->GetNbinsX(); i++) {
904  hHitComposition->GetXaxis()->SetBinLabel(i, (dets[i - 1]).c_str());
905  }
906 
907  vTrackHistos_.push_back(
908  fs->make<TH1F>("h_tracketa", "Track #eta;#eta_{Track};Number of Tracks", 90, -etaMax_, etaMax_));
909  vTrackHistos_.push_back(fs->make<TH1F>("h_trackphi", "Track #phi;#phi_{Track};Number of Tracks", 90, -M_PI, M_PI));
910  vTrackHistos_.push_back(fs->make<TH1F>(
911  "h_trackNumberOfValidHits", "Track # of valid hits;# of valid hits _{Track};Number of Tracks", 40, 0., 40.));
912  vTrackHistos_.push_back(fs->make<TH1F>(
913  "h_trackNumberOfLostHits", "Track # of lost hits;# of lost hits _{Track};Number of Tracks", 10, 0., 10.));
914  vTrackHistos_.push_back(
915  fs->make<TH1F>("h_curvature", "Curvature #kappa;#kappa_{Track};Number of Tracks", 100, -.05, .05));
916  vTrackHistos_.push_back(fs->make<TH1F>(
917  "h_curvature_pos", "Curvature |#kappa| Positive Tracks;|#kappa_{pos Track}|;Number of Tracks", 100, .0, .05));
918  vTrackHistos_.push_back(fs->make<TH1F>(
919  "h_curvature_neg", "Curvature |#kappa| Negative Tracks;|#kappa_{neg Track}|;Number of Tracks", 100, .0, .05));
920  vTrackHistos_.push_back(
921  fs->make<TH1F>("h_diff_curvature",
922  "Curvature |#kappa| Tracks Difference;|#kappa_{Track}|;# Pos Tracks - # Neg Tracks",
923  100,
924  .0,
925  .05));
926 
927  vTrackHistos_.push_back(
928  fs->make<TH1F>("h_chi2", "Track #chi^{2};#chi^{2}_{Track};Number of Tracks", 500, -0.01, 500.));
929  vTrackHistos_.push_back(
930  fs->make<TH1F>("h_chi2Prob", "#chi^{2} probability;Track Prob(#chi^{2},ndof);Number of Tracks", 100, 0.0, 1.));
931  vTrackHistos_.push_back(
932  fs->make<TH1F>("h_normchi2", "#chi^{2}/ndof;#chi^{2}/ndof;Number of Tracks", 100, -0.01, 10.));
933  //variable binning for chi2/ndof vs. pT
934  double xBins[19] = {0., 0.15, 0.5, 1., 1.5, 2., 2.5, 3., 3.5, 4., 4.5, 5., 7., 10., 15., 25., 40., 100., 200.};
935  vTrackHistos_.push_back(fs->make<TH1F>("h_pt", "Track p_{T};p_{T}^{track} [GeV];Number of Tracks", 250, 0., 250));
936  vTrackHistos_.push_back(fs->make<TH1F>("h_ptrebin", "Track p_{T};p_{T}^{track} [GeV];Number of Tracks", 18, xBins));
937  vTrackHistos_.push_back(fs->make<TH1F>(
938  "h_ptResolution", "#delta_{p_{T}}/p_{T}^{track};#delta_{p_{T}}/p_{T}^{track};Number of Tracks", 100, 0., 0.5));
939  vTrackProfiles_.push_back(fs->make<TProfile>(
940  "p_d0_vs_phi", "Transverse Impact Parameter vs. #phi;#phi_{Track};#LT d_{0} #GT [cm]", 100, -M_PI, M_PI));
941  vTrackProfiles_.push_back(fs->make<TProfile>(
942  "p_dz_vs_phi", "Longitudinal Impact Parameter vs. #phi;#phi_{Track};#LT d_{z} #GT [cm]", 100, -M_PI, M_PI));
943  vTrackProfiles_.push_back(fs->make<TProfile>(
944  "p_d0_vs_eta", "Transverse Impact Parameter vs. #eta;#eta_{Track};#LT d_{0} #GT [cm]", 100, -etaMax_, etaMax_));
945  vTrackProfiles_.push_back(
946  fs->make<TProfile>("p_dz_vs_eta",
947  "Longitudinal Impact Parameter vs. #eta;#eta_{Track};#LT d_{z} #GT [cm]",
948  100,
949  -etaMax_,
950  etaMax_));
951  vTrackProfiles_.push_back(
952  fs->make<TProfile>("p_chi2_vs_phi", "#chi^{2} vs. #phi;#phi_{Track};#LT #chi^{2} #GT", 100, -M_PI, M_PI));
953  vTrackProfiles_.push_back(
954  fs->make<TProfile>("p_chi2Prob_vs_phi",
955  "#chi^{2} probablility vs. #phi;#phi_{Track};#LT #chi^{2} probability#GT",
956  100,
957  -M_PI,
958  M_PI));
959 
960  vTrackProfiles_.push_back(fs->make<TProfile>(
961  "p_chi2Prob_vs_d0", "#chi^{2} probablility vs. |d_{0}|;|d_{0}|[cm];#LT #chi^{2} probability#GT", 100, 0, 80));
962  vTrackProfiles_.push_back(fs->make<TProfile>(
963  "p_chi2Prob_vs_dz", "#chi^{2} probablility vs. dz;d_{z} [cm];#LT #chi^{2} probability#GT", 100, -30, 30));
964  vTrackProfiles_.push_back(fs->make<TProfile>(
965  "p_normchi2_vs_phi", "#chi^{2}/ndof vs. #phi;#phi_{Track};#LT #chi^{2}/ndof #GT", 100, -M_PI, M_PI));
966  vTrackProfiles_.push_back(
967  fs->make<TProfile>("p_chi2_vs_eta", "#chi^{2} vs. #eta;#eta_{Track};#LT #chi^{2} #GT", 100, -etaMax_, etaMax_));
968  vTrackProfiles_.push_back(fs->make<TProfile>(
969  "p_normchi2_vs_pt", "norm #chi^{2} vs. p_{T}_{Track}; p_{T}_{Track};#LT #chi^{2}/ndof #GT", 18, xBins));
970  vTrackProfiles_.push_back(fs->make<TProfile>(
971  "p_normchi2_vs_p", "#chi^{2}/ndof vs. p_{Track};p_{Track};#LT #chi^{2}/ndof #GT", 18, xBins));
972  vTrackProfiles_.push_back(
973  fs->make<TProfile>("p_chi2Prob_vs_eta",
974  "#chi^{2} probability vs. #eta;#eta_{Track};#LT #chi^{2} probability #GT",
975  100,
976  -etaMax_,
977  etaMax_));
978 
979  vTrackProfiles_.push_back(fs->make<TProfile>(
980  "p_normchi2_vs_eta", "#chi^{2}/ndof vs. #eta;#eta_{Track};#LT #chi^{2}/ndof #GT", 100, -etaMax_, etaMax_));
981  vTrackProfiles_.push_back(
982  fs->make<TProfile>("p_kappa_vs_phi", "#kappa vs. #phi;#phi_{Track};#kappa", 100, -M_PI, M_PI));
983  vTrackProfiles_.push_back(
984  fs->make<TProfile>("p_kappa_vs_eta", "#kappa vs. #eta;#eta_{Track};#kappa", 100, -etaMax_, etaMax_));
985  vTrackProfiles_.push_back(
986  fs->make<TProfile>("p_ptResolution_vs_phi",
987  "#delta_{p_{T}}/p_{T}^{track};#phi^{track};#delta_{p_{T}}/p_{T}^{track}",
988  100,
989  -M_PI,
990  M_PI));
991 
992  vTrackProfiles_.push_back(
993  fs->make<TProfile>("p_ptResolution_vs_eta",
994  "#delta_{p_{T}}/p_{T}^{track};#eta^{track};#delta_{p_{T}}/p_{T}^{track}",
995  100,
996  -etaMax_,
997  etaMax_));
998 
999  vTrack2DHistos_.push_back(fs->make<TH2F>(
1000  "h2_d0_vs_phi", "Transverse Impact Parameter vs. #phi;#phi_{Track};d_{0} [cm]", 100, -M_PI, M_PI, 100, -1., 1.));
1001  vTrack2DHistos_.push_back(fs->make<TH2F>("h2_dz_vs_phi",
1002  "Longitudinal Impact Parameter vs. #phi;#phi_{Track};d_{z} [cm]",
1003  100,
1004  -M_PI,
1005  M_PI,
1006  100,
1007  -100.,
1008  100.));
1009 
1010  vTrack2DHistos_.push_back(fs->make<TH2F>("h2_d0_vs_eta",
1011  "Transverse Impact Parameter vs. #eta;#eta_{Track};d_{0} [cm]",
1012  100,
1013  -etaMax_,
1014  etaMax_,
1015  100,
1016  -1.,
1017  1.));
1018  vTrack2DHistos_.push_back(fs->make<TH2F>("h2_dz_vs_eta",
1019  "Longitudinal Impact Parameter vs. #eta;#eta_{Track};d_{z} [cm]",
1020  100,
1021  -etaMax_,
1022  etaMax_,
1023  100,
1024  -100.,
1025  100.));
1026 
1027  vTrack2DHistos_.push_back(
1028  fs->make<TH2F>("h2_chi2_vs_phi", "#chi^{2} vs. #phi;#phi_{Track};#chi^{2}", 100, -M_PI, M_PI, 500, 0., 500.));
1029  vTrack2DHistos_.push_back(fs->make<TH2F>("h2_chi2Prob_vs_phi",
1030  "#chi^{2} probability vs. #phi;#phi_{Track};#chi^{2} probability",
1031  100,
1032  -M_PI,
1033  M_PI,
1034  100,
1035  0.,
1036  1.));
1037 
1038  vTrack2DHistos_.push_back(fs->make<TH2F>("h2_chi2Prob_vs_d0",
1039  "#chi^{2} probability vs. |d_{0}|;|d_{0}| [cm];#chi^{2} probability",
1040  100,
1041  0,
1042  80,
1043  100,
1044  0.,
1045  1.));
1046 
1047  vTrack2DHistos_.push_back(fs->make<TH2F>(
1048  "h2_normchi2_vs_phi", "#chi^{2}/ndof vs. #phi;#phi_{Track};#chi^{2}/ndof", 100, -M_PI, M_PI, 100, 0., 10.));
1049  vTrack2DHistos_.push_back(fs->make<TH2F>(
1050  "h2_chi2_vs_eta", "#chi^{2} vs. #eta;#eta_{Track};#chi^{2}", 100, -etaMax_, etaMax_, 500, 0., 500.));
1051  vTrack2DHistos_.push_back(fs->make<TH2F>("h2_chi2Prob_vs_eta",
1052  "#chi^{2} probaility vs. #eta;#eta_{Track};#chi^{2} probability",
1053  100,
1054  -etaMax_,
1055  etaMax_,
1056  100,
1057  0.,
1058  1.));
1059 
1060  vTrack2DHistos_.push_back(fs->make<TH2F>(
1061  "h2_normchi2_vs_eta", "#chi^{2}/ndof vs. #eta;#eta_{Track};#chi^{2}/ndof", 100, -etaMax_, etaMax_, 100, 0., 10.));
1062  vTrack2DHistos_.push_back(
1063  fs->make<TH2F>("h2_kappa_vs_phi", "#kappa vs. #phi;#phi_{Track};#kappa", 100, -M_PI, M_PI, 100, .0, .05));
1064  vTrack2DHistos_.push_back(
1065  fs->make<TH2F>("h2_kappa_vs_eta", "#kappa vs. #eta;#eta_{Track};#kappa", 100, -etaMax_, etaMax_, 100, .0, .05));
1066  vTrack2DHistos_.push_back(fs->make<TH2F>(
1067  "h2_normchi2_vs_kappa", "#kappa vs. #chi^{2}/ndof;#chi^{2}/ndof;#kappa", 100, 0., 10, 100, -.03, .03));
1068 
1069  firstEvent_ = true;
1070 
1071  } //beginJob
1072 
1073  //*************************************************************
1074  void endJob() override
1075  //*************************************************************
1076  {
1077  edm::LogPrint("GeneralPurposeTrackAnalyzer") << "*******************************" << std::endl;
1078  edm::LogPrint("GeneralPurposeTrackAnalyzer") << "Events run in total: " << ievt << std::endl;
1079  edm::LogPrint("GeneralPurposeTrackAnalyzer") << "n. tracks: " << itrks << std::endl;
1080  edm::LogPrint("GeneralPurposeTrackAnalyzer") << "*******************************" << std::endl;
1081 
1082  int nFiringTriggers = triggerMap_.size();
1083  edm::LogPrint("GeneralPurposeTrackAnalyzer") << "firing triggers: " << nFiringTriggers << std::endl;
1084  edm::LogPrint("GeneralPurposeTrackAnalyzer") << "*******************************" << std::endl;
1085 
1086  tksByTrigger_ = fs->make<TH1D>(
1087  "tksByTrigger", "tracks by HLT path;;% of # traks", nFiringTriggers, -0.5, nFiringTriggers - 0.5);
1088  evtsByTrigger_ = fs->make<TH1D>(
1089  "evtsByTrigger", "events by HLT path;;% of # events", nFiringTriggers, -0.5, nFiringTriggers - 0.5);
1090 
1091  int i = 0;
1092  for (const auto &it : triggerMap_) {
1093  i++;
1094 
1095  double trkpercent = ((it.second).second) * 100. / double(itrks);
1096  double evtpercent = ((it.second).first) * 100. / double(ievt);
1097 
1098  std::cout.precision(4);
1099 
1100  edm::LogPrint("GeneralPurposeTrackAnalyzer")
1101  << "HLT path: " << std::setw(60) << std::left << it.first << " | events firing: " << std::right
1102  << std::setw(8) << (it.second).first << " (" << std::setw(8) << std::fixed << std::setprecision(4)
1103  << evtpercent << "%)"
1104  << " | tracks collected: " << std::setw(10) << (it.second).second << " (" << std::setw(8) << std::fixed
1105  << std::setprecision(4) << trkpercent << "%)";
1106 
1107  tksByTrigger_->SetBinContent(i, trkpercent);
1108  tksByTrigger_->GetXaxis()->SetBinLabel(i, (it.first).c_str());
1109 
1110  evtsByTrigger_->SetBinContent(i, evtpercent);
1111  evtsByTrigger_->GetXaxis()->SetBinLabel(i, (it.first).c_str());
1112  }
1113 
1114  int nRuns = conditionsMap_.size();
1115 
1116  std::vector<int> theRuns_;
1117  for (const auto &it : conditionsMap_) {
1118  theRuns_.push_back(it.first);
1119  }
1120 
1121  sort(theRuns_.begin(), theRuns_.end());
1122  int runRange = theRuns_.back() - theRuns_.front() + 1;
1123 
1124  edm::LogPrint("GeneralPurposeTrackAnalyzer") << "*******************************" << std::endl;
1125  edm::LogPrint("GeneralPurposeTrackAnalyzer") << "first run: " << theRuns_.front() << std::endl;
1126  edm::LogPrint("GeneralPurposeTrackAnalyzer") << "last run: " << theRuns_.back() << std::endl;
1127  edm::LogPrint("GeneralPurposeTrackAnalyzer") << "considered runs: " << nRuns << std::endl;
1128  edm::LogPrint("GeneralPurposeTrackAnalyzer") << "*******************************" << std::endl;
1129 
1130  modeByRun_ = fs->make<TH1D>("modeByRun",
1131  "Strip APV mode by run number;;APV mode (-1=deco,+1=peak)",
1132  runRange,
1133  theRuns_.front() - 0.5,
1134  theRuns_.back() + 0.5);
1135 
1136  fieldByRun_ = fs->make<TH1D>("fieldByRun",
1137  "CMS B-field intensity by run number;;B-field intensity [T]",
1138  runRange,
1139  theRuns_.front() - 0.5,
1140  theRuns_.back() + 0.5);
1141 
1142  for (const auto &the_r : theRuns_) {
1143  if (conditionsMap_.find(the_r)->second.first != 0) {
1144  edm::LogPrint("GeneralPurposeTrackAnalyzer")
1145  << "run:" << the_r << " | isPeak: " << std::setw(4) << conditionsMap_.find(the_r)->second.first
1146  << "| B-field: " << conditionsMap_.find(the_r)->second.second << " [T]"
1147  << "| events: " << std::setw(10) << runInfoMap_.find(the_r)->second.first << ", tracks " << std::setw(10)
1148  << runInfoMap_.find(the_r)->second.second << std::endl;
1149  }
1150 
1151  modeByRun_->SetBinContent((the_r - theRuns_.front()) + 1, conditionsMap_.find(the_r)->second.first);
1152  fieldByRun_->SetBinContent((the_r - theRuns_.front()) + 1, conditionsMap_.find(the_r)->second.second);
1153  modeByRun_->GetXaxis()->SetBinLabel((the_r - theRuns_.front()) + 1, std::to_string(the_r).c_str());
1154  fieldByRun_->GetXaxis()->SetBinLabel((the_r - theRuns_.front()) + 1, std::to_string(the_r).c_str());
1155  }
1156 
1157  if (!isPhase1_) {
1158  pmap->save(true, 0, 0, "PixelHitMap.pdf", 600, 800);
1159  pmap->save(true, 0, 0, "PixelHitMap.png", 500, 750);
1160  }
1161 
1162  tmap->save(true, 0, 0, "StripHitMap.pdf");
1163  tmap->save(true, 0, 0, "StripHitMap.png");
1164 
1165  gStyle->SetPalette(kRainBow);
1166  pixelmap->beautifyAllHistograms();
1167 
1168  TCanvas cB("CanvBarrel", "CanvBarrel", 1200, 1000);
1169  pixelmap->drawBarrelMaps("entriesBarrel", cB);
1170  cB.SaveAs("pixelBarrelEntries.png");
1171 
1172  TCanvas cF("CanvForward", "CanvForward", 1600, 1000);
1173  pixelmap->drawForwardMaps("entriesForward", cF);
1174  cF.SaveAs("pixelForwardEntries.png");
1175 
1176  TCanvas cRocs = TCanvas("cRocs", "cRocs", 1200, 1600);
1177  pixelrocsmap_->drawMaps(cRocs, "Pixel on-track clusters occupancy");
1178  cRocs.SaveAs("Phase1PixelROCMaps_fullROCs.png");
1179  }
1180 
1181  //*************************************************************
1183  //*************************************************************
1184  {
1185  bool countStereoHitAs2D_ = true;
1186  // we count SiStrip stereo modules as 2D if selected via countStereoHitAs2D_
1187  // (since they provide theta information)
1188  if (!hit.isValid() ||
1189  (hit.dimension() < 2 && !countStereoHitAs2D_ && !dynamic_cast<const SiStripRecHit1D *>(&hit))) {
1190  return false; // real RecHit1D - but SiStripRecHit1D depends on countStereoHitAs2D_
1191  } else {
1192  const DetId detId(hit.geographicalId());
1193  if (detId.det() == DetId::Tracker) {
1194  if (detId.subdetId() == kBPIX || detId.subdetId() == kFPIX) {
1195  return true; // pixel is always 2D
1196  } else { // should be SiStrip now
1197  const SiStripDetId stripId(detId);
1198  if (stripId.stereo())
1199  return countStereoHitAs2D_; // stereo modules
1200  else if (dynamic_cast<const SiStripRecHit1D *>(&hit) || dynamic_cast<const SiStripRecHit2D *>(&hit))
1201  return false; // rphi modules hit
1202  //the following two are not used any more since ages...
1203  else if (dynamic_cast<const SiStripMatchedRecHit2D *>(&hit))
1204  return true; // matched is 2D
1205  else if (dynamic_cast<const ProjectedSiStripRecHit2D *>(&hit)) {
1206  const ProjectedSiStripRecHit2D *pH = static_cast<const ProjectedSiStripRecHit2D *>(&hit);
1207  return (countStereoHitAs2D_ && this->isHit2D(pH->originalHit())); // depends on original...
1208  } else {
1209  edm::LogError("UnknownType") << "@SUB=GeneralPurposeTrackAnalyzer::isHit2D"
1210  << "Tracker hit not in pixel, neither SiStripRecHit[12]D nor "
1211  << "SiStripMatchedRecHit2D nor ProjectedSiStripRecHit2D.";
1212  return false;
1213  }
1214  }
1215  } else { // not tracker??
1216  edm::LogWarning("DetectorMismatch") << "@SUB=GeneralPurposeTrackAnalyzer::isHit2D"
1217  << "Hit not in tracker with 'official' dimension >=2.";
1218  return true; // dimension() >= 2 so accept that...
1219  }
1220  }
1221  // never reached...
1222  }
1223 };
1224 
1225 //*************************************************************
1227 //*************************************************************
1228 {
1230  desc.setComment("Generic track analyzer to check ALCARECO sample quantities");
1231  desc.add<edm::InputTag>("TkTag", edm::InputTag("generalTracks"));
1232  desc.add<edm::InputTag>("TriggerResultsTag", edm::InputTag("TriggerResults", "", "HLT"));
1233  desc.add<edm::InputTag>("BeamSpotTag", edm::InputTag("offlineBeamSpot"));
1234  desc.add<edm::InputTag>("VerticesTag", edm::InputTag("offlinePrimaryVertices"));
1235  desc.add<bool>("isCosmics", false);
1236  descriptions.add("GeneralPurposeTrackAnalyzer", desc);
1237 }
1238 
static const std::string kSharedResource
Definition: TFileService.h:76
void analyze(const edm::Event &event, const edm::EventSetup &setup) override
double z0() const
z coordinate
Definition: BeamSpot.h:65
const edm::ESGetToken< SiStripLatency, SiStripLatencyRcd > latencyToken_
virtual int dimension() const =0
std::map< int, std::pair< int, float > > conditionsMap_
edm::EDGetTokenT< edm::TriggerResults > hltresultsToken_
edm::EDGetTokenT< reco::TrackCollection > theTrackCollectionToken_
uint32_t stereo() const
Definition: SiStripDetId.h:168
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
double y() const
y coordinate
Definition: Vertex.h:131
std::unique_ptr< TrackerMap > tmap
Geom::Phi< T > phi() const
Definition: PV3DBase.h:66
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
constexpr uint32_t rawId() const
get the raw id
Definition: DetId.h:57
std::vector< Track > TrackCollection
collection of Tracks
Definition: TrackFwd.h:14
T y() const
Definition: PV3DBase.h:60
static void fillDescriptions(edm::ConfigurationDescriptions &)
std::unique_ptr< TrackerMap > pmap
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
Log< level::Error, false > LogError
bool accept(const edm::Event &event, const edm::TriggerResults &triggerTable, const std::string &triggerPath)
Definition: TopDQMHelpers.h:31
Strings const & triggerNames() const
Definition: TriggerNames.cc:48
void beginRun(edm::Run const &run, edm::EventSetup const &setup) override
const edm::ESGetToken< TrackerGeometry, TrackerDigiGeometryRecord > geomTokenBR_
Geom::Theta< T > theta() const
Definition: PV3DBase.h:72
bool isThere(GeomDetEnumerators::SubDetector subdet) const
const int kFPIX
U second(std::pair< T, U > const &p)
void setComment(std::string const &value)
edm::ESHandle< MagneticField > magneticField_
std::map< std::string, std::pair< int, int > > triggerMap_
edm::EDGetTokenT< reco::VertexCollection > vertexToken_
void Fill(HcalDetId &id, double val, std::vector< TH2F > &depth)
T z() const
Definition: PV3DBase.h:61
GeneralPurposeTrackAnalyzer(const edm::ParameterSet &pset)
constexpr int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:48
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
double z() const
z coordinate
Definition: Vertex.h:133
Transition
Definition: Transition.h:12
void init(const TrackerTopology *, const TrackerGeometry *, const SiPixelFedCablingMap *)
int roc(const DetId &, const std::pair< int, int > &)
std::map< int, std::pair< int, int > > runInfoMap_
ParameterDescriptionBase * add(U const &iLabel, T const &value)
bool isValid() const
Definition: HandleBase.h:70
Log< level::Warning, true > LogPrint
const TrackerGeomDet * idToDet(DetId) const override
#define M_PI
edm::EDGetTokenT< reco::BeamSpot > beamspotToken_
SiStripRecHit2D originalHit() const
Log< level::Info, false > LogInfo
Detector identifier class for the strip tracker.
Definition: SiStripDetId.h:18
const edm::ESGetToken< MagneticField, IdealMagneticFieldRecord > magFieldToken_
Definition: DetId.h:17
double x() const
x coordinate
Definition: Vertex.h:129
tuple trackCollection
T const * product() const
Definition: Handle.h:70
std::string const & triggerName(unsigned int index) const
Definition: TriggerNames.cc:50
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
const edm::ESGetToken< TrackerTopology, TrackerTopologyRcd > trackerTopologyTokenBR_
T const * product() const
Definition: ESHandle.h:86
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
bool isValid() const
void add(std::string const &label, ParameterSetDescription const &psetDescription)
static const std::string algoNames[]
Definition: TrackBase.h:147
void endRun(edm::Run const &, edm::EventSetup const &) override
int GetIndex(const std::vector< OBJECT_TYPE * > &vec, const TString &name)
std::unique_ptr< Phase1PixelROCMaps > pixelrocsmap_
double y0() const
y coordinate
Definition: BeamSpot.h:63
tuple cout
Definition: gather_cfg.py:144
DetId geographicalId() const
Log< level::Warning, false > LogWarning
std::unique_ptr< Phase1PixelMaps > pixelmap
int col
Definition: cuy.py:1009
T x() const
Definition: PV3DBase.h:59
bool isHit2D(const TrackingRecHit &hit)
~GeneralPurposeTrackAnalyzer() override=default
*vegas h *****************************************************used in the default bin number in original ***version of VEGAS is ***a higher bin number might help to derive a more precise ***grade subtle point
Definition: invegas.h:5
const edm::ESGetToken< SiPixelFedCablingMap, SiPixelFedCablingMapRcd > siPixelFedCablingMapTokenBR_
Definition: Run.h:45
Our base class.
Definition: SiPixelRecHit.h:23
const edm::ESGetToken< TrackerGeometry, TrackerDigiGeometryRecord > geomToken_
double x0() const
x coordinate
Definition: BeamSpot.h:61
const int kBPIX