CMS 3D CMS Logo

PrimaryVertexValidation.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: Alignment/OfflineValidation
4 // Class: PrimaryVertexValidation
5 //
13 //
14 // Original Author: Marco Musich
15 // Created: Tue Mar 02 10:39:34 CDT 2010
16 //
17 
18 // system include files
19 #include <memory>
20 #include <vector>
21 #include <regex>
22 #include <cassert>
23 #include <chrono>
24 #include <iomanip>
25 #include <boost/range/adaptor/indexed.hpp>
26 
27 // user include files
29 
30 // ROOT includes
31 #include "TH1F.h"
32 #include "TH2F.h"
33 #include "TF1.h"
34 #include "TVector3.h"
35 #include "TFile.h"
36 #include "TMath.h"
37 #include "TROOT.h"
38 #include "TChain.h"
39 #include "TNtuple.h"
40 #include "TMatrixD.h"
41 #include "TVectorD.h"
42 
43 // CMSSW includes
65 
67 
68 // Constructor
72  ttkToken_(esConsumes<TransientTrackBuilder, TransientTrackRecord>(edm::ESInputTag("", "TransientTrackBuilder"))),
74  topoTokenBR_(esConsumes<TrackerTopology, TrackerTopologyRcd, edm::Transition::BeginRun>()),
76  runInfoTokenBR_(esConsumes<RunInfo, RunInfoRcd, edm::Transition::BeginRun>()),
77  compressionSettings_(iConfig.getUntrackedParameter<int>("compressionSettings", -1)),
78  storeNtuple_(iConfig.getParameter<bool>("storeNtuple")),
79  lightNtupleSwitch_(iConfig.getParameter<bool>("isLightNtuple")),
80  useTracksFromRecoVtx_(iConfig.getParameter<bool>("useTracksFromRecoVtx")),
81  vertexZMax_(iConfig.getUntrackedParameter<double>("vertexZMax", 99.)),
82  intLumi_(iConfig.getUntrackedParameter<double>("intLumi", 0.)),
83  askFirstLayerHit_(iConfig.getParameter<bool>("askFirstLayerHit")),
84  doBPix_(iConfig.getUntrackedParameter<bool>("doBPix", true)),
85  doFPix_(iConfig.getUntrackedParameter<bool>("doFPix", true)),
86  ptOfProbe_(iConfig.getUntrackedParameter<double>("probePt", 0.)),
87  pOfProbe_(iConfig.getUntrackedParameter<double>("probeP", 0.)),
88  etaOfProbe_(iConfig.getUntrackedParameter<double>("probeEta", 2.4)),
89  nHitsOfProbe_(iConfig.getUntrackedParameter<double>("probeNHits", 0.)),
90  nBins_(iConfig.getUntrackedParameter<int>("numberOfBins", 24)),
91  minPt_(iConfig.getUntrackedParameter<double>("minPt", 1.)),
92  maxPt_(iConfig.getUntrackedParameter<double>("maxPt", 20.)),
93  debug_(iConfig.getParameter<bool>("Debug")),
94  runControl_(iConfig.getUntrackedParameter<bool>("runControl", false)),
95  forceBeamSpotContraint_(iConfig.getUntrackedParameter<bool>("forceBeamSpot", false)) {
96  // now do what ever initialization is needed
97  // initialize phase space boundaries
98 
99  usesResource(TFileService::kSharedResource);
100 
101  std::vector<unsigned int> defaultRuns;
102  defaultRuns.push_back(0);
103  runControlNumbers_ = iConfig.getUntrackedParameter<std::vector<unsigned int>>("runControlNumber", defaultRuns);
104 
105  edm::InputTag TrackCollectionTag_ = iConfig.getParameter<edm::InputTag>("TrackCollectionTag");
106  theTrackCollectionToken_ = consumes<reco::TrackCollection>(TrackCollectionTag_);
107 
108  edm::InputTag VertexCollectionTag_ = iConfig.getParameter<edm::InputTag>("VertexCollectionTag");
109  theVertexCollectionToken_ = consumes<reco::VertexCollection>(VertexCollectionTag_);
110 
111  edm::InputTag BeamspotTag_ = iConfig.getParameter<edm::InputTag>("BeamSpotTag");
112  theBeamspotToken_ = consumes<reco::BeamSpot>(BeamspotTag_);
113 
114  // select and configure the track filter
116  std::make_unique<TrackFilterForPVFinding>(iConfig.getParameter<edm::ParameterSet>("TkFilterParameters"));
117  // select and configure the track clusterizer
118  std::string clusteringAlgorithm =
119  iConfig.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<std::string>("algorithm");
120  if (clusteringAlgorithm == "gap") {
122  std::make_unique<GapClusterizerInZ>(iConfig.getParameter<edm::ParameterSet>("TkClusParameters")
123  .getParameter<edm::ParameterSet>("TkGapClusParameters"));
124  } else if (clusteringAlgorithm == "DA_vect") {
126  std::make_unique<DAClusterizerInZ_vect>(iConfig.getParameter<edm::ParameterSet>("TkClusParameters")
127  .getParameter<edm::ParameterSet>("TkDAClusParameters"));
128  } else {
129  throw VertexException("PrimaryVertexProducerAlgorithm: unknown clustering algorithm: " + clusteringAlgorithm);
130  }
131 
132  theDetails_.histobins = 500;
139 
140  for (int i = PVValHelper::phi; i < PVValHelper::END_OF_PLOTS; i++) {
141  for (int j = PVValHelper::dx; j < PVValHelper::END_OF_TYPES; j++) {
142  auto plot_index = static_cast<PVValHelper::plotVariable>(i);
143  auto res_index = static_cast<PVValHelper::residualType>(j);
144 
145  if (debug_) {
146  edm::LogInfo("PrimaryVertexValidation")
147  << "==> " << std::get<0>(PVValHelper::getTypeString(res_index)) << " " << std::setw(10)
148  << std::get<0>(PVValHelper::getVarString(plot_index)) << std::endl;
149  }
150  if (res_index != PVValHelper::d3D && res_index != PVValHelper::norm_d3D)
151  theDetails_.setMap(res_index,
152  plot_index,
153  theDetails_.getLow(PVValHelper::dxy, plot_index),
154  theDetails_.getHigh(PVValHelper::dxy, plot_index));
155  else
156  theDetails_.setMap(res_index, plot_index, 0., theDetails_.getHigh(PVValHelper::dxy, plot_index));
157  }
158  }
159 
160  edm::LogVerbatim("PrimaryVertexValidation") << "######################################";
161  for (const auto& it : theDetails_.range) {
162  edm::LogVerbatim("PrimaryVertexValidation")
163  << "|" << std::setw(10) << std::get<0>(PVValHelper::getTypeString(it.first.first)) << "|" << std::setw(10)
164  << std::get<0>(PVValHelper::getVarString(it.first.second)) << "| (" << std::setw(5) << it.second.first << ";"
165  << std::setw(5) << it.second.second << ") |" << std::endl;
166  }
167 
170 
171  if (debug_) {
172  edm::LogVerbatim("PrimaryVertexValidation") << "etaBins: ";
174  edm::LogVerbatim("PrimaryVertexValidation") << ieta << " ";
175  }
176  edm::LogVerbatim("PrimaryVertexValidation") << "\n";
177 
178  edm::LogVerbatim("PrimaryVertexValidation") << "phiBins: ";
180  edm::LogVerbatim("PrimaryVertexValidation") << iphi << " ";
181  }
182  edm::LogVerbatim("PrimaryVertexValidation") << "\n";
183  }
184 
185  // create the bins of the pT-binned distributions
186 
187  mypT_bins_ = PVValHelper::makeLogBins<float, nPtBins_>(minPt_, maxPt_);
188 
189  std::string toOutput = "";
190  for (auto ptbin : mypT_bins_) {
191  toOutput += " ";
192  toOutput += std::to_string(ptbin);
193  toOutput += ",";
194  }
195 
196  edm::LogVerbatim("PrimaryVertexValidation") << "######################################\n";
197  edm::LogVerbatim("PrimaryVertexValidation") << "The pT binning is: [" << toOutput << "] \n";
198 }
199 
200 // Destructor
202 //
203 // member functions
204 //
205 
206 // ------------ method called to for each event ------------
208  using namespace std;
209  using namespace reco;
210  using namespace IPTools;
211 
212  if (nBins_ != 24 && debug_) {
213  edm::LogInfo("PrimaryVertexValidation") << "Using: " << nBins_ << " bins plots";
214  }
215 
216  bool passesRunControl = false;
217 
218  if (runControl_) {
219  for (const auto& runControlNumber : runControlNumbers_) {
220  if (iEvent.eventAuxiliary().run() == runControlNumber) {
221  if (debug_) {
222  edm::LogInfo("PrimaryVertexValidation")
223  << " run number: " << iEvent.eventAuxiliary().run() << " keeping run:" << runControlNumber;
224  }
225  passesRunControl = true;
226  break;
227  }
228  }
229  if (!passesRunControl)
230  return;
231  }
232 
233  Nevt_++;
234 
235  //=======================================================
236  // Initialize Root-tuple variables
237  //=======================================================
238 
239  SetVarToZero();
240 
241  //=======================================================
242  // Retrieve tracker topology from geometry
243  //=======================================================
244 
245  const TrackerTopology* const tTopo = &iSetup.getData(topoToken_);
246 
247  //=======================================================
248  // Retrieve the Magnetic Field information
249  //=======================================================
250 
252 
253  //=======================================================
254  // Retrieve the Tracking Geometry information
255  //=======================================================
256 
258 
259  //=======================================================
260  // Retrieve the Transient Track Builder information
261  //=======================================================
262 
264  double fBfield_ = ((*theB_).field()->inTesla(GlobalPoint(0., 0., 0.))).z();
265 
266  //=======================================================
267  // Retrieve the Track information
268  //=======================================================
269 
270  edm::Handle<TrackCollection> trackCollectionHandle = iEvent.getHandle(theTrackCollectionToken_);
271  if (!trackCollectionHandle.isValid())
272  return;
273  auto const& tracks = *trackCollectionHandle;
274 
275  //=======================================================
276  // Retrieve offline vartex information (only for reco)
277  //=======================================================
278 
279  //edm::Handle<VertexCollection> vertices;
282  if (!vertices.isValid()) {
283  edm::LogError("PrimaryVertexValidation") << "Vertex collection handle is not valid. Aborting!" << std::endl;
284  return;
285  }
286 
287  std::vector<Vertex> vsorted = *(vertices);
288  // sort the vertices by number of tracks in descending order
289  // use chi2 as tiebreaker
290  std::sort(vsorted.begin(), vsorted.end(), PrimaryVertexValidation::vtxSort);
291 
292  // skip events with no PV, this should not happen
293  if (vsorted.empty())
294  return;
295 
296  // skip events failing vertex cut
297  if (std::abs(vsorted[0].z()) > vertexZMax_)
298  return;
299 
300  if (vsorted[0].isValid()) {
301  xOfflineVertex_ = (vsorted)[0].x();
302  yOfflineVertex_ = (vsorted)[0].y();
303  zOfflineVertex_ = (vsorted)[0].z();
304 
305  xErrOfflineVertex_ = (vsorted)[0].xError();
306  yErrOfflineVertex_ = (vsorted)[0].yError();
307  zErrOfflineVertex_ = (vsorted)[0].zError();
308  }
309 
316 
317  unsigned int vertexCollectionSize = vsorted.size();
318  int nvvertex = 0;
319 
320  for (unsigned int i = 0; i < vertexCollectionSize; i++) {
321  const Vertex& vertex = vsorted.at(i);
322  if (vertex.isValid())
323  nvvertex++;
324  }
325 
326  nOfflineVertices_ = nvvertex;
327  h_nOfflineVertices->Fill(nvvertex);
328 
329  if (!vsorted.empty() && useTracksFromRecoVtx_) {
330  double sumpt = 0;
331  size_t ntracks = 0;
332  double chi2ndf = 0.;
333  double chi2prob = 0.;
334 
335  if (!vsorted.at(0).isFake()) {
336  Vertex pv = vsorted.at(0);
337 
338  ntracks = pv.tracksSize();
339  chi2ndf = pv.normalizedChi2();
340  chi2prob = TMath::Prob(pv.chi2(), (int)pv.ndof());
341 
342  h_recoVtxNtracks_->Fill(ntracks);
343  h_recoVtxChi2ndf_->Fill(chi2ndf);
344  h_recoVtxChi2Prob_->Fill(chi2prob);
345 
346  for (Vertex::trackRef_iterator itrk = pv.tracks_begin(); itrk != pv.tracks_end(); ++itrk) {
347  double pt = (**itrk).pt();
348  sumpt += pt * pt;
349 
350  const math::XYZPoint myVertex(pv.position().x(), pv.position().y(), pv.position().z());
351 
352  double dxyRes = (**itrk).dxy(myVertex);
353  double dzRes = (**itrk).dz(myVertex);
354 
355  double dxy_err = (**itrk).dxyError();
356  double dz_err = (**itrk).dzError();
357 
358  float trackphi = ((**itrk).phi()) * (180 / M_PI);
359  float tracketa = (**itrk).eta();
360 
361  for (int i = 0; i < nBins_; i++) {
362  float phiF = theDetails_.trendbins[PVValHelper::phi][i];
363  float phiL = theDetails_.trendbins[PVValHelper::phi][i + 1];
364 
365  float etaF = theDetails_.trendbins[PVValHelper::eta][i];
366  float etaL = theDetails_.trendbins[PVValHelper::eta][i + 1];
367 
368  if (tracketa >= etaF && tracketa < etaL) {
371  PVValHelper::fillByIndex(n_dxyEtaBiasResiduals, i, (dxyRes) / dxy_err, "3");
372  PVValHelper::fillByIndex(n_dzEtaBiasResiduals, i, (dzRes) / dz_err, "4");
373  }
374 
375  if (trackphi >= phiF && trackphi < phiL) {
378  PVValHelper::fillByIndex(n_dxyPhiBiasResiduals, i, (dxyRes) / dxy_err, "7");
379  PVValHelper::fillByIndex(n_dzPhiBiasResiduals, i, (dzRes) / dz_err, "8");
380 
381  for (int j = 0; j < nBins_; j++) {
382  float etaJ = theDetails_.trendbins[PVValHelper::eta][j];
383  float etaK = theDetails_.trendbins[PVValHelper::eta][j + 1];
384 
385  if (tracketa >= etaJ && tracketa < etaK) {
386  a_dxyBiasResidualsMap[i][j]->Fill(dxyRes * cmToum);
387  a_dzBiasResidualsMap[i][j]->Fill(dzRes * cmToum);
388 
389  n_dxyBiasResidualsMap[i][j]->Fill((dxyRes) / dxy_err);
390  n_dzBiasResidualsMap[i][j]->Fill((dzRes) / dz_err);
391  }
392  }
393  }
394  }
395  }
396 
397  h_recoVtxSumPt_->Fill(sumpt);
398  }
399  }
400 
401  //=======================================================
402  // Retrieve Beamspot information
403  //=======================================================
404 
406  edm::Handle<BeamSpot> beamSpotHandle = iEvent.getHandle(theBeamspotToken_);
407 
408  if (beamSpotHandle.isValid()) {
409  beamSpot = *beamSpotHandle;
410  BSx0_ = beamSpot.x0();
411  BSy0_ = beamSpot.y0();
412  BSz0_ = beamSpot.z0();
413  Beamsigmaz_ = beamSpot.sigmaZ();
414  Beamdxdz_ = beamSpot.dxdz();
415  BeamWidthX_ = beamSpot.BeamWidthX();
416  BeamWidthY_ = beamSpot.BeamWidthY();
417 
418  wxy2_ = TMath::Power(BeamWidthX_, 2) + TMath::Power(BeamWidthY_, 2);
419 
420  } else {
421  edm::LogWarning("PrimaryVertexValidation") << "No BeamSpot found!";
422  }
423 
424  h_BSx0->Fill(BSx0_);
425  h_BSy0->Fill(BSy0_);
426  h_BSz0->Fill(BSz0_);
427  h_Beamsigmaz->Fill(Beamsigmaz_);
428  h_BeamWidthX->Fill(BeamWidthX_);
429  h_BeamWidthY->Fill(BeamWidthY_);
430 
431  if (debug_)
432  edm::LogInfo("PrimaryVertexValidation") << "Beamspot x:" << BSx0_ << " y:" << BSy0_ << " z:" << BSz0_;
433 
434  //=======================================================
435  // Starts here ananlysis
436  //=======================================================
437 
438  RunNumber_ = iEvent.eventAuxiliary().run();
439  LuminosityBlockNumber_ = iEvent.eventAuxiliary().luminosityBlock();
440  EventNumber_ = iEvent.eventAuxiliary().id().event();
441 
442  if (debug_)
443  edm::LogInfo("PrimaryVertexValidation") << " looping over " << trackCollectionHandle->size() << "tracks";
444 
445  h_nTracks->Fill(trackCollectionHandle->size());
446 
447  //======================================================
448  // Interface RECO tracks to vertex reconstruction
449  //======================================================
450 
451  std::vector<TransientTrack> t_tks;
452  for (const auto& track : tracks) {
453  TransientTrack tt = theB_->build(&(track));
454  tt.setBeamSpot(beamSpot);
455  t_tks.push_back(tt);
456  }
457 
458  if (debug_) {
459  edm::LogInfo("PrimaryVertexValidation") << "Found: " << t_tks.size() << " reconstructed tracks";
460  }
461 
462  //======================================================
463  // select the tracks
464  //======================================================
465 
466  std::vector<TransientTrack> seltks = theTrackFilter_->select(t_tks);
467 
468  //======================================================
469  // clusterize tracks in Z
470  //======================================================
471 
472  vector<vector<TransientTrack>> clusters = theTrackClusterizer_->clusterize(seltks);
473 
474  if (debug_) {
475  edm::LogInfo("PrimaryVertexValidation")
476  << " looping over: " << clusters.size() << " clusters from " << t_tks.size() << " selected tracks";
477  }
478 
479  nClus_ = clusters.size();
480  h_nClus->Fill(nClus_);
481 
482  //======================================================
483  // Starts loop on clusters
484  //======================================================
485  for (const auto& iclus : clusters) {
486  nTracksPerClus_ = 0;
487 
488  unsigned int i = 0;
489  for (const auto& theTTrack : iclus) {
490  i++;
491 
492  if (nTracks_ >= nMaxtracks_) {
493  edm::LogError("PrimaryVertexValidation")
494  << " Warning - Number of tracks: " << nTracks_ << " , greater than " << nMaxtracks_;
495  continue;
496  }
497 
498  const Track& theTrack = theTTrack.track();
499 
500  pt_[nTracks_] = theTrack.pt();
501  p_[nTracks_] = theTrack.p();
502  nhits_[nTracks_] = theTrack.numberOfValidHits();
503  eta_[nTracks_] = theTrack.eta();
504  theta_[nTracks_] = theTrack.theta();
505  phi_[nTracks_] = theTrack.phi();
506  chi2_[nTracks_] = theTrack.chi2();
507  chi2ndof_[nTracks_] = theTrack.normalizedChi2();
508  charge_[nTracks_] = theTrack.charge();
509  qoverp_[nTracks_] = theTrack.qoverp();
510  dz_[nTracks_] = theTrack.dz();
511  dxy_[nTracks_] = theTrack.dxy();
512 
513  TrackBase::TrackQuality _trackQuality = TrackBase::qualityByName("highPurity");
514  isHighPurity_[nTracks_] = theTrack.quality(_trackQuality);
515 
517  dxyBs_[nTracks_] = theTrack.dxy(point);
518  dzBs_[nTracks_] = theTrack.dz(point);
519 
520  xPCA_[nTracks_] = theTrack.vertex().x();
521  yPCA_[nTracks_] = theTrack.vertex().y();
522  zPCA_[nTracks_] = theTrack.vertex().z();
523 
524  //=======================================================
525  // Retrieve rechit information
526  //=======================================================
527 
528  const reco::HitPattern& hits = theTrack.hitPattern();
529 
530  int nRecHit1D = 0;
531  int nRecHit2D = 0;
532  int nhitinTIB = hits.numberOfValidStripTIBHits();
533  int nhitinTOB = hits.numberOfValidStripTOBHits();
534  int nhitinTID = hits.numberOfValidStripTIDHits();
535  int nhitinTEC = hits.numberOfValidStripTECHits();
536  int nhitinBPIX = hits.numberOfValidPixelBarrelHits();
537  int nhitinFPIX = hits.numberOfValidPixelEndcapHits();
538  for (trackingRecHit_iterator iHit = theTTrack.recHitsBegin(); iHit != theTTrack.recHitsEnd(); ++iHit) {
539  if ((*iHit)->isValid()) {
540  if (this->isHit2D(**iHit, phase_)) {
541  ++nRecHit2D;
542  } else {
543  ++nRecHit1D;
544  }
545  }
546  }
547 
548  nhits1D_[nTracks_] = nRecHit1D;
549  nhits2D_[nTracks_] = nRecHit2D;
550  nhitsBPIX_[nTracks_] = nhitinBPIX;
551  nhitsFPIX_[nTracks_] = nhitinFPIX;
552  nhitsTIB_[nTracks_] = nhitinTIB;
553  nhitsTID_[nTracks_] = nhitinTID;
554  nhitsTOB_[nTracks_] = nhitinTOB;
555  nhitsTEC_[nTracks_] = nhitinTEC;
556 
557  //=======================================================
558  // Good tracks for vertexing selection
559  //=======================================================
560 
561  bool pass = true;
562  if (askFirstLayerHit_)
563  pass = this->hasFirstLayerPixelHits(theTTrack);
564  if (pass && (theTrack.pt() >= ptOfProbe_) && std::abs(theTrack.eta()) <= etaOfProbe_ &&
565  (theTrack.numberOfValidHits()) >= nHitsOfProbe_ && (theTrack.p()) >= pOfProbe_) {
566  isGoodTrack_[nTracks_] = 1;
567  }
568 
569  //=======================================================
570  // Fit unbiased vertex
571  //=======================================================
572 
573  vector<TransientTrack> theFinalTracks;
574  theFinalTracks.clear();
575 
576  for (const auto& tk : iclus) {
577  pass = this->hasFirstLayerPixelHits(tk);
578  if (pass) {
579  if (tk == theTTrack)
580  continue;
581  else {
582  theFinalTracks.push_back(tk);
583  }
584  }
585  }
586 
587  if (theFinalTracks.size() > 1) {
588  if (debug_)
589  edm::LogInfo("PrimaryVertexValidation") << "Transient Track Collection size: " << theFinalTracks.size();
590  try {
591  //AdaptiveVertexFitter* theFitter = new AdaptiveVertexFitter;
592  auto theFitter = std::unique_ptr<VertexFitter<5>>(new AdaptiveVertexFitter());
593  TransientVertex theFittedVertex;
594 
596  theFittedVertex = theFitter->vertex(theFinalTracks, beamSpot); // if you want the beam constraint
597  } else {
598  theFittedVertex = theFitter->vertex(theFinalTracks);
599  }
600 
601  double totalTrackWeights = 0;
602  if (theFittedVertex.isValid()) {
603  if (theFittedVertex.hasTrackWeight()) {
604  for (const auto& theFinalTrack : theFinalTracks) {
605  sumOfWeightsUnbiasedVertex_[nTracks_] += theFittedVertex.trackWeight(theFinalTrack);
606  totalTrackWeights += theFittedVertex.trackWeight(theFinalTrack);
607  h_fitVtxTrackWeights_->Fill(theFittedVertex.trackWeight(theFinalTrack));
608  }
609  }
610 
611  h_fitVtxTrackAverageWeight_->Fill(totalTrackWeights / theFinalTracks.size());
612 
614  const math::XYZPoint myVertex(
615  theFittedVertex.position().x(), theFittedVertex.position().y(), theFittedVertex.position().z());
616 
617  const Vertex vertex = theFittedVertex;
618  fillTrackHistos(hDA, "all", &theTTrack, vertex, beamSpot, fBfield_);
619 
620  hasRecVertex_[nTracks_] = 1;
621  xUnbiasedVertex_[nTracks_] = theFittedVertex.position().x();
622  yUnbiasedVertex_[nTracks_] = theFittedVertex.position().y();
623  zUnbiasedVertex_[nTracks_] = theFittedVertex.position().z();
624 
625  chi2normUnbiasedVertex_[nTracks_] = theFittedVertex.normalisedChiSquared();
626  chi2UnbiasedVertex_[nTracks_] = theFittedVertex.totalChiSquared();
627  DOFUnbiasedVertex_[nTracks_] = theFittedVertex.degreesOfFreedom();
629  TMath::Prob(theFittedVertex.totalChiSquared(), (int)theFittedVertex.degreesOfFreedom());
630  tracksUsedForVertexing_[nTracks_] = theFinalTracks.size();
631 
632  h_fitVtxNtracks_->Fill(theFinalTracks.size());
633  h_fitVtxChi2_->Fill(theFittedVertex.totalChiSquared());
634  h_fitVtxNdof_->Fill(theFittedVertex.degreesOfFreedom());
635  h_fitVtxChi2ndf_->Fill(theFittedVertex.normalisedChiSquared());
636  h_fitVtxChi2Prob_->Fill(
637  TMath::Prob(theFittedVertex.totalChiSquared(), (int)theFittedVertex.degreesOfFreedom()));
638 
639  // from my Vertex
640  double dxyFromMyVertex = theTrack.dxy(myVertex);
641  double dzFromMyVertex = theTrack.dz(myVertex);
642 
643  GlobalPoint vert(
644  theFittedVertex.position().x(), theFittedVertex.position().y(), theFittedVertex.position().z());
645 
646  //FreeTrajectoryState theTrackNearVertex = theTTrack.trajectoryStateClosestToPoint(vert).theState();
647  //double dz_err = sqrt(theFittedVertex.positionError().czz() + theTrackNearVertex.cartesianError().position().czz());
648  //double dz_err = hypot(theTrack.dzError(),theFittedVertex.positionError().czz());
649 
650  double dz_err = sqrt(std::pow(theTrack.dzError(), 2) + theFittedVertex.positionError().czz());
651 
652  // PV2D
653  std::pair<bool, Measurement1D> s_ip2dpv = signedTransverseImpactParameter(
654  theTTrack, GlobalVector(theTrack.px(), theTrack.py(), theTrack.pz()), theFittedVertex);
655 
656  double s_ip2dpv_corr = s_ip2dpv.second.value();
657  double s_ip2dpv_err = s_ip2dpv.second.error();
658 
659  // PV3D
660  std::pair<bool, Measurement1D> s_ip3dpv = signedImpactParameter3D(
661  theTTrack, GlobalVector(theTrack.px(), theTrack.py(), theTrack.pz()), theFittedVertex);
662 
663  double s_ip3dpv_corr = s_ip3dpv.second.value();
664  double s_ip3dpv_err = s_ip3dpv.second.error();
665 
666  // PV3D absolute
667  std::pair<bool, Measurement1D> ip3dpv = absoluteImpactParameter3D(theTTrack, theFittedVertex);
668  double ip3d_corr = ip3dpv.second.value();
669  double ip3d_err = ip3dpv.second.error();
670 
671  // with respect to any specified vertex, such as primary vertex
673 
674  GlobalPoint refPoint = traj.position();
675  GlobalPoint cPToVtx = traj.theState().position();
676 
677  float my_dx = refPoint.x() - myVertex.x();
678  float my_dy = refPoint.y() - myVertex.y();
679 
680  float my_dx2 = cPToVtx.x() - myVertex.x();
681  float my_dy2 = cPToVtx.y() - myVertex.y();
682 
683  float my_dxy = std::sqrt(my_dx * my_dx + my_dy * my_dy);
684 
686  //double d0_error = traj.perigeeError().transverseImpactParameterError();
688  double z0_error = traj.perigeeError().longitudinalImpactParameterError();
689 
690  if (debug_) {
691  edm::LogInfo("PrimaryVertexValidation")
692  << "my_dx:" << my_dx << " my_dy:" << my_dy << " my_dxy:" << my_dxy << " my_dx2:" << my_dx2
693  << " my_dy2:" << my_dy2 << " d0: " << d0 << " dxyFromVtx:" << dxyFromMyVertex << "\n"
694  << " ============================== "
695  << "\n"
696  << "diff1:" << std::abs(d0) - std::abs(my_dxy) << "\n"
697  << "diff2:" << std::abs(d0) - std::abs(dxyFromMyVertex) << "\n"
698  << "diff3:" << (my_dx - my_dx2) << " " << (my_dy - my_dy2) << "\n"
699  << std::endl;
700  }
701 
702  // define IPs
703 
704  dxyFromMyVertex_[nTracks_] = dxyFromMyVertex;
705  dxyErrorFromMyVertex_[nTracks_] = s_ip2dpv_err;
706  IPTsigFromMyVertex_[nTracks_] = dxyFromMyVertex / s_ip2dpv_err;
707 
708  dzFromMyVertex_[nTracks_] = dzFromMyVertex;
709  dzErrorFromMyVertex_[nTracks_] = dz_err;
710  IPLsigFromMyVertex_[nTracks_] = dzFromMyVertex / dz_err;
711 
712  d3DFromMyVertex_[nTracks_] = ip3d_corr;
713  d3DErrorFromMyVertex_[nTracks_] = ip3d_err;
714  IP3DsigFromMyVertex_[nTracks_] = (ip3d_corr / ip3d_err);
715 
716  // fill directly the histograms of residuals
717 
718  float trackphi = (theTrack.phi()) * (180. / M_PI);
719  float tracketa = theTrack.eta();
720  float trackpt = theTrack.pt();
721  float trackp = theTrack.p();
722  float tracknhits = theTrack.numberOfValidHits();
723 
724  // determine the module number and ladder
725 
726  int ladder_num = -1.;
727  int module_num = -1.;
728  int L1BPixHitCount = 0;
729 
730  for (auto const& hit : theTrack.recHits()) {
731  const DetId& detId = hit->geographicalId();
732  unsigned int subid = detId.subdetId();
733 
734  if (hit->isValid() && (subid == PixelSubdetector::PixelBarrel)) {
735  int layer = tTopo->pxbLayer(detId);
736  if (layer == 1) {
737  const SiPixelRecHit* prechit = dynamic_cast<const SiPixelRecHit*>(
738  hit); //to be used to get the associated cluster and the cluster probability
739  double clusterProbability = prechit->clusterProbability(0);
740  if (clusterProbability > 0) {
741  h_probeL1ClusterProb_->Fill(log10(clusterProbability));
742  }
743 
744  L1BPixHitCount += 1;
745  ladder_num = tTopo->pxbLadder(detId);
746  module_num = tTopo->pxbModule(detId);
747  }
748  }
749  }
750 
751  h_probeL1Ladder_->Fill(ladder_num);
752  h_probeL1Module_->Fill(module_num);
753  h2_probeLayer1Map_->Fill(module_num, ladder_num);
754  h_probeHasBPixL1Overlap_->Fill(L1BPixHitCount);
755 
756  // residuals vs ladder and module number for map
757  if (module_num > 0 && ladder_num > 0) { // only if we are on BPix Layer 1
758  a_dxyL1ResidualsMap[ladder_num - 1][module_num - 1]->Fill(dxyFromMyVertex * cmToum);
759  a_dzL1ResidualsMap[ladder_num - 1][module_num - 1]->Fill(dzFromMyVertex * cmToum);
760  n_dxyL1ResidualsMap[ladder_num - 1][module_num - 1]->Fill(dxyFromMyVertex / s_ip2dpv_err);
761  n_dzL1ResidualsMap[ladder_num - 1][module_num - 1]->Fill(dzFromMyVertex / dz_err);
762  }
763 
764  // filling the pT-binned distributions
765 
766  for (int ipTBin = 0; ipTBin < nPtBins_; ipTBin++) {
767  float pTF = mypT_bins_[ipTBin];
768  float pTL = mypT_bins_[ipTBin + 1];
769 
770  if (debug_)
771  edm::LogInfo("PrimaryVertexValidation") << "ipTBin:" << ipTBin << " " << mypT_bins_[ipTBin]
772  << " < pT < " << mypT_bins_[ipTBin + 1] << std::endl;
773 
774  if (std::abs(tracketa) < 1.5 && (trackpt >= pTF && trackpt < pTL)) {
775  if (debug_)
776  edm::LogInfo("PrimaryVertexValidation") << "passes this cut: " << mypT_bins_[ipTBin] << std::endl;
777  PVValHelper::fillByIndex(h_dxy_pT_, ipTBin, dxyFromMyVertex * cmToum, "9");
778  PVValHelper::fillByIndex(h_dz_pT_, ipTBin, dzFromMyVertex * cmToum, "10");
779  PVValHelper::fillByIndex(h_norm_dxy_pT_, ipTBin, dxyFromMyVertex / s_ip2dpv_err, "11");
780  PVValHelper::fillByIndex(h_norm_dz_pT_, ipTBin, dzFromMyVertex / dz_err, "12");
781 
782  if (std::abs(tracketa) < 1.) {
783  if (debug_)
784  edm::LogInfo("PrimaryVertexValidation")
785  << "passes tight eta cut: " << mypT_bins_[ipTBin] << std::endl;
786  PVValHelper::fillByIndex(h_dxy_Central_pT_, ipTBin, dxyFromMyVertex * cmToum, "13");
787  PVValHelper::fillByIndex(h_dz_Central_pT_, ipTBin, dzFromMyVertex * cmToum, "14");
788  PVValHelper::fillByIndex(h_norm_dxy_Central_pT_, ipTBin, dxyFromMyVertex / s_ip2dpv_err, "15");
789  PVValHelper::fillByIndex(h_norm_dz_Central_pT_, ipTBin, dzFromMyVertex / dz_err, "16");
790  }
791  }
792  }
793 
794  // checks on the probe track quality
795  if (trackpt >= ptOfProbe_ && std::abs(tracketa) <= etaOfProbe_ && tracknhits >= nHitsOfProbe_ &&
796  trackp >= pOfProbe_) {
797  std::pair<bool, bool> pixelOcc = pixelHitsCheck((theTTrack));
798 
799  if (debug_) {
800  if (pixelOcc.first == true)
801  edm::LogInfo("PrimaryVertexValidation") << "has BPIx hits" << std::endl;
802  if (pixelOcc.second == true)
803  edm::LogInfo("PrimaryVertexValidation") << "has FPix hits" << std::endl;
804  }
805 
806  if (!doBPix_ && (pixelOcc.first == true))
807  continue;
808  if (!doFPix_ && (pixelOcc.second == true))
809  continue;
810 
811  fillTrackHistos(hDA, "sel", &(theTTrack), vertex, beamSpot, fBfield_);
812 
813  // probe checks
814  h_probePt_->Fill(theTrack.pt());
815  h_probePtRebin_->Fill(theTrack.pt());
816  h_probeP_->Fill(theTrack.p());
817  h_probeEta_->Fill(theTrack.eta());
818  h_probePhi_->Fill(theTrack.phi());
819  h2_probeEtaPhi_->Fill(theTrack.eta(), theTrack.phi());
820  h2_probeEtaPt_->Fill(theTrack.eta(), theTrack.pt());
821 
822  h_probeChi2_->Fill(theTrack.chi2());
823  h_probeNormChi2_->Fill(theTrack.normalizedChi2());
824  h_probeCharge_->Fill(theTrack.charge());
825  h_probeQoverP_->Fill(theTrack.qoverp());
826  h_probeHits_->Fill(theTrack.numberOfValidHits());
827  h_probeHits1D_->Fill(nRecHit1D);
828  h_probeHits2D_->Fill(nRecHit2D);
829  h_probeHitsInTIB_->Fill(nhitinTIB);
830  h_probeHitsInTOB_->Fill(nhitinTOB);
831  h_probeHitsInTID_->Fill(nhitinTID);
832  h_probeHitsInTEC_->Fill(nhitinTEC);
833  h_probeHitsInBPIX_->Fill(nhitinBPIX);
834  h_probeHitsInFPIX_->Fill(nhitinFPIX);
835 
836  float dxyRecoV = theTrack.dz(theRecoVertex);
837  float dzRecoV = theTrack.dxy(theRecoVertex);
838  float dxysigmaRecoV =
839  TMath::Sqrt(theTrack.d0Error() * theTrack.d0Error() + xErrOfflineVertex_ * yErrOfflineVertex_);
840  float dzsigmaRecoV =
841  TMath::Sqrt(theTrack.dzError() * theTrack.dzError() + zErrOfflineVertex_ * zErrOfflineVertex_);
842 
843  double zTrack = (theTTrack.stateAtBeamLine().trackStateAtPCA()).position().z();
844  double zVertex = theFittedVertex.position().z();
845  double tantheta = tan((theTTrack.stateAtBeamLine().trackStateAtPCA()).momentum().theta());
846 
847  double dz2 = pow(theTrack.dzError(), 2) + wxy2_ / pow(tantheta, 2);
848  double restrkz = zTrack - zVertex;
849  double pulltrkz = (zTrack - zVertex) / TMath::Sqrt(dz2);
850 
851  h_probedxyRecoV_->Fill(dxyRecoV);
852  h_probedzRecoV_->Fill(dzRecoV);
853 
854  h_probedzRefitV_->Fill(dxyFromMyVertex);
855  h_probedxyRefitV_->Fill(dzFromMyVertex);
856 
857  h_probed0RefitV_->Fill(d0);
858  h_probez0RefitV_->Fill(z0);
859 
860  h_probesignIP2DRefitV_->Fill(s_ip2dpv_corr);
861  h_probed3DRefitV_->Fill(ip3d_corr);
862  h_probereszRefitV_->Fill(restrkz);
863 
864  h_probeRecoVSigZ_->Fill(dzRecoV / dzsigmaRecoV);
865  h_probeRecoVSigXY_->Fill(dxyRecoV / dxysigmaRecoV);
866  h_probeRefitVSigZ_->Fill(dzFromMyVertex / dz_err);
867  h_probeRefitVSigXY_->Fill(dxyFromMyVertex / s_ip2dpv_err);
868  h_probeRefitVSig3D_->Fill(ip3d_corr / ip3d_err);
869  h_probeRefitVLogSig3D_->Fill(log10(ip3d_corr / ip3d_err));
870  h_probeRefitVSigResZ_->Fill(pulltrkz);
871 
872  a_dxyVsPhi->Fill(trackphi, dxyFromMyVertex * cmToum);
873  a_dzVsPhi->Fill(trackphi, z0 * cmToum);
874  n_dxyVsPhi->Fill(trackphi, dxyFromMyVertex / s_ip2dpv_err);
875  n_dzVsPhi->Fill(trackphi, z0 / z0_error);
876 
877  a_dxyVsEta->Fill(tracketa, dxyFromMyVertex * cmToum);
878  a_dzVsEta->Fill(tracketa, z0 * cmToum);
879  n_dxyVsEta->Fill(tracketa, dxyFromMyVertex / s_ip2dpv_err);
880  n_dzVsEta->Fill(tracketa, z0 / z0_error);
881 
882  if (ladder_num > 0 && module_num > 0) {
883  LogDebug("PrimaryVertexValidation")
884  << " ladder_num: " << ladder_num << " module_num: " << module_num << std::endl;
885 
886  PVValHelper::fillByIndex(h_dxy_modZ_, module_num - 1, dxyFromMyVertex * cmToum, "17");
887  PVValHelper::fillByIndex(h_dz_modZ_, module_num - 1, dzFromMyVertex * cmToum, "18");
888  PVValHelper::fillByIndex(h_norm_dxy_modZ_, module_num - 1, dxyFromMyVertex / s_ip2dpv_err, "19");
889  PVValHelper::fillByIndex(h_norm_dz_modZ_, module_num - 1, dzFromMyVertex / dz_err, "20");
890 
891  PVValHelper::fillByIndex(h_dxy_ladder_, ladder_num - 1, dxyFromMyVertex * cmToum, "21");
892 
893  LogDebug("PrimaryVertexValidation") << "h_dxy_ladder size:" << h_dxy_ladder_.size() << std::endl;
894 
895  if (L1BPixHitCount == 1) {
896  PVValHelper::fillByIndex(h_dxy_ladderNoOverlap_, ladder_num - 1, dxyFromMyVertex * cmToum);
897  } else {
898  PVValHelper::fillByIndex(h_dxy_ladderOverlap_, ladder_num - 1, dxyFromMyVertex * cmToum);
899  }
900 
901  h2_probePassingLayer1Map_->Fill(module_num, ladder_num);
902 
903  PVValHelper::fillByIndex(h_dz_ladder_, ladder_num - 1, dzFromMyVertex * cmToum, "22");
904  PVValHelper::fillByIndex(h_norm_dxy_ladder_, ladder_num - 1, dxyFromMyVertex / s_ip2dpv_err, "23");
905  PVValHelper::fillByIndex(h_norm_dz_ladder_, ladder_num - 1, dzFromMyVertex / dz_err, "24");
906  }
907 
908  // filling the binned distributions
909  for (int i = 0; i < nBins_; i++) {
910  float phiF = theDetails_.trendbins[PVValHelper::phi][i];
911  float phiL = theDetails_.trendbins[PVValHelper::phi][i + 1];
912 
913  float etaF = theDetails_.trendbins[PVValHelper::eta][i];
914  float etaL = theDetails_.trendbins[PVValHelper::eta][i + 1];
915 
916  if (tracketa >= etaF && tracketa < etaL) {
917  PVValHelper::fillByIndex(a_dxyEtaResiduals, i, dxyFromMyVertex * cmToum, "25");
920  PVValHelper::fillByIndex(a_dzEtaResiduals, i, dzFromMyVertex * cmToum, "28");
921  PVValHelper::fillByIndex(n_dxyEtaResiduals, i, dxyFromMyVertex / s_ip2dpv_err, "29");
922  PVValHelper::fillByIndex(n_dzEtaResiduals, i, dzFromMyVertex / dz_err, "30");
923  PVValHelper::fillByIndex(a_IP2DEtaResiduals, i, s_ip2dpv_corr * cmToum, "31");
924  PVValHelper::fillByIndex(n_IP2DEtaResiduals, i, s_ip2dpv_corr / s_ip2dpv_err, "32");
927  PVValHelper::fillByIndex(a_d3DEtaResiduals, i, ip3d_corr * cmToum, "35");
928  PVValHelper::fillByIndex(n_d3DEtaResiduals, i, ip3d_corr / ip3d_err, "36");
929  PVValHelper::fillByIndex(a_IP3DEtaResiduals, i, s_ip3dpv_corr * cmToum, "37");
930  PVValHelper::fillByIndex(n_IP3DEtaResiduals, i, s_ip3dpv_corr / s_ip3dpv_err, "38");
931  }
932 
933  if (trackphi >= phiF && trackphi < phiL) {
934  PVValHelper::fillByIndex(a_dxyPhiResiduals, i, dxyFromMyVertex * cmToum, "39");
937  PVValHelper::fillByIndex(a_dzPhiResiduals, i, dzFromMyVertex * cmToum, "42");
938  PVValHelper::fillByIndex(n_dxyPhiResiduals, i, dxyFromMyVertex / s_ip2dpv_err, "43");
939  PVValHelper::fillByIndex(n_dzPhiResiduals, i, dzFromMyVertex / dz_err, "44");
940  PVValHelper::fillByIndex(a_IP2DPhiResiduals, i, s_ip2dpv_corr * cmToum, "45");
941  PVValHelper::fillByIndex(n_IP2DPhiResiduals, i, s_ip2dpv_corr / s_ip2dpv_err, "46");
944  PVValHelper::fillByIndex(a_d3DPhiResiduals, i, ip3d_corr * cmToum, "49");
945  PVValHelper::fillByIndex(n_d3DPhiResiduals, i, ip3d_corr / ip3d_err, "50");
946  PVValHelper::fillByIndex(a_IP3DPhiResiduals, i, s_ip3dpv_corr * cmToum, "51");
947  PVValHelper::fillByIndex(n_IP3DPhiResiduals, i, s_ip3dpv_corr / s_ip3dpv_err, "52");
948 
949  for (int j = 0; j < nBins_; j++) {
950  float etaJ = theDetails_.trendbins[PVValHelper::eta][j];
951  float etaK = theDetails_.trendbins[PVValHelper::eta][j + 1];
952 
953  if (tracketa >= etaJ && tracketa < etaK) {
954  a_dxyResidualsMap[i][j]->Fill(dxyFromMyVertex * cmToum);
955  a_dzResidualsMap[i][j]->Fill(dzFromMyVertex * cmToum);
956  n_dxyResidualsMap[i][j]->Fill(dxyFromMyVertex / s_ip2dpv_err);
957  n_dzResidualsMap[i][j]->Fill(dzFromMyVertex / dz_err);
958  a_d3DResidualsMap[i][j]->Fill(ip3d_corr * cmToum);
959  n_d3DResidualsMap[i][j]->Fill(ip3d_corr / ip3d_err);
960  }
961  }
962  }
963  }
964  }
965 
966  if (debug_) {
967  edm::LogInfo("PrimaryVertexValidation")
968  << " myVertex.x()= " << myVertex.x() << "\n"
969  << " myVertex.y()= " << myVertex.y() << " \n"
970  << " myVertex.z()= " << myVertex.z() << " \n"
971  << " theTrack.dz(myVertex)= " << theTrack.dz(myVertex) << " \n"
972  << " zPCA -myVertex.z() = " << (theTrack.vertex().z() - myVertex.z());
973 
974  } // ends if debug_
975  } // ends if the fitted vertex is Valid
976 
977  //delete theFitter;
978 
979  } catch (cms::Exception& er) {
980  LogTrace("PrimaryVertexValidation") << "caught std::exception " << er.what() << std::endl;
981  }
982 
983  } //ends if theFinalTracks.size() > 2
984 
985  else {
986  if (debug_)
987  edm::LogInfo("PrimaryVertexValidation") << "Not enough tracks to make a vertex. Returns no vertex info";
988  }
989 
990  ++nTracks_;
991  ++nTracksPerClus_;
992 
993  if (debug_)
994  edm::LogInfo("PrimaryVertexValidation") << "Track " << i << " : pT = " << theTrack.pt();
995 
996  } // for loop on tracks
997  } // for loop on track clusters
998 
999  // Fill the TTree if needed
1000 
1001  if (storeNtuple_) {
1002  rootTree_->Fill();
1003  }
1004 }
1005 
1006 // ------------ method called to discriminate 1D from 2D hits ------------
1008  if (hit.dimension() < 2) {
1009  return false; // some (muon...) stuff really has RecHit1D
1010  } else {
1011  const DetId detId(hit.geographicalId());
1012  if (detId.det() == DetId::Tracker) {
1013  if (detId.subdetId() == PixelSubdetector::PixelBarrel || detId.subdetId() == PixelSubdetector::PixelEndcap) {
1014  return true; // pixel is always 2D
1015  } else if (thePhase != PVValHelper::phase2) { // should be SiStrip now
1016  if (dynamic_cast<const SiStripRecHit2D*>(&hit))
1017  return false; // normal hit
1018  else if (dynamic_cast<const SiStripMatchedRecHit2D*>(&hit))
1019  return true; // matched is 2D
1020  else if (dynamic_cast<const ProjectedSiStripRecHit2D*>(&hit))
1021  return false; // crazy hit...
1022  else {
1023  edm::LogError("UnknownType") << "@SUB=PrimaryVertexValidation::isHit2D"
1024  << "Tracker hit not in pixel and neither SiStripRecHit2D nor "
1025  << "SiStripMatchedRecHit2D nor ProjectedSiStripRecHit2D.";
1026  return false;
1027  }
1028  } else {
1029  return false;
1030  }
1031  } else { // not tracker??
1032  edm::LogWarning("DetectorMismatch") << "@SUB=PrimaryVertexValidation::isHit2D"
1033  << "Hit not in tracker with 'official' dimension >=2.";
1034  return true; // dimension() >= 2 so accept that...
1035  }
1036  }
1037  // never reached...
1038 }
1039 
1040 // ------------ method to check the presence of pixel hits ------------
1042  bool hasBPixHits = false;
1043  bool hasFPixHits = false;
1044 
1045  const reco::HitPattern& p = track.hitPattern();
1046  if (p.numberOfValidPixelEndcapHits() != 0) {
1047  hasFPixHits = true;
1048  }
1049  if (p.numberOfValidPixelBarrelHits() != 0) {
1050  hasBPixHits = true;
1051  }
1052 
1053  return std::make_pair(hasBPixHits, hasFPixHits);
1054 }
1055 
1056 // ------------ method to check the presence of pixel hits ------------
1058  using namespace reco;
1059  const HitPattern& p = track.hitPattern();
1060  for (int i = 0; i < p.numberOfAllHits(HitPattern::TRACK_HITS); i++) {
1061  uint32_t pattern = p.getHitPattern(HitPattern::TRACK_HITS, i);
1062  if (p.pixelBarrelHitFilter(pattern) || p.pixelEndcapHitFilter(pattern)) {
1063  if (p.getLayer(pattern) == 1) {
1064  if (p.validHitFilter(pattern)) {
1065  return true;
1066  }
1067  }
1068  }
1069  }
1070  return false;
1071 }
1072 
1073 // ------------ method called once each job before begining the event loop ------------
1075  edm::LogInfo("PrimaryVertexValidation") << "######################################\n"
1076  << "Begin Job \n"
1077  << "######################################";
1078 
1079  // Define TTree for output
1080  Nevt_ = 0;
1081  if (compressionSettings_ > 0) {
1082  fs->file().SetCompressionSettings(compressionSettings_);
1083  }
1084 
1085  // rootFile_ = new TFile(filename_.c_str(),"recreate");
1086  rootTree_ = fs->make<TTree>("tree", "PV Validation tree");
1087 
1088  // Track Paramters
1089 
1090  if (lightNtupleSwitch_) {
1091  rootTree_->Branch("EventNumber", &EventNumber_, "EventNumber/i");
1092  rootTree_->Branch("RunNumber", &RunNumber_, "RunNumber/i");
1093  rootTree_->Branch("LuminosityBlockNumber", &LuminosityBlockNumber_, "LuminosityBlockNumber/i");
1094  rootTree_->Branch("nOfflineVertices", &nOfflineVertices_, "nOfflineVertices/I");
1095  rootTree_->Branch("nTracks", &nTracks_, "nTracks/I");
1096  rootTree_->Branch("phi", &phi_, "phi[nTracks]/D");
1097  rootTree_->Branch("eta", &eta_, "eta[nTracks]/D");
1098  rootTree_->Branch("pt", &pt_, "pt[nTracks]/D");
1099  rootTree_->Branch("dxyFromMyVertex", &dxyFromMyVertex_, "dxyFromMyVertex[nTracks]/D");
1100  rootTree_->Branch("dzFromMyVertex", &dzFromMyVertex_, "dzFromMyVertex[nTracks]/D");
1101  rootTree_->Branch("d3DFromMyVertex", &d3DFromMyVertex_, "d3DFromMyVertex[nTracks]/D");
1102  rootTree_->Branch("IPTsigFromMyVertex", &IPTsigFromMyVertex_, "IPTsigFromMyVertex_[nTracks]/D");
1103  rootTree_->Branch("IPLsigFromMyVertex", &IPLsigFromMyVertex_, "IPLsigFromMyVertex_[nTracks]/D");
1104  rootTree_->Branch("IP3DsigFromMyVertex", &IP3DsigFromMyVertex_, "IP3DsigFromMyVertex_[nTracks]/D");
1105  rootTree_->Branch("hasRecVertex", &hasRecVertex_, "hasRecVertex[nTracks]/I");
1106  rootTree_->Branch("isGoodTrack", &isGoodTrack_, "isGoodTrack[nTracks]/I");
1107  rootTree_->Branch("isHighPurity", &isHighPurity_, "isHighPurity_[nTracks]/I");
1108 
1109  } else {
1110  rootTree_->Branch("nTracks", &nTracks_, "nTracks/I");
1111  rootTree_->Branch("nTracksPerClus", &nTracksPerClus_, "nTracksPerClus/I");
1112  rootTree_->Branch("nClus", &nClus_, "nClus/I");
1113  rootTree_->Branch("xOfflineVertex", &xOfflineVertex_, "xOfflineVertex/D");
1114  rootTree_->Branch("yOfflineVertex", &yOfflineVertex_, "yOfflineVertex/D");
1115  rootTree_->Branch("zOfflineVertex", &zOfflineVertex_, "zOfflineVertex/D");
1116  rootTree_->Branch("BSx0", &BSx0_, "BSx0/D");
1117  rootTree_->Branch("BSy0", &BSy0_, "BSy0/D");
1118  rootTree_->Branch("BSz0", &BSz0_, "BSz0/D");
1119  rootTree_->Branch("Beamsigmaz", &Beamsigmaz_, "Beamsigmaz/D");
1120  rootTree_->Branch("Beamdxdz", &Beamdxdz_, "Beamdxdz/D");
1121  rootTree_->Branch("BeamWidthX", &BeamWidthX_, "BeamWidthX/D");
1122  rootTree_->Branch("BeamWidthY", &BeamWidthY_, "BeamWidthY/D");
1123  rootTree_->Branch("pt", &pt_, "pt[nTracks]/D");
1124  rootTree_->Branch("p", &p_, "p[nTracks]/D");
1125  rootTree_->Branch("nhits", &nhits_, "nhits[nTracks]/I");
1126  rootTree_->Branch("nhits1D", &nhits1D_, "nhits1D[nTracks]/I");
1127  rootTree_->Branch("nhits2D", &nhits2D_, "nhits2D[nTracks]/I");
1128  rootTree_->Branch("nhitsBPIX", &nhitsBPIX_, "nhitsBPIX[nTracks]/I");
1129  rootTree_->Branch("nhitsFPIX", &nhitsFPIX_, "nhitsFPIX[nTracks]/I");
1130  rootTree_->Branch("nhitsTIB", &nhitsTIB_, "nhitsTIB[nTracks]/I");
1131  rootTree_->Branch("nhitsTID", &nhitsTID_, "nhitsTID[nTracks]/I");
1132  rootTree_->Branch("nhitsTOB", &nhitsTOB_, "nhitsTOB[nTracks]/I");
1133  rootTree_->Branch("nhitsTEC", &nhitsTEC_, "nhitsTEC[nTracks]/I");
1134  rootTree_->Branch("eta", &eta_, "eta[nTracks]/D");
1135  rootTree_->Branch("theta", &theta_, "theta[nTracks]/D");
1136  rootTree_->Branch("phi", &phi_, "phi[nTracks]/D");
1137  rootTree_->Branch("chi2", &chi2_, "chi2[nTracks]/D");
1138  rootTree_->Branch("chi2ndof", &chi2ndof_, "chi2ndof[nTracks]/D");
1139  rootTree_->Branch("charge", &charge_, "charge[nTracks]/I");
1140  rootTree_->Branch("qoverp", &qoverp_, "qoverp[nTracks]/D");
1141  rootTree_->Branch("dz", &dz_, "dz[nTracks]/D");
1142  rootTree_->Branch("dxy", &dxy_, "dxy[nTracks]/D");
1143  rootTree_->Branch("dzBs", &dzBs_, "dzBs[nTracks]/D");
1144  rootTree_->Branch("dxyBs", &dxyBs_, "dxyBs[nTracks]/D");
1145  rootTree_->Branch("xPCA", &xPCA_, "xPCA[nTracks]/D");
1146  rootTree_->Branch("yPCA", &yPCA_, "yPCA[nTracks]/D");
1147  rootTree_->Branch("zPCA", &zPCA_, "zPCA[nTracks]/D");
1148  rootTree_->Branch("xUnbiasedVertex", &xUnbiasedVertex_, "xUnbiasedVertex[nTracks]/D");
1149  rootTree_->Branch("yUnbiasedVertex", &yUnbiasedVertex_, "yUnbiasedVertex[nTracks]/D");
1150  rootTree_->Branch("zUnbiasedVertex", &zUnbiasedVertex_, "zUnbiasedVertex[nTracks]/D");
1151  rootTree_->Branch("chi2normUnbiasedVertex", &chi2normUnbiasedVertex_, "chi2normUnbiasedVertex[nTracks]/F");
1152  rootTree_->Branch("chi2UnbiasedVertex", &chi2UnbiasedVertex_, "chi2UnbiasedVertex[nTracks]/F");
1153  rootTree_->Branch("DOFUnbiasedVertex", &DOFUnbiasedVertex_, " DOFUnbiasedVertex[nTracks]/F");
1154  rootTree_->Branch("chi2ProbUnbiasedVertex", &chi2ProbUnbiasedVertex_, "chi2ProbUnbiasedVertex[nTracks]/F");
1155  rootTree_->Branch(
1156  "sumOfWeightsUnbiasedVertex", &sumOfWeightsUnbiasedVertex_, "sumOfWeightsUnbiasedVertex[nTracks]/F");
1157  rootTree_->Branch("tracksUsedForVertexing", &tracksUsedForVertexing_, "tracksUsedForVertexing[nTracks]/I");
1158  rootTree_->Branch("dxyFromMyVertex", &dxyFromMyVertex_, "dxyFromMyVertex[nTracks]/D");
1159  rootTree_->Branch("dzFromMyVertex", &dzFromMyVertex_, "dzFromMyVertex[nTracks]/D");
1160  rootTree_->Branch("dxyErrorFromMyVertex", &dxyErrorFromMyVertex_, "dxyErrorFromMyVertex_[nTracks]/D");
1161  rootTree_->Branch("dzErrorFromMyVertex", &dzErrorFromMyVertex_, "dzErrorFromMyVertex_[nTracks]/D");
1162  rootTree_->Branch("IPTsigFromMyVertex", &IPTsigFromMyVertex_, "IPTsigFromMyVertex_[nTracks]/D");
1163  rootTree_->Branch("IPLsigFromMyVertex", &IPLsigFromMyVertex_, "IPLsigFromMyVertex_[nTracks]/D");
1164  rootTree_->Branch("hasRecVertex", &hasRecVertex_, "hasRecVertex[nTracks]/I");
1165  rootTree_->Branch("isGoodTrack", &isGoodTrack_, "isGoodTrack[nTracks]/I");
1166  }
1167 
1168  // event histograms
1169  TFileDirectory EventFeatures = fs->mkdir("EventFeatures");
1170 
1171  TH1F::SetDefaultSumw2(kTRUE);
1172 
1174  EventFeatures.make<TH1F>("h_lumiFromConfig", "luminosity from config;;luminosity of present run", 1, -0.5, 0.5);
1175  h_lumiFromConfig->SetBinContent(1, intLumi_);
1176 
1177  h_runFromConfig = EventFeatures.make<TH1I>("h_runFromConfig",
1178  "run number from config;;run number (from configuration)",
1179  runControlNumbers_.size(),
1180  0.,
1181  runControlNumbers_.size());
1182 
1183  for (const auto& run : runControlNumbers_ | boost::adaptors::indexed(1)) {
1184  h_runFromConfig->SetBinContent(run.index(), run.value());
1185  }
1186 
1187  h_runFromEvent =
1188  EventFeatures.make<TH1I>("h_runFromEvent", "run number from event;;run number (from event)", 1, -0.5, 0.5);
1189  h_nTracks =
1190  EventFeatures.make<TH1F>("h_nTracks", "number of tracks per event;n_{tracks}/event;n_{events}", 300, -0.5, 299.5);
1191  h_nClus =
1192  EventFeatures.make<TH1F>("h_nClus", "number of track clusters;n_{clusters}/event;n_{events}", 50, -0.5, 49.5);
1193  h_nOfflineVertices = EventFeatures.make<TH1F>(
1194  "h_nOfflineVertices", "number of offline reconstructed vertices;n_{vertices}/event;n_{events}", 50, -0.5, 49.5);
1195  h_runNumber = EventFeatures.make<TH1F>("h_runNumber", "run number;run number;n_{events}", 100000, 250000., 350000.);
1196  h_xOfflineVertex = EventFeatures.make<TH1F>(
1197  "h_xOfflineVertex", "x-coordinate of offline vertex;x_{vertex};n_{events}", 100, -0.1, 0.1);
1198  h_yOfflineVertex = EventFeatures.make<TH1F>(
1199  "h_yOfflineVertex", "y-coordinate of offline vertex;y_{vertex};n_{events}", 100, -0.1, 0.1);
1200  h_zOfflineVertex = EventFeatures.make<TH1F>(
1201  "h_zOfflineVertex", "z-coordinate of offline vertex;z_{vertex};n_{events}", 100, -30., 30.);
1202  h_xErrOfflineVertex = EventFeatures.make<TH1F>(
1203  "h_xErrOfflineVertex", "x-coordinate error of offline vertex;err_{x}^{vtx};n_{events}", 100, 0., 0.01);
1204  h_yErrOfflineVertex = EventFeatures.make<TH1F>(
1205  "h_yErrOfflineVertex", "y-coordinate error of offline vertex;err_{y}^{vtx};n_{events}", 100, 0., 0.01);
1206  h_zErrOfflineVertex = EventFeatures.make<TH1F>(
1207  "h_zErrOfflineVertex", "z-coordinate error of offline vertex;err_{z}^{vtx};n_{events}", 100, 0., 10.);
1208  h_BSx0 = EventFeatures.make<TH1F>("h_BSx0", "x-coordinate of reco beamspot;x^{BS}_{0};n_{events}", 100, -0.1, 0.1);
1209  h_BSy0 = EventFeatures.make<TH1F>("h_BSy0", "y-coordinate of reco beamspot;y^{BS}_{0};n_{events}", 100, -0.1, 0.1);
1210  h_BSz0 = EventFeatures.make<TH1F>("h_BSz0", "z-coordinate of reco beamspot;z^{BS}_{0};n_{events}", 100, -1., 1.);
1211  h_Beamsigmaz =
1212  EventFeatures.make<TH1F>("h_Beamsigmaz", "z-coordinate beam width;#sigma_{Z}^{beam};n_{events}", 100, 0., 1.);
1213  h_BeamWidthX =
1214  EventFeatures.make<TH1F>("h_BeamWidthX", "x-coordinate beam width;#sigma_{X}^{beam};n_{events}", 100, 0., 0.01);
1215  h_BeamWidthY =
1216  EventFeatures.make<TH1F>("h_BeamWidthY", "y-coordinate beam width;#sigma_{Y}^{beam};n_{events}", 100, 0., 0.01);
1217 
1218  h_etaMax = EventFeatures.make<TH1F>("etaMax", "etaMax", 1, -0.5, 0.5);
1219  h_pTinfo = EventFeatures.make<TH1F>("pTinfo", "pTinfo", 3, -1.5, 1.5);
1220  h_pTinfo->GetXaxis()->SetBinLabel(1, "n. bins");
1221  h_pTinfo->GetXaxis()->SetBinLabel(2, "pT min");
1222  h_pTinfo->GetXaxis()->SetBinLabel(3, "pT max");
1223 
1224  h_nbins = EventFeatures.make<TH1F>("nbins", "nbins", 1, -0.5, 0.5);
1225  h_nLadders = EventFeatures.make<TH1F>("nladders", "n. ladders", 1, -0.5, 0.5);
1226  h_nModZ = EventFeatures.make<TH1F>("nModZ", "n. modules along z", 1, -0.5, 0.5);
1227 
1228  // probe track histograms
1229  TFileDirectory ProbeFeatures = fs->mkdir("ProbeTrackFeatures");
1230 
1231  h_probePt_ = ProbeFeatures.make<TH1F>("h_probePt", "p_{T} of probe track;track p_{T} (GeV); tracks", 100, 0., 50.);
1232  h_probePtRebin_ = ProbeFeatures.make<TH1F>(
1233  "h_probePtRebin", "p_{T} of probe track;track p_{T} (GeV); tracks", mypT_bins_.size() - 1, mypT_bins_.data());
1234  h_probeP_ = ProbeFeatures.make<TH1F>("h_probeP", "momentum of probe track;track p (GeV); tracks", 100, 0., 100.);
1235  h_probeEta_ = ProbeFeatures.make<TH1F>(
1236  "h_probeEta", "#eta of the probe track;track #eta;tracks", 54, -etaOfProbe_, etaOfProbe_);
1237  h_probePhi_ = ProbeFeatures.make<TH1F>("h_probePhi", "#phi of probe track;track #phi (rad);tracks", 100, -3.15, 3.15);
1238 
1239  h2_probeEtaPhi_ =
1240  ProbeFeatures.make<TH2F>("h2_probeEtaPhi",
1241  "probe track #phi vs #eta;#eta of probe track;track #phi of probe track (rad); tracks",
1242  54,
1243  -etaOfProbe_,
1244  etaOfProbe_,
1245  100,
1246  -M_PI,
1247  M_PI);
1248  h2_probeEtaPt_ = ProbeFeatures.make<TH2F>("h2_probeEtaPt",
1249  "probe track p_{T} vs #eta;#eta of probe track;track p_{T} (GeV); tracks",
1250  54,
1251  -etaOfProbe_,
1252  etaOfProbe_,
1253  100,
1254  0.,
1255  50.);
1256 
1257  h_probeChi2_ =
1258  ProbeFeatures.make<TH1F>("h_probeChi2", "#chi^{2} of probe track;track #chi^{2}; tracks", 100, 0., 100.);
1259  h_probeNormChi2_ = ProbeFeatures.make<TH1F>(
1260  "h_probeNormChi2", " normalized #chi^{2} of probe track;track #chi^{2}/ndof; tracks", 100, 0., 10.);
1261  h_probeCharge_ =
1262  ProbeFeatures.make<TH1F>("h_probeCharge", "charge of probe track;track charge Q;tracks", 3, -1.5, 1.5);
1263  h_probeQoverP_ =
1264  ProbeFeatures.make<TH1F>("h_probeQoverP", "q/p of probe track; track Q/p (GeV^{-1});tracks", 200, -1., 1.);
1265  h_probedzRecoV_ = ProbeFeatures.make<TH1F>(
1266  "h_probedzRecoV", "d_{z}(V_{offline}) of probe track;track d_{z}(V_{off}) (cm);tracks", 200, -1., 1.);
1267  h_probedxyRecoV_ = ProbeFeatures.make<TH1F>(
1268  "h_probedxyRecoV", "d_{xy}(V_{offline}) of probe track;track d_{xy}(V_{off}) (cm);tracks", 200, -1., 1.);
1269  h_probedzRefitV_ = ProbeFeatures.make<TH1F>(
1270  "h_probedzRefitV", "d_{z}(V_{refit}) of probe track;track d_{z}(V_{fit}) (cm);tracks", 200, -0.5, 0.5);
1271  h_probesignIP2DRefitV_ = ProbeFeatures.make<TH1F>(
1272  "h_probesignIPRefitV", "ip_{2D}(V_{refit}) of probe track;track ip_{2D}(V_{fit}) (cm);tracks", 200, -1., 1.);
1273  h_probedxyRefitV_ = ProbeFeatures.make<TH1F>(
1274  "h_probedxyRefitV", "d_{xy}(V_{refit}) of probe track;track d_{xy}(V_{fit}) (cm);tracks", 200, -0.5, 0.5);
1275 
1276  h_probez0RefitV_ = ProbeFeatures.make<TH1F>(
1277  "h_probez0RefitV", "z_{0}(V_{refit}) of probe track;track z_{0}(V_{fit}) (cm);tracks", 200, -1., 1.);
1278  h_probed0RefitV_ = ProbeFeatures.make<TH1F>(
1279  "h_probed0RefitV", "d_{0}(V_{refit}) of probe track;track d_{0}(V_{fit}) (cm);tracks", 200, -1., 1.);
1280 
1281  h_probed3DRefitV_ = ProbeFeatures.make<TH1F>(
1282  "h_probed3DRefitV", "d_{3D}(V_{refit}) of probe track;track d_{3D}(V_{fit}) (cm);tracks", 200, 0., 1.);
1283  h_probereszRefitV_ = ProbeFeatures.make<TH1F>(
1284  "h_probeReszRefitV", "z_{track} -z_{V_{refit}};track res_{z}(V_{refit}) (cm);tracks", 200, -1., 1.);
1285 
1286  h_probeRecoVSigZ_ = ProbeFeatures.make<TH1F>(
1287  "h_probeRecoVSigZ", "Longitudinal DCA Significance (reco);d_{z}(V_{off})/#sigma_{dz};tracks", 100, -8, 8);
1288  h_probeRecoVSigXY_ = ProbeFeatures.make<TH1F>(
1289  "h_probeRecoVSigXY", "Transverse DCA Significance (reco);d_{xy}(V_{off})/#sigma_{dxy};tracks", 100, -8, 8);
1290  h_probeRefitVSigZ_ = ProbeFeatures.make<TH1F>(
1291  "h_probeRefitVSigZ", "Longitudinal DCA Significance (refit);d_{z}(V_{fit})/#sigma_{dz};tracks", 100, -8, 8);
1292  h_probeRefitVSigXY_ = ProbeFeatures.make<TH1F>(
1293  "h_probeRefitVSigXY", "Transverse DCA Significance (refit);d_{xy}(V_{fit})/#sigma_{dxy};tracks", 100, -8, 8);
1294  h_probeRefitVSig3D_ = ProbeFeatures.make<TH1F>(
1295  "h_probeRefitVSig3D", "3D DCA Significance (refit);d_{3D}/#sigma_{3D};tracks", 100, 0., 20.);
1297  ProbeFeatures.make<TH1F>("h_probeRefitVLogSig3D",
1298  "log_{10}(3D DCA-Significance) (refit);log_{10}(d_{3D}/#sigma_{3D});tracks",
1299  100,
1300  -5.,
1301  4.);
1302  h_probeRefitVSigResZ_ = ProbeFeatures.make<TH1F>(
1303  "h_probeRefitVSigResZ",
1304  "Longitudinal residual significance (refit);(z_{track} -z_{V_{fit}})/#sigma_{res_{z}};tracks",
1305  100,
1306  -8,
1307  8);
1308 
1309  h_probeHits_ = ProbeFeatures.make<TH1F>("h_probeNRechits", "N_{hits} ;N_{hits} ;tracks", 40, -0.5, 39.5);
1310  h_probeHits1D_ = ProbeFeatures.make<TH1F>("h_probeNRechits1D", "N_{hits} 1D ;N_{hits} 1D ;tracks", 40, -0.5, 39.5);
1311  h_probeHits2D_ = ProbeFeatures.make<TH1F>("h_probeNRechits2D", "N_{hits} 2D ;N_{hits} 2D ;tracks", 40, -0.5, 39.5);
1313  ProbeFeatures.make<TH1F>("h_probeNRechitsTIB", "N_{hits} TIB ;N_{hits} TIB;tracks", 40, -0.5, 39.5);
1315  ProbeFeatures.make<TH1F>("h_probeNRechitsTOB", "N_{hits} TOB ;N_{hits} TOB;tracks", 40, -0.5, 39.5);
1317  ProbeFeatures.make<TH1F>("h_probeNRechitsTID", "N_{hits} TID ;N_{hits} TID;tracks", 40, -0.5, 39.5);
1319  ProbeFeatures.make<TH1F>("h_probeNRechitsTEC", "N_{hits} TEC ;N_{hits} TEC;tracks", 40, -0.5, 39.5);
1321  ProbeFeatures.make<TH1F>("h_probeNRechitsBPIX", "N_{hits} BPIX;N_{hits} BPIX;tracks", 40, -0.5, 39.5);
1323  ProbeFeatures.make<TH1F>("h_probeNRechitsFPIX", "N_{hits} FPIX;N_{hits} FPIX;tracks", 40, -0.5, 39.5);
1324 
1325  h_probeL1Ladder_ = ProbeFeatures.make<TH1F>(
1326  "h_probeL1Ladder", "Ladder number (L1 hit); ladder number", nLadders_ + 2, -1.5, nLadders_ + 0.5);
1327  h_probeL1Module_ = ProbeFeatures.make<TH1F>(
1328  "h_probeL1Module", "Module number (L1 hit); module number", nModZ_ + 2, -1.5, nModZ_ + 0.5);
1329 
1330  h2_probeLayer1Map_ = ProbeFeatures.make<TH2F>("h2_probeLayer1Map",
1331  "Position in Layer 1 of first hit;module number;ladder number",
1332  nModZ_,
1333  0.5,
1334  nModZ_ + 0.5,
1335  nLadders_,
1336  0.5,
1337  nLadders_ + 0.5);
1338 
1339  h2_probePassingLayer1Map_ = ProbeFeatures.make<TH2F>("h2_probePassingLayer1Map",
1340  "Position in Layer 1 of first hit;module number;ladder number",
1341  nModZ_,
1342  0.5,
1343  nModZ_ + 0.5,
1344  nLadders_,
1345  0.5,
1346  nLadders_ + 0.5);
1348  ProbeFeatures.make<TH1I>("h_probeHasBPixL1Overlap", "n. hits in L1;n. L1-BPix hits;tracks", 5, -0.5, 4.5);
1349  h_probeL1ClusterProb_ = ProbeFeatures.make<TH1F>(
1350  "h_probeL1ClusterProb",
1351  "log_{10}(Cluster Probability) for Layer1 hits;log_{10}(cluster probability); n. Layer1 hits",
1352  100,
1353  -10.,
1354  0.);
1355 
1356  // refit vertex features
1357  TFileDirectory RefitVertexFeatures = fs->mkdir("RefitVertexFeatures");
1358  h_fitVtxNtracks_ = RefitVertexFeatures.make<TH1F>(
1359  "h_fitVtxNtracks", "N_{trks} used in vertex fit;N^{fit}_{tracks};vertices", 100, -0.5, 99.5);
1360  h_fitVtxNdof_ = RefitVertexFeatures.make<TH1F>(
1361  "h_fitVtxNdof", "N_{DOF} of vertex fit;N_{DOF} of refit vertex;vertices", 100, -0.5, 99.5);
1362  h_fitVtxChi2_ = RefitVertexFeatures.make<TH1F>(
1363  "h_fitVtxChi2", "#chi^{2} of vertex fit;vertex #chi^{2};vertices", 100, -0.5, 99.5);
1364  h_fitVtxChi2ndf_ = RefitVertexFeatures.make<TH1F>(
1365  "h_fitVtxChi2ndf", "#chi^{2}/ndf of vertex fit;vertex #chi^{2}/ndf;vertices", 100, -0.5, 9.5);
1366  h_fitVtxChi2Prob_ = RefitVertexFeatures.make<TH1F>(
1367  "h_fitVtxChi2Prob", "Prob(#chi^{2},ndf) of vertex fit;Prob(#chi^{2},ndf);vertices", 40, 0., 1.);
1368  h_fitVtxTrackWeights_ = RefitVertexFeatures.make<TH1F>(
1369  "h_fitVtxTrackWeights", "track weights associated to track;track weights;tracks", 40, 0., 1.);
1370  h_fitVtxTrackAverageWeight_ = RefitVertexFeatures.make<TH1F>(
1371  "h_fitVtxTrackAverageWeight_", "average track weight per vertex;#LT track weight #GT;vertices", 40, 0., 1.);
1372 
1373  if (useTracksFromRecoVtx_) {
1374  TFileDirectory RecoVertexFeatures = fs->mkdir("RecoVertexFeatures");
1376  RecoVertexFeatures.make<TH1F>("h_recoVtxNtracks", "N^{vtx}_{trks};N^{vtx}_{trks};vertices", 100, -0.5, 99.5);
1378  RecoVertexFeatures.make<TH1F>("h_recoVtxChi2ndf", "#chi^{2}/ndf vtx;#chi^{2}/ndf vtx;vertices", 10, -0.5, 9.5);
1379  h_recoVtxChi2Prob_ = RecoVertexFeatures.make<TH1F>(
1380  "h_recoVtxChi2Prob", "Prob(#chi^{2},ndf);Prob(#chi^{2},ndf);vertices", 40, 0., 1.);
1381  h_recoVtxSumPt_ =
1382  RecoVertexFeatures.make<TH1F>("h_recoVtxSumPt", "Sum(p^{trks}_{T});Sum(p^{trks}_{T});vertices", 100, 0., 200.);
1383  }
1384 
1385  TFileDirectory DA = fs->mkdir("DA");
1386  //DA.cd();
1387  hDA = bookVertexHistograms(DA);
1388  //for(std::map<std::string,TH1*>::const_iterator hist=hDA.begin(); hist!=hDA.end(); hist++){
1389  //hist->second->SetDirectory(DA);
1390  // DA.make<TH1F>(hist->second);
1391  // }
1392 
1393  // initialize the residuals histograms
1394 
1395  const float dxymax_phi = theDetails_.getHigh(PVValHelper::dxy, PVValHelper::phi);
1396  const float dzmax_phi = theDetails_.getHigh(PVValHelper::dz, PVValHelper::eta);
1397  const float dxymax_eta = theDetails_.getHigh(PVValHelper::dxy, PVValHelper::phi);
1398  const float dzmax_eta = theDetails_.getHigh(PVValHelper::dz, PVValHelper::eta);
1399  //const float d3Dmax_phi = theDetails_.getHigh(PVValHelper::d3D,PVValHelper::phi);
1400  const float d3Dmax_eta = theDetails_.getHigh(PVValHelper::d3D, PVValHelper::eta);
1401 
1403  //
1404  // Unbiased track-to-vertex residuals
1405  // The vertex is refit without the probe track
1406  //
1408 
1409  // _ _ _ _ ___ _ _ _
1410  // /_\ | |__ ___ ___| |_ _| |_ ___ | _ \___ __(_)__| |_ _ __ _| |___
1411  // / _ \| '_ (_-</ _ \ | || | _/ -_) | / -_|_-< / _` | || / _` | (_-<
1412  // /_/ \_\_.__/__/\___/_|\_,_|\__\___| |_|_\___/__/_\__,_|\_,_\__,_|_/__/
1413  //
1414 
1415  TFileDirectory AbsTransPhiRes = fs->mkdir("Abs_Transv_Phi_Residuals");
1420 
1421  TFileDirectory AbsTransEtaRes = fs->mkdir("Abs_Transv_Eta_Residuals");
1426 
1427  TFileDirectory AbsLongPhiRes = fs->mkdir("Abs_Long_Phi_Residuals");
1430 
1431  TFileDirectory AbsLongEtaRes = fs->mkdir("Abs_Long_Eta_Residuals");
1434 
1435  TFileDirectory Abs3DPhiRes = fs->mkdir("Abs_3D_Phi_Residuals");
1438 
1439  TFileDirectory Abs3DEtaRes = fs->mkdir("Abs_3D_Eta_Residuals");
1442 
1443  TFileDirectory NormTransPhiRes = fs->mkdir("Norm_Transv_Phi_Residuals");
1446 
1447  TFileDirectory NormTransEtaRes = fs->mkdir("Norm_Transv_Eta_Residuals");
1450 
1451  TFileDirectory NormLongPhiRes = fs->mkdir("Norm_Long_Phi_Residuals");
1454 
1455  TFileDirectory NormLongEtaRes = fs->mkdir("Norm_Long_Eta_Residuals");
1458 
1459  TFileDirectory Norm3DPhiRes = fs->mkdir("Norm_3D_Phi_Residuals");
1462 
1463  TFileDirectory Norm3DEtaRes = fs->mkdir("Norm_3D_Eta_Residuals");
1466 
1467  TFileDirectory AbsDoubleDiffRes = fs->mkdir("Abs_DoubleDiffResiduals");
1468  TFileDirectory NormDoubleDiffRes = fs->mkdir("Norm_DoubleDiffResiduals");
1469 
1470  TFileDirectory AbsL1Map = fs->mkdir("Abs_L1Residuals");
1471  TFileDirectory NormL1Map = fs->mkdir("Norm_L1Residuals");
1472 
1473  // book residuals vs pT histograms
1474 
1475  TFileDirectory AbsTranspTRes = fs->mkdir("Abs_Transv_pT_Residuals");
1477 
1478  TFileDirectory AbsLongpTRes = fs->mkdir("Abs_Long_pT_Residuals");
1480 
1481  TFileDirectory NormTranspTRes = fs->mkdir("Norm_Transv_pT_Residuals");
1483 
1484  TFileDirectory NormLongpTRes = fs->mkdir("Norm_Long_pT_Residuals");
1486 
1487  // book residuals vs pT histograms in central region (|eta|<1.0)
1488 
1489  TFileDirectory AbsTranspTCentralRes = fs->mkdir("Abs_Transv_pTCentral_Residuals");
1491 
1492  TFileDirectory AbsLongpTCentralRes = fs->mkdir("Abs_Long_pTCentral_Residuals");
1494 
1495  TFileDirectory NormTranspTCentralRes = fs->mkdir("Norm_Transv_pTCentral_Residuals");
1498 
1499  TFileDirectory NormLongpTCentralRes = fs->mkdir("Norm_Long_pTCentral_Residuals");
1502 
1503  // book residuals vs module number
1504 
1505  TFileDirectory AbsTransModZRes = fs->mkdir("Abs_Transv_modZ_Residuals");
1507 
1508  TFileDirectory AbsLongModZRes = fs->mkdir("Abs_Long_modZ_Residuals");
1510 
1511  // _ _ _ _ _ ___ _ _ _
1512  // | \| |___ _ _ _ __ __ _| (_)______ __| | | _ \___ __(_)__| |_ _ __ _| |___
1513  // | .` / _ \ '_| ' \/ _` | | |_ / -_) _` | | / -_|_-< / _` | || / _` | (_-<
1514  // |_|\_\___/_| |_|_|_\__,_|_|_/__\___\__,_| |_|_\___/__/_\__,_|\_,_\__,_|_/__/
1515  //
1516 
1517  TFileDirectory NormTransModZRes = fs->mkdir("Norm_Transv_modZ_Residuals");
1519 
1520  TFileDirectory NormLongModZRes = fs->mkdir("Norm_Long_modZ_Residuals");
1522 
1523  TFileDirectory AbsTransLadderRes = fs->mkdir("Abs_Transv_ladder_Residuals");
1525 
1526  TFileDirectory AbsTransLadderResOverlap = fs->mkdir("Abs_Transv_ladderOverlap_Residuals");
1529 
1530  TFileDirectory AbsTransLadderResNoOverlap = fs->mkdir("Abs_Transv_ladderNoOverlap_Residuals");
1533 
1534  TFileDirectory AbsLongLadderRes = fs->mkdir("Abs_Long_ladder_Residuals");
1536 
1537  TFileDirectory NormTransLadderRes = fs->mkdir("Norm_Transv_ladder_Residuals");
1540 
1541  TFileDirectory NormLongLadderRes = fs->mkdir("Norm_Long_ladder_Residuals");
1544 
1545  // book residuals as function of nLadders and nModules
1546 
1547  for (unsigned int iLadder = 0; iLadder < nLadders_; iLadder++) {
1548  for (unsigned int iModule = 0; iModule < nModZ_; iModule++) {
1549  a_dxyL1ResidualsMap[iLadder][iModule] =
1550  AbsL1Map.make<TH1F>(Form("histo_dxy_ladder%i_module%i", iLadder, iModule),
1551  Form("d_{xy} ladder=%i module=%i;d_{xy} [#mum];tracks", iLadder, iModule),
1553  -dzmax_eta,
1554  dzmax_eta);
1555 
1556  a_dzL1ResidualsMap[iLadder][iModule] =
1557  AbsL1Map.make<TH1F>(Form("histo_dz_ladder%i_module%i", iLadder, iModule),
1558  Form("d_{z} ladder=%i module=%i;d_{z} [#mum];tracks", iLadder, iModule),
1560  -dzmax_eta,
1561  dzmax_eta);
1562 
1563  n_dxyL1ResidualsMap[iLadder][iModule] =
1564  NormL1Map.make<TH1F>(Form("histo_norm_dxy_ladder%i_module%i", iLadder, iModule),
1565  Form("d_{xy} ladder=%i module=%i;d_{xy}/#sigma_{d_{xy}};tracks", iLadder, iModule),
1567  -dzmax_eta / 100,
1568  dzmax_eta / 100);
1569 
1570  n_dzL1ResidualsMap[iLadder][iModule] =
1571  NormL1Map.make<TH1F>(Form("histo_norm_dz_ladder%i_module%i", iLadder, iModule),
1572  Form("d_{z} ladder=%i module=%i;d_{z}/#sigma_{d_{z}};tracks", iLadder, iModule),
1574  -dzmax_eta / 100,
1575  dzmax_eta / 100);
1576  }
1577  }
1578 
1579  // book residuals as function of phi and eta
1580 
1581  for (int i = 0; i < nBins_; ++i) {
1582  float phiF = theDetails_.trendbins[PVValHelper::phi][i];
1583  float phiL = theDetails_.trendbins[PVValHelper::phi][i + 1];
1584 
1585  // ___ _ _ ___ _ __ __ ___ _ _ _
1586  // | \ ___ _ _| |__| |___| \(_)/ _|/ _| | _ \___ __(_)__| |_ _ __ _| |___
1587  // | |) / _ \ || | '_ \ / -_) |) | | _| _| | / -_|_-< / _` | || / _` | (_-<
1588  // |___/\___/\_,_|_.__/_\___|___/|_|_| |_| |_|_\___/__/_\__,_|\_,_\__,_|_/__/
1589 
1590  for (int j = 0; j < nBins_; ++j) {
1591  float etaF = theDetails_.trendbins[PVValHelper::eta][j];
1592  float etaL = theDetails_.trendbins[PVValHelper::eta][j + 1];
1593 
1594  a_dxyResidualsMap[i][j] = AbsDoubleDiffRes.make<TH1F>(
1595  Form("histo_dxy_eta_plot%i_phi_plot%i", i, j),
1596  Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{xy};tracks", etaF, etaL, phiF, phiL),
1598  -dzmax_eta,
1599  dzmax_eta);
1600 
1601  a_dzResidualsMap[i][j] = AbsDoubleDiffRes.make<TH1F>(
1602  Form("histo_dz_eta_plot%i_phi_plot%i", i, j),
1603  Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{z};tracks", etaF, etaL, phiF, phiL),
1605  -dzmax_eta,
1606  dzmax_eta);
1607 
1608  a_d3DResidualsMap[i][j] = AbsDoubleDiffRes.make<TH1F>(
1609  Form("histo_d3D_eta_plot%i_phi_plot%i", i, j),
1610  Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{3D};tracks", etaF, etaL, phiF, phiL),
1612  0.,
1613  d3Dmax_eta);
1614 
1615  n_dxyResidualsMap[i][j] = NormDoubleDiffRes.make<TH1F>(
1616  Form("histo_norm_dxy_eta_plot%i_phi_plot%i", i, j),
1617  Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{xy}/#sigma_{d_{xy}};tracks",
1618  etaF,
1619  etaL,
1620  phiF,
1621  phiL),
1623  -dzmax_eta / 100,
1624  dzmax_eta / 100);
1625 
1626  n_dzResidualsMap[i][j] = NormDoubleDiffRes.make<TH1F>(
1627  Form("histo_norm_dz_eta_plot%i_phi_plot%i", i, j),
1628  Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{z}/#sigma_{d_{z}};tracks",
1629  etaF,
1630  etaL,
1631  phiF,
1632  phiL),
1634  -dzmax_eta / 100,
1635  dzmax_eta / 100);
1636 
1637  n_d3DResidualsMap[i][j] = NormDoubleDiffRes.make<TH1F>(
1638  Form("histo_norm_d3D_eta_plot%i_phi_plot%i", i, j),
1639  Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{3D}/#sigma_{d_{3D}};tracks",
1640  etaF,
1641  etaL,
1642  phiF,
1643  phiL),
1645  0.,
1646  d3Dmax_eta);
1647  }
1648  }
1649 
1650  // declaration of the directories
1651 
1652  TFileDirectory BiasVsParameter = fs->mkdir("BiasVsParameter");
1653 
1654  a_dxyVsPhi = BiasVsParameter.make<TH2F>("h2_dxy_vs_phi",
1655  "d_{xy} vs track #phi;track #phi [rad];track d_{xy}(PV) [#mum]",
1656  nBins_,
1657  -M_PI,
1658  M_PI,
1660  -dxymax_phi,
1661  dxymax_phi);
1662 
1663  a_dzVsPhi = BiasVsParameter.make<TH2F>("h2_dz_vs_phi",
1664  "d_{z} vs track #phi;track #phi [rad];track d_{z}(PV) [#mum]",
1665  nBins_,
1666  -M_PI,
1667  M_PI,
1669  -dzmax_phi,
1670  dzmax_phi);
1671 
1672  n_dxyVsPhi = BiasVsParameter.make<TH2F>(
1673  "h2_n_dxy_vs_phi",
1674  "d_{xy}/#sigma_{d_{xy}} vs track #phi;track #phi [rad];track d_{xy}(PV)/#sigma_{d_{xy}}",
1675  nBins_,
1676  -M_PI,
1677  M_PI,
1679  -dxymax_phi / 100.,
1680  dxymax_phi / 100.);
1681 
1682  n_dzVsPhi =
1683  BiasVsParameter.make<TH2F>("h2_n_dz_vs_phi",
1684  "d_{z}/#sigma_{d_{z}} vs track #phi;track #phi [rad];track d_{z}(PV)/#sigma_{d_{z}}",
1685  nBins_,
1686  -M_PI,
1687  M_PI,
1689  -dzmax_phi / 100.,
1690  dzmax_phi / 100.);
1691 
1692  a_dxyVsEta = BiasVsParameter.make<TH2F>("h2_dxy_vs_eta",
1693  "d_{xy} vs track #eta;track #eta;track d_{xy}(PV) [#mum]",
1694  nBins_,
1695  -etaOfProbe_,
1696  etaOfProbe_,
1698  -dxymax_eta,
1699  dzmax_eta);
1700 
1701  a_dzVsEta = BiasVsParameter.make<TH2F>("h2_dz_vs_eta",
1702  "d_{z} vs track #eta;track #eta;track d_{z}(PV) [#mum]",
1703  nBins_,
1704  -etaOfProbe_,
1705  etaOfProbe_,
1707  -dzmax_eta,
1708  dzmax_eta);
1709 
1710  n_dxyVsEta =
1711  BiasVsParameter.make<TH2F>("h2_n_dxy_vs_eta",
1712  "d_{xy}/#sigma_{d_{xy}} vs track #eta;track #eta;track d_{xy}(PV)/#sigma_{d_{xy}}",
1713  nBins_,
1714  -etaOfProbe_,
1715  etaOfProbe_,
1717  -dxymax_eta / 100.,
1718  dxymax_eta / 100.);
1719 
1720  n_dzVsEta = BiasVsParameter.make<TH2F>("h2_n_dz_vs_eta",
1721  "d_{z}/#sigma_{d_{z}} vs track #eta;track #eta;track d_{z}(PV)/#sigma_{d_{z}}",
1722  nBins_,
1723  -etaOfProbe_,
1724  etaOfProbe_,
1726  -dzmax_eta / 100.,
1727  dzmax_eta / 100.);
1728 
1729  MeanTrendsDir = fs->mkdir("MeanTrends");
1730  WidthTrendsDir = fs->mkdir("WidthTrends");
1731  MedianTrendsDir = fs->mkdir("MedianTrends");
1732  MADTrendsDir = fs->mkdir("MADTrends");
1733 
1734  Mean2DMapsDir = fs->mkdir("MeanMaps");
1735  Width2DMapsDir = fs->mkdir("WidthMaps");
1736 
1737  double highedge = nBins_ - 0.5;
1738  double lowedge = -0.5;
1739 
1740  // means and widths from the fit
1741 
1743  MeanTrendsDir.make<TH1F>("means_dxy_phi",
1744  "#LT d_{xy} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{xy} #GT [#mum]",
1745  nBins_,
1746  lowedge,
1747  highedge);
1748 
1750  WidthTrendsDir.make<TH1F>("widths_dxy_phi",
1751  "#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees];#sigma_{d_{xy}} [#mum]",
1752  nBins_,
1753  lowedge,
1754  highedge);
1755 
1757  MeanTrendsDir.make<TH1F>("means_dz_phi",
1758  "#LT d_{z} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{z} #GT [#mum]",
1759  nBins_,
1760  lowedge,
1761  highedge);
1762 
1764  WidthTrendsDir.make<TH1F>("widths_dz_phi",
1765  "#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];#sigma_{d_{z}} [#mum]",
1766  nBins_,
1767  lowedge,
1768  highedge);
1769 
1771  "means_dxy_eta", "#LT d_{xy} #GT vs #eta sector;#eta (sector);#LT d_{xy} #GT [#mum]", nBins_, lowedge, highedge);
1772 
1773  a_dxyEtaWidthTrend = WidthTrendsDir.make<TH1F>("widths_dxy_eta",
1774  "#sigma_{d_{xy}} vs #eta sector;#eta (sector);#sigma_{d_{xy}} [#mum]",
1775  nBins_,
1776  lowedge,
1777  highedge);
1778 
1780  "means_dz_eta", "#LT d_{z} #GT vs #eta sector;#eta (sector);#LT d_{z} #GT [#mum]", nBins_, lowedge, highedge);
1781 
1783  "widths_dz_eta", "#sigma_{d_{xy}} vs #eta sector;#eta (sector);#sigma_{d_{z}} [#mum]", nBins_, lowedge, highedge);
1784 
1786  "norm_means_dxy_phi",
1787  "#LT d_{xy}/#sigma_{d_{xy}} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{xy}/#sigma_{d_{xy}} #GT",
1788  nBins_,
1789  lowedge,
1790  highedge);
1791 
1793  "norm_widths_dxy_phi",
1794  "width(d_{xy}/#sigma_{d_{xy}}) vs #phi sector;#varphi (sector) [degrees]; width(d_{xy}/#sigma_{d_{xy}})",
1795  nBins_,
1796  lowedge,
1797  highedge);
1798 
1800  "norm_means_dz_phi",
1801  "#LT d_{z}/#sigma_{d_{z}} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{z}/#sigma_{d_{z}} #GT",
1802  nBins_,
1803  lowedge,
1804  highedge);
1805 
1807  "norm_widths_dz_phi",
1808  "width(d_{z}/#sigma_{d_{z}}) vs #phi sector;#varphi (sector) [degrees];width(d_{z}/#sigma_{d_{z}})",
1809  nBins_,
1810  lowedge,
1811  highedge);
1812 
1814  "norm_means_dxy_eta",
1815  "#LT d_{xy}/#sigma_{d_{xy}} #GT vs #eta sector;#eta (sector);#LT d_{xy}/#sigma_{d_{z}} #GT",
1816  nBins_,
1817  lowedge,
1818  highedge);
1819 
1821  "norm_widths_dxy_eta",
1822  "width(d_{xy}/#sigma_{d_{xy}}) vs #eta sector;#eta (sector);width(d_{xy}/#sigma_{d_{z}})",
1823  nBins_,
1824  lowedge,
1825  highedge);
1826 
1828  MeanTrendsDir.make<TH1F>("norm_means_dz_eta",
1829  "#LT d_{z}/#sigma_{d_{z}} #GT vs #eta sector;#eta (sector);#LT d_{z}/#sigma_{d_{z}} #GT",
1830  nBins_,
1831  lowedge,
1832  highedge);
1833 
1835  WidthTrendsDir.make<TH1F>("norm_widths_dz_eta",
1836  "width(d_{z}/#sigma_{d_{z}}) vs #eta sector;#eta (sector);width(d_{z}/#sigma_{d_{z}})",
1837  nBins_,
1838  lowedge,
1839  highedge);
1840 
1841  // means and widhts vs pT and pTCentral
1842 
1843  a_dxypTMeanTrend = MeanTrendsDir.make<TH1F>("means_dxy_pT",
1844  "#LT d_{xy} #GT vs pT;p_{T} [GeV];#LT d_{xy} #GT [#mum]",
1845  mypT_bins_.size() - 1,
1846  mypT_bins_.data());
1847 
1848  a_dxypTWidthTrend = WidthTrendsDir.make<TH1F>("widths_dxy_pT",
1849  "#sigma_{d_{xy}} vs pT;p_{T} [GeV];#sigma_{d_{xy}} [#mum]",
1850  mypT_bins_.size() - 1,
1851  mypT_bins_.data());
1852 
1854  "means_dz_pT", "#LT d_{z} #GT vs pT;p_{T} [GeV];#LT d_{z} #GT [#mum]", mypT_bins_.size() - 1, mypT_bins_.data());
1855 
1856  a_dzpTWidthTrend = WidthTrendsDir.make<TH1F>("widths_dz_pT",
1857  "#sigma_{d_{z}} vs pT;p_{T} [GeV];#sigma_{d_{z}} [#mum]",
1858  mypT_bins_.size() - 1,
1859  mypT_bins_.data());
1860 
1862  MeanTrendsDir.make<TH1F>("norm_means_dxy_pT",
1863  "#LT d_{xy}/#sigma_{d_{xy}} #GT vs pT;p_{T} [GeV];#LT d_{xy}/#sigma_{d_{xy}} #GT",
1864  mypT_bins_.size() - 1,
1865  mypT_bins_.data());
1866 
1868  WidthTrendsDir.make<TH1F>("norm_widths_dxy_pT",
1869  "width(d_{xy}/#sigma_{d_{xy}}) vs pT;p_{T} [GeV]; width(d_{xy}/#sigma_{d_{xy}})",
1870  mypT_bins_.size() - 1,
1871  mypT_bins_.data());
1872 
1873  n_dzpTMeanTrend =
1874  MeanTrendsDir.make<TH1F>("norm_means_dz_pT",
1875  "#LT d_{z}/#sigma_{d_{z}} #GT vs pT;p_{T} [GeV];#LT d_{z}/#sigma_{d_{z}} #GT",
1876  mypT_bins_.size() - 1,
1877  mypT_bins_.data());
1878 
1880  WidthTrendsDir.make<TH1F>("norm_widths_dz_pT",
1881  "width(d_{z}/#sigma_{d_{z}}) vs pT;p_{T} [GeV];width(d_{z}/#sigma_{d_{z}})",
1882  mypT_bins_.size() - 1,
1883  mypT_bins_.data());
1884 
1886  MeanTrendsDir.make<TH1F>("means_dxy_pTCentral",
1887  "#LT d_{xy} #GT vs p_{T};p_{T}(|#eta|<1.) [GeV];#LT d_{xy} #GT [#mum]",
1888  mypT_bins_.size() - 1,
1889  mypT_bins_.data());
1890 
1892  WidthTrendsDir.make<TH1F>("widths_dxy_pTCentral",
1893  "#sigma_{d_{xy}} vs p_{T};p_{T}(|#eta|<1.) [GeV];#sigma_{d_{xy}} [#mum]",
1894  mypT_bins_.size() - 1,
1895  mypT_bins_.data());
1896 
1898  MeanTrendsDir.make<TH1F>("means_dz_pTCentral",
1899  "#LT d_{z} #GT vs p_{T};p_{T}(|#eta|<1.) [GeV];#LT d_{z} #GT [#mum]",
1900  mypT_bins_.size() - 1,
1901  mypT_bins_.data());
1902 
1904  WidthTrendsDir.make<TH1F>("widths_dz_pTCentral",
1905  "#sigma_{d_{z}} vs p_{T};p_{T}(|#eta|<1.) [GeV];#sigma_{d_{z}} [#mum]",
1906  mypT_bins_.size() - 1,
1907  mypT_bins_.data());
1908 
1910  "norm_means_dxy_pTCentral",
1911  "#LT d_{xy}/#sigma_{d_{xy}} #GT vs p_{T};p_{T}(|#eta|<1.) [GeV];#LT d_{xy}/#sigma_{d_{z}} #GT",
1912  mypT_bins_.size() - 1,
1913  mypT_bins_.data());
1914 
1916  "norm_widths_dxy_pTCentral",
1917  "width(d_{xy}/#sigma_{d_{xy}}) vs p_{T};p_{T}(|#eta|<1.) [GeV];width(d_{xy}/#sigma_{d_{z}})",
1918  mypT_bins_.size() - 1,
1919  mypT_bins_.data());
1920 
1922  "norm_means_dz_pTCentral",
1923  "#LT d_{z}/#sigma_{d_{z}} #GT vs p_{T};p_{T}(|#eta|<1.) [GeV];#LT d_{z}/#sigma_{d_{z}} #GT",
1924  mypT_bins_.size() - 1,
1925  mypT_bins_.data());
1926 
1928  "norm_widths_dz_pTCentral",
1929  "width(d_{z}/#sigma_{d_{z}}) vs p_{T};p_{T}(|#eta|<1.) [GeV];width(d_{z}/#sigma_{d_{z}})",
1930  mypT_bins_.size() - 1,
1931  mypT_bins_.data());
1932 
1933  // 2D maps
1934 
1935  a_dxyMeanMap = Mean2DMapsDir.make<TH2F>("means_dxy_map",
1936  "#LT d_{xy} #GT map;#eta (sector);#varphi (sector) [degrees]",
1937  nBins_,
1938  lowedge,
1939  highedge,
1940  nBins_,
1941  lowedge,
1942  highedge);
1943 
1944  a_dzMeanMap = Mean2DMapsDir.make<TH2F>("means_dz_map",
1945  "#LT d_{z} #GT map;#eta (sector);#varphi (sector) [degrees]",
1946  nBins_,
1947  lowedge,
1948  highedge,
1949  nBins_,
1950  lowedge,
1951  highedge);
1952 
1953  n_dxyMeanMap = Mean2DMapsDir.make<TH2F>("norm_means_dxy_map",
1954  "#LT d_{xy}/#sigma_{d_{xy}} #GT map;#eta (sector);#varphi (sector) [degrees]",
1955  nBins_,
1956  lowedge,
1957  highedge,
1958  nBins_,
1959  lowedge,
1960  highedge);
1961 
1962  n_dzMeanMap = Mean2DMapsDir.make<TH2F>("norm_means_dz_map",
1963  "#LT d_{z}/#sigma_{d_{z}} #GT map;#eta (sector);#varphi (sector) [degrees]",
1964  nBins_,
1965  lowedge,
1966  highedge,
1967  nBins_,
1968  lowedge,
1969  highedge);
1970 
1971  a_dxyWidthMap = Width2DMapsDir.make<TH2F>("widths_dxy_map",
1972  "#sigma_{d_{xy}} map;#eta (sector);#varphi (sector) [degrees]",
1973  nBins_,
1974  lowedge,
1975  highedge,
1976  nBins_,
1977  lowedge,
1978  highedge);
1979 
1980  a_dzWidthMap = Width2DMapsDir.make<TH2F>("widths_dz_map",
1981  "#sigma_{d_{z}} map;#eta (sector);#varphi (sector) [degrees]",
1982  nBins_,
1983  lowedge,
1984  highedge,
1985  nBins_,
1986  lowedge,
1987  highedge);
1988 
1989  n_dxyWidthMap =
1990  Width2DMapsDir.make<TH2F>("norm_widths_dxy_map",
1991  "width(d_{xy}/#sigma_{d_{xy}}) map;#eta (sector);#varphi (sector) [degrees]",
1992  nBins_,
1993  lowedge,
1994  highedge,
1995  nBins_,
1996  lowedge,
1997  highedge);
1998 
1999  n_dzWidthMap = Width2DMapsDir.make<TH2F>("norm_widths_dz_map",
2000  "width(d_{z}/#sigma_{d_{z}}) map;#eta (sector);#varphi (sector) [degrees]",
2001  nBins_,
2002  lowedge,
2003  highedge,
2004  nBins_,
2005  lowedge,
2006  highedge);
2007 
2008  // medians and MADs
2009 
2011  MedianTrendsDir.make<TH1F>("medians_dxy_phi",
2012  "Median of d_{xy} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{xy}) [#mum]",
2013  nBins_,
2014  lowedge,
2015  highedge);
2016 
2018  "MADs_dxy_phi",
2019  "Median absolute deviation of d_{xy} vs #phi sector;#varphi (sector) [degrees];MAD(d_{xy}) [#mum]",
2020  nBins_,
2021  lowedge,
2022  highedge);
2023 
2025  MedianTrendsDir.make<TH1F>("medians_dz_phi",
2026  "Median of d_{z} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{z}) [#mum]",
2027  nBins_,
2028  lowedge,
2029  highedge);
2030 
2032  "MADs_dz_phi",
2033  "Median absolute deviation of d_{z} vs #phi sector;#varphi (sector) [degrees];MAD(d_{z}) [#mum]",
2034  nBins_,
2035  lowedge,
2036  highedge);
2037 
2039  MedianTrendsDir.make<TH1F>("medians_dxy_eta",
2040  "Median of d_{xy} vs #eta sector;#eta (sector);#mu_{1/2}(d_{xy}) [#mum]",
2041  nBins_,
2042  lowedge,
2043  highedge);
2044 
2046  MADTrendsDir.make<TH1F>("MADs_dxy_eta",
2047  "Median absolute deviation of d_{xy} vs #eta sector;#eta (sector);MAD(d_{xy}) [#mum]",
2048  nBins_,
2049  lowedge,
2050  highedge);
2051 
2053  MedianTrendsDir.make<TH1F>("medians_dz_eta",
2054  "Median of d_{z} vs #eta sector;#eta (sector);#mu_{1/2}(d_{z}) [#mum]",
2055  nBins_,
2056  lowedge,
2057  highedge);
2058 
2059  a_dzEtaMADTrend =
2060  MADTrendsDir.make<TH1F>("MADs_dz_eta",
2061  "Median absolute deviation of d_{z} vs #eta sector;#eta (sector);MAD(d_{z}) [#mum]",
2062  nBins_,
2063  lowedge,
2064  highedge);
2065 
2067  "norm_medians_dxy_phi",
2068  "Median of d_{xy}/#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{xy}/#sigma_{d_{xy}})",
2069  nBins_,
2070  lowedge,
2071  highedge);
2072 
2073  n_dxyPhiMADTrend = MADTrendsDir.make<TH1F>("norm_MADs_dxy_phi",
2074  "Median absolute deviation of d_{xy}/#sigma_{d_{xy}} vs #phi "
2075  "sector;#varphi (sector) [degrees]; MAD(d_{xy}/#sigma_{d_{xy}})",
2076  nBins_,
2077  lowedge,
2078  highedge);
2079 
2081  "norm_medians_dz_phi",
2082  "Median of d_{z}/#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{z}/#sigma_{d_{z}})",
2083  nBins_,
2084  lowedge,
2085  highedge);
2086 
2087  n_dzPhiMADTrend = MADTrendsDir.make<TH1F>("norm_MADs_dz_phi",
2088  "Median absolute deviation of d_{z}/#sigma_{d_{z}} vs #phi sector;#varphi "
2089  "(sector) [degrees];MAD(d_{z}/#sigma_{d_{z}})",
2090  nBins_,
2091  lowedge,
2092  highedge);
2093 
2095  "norm_medians_dxy_eta",
2096  "Median of d_{xy}/#sigma_{d_{xy}} vs #eta sector;#eta (sector);#mu_{1/2}(d_{xy}/#sigma_{d_{z}})",
2097  nBins_,
2098  lowedge,
2099  highedge);
2100 
2102  "norm_MADs_dxy_eta",
2103  "Median absolute deviation of d_{xy}/#sigma_{d_{xy}} vs #eta sector;#eta (sector);MAD(d_{xy}/#sigma_{d_{z}})",
2104  nBins_,
2105  lowedge,
2106  highedge);
2107 
2109  "norm_medians_dz_eta",
2110  "Median of d_{z}/#sigma_{d_{z}} vs #eta sector;#eta (sector);#mu_{1/2}(d_{z}/#sigma_{d_{z}})",
2111  nBins_,
2112  lowedge,
2113  highedge);
2114 
2116  "norm_MADs_dz_eta",
2117  "Median absolute deviation of d_{z}/#sigma_{d_{z}} vs #eta sector;#eta (sector);MAD(d_{z}/#sigma_{d_{z}})",
2118  nBins_,
2119  lowedge,
2120  highedge);
2121 
2123  //
2124  // plots of biased residuals
2125  // The vertex includes the probe track
2126  //
2128 
2129  if (useTracksFromRecoVtx_) {
2130  TFileDirectory AbsTransPhiBiasRes = fs->mkdir("Abs_Transv_Phi_BiasResiduals");
2132 
2133  TFileDirectory AbsTransEtaBiasRes = fs->mkdir("Abs_Transv_Eta_BiasResiduals");
2135 
2136  TFileDirectory AbsLongPhiBiasRes = fs->mkdir("Abs_Long_Phi_BiasResiduals");
2138 
2139  TFileDirectory AbsLongEtaBiasRes = fs->mkdir("Abs_Long_Eta_BiasResiduals");
2141 
2142  TFileDirectory NormTransPhiBiasRes = fs->mkdir("Norm_Transv_Phi_BiasResiduals");
2144 
2145  TFileDirectory NormTransEtaBiasRes = fs->mkdir("Norm_Transv_Eta_BiasResiduals");
2147 
2148  TFileDirectory NormLongPhiBiasRes = fs->mkdir("Norm_Long_Phi_BiasResiduals");
2150 
2151  TFileDirectory NormLongEtaBiasRes = fs->mkdir("Norm_Long_Eta_BiasResiduals");
2153 
2154  TFileDirectory AbsDoubleDiffBiasRes = fs->mkdir("Abs_DoubleDiffBiasResiduals");
2155  TFileDirectory NormDoubleDiffBiasRes = fs->mkdir("Norm_DoubleDiffBiasResiduals");
2156 
2157  for (int i = 0; i < nBins_; ++i) {
2158  float phiF = theDetails_.trendbins[PVValHelper::phi][i];
2159  float phiL = theDetails_.trendbins[PVValHelper::phi][i + 1];
2160 
2161  for (int j = 0; j < nBins_; ++j) {
2162  float etaF = theDetails_.trendbins[PVValHelper::eta][j];
2163  float etaL = theDetails_.trendbins[PVValHelper::eta][j + 1];
2164 
2165  a_dxyBiasResidualsMap[i][j] = AbsDoubleDiffBiasRes.make<TH1F>(
2166  Form("histo_dxy_eta_plot%i_phi_plot%i", i, j),
2167  Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{xy} [#mum];tracks", etaF, etaL, phiF, phiL),
2169  -dzmax_eta,
2170  dzmax_eta);
2171 
2172  a_dzBiasResidualsMap[i][j] = AbsDoubleDiffBiasRes.make<TH1F>(
2173  Form("histo_dxy_eta_plot%i_phi_plot%i", i, j),
2174  Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{z} [#mum];tracks", etaF, etaL, phiF, phiL),
2176  -dzmax_eta,
2177  dzmax_eta);
2178 
2179  n_dxyBiasResidualsMap[i][j] = NormDoubleDiffBiasRes.make<TH1F>(
2180  Form("histo_norm_dxy_eta_plot%i_phi_plot%i", i, j),
2181  Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{xy}/#sigma_{d_{xy}};tracks",
2182  etaF,
2183  etaL,
2184  phiF,
2185  phiL),
2187  -dzmax_eta / 100,
2188  dzmax_eta / 100);
2189 
2190  n_dzBiasResidualsMap[i][j] = NormDoubleDiffBiasRes.make<TH1F>(
2191  Form("histo_norm_dxy_eta_plot%i_phi_plot%i", i, j),
2192  Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{z}/#sigma_{d_{z}};tracks",
2193  etaF,
2194  etaL,
2195  phiF,
2196  phiL),
2198  -dzmax_eta / 100,
2199  dzmax_eta / 100);
2200  }
2201  }
2202 
2203  // declaration of the directories
2204 
2205  TFileDirectory MeanBiasTrendsDir = fs->mkdir("MeanBiasTrends");
2206  TFileDirectory WidthBiasTrendsDir = fs->mkdir("WidthBiasTrends");
2207  TFileDirectory MedianBiasTrendsDir = fs->mkdir("MedianBiasTrends");
2208  TFileDirectory MADBiasTrendsDir = fs->mkdir("MADBiasTrends");
2209 
2210  TFileDirectory Mean2DBiasMapsDir = fs->mkdir("MeanBiasMaps");
2211  TFileDirectory Width2DBiasMapsDir = fs->mkdir("WidthBiasMaps");
2212 
2213  // means and widths from the fit
2214 
2216  MeanBiasTrendsDir.make<TH1F>("means_dxy_phi",
2217  "#LT d_{xy} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{xy} #GT [#mum]",
2218  nBins_,
2219  lowedge,
2220  highedge);
2221 
2222  a_dxyPhiWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>(
2223  "widths_dxy_phi",
2224  "#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees];#sigma_{d_{xy}} [#mum]",
2225  nBins_,
2226  lowedge,
2227  highedge);
2228 
2230  MeanBiasTrendsDir.make<TH1F>("means_dz_phi",
2231  "#LT d_{z} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{z} #GT [#mum]",
2232  nBins_,
2233  lowedge,
2234  highedge);
2235 
2237  WidthBiasTrendsDir.make<TH1F>("widths_dz_phi",
2238  "#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];#sigma_{d_{z}} [#mum]",
2239  nBins_,
2240  lowedge,
2241  highedge);
2242 
2243  a_dxyEtaMeanBiasTrend = MeanBiasTrendsDir.make<TH1F>(
2244  "means_dxy_eta", "#LT d_{xy} #GT vs #eta sector;#eta (sector);#LT d_{xy} #GT [#mum]", nBins_, lowedge, highedge);
2245 
2247  WidthBiasTrendsDir.make<TH1F>("widths_dxy_eta",
2248  "#sigma_{d_{xy}} vs #eta sector;#eta (sector);#sigma_{d_{xy}} [#mum]",
2249  nBins_,
2250  lowedge,
2251  highedge);
2252 
2253  a_dzEtaMeanBiasTrend = MeanBiasTrendsDir.make<TH1F>(
2254  "means_dz_eta", "#LT d_{z} #GT vs #eta sector;#eta (sector);#LT d_{z} #GT [#mum]", nBins_, lowedge, highedge);
2255 
2257  WidthBiasTrendsDir.make<TH1F>("widths_dz_eta",
2258  "#sigma_{d_{xy}} vs #eta sector;#eta (sector);#sigma_{d_{z}} [#mum]",
2259  nBins_,
2260  lowedge,
2261  highedge);
2262 
2263  n_dxyPhiMeanBiasTrend = MeanBiasTrendsDir.make<TH1F>(
2264  "norm_means_dxy_phi",
2265  "#LT d_{xy}/#sigma_{d_{xy}} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{xy}/#sigma_{d_{xy}} #GT",
2266  nBins_,
2267  lowedge,
2268  highedge);
2269 
2270  n_dxyPhiWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>(
2271  "norm_widths_dxy_phi",
2272  "width(d_{xy}/#sigma_{d_{xy}}) vs #phi sector;#varphi (sector) [degrees]; width(d_{xy}/#sigma_{d_{xy}})",
2273  nBins_,
2274  lowedge,
2275  highedge);
2276 
2277  n_dzPhiMeanBiasTrend = MeanBiasTrendsDir.make<TH1F>(
2278  "norm_means_dz_phi",
2279  "#LT d_{z}/#sigma_{d_{z}} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{z}/#sigma_{d_{z}} #GT",
2280  nBins_,
2281  lowedge,
2282  highedge);
2283 
2284  n_dzPhiWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>(
2285  "norm_widths_dz_phi",
2286  "width(d_{z}/#sigma_{d_{z}}) vs #phi sector;#varphi (sector) [degrees];width(d_{z}/#sigma_{d_{z}})",
2287  nBins_,
2288  lowedge,
2289  highedge);
2290 
2291  n_dxyEtaMeanBiasTrend = MeanBiasTrendsDir.make<TH1F>(
2292  "norm_means_dxy_eta",
2293  "#LT d_{xy}/#sigma_{d_{xy}} #GT vs #eta sector;#eta (sector);#LT d_{xy}/#sigma_{d_{z}} #GT",
2294  nBins_,
2295  lowedge,
2296  highedge);
2297 
2298  n_dxyEtaWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>(
2299  "norm_widths_dxy_eta",
2300  "width(d_{xy}/#sigma_{d_{xy}}) vs #eta sector;#eta (sector);width(d_{xy}/#sigma_{d_{z}})",
2301  nBins_,
2302  lowedge,
2303  highedge);
2304 
2305  n_dzEtaMeanBiasTrend = MeanBiasTrendsDir.make<TH1F>(
2306  "norm_means_dz_eta",
2307  "#LT d_{z}/#sigma_{d_{z}} #GT vs #eta sector;#eta (sector);#LT d_{z}/#sigma_{d_{z}} #GT",
2308  nBins_,
2309  lowedge,
2310  highedge);
2311 
2312  n_dzEtaWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>(
2313  "norm_widths_dz_eta",
2314  "width(d_{z}/#sigma_{d_{z}}) vs #eta sector;#eta (sector);width(d_{z}/#sigma_{d_{z}})",
2315  nBins_,
2316  lowedge,
2317  highedge);
2318 
2319  // 2D maps
2320 
2321  a_dxyMeanBiasMap = Mean2DBiasMapsDir.make<TH2F>("means_dxy_map",
2322  "#LT d_{xy} #GT map;#eta (sector);#varphi (sector) [degrees]",
2323  nBins_,
2324  lowedge,
2325  highedge,
2326  nBins_,
2327  lowedge,
2328  highedge);
2329 
2330  a_dzMeanBiasMap = Mean2DBiasMapsDir.make<TH2F>("means_dz_map",
2331  "#LT d_{z} #GT map;#eta (sector);#varphi (sector) [degrees]",
2332  nBins_,
2333  lowedge,
2334  highedge,
2335  nBins_,
2336  lowedge,
2337  highedge);
2338 
2340  Mean2DBiasMapsDir.make<TH2F>("norm_means_dxy_map",
2341  "#LT d_{xy}/#sigma_{d_{xy}} #GT map;#eta (sector);#varphi (sector) [degrees]",
2342  nBins_,
2343  lowedge,
2344  highedge,
2345  nBins_,
2346  lowedge,
2347  highedge);
2348 
2349  n_dzMeanBiasMap =
2350  Mean2DBiasMapsDir.make<TH2F>("norm_means_dz_map",
2351  "#LT d_{z}/#sigma_{d_{z}} #GT map;#eta (sector);#varphi (sector) [degrees]",
2352  nBins_,
2353  lowedge,
2354  highedge,
2355  nBins_,
2356  lowedge,
2357  highedge);
2358 
2359  a_dxyWidthBiasMap = Width2DBiasMapsDir.make<TH2F>("widths_dxy_map",
2360  "#sigma_{d_{xy}} map;#eta (sector);#varphi (sector) [degrees]",
2361  nBins_,
2362  lowedge,
2363  highedge,
2364  nBins_,
2365  lowedge,
2366  highedge);
2367 
2368  a_dzWidthBiasMap = Width2DBiasMapsDir.make<TH2F>("widths_dz_map",
2369  "#sigma_{d_{z}} map;#eta (sector);#varphi (sector) [degrees]",
2370  nBins_,
2371  lowedge,
2372  highedge,
2373  nBins_,
2374  lowedge,
2375  highedge);
2376 
2378  Width2DBiasMapsDir.make<TH2F>("norm_widths_dxy_map",
2379  "width(d_{xy}/#sigma_{d_{xy}}) map;#eta (sector);#varphi (sector) [degrees]",
2380  nBins_,
2381  lowedge,
2382  highedge,
2383  nBins_,
2384  lowedge,
2385  highedge);
2386 
2388  Width2DBiasMapsDir.make<TH2F>("norm_widths_dz_map",
2389  "width(d_{z}/#sigma_{d_{z}}) map;#eta (sector);#varphi (sector) [degrees]",
2390  nBins_,
2391  lowedge,
2392  highedge,
2393  nBins_,
2394  lowedge,
2395  highedge);
2396 
2397  // medians and MADs
2398 
2399  a_dxyPhiMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>(
2400  "medians_dxy_phi",
2401  "Median of d_{xy} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{xy}) [#mum]",
2402  nBins_,
2403  lowedge,
2404  highedge);
2405 
2406  a_dxyPhiMADBiasTrend = MADBiasTrendsDir.make<TH1F>(
2407  "MADs_dxy_phi",
2408  "Median absolute deviation of d_{xy} vs #phi sector;#varphi (sector) [degrees];MAD(d_{xy}) [#mum]",
2409  nBins_,
2410  lowedge,
2411  highedge);
2412 
2413  a_dzPhiMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>(
2414  "medians_dz_phi",
2415  "Median of d_{z} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{z}) [#mum]",
2416  nBins_,
2417  lowedge,
2418  highedge);
2419 
2420  a_dzPhiMADBiasTrend = MADBiasTrendsDir.make<TH1F>(
2421  "MADs_dz_phi",
2422  "Median absolute deviation of d_{z} vs #phi sector;#varphi (sector) [degrees];MAD(d_{z}) [#mum]",
2423  nBins_,
2424  lowedge,
2425  highedge);
2426 
2428  MedianBiasTrendsDir.make<TH1F>("medians_dxy_eta",
2429  "Median of d_{xy} vs #eta sector;#eta (sector);#mu_{1/2}(d_{xy}) [#mum]",
2430  nBins_,
2431  lowedge,
2432  highedge);
2433 
2434  a_dxyEtaMADBiasTrend = MADBiasTrendsDir.make<TH1F>(
2435  "MADs_dxy_eta",
2436  "Median absolute deviation of d_{xy} vs #eta sector;#eta (sector);MAD(d_{xy}) [#mum]",
2437  nBins_,
2438  lowedge,
2439  highedge);
2440 
2442  MedianBiasTrendsDir.make<TH1F>("medians_dz_eta",
2443  "Median of d_{z} vs #eta sector;#eta (sector);#mu_{1/2}(d_{z}) [#mum]",
2444  nBins_,
2445  lowedge,
2446  highedge);
2447 
2449  MADBiasTrendsDir.make<TH1F>("MADs_dz_eta",
2450  "Median absolute deviation of d_{z} vs #eta sector;#eta (sector);MAD(d_{z}) [#mum]",
2451  nBins_,
2452  lowedge,
2453  highedge);
2454 
2455  n_dxyPhiMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>(
2456  "norm_medians_dxy_phi",
2457  "Median of d_{xy}/#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{xy}/#sigma_{d_{xy}})",
2458  nBins_,
2459  lowedge,
2460  highedge);
2461 
2462  n_dxyPhiMADBiasTrend = MADBiasTrendsDir.make<TH1F>("norm_MADs_dxy_phi",
2463  "Median absolute deviation of d_{xy}/#sigma_{d_{xy}} vs #phi "
2464  "sector;#varphi (sector) [degrees]; MAD(d_{xy}/#sigma_{d_{xy}})",
2465  nBins_,
2466  lowedge,
2467  highedge);
2468 
2469  n_dzPhiMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>(
2470  "norm_medians_dz_phi",
2471  "Median of d_{z}/#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{z}/#sigma_{d_{z}})",
2472  nBins_,
2473  lowedge,
2474  highedge);
2475 
2476  n_dzPhiMADBiasTrend = MADBiasTrendsDir.make<TH1F>("norm_MADs_dz_phi",
2477  "Median absolute deviation of d_{z}/#sigma_{d_{z}} vs #phi "
2478  "sector;#varphi (sector) [degrees];MAD(d_{z}/#sigma_{d_{z}})",
2479  nBins_,
2480  lowedge,
2481  highedge);
2482 
2483  n_dxyEtaMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>(
2484  "norm_medians_dxy_eta",
2485  "Median of d_{xy}/#sigma_{d_{xy}} vs #eta sector;#eta (sector);#mu_{1/2}(d_{xy}/#sigma_{d_{z}})",
2486  nBins_,
2487  lowedge,
2488  highedge);
2489 
2490  n_dxyEtaMADBiasTrend = MADBiasTrendsDir.make<TH1F>(
2491  "norm_MADs_dxy_eta",
2492  "Median absolute deviation of d_{xy}/#sigma_{d_{xy}} vs #eta sector;#eta (sector);MAD(d_{xy}/#sigma_{d_{z}})",
2493  nBins_,
2494  lowedge,
2495  highedge);
2496 
2497  n_dzEtaMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>(
2498  "norm_medians_dz_eta",
2499  "Median of d_{z}/#sigma_{d_{z}} vs #eta sector;#eta (sector);#mu_{1/2}(d_{z}/#sigma_{d_{z}})",
2500  nBins_,
2501  lowedge,
2502  highedge);
2503 
2504  n_dzEtaMADBiasTrend = MADBiasTrendsDir.make<TH1F>(
2505  "norm_MADs_dz_eta",
2506  "Median absolute deviation of d_{z}/#sigma_{d_{z}} vs #eta sector;#eta (sector);MAD(d_{z}/#sigma_{d_{z}})",
2507  nBins_,
2508  lowedge,
2509  highedge);
2510  }
2511 }
2512 
2513 // ------------ method called once every run before doing the event loop ----------------
2515  //=======================================================
2516  // check on the b-field consistency
2517  //=======================================================
2518 
2519  if (!isBFieldConsistentWithMode(iSetup)) {
2520  edm::LogWarning("PrimaryVertexValidation")
2521  << "*********************************************************************************\n"
2522  << "* The configuration (ptOfProbe > " << ptOfProbe_
2523  << "GeV) is not correctly set for current value of magnetic field \n"
2524  << "* Switching it to 0. !!! \n"
2525  << "*********************************************************************************" << std::endl;
2526  ptOfProbe_ = 0.;
2527  }
2528 
2529  //=======================================================
2530  // Run numbers analysis
2531  //=======================================================
2532 
2533  const auto runNumber_ = iRun.run();
2534  h_runNumber->Fill(runNumber_);
2535 
2536  if (!runNumbersTimesLog_.count(runNumber_)) {
2537  auto times = getRunTime(iSetup);
2538 
2539  if (debug_) {
2540  const time_t start_time = times.first / 1000000;
2541  edm::LogInfo("PrimaryVertexValidation")
2542  << runNumber_ << " has start time: " << times.first << " - " << times.second << std::endl;
2543  edm::LogInfo("PrimaryVertexValidation")
2544  << "human readable time: " << std::asctime(std::gmtime(&start_time)) << std::endl;
2545  }
2546 
2548  }
2549 
2550  if (h_runFromEvent->GetEntries() == 0) {
2551  h_runFromEvent->SetBinContent(1, RunNumber_);
2552  }
2553 
2554  //=======================================================
2555  // Retrieve tracker topology from geometry
2556  //=======================================================
2557 
2558  const TrackerTopology* const tTopo = &iSetup.getData(topoTokenBR_);
2559 
2560  //=======================================================
2561  // Retrieve geometry information
2562  //=======================================================
2563 
2564  edm::LogInfo("read tracker geometry...");
2566  edm::LogInfo("tracker geometry read") << "There are: " << pDD->dets().size() << " detectors";
2567 
2568  PixelTopologyMap PTMap = PixelTopologyMap(pDD.product(), tTopo);
2569  nLadders_ = PTMap.getPXBLadders(1);
2570  nModZ_ = PTMap.getPXBModules(1);
2571 
2572  //=======================================================
2573  // shrink to fit
2574  //=======================================================
2575 
2576  if (h_dxy_ladderOverlap_.size() != nLadders_) {
2583 
2584  if (debug_) {
2585  edm::LogInfo("PrimaryVertexValidation") << "checking size:" << h_dxy_ladder_.size() << std::endl;
2586  }
2587  }
2588 
2589  if (h_dxy_modZ_.size() != nModZ_) {
2594 
2595  if (debug_) {
2596  edm::LogInfo("PrimaryVertexValidation") << "checking size:" << h_dxy_modZ_.size() << std::endl;
2597  }
2598  }
2599 
2601  // switch on the phase-2
2603  if (debug_) {
2604  edm::LogInfo("PrimaryVertexValidation")
2605  << " pixel phase2 setup, nLadders: " << nLadders_ << " nModules:" << nModZ_;
2606  }
2607 
2608  } else if ((pDD->isThere(GeomDetEnumerators::P1PXB)) || (pDD->isThere(GeomDetEnumerators::P1PXEC))) {
2609  // switch on the phase-1
2611  if (debug_) {
2612  edm::LogInfo("PrimaryVertexValidation")
2613  << " pixel phase1 setup, nLadders: " << nLadders_ << " nModules:" << nModZ_;
2614  }
2615  } else {
2616  // switch on the phase-0
2618  if (debug_) {
2619  edm::LogInfo("PrimaryVertexValidation")
2620  << " pixel phase0 setup, nLadders: " << nLadders_ << " nModules:" << nModZ_;
2621  }
2622  }
2623 
2624  switch (phase_) {
2625  case PVValHelper::phase0:
2627  break;
2628  case PVValHelper::phase1:
2630  break;
2631  case PVValHelper::phase2:
2633  break;
2634  default:
2635  throw cms::Exception("LogicError") << "Unknown detector phase: " << phase_;
2636  }
2637 
2638  if (h_etaMax->GetEntries() == 0.) {
2639  h_etaMax->SetBinContent(1., etaOfProbe_);
2640  h_nbins->SetBinContent(1., nBins_);
2641  h_nLadders->SetBinContent(1., nLadders_);
2642  h_nModZ->SetBinContent(1., nModZ_);
2643  h_pTinfo->SetBinContent(1., mypT_bins_.size());
2644  h_pTinfo->SetBinContent(2., minPt_);
2645  h_pTinfo->SetBinContent(3., maxPt_);
2646  }
2647 }
2648 
2649 // ------------ method called once each job just after ending the event loop ------------
2651  // shring the histograms to fit
2652  h_probeL1Ladder_->GetXaxis()->SetRangeUser(-1.5, nLadders_ + 0.5);
2653  h_probeL1Module_->GetXaxis()->SetRangeUser(-1.5, nModZ_ + 0.5);
2654  h2_probeLayer1Map_->GetXaxis()->SetRangeUser(0.5, nModZ_ + 0.5);
2655  h2_probeLayer1Map_->GetYaxis()->SetRangeUser(0.5, nLadders_ + 0.5);
2656  h2_probePassingLayer1Map_->GetXaxis()->SetRangeUser(0.5, nModZ_ + 0.5);
2657  h2_probePassingLayer1Map_->GetYaxis()->SetRangeUser(0.5, nLadders_ + 0.5);
2658 
2659  TFileDirectory RunFeatures = fs->mkdir("RunFeatures");
2660  h_runStartTimes = RunFeatures.make<TH1I>(
2661  "runStartTimes", "run start times", runNumbersTimesLog_.size(), 0, runNumbersTimesLog_.size());
2662  h_runEndTimes =
2663  RunFeatures.make<TH1I>("runEndTimes", "run end times", runNumbersTimesLog_.size(), 0, runNumbersTimesLog_.size());
2664 
2665  unsigned int count = 1;
2666  for (const auto& run : runNumbersTimesLog_) {
2667  // strip down the microseconds
2668  h_runStartTimes->SetBinContent(count, run.second.first / 10e6);
2669  h_runStartTimes->GetXaxis()->SetBinLabel(count, (std::to_string(run.first)).c_str());
2670 
2671  h_runEndTimes->SetBinContent(count, run.second.second / 10e6);
2672  h_runEndTimes->GetXaxis()->SetBinLabel(count, (std::to_string(run.first)).c_str());
2673 
2674  count++;
2675  }
2676 
2677  edm::LogInfo("PrimaryVertexValidation") << "######################################\n"
2678  << "# PrimaryVertexValidation::endJob()\n"
2679  << "# Number of analyzed events: " << Nevt_ << "\n"
2680  << "######################################";
2681 
2682  // means and widhts vs ladder and module number
2683 
2685  "means_dxy_modZ", "#LT d_{xy} #GT vs modZ;module number (Z);#LT d_{xy} #GT [#mum]", nModZ_, 0., nModZ_);
2686 
2688  "widths_dxy_modZ", "#sigma_{d_{xy}} vs modZ;module number (Z);#sigma_{d_{xy}} [#mum]", nModZ_, 0., nModZ_);
2689 
2691  "means_dz_modZ", "#LT d_{z} #GT vs modZ;module number (Z);#LT d_{z} #GT [#mum]", nModZ_, 0., nModZ_);
2692 
2694  "widths_dz_modZ", "#sigma_{d_{z}} vs modZ;module number (Z);#sigma_{d_{z}} [#mum]", nModZ_, 0., nModZ_);
2695 
2696  a_dxyladderMeanTrend = MeanTrendsDir.make<TH1F>("means_dxy_ladder",
2697  "#LT d_{xy} #GT vs ladder;ladder number (#phi);#LT d_{xy} #GT [#mum]",
2698  nLadders_,
2699  0.,
2700  nLadders_);
2701 
2703  WidthTrendsDir.make<TH1F>("widths_dxy_ladder",
2704  "#sigma_{d_{xy}} vs ladder;ladder number (#phi);#sigma_{d_{xy}} [#mum]",
2705  nLadders_,
2706  0.,
2707  nLadders_);
2708 
2710  "means_dz_ladder", "#LT d_{z} #GT vs ladder;ladder number (#phi);#LT d_{z} #GT [#mum]", nLadders_, 0., nLadders_);
2711 
2713  WidthTrendsDir.make<TH1F>("widths_dz_ladder",
2714  "#sigma_{d_{z}} vs ladder;ladder number (#phi);#sigma_{d_{z}} [#mum]",
2715  nLadders_,
2716  0.,
2717  nLadders_);
2718 
2720  "norm_means_dxy_modZ",
2721  "#LT d_{xy}/#sigma_{d_{xy}} #GT vs modZ;module number (Z);#LT d_{xy}/#sigma_{d_{xy}} #GT",
2722  nModZ_,
2723  0.,
2724  nModZ_);
2725 
2727  "norm_widths_dxy_modZ",
2728  "width(d_{xy}/#sigma_{d_{xy}}) vs modZ;module number (Z); width(d_{xy}/#sigma_{d_{xy}})",
2729  nModZ_,
2730  0.,
2731  nModZ_);
2732 
2734  MeanTrendsDir.make<TH1F>("norm_means_dz_modZ",
2735  "#LT d_{z}/#sigma_{d_{z}} #GT vs modZ;module number (Z);#LT d_{z}/#sigma_{d_{z}} #GT",
2736  nModZ_,
2737  0.,
2738  nModZ_);
2739 
2741  WidthTrendsDir.make<TH1F>("norm_widths_dz_modZ",
2742  "width(d_{z}/#sigma_{d_{z}}) vs pT;module number (Z);width(d_{z}/#sigma_{d_{z}})",
2743  nModZ_,
2744  0.,
2745  nModZ_);
2746 
2748  "norm_means_dxy_ladder",
2749  "#LT d_{xy}/#sigma_{d_{xy}} #GT vs ladder;ladder number (#phi);#LT d_{xy}/#sigma_{d_{z}} #GT",
2750  nLadders_,
2751  0.,
2752  nLadders_);
2753 
2755  "norm_widths_dxy_ladder",
2756  "width(d_{xy}/#sigma_{d_{xy}}) vs ladder;ladder number (#phi);width(d_{xy}/#sigma_{d_{z}})",
2757  nLadders_,
2758  0.,
2759  nLadders_);
2760 
2762  "norm_means_dz_ladder",
2763  "#LT d_{z}/#sigma_{d_{z}} #GT vs ladder;ladder number (#phi);#LT d_{z}/#sigma_{d_{z}} #GT",
2764  nLadders_,
2765  0.,
2766  nLadders_);
2767 
2769  "norm_widths_dz_ladder",
2770  "width(d_{z}/#sigma_{d_{z}}) vs ladder;ladder number (#phi);width(d_{z}/#sigma_{d_{z}})",
2771  nLadders_,
2772  0.,
2773  nLadders_);
2774 
2775  // 2D maps of residuals in bins of L1 modules
2776 
2777  a_dxyL1MeanMap = Mean2DMapsDir.make<TH2F>("means_dxy_l1map",
2778  "#LT d_{xy} #GT map;module number [z];ladder number [#varphi]",
2779  nModZ_,
2780  0.,
2781  nModZ_,
2782  nLadders_,
2783  0.,
2784  nLadders_);
2785 
2786  a_dzL1MeanMap = Mean2DMapsDir.make<TH2F>("means_dz_l1map",
2787  "#LT d_{z} #GT map;module number [z];ladder number [#varphi]",
2788  nModZ_,
2789  0.,
2790  nModZ_,
2791  nLadders_,
2792  0.,
2793  nLadders_);
2794 
2795  n_dxyL1MeanMap =
2796  Mean2DMapsDir.make<TH2F>("norm_means_dxy_l1map",
2797  "#LT d_{xy}/#sigma_{d_{xy}} #GT map;module number [z];ladder number [#varphi]",
2798  nModZ_,
2799  0.,
2800  nModZ_,
2801  nLadders_,
2802  0.,
2803  nLadders_);
2804 
2805  n_dzL1MeanMap = Mean2DMapsDir.make<TH2F>("norm_means_dz_l1map",
2806  "#LT d_{z}/#sigma_{d_{z}} #GT map;module number [z];ladder number [#varphi]",
2807  nModZ_,
2808  0.,
2809  nModZ_,
2810  nLadders_,
2811  0.,
2812  nLadders_);
2813 
2814  a_dxyL1WidthMap = Width2DMapsDir.make<TH2F>("widths_dxy_l1map",
2815  "#sigma_{d_{xy}} map;module number [z];ladder number [#varphi]",
2816  nModZ_,
2817  0.,
2818  nModZ_,
2819  nLadders_,
2820  0.,
2821  nLadders_);
2822 
2823  a_dzL1WidthMap = Width2DMapsDir.make<TH2F>("widths_dz_l1map",
2824  "#sigma_{d_{z}} map;module number [z];ladder number [#varphi]",
2825  nModZ_,
2826  0.,
2827  nModZ_,
2828  nLadders_,
2829  0.,
2830  nLadders_);
2831 
2832  n_dxyL1WidthMap =
2833  Width2DMapsDir.make<TH2F>("norm_widths_dxy_l1map",
2834  "width(d_{xy}/#sigma_{d_{xy}}) map;module number [z];ladder number [#varphi]",
2835  nModZ_,
2836  0.,
2837  nModZ_,
2838  nLadders_,
2839  0.,
2840  nLadders_);
2841 
2842  n_dzL1WidthMap =
2843  Width2DMapsDir.make<TH2F>("norm_widths_dz_l1map",
2844  "width(d_{z}/#sigma_{d_{z}}) map;module number [z];ladder number [#varphi]",
2845  nModZ_,
2846  0.,
2847  nModZ_,
2848  nLadders_,
2849  0.,
2850  nLadders_);
2851 
2852  if (useTracksFromRecoVtx_) {
2857 
2862 
2867 
2872 
2873  // medians and MADs
2874 
2879 
2884 
2889 
2894 
2895  // 2d Maps
2896 
2901 
2906  }
2907 
2908  // do profiles
2909 
2914 
2919 
2924 
2929 
2930  // vs transverse momentum
2931 
2936 
2941 
2946 
2951 
2952  // vs ladder and module number
2953 
2958 
2963 
2968 
2973 
2974  // medians and MADs
2975 
2978 
2981 
2986 
2991 
2996 
2997  // 2D Maps
2998 
3003 
3008 
3009  // 2D Maps of residuals in bins of L1 modules
3010 
3015 
3020 }
3021 
3022 //*************************************************************
3023 std::pair<long long, long long> PrimaryVertexValidation::getRunTime(const edm::EventSetup& iSetup) const
3024 //*************************************************************
3025 {
3026  edm::ESHandle<RunInfo> runInfo = iSetup.getHandle(runInfoTokenBR_);
3027  if (debug_) {
3028  edm::LogInfo("PrimaryVertexValidation")
3029  << runInfo.product()->m_start_time_str << " " << runInfo.product()->m_stop_time_str << std::endl;
3030  }
3031  return std::make_pair(runInfo.product()->m_start_time_ll, runInfo.product()->m_stop_time_ll);
3032 }
3033 
3034 //*************************************************************
3036 //*************************************************************
3037 {
3038  edm::ESHandle<RunInfo> runInfo = iSetup.getHandle(runInfoTokenBR_);
3039  double average_current = runInfo.product()->m_avg_current;
3040  bool isOn = (average_current > 2000.);
3041  bool is0T = (ptOfProbe_ == 0.);
3042 
3043  return ((isOn && !is0T) || (!isOn && is0T));
3044 }
3045 
3046 //*************************************************************
3048 //*************************************************************
3049 {
3050  nTracks_ = 0;
3051  nClus_ = 0;
3052  nOfflineVertices_ = 0;
3053  RunNumber_ = 0;
3055  xOfflineVertex_ = -999.;
3056  yOfflineVertex_ = -999.;
3057  zOfflineVertex_ = -999.;
3058  xErrOfflineVertex_ = 0.;
3059  yErrOfflineVertex_ = 0.;
3060  zErrOfflineVertex_ = 0.;
3061  BSx0_ = -999.;
3062  BSy0_ = -999.;
3063  BSz0_ = -999.;
3064  Beamsigmaz_ = -999.;
3065  Beamdxdz_ = -999.;
3066  BeamWidthX_ = -999.;
3067  BeamWidthY_ = -999.;
3068  wxy2_ = -999.;
3069 
3070  for (int i = 0; i < nMaxtracks_; ++i) {
3071  pt_[i] = 0;
3072  p_[i] = 0;
3073  nhits_[i] = 0;
3074  nhits1D_[i] = 0;
3075  nhits2D_[i] = 0;
3076  nhitsBPIX_[i] = 0;
3077  nhitsFPIX_[i] = 0;
3078  nhitsTIB_[i] = 0;
3079  nhitsTID_[i] = 0;
3080  nhitsTOB_[i] = 0;
3081  nhitsTEC_[i] = 0;
3082  isHighPurity_[i] = 0;
3083  eta_[i] = 0;
3084  theta_[i] = 0;
3085  phi_[i] = 0;
3086  chi2_[i] = 0;
3087  chi2ndof_[i] = 0;
3088  charge_[i] = 0;
3089  qoverp_[i] = 0;
3090  dz_[i] = 0;
3091  dxy_[i] = 0;
3092  dzBs_[i] = 0;
3093  dxyBs_[i] = 0;
3094  xPCA_[i] = 0;
3095  yPCA_[i] = 0;
3096  zPCA_[i] = 0;
3097  xUnbiasedVertex_[i] = 0;
3098  yUnbiasedVertex_[i] = 0;
3099  zUnbiasedVertex_[i] = 0;
3101  chi2UnbiasedVertex_[i] = 0;
3103  DOFUnbiasedVertex_[i] = 0;
3106  dxyFromMyVertex_[i] = 0;
3107  dzFromMyVertex_[i] = 0;
3108  d3DFromMyVertex_[i] = 0;
3109  dxyErrorFromMyVertex_[i] = 0;
3110  dzErrorFromMyVertex_[i] = 0;
3111  d3DErrorFromMyVertex_[i] = 0;
3112  IPTsigFromMyVertex_[i] = 0;
3113  IPLsigFromMyVertex_[i] = 0;
3114  IP3DsigFromMyVertex_[i] = 0;
3115  hasRecVertex_[i] = 0;
3116  isGoodTrack_[i] = 0;
3117  }
3118 }
3119 
3120 //*************************************************************
3122  TH1F* residualsPlot[100],
3123  PVValHelper::estimator fitPar_,
3124  const std::string& var_)
3125 //*************************************************************
3126 {
3127  for (int i = 0; i < nBins_; ++i) {
3128  char phibincenter[129];
3130  sprintf(phibincenter, "%.f", (phiBins[i] + phiBins[i + 1]) / 2.);
3131 
3132  char etabincenter[129];
3134  sprintf(etabincenter, "%.1f", (etaBins[i] + etaBins[i + 1]) / 2.);
3135 
3136  switch (fitPar_) {
3137  case PVValHelper::MEAN: {
3138  float mean_ = PVValHelper::fitResiduals(residualsPlot[i]).first.value();
3139  float meanErr_ = PVValHelper::fitResiduals(residualsPlot[i]).first.error();
3140  trendPlot->SetBinContent(i + 1, mean_);
3141  trendPlot->SetBinError(i + 1, meanErr_);
3142  break;
3143  }
3144  case PVValHelper::WIDTH: {
3145  float width_ = PVValHelper::fitResiduals(residualsPlot[i]).second.value();
3146  float widthErr_ = PVValHelper::fitResiduals(residualsPlot[i]).second.error();
3147  trendPlot->SetBinContent(i + 1, width_);
3148  trendPlot->SetBinError(i + 1, widthErr_);
3149  break;
3150  }
3151  case PVValHelper::MEDIAN: {
3152  float median_ = PVValHelper::getMedian(residualsPlot[i]).value();
3153  float medianErr_ = PVValHelper::getMedian(residualsPlot[i]).error();
3154  trendPlot->SetBinContent(i + 1, median_);
3155  trendPlot->SetBinError(i + 1, medianErr_);
3156  break;
3157  }
3158  case PVValHelper::MAD: {
3159  float mad_ = PVValHelper::getMAD(residualsPlot[i]).value();
3160  float madErr_ = PVValHelper::getMAD(residualsPlot[i]).error();
3161  trendPlot->SetBinContent(i + 1, mad_);
3162  trendPlot->SetBinError(i + 1, madErr_);
3163  break;
3164  }
3165  default:
3166  edm::LogWarning("PrimaryVertexValidation")
3167  << "fillTrendPlot() " << fitPar_ << " unknown estimator!" << std::endl;
3168  break;
3169  }
3170 
3171  if (var_.find("eta") != std::string::npos) {
3172  trendPlot->GetXaxis()->SetBinLabel(i + 1, etabincenter);
3173  } else if (var_.find("phi") != std::string::npos) {
3174  trendPlot->GetXaxis()->SetBinLabel(i + 1, phibincenter);
3175  } else {
3176  edm::LogWarning("PrimaryVertexValidation")
3177  << "fillTrendPlot() " << var_ << " unknown track parameter!" << std::endl;
3178  }
3179  }
3180 }
3181 
3182 //*************************************************************
3184  std::vector<TH1F*>& h,
3185  PVValHelper::estimator fitPar_,
3186  PVValHelper::plotVariable plotVar)
3187 //*************************************************************
3188 {
3189  for (auto iterator = h.begin(); iterator != h.end(); iterator++) {
3190  unsigned int bin = std::distance(h.begin(), iterator) + 1;
3191  std::pair<Measurement1D, Measurement1D> myFit = PVValHelper::fitResiduals((*iterator));
3192 
3193  switch (fitPar_) {
3194  case PVValHelper::MEAN: {
3195  float mean_ = myFit.first.value();
3196  float meanErr_ = myFit.first.error();
3197  trendPlot->SetBinContent(bin, mean_);
3198  trendPlot->SetBinError(bin, meanErr_);
3199  break;
3200  }
3201  case PVValHelper::WIDTH: {
3202  float width_ = myFit.second.value();
3203  float widthErr_ = myFit.second.error();
3204  trendPlot->SetBinContent(bin, width_);
3205  trendPlot->SetBinError(bin, widthErr_);
3206  break;
3207  }
3208  case PVValHelper::MEDIAN: {
3209  float median_ = PVValHelper::getMedian(*iterator).value();
3210  float medianErr_ = PVValHelper::getMedian(*iterator).error();
3211  trendPlot->SetBinContent(bin, median_);
3212  trendPlot->SetBinError(bin, medianErr_);
3213  break;
3214  }
3215  case PVValHelper::MAD: {
3216  float mad_ = PVValHelper::getMAD(*iterator).value();
3217  float madErr_ = PVValHelper::getMAD(*iterator).error();
3218  trendPlot->SetBinContent(bin, mad_);
3219  trendPlot->SetBinError(bin, madErr_);
3220  break;
3221  }
3222  default:
3223  edm::LogWarning("PrimaryVertexValidation")
3224  << "fillTrendPlotByIndex() " << fitPar_ << " unknown estimator!" << std::endl;
3225  break;
3226  }
3227 
3228  char bincenter[129];
3229  if (plotVar == PVValHelper::eta) {
3231  sprintf(bincenter, "%.1f", (etaBins[bin - 1] + etaBins[bin]) / 2.);
3232  trendPlot->GetXaxis()->SetBinLabel(bin, bincenter);
3233  } else if (plotVar == PVValHelper::phi) {
3235  sprintf(bincenter, "%.f", (phiBins[bin - 1] + phiBins[bin]) / 2.);
3236  trendPlot->GetXaxis()->SetBinLabel(bin, bincenter);
3237  } else {
3239  //edm::LogWarning("PrimaryVertexValidation")<<"fillTrendPlotByIndex() "<< plotVar <<" unknown track parameter!"<<std::endl;
3240  }
3241  }
3242 }
3243 
3244 //*************************************************************
3246  TH1F* residualsMapPlot[100][100],
3247  PVValHelper::estimator fitPar_,
3248  const int nXBins_,
3249  const int nYBins_)
3250 //*************************************************************
3251 {
3252  for (int i = 0; i < nYBins_; ++i) {
3253  char phibincenter[129];
3255  sprintf(phibincenter, "%.f", (phiBins[i] + phiBins[i + 1]) / 2.);
3256 
3257  if (nXBins_ == nYBins_) {
3258  trendMap->GetYaxis()->SetBinLabel(i + 1, phibincenter);
3259  }
3260 
3261  for (int j = 0; j < nXBins_; ++j) {
3262  char etabincenter[129];
3264  sprintf(etabincenter, "%.1f", (etaBins[j] + etaBins[j + 1]) / 2.);
3265 
3266  if (i == 0) {
3267  if (nXBins_ == nYBins_) {
3268  trendMap->GetXaxis()->SetBinLabel(j + 1, etabincenter);
3269  }
3270  }
3271 
3272  switch (fitPar_) {
3273  case PVValHelper::MEAN: {
3274  float mean_ = PVValHelper::fitResiduals(residualsMapPlot[i][j]).first.value();
3275  float meanErr_ = PVValHelper::fitResiduals(residualsMapPlot[i][j]).first.error();
3276  trendMap->SetBinContent(j + 1, i + 1, mean_);
3277  trendMap->SetBinError(j + 1, i + 1, meanErr_);
3278  break;
3279  }
3280  case PVValHelper::WIDTH: {
3281  float width_ = PVValHelper::fitResiduals(residualsMapPlot[i][j]).second.value();
3282  float widthErr_ = PVValHelper::fitResiduals(residualsMapPlot[i][j]).second.error();
3283  trendMap->SetBinContent(j + 1, i + 1, width_);
3284  trendMap->SetBinError(j + 1, i + 1, widthErr_);
3285  break;
3286  }
3287  case PVValHelper::MEDIAN: {
3288  float median_ = PVValHelper::getMedian(residualsMapPlot[i][j]).value();
3289  float medianErr_ = PVValHelper::getMedian(residualsMapPlot[i][j]).error();
3290  trendMap->SetBinContent(j + 1, i + 1, median_);
3291  trendMap->SetBinError(j + 1, i + 1, medianErr_);
3292  break;
3293  }
3294  case PVValHelper::MAD: {
3295  float mad_ = PVValHelper::getMAD(residualsMapPlot[i][j]).value();
3296  float madErr_ = PVValHelper::getMAD(residualsMapPlot[i][j]).error();
3297  trendMap->SetBinContent(j + 1, i + 1, mad_);
3298  trendMap->SetBinError(j + 1, i + 1, madErr_);
3299  break;
3300  }
3301  default:
3302  edm::LogWarning("PrimaryVertexValidation:") << " fillMap() " << fitPar_ << " unknown estimator!" << std::endl;
3303  }
3304  } // closes loop on eta bins
3305  } // cloeses loop on phi bins
3306 }
3307 
3308 //*************************************************************
3310 //*************************************************************
3311 {
3312  if (a.tracksSize() != b.tracksSize())
3313  return a.tracksSize() > b.tracksSize() ? true : false;
3314  else
3315  return a.chi2() < b.chi2() ? true : false;
3316 }
3317 
3318 //*************************************************************
3320  const reco::Vertex& vertex,
3321  const std::string& qualityString_,
3322  double dxyErrMax_,
3323  double dzErrMax_,
3324  double ptErrMax_)
3325 //*************************************************************
3326 {
3327  math::XYZPoint vtxPoint(0.0, 0.0, 0.0);
3328  double vzErr = 0.0, vxErr = 0.0, vyErr = 0.0;
3329  vtxPoint = vertex.position();
3330  vzErr = vertex.zError();
3331  vxErr = vertex.xError();
3332  vyErr = vertex.yError();
3333 
3334  double dxy = 0.0, dz = 0.0, dxysigma = 0.0, dzsigma = 0.0;
3335  dxy = track.dxy(vtxPoint);
3336  dz = track.dz(vtxPoint);
3337  dxysigma = sqrt(track.d0Error() * track.d0Error() + vxErr * vyErr);
3338  dzsigma = sqrt(track.dzError() * track.dzError() + vzErr * vzErr);
3339 
3340  if (track.quality(reco::TrackBase::qualityByName(qualityString_)) != 1)
3341  return false;
3342  if (std::abs(dxy / dxysigma) > dxyErrMax_)
3343  return false;
3344  if (std::abs(dz / dzsigma) > dzErrMax_)
3345  return false;
3346  if (track.ptError() / track.pt() > ptErrMax_)
3347  return false;
3348 
3349  return true;
3350 }
3351 
3352 //*************************************************************
3354 //*************************************************************
3355 {
3356  TH1F::SetDefaultSumw2(kTRUE);
3357 
3358  std::map<std::string, TH1*> h;
3359 
3360  // histograms of track quality (Data and MC)
3361  std::string types[] = {"all", "sel"};
3362  for (const auto& type : types) {
3363  h["pseudorapidity_" + type] =
3364  dir.make<TH1F>(("rapidity_" + type).c_str(), "track pseudorapidity; track #eta; tracks", 100, -3., 3.);
3365  h["z0_" + type] = dir.make<TH1F>(("z0_" + type).c_str(), "track z_{0};track z_{0} (cm);tracks", 80, -40., 40.);
3366  h["phi_" + type] = dir.make<TH1F>(("phi_" + type).c_str(), "track #phi; track #phi;tracks", 80, -M_PI, M_PI);
3367  h["eta_" + type] = dir.make<TH1F>(("eta_" + type).c_str(), "track #eta; track #eta;tracks", 80, -4., 4.);
3368  h["pt_" + type] = dir.make<TH1F>(("pt_" + type).c_str(), "track p_{T}; track p_{T} [GeV];tracks", 100, 0., 20.);
3369  h["p_" + type] = dir.make<TH1F>(("p_" + type).c_str(), "track p; track p [GeV];tracks", 100, 0., 20.);
3370  h["found_" + type] =
3371  dir.make<TH1F>(("found_" + type).c_str(), "n. found hits;n^{found}_{hits};tracks", 30, 0., 30.);
3372  h["lost_" + type] = dir.make<TH1F>(("lost_" + type).c_str(), "n. lost hits;n^{lost}_{hits};tracks", 20, 0., 20.);
3373  h["nchi2_" + type] =
3374  dir.make<TH1F>(("nchi2_" + type).c_str(), "normalized track #chi^{2};track #chi^{2}/ndf;tracks", 100, 0., 20.);
3375  h["rstart_" + type] = dir.make<TH1F>(
3376  ("rstart_" + type).c_str(), "track start radius; track innermost radius r (cm);tracks", 100, 0., 20.);
3377  h["expectedInner_" + type] = dir.make<TH1F>(
3378  ("expectedInner_" + type).c_str(), "n. expected inner hits;n^{expected}_{inner};tracks", 10, 0., 10.);
3379  h["expectedOuter_" + type] = dir.make<TH1F>(
3380  ("expectedOuter_" + type).c_str(), "n. expected outer hits;n^{expected}_{outer};tracks ", 10, 0., 10.);
3381  h["logtresxy_" + type] =
3382  dir.make<TH1F>(("logtresxy_" + type).c_str(),
3383  "log10(track r-#phi resolution/#mum);log10(track r-#phi resolution/#mum);tracks",
3384  100,
3385  0.,
3386  5.);
3387  h["logtresz_" + type] = dir.make<TH1F>(("logtresz_" + type).c_str(),
3388  "log10(track z resolution/#mum);log10(track z resolution/#mum);tracks",
3389  100,
3390  0.,
3391  5.);
3392  h["tpullxy_" + type] =
3393  dir.make<TH1F>(("tpullxy_" + type).c_str(), "track r-#phi pull;pull_{r-#phi};tracks", 100, -10., 10.);
3394  h["tpullz_" + type] =
3395  dir.make<TH1F>(("tpullz_" + type).c_str(), "track r-z pull;pull_{r-z};tracks", 100, -50., 50.);
3396  h["tlogDCAxy_" + type] = dir.make<TH1F>(
3397  ("tlogDCAxy_" + type).c_str(), "track log_{10}(DCA_{r-#phi});track log_{10}(DCA_{r-#phi});tracks", 200, -5., 3.);
3398  h["tlogDCAz_" + type] = dir.make<TH1F>(
3399  ("tlogDCAz_" + type).c_str(), "track log_{10}(DCA_{r-z});track log_{10}(DCA_{r-z});tracks", 200, -5., 5.);
3400  h["lvseta_" + type] = dir.make<TH2F>(
3401  ("lvseta_" + type).c_str(), "cluster length vs #eta;track #eta;cluster length", 60, -3., 3., 20, 0., 20);
3402  h["lvstanlambda_" + type] = dir.make<TH2F>(("lvstanlambda_" + type).c_str(),
3403  "cluster length vs tan #lambda; tan#lambda;cluster length",
3404  60,
3405  -6.,
3406  6.,
3407  20,
3408  0.,
3409  20);
3410  h["restrkz_" + type] =
3411  dir.make<TH1F>(("restrkz_" + type).c_str(), "z-residuals (track vs vertex);res_{z} (cm);tracks", 200, -5., 5.);
3412  h["restrkzvsphi_" + type] = dir.make<TH2F>(("restrkzvsphi_" + type).c_str(),
3413  "z-residuals (track - vertex) vs track #phi;track #phi;res_{z} (cm)",
3414  12,
3415  -M_PI,
3416  M_PI,
3417  100,
3418  -0.5,
3419  0.5);
3420  h["restrkzvseta_" + type] = dir.make<TH2F>(("restrkzvseta_" + type).c_str(),
3421  "z-residuals (track - vertex) vs track #eta;track #eta;res_{z} (cm)",
3422  12,
3423  -3.,
3424  3.,
3425  200,
3426  -0.5,
3427  0.5);
3428  h["pulltrkzvsphi_" + type] =
3429  dir.make<TH2F>(("pulltrkzvsphi_" + type).c_str(),
3430  "normalized z-residuals (track - vertex) vs track #phi;track #phi;res_{z}/#sigma_{res_{z}}",
3431  12,
3432  -M_PI,
3433  M_PI,
3434  100,
3435  -5.,
3436  5.);
3437  h["pulltrkzvseta_" + type] =
3438  dir.make<TH2F>(("pulltrkzvseta_" + type).c_str(),
3439  "normalized z-residuals (track - vertex) vs track #eta;track #eta;res_{z}/#sigma_{res_{z}}",
3440  12,
3441  -3.,
3442  3.,
3443  100,
3444  -5.,
3445  5.);
3446  h["pulltrkz_" + type] = dir.make<TH1F>(("pulltrkz_" + type).c_str(),
3447  "normalized z-residuals (track vs vertex);res_{z}/#sigma_{res_{z}};tracks",
3448  100,
3449  -5.,
3450  5.);
3451  h["sigmatrkz0_" + type] = dir.make<TH1F>(
3452  ("sigmatrkz0_" + type).c_str(), "z-resolution (excluding beam);#sigma^{trk}_{z_{0}} (cm);tracks", 100, 0., 5.);
3453  h["sigmatrkz_" + type] = dir.make<TH1F>(
3454  ("sigmatrkz_" + type).c_str(), "z-resolution (including beam);#sigma^{trk}_{z} (cm);tracks", 100, 0., 5.);
3455  h["nbarrelhits_" + type] = dir.make<TH1F>(
3456  ("nbarrelhits_" + type).c_str(), "number of pixel barrel hits;n. hits Barrel Pixel;tracks", 10, 0., 10.);
3457  h["nbarrelLayers_" + type] = dir.make<TH1F>(
3458  ("nbarrelLayers_" + type).c_str(), "number of pixel barrel layers;n. layers Barrel Pixel;tracks", 10, 0., 10.);
3459  h["nPxLayers_" + type] = dir.make<TH1F>(
3460  ("nPxLayers_" + type).c_str(), "number of pixel layers (barrel+endcap);n. Pixel layers;tracks", 10, 0., 10.);
3461  h["nSiLayers_" + type] =
3462  dir.make<TH1F>(("nSiLayers_" + type).c_str(), "number of Tracker layers;n. Tracker layers;tracks", 20, 0., 20.);
3463  h["trackAlgo_" + type] =
3464  dir.make<TH1F>(("trackAlgo_" + type).c_str(), "track algorithm;track algo;tracks", 30, 0., 30.);
3465  h["trackQuality_" + type] =
3466  dir.make<TH1F>(("trackQuality_" + type).c_str(), "track quality;track quality;tracks", 7, -1., 6.);
3467  }
3468 
3469  return h;
3470 }
3471 
3472 //*************************************************************
3473 // Generic booker function
3474 //*************************************************************
3476  unsigned int theNOfBins,
3477  PVValHelper::residualType resType,
3479  bool isNormalized) {
3480  TH1F::SetDefaultSumw2(kTRUE);
3481 
3482  auto hash = std::make_pair(resType, varType);
3483 
3484  double down = theDetails_.range[hash].first;
3485  double up = theDetails_.range[hash].second;
3486 
3487  if (isNormalized) {
3488  up = up / 100.;
3489  down = down / 100.;
3490  }
3491 
3492  std::vector<TH1F*> h;
3493  h.reserve(theNOfBins);
3494 
3495  if (theNOfBins == 0) {
3496  edm::LogError("PrimaryVertexValidation")
3497  << "bookResidualsHistogram() The number of bins cannot be identically 0" << std::endl;
3498  assert(false);
3499  }
3500 
3501  std::string s_resType = std::get<0>(PVValHelper::getTypeString(resType));
3502  std::string s_varType = std::get<0>(PVValHelper::getVarString(varType));
3503 
3504  std::string t_resType = std::get<1>(PVValHelper::getTypeString(resType));
3505  std::string t_varType = std::get<1>(PVValHelper::getVarString(varType));
3506  std::string units = std::get<2>(PVValHelper::getTypeString(resType));
3507 
3508  for (unsigned int i = 0; i < theNOfBins; i++) {
3510  ? Form("%s vs %s - bin %i (%f < %s < %f);%s %s;tracks",
3511  t_resType.c_str(),
3512  t_varType.c_str(),
3513  i,
3515  t_varType.c_str(),
3517  t_resType.c_str(),
3518  units.c_str())
3519  : Form("%s vs %s - bin %i;%s %s;tracks",
3520  t_resType.c_str(),
3521  t_varType.c_str(),
3522  i,
3523  t_resType.c_str(),
3524  units.c_str());
3525 
3526  TH1F* htemp = dir.make<TH1F>(
3527  Form("histo_%s_%s_plot%i", s_resType.c_str(), s_varType.c_str(), i),
3528  //Form("%s vs %s - bin %i;%s %s;tracks",t_resType.c_str(),t_varType.c_str(),i,t_resType.c_str(),units.c_str()),
3529  title.Data(),
3531  down,
3532  up);
3533  h.push_back(htemp);
3534  }
3535 
3536  return h;
3537 }
3538 
3539 //*************************************************************
3540 void PrimaryVertexValidation::fillTrackHistos(std::map<std::string, TH1*>& h,
3541  const std::string& ttype,
3542  const reco::TransientTrack* tt,
3543  const reco::Vertex& v,
3544  const reco::BeamSpot& beamSpot,
3545  double fBfield_)
3546 //*************************************************************
3547 {
3548  using namespace reco;
3549 
3550  PVValHelper::fill(h, "pseudorapidity_" + ttype, tt->track().eta());
3551  PVValHelper::fill(h, "z0_" + ttype, tt->track().vz());
3552  PVValHelper::fill(h, "phi_" + ttype, tt->track().phi());
3553  PVValHelper::fill(h, "eta_" + ttype, tt->track().eta());
3554  PVValHelper::fill(h, "pt_" + ttype, tt->track().pt());
3555  PVValHelper::fill(h, "p_" + ttype, tt->track().p());
3556  PVValHelper::fill(h, "found_" + ttype, tt->track().found());
3557  PVValHelper::fill(h, "lost_" + ttype, tt->track().lost());
3558  PVValHelper::fill(h, "nchi2_" + ttype, tt->track().normalizedChi2());
3559  PVValHelper::fill(h, "rstart_" + ttype, (tt->track().innerPosition()).Rho());
3560 
3561  double d0Error = tt->track().d0Error();
3562  double d0 = tt->track().dxy(beamSpot.position());
3563  double dz = tt->track().dz(beamSpot.position());
3564  if (d0Error > 0) {
3565  PVValHelper::fill(h, "logtresxy_" + ttype, log(d0Error / 0.0001) / log(10.));
3566  PVValHelper::fill(h, "tpullxy_" + ttype, d0 / d0Error);
3567  PVValHelper::fill(h, "tlogDCAxy_" + ttype, log(std::abs(d0 / d0Error)));
3568  }
3569  //double z0=tt->track().vz();
3570  double dzError = tt->track().dzError();
3571  if (dzError > 0) {
3572  PVValHelper::fill(h, "logtresz_" + ttype, log(dzError / 0.0001) / log(10.));
3573  PVValHelper::fill(h, "tpullz_" + ttype, dz / dzError);
3574  PVValHelper::fill(h, "tlogDCAz_" + ttype, log(std::abs(dz / dzError)));
3575  }
3576 
3577  //
3578  double wxy2_ = pow(beamSpot.BeamWidthX(), 2) + pow(beamSpot.BeamWidthY(), 2);
3579 
3581  h, "sigmatrkz_" + ttype, sqrt(pow(tt->track().dzError(), 2) + wxy2_ / pow(tan(tt->track().theta()), 2)));
3582  PVValHelper::fill(h, "sigmatrkz0_" + ttype, tt->track().dzError());
3583 
3584  // track vs vertex
3585  if (v.isValid()) { // && (v.ndof()<10.)) {
3586  // emulate clusterizer input
3587  //const TransientTrack & tt = theB_->build(&t); wrong !!!!
3588  //reco::TransientTrack tt = theB_->build(&t);
3589  //ttt->track().setBeamSpot(beamSpot); // need the setBeamSpot !
3590  double z = (tt->stateAtBeamLine().trackStateAtPCA()).position().z();
3591  double tantheta = tan((tt->stateAtBeamLine().trackStateAtPCA()).momentum().theta());
3592  double dz2 = pow(tt->track().dzError(), 2) + wxy2_ / pow(tantheta, 2);
3593 
3594  PVValHelper::fill(h, "restrkz_" + ttype, z - v.position().z());
3595  PVValHelper::fill(h, "restrkzvsphi_" + ttype, tt->track().phi(), z - v.position().z());
3596  PVValHelper::fill(h, "restrkzvseta_" + ttype, tt->track().eta(), z - v.position().z());
3597  PVValHelper::fill(h, "pulltrkzvsphi_" + ttype, tt->track().phi(), (z - v.position().z()) / sqrt(dz2));
3598  PVValHelper::fill(h, "pulltrkzvseta_" + ttype, tt->track().eta(), (z - v.position().z()) / sqrt(dz2));
3599 
3600  PVValHelper::fill(h, "pulltrkz_" + ttype, (z - v.position().z()) / sqrt(dz2));
3601 
3602  double x1 = tt->track().vx() - beamSpot.x0();
3603  double y1 = tt->track().vy() - beamSpot.y0();
3604 
3605  double kappa = -0.002998 * fBfield_ * tt->track().qoverp() / cos(tt->track().theta());
3606  double D0 = x1 * sin(tt->track().phi()) - y1 * cos(tt->track().phi()) - 0.5 * kappa * (x1 * x1 + y1 * y1);
3607  double q = sqrt(1. - 2. * kappa * D0);
3608  double s0 = (x1 * cos(tt->track().phi()) + y1 * sin(tt->track().phi())) / q;
3609  // double s1;
3610  if (std::abs(kappa * s0) > 0.001) {
3611  //s1=asin(kappa*s0)/kappa;
3612  } else {
3613  //double ks02=(kappa*s0)*(kappa*s0);
3614  //s1=s0*(1.+ks02/6.+3./40.*ks02*ks02+5./112.*pow(ks02,3));
3615  }
3616  // sp.ddcap=-2.*D0/(1.+q);
3617  //double zdcap=tt->track().vz()-s1/tan(tt->track().theta());
3618  }
3619  //
3620 
3621  // collect some info on hits and clusters
3622  PVValHelper::fill(h, "nbarrelLayers_" + ttype, tt->track().hitPattern().pixelBarrelLayersWithMeasurement());
3623  PVValHelper::fill(h, "nPxLayers_" + ttype, tt->track().hitPattern().pixelLayersWithMeasurement());
3624  PVValHelper::fill(h, "nSiLayers_" + ttype, tt->track().hitPattern().trackerLayersWithMeasurement());
3626  h, "expectedInner_" + ttype, tt->track().hitPattern().numberOfLostHits(HitPattern::MISSING_INNER_HITS));
3628  h, "expectedOuter_" + ttype, tt->track().hitPattern().numberOfLostHits(HitPattern::MISSING_OUTER_HITS));
3629  PVValHelper::fill(h, "trackAlgo_" + ttype, tt->track().algo());
3630  PVValHelper::fill(h, "trackQuality_" + ttype, tt->track().qualityMask());
3631 
3632  //
3633  int longesthit = 0, nbarrel = 0;
3634  for (auto const& hit : tt->track().recHits()) {
3635  if (hit->isValid() && hit->geographicalId().det() == DetId::Tracker) {
3636  bool barrel = DetId(hit->geographicalId()).subdetId() == static_cast<int>(PixelSubdetector::PixelBarrel);
3637  //bool endcap = DetId::DetId(hit->geographicalId()).subdetId() == static_cast<int>(PixelSubdetector::PixelEndcap);
3638  if (barrel) {
3639  const SiPixelRecHit* pixhit = dynamic_cast<const SiPixelRecHit*>(&(*hit));
3640  edm::Ref<edmNew::DetSetVector<SiPixelCluster>, SiPixelCluster> const& clust = (*pixhit).cluster();
3641  if (clust.isNonnull()) {
3642  nbarrel++;
3643  if (clust->sizeY() - longesthit > 0)
3644  longesthit = clust->sizeY();
3645  if (clust->sizeY() > 20.) {
3646  PVValHelper::fill(h, "lvseta_" + ttype, tt->track().eta(), 19.9);
3647  PVValHelper::fill(h, "lvstanlambda_" + ttype, tan(tt->track().lambda()), 19.9);
3648  } else {
3649  PVValHelper::fill(h, "lvseta_" + ttype, tt->track().eta(), float(clust->sizeY()));
3650  PVValHelper::fill(h, "lvstanlambda_" + ttype, tan(tt->track().lambda()), float(clust->sizeY()));
3651  }
3652  }
3653  }
3654  }
3655  }
3656  PVValHelper::fill(h, "nbarrelhits_" + ttype, float(nbarrel));
3657  //-------------------------------------------------------------------
3658 }
3659 
3662  desc.setComment("Validates alignment payloads by evaluating unbiased track paramter resisuals to vertices");
3663 
3664  // PV Validation specific
3665 
3666  desc.addUntracked<int>("compressionSettings", -1);
3667  desc.add<bool>("storeNtuple", false);
3668  desc.add<bool>("isLightNtuple", true);
3669  desc.add<bool>("useTracksFromRecoVtx", false);
3670  desc.addUntracked<double>("vertexZMax", 99);
3671  desc.addUntracked<double>("intLumi", 0.);
3672  desc.add<bool>("askFirstLayerHit", false);
3673  desc.addUntracked<bool>("doBPix", true);
3674  desc.addUntracked<bool>("doFPix", true);
3675  desc.addUntracked<double>("probePt", 0.);
3676  desc.addUntracked<double>("probeP", 0.);
3677  desc.addUntracked<double>("probeEta", 2.4);
3678  desc.addUntracked<double>("probeNHits", 0.);
3679  desc.addUntracked<int>("numberOfBins", 24);
3680  desc.addUntracked<double>("minPt", 1.);
3681  desc.addUntracked<double>("maxPt", 20.);
3682  desc.add<bool>("Debug", false);
3683  desc.addUntracked<bool>("runControl", false);
3684  desc.addUntracked<bool>("forceBeamSpot", false);
3685 
3686  std::vector<unsigned int> defaultRuns;
3687  defaultRuns.push_back(0);
3688  desc.addUntracked<std::vector<unsigned int>>("runControlNumber", defaultRuns);
3689 
3690  // event sources
3691 
3692  desc.add<edm::InputTag>("TrackCollectionTag", edm::InputTag("ALCARECOTkAlMinBias"));
3693  desc.add<edm::InputTag>("VertexCollectionTag", edm::InputTag("offlinePrimaryVertices"));
3694  desc.add<edm::InputTag>("BeamSpotTag", edm::InputTag("offlineBeamSpot"));
3695 
3696  // track filtering
3700  desc.add<edm::ParameterSetDescription>("TkFilterParameters", psd0);
3701 
3702  // PV Clusterization
3703  {
3705  {
3708  psd0.add<edm::ParameterSetDescription>("TkDAClusParameters", psd1);
3709 
3712  psd0.add<edm::ParameterSetDescription>("TkGapClusParameters", psd2);
3713  }
3714  psd0.add<std::string>("algorithm", "DA_vect");
3715  desc.add<edm::ParameterSetDescription>("TkClusParameters", psd0);
3716  }
3717 
3718  descriptions.add("primaryVertexValidation", desc);
3719 }
3720 
3721 //define this as a plug-in
bool quality(const TrackQuality) const
Track quality.
Definition: TrackBase.h:552
std::vector< TH1F * > h_norm_dxy_modZ_
TrajectoryStateClosestToPoint trajectoryStateClosestToPoint(const AlgebraicVector3 &momentum, const GlobalPoint &referencePoint, const TrackCharge &charge, const AlgebraicSymMatrix66 &theCovarianceMatrix, const MagneticField *field)
static const char runNumber_[]
static const std::string kSharedResource
Definition: TFileService.h:76
Definition: BitonicSort.h:7
Log< level::Info, true > LogVerbatim
void fillTrackHistos(std::map< std::string, TH1 *> &h, const std::string &ttype, const reco::TransientTrack *tt, const reco::Vertex &v, const reco::BeamSpot &beamSpot, double fBfield)
double qoverp() const
q / p
Definition: TrackBase.h:599
edm::EDGetTokenT< reco::TrackCollection > theTrackCollectionToken_
std::vector< TH1F * > n_IP3DPhiResiduals
Measurement1D getMedian(TH1F *histo)
std::vector< TH1F * > a_dxyEtaResiduals
TH1F * n_dzResidualsMap[nMaxBins_][nMaxBins_]
std::vector< TH1F * > a_d3DEtaResiduals
ESGetTokenH3DDVariant esConsumes(std::string const &Record, edm::ConsumesCollector &)
Definition: DeDxTools.cc:283
T getParameter(std::string const &) const
Definition: ParameterSet.h:307
unsigned int pxbLayer(const DetId &id) const
std::vector< TH1F * > n_reszPhiResiduals
std::vector< unsigned int > runControlNumbers_
T const & getData(const ESGetToken< T, R > &iToken) const noexcept(false)
Definition: EventSetup.h:119
const bool isValid(const Frame &aFrame, const FrameQuality &aQuality, const uint16_t aExpectedPos)
static bool vtxSort(const reco::Vertex &a, const reco::Vertex &b)
static void fillPSetDescription(edm::ParameterSetDescription &desc)
unsigned short numberOfValidHits() const
number of valid hits found
Definition: TrackBase.h:798
std::pair< Measurement1D, Measurement1D > fitResiduals(TH1 *hist)
Quality qualityByName(std::string const &name)
edm::EDGetTokenT< reco::BeamSpot > theBeamspotToken_
void fillTrendPlot(TH1F *trendPlot, TH1F *residualsPlot[100], PVValHelper::estimator fitPar_, const std::string &var_)
std::vector< TH1F * > h_dxy_pT_
TH1F * a_dxyResidualsMap[nMaxBins_][nMaxBins_]
std::map< std::string, TH1 * > hDA
TrackQuality
track quality
Definition: TrackBase.h:150
std::vector< TH1F * > a_IP3DEtaResiduals
Common base class.
std::vector< float > generateBins(int n, float start, float range)
void beginRun(edm::Run const &iRun, edm::EventSetup const &iSetup) override
double px() const
x coordinate of momentum vector
Definition: TrackBase.h:640
std::vector< TH1F * > n_dzPhiResiduals
double p() const
momentum vector magnitude
Definition: TrackBase.h:631
Divides< arg, void > D0
Definition: Factorize.h:135
std::vector< TH1F * > a_d3DPhiResiduals
std::vector< TH1F * > n_d3DPhiResiduals
PrimaryVertexValidation(const edm::ParameterSet &)
std::vector< TH1F * > n_dxyPhiBiasResiduals
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
std::pair< bool, Measurement1D > signedTransverseImpactParameter(const reco::TransientTrack &track, const GlobalVector &direction, const reco::Vertex &vertex)
Definition: IPTools.cc:57
TH1F * n_dxyBiasResidualsMap[nMaxBins_][nMaxBins_]
std::vector< TH1F * > n_d3DEtaResiduals
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
const edm::ESGetToken< GlobalTrackingGeometry, GlobalTrackingGeometryRecord > trackingGeomToken_
std::vector< TH1F * > h_norm_dz_Central_pT_
std::pair< bool, Measurement1D > absoluteImpactParameter3D(const reco::TransientTrack &transientTrack, const reco::Vertex &vertex)
Definition: IPTools.cc:38
const unsigned getPXBModules(unsigned int lay) const
double py() const
y coordinate of momentum vector
Definition: TrackBase.h:643
std::pair< bool, Measurement1D > signedImpactParameter3D(const reco::TransientTrack &track, const GlobalVector &direction, const reco::Vertex &vertex)
Definition: IPTools.cc:81
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:238
double longitudinalImpactParameterError() const
float DOFUnbiasedVertex_[nMaxtracks_]
std::map< std::pair< residualType, plotVariable >, std::pair< float, float > > range
double d3DFromMyVertex_[nMaxtracks_]
const DetContainer & dets() const override
Returm a vector of all GeomDet (including all GeomDetUnits)
std::vector< TH1F * > a_reszPhiResiduals
std::vector< TH1F * > a_IP2DPhiResiduals
float chi2ProbUnbiasedVertex_[nMaxtracks_]
unsigned int pxbLadder(const DetId &id) const
const PerigeeTrajectoryError & perigeeError() const
Log< level::Error, false > LogError
std::vector< TH1F * > n_reszEtaResiduals
double dxyFromMyVertex_[nMaxtracks_]
assert(be >=bs)
constexpr double max_eta_phase2
Measurement1D getMAD(TH1F *histo)
std::map< std::string, TH1 * > bookVertexHistograms(const TFileDirectory &dir)
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
std::vector< TH1F * > h_dz_modZ_
std::vector< TH1F * > h_dxy_modZ_
std::vector< TH1F * > a_dxEtaResiduals
std::pair< bool, bool > pixelHitsCheck(const reco::TransientTrack &track)
plotLabels getVarString(plotVariable var)
double dzErrorFromMyVertex_[nMaxtracks_]
static std::string to_string(const XMLCh *ch)
std::vector< TH1F * > n_dxyEtaBiasResiduals
#define LogTrace(id)
std::vector< TH1F * > n_dxyEtaResiduals
double IPTsigFromMyVertex_[nMaxtracks_]
T getUntrackedParameter(std::string const &, T const &) const
reco::TransientTrack build(const reco::Track *p) const
GlobalPoint position() const
std::vector< TH1F * > a_dzPhiBiasResiduals
std::map< plotVariable, std::vector< float > > trendbins
double pt() const
track transverse momentum
Definition: TrackBase.h:637
T * make(const Args &...args) const
make new ROOT object
constexpr double max_eta_phase1
auto recHits() const
Access to reconstructed hits on the track.
Definition: Track.h:85
std::vector< TH1F * > h_norm_dxy_pT_
std::vector< TH1F * > a_dzEtaResiduals
T x() const
Definition: PV3DBase.h:59
T y() const
Definition: PV3DBase.h:60
edm::EDGetTokenT< reco::VertexCollection > theVertexCollectionToken_
float getHigh(residualType type, plotVariable plot)
int tracksUsedForVertexing_[nMaxtracks_]
float getLow(residualType type, plotVariable plot)
int charge() const
track electric charge
Definition: TrackBase.h:596
double dz() const
dz parameter (= dsz/cos(lambda)). This is the track z0 w.r.t (0,0,0) only if the refPoint is close to...
Definition: TrackBase.h:622
int iEvent
Definition: GenABIO.cc:224
T const * product() const
Definition: ESHandle.h:86
Definition: TTTypes.h:54
RunNumber_t run() const
Definition: RunBase.h:40
runControlNumber
Definition: PV_cfg.py:170
std::vector< TH1F * > a_dxyEtaBiasResiduals
const edm::ESGetToken< RunInfo, RunInfoRcd > runInfoTokenBR_
std::vector< TH1F * > h_norm_dxy_ladder_
std::vector< TH1F * > a_dzPhiResiduals
double dzError() const
error on dz
Definition: TrackBase.h:778
double yUnbiasedVertex_[nMaxtracks_]
bool isThere(GeomDetEnumerators::SubDetector subdet) const
std::unique_ptr< TrackClusterizerInZ > theTrackClusterizer_
T sqrt(T t)
Definition: SSEVec.h:19
static void fillPSetDescription(edm::ParameterSetDescription &desc)
std::vector< TH1F * > n_dxyPhiResiduals
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
TH1F * a_dzL1ResidualsMap[nMaxBins_][nMaxBins_]
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
double zUnbiasedVertex_[nMaxtracks_]
Transition
Definition: Transition.h:12
double phi() const
azimuthal angle of momentum vector
Definition: TrackBase.h:649
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
std::vector< TH1F * > h_dxy_ladder_
void fill(std::map< std::string, TH1 *> &h, const std::string &s, double x)
float chi2normUnbiasedVertex_[nMaxtracks_]
const Point & vertex() const
reference point on the track. This method is DEPRECATED, please use referencePoint() instead ...
Definition: TrackBase.h:676
TH1F * n_d3DResidualsMap[nMaxBins_][nMaxBins_]
std::vector< TH1F * > a_dxyPhiResiduals
std::vector< TH1F * > a_dxyPhiBiasResiduals
void analyze(const edm::Event &, const edm::EventSetup &) override
TH1F * a_d3DResidualsMap[nMaxBins_][nMaxBins_]
ParameterDescriptionBase * add(U const &iLabel, T const &value)
float clusterProbability(unsigned int flags=0) const
Definition: SiPixelRecHit.cc:9
std::vector< TH1F * > h_dxy_ladderNoOverlap_
ESHandle< T > getHandle(const ESGetToken< T, R > &iToken) const
Definition: EventSetup.h:130
TH1F * n_dxyL1ResidualsMap[nMaxBins_][nMaxBins_]
double IPLsigFromMyVertex_[nMaxtracks_]
void fillMap(TH2F *trendMap, TH1F *residualsMapPlot[100][100], PVValHelper::estimator fitPar_, const int nXBins_, const int nYBins_)
const PerigeeTrajectoryParameters & perigeeParameters() const
const edm::ESGetToken< TrackerTopology, TrackerTopologyRcd > topoToken_
bool hasFirstLayerPixelHits(const reco::TransientTrack &track)
std::vector< TH1F * > n_dzEtaResiduals
void setMap(residualType type, plotVariable plot, float low, float high)
std::unique_ptr< TrackFilterForPVFindingBase > theTrackFilter_
#define M_PI
std::vector< TH1F * > n_IP3DEtaResiduals
Log< level::Info, false > LogInfo
std::vector< TH1F * > h_norm_dz_pT_
bool isHit2D(const TrackingRecHit &hit, const PVValHelper::detectorPhase &thePhase) const
double eta() const
pseudorapidity of momentum vector
Definition: TrackBase.h:652
Definition: DetId.h:17
bool isBFieldConsistentWithMode(const edm::EventSetup &iSetup) const
static constexpr float d0
TH1F * a_dxyL1ResidualsMap[nMaxBins_][nMaxBins_]
void shrinkHistVectorToFit(std::vector< TH1F *> &h, unsigned int desired_size)
static TrackQuality qualityByName(const std::string &name)
Definition: TrackBase.cc:126
std::vector< TH1F * > h_dz_pT_
const edm::ESGetToken< TransientTrackBuilder, TransientTrackRecord > ttkToken_
plotLabels getTypeString(residualType type)
TFileDirectory mkdir(const std::string &dir, const std::string &descr="")
create a new subdirectory
Definition: TFileService.h:69
const HitPattern & hitPattern() const
Access the hit pattern, indicating in which Tracker layers the track has hits.
Definition: TrackBase.h:504
std::vector< TH1F * > n_dzEtaBiasResiduals
edm::Service< TFileService > fs
std::vector< TH1F * > h_dz_Central_pT_
static void fillPSetDescription(edm::ParameterSetDescription &desc)
std::vector< TH1F * > h_dz_ladder_
std::vector< TH1F * > a_IP2DEtaResiduals
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
double pz() const
z coordinate of momentum vector
Definition: TrackBase.h:646
double chi2() const
chi-squared of the fit
Definition: TrackBase.h:587
double normalizedChi2() const
chi-squared divided by n.d.o.f. (or chi-squared * 1e6 if n.d.o.f. is zero)
Definition: TrackBase.h:593
constexpr double max_eta_phase0
double b
Definition: hdecay.h:120
const FreeTrajectoryState & theState() const
std::vector< TH1F * > n_IP2DPhiResiduals
void add(std::string const &label, ParameterSetDescription const &psetDescription)
double xUnbiasedVertex_[nMaxtracks_]
std::vector< TH1F * > h_dxy_Central_pT_
PVValHelper::detectorPhase phase_
bool isValid() const
Definition: HandleBase.h:70
double value() const
Definition: Measurement1D.h:25
TH1F * n_dzL1ResidualsMap[nMaxBins_][nMaxBins_]
double theta() const
polar angle
Definition: TrackBase.h:602
std::map< unsigned int, std::pair< long long, long long > > runNumbersTimesLog_
~PrimaryVertexValidation() override
Pixel cluster – collection of neighboring pixels above threshold.
TString units(TString variable, Char_t axis)
double dzFromMyVertex_[nMaxtracks_]
fixed size matrix
HLT enums.
double a
Definition: hdecay.h:121
static int position[264][3]
Definition: ReadPGInfo.cc:289
std::vector< TH1F * > a_reszEtaResiduals
std::vector< TH1F * > bookResidualsHistogram(const TFileDirectory &dir, unsigned int theNOfBins, PVValHelper::residualType resType, PVValHelper::plotVariable varType, bool isNormalized=false)
double error() const
Definition: Measurement1D.h:27
TH1F * a_dxyBiasResidualsMap[nMaxBins_][nMaxBins_]
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
static void fillPSetDescription(edm::ParameterSetDescription &desc)
std::vector< TH1F * > a_dxPhiResiduals
std::vector< TH1F * > h_norm_dz_modZ_
TH1F * n_dxyResidualsMap[nMaxBins_][nMaxBins_]
bool passesTrackCuts(const reco::Track &track, const reco::Vertex &vertex, const std::string &qualityString_, double dxyErrMax_, double dzErrMax_, double ptErrMax_)
const edm::ESGetToken< TrackerTopology, TrackerTopologyRcd > topoTokenBR_
double dxyErrorFromMyVertex_[nMaxtracks_]
std::vector< TH1F * > a_dyEtaResiduals
unsigned int pxbModule(const DetId &id) const
Log< level::Warning, false > LogWarning
std::vector< TH1F * > n_IP2DEtaResiduals
std::vector< TH1F * > h_norm_dz_ladder_
char const * what() const noexcept override
Definition: Exception.cc:107
float sumOfWeightsUnbiasedVertex_[nMaxtracks_]
void fillByIndex(std::vector< TH1F *> &h, unsigned int index, double x, std::string tag="")
TH1F * n_dzBiasResidualsMap[nMaxBins_][nMaxBins_]
double d3DErrorFromMyVertex_[nMaxtracks_]
std::vector< TH1F * > a_dyPhiResiduals
double d0Error() const
error on d0
Definition: TrackBase.h:772
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
std::vector< TH1F * > h_dxy_ladderOverlap_
const edm::ESGetToken< MagneticField, IdealMagneticFieldRecord > magFieldToken_
TH1F * a_dzBiasResidualsMap[nMaxBins_][nMaxBins_]
PVValHelper::histodetails theDetails_
std::array< float, nPtBins_+1 > mypT_bins_
TH1F * a_dzResidualsMap[nMaxBins_][nMaxBins_]
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
std::vector< TrackBaseRef >::const_iterator trackRef_iterator
The iteratator for the vector<TrackRef>
Definition: Vertex.h:38
float chi2UnbiasedVertex_[nMaxtracks_]
std::pair< long long, long long > getRunTime(const edm::EventSetup &iSetup) const
*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 unsigned getPXBLadders(unsigned int lay) const
std::vector< TH1F * > h_norm_dxy_Central_pT_
Definition: Run.h:45
Global3DVector GlobalVector
Definition: GlobalVector.h:10
void fillTrendPlotByIndex(TH1F *trendPlot, std::vector< TH1F *> &h, PVValHelper::estimator fitPar_, PVValHelper::plotVariable plotVar=PVValHelper::END_OF_PLOTS)
std::vector< TH1F * > a_dzEtaBiasResiduals
Our base class.
Definition: SiPixelRecHit.h:23
std::vector< TH1F * > a_IP3DPhiResiduals
double dxy() const
dxy parameter. (This is the transverse impact parameter w.r.t. to (0,0,0) ONLY if refPoint is close t...
Definition: TrackBase.h:608
#define LogDebug(id)
std::vector< TH1F * > n_dzPhiBiasResiduals
double IP3DsigFromMyVertex_[nMaxtracks_]
TFile & file() const
return opened TFile
Definition: TFileService.h:37
const edm::ESGetToken< TrackerGeometry, TrackerDigiGeometryRecord > geomTokenBR_