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