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