CMS 3D CMS Logo
PrimaryVertexValidation.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: Alignment/OfflineValidation
4 // Class: PrimaryVertexValidation
5 //
13 //
14 // Original Author: Marco Musich
15 // Created: Tue Mar 02 10:39:34 CDT 2010
16 //
17 
18 // system include files
19 #include <memory>
20 #include <vector>
21 #include <regex>
22 #include <cassert>
23 
24 // user include files
25 #include "Alignment/OfflineValidation/plugins/PrimaryVertexValidation.h"
26 
27 // ROOT includes
28 #include "TH1F.h"
29 #include "TH2F.h"
30 #include "TF1.h"
31 #include "TVector3.h"
32 #include "TFile.h"
33 #include "TMath.h"
34 #include "TROOT.h"
35 #include "TChain.h"
36 #include "TNtuple.h"
37 #include "TMatrixD.h"
38 #include "TVectorD.h"
39 
40 // CMSSW includes
41 #include "CondFormats/RunInfo/interface/RunInfo.h"
42 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
43 #include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
44 #include "DataFormats/TrackReco/interface/Track.h"
45 #include "DataFormats/TrackReco/interface/TrackFwd.h"
46 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
47 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
48 #include "FWCore/MessageLogger/interface/MessageLogger.h"
49 #include "FWCore/ParameterSet/interface/ParameterSet.h"
50 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
51 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
52 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
53 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
54 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
55 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
56 #include "MagneticField/Engine/interface/MagneticField.h"
57 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.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/Records/interface/TransientTrackRecord.h"
63 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
64 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateClosestToPoint.h"
65 
66 const int PrimaryVertexValidation::nMaxtracks_;
67 
68 // Constructor
69 PrimaryVertexValidation::PrimaryVertexValidation(const edm::ParameterSet& iConfig):
70  storeNtuple_(iConfig.getParameter<bool>("storeNtuple")),
71  lightNtupleSwitch_(iConfig.getParameter<bool>("isLightNtuple")),
72  useTracksFromRecoVtx_(iConfig.getParameter<bool>("useTracksFromRecoVtx")),
73  vertexZMax_(iConfig.getUntrackedParameter<double>("vertexZMax",99.)),
74  intLumi_(iConfig.getUntrackedParameter<double>("intLumi",0.)),
75  askFirstLayerHit_(iConfig.getParameter<bool>("askFirstLayerHit")),
76  doBPix_(iConfig.getUntrackedParameter<bool>("doBPix",true)),
77  doFPix_(iConfig.getUntrackedParameter<bool>("doFPix",true)),
78  ptOfProbe_(iConfig.getUntrackedParameter<double>("probePt",0.)),
79  pOfProbe_(iConfig.getUntrackedParameter<double>("probeP",0.)),
80  etaOfProbe_(iConfig.getUntrackedParameter<double>("probeEta",2.4)),
81  nHitsOfProbe_(iConfig.getUntrackedParameter<double>("probeNHits",0.)),
82  nBins_(iConfig.getUntrackedParameter<int>("numberOfBins",24)),
83  minPt_(iConfig.getUntrackedParameter<double>("minPt",1.)),
84  maxPt_(iConfig.getUntrackedParameter<double>("maxPt",20.)),
85  debug_(iConfig.getParameter<bool>("Debug")),
86  runControl_(iConfig.getUntrackedParameter<bool>("runControl",false))
87 {
88 
89  // now do what ever initialization is needed
90  // initialize phase space boundaries
91 
92  usesResource(TFileService::kSharedResource);
93 
94  std::vector<unsigned int> defaultRuns;
95  defaultRuns.push_back(0);
96  runControlNumbers_ = iConfig.getUntrackedParameter<std::vector<unsigned int> >("runControlNumber",defaultRuns);
97 
98  edm::InputTag TrackCollectionTag_ = iConfig.getParameter<edm::InputTag>("TrackCollectionTag");
99  theTrackCollectionToken = consumes<reco::TrackCollection>(TrackCollectionTag_);
100 
101  edm::InputTag VertexCollectionTag_ = iConfig.getParameter<edm::InputTag>("VertexCollectionTag");
102  theVertexCollectionToken = consumes<reco::VertexCollection>(VertexCollectionTag_);
103 
104  edm::InputTag BeamspotTag_ = edm::InputTag("offlineBeamSpot");
105  theBeamspotToken = consumes<reco::BeamSpot>(BeamspotTag_);
106 
107  // select and configure the track filter
108  theTrackFilter_= new TrackFilterForPVFinding(iConfig.getParameter<edm::ParameterSet>("TkFilterParameters") );
109  // select and configure the track clusterizer
110  std::string clusteringAlgorithm=iConfig.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<std::string>("algorithm");
111  if (clusteringAlgorithm=="gap"){
112  theTrackClusterizer_ = new GapClusterizerInZ(iConfig.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<edm::ParameterSet>("TkGapClusParameters"));
113  }else if(clusteringAlgorithm=="DA"){
114  theTrackClusterizer_ = new DAClusterizerInZ(iConfig.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<edm::ParameterSet>("TkDAClusParameters"));
115  // provide the vectorized version of the clusterizer, if supported by the build
116  } else if(clusteringAlgorithm=="DA_vect") {
117  theTrackClusterizer_ = new DAClusterizerInZ_vect(iConfig.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<edm::ParameterSet>("TkDAClusParameters"));
118  } else {
119  throw VertexException("PrimaryVertexProducerAlgorithm: unknown clustering algorithm: " + clusteringAlgorithm);
120  }
121 
122  theDetails_.histobins = 500;
123  theDetails_.setMap(PVValHelper::dxy,PVValHelper::phi ,-2000.,2000.);
124  theDetails_.setMap(PVValHelper::dxy,PVValHelper::eta ,-3000.,3000.);
125  theDetails_.setMap(PVValHelper::dxy,PVValHelper::pT ,-1000.,1000.);
126  theDetails_.setMap(PVValHelper::dxy,PVValHelper::pTCentral,-1000.,1000.);
127  theDetails_.setMap(PVValHelper::dxy,PVValHelper::ladder ,-1000.,1000.);
128  theDetails_.setMap(PVValHelper::dxy,PVValHelper::modZ ,-1000.,1000.);
129 
130  for (int i = PVValHelper::phi; i < PVValHelper::END_OF_PLOTS; i++ ){
131  for (int j = PVValHelper::dx; j < PVValHelper::END_OF_TYPES; j++ ){
132 
133  auto plot_index = static_cast<PVValHelper::plotVariable>(i);
134  auto res_index = static_cast<PVValHelper::residualType>(j);
135 
136  if(debug_){
137  edm::LogInfo("PrimaryVertexValidation")<<"==> "<<std::get<0>(PVValHelper::getTypeString(res_index)) << " "<< std::setw(10)<< std::get<0>(PVValHelper::getVarString(plot_index))<<std::endl;
138  }
139  if(res_index!=PVValHelper::d3D && res_index!=PVValHelper::norm_d3D)
140  theDetails_.setMap(res_index,plot_index,theDetails_.getLow(PVValHelper::dxy,plot_index),theDetails_.getHigh(PVValHelper::dxy,plot_index));
141  else
142  theDetails_.setMap(res_index,plot_index,0.,theDetails_.getHigh(PVValHelper::dxy,plot_index));
143  }
144  }
145 
146  edm::LogVerbatim("PrimaryVertexValidation") <<"######################################";
147  for (const auto & it : theDetails_.range){
148  edm::LogVerbatim("PrimaryVertexValidation")<< "|" <<std::setw(10) << std::get<0>(PVValHelper::getTypeString(it.first.first)) << "|" << std::setw(10)<< std::get<0>(PVValHelper::getVarString(it.first.second)) << "| (" << std::setw(5)<< it.second.first << ";" <<std::setw(5)<< it.second.second << ") |"<<std::endl;
149  }
150 
151  theDetails_.trendbins[PVValHelper::phi] = PVValHelper::generateBins(nBins_+1,-180.,360.);
152  theDetails_.trendbins[PVValHelper::eta] = PVValHelper::generateBins(nBins_+1,-etaOfProbe_,2*etaOfProbe_);
153 
154  if(debug_){
155  edm::LogVerbatim("PrimaryVertexValidation") << "etaBins: ";
156  for (auto ieta: theDetails_.trendbins[PVValHelper::eta]) {
157  edm::LogVerbatim("PrimaryVertexValidation") << ieta << " ";
158  }
159  edm::LogVerbatim("PrimaryVertexValidation") << "\n";
160 
161  edm::LogVerbatim("PrimaryVertexValidation") << "phiBins: ";
162  for (auto iphi: theDetails_.trendbins[PVValHelper::phi]) {
163  edm::LogVerbatim("PrimaryVertexValidation") << iphi << " ";
164  }
165  edm::LogVerbatim("PrimaryVertexValidation") << "\n";
166  }
167 
168  // create the bins of the pT-binned distributions
169 
170  mypT_bins_ = PVValHelper::makeLogBins<float,nPtBins_>(minPt_,maxPt_);
171 
172  std::string toOutput="";
173  for (auto ptbin: mypT_bins_){
174  toOutput+=" ";
175  toOutput+=std::to_string(ptbin);
176  toOutput+=",";
177  }
178 
179  edm::LogVerbatim("PrimaryVertexValidation") <<"######################################\n";
180  edm::LogVerbatim("PrimaryVertexValidation") <<"The pT binning is: [" << toOutput << "] \n";
181 
182 }
183 
184 // Destructor
185 PrimaryVertexValidation::~PrimaryVertexValidation()
186 {
187  // do anything here that needs to be done at desctruction time
188  // (e.g. close files, deallocate resources etc.)
189  if (theTrackFilter_) delete theTrackFilter_;
190  if (theTrackClusterizer_) delete theTrackClusterizer_;
191 }
192 
193 
194 //
195 // member functions
196 //
197 
198 // ------------ method called to for each event ------------
199 void
200 PrimaryVertexValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
201 {
202 
203  using namespace std;
204  using namespace reco;
205  using namespace IPTools;
206 
207  if (!isBFieldConsistentWithMode(iSetup)) {
208  edm::LogWarning("PrimaryVertexValidation") << "*********************************************************************************\n"
209  << "* The configuration (ptOfProbe > " << ptOfProbe_ << "GeV) is not correctly set for current value of magnetic field \n"
210  << "* Switching it to 0. !!! \n"
211  << "*********************************************************************************"<< std::endl;
212  ptOfProbe_=0.;
213  }
214 
215  if(nBins_!=24 && debug_){
216  edm::LogInfo("PrimaryVertexValidation")<<"Using: "<<nBins_<<" bins plots";
217  }
218 
219  bool passesRunControl = false;
220 
221  if(runControl_){
222  for(const auto & runControlNumber : runControlNumbers_){
223  if(iEvent.eventAuxiliary().run() == runControlNumber){
224  if (debug_){
225  edm::LogInfo("PrimaryVertexValidation")<<" run number: "<<iEvent.eventAuxiliary().run()<<" keeping run:"<<runControlNumber;
226  }
227  passesRunControl = true;
228  break;
229  }
230  }
231  if (!passesRunControl) return;
232  }
233 
234  Nevt_++;
235 
236  //=======================================================
237  // Initialize Root-tuple variables
238  //=======================================================
239 
240  SetVarToZero();
241 
242  //=======================================================
243  // Retrieve the Magnetic Field information
244  //=======================================================
245 
246  edm::ESHandle<MagneticField> theMGField;
247  iSetup.get<IdealMagneticFieldRecord>().get( theMGField );
248 
249  //=======================================================
250  // Retrieve the Tracking Geometry information
251  //=======================================================
252 
253  edm::ESHandle<GlobalTrackingGeometry> theTrackingGeometry;
254  iSetup.get<GlobalTrackingGeometryRecord>().get( theTrackingGeometry );
255 
256  //=======================================================
257  // Retrieve geometry information
258  //=======================================================
259 
260  edm::LogInfo("read tracker geometry...");
261  edm::ESHandle<TrackerGeometry> pDD;
262  iSetup.get<TrackerDigiGeometryRecord>().get( pDD );
263  edm::LogInfo("tracker geometry read")<<"There are: "<< pDD->dets().size() <<" detectors";
264 
265  // switch on the phase1
266  if( (pDD->isThere(GeomDetEnumerators::P1PXB)) ||
267  (pDD->isThere(GeomDetEnumerators::P1PXEC)) ) {
268  isPhase1_ = true;
269  nLadders_ = 12;
270 
271  if(h_dxy_ladderOverlap_.size()!=nLadders_){
272 
273  PVValHelper::shrinkHistVectorToFit(h_dxy_ladderOverlap_,nLadders_);
274  PVValHelper::shrinkHistVectorToFit(h_dxy_ladderNoOverlap_,nLadders_);
275  PVValHelper::shrinkHistVectorToFit(h_dxy_ladder_,nLadders_);
276  PVValHelper::shrinkHistVectorToFit(h_dz_ladder_,nLadders_);
277  PVValHelper::shrinkHistVectorToFit(h_norm_dxy_ladder_,nLadders_);
278  PVValHelper::shrinkHistVectorToFit(h_norm_dz_ladder_,nLadders_);
279 
280  if (debug_){
281  edm::LogInfo("PrimaryVertexValidation")<<"checking size:"<<h_dxy_ladder_.size()<<std::endl;
282  }
283  }
284 
285  if (debug_){
286  edm::LogInfo("PrimaryVertexValidation")<<" pixel phase1 setup, nLadders: "<<nLadders_;
287  }
288 
289  } else {
290  isPhase1_ = false;
291  nLadders_ = 20;
292  if (debug_){
293  edm::LogInfo("PrimaryVertexValidation")<<" pixel phase0 setup, nLadders: "<<nLadders_;
294  }
295  }
296 
297  if(isPhase1_){
298  etaOfProbe_ = std::min(etaOfProbe_,2.7);
299  } else {
300  etaOfProbe_ = std::min(etaOfProbe_,2.5);
301  }
302 
303  if(h_etaMax->GetEntries()==0.){
304  h_etaMax->SetBinContent(1.,etaOfProbe_);
305  h_nbins->SetBinContent(1.,nBins_);
306  h_nLadders->SetBinContent(1.,nLadders_);
307  h_pTinfo->SetBinContent(1.,mypT_bins_.size());
308  h_pTinfo->SetBinContent(2.,minPt_);
309  h_pTinfo->SetBinContent(3.,maxPt_);
310  }
311 
312  //=======================================================
313  // Retrieve the Transient Track Builder information
314  //=======================================================
315 
316  edm::ESHandle<TransientTrackBuilder> theB_;
317  iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder",theB_);
318  double fBfield_=((*theB_).field()->inTesla(GlobalPoint(0.,0.,0.))).z();
319 
320  //=======================================================
321  // Retrieve the Track information
322  //=======================================================
323 
324  edm::Handle<TrackCollection> trackCollectionHandle;
325  iEvent.getByToken(theTrackCollectionToken, trackCollectionHandle);
326  if(!trackCollectionHandle.isValid()) return;
327  auto const & tracks = *trackCollectionHandle;
328 
329  //=======================================================
330  // Retrieve tracker topology from geometry
331  //=======================================================
332 
333  edm::ESHandle<TrackerTopology> tTopoHandle;
334  iSetup.get<TrackerTopologyRcd>().get(tTopoHandle);
335  const TrackerTopology* const tTopo = tTopoHandle.product();
336 
337  //=======================================================
338  // Retrieve offline vartex information (only for reco)
339  //=======================================================
340 
341  //edm::Handle<VertexCollection> vertices;
342  edm::Handle<std::vector<Vertex> > vertices;
343 
344  try {
345  iEvent.getByToken(theVertexCollectionToken, vertices);
346  } catch ( cms::Exception& er ) {
347  LogTrace("PrimaryVertexValidation")<<"caught std::exception "<<er.what()<<std::endl;
348  }
349 
350  std::vector<Vertex> vsorted = *(vertices);
351  // sort the vertices by number of tracks in descending order
352  // use chi2 as tiebreaker
353  std::sort( vsorted.begin(), vsorted.end(), PrimaryVertexValidation::vtxSort );
354 
355  // skip events with no PV, this should not happen
356 
357  if( vsorted.empty()) return;
358 
359  // skip events failing vertex cut
360  if( std::abs(vsorted[0].z()) > vertexZMax_ ) return;
361 
362  if ( vsorted[0].isValid() ) {
363  xOfflineVertex_ = (vsorted)[0].x();
364  yOfflineVertex_ = (vsorted)[0].y();
365  zOfflineVertex_ = (vsorted)[0].z();
366 
367  xErrOfflineVertex_ = (vsorted)[0].xError();
368  yErrOfflineVertex_ = (vsorted)[0].yError();
369  zErrOfflineVertex_ = (vsorted)[0].zError();
370  }
371 
372  h_xOfflineVertex->Fill(xOfflineVertex_);
373  h_yOfflineVertex->Fill(yOfflineVertex_);
374  h_zOfflineVertex->Fill(zOfflineVertex_);
375  h_xErrOfflineVertex->Fill(xErrOfflineVertex_);
376  h_yErrOfflineVertex->Fill(yErrOfflineVertex_);
377  h_zErrOfflineVertex->Fill(zErrOfflineVertex_);
378 
379  unsigned int vertexCollectionSize = vsorted.size();
380  int nvvertex = 0;
381 
382  for (unsigned int i=0; i<vertexCollectionSize; i++) {
383  const Vertex& vertex = vsorted.at(i);
384  if (vertex.isValid()) nvvertex++;
385  }
386 
387  nOfflineVertices_ = nvvertex;
388  h_nOfflineVertices->Fill(nvvertex);
389 
390  if ( !vsorted.empty() && useTracksFromRecoVtx_ ) {
391 
392  double sumpt = 0;
393  size_t ntracks = 0;
394  double chi2ndf = 0.;
395  double chi2prob = 0.;
396 
397  if (!vsorted.at(0).isFake()) {
398 
399  Vertex pv = vsorted.at(0);
400 
401  ntracks = pv.tracksSize();
402  chi2ndf = pv.normalizedChi2();
403  chi2prob = TMath::Prob(pv.chi2(),(int)pv.ndof());
404 
405  h_recoVtxNtracks_->Fill(ntracks);
406  h_recoVtxChi2ndf_->Fill(chi2ndf);
407  h_recoVtxChi2Prob_->Fill(chi2prob);
408 
409  for (Vertex::trackRef_iterator itrk = pv.tracks_begin();itrk != pv.tracks_end(); ++itrk) {
410  double pt = (**itrk).pt();
411  sumpt += pt*pt;
412 
413  const math::XYZPoint myVertex(pv.position().x(),pv.position().y(),pv.position().z());
414 
415  double dxyRes = (**itrk).dxy(myVertex);
416  double dzRes = (**itrk).dz(myVertex);
417 
418  double dxy_err = (**itrk).dxyError();
419  double dz_err = (**itrk).dzError();
420 
421  float trackphi = ((**itrk).phi())*(180/M_PI);
422  float tracketa = (**itrk).eta();
423 
424  for(int i=0; i<nBins_; i++){
425 
426  float phiF = theDetails_.trendbins[PVValHelper::phi][i];
427  float phiL = theDetails_.trendbins[PVValHelper::phi][i+1];
428 
429  float etaF = theDetails_.trendbins[PVValHelper::eta][i];
430  float etaL = theDetails_.trendbins[PVValHelper::eta][i+1];
431 
432  if(tracketa >= etaF && tracketa < etaL ){
433 
434  PVValHelper::fillByIndex(a_dxyEtaBiasResiduals,i,dxyRes*cmToum);
435  PVValHelper::fillByIndex(a_dzEtaBiasResiduals,i,dzRes*cmToum);
436  PVValHelper::fillByIndex(n_dxyEtaBiasResiduals,i,(dxyRes)/dxy_err);
437  PVValHelper::fillByIndex(n_dzEtaBiasResiduals,i,(dzRes)/dz_err);
438 
439  }
440 
441  if(trackphi >= phiF && trackphi < phiL ){
442 
443  PVValHelper::fillByIndex(a_dxyPhiBiasResiduals,i,dxyRes*cmToum);
444  PVValHelper::fillByIndex(a_dzPhiBiasResiduals,i,dzRes*cmToum);
445  PVValHelper::fillByIndex(n_dxyPhiBiasResiduals,i,(dxyRes)/dxy_err);
446  PVValHelper::fillByIndex(n_dzPhiBiasResiduals,i,(dzRes)/dz_err);
447 
448  for(int j=0; j<nBins_; j++){
449 
450  float etaJ = theDetails_.trendbins[PVValHelper::eta][j];
451  float etaK = theDetails_.trendbins[PVValHelper::eta][j+1];
452 
453  if(tracketa >= etaJ && tracketa < etaK ){
454 
455  a_dxyBiasResidualsMap[i][j]->Fill(dxyRes*cmToum);
456  a_dzBiasResidualsMap[i][j]->Fill(dzRes*cmToum);
457 
458  n_dxyBiasResidualsMap[i][j]->Fill((dxyRes)/dxy_err);
459  n_dzBiasResidualsMap[i][j]->Fill((dzRes)/dz_err);
460 
461  }
462  }
463  }
464  }
465  }
466 
467  h_recoVtxSumPt_->Fill(sumpt);
468 
469  }
470  }
471 
472  //=======================================================
473  // Retrieve Beamspot information
474  //=======================================================
475 
476  BeamSpot beamSpot;
477  edm::Handle<BeamSpot> beamSpotHandle;
478  iEvent.getByToken(theBeamspotToken, beamSpotHandle);
479 
480  if ( beamSpotHandle.isValid() ) {
481  beamSpot = *beamSpotHandle;
482  BSx0_ = beamSpot.x0();
483  BSy0_ = beamSpot.y0();
484  BSz0_ = beamSpot.z0();
485  Beamsigmaz_ = beamSpot.sigmaZ();
486  Beamdxdz_ = beamSpot.dxdz();
487  BeamWidthX_ = beamSpot.BeamWidthX();
488  BeamWidthY_ = beamSpot.BeamWidthY();
489 
490  wxy2_=TMath::Power(BeamWidthX_,2)+TMath::Power(BeamWidthY_,2);
491 
492  } else {
493  edm::LogWarning("PrimaryVertexValidation")<<"No BeamSpot found!";
494  }
495 
496  h_BSx0->Fill(BSx0_);
497  h_BSy0->Fill(BSy0_);
498  h_BSz0->Fill(BSz0_);
499  h_Beamsigmaz->Fill(Beamsigmaz_);
500  h_BeamWidthX->Fill(BeamWidthX_);
501  h_BeamWidthY->Fill(BeamWidthY_);
502 
503  if(debug_)
504  edm::LogInfo("PrimaryVertexValidation")<<"Beamspot x:" <<BSx0_<<" y:"<<BSy0_<<" z:"<<BSz0_;
505 
506  //=======================================================
507  // Starts here ananlysis
508  //=======================================================
509 
510  RunNumber_=iEvent.eventAuxiliary().run();
511  h_runNumber->Fill(RunNumber_);
512 
513  if(h_runFromEvent->GetEntries()==0){
514  h_runFromEvent->SetBinContent(1,RunNumber_);
515  }
516 
517  LuminosityBlockNumber_=iEvent.eventAuxiliary().luminosityBlock();
518  EventNumber_=iEvent.eventAuxiliary().id().event();
519 
520  if(debug_)
521  edm::LogInfo("PrimaryVertexValidation")<<" looping over "<<trackCollectionHandle->size()<< "tracks";
522 
523  h_nTracks->Fill(trackCollectionHandle->size());
524 
525  //======================================================
526  // Interface RECO tracks to vertex reconstruction
527  //======================================================
528 
529  std::vector<TransientTrack> t_tks;
530  for (const auto & track : tracks){
531  TransientTrack tt = theB_->build(&(track));
532  tt.setBeamSpot(beamSpot);
533  t_tks.push_back(tt);
534 
535  }
536 
537  if(debug_) {
538  edm::LogInfo("PrimaryVertexValidation") << "Found: " << t_tks.size() << " reconstructed tracks";
539  }
540 
541  //======================================================
542  // select the tracks
543  //======================================================
544 
545  std::vector<TransientTrack> seltks = theTrackFilter_->select(t_tks);
546 
547  //======================================================
548  // clusterize tracks in Z
549  //======================================================
550 
551  vector< vector<TransientTrack> > clusters = theTrackClusterizer_->clusterize(seltks);
552 
553  if (debug_){
554  edm::LogInfo("PrimaryVertexValidation")<<" looping over: "<< clusters.size() << " clusters from " << t_tks.size() << " selected tracks";
555  }
556 
557  nClus_=clusters.size();
558  h_nClus->Fill(nClus_);
559 
560  //======================================================
561  // Starts loop on clusters
562  //======================================================
563  for (const auto & iclus : clusters){
564 
565  nTracksPerClus_=0;
566 
567  unsigned int i=0;
568  for(const auto & theTTrack : iclus)
569  {
570  i++;
571 
572  if ( nTracks_ >= nMaxtracks_ ) {
573  edm::LogError("PrimaryVertexValidation")<<" Warning - Number of tracks: " << nTracks_ << " , greater than " << nMaxtracks_;
574  continue;
575  }
576 
577  const Track & theTrack = theTTrack.track();
578 
579  pt_[nTracks_] = theTrack.pt();
580  p_[nTracks_] = theTrack.p();
581  nhits_[nTracks_] = theTrack.numberOfValidHits();
582  eta_[nTracks_] = theTrack.eta();
583  theta_[nTracks_] = theTrack.theta();
584  phi_[nTracks_] = theTrack.phi();
585  chi2_[nTracks_] = theTrack.chi2();
586  chi2ndof_[nTracks_] = theTrack.normalizedChi2();
587  charge_[nTracks_] = theTrack.charge();
588  qoverp_[nTracks_] = theTrack.qoverp();
589  dz_[nTracks_] = theTrack.dz();
590  dxy_[nTracks_] = theTrack.dxy();
591 
592  TrackBase::TrackQuality _trackQuality = TrackBase::qualityByName("highPurity");
593  isHighPurity_[nTracks_] = theTrack.quality(_trackQuality);
594 
595  math::XYZPoint point(BSx0_,BSy0_,BSz0_);
596  dxyBs_[nTracks_] = theTrack.dxy(point);
597  dzBs_[nTracks_] = theTrack.dz(point);
598 
599  xPCA_[nTracks_] = theTrack.vertex().x();
600  yPCA_[nTracks_] = theTrack.vertex().y();
601  zPCA_[nTracks_] = theTrack.vertex().z();
602 
603  //=======================================================
604  // Retrieve rechit information
605  //=======================================================
606 
607  const reco::HitPattern& hits = theTrack.hitPattern();
608 
609  int nRecHit1D=0;
610  int nRecHit2D=0;
611  int nhitinTIB = hits.numberOfValidStripTIBHits();
612  int nhitinTOB = hits.numberOfValidStripTOBHits();
613  int nhitinTID = hits.numberOfValidStripTIDHits();
614  int nhitinTEC = hits.numberOfValidStripTECHits();
615  int nhitinBPIX = hits.numberOfValidPixelBarrelHits();
616  int nhitinFPIX = hits.numberOfValidPixelEndcapHits();
617 
618  for (trackingRecHit_iterator iHit = theTTrack.recHitsBegin(); iHit != theTTrack.recHitsEnd(); ++iHit) {
619  if((*iHit)->isValid()) {
620 
621  if (this->isHit2D(**iHit)) {++nRecHit2D;}
622  else {++nRecHit1D; }
623  }
624  }
625 
626  nhits1D_[nTracks_] = nRecHit1D;
627  nhits2D_[nTracks_] = nRecHit2D;
628  nhitsBPIX_[nTracks_] = nhitinBPIX;
629  nhitsFPIX_[nTracks_] = nhitinFPIX;
630  nhitsTIB_[nTracks_] = nhitinTIB;
631  nhitsTID_[nTracks_] = nhitinTID;
632  nhitsTOB_[nTracks_] = nhitinTOB;
633  nhitsTEC_[nTracks_] = nhitinTEC;
634 
635  //=======================================================
636  // Good tracks for vertexing selection
637  //=======================================================
638 
639  bool pass = true;
640  if(askFirstLayerHit_) pass = this->hasFirstLayerPixelHits(theTTrack);
641  if (pass
642  && (theTrack.pt() >=ptOfProbe_)
643  && std::abs(theTrack.eta()) <= etaOfProbe_
644  && (theTrack.numberOfValidHits())>=nHitsOfProbe_
645  && (theTrack.p()) >= pOfProbe_ ){
646  isGoodTrack_[nTracks_]=1;
647  }
648 
649  //=======================================================
650  // Fit unbiased vertex
651  //=======================================================
652 
653  vector<TransientTrack> theFinalTracks;
654  theFinalTracks.clear();
655 
656  for (const auto & tk : iclus) {
657 
658  pass = this->hasFirstLayerPixelHits(tk);
659  if (pass){
660  if( tk == theTTrack ) continue;
661  else {
662  theFinalTracks.push_back(tk);
663  }
664  }
665  }
666 
667  if(theFinalTracks.size() > 1){
668 
669  if(debug_)
670  edm::LogInfo("PrimaryVertexValidation")<<"Transient Track Collection size: "<<theFinalTracks.size();
671  try{
672 
673  auto theFitter = std::unique_ptr<VertexFitter<5> >( new AdaptiveVertexFitter());
674  TransientVertex theFittedVertex = theFitter->vertex(theFinalTracks);
675 
676  //AdaptiveVertexFitter* theFitter = new AdaptiveVertexFitter;
677  //TransientVertex theFittedVertex = theFitter->vertex(theFinalTracks,beamSpot); // if you want the beam constraint
678 
679  double totalTrackWeights=0;
680  if(theFittedVertex.isValid ()){
681 
682 
683  if(theFittedVertex.hasTrackWeight()){
684  for(const auto & theFinalTrack : theFinalTracks){
685  sumOfWeightsUnbiasedVertex_[nTracks_] += theFittedVertex.trackWeight(theFinalTrack);
686  totalTrackWeights+= theFittedVertex.trackWeight(theFinalTrack);
687  h_fitVtxTrackWeights_->Fill(theFittedVertex.trackWeight(theFinalTrack));
688  }
689  }
690 
691  h_fitVtxTrackAverageWeight_->Fill(totalTrackWeights/theFinalTracks.size());
692 
693  const math::XYZPoint theRecoVertex(xOfflineVertex_,yOfflineVertex_,zOfflineVertex_);
694  const math::XYZPoint myVertex(theFittedVertex.position().x(),theFittedVertex.position().y(),theFittedVertex.position().z());
695 
696  const Vertex vertex = theFittedVertex;
697  fillTrackHistos(hDA,"all",&theTTrack,vertex,beamSpot,fBfield_);
698 
699  hasRecVertex_[nTracks_] = 1;
700  xUnbiasedVertex_[nTracks_] = theFittedVertex.position().x();
701  yUnbiasedVertex_[nTracks_] = theFittedVertex.position().y();
702  zUnbiasedVertex_[nTracks_] = theFittedVertex.position().z();
703 
704  chi2normUnbiasedVertex_[nTracks_] = theFittedVertex.normalisedChiSquared();
705  chi2UnbiasedVertex_[nTracks_] = theFittedVertex.totalChiSquared();
706  DOFUnbiasedVertex_[nTracks_] = theFittedVertex.degreesOfFreedom();
707  chi2ProbUnbiasedVertex_[nTracks_] = TMath::Prob(theFittedVertex.totalChiSquared(),(int)theFittedVertex.degreesOfFreedom());
708  tracksUsedForVertexing_[nTracks_] = theFinalTracks.size();
709 
710  h_fitVtxNtracks_->Fill(theFinalTracks.size());
711  h_fitVtxChi2_->Fill(theFittedVertex.totalChiSquared());
712  h_fitVtxNdof_->Fill(theFittedVertex.degreesOfFreedom());
713  h_fitVtxChi2ndf_->Fill(theFittedVertex.normalisedChiSquared());
714  h_fitVtxChi2Prob_->Fill(TMath::Prob(theFittedVertex.totalChiSquared(),(int)theFittedVertex.degreesOfFreedom()));
715 
716  // from my Vertex
717  double dxyFromMyVertex = theTrack.dxy(myVertex);
718  double dzFromMyVertex = theTrack.dz(myVertex);
719 
720  GlobalPoint vert(theFittedVertex.position().x(),theFittedVertex.position().y(),theFittedVertex.position().z());
721 
722  //FreeTrajectoryState theTrackNearVertex = theTTrack.trajectoryStateClosestToPoint(vert).theState();
723  //double dz_err = sqrt(theFittedVertex.positionError().czz() + theTrackNearVertex.cartesianError().position().czz());
724  //double dz_err = hypot(theTrack.dzError(),theFittedVertex.positionError().czz());
725 
726  double dz_err = sqrt(std::pow(theTrack.dzError(), 2) + theFittedVertex.positionError().czz());
727 
728 
729  // PV2D
730  std::pair<bool,Measurement1D> s_ip2dpv = signedTransverseImpactParameter(theTTrack,
731  GlobalVector(theTrack.px(),
732  theTrack.py(),
733  theTrack.pz()),
734  theFittedVertex);
735 
736  double s_ip2dpv_corr = s_ip2dpv.second.value();
737  double s_ip2dpv_err = s_ip2dpv.second.error();
738 
739  // PV3D
740  std::pair<bool, Measurement1D> s_ip3dpv = signedImpactParameter3D(theTTrack,
741  GlobalVector(theTrack.px(),
742  theTrack.py(),
743  theTrack.pz()),
744  theFittedVertex);
745 
746  double s_ip3dpv_corr = s_ip3dpv.second.value();
747  double s_ip3dpv_err = s_ip3dpv.second.error();
748 
749  // PV3D absolute
750  std::pair<bool,Measurement1D> ip3dpv = absoluteImpactParameter3D(theTTrack,theFittedVertex);
751  double ip3d_corr = ip3dpv.second.value();
752  double ip3d_err = ip3dpv.second.error();
753 
754  // with respect to any specified vertex, such as primary vertex
755  TrajectoryStateClosestToPoint traj = (theTTrack).trajectoryStateClosestToPoint(vert);
756 
757  GlobalPoint refPoint = traj.position();
758  GlobalPoint cPToVtx = traj.theState().position();
759 
760  float my_dx = refPoint.x() - myVertex.x();
761  float my_dy = refPoint.y() - myVertex.y();
762 
763  float my_dx2 = cPToVtx.x() - myVertex.x();
764  float my_dy2 = cPToVtx.y() - myVertex.y();
765 
766  float my_dxy = std::sqrt(my_dx*my_dx + my_dy*my_dy);
767 
768  double d0 = traj.perigeeParameters().transverseImpactParameter();
769  //double d0_error = traj.perigeeError().transverseImpactParameterError();
770  double z0 = traj.perigeeParameters().longitudinalImpactParameter();
771  double z0_error = traj.perigeeError().longitudinalImpactParameterError();
772 
773  if(debug_){
774  edm::LogInfo("PrimaryVertexValidation")<< "my_dx:" << my_dx
775  << " my_dy:" << my_dy
776  << " my_dxy:" << my_dxy
777  << " my_dx2:" << my_dx2
778  << " my_dy2:" << my_dy2
779  << " d0: " << d0
780  << " dxyFromVtx:" << dxyFromMyVertex << "\n"
781  << " ============================== "<< "\n"
782  << "diff1:" << std::abs(d0) - std::abs(my_dxy) << "\n"
783  << "diff2:" << std::abs(d0) - std::abs(dxyFromMyVertex) << "\n"
784  << "diff3:" << (my_dx - my_dx2) << " " << (my_dy - my_dy2) << "\n"
785  << std::endl;
786  }
787 
788  // define IPs
789 
790  dxyFromMyVertex_[nTracks_] = dxyFromMyVertex;
791  dxyErrorFromMyVertex_[nTracks_] = s_ip2dpv_err;
792  IPTsigFromMyVertex_[nTracks_] = dxyFromMyVertex/s_ip2dpv_err;
793 
794  dzFromMyVertex_[nTracks_] = dzFromMyVertex;
795  dzErrorFromMyVertex_[nTracks_] = dz_err;
796  IPLsigFromMyVertex_[nTracks_] = dzFromMyVertex/dz_err;
797 
798  d3DFromMyVertex_[nTracks_] = ip3d_corr;
799  d3DErrorFromMyVertex_[nTracks_] = ip3d_err;
800  IP3DsigFromMyVertex_[nTracks_] = (ip3d_corr/ip3d_err);
801 
802  // fill directly the histograms of residuals
803 
804  float trackphi = (theTrack.phi())*(180./M_PI);
805  float tracketa = theTrack.eta();
806  float trackpt = theTrack.pt();
807  float trackp = theTrack.p();
808  float tracknhits = theTrack.numberOfValidHits();
809 
810  // determine the module number and ladder
811 
812  int ladder_num = -1.;
813  int module_num = -1.;
814  int L1BPixHitCount = 0;
815 
816  for (trackingRecHit_iterator iHit = theTrack.recHitsBegin(); iHit != theTrack.recHitsEnd(); ++iHit) {
817  const DetId& detId = (*iHit)->geographicalId();
818  unsigned int subid = detId.subdetId();
819 
820  if((*iHit)->isValid() && ( subid == PixelSubdetector::PixelBarrel ) ) {
821  int layer = tTopo->pxbLayer(detId);
822  if(layer==1){
823  L1BPixHitCount+=1;
824  ladder_num = tTopo->pxbLadder(detId);
825  module_num = tTopo->pxbModule(detId);
826  }
827  }
828  }
829 
830  h_probeL1Ladder_->Fill(ladder_num);
831  h_probeL1Module_->Fill(module_num);
832  h_probeHasBPixL1Overlap_->Fill(L1BPixHitCount);
833 
834  // filling the pT-binned distributions
835 
836  for(int ipTBin=0; ipTBin<nPtBins_; ipTBin++){
837 
838  float pTF = mypT_bins_[ipTBin];
839  float pTL = mypT_bins_[ipTBin+1];
840 
841  if(debug_)
842  edm::LogInfo("PrimaryVertexValidation")<<"ipTBin:"<<ipTBin<< " "<<mypT_bins_[ipTBin]<< " < pT < "<<mypT_bins_[ipTBin+1]<<std::endl;
843 
844  if( std::abs(tracketa)<1.5 && (trackpt >= pTF && trackpt < pTL) ){
845 
846  if(debug_)
847  edm::LogInfo("PrimaryVertexValidation")<<"passes this cut: "<<mypT_bins_[ipTBin]<<std::endl;
848  PVValHelper::fillByIndex(h_dxy_pT_,ipTBin,dxyFromMyVertex*cmToum);
849  PVValHelper::fillByIndex(h_dz_pT_,ipTBin,dzFromMyVertex*cmToum);
850  PVValHelper::fillByIndex(h_norm_dxy_pT_,ipTBin,dxyFromMyVertex/s_ip2dpv_err);
851  PVValHelper::fillByIndex(h_norm_dz_pT_,ipTBin,dzFromMyVertex/dz_err);
852 
853  if(std::abs(tracketa)<1.){
854 
855  if(debug_)
856  edm::LogInfo("PrimaryVertexValidation")<<"passes tight eta cut: "<<mypT_bins_[ipTBin]<<std::endl;
857  PVValHelper::fillByIndex(h_dxy_Central_pT_,ipTBin,dxyFromMyVertex*cmToum);
858  PVValHelper::fillByIndex(h_dz_Central_pT_,ipTBin,dzFromMyVertex*cmToum);
859  PVValHelper::fillByIndex(h_norm_dxy_Central_pT_,ipTBin,dxyFromMyVertex/s_ip2dpv_err);
860  PVValHelper::fillByIndex(h_norm_dz_Central_pT_,ipTBin,dzFromMyVertex/dz_err);
861  }
862  }
863  }
864 
865  // checks on the probe track quality
866  if(trackpt >= ptOfProbe_
867  && std::abs(tracketa)<= etaOfProbe_
868  && tracknhits>=nHitsOfProbe_
869  && trackp >= pOfProbe_){
870 
871  std::pair<bool,bool> pixelOcc = pixelHitsCheck((theTTrack));
872 
873  if(debug_){
874  if(pixelOcc.first == true)
875  edm::LogInfo("PrimaryVertexValidation")<<"has BPIx hits"<<std::endl;
876  if(pixelOcc.second == true)
877  edm::LogInfo("PrimaryVertexValidation")<<"has FPix hits"<<std::endl;
878  }
879 
880  if(!doBPix_ && (pixelOcc.first == true)) continue;
881  if(!doFPix_ && (pixelOcc.second == true)) continue;
882 
883  fillTrackHistos(hDA,"sel",&(theTTrack),vertex,beamSpot,fBfield_);
884 
885  // probe checks
886  h_probePt_->Fill(theTrack.pt());
887  h_probePtRebin_->Fill(theTrack.pt());
888  h_probeP_->Fill(theTrack.p());
889  h_probeEta_->Fill(theTrack.eta());
890  h_probePhi_->Fill(theTrack.phi());
891  h2_probeEtaPhi_->Fill(theTrack.eta(),theTrack.phi());
892  h2_probeEtaPt_->Fill(theTrack.eta(),theTrack.pt());
893 
894  h_probeChi2_->Fill(theTrack.chi2());
895  h_probeNormChi2_->Fill(theTrack.normalizedChi2());
896  h_probeCharge_->Fill(theTrack.charge());
897  h_probeQoverP_->Fill(theTrack.qoverp());
898  h_probeHits_->Fill(theTrack.numberOfValidHits());
899  h_probeHits1D_->Fill(nRecHit1D);
900  h_probeHits2D_->Fill(nRecHit2D);
901  h_probeHitsInTIB_->Fill(nhitinBPIX);
902  h_probeHitsInTOB_->Fill(nhitinFPIX);
903  h_probeHitsInTID_->Fill(nhitinTIB);
904  h_probeHitsInTEC_->Fill(nhitinTID);
905  h_probeHitsInBPIX_->Fill(nhitinTOB);
906  h_probeHitsInFPIX_->Fill(nhitinTEC);
907 
908  float dxyRecoV = theTrack.dz(theRecoVertex);
909  float dzRecoV = theTrack.dxy(theRecoVertex);
910  float dxysigmaRecoV = TMath::Sqrt(theTrack.d0Error()*theTrack.d0Error()+xErrOfflineVertex_*yErrOfflineVertex_);
911  float dzsigmaRecoV = TMath::Sqrt(theTrack.dzError()*theTrack.dzError()+zErrOfflineVertex_*zErrOfflineVertex_);
912 
913  double zTrack=(theTTrack.stateAtBeamLine().trackStateAtPCA()).position().z();
914  double zVertex=theFittedVertex.position().z();
915  double tantheta=tan((theTTrack.stateAtBeamLine().trackStateAtPCA()).momentum().theta());
916 
917  double dz2= pow(theTrack.dzError(),2)+wxy2_/pow(tantheta,2);
918  double restrkz = zTrack-zVertex;
919  double pulltrkz = (zTrack-zVertex)/TMath::Sqrt(dz2);
920 
921  h_probedxyRecoV_->Fill(dxyRecoV);
922  h_probedzRecoV_->Fill(dzRecoV);
923 
924  h_probedzRefitV_->Fill(dxyFromMyVertex);
925  h_probedxyRefitV_->Fill(dzFromMyVertex);
926 
927  h_probed0RefitV_->Fill(d0);
928  h_probez0RefitV_->Fill(z0);
929 
930  h_probesignIP2DRefitV_->Fill(s_ip2dpv_corr);
931  h_probed3DRefitV_->Fill(ip3d_corr);
932  h_probereszRefitV_->Fill(restrkz);
933 
934  h_probeRecoVSigZ_->Fill(dzRecoV/dzsigmaRecoV);
935  h_probeRecoVSigXY_->Fill(dxyRecoV/dxysigmaRecoV);
936  h_probeRefitVSigZ_->Fill(dzFromMyVertex/dz_err);
937  h_probeRefitVSigXY_->Fill(dxyFromMyVertex/s_ip2dpv_err);
938  h_probeRefitVSig3D_->Fill(ip3d_corr/ip3d_err);
939  h_probeRefitVLogSig3D_->Fill(log10(ip3d_corr/ip3d_err));
940  h_probeRefitVSigResZ_->Fill(pulltrkz);
941 
942  a_dxyVsPhi->Fill(trackphi,dxyFromMyVertex*cmToum);
943  a_dzVsPhi->Fill(trackphi,z0*cmToum);
944  n_dxyVsPhi->Fill(trackphi,dxyFromMyVertex/s_ip2dpv_err);
945  n_dzVsPhi->Fill(trackphi,z0/z0_error);
946 
947  a_dxyVsEta->Fill(tracketa,dxyFromMyVertex*cmToum);
948  a_dzVsEta->Fill(tracketa,z0*cmToum);
949  n_dxyVsEta->Fill(tracketa,dxyFromMyVertex/s_ip2dpv_err);
950  n_dzVsEta->Fill(tracketa,z0/z0_error);
951 
952  if( ladder_num > 0 && module_num > 0 ) {
953 
954  LogDebug("PrimaryVertexValidation")<<" ladder_num"<<ladder_num <<" module_num"<<module_num <<std::endl;
955 
956  PVValHelper::fillByIndex(h_dxy_modZ_,module_num-1,dxyFromMyVertex*cmToum);
957  PVValHelper::fillByIndex(h_dz_modZ_,module_num-1,dzFromMyVertex*cmToum);
958  PVValHelper::fillByIndex(h_norm_dxy_modZ_,module_num-1,dxyFromMyVertex/s_ip2dpv_err);
959  PVValHelper::fillByIndex(h_norm_dz_modZ_,module_num-1,dzFromMyVertex/dz_err);
960 
961  PVValHelper::fillByIndex(h_dxy_ladder_,ladder_num-1,dxyFromMyVertex*cmToum);
962 
963  LogDebug("PrimaryVertexValidation")<<"h_dxy_ladder size:" <<h_dxy_ladder_.size() << std::endl;
964 
965  if(L1BPixHitCount==1){
966  PVValHelper::fillByIndex(h_dxy_ladderNoOverlap_,ladder_num-1,dxyFromMyVertex*cmToum);
967  } else {
968  PVValHelper::fillByIndex(h_dxy_ladderOverlap_,ladder_num-1,dxyFromMyVertex*cmToum);
969  }
970 
971  PVValHelper::fillByIndex(h_dz_ladder_,ladder_num-1,dzFromMyVertex*cmToum);
972  PVValHelper::fillByIndex(h_norm_dxy_ladder_,ladder_num-1,dxyFromMyVertex/s_ip2dpv_err);
973  PVValHelper::fillByIndex(h_norm_dz_ladder_,ladder_num-1,dzFromMyVertex/dz_err);
974 
975  }
976 
977  // filling the binned distributions
978  for(int i=0; i<nBins_; i++){
979 
980  float phiF = theDetails_.trendbins[PVValHelper::phi][i];
981  float phiL = theDetails_.trendbins[PVValHelper::phi][i+1];
982 
983  float etaF = theDetails_.trendbins[PVValHelper::eta][i];
984  float etaL = theDetails_.trendbins[PVValHelper::eta][i+1];
985 
986  if(tracketa >= etaF && tracketa < etaL ){
987 
988  PVValHelper::fillByIndex(a_dxyEtaResiduals,i,dxyFromMyVertex*cmToum,"1");
989  PVValHelper::fillByIndex(a_dxEtaResiduals,i,my_dx*cmToum,"2");
990  PVValHelper::fillByIndex(a_dyEtaResiduals,i,my_dy*cmToum,"3");
991  PVValHelper::fillByIndex(a_dzEtaResiduals,i,dzFromMyVertex*cmToum,"4");
992  PVValHelper::fillByIndex(n_dxyEtaResiduals,i,dxyFromMyVertex/s_ip2dpv_err,"5");
993  PVValHelper::fillByIndex(n_dzEtaResiduals,i,dzFromMyVertex/dz_err,"6");
994  PVValHelper::fillByIndex(a_IP2DEtaResiduals,i,s_ip2dpv_corr*cmToum,"7");
995  PVValHelper::fillByIndex(n_IP2DEtaResiduals,i,s_ip2dpv_corr/s_ip2dpv_err,"8");
996  PVValHelper::fillByIndex(a_reszEtaResiduals,i,restrkz*cmToum,"9");
997  PVValHelper::fillByIndex(n_reszEtaResiduals,i,pulltrkz,"10");
998  PVValHelper::fillByIndex(a_d3DEtaResiduals,i,ip3d_corr*cmToum,"11");
999  PVValHelper::fillByIndex(n_d3DEtaResiduals,i,ip3d_corr/ip3d_err,"12");
1000  PVValHelper::fillByIndex(a_IP3DEtaResiduals,i,s_ip3dpv_corr*cmToum,"13");
1001  PVValHelper::fillByIndex(n_IP3DEtaResiduals,i,s_ip3dpv_corr/s_ip3dpv_err,"14");
1002 
1003  }
1004 
1005  if(trackphi >= phiF && trackphi < phiL ){
1006 
1007  PVValHelper::fillByIndex(a_dxyPhiResiduals,i,dxyFromMyVertex*cmToum,"15");
1008  PVValHelper::fillByIndex(a_dxPhiResiduals,i,my_dx*cmToum,"16");
1009  PVValHelper::fillByIndex(a_dyPhiResiduals,i,my_dy*cmToum,"17");
1010  PVValHelper::fillByIndex(a_dzPhiResiduals,i,dzFromMyVertex*cmToum,"18");
1011  PVValHelper::fillByIndex(n_dxyPhiResiduals,i,dxyFromMyVertex/s_ip2dpv_err,"19");
1012  PVValHelper::fillByIndex(n_dzPhiResiduals,i,dzFromMyVertex/dz_err,"20");
1013  PVValHelper::fillByIndex(a_IP2DPhiResiduals,i,s_ip2dpv_corr*cmToum,"21");
1014  PVValHelper::fillByIndex(n_IP2DPhiResiduals,i,s_ip2dpv_corr/s_ip2dpv_err,"22");
1015  PVValHelper::fillByIndex(a_reszPhiResiduals,i,restrkz*cmToum,"23");
1016  PVValHelper::fillByIndex(n_reszPhiResiduals,i,pulltrkz,"24");
1017  PVValHelper::fillByIndex(a_d3DPhiResiduals,i,ip3d_corr*cmToum,"25");
1018  PVValHelper::fillByIndex(n_d3DPhiResiduals,i,ip3d_corr/ip3d_err,"26");
1019  PVValHelper::fillByIndex(a_IP3DPhiResiduals,i,s_ip3dpv_corr*cmToum,"27");
1020  PVValHelper::fillByIndex(n_IP3DPhiResiduals,i,s_ip3dpv_corr/s_ip3dpv_err,"28");
1021 
1022  for(int j=0; j<nBins_; j++){
1023 
1024  float etaJ = theDetails_.trendbins[PVValHelper::eta][j];
1025  float etaK = theDetails_.trendbins[PVValHelper::eta][j+1];
1026 
1027  if(tracketa >= etaJ && tracketa < etaK ){
1028  a_dxyResidualsMap[i][j]->Fill(dxyFromMyVertex*cmToum);
1029  a_dzResidualsMap[i][j]->Fill(dzFromMyVertex*cmToum);
1030  n_dxyResidualsMap[i][j]->Fill(dxyFromMyVertex/s_ip2dpv_err);
1031  n_dzResidualsMap[i][j]->Fill(dzFromMyVertex/dz_err);
1032  a_d3DResidualsMap[i][j]->Fill(ip3d_corr*cmToum);
1033  n_d3DResidualsMap[i][j]->Fill(ip3d_corr/ip3d_err);
1034 
1035  }
1036  }
1037  }
1038  }
1039  }
1040 
1041  if(debug_){
1042  edm::LogInfo("PrimaryVertexValidation")<<" myVertex.x()= "<<myVertex.x()<<"\n"
1043  <<" myVertex.y()= "<<myVertex.y()<<" \n"
1044  <<" myVertex.z()= "<<myVertex.z()<<" \n"
1045  <<" theTrack.dz(myVertex)= "<<theTrack.dz(myVertex)<<" \n"
1046  <<" zPCA -myVertex.z() = "<<(theTrack.vertex().z() -myVertex.z());
1047 
1048  }// ends if debug_
1049  } // ends if the fitted vertex is Valid
1050 
1051  //delete theFitter;
1052 
1053  } catch ( cms::Exception& er ) {
1054  LogTrace("PrimaryVertexValidation")<<"caught std::exception "<<er.what()<<std::endl;
1055  }
1056 
1057  } //ends if theFinalTracks.size() > 2
1058 
1059  else {
1060  if(debug_)
1061  edm::LogInfo("PrimaryVertexValidation")<<"Not enough tracks to make a vertex. Returns no vertex info";
1062  }
1063 
1064  ++nTracks_;
1065  ++nTracksPerClus_;
1066 
1067  if(debug_)
1068  edm::LogInfo("PrimaryVertexValidation")<<"Track "<<i<<" : pT = "<<theTrack.pt();
1069 
1070  }// for loop on tracks
1071  } // for loop on track clusters
1072 
1073  // Fill the TTree if needed
1074 
1075  if(storeNtuple_){
1076  rootTree_->Fill();
1077  }
1078 
1079 
1080 }
1081 
1082 // ------------ method called to discriminate 1D from 2D hits ------------
1083 bool PrimaryVertexValidation::isHit2D(const TrackingRecHit &hit) const
1084 {
1085  if (hit.dimension() < 2) {
1086  return false; // some (muon...) stuff really has RecHit1D
1087  } else {
1088  const DetId detId(hit.geographicalId());
1089  if (detId.det() == DetId::Tracker) {
1090  if (detId.subdetId() == PixelSubdetector::PixelBarrel || detId.subdetId() == PixelSubdetector::PixelEndcap) {
1091  return true; // pixel is always 2D
1092  } else { // should be SiStrip now
1093  if (dynamic_cast<const SiStripRecHit2D*>(&hit)) return false; // normal hit
1094  else if (dynamic_cast<const SiStripMatchedRecHit2D*>(&hit)) return true; // matched is 2D
1095  else if (dynamic_cast<const ProjectedSiStripRecHit2D*>(&hit)) return false; // crazy hit...
1096  else {
1097  edm::LogError("UnkownType") << "@SUB=AlignmentTrackSelector::isHit2D"
1098  << "Tracker hit not in pixel and neither SiStripRecHit2D nor "
1099  << "SiStripMatchedRecHit2D nor ProjectedSiStripRecHit2D.";
1100  return false;
1101  }
1102  }
1103  } else { // not tracker??
1104  edm::LogWarning("DetectorMismatch") << "@SUB=AlignmentTrackSelector::isHit2D"
1105  << "Hit not in tracker with 'official' dimension >=2.";
1106  return true; // dimension() >= 2 so accept that...
1107  }
1108  }
1109  // never reached...
1110 }
1111 
1112 
1113 // ------------ method to check the presence of pixel hits ------------
1114 std::pair<bool,bool> PrimaryVertexValidation::pixelHitsCheck(const reco::TransientTrack& track){
1115 
1116  bool hasBPixHits = false;
1117  bool hasFPixHits = false;
1118 
1119  const reco::HitPattern& p = track.hitPattern();
1120  if(p.numberOfValidPixelEndcapHits()!=0){
1121  hasFPixHits = true;
1122  }
1123  if(p.numberOfValidPixelBarrelHits()!=0){
1124  hasBPixHits = true;
1125  }
1126 
1127  return std::make_pair(hasBPixHits,hasFPixHits);
1128 }
1129 
1130 
1131 // ------------ method to check the presence of pixel hits ------------
1132 bool PrimaryVertexValidation::hasFirstLayerPixelHits(const reco::TransientTrack& track)
1133 {
1134  using namespace reco;
1135  const HitPattern& p = track.hitPattern();
1136  for (int i=0; i<p.numberOfAllHits(HitPattern::TRACK_HITS); i++) {
1137  uint32_t pattern = p.getHitPattern(HitPattern::TRACK_HITS, i);
1138  if (p.pixelBarrelHitFilter(pattern) || p.pixelEndcapHitFilter(pattern) ) {
1139  if (p.getLayer(pattern) == 1) {
1140  if (p.validHitFilter(pattern)) {
1141  return true;
1142  }
1143  }
1144  }
1145  }
1146  return false;
1147 }
1148 
1149 // ------------ method called once each job before begining the event loop ------------
1150 void PrimaryVertexValidation::beginJob()
1151 {
1152  edm::LogInfo("PrimaryVertexValidation")
1153  <<"######################################\n"
1154  <<"Begin Job \n"
1155  <<"######################################";
1156 
1157  // Define TTree for output
1158  Nevt_ = 0;
1159 
1160  // rootFile_ = new TFile(filename_.c_str(),"recreate");
1161  rootTree_ = fs->make<TTree>("tree","PV Validation tree");
1162 
1163  // Track Paramters
1164 
1165  if(lightNtupleSwitch_){
1166 
1167  rootTree_->Branch("EventNumber",&EventNumber_,"EventNumber/i");
1168  rootTree_->Branch("RunNumber",&RunNumber_,"RunNumber/i");
1169  rootTree_->Branch("LuminosityBlockNumber",&LuminosityBlockNumber_,"LuminosityBlockNumber/i");
1170  rootTree_->Branch("nOfflineVertices",&nOfflineVertices_,"nOfflineVertices/I");
1171  rootTree_->Branch("nTracks",&nTracks_,"nTracks/I");
1172  rootTree_->Branch("phi",&phi_,"phi[nTracks]/D");
1173  rootTree_->Branch("eta",&eta_,"eta[nTracks]/D");
1174  rootTree_->Branch("pt",&pt_,"pt[nTracks]/D");
1175  rootTree_->Branch("dxyFromMyVertex",&dxyFromMyVertex_,"dxyFromMyVertex[nTracks]/D");
1176  rootTree_->Branch("dzFromMyVertex",&dzFromMyVertex_,"dzFromMyVertex[nTracks]/D");
1177  rootTree_->Branch("d3DFromMyVertex",&d3DFromMyVertex_,"d3DFromMyVertex[nTracks]/D");
1178  rootTree_->Branch("IPTsigFromMyVertex",&IPTsigFromMyVertex_,"IPTsigFromMyVertex_[nTracks]/D");
1179  rootTree_->Branch("IPLsigFromMyVertex",&IPLsigFromMyVertex_,"IPLsigFromMyVertex_[nTracks]/D");
1180  rootTree_->Branch("IP3DsigFromMyVertex",&IP3DsigFromMyVertex_,"IP3DsigFromMyVertex_[nTracks]/D");
1181  rootTree_->Branch("hasRecVertex",&hasRecVertex_,"hasRecVertex[nTracks]/I");
1182  rootTree_->Branch("isGoodTrack",&isGoodTrack_,"isGoodTrack[nTracks]/I");
1183  rootTree_->Branch("isHighPurity",&isHighPurity_,"isHighPurity_[nTracks]/I");
1184 
1185  } else {
1186 
1187  rootTree_->Branch("nTracks",&nTracks_,"nTracks/I");
1188  rootTree_->Branch("nTracksPerClus",&nTracksPerClus_,"nTracksPerClus/I");
1189  rootTree_->Branch("nClus",&nClus_,"nClus/I");
1190  rootTree_->Branch("xOfflineVertex",&xOfflineVertex_,"xOfflineVertex/D");
1191  rootTree_->Branch("yOfflineVertex",&yOfflineVertex_,"yOfflineVertex/D");
1192  rootTree_->Branch("zOfflineVertex",&zOfflineVertex_,"zOfflineVertex/D");
1193  rootTree_->Branch("BSx0",&BSx0_,"BSx0/D");
1194  rootTree_->Branch("BSy0",&BSy0_,"BSy0/D");
1195  rootTree_->Branch("BSz0",&BSz0_,"BSz0/D");
1196  rootTree_->Branch("Beamsigmaz",&Beamsigmaz_,"Beamsigmaz/D");
1197  rootTree_->Branch("Beamdxdz",&Beamdxdz_,"Beamdxdz/D");
1198  rootTree_->Branch("BeamWidthX",&BeamWidthX_,"BeamWidthX/D");
1199  rootTree_->Branch("BeamWidthY",&BeamWidthY_,"BeamWidthY/D");
1200  rootTree_->Branch("pt",&pt_,"pt[nTracks]/D");
1201  rootTree_->Branch("p",&p_,"p[nTracks]/D");
1202  rootTree_->Branch("nhits",&nhits_,"nhits[nTracks]/I");
1203  rootTree_->Branch("nhits1D",&nhits1D_,"nhits1D[nTracks]/I");
1204  rootTree_->Branch("nhits2D",&nhits2D_,"nhits2D[nTracks]/I");
1205  rootTree_->Branch("nhitsBPIX",&nhitsBPIX_,"nhitsBPIX[nTracks]/I");
1206  rootTree_->Branch("nhitsFPIX",&nhitsFPIX_,"nhitsFPIX[nTracks]/I");
1207  rootTree_->Branch("nhitsTIB",&nhitsTIB_,"nhitsTIB[nTracks]/I");
1208  rootTree_->Branch("nhitsTID",&nhitsTID_,"nhitsTID[nTracks]/I");
1209  rootTree_->Branch("nhitsTOB",&nhitsTOB_,"nhitsTOB[nTracks]/I");
1210  rootTree_->Branch("nhitsTEC",&nhitsTEC_,"nhitsTEC[nTracks]/I");
1211  rootTree_->Branch("eta",&eta_,"eta[nTracks]/D");
1212  rootTree_->Branch("theta",&theta_,"theta[nTracks]/D");
1213  rootTree_->Branch("phi",&phi_,"phi[nTracks]/D");
1214  rootTree_->Branch("chi2",&chi2_,"chi2[nTracks]/D");
1215  rootTree_->Branch("chi2ndof",&chi2ndof_,"chi2ndof[nTracks]/D");
1216  rootTree_->Branch("charge",&charge_,"charge[nTracks]/I");
1217  rootTree_->Branch("qoverp",&qoverp_,"qoverp[nTracks]/D");
1218  rootTree_->Branch("dz",&dz_,"dz[nTracks]/D");
1219  rootTree_->Branch("dxy",&dxy_,"dxy[nTracks]/D");
1220  rootTree_->Branch("dzBs",&dzBs_,"dzBs[nTracks]/D");
1221  rootTree_->Branch("dxyBs",&dxyBs_,"dxyBs[nTracks]/D");
1222  rootTree_->Branch("xPCA",&xPCA_,"xPCA[nTracks]/D");
1223  rootTree_->Branch("yPCA",&yPCA_,"yPCA[nTracks]/D");
1224  rootTree_->Branch("zPCA",&zPCA_,"zPCA[nTracks]/D");
1225  rootTree_->Branch("xUnbiasedVertex",&xUnbiasedVertex_,"xUnbiasedVertex[nTracks]/D");
1226  rootTree_->Branch("yUnbiasedVertex",&yUnbiasedVertex_,"yUnbiasedVertex[nTracks]/D");
1227  rootTree_->Branch("zUnbiasedVertex",&zUnbiasedVertex_,"zUnbiasedVertex[nTracks]/D");
1228  rootTree_->Branch("chi2normUnbiasedVertex",&chi2normUnbiasedVertex_,"chi2normUnbiasedVertex[nTracks]/F");
1229  rootTree_->Branch("chi2UnbiasedVertex",&chi2UnbiasedVertex_,"chi2UnbiasedVertex[nTracks]/F");
1230  rootTree_->Branch("DOFUnbiasedVertex",&DOFUnbiasedVertex_," DOFUnbiasedVertex[nTracks]/F");
1231  rootTree_->Branch("chi2ProbUnbiasedVertex",&chi2ProbUnbiasedVertex_,"chi2ProbUnbiasedVertex[nTracks]/F");
1232  rootTree_->Branch("sumOfWeightsUnbiasedVertex",&sumOfWeightsUnbiasedVertex_,"sumOfWeightsUnbiasedVertex[nTracks]/F");
1233  rootTree_->Branch("tracksUsedForVertexing",&tracksUsedForVertexing_,"tracksUsedForVertexing[nTracks]/I");
1234  rootTree_->Branch("dxyFromMyVertex",&dxyFromMyVertex_,"dxyFromMyVertex[nTracks]/D");
1235  rootTree_->Branch("dzFromMyVertex",&dzFromMyVertex_,"dzFromMyVertex[nTracks]/D");
1236  rootTree_->Branch("dxyErrorFromMyVertex",&dxyErrorFromMyVertex_,"dxyErrorFromMyVertex_[nTracks]/D");
1237  rootTree_->Branch("dzErrorFromMyVertex",&dzErrorFromMyVertex_,"dzErrorFromMyVertex_[nTracks]/D");
1238  rootTree_->Branch("IPTsigFromMyVertex",&IPTsigFromMyVertex_,"IPTsigFromMyVertex_[nTracks]/D");
1239  rootTree_->Branch("IPLsigFromMyVertex",&IPLsigFromMyVertex_,"IPLsigFromMyVertex_[nTracks]/D");
1240  rootTree_->Branch("hasRecVertex",&hasRecVertex_,"hasRecVertex[nTracks]/I");
1241  rootTree_->Branch("isGoodTrack",&isGoodTrack_,"isGoodTrack[nTracks]/I");
1242  }
1243 
1244  // event histograms
1245  TFileDirectory EventFeatures = fs->mkdir("EventFeatures");
1246 
1247  TH1F::SetDefaultSumw2(kTRUE);
1248 
1249  h_lumiFromConfig = EventFeatures.make<TH1F>("h_lumiFromConfig","luminosity from config;;luminosity of present run",1,-0.5,0.5);
1250  h_lumiFromConfig->SetBinContent(1,intLumi_);
1251 
1252  h_runFromConfig = EventFeatures.make<TH1I>("h_runFromConfig","run number from config;;run number (from configuration)",
1253  runControlNumbers_.size(),0.,runControlNumbers_.size());
1254  for(const auto r : runControlNumbers_){
1255  h_runFromConfig->SetBinContent(r+1, r);
1256  }
1257 
1258  h_runFromEvent = EventFeatures.make<TH1I>("h_runFromEvent","run number from config;;run number (from event)",1,-0.5,0.5);
1259  h_nTracks = EventFeatures.make<TH1F>("h_nTracks","number of tracks per event;n_{tracks}/event;n_{events}",300,-0.5,299.5);
1260  h_nClus = EventFeatures.make<TH1F>("h_nClus","number of track clusters;n_{clusters}/event;n_{events}",50,-0.5,49.5);
1261  h_nOfflineVertices = EventFeatures.make<TH1F>("h_nOfflineVertices","number of offline reconstructed vertices;n_{vertices}/event;n_{events}",50,-0.5,49.5);
1262  h_runNumber = EventFeatures.make<TH1F>("h_runNumber","run number;run number;n_{events}",100000,250000.,350000.);
1263  h_xOfflineVertex = EventFeatures.make<TH1F>("h_xOfflineVertex","x-coordinate of offline vertex;x_{vertex};n_{events}",100,-0.1,0.1);
1264  h_yOfflineVertex = EventFeatures.make<TH1F>("h_yOfflineVertex","y-coordinate of offline vertex;y_{vertex};n_{events}",100,-0.1,0.1);
1265  h_zOfflineVertex = EventFeatures.make<TH1F>("h_zOfflineVertex","z-coordinate of offline vertex;z_{vertex};n_{events}",100,-30.,30.);
1266  h_xErrOfflineVertex = EventFeatures.make<TH1F>("h_xErrOfflineVertex","x-coordinate error of offline vertex;err_{x}^{vtx};n_{events}",100,0.,0.01);
1267  h_yErrOfflineVertex = EventFeatures.make<TH1F>("h_yErrOfflineVertex","y-coordinate error of offline vertex;err_{y}^{vtx};n_{events}",100,0.,0.01);
1268  h_zErrOfflineVertex = EventFeatures.make<TH1F>("h_zErrOfflineVertex","z-coordinate error of offline vertex;err_{z}^{vtx};n_{events}",100,0.,10.);
1269  h_BSx0 = EventFeatures.make<TH1F>("h_BSx0","x-coordinate of reco beamspot;x^{BS}_{0};n_{events}",100,-0.1,0.1);
1270  h_BSy0 = EventFeatures.make<TH1F>("h_BSy0","y-coordinate of reco beamspot;y^{BS}_{0};n_{events}",100,-0.1,0.1);
1271  h_BSz0 = EventFeatures.make<TH1F>("h_BSz0","z-coordinate of reco beamspot;z^{BS}_{0};n_{events}",100,-1.,1.);
1272  h_Beamsigmaz = EventFeatures.make<TH1F>("h_Beamsigmaz","z-coordinate beam width;#sigma_{Z}^{beam};n_{events}",100,0.,1.);
1273  h_BeamWidthX = EventFeatures.make<TH1F>("h_BeamWidthX","x-coordinate beam width;#sigma_{X}^{beam};n_{events}",100,0.,0.01);
1274  h_BeamWidthY = EventFeatures.make<TH1F>("h_BeamWidthY","y-coordinate beam width;#sigma_{Y}^{beam};n_{events}",100,0.,0.01);
1275 
1276  h_etaMax = EventFeatures.make<TH1F>("etaMax","etaMax",1,-0.5,0.5);
1277  h_pTinfo = EventFeatures.make<TH1F>("pTinfo","pTinfo",3,-1.5,1.5);
1278  h_pTinfo->GetXaxis()->SetBinLabel(1,"n. bins");
1279  h_pTinfo->GetXaxis()->SetBinLabel(2,"pT min");
1280  h_pTinfo->GetXaxis()->SetBinLabel(3,"pT max");
1281 
1282  h_nbins = EventFeatures.make<TH1F>("nbins","nbins",1,-0.5,0.5);
1283  h_nLadders = EventFeatures.make<TH1F>("nladders","n. ladders",1,-0.5,0.5);
1284 
1285  // probe track histograms
1286  TFileDirectory ProbeFeatures = fs->mkdir("ProbeTrackFeatures");
1287 
1288  h_probePt_ = ProbeFeatures.make<TH1F>("h_probePt","p_{T} of probe track;track p_{T} (GeV); tracks",100,0.,50.);
1289  h_probePtRebin_ = ProbeFeatures.make<TH1F>("h_probePtRebin","p_{T} of probe track;track p_{T} (GeV); tracks",mypT_bins_.size()-1,mypT_bins_.data());
1290  h_probeP_ = ProbeFeatures.make<TH1F>("h_probeP","momentum of probe track;track p (GeV); tracks",100,0.,100.);
1291  h_probeEta_ = ProbeFeatures.make<TH1F>("h_probeEta","#eta of the probe track;track #eta;tracks",54,-2.8,2.8);
1292  h_probePhi_ = ProbeFeatures.make<TH1F>("h_probePhi","#phi of probe track;track #phi (rad);tracks",100,-3.15,3.15);
1293 
1294  h2_probeEtaPhi_ = ProbeFeatures.make<TH2F>("h2_probeEtaPhi","probe track #phi vs #eta;#eta of probe track;track #phi of probe track (rad); tracks",54,-2.8,2.8,100,-3.15,3.15);
1295  h2_probeEtaPt_ = ProbeFeatures.make<TH2F>("h2_probeEtaPt","probe track p_{T} vs #eta;#eta of probe track;track p_{T} (GeV); tracks",54,-2.8,2.8,100,0.,50.);
1296 
1297  h_probeChi2_ = ProbeFeatures.make<TH1F>("h_probeChi2","#chi^{2} of probe track;track #chi^{2}; tracks",100,0.,100.);
1298  h_probeNormChi2_ = ProbeFeatures.make<TH1F>("h_probeNormChi2"," normalized #chi^{2} of probe track;track #chi^{2}/ndof; tracks",100,0.,10.);
1299  h_probeCharge_ = ProbeFeatures.make<TH1F>("h_probeCharge","charge of profe track;track charge Q;tracks",3,-1.5,1.5);
1300  h_probeQoverP_ = ProbeFeatures.make<TH1F>("h_probeQoverP","q/p of probe track; track Q/p (GeV^{-1});tracks",200,-1.,1.);
1301  h_probedzRecoV_ = ProbeFeatures.make<TH1F>("h_probedzRecoV","d_{z}(V_{offline}) of probe track;track d_{z}(V_{off}) (cm);tracks",200,-1.,1.);
1302  h_probedxyRecoV_ = ProbeFeatures.make<TH1F>("h_probedxyRecoV","d_{xy}(V_{offline}) of probe track;track d_{xy}(V_{off}) (cm);tracks",200,-1.,1.);
1303  h_probedzRefitV_ = ProbeFeatures.make<TH1F>("h_probedzRefitV","d_{z}(V_{refit}) of probe track;track d_{z}(V_{fit}) (cm);tracks",200,-0.5,0.5);
1304  h_probesignIP2DRefitV_ = ProbeFeatures.make<TH1F>("h_probesignIPRefitV","ip_{2D}(V_{refit}) of probe track;track ip_{2D}(V_{fit}) (cm);tracks",200,-1.,1.);
1305  h_probedxyRefitV_ = ProbeFeatures.make<TH1F>("h_probedxyRefitV","d_{xy}(V_{refit}) of probe track;track d_{xy}(V_{fit}) (cm);tracks",200,-0.5,0.5);
1306 
1307  h_probez0RefitV_ = ProbeFeatures.make<TH1F>("h_probez0RefitV","z_{0}(V_{refit}) of probe track;track z_{0}(V_{fit}) (cm);tracks",200,-1.,1.);
1308  h_probed0RefitV_ = ProbeFeatures.make<TH1F>("h_probed0RefitV","d_{0}(V_{refit}) of probe track;track d_{0}(V_{fit}) (cm);tracks",200,-1.,1.);
1309 
1310  h_probed3DRefitV_ = ProbeFeatures.make<TH1F>("h_probed3DRefitV","d_{3D}(V_{refit}) of probe track;track d_{3D}(V_{fit}) (cm);tracks",200,0.,1.);
1311  h_probereszRefitV_ = ProbeFeatures.make<TH1F>("h_probeReszRefitV","z_{track} -z_{V_{refit}};track res_{z}(V_{refit}) (cm);tracks",200,-1.,1.);
1312 
1313  h_probeRecoVSigZ_ = ProbeFeatures.make<TH1F>("h_probeRecoVSigZ" ,"Longitudinal DCA Significance (reco);d_{z}(V_{off})/#sigma_{dz};tracks",100,-8,8);
1314  h_probeRecoVSigXY_ = ProbeFeatures.make<TH1F>("h_probeRecoVSigXY" ,"Transverse DCA Significance (reco);d_{xy}(V_{off})/#sigma_{dxy};tracks",100,-8,8);
1315  h_probeRefitVSigZ_ = ProbeFeatures.make<TH1F>("h_probeRefitVSigZ" ,"Longitudinal DCA Significance (refit);d_{z}(V_{fit})/#sigma_{dz};tracks",100,-8,8);
1316  h_probeRefitVSigXY_= ProbeFeatures.make<TH1F>("h_probeRefitVSigXY","Transverse DCA Significance (refit);d_{xy}(V_{fit})/#sigma_{dxy};tracks",100,-8,8);
1317  h_probeRefitVSig3D_= ProbeFeatures.make<TH1F>("h_probeRefitVSig3D","3D DCA Significance (refit);d_{3D}/#sigma_{3D};tracks",100,0.,20.);
1318  h_probeRefitVLogSig3D_ = ProbeFeatures.make<TH1F>("h_probeRefitVLogSig3D","log_{10}(3D DCA-Significance) (refit);log_{10}(d_{3D}/#sigma_{3D});tracks",100,-5.,4.);
1319  h_probeRefitVSigResZ_ = ProbeFeatures.make<TH1F>("h_probeRefitVSigResZ" ,"Longitudinal residual significance (refit);(z_{track} -z_{V_{fit}})/#sigma_{res_{z}};tracks",100,-8,8);
1320 
1321  h_probeHits_ = ProbeFeatures.make<TH1F>("h_probeNRechits" ,"N_{hits} ;N_{hits} ;tracks",40,-0.5,39.5);
1322  h_probeHits1D_ = ProbeFeatures.make<TH1F>("h_probeNRechits1D" ,"N_{hits} 1D ;N_{hits} 1D ;tracks",40,-0.5,39.5);
1323  h_probeHits2D_ = ProbeFeatures.make<TH1F>("h_probeNRechits2D" ,"N_{hits} 2D ;N_{hits} 2D ;tracks",40,-0.5,39.5);
1324  h_probeHitsInTIB_ = ProbeFeatures.make<TH1F>("h_probeNRechitsTIB" ,"N_{hits} TIB ;N_{hits} TIB;tracks",40,-0.5,39.5);
1325  h_probeHitsInTOB_ = ProbeFeatures.make<TH1F>("h_probeNRechitsTOB" ,"N_{hits} TOB ;N_{hits} TOB;tracks",40,-0.5,39.5);
1326  h_probeHitsInTID_ = ProbeFeatures.make<TH1F>("h_probeNRechitsTID" ,"N_{hits} TID ;N_{hits} TID;tracks",40,-0.5,39.5);
1327  h_probeHitsInTEC_ = ProbeFeatures.make<TH1F>("h_probeNRechitsTEC" ,"N_{hits} TEC ;N_{hits} TEC;tracks",40,-0.5,39.5);
1328  h_probeHitsInBPIX_ = ProbeFeatures.make<TH1F>("h_probeNRechitsBPIX","N_{hits} BPIX;N_{hits} BPIX;tracks",40,-0.5,39.5);
1329  h_probeHitsInFPIX_ = ProbeFeatures.make<TH1F>("h_probeNRechitsFPIX","N_{hits} FPIX;N_{hits} FPIX;tracks",40,-0.5,39.5);
1330 
1331  h_probeL1Ladder_ = ProbeFeatures.make<TH1F>("h_probeL1Ladder","Ladder number (L1 hit); ladder number",22,-1.5,20.5);
1332  h_probeL1Module_ = ProbeFeatures.make<TH1F>("h_probeL1Module","Module number (L1 hit); module number",10,-1.5,8.5);
1333  h_probeHasBPixL1Overlap_ = ProbeFeatures.make<TH1I>("h_probeHasBPixL1Overlap","n. hits in L1;n. L1-BPix hits;tracks",5,0,5);
1334 
1335  // refit vertex features
1336  TFileDirectory RefitVertexFeatures = fs->mkdir("RefitVertexFeatures");
1337  h_fitVtxNtracks_ = RefitVertexFeatures.make<TH1F>("h_fitVtxNtracks" ,"N_{trks} used in vertex fit;N^{fit}_{tracks};vertices" ,100,-0.5,99.5);
1338  h_fitVtxNdof_ = RefitVertexFeatures.make<TH1F>("h_fitVtxNdof" ,"N_{DOF} of vertex fit;N_{DOF} of refit vertex;vertices" ,100,-0.5,99.5);
1339  h_fitVtxChi2_ = RefitVertexFeatures.make<TH1F>("h_fitVtxChi2" ,"#chi^{2} of vertex fit;vertex #chi^{2};vertices" ,100,-0.5,99.5);
1340  h_fitVtxChi2ndf_ = RefitVertexFeatures.make<TH1F>("h_fitVtxChi2ndf" ,"#chi^{2}/ndf of vertex fit;vertex #chi^{2}/ndf;vertices" ,100,-0.5,9.5);
1341  h_fitVtxChi2Prob_ = RefitVertexFeatures.make<TH1F>("h_fitVtxChi2Prob" ,"Prob(#chi^{2},ndf) of vertex fit;Prob(#chi^{2},ndf);vertices",40,0.,1.);
1342  h_fitVtxTrackWeights_ = RefitVertexFeatures.make<TH1F>("h_fitVtxTrackWeights","track weights associated to track;track weights;tracks",40,0.,1.);
1343  h_fitVtxTrackAverageWeight_ = RefitVertexFeatures.make<TH1F>("h_fitVtxTrackAverageWeight_","average track weight per vertex;#LT track weight #GT;vertices",40,0.,1.);
1344 
1345  if(useTracksFromRecoVtx_) {
1346 
1347  TFileDirectory RecoVertexFeatures = fs->mkdir("RecoVertexFeatures");
1348  h_recoVtxNtracks_ = RecoVertexFeatures.make<TH1F>("h_recoVtxNtracks" ,"N^{vtx}_{trks};N^{vtx}_{trks};vertices" ,100,-0.5,99.5);
1349  h_recoVtxChi2ndf_ = RecoVertexFeatures.make<TH1F>("h_recoVtxChi2ndf" ,"#chi^{2}/ndf vtx;#chi^{2}/ndf vtx;vertices" ,10,-0.5,9.5);
1350  h_recoVtxChi2Prob_ = RecoVertexFeatures.make<TH1F>("h_recoVtxChi2Prob" ,"Prob(#chi^{2},ndf);Prob(#chi^{2},ndf);vertices",40,0.,1.);
1351  h_recoVtxSumPt_ = RecoVertexFeatures.make<TH1F>("h_recoVtxSumPt" ,"Sum(p^{trks}_{T});Sum(p^{trks}_{T});vertices" ,100,0.,200.);
1352 
1353  }
1354 
1355 
1356  TFileDirectory DA = fs->mkdir("DA");
1357  //DA.cd();
1358  hDA=bookVertexHistograms(DA);
1359  //for(std::map<std::string,TH1*>::const_iterator hist=hDA.begin(); hist!=hDA.end(); hist++){
1360  //hist->second->SetDirectory(DA);
1361  // DA.make<TH1F>(hist->second);
1362  // }
1363 
1364  // initialize the residuals histograms
1365 
1366  const float dxymax_phi = theDetails_.getHigh(PVValHelper::dxy,PVValHelper::phi);
1367  const float dzmax_phi = theDetails_.getHigh(PVValHelper::dz ,PVValHelper::eta);
1368  const float dxymax_eta = theDetails_.getHigh(PVValHelper::dxy,PVValHelper::phi);
1369  const float dzmax_eta = theDetails_.getHigh(PVValHelper::dz, PVValHelper::eta);
1370  //const float d3Dmax_phi = theDetails_.getHigh(PVValHelper::d3D,PVValHelper::phi);
1371  const float d3Dmax_eta = theDetails_.getHigh(PVValHelper::d3D,PVValHelper::eta);
1372 
1374  //
1375  // Unbiased track-to-vertex residuals
1376  // The vertex is refit without the probe track
1377  //
1379 
1380  // _ _ _ _ ___ _ _ _
1381  // /_\ | |__ ___ ___| |_ _| |_ ___ | _ \___ __(_)__| |_ _ __ _| |___
1382  // / _ \| '_ (_-</ _ \ | || | _/ -_) | / -_|_-< / _` | || / _` | (_-<
1383  // /_/ \_\_.__/__/\___/_|\_,_|\__\___| |_|_\___/__/_\__,_|\_,_\__,_|_/__/
1384  //
1385 
1386  TFileDirectory AbsTransPhiRes = fs->mkdir("Abs_Transv_Phi_Residuals");
1387  a_dxyPhiResiduals = bookResidualsHistogram(AbsTransPhiRes,nBins_,PVValHelper::dxy,PVValHelper::phi);
1388  a_dxPhiResiduals = bookResidualsHistogram(AbsTransPhiRes,nBins_,PVValHelper::dx,PVValHelper::phi);
1389  a_dyPhiResiduals = bookResidualsHistogram(AbsTransPhiRes,nBins_,PVValHelper::dy,PVValHelper::phi);
1390  a_IP2DPhiResiduals = bookResidualsHistogram(AbsTransPhiRes,nBins_,PVValHelper::IP2D,PVValHelper::phi);
1391 
1392  TFileDirectory AbsTransEtaRes = fs->mkdir("Abs_Transv_Eta_Residuals");
1393  a_dxyEtaResiduals = bookResidualsHistogram(AbsTransEtaRes,nBins_,PVValHelper::dxy,PVValHelper::eta);
1394  a_dxEtaResiduals = bookResidualsHistogram(AbsTransEtaRes,nBins_,PVValHelper::dx,PVValHelper::eta);
1395  a_dyEtaResiduals = bookResidualsHistogram(AbsTransEtaRes,nBins_,PVValHelper::dy,PVValHelper::eta);
1396  a_IP2DEtaResiduals = bookResidualsHistogram(AbsTransEtaRes,nBins_,PVValHelper::IP2D,PVValHelper::eta);
1397 
1398  TFileDirectory AbsLongPhiRes = fs->mkdir("Abs_Long_Phi_Residuals");
1399  a_dzPhiResiduals = bookResidualsHistogram(AbsLongPhiRes,nBins_,PVValHelper::dz,PVValHelper::phi);
1400  a_reszPhiResiduals = bookResidualsHistogram(AbsLongPhiRes,nBins_,PVValHelper::resz,PVValHelper::phi);
1401 
1402  TFileDirectory AbsLongEtaRes = fs->mkdir("Abs_Long_Eta_Residuals");
1403  a_dzEtaResiduals = bookResidualsHistogram(AbsLongEtaRes,nBins_,PVValHelper::dz,PVValHelper::eta);
1404  a_reszEtaResiduals = bookResidualsHistogram(AbsLongEtaRes,nBins_,PVValHelper::resz,PVValHelper::eta);
1405 
1406  TFileDirectory Abs3DPhiRes = fs->mkdir("Abs_3D_Phi_Residuals");
1407  a_d3DPhiResiduals = bookResidualsHistogram(Abs3DPhiRes,nBins_,PVValHelper::d3D,PVValHelper::phi);
1408  a_IP3DPhiResiduals = bookResidualsHistogram(Abs3DPhiRes,nBins_,PVValHelper::IP3D,PVValHelper::phi);
1409 
1410  TFileDirectory Abs3DEtaRes = fs->mkdir("Abs_3D_Eta_Residuals");
1411  a_d3DEtaResiduals = bookResidualsHistogram(Abs3DEtaRes ,nBins_,PVValHelper::d3D,PVValHelper::eta);
1412  a_IP3DEtaResiduals = bookResidualsHistogram(Abs3DEtaRes ,nBins_,PVValHelper::IP3D,PVValHelper::eta);
1413 
1414  TFileDirectory NormTransPhiRes = fs->mkdir("Norm_Transv_Phi_Residuals");
1415  n_dxyPhiResiduals = bookResidualsHistogram(NormTransPhiRes,nBins_,PVValHelper::norm_dxy,PVValHelper::phi,true);
1416  n_IP2DPhiResiduals = bookResidualsHistogram(NormTransPhiRes,nBins_,PVValHelper::norm_IP2D,PVValHelper::phi,true);
1417 
1418  TFileDirectory NormTransEtaRes = fs->mkdir("Norm_Transv_Eta_Residuals");
1419  n_dxyEtaResiduals = bookResidualsHistogram(NormTransEtaRes,nBins_,PVValHelper::norm_dxy,PVValHelper::eta,true);
1420  n_IP2DEtaResiduals = bookResidualsHistogram(NormTransEtaRes,nBins_,PVValHelper::norm_IP2D,PVValHelper::eta,true);
1421 
1422  TFileDirectory NormLongPhiRes = fs->mkdir("Norm_Long_Phi_Residuals");
1423  n_dzPhiResiduals = bookResidualsHistogram(NormLongPhiRes,nBins_,PVValHelper::norm_dz,PVValHelper::phi,true);
1424  n_reszPhiResiduals = bookResidualsHistogram(NormLongPhiRes,nBins_,PVValHelper::norm_resz,PVValHelper::phi,true);
1425 
1426  TFileDirectory NormLongEtaRes = fs->mkdir("Norm_Long_Eta_Residuals");
1427  n_dzEtaResiduals = bookResidualsHistogram(NormLongEtaRes,nBins_,PVValHelper::norm_dz,PVValHelper::eta,true);
1428  n_reszEtaResiduals = bookResidualsHistogram(NormLongEtaRes,nBins_,PVValHelper::norm_resz,PVValHelper::eta,true);
1429 
1430  TFileDirectory Norm3DPhiRes = fs->mkdir("Norm_3D_Phi_Residuals");
1431  n_d3DPhiResiduals = bookResidualsHistogram(Norm3DPhiRes,nBins_,PVValHelper::norm_d3D,PVValHelper::phi,true);
1432  n_IP3DPhiResiduals = bookResidualsHistogram(Norm3DPhiRes,nBins_,PVValHelper::norm_IP3D,PVValHelper::phi,true);
1433 
1434  TFileDirectory Norm3DEtaRes = fs->mkdir("Norm_3D_Eta_Residuals");
1435  n_d3DEtaResiduals = bookResidualsHistogram(Norm3DEtaRes,nBins_,PVValHelper::norm_d3D,PVValHelper::eta,true);
1436  n_IP3DEtaResiduals = bookResidualsHistogram(Norm3DEtaRes,nBins_,PVValHelper::norm_IP3D,PVValHelper::eta,true);
1437 
1438  TFileDirectory AbsDoubleDiffRes = fs->mkdir("Abs_DoubleDiffResiduals");
1439  TFileDirectory NormDoubleDiffRes = fs->mkdir("Norm_DoubleDiffResiduals");
1440 
1441  // book residuals vs pT histograms
1442 
1443  TFileDirectory AbsTranspTRes = fs->mkdir("Abs_Transv_pT_Residuals");
1444  h_dxy_pT_ = bookResidualsHistogram(AbsTranspTRes,nPtBins_,PVValHelper::dxy,PVValHelper::pT);
1445 
1446  TFileDirectory AbsLongpTRes = fs->mkdir("Abs_Long_pT_Residuals");
1447  h_dz_pT_ = bookResidualsHistogram(AbsLongpTRes,nPtBins_,PVValHelper::dz,PVValHelper::pT);
1448 
1449  TFileDirectory NormTranspTRes = fs->mkdir("Norm_Transv_pT_Residuals");
1450  h_norm_dxy_pT_ = bookResidualsHistogram(NormTranspTRes,nPtBins_,PVValHelper::norm_dxy,PVValHelper::pT,true);
1451 
1452  TFileDirectory NormLongpTRes = fs->mkdir("Norm_Long_pT_Residuals");
1453  h_norm_dz_pT_ = bookResidualsHistogram(NormLongpTRes,nPtBins_,PVValHelper::norm_dz,PVValHelper::pT,true);
1454 
1455  // book residuals vs pT histograms in central region (|eta|<1.0)
1456 
1457  TFileDirectory AbsTranspTCentralRes = fs->mkdir("Abs_Transv_pTCentral_Residuals");
1458  h_dxy_Central_pT_ = bookResidualsHistogram(AbsTranspTCentralRes,nPtBins_,PVValHelper::dxy,PVValHelper::pTCentral);
1459 
1460  TFileDirectory AbsLongpTCentralRes = fs->mkdir("Abs_Long_pTCentral_Residuals");
1461  h_dz_Central_pT_ = bookResidualsHistogram(AbsLongpTCentralRes,nPtBins_,PVValHelper::dz,PVValHelper::pTCentral);
1462 
1463  TFileDirectory NormTranspTCentralRes = fs->mkdir("Norm_Transv_pTCentral_Residuals");
1464  h_norm_dxy_Central_pT_ = bookResidualsHistogram(NormTranspTCentralRes,nPtBins_,PVValHelper::norm_dxy,PVValHelper::pTCentral,true);
1465 
1466  TFileDirectory NormLongpTCentralRes = fs->mkdir("Norm_Long_pTCentral_Residuals");
1467  h_norm_dz_Central_pT_ = bookResidualsHistogram(NormLongpTCentralRes,nPtBins_,PVValHelper::norm_dz,PVValHelper::pTCentral,true);
1468 
1469  // book residuals vs module number
1470 
1471  TFileDirectory AbsTransModZRes = fs->mkdir("Abs_Transv_modZ_Residuals");
1472  h_dxy_modZ_ = bookResidualsHistogram(AbsTransModZRes,8,PVValHelper::dxy,PVValHelper::modZ);
1473 
1474  TFileDirectory AbsLongModZRes = fs->mkdir("Abs_Long_modZ_Residuals");
1475  h_dz_modZ_ = bookResidualsHistogram(AbsLongModZRes,8,PVValHelper::dz,PVValHelper::modZ);
1476 
1477 
1478  // _ _ _ _ _ ___ _ _ _
1479  // | \| |___ _ _ _ __ __ _| (_)______ __| | | _ \___ __(_)__| |_ _ __ _| |___
1480  // | .` / _ \ '_| ' \/ _` | | |_ / -_) _` | | / -_|_-< / _` | || / _` | (_-<
1481  // |_|\_\___/_| |_|_|_\__,_|_|_/__\___\__,_| |_|_\___/__/_\__,_|\_,_\__,_|_/__/
1482  //
1483 
1484  TFileDirectory NormTransModZRes = fs->mkdir("Norm_Transv_modZ_Residuals");
1485  h_norm_dxy_modZ_ = bookResidualsHistogram(NormTransModZRes,8,PVValHelper::norm_dxy,PVValHelper::modZ,true);
1486 
1487  TFileDirectory NormLongModZRes = fs->mkdir("Norm_Long_modZ_Residuals");
1488  h_norm_dz_modZ_ = bookResidualsHistogram(NormLongModZRes,8,PVValHelper::norm_dz,PVValHelper::modZ,true);
1489 
1490  TFileDirectory AbsTransLadderRes = fs->mkdir("Abs_Transv_ladder_Residuals");
1491  h_dxy_ladder_ = bookResidualsHistogram(AbsTransLadderRes,nLadders_,PVValHelper::dxy,PVValHelper::ladder);
1492 
1493  TFileDirectory AbsTransLadderResOverlap = fs->mkdir("Abs_Transv_ladderOverlap_Residuals");
1494  h_dxy_ladderOverlap_ = bookResidualsHistogram(AbsTransLadderResOverlap,nLadders_,PVValHelper::dxy,PVValHelper::ladder);
1495 
1496  TFileDirectory AbsTransLadderResNoOverlap = fs->mkdir("Abs_Transv_ladderNoOverlap_Residuals");
1497  h_dxy_ladderNoOverlap_ = bookResidualsHistogram(AbsTransLadderResNoOverlap,nLadders_,PVValHelper::dxy,PVValHelper::ladder);
1498 
1499  TFileDirectory AbsLongLadderRes = fs->mkdir("Abs_Long_ladder_Residuals");
1500  h_dz_ladder_ = bookResidualsHistogram(AbsLongLadderRes,nLadders_,PVValHelper::dz,PVValHelper::ladder);
1501 
1502  TFileDirectory NormTransLadderRes = fs->mkdir("Norm_Transv_ladder_Residuals");
1503  h_norm_dxy_ladder_ = bookResidualsHistogram(NormTransLadderRes,nLadders_,PVValHelper::norm_dxy,PVValHelper::ladder,true);
1504 
1505  TFileDirectory NormLongLadderRes = fs->mkdir("Norm_Long_ladder_Residuals");
1506  h_norm_dz_ladder_ = bookResidualsHistogram(NormLongLadderRes,nLadders_,PVValHelper::norm_dz,PVValHelper::ladder,true);
1507 
1508 
1509  // book residuals as function of phi and eta
1510 
1511  for (int i=0; i<nBins_; ++i ) {
1512 
1513  float phiF = theDetails_.trendbins[PVValHelper::phi][i];
1514  float phiL = theDetails_.trendbins[PVValHelper::phi][i+1];
1515 
1516  // ___ _ _ ___ _ __ __ ___ _ _ _
1517  // | \ ___ _ _| |__| |___| \(_)/ _|/ _| | _ \___ __(_)__| |_ _ __ _| |___
1518  // | |) / _ \ || | '_ \ / -_) |) | | _| _| | / -_|_-< / _` | || / _` | (_-<
1519  // |___/\___/\_,_|_.__/_\___|___/|_|_| |_| |_|_\___/__/_\__,_|\_,_\__,_|_/__/
1520 
1521  for ( int j=0; j<nBins_; ++j ) {
1522 
1523  float etaF = theDetails_.trendbins[PVValHelper::eta][j];
1524  float etaL = theDetails_.trendbins[PVValHelper::eta][j+1];
1525 
1526  a_dxyResidualsMap[i][j] = AbsDoubleDiffRes.make<TH1F>(Form("histo_dxy_eta_plot%i_phi_plot%i",i,j),
1527  Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{xy};tracks",etaF,etaL,phiF,phiL),
1528  theDetails_.histobins,-dzmax_eta,dzmax_eta);
1529 
1530  a_dzResidualsMap[i][j] = AbsDoubleDiffRes.make<TH1F>(Form("histo_dz_eta_plot%i_phi_plot%i",i,j),
1531  Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{z};tracks",etaF,etaL,phiF,phiL),
1532  theDetails_.histobins,-dzmax_eta,dzmax_eta);
1533 
1534  a_d3DResidualsMap[i][j] = AbsDoubleDiffRes.make<TH1F>(Form("histo_d3D_eta_plot%i_phi_plot%i",i,j),
1535  Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{3D};tracks",etaF,etaL,phiF,phiL),
1536  theDetails_.histobins,0.,d3Dmax_eta);
1537 
1538  n_dxyResidualsMap[i][j] = NormDoubleDiffRes.make<TH1F>(Form("histo_norm_dxy_eta_plot%i_phi_plot%i",i,j),
1539  Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{xy}/#sigma_{d_{xy}};tracks",etaF,etaL,phiF,phiL),
1540  theDetails_.histobins,-dzmax_eta/100,dzmax_eta/100);
1541 
1542  n_dzResidualsMap[i][j] = NormDoubleDiffRes.make<TH1F>(Form("histo_norm_dz_eta_plot%i_phi_plot%i",i,j),
1543  Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{z}/#sigma_{d_{z}};tracks",etaF,etaL,phiF,phiL),
1544  theDetails_.histobins,-dzmax_eta/100,dzmax_eta/100);
1545 
1546  n_d3DResidualsMap[i][j] = NormDoubleDiffRes.make<TH1F>(Form("histo_norm_d3D_eta_plot%i_phi_plot%i",i,j),
1547  Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{3D}/#sigma_{d_{3D}};tracks",etaF,etaL,phiF,phiL),
1548  theDetails_.histobins,0.,d3Dmax_eta);
1549 
1550  }
1551  }
1552 
1553  // declaration of the directories
1554 
1555  TFileDirectory BiasVsParameter = fs->mkdir("BiasVsParameter");
1556 
1557  a_dxyVsPhi = BiasVsParameter.make<TH2F>("h2_dxy_vs_phi","d_{xy} vs track #phi;track #phi [rad];track d_{xy}(PV) [#mum]",
1558  nBins_,-M_PI,M_PI,theDetails_.histobins,-dxymax_phi,dxymax_phi);
1559 
1560  a_dzVsPhi = BiasVsParameter.make<TH2F>("h2_dz_vs_phi","d_{z} vs track #phi;track #phi [rad];track d_{z}(PV) [#mum]",
1561  nBins_,-M_PI,M_PI,theDetails_.histobins,-dzmax_phi,dzmax_phi);
1562 
1563  n_dxyVsPhi = BiasVsParameter.make<TH2F>("h2_n_dxy_vs_phi","d_{xy}/#sigma_{d_{xy}} vs track #phi;track #phi [rad];track d_{xy}(PV)/#sigma_{d_{xy}}",
1564  nBins_,-M_PI,M_PI,theDetails_.histobins,-dxymax_phi/100.,dxymax_phi/100.);
1565 
1566  n_dzVsPhi = BiasVsParameter.make<TH2F>("h2_n_dz_vs_phi","d_{z}/#sigma_{d_{z}} vs track #phi;track #phi [rad];track d_{z}(PV)/#sigma_{d_{z}}",
1567  nBins_,-M_PI,M_PI,theDetails_.histobins,-dzmax_phi/100.,dzmax_phi/100.);
1568 
1569  a_dxyVsEta = BiasVsParameter.make<TH2F>("h2_dxy_vs_eta","d_{xy} vs track #eta;track #eta;track d_{xy}(PV) [#mum]",
1570  nBins_,-etaOfProbe_,etaOfProbe_,theDetails_.histobins,-dxymax_eta,dzmax_eta);
1571 
1572  a_dzVsEta = BiasVsParameter.make<TH2F>("h2_dz_vs_eta","d_{z} vs track #eta;track #eta;track d_{z}(PV) [#mum]",
1573  nBins_,-etaOfProbe_,etaOfProbe_,theDetails_.histobins,-dzmax_eta,dzmax_eta);
1574 
1575  n_dxyVsEta = BiasVsParameter.make<TH2F>("h2_n_dxy_vs_eta","d_{xy}/#sigma_{d_{xy}} vs track #eta;track #eta;track d_{xy}(PV)/#sigma_{d_{xy}}",
1576  nBins_,-etaOfProbe_,etaOfProbe_,theDetails_.histobins,-dxymax_eta/100.,dxymax_eta/100.);
1577 
1578  n_dzVsEta = BiasVsParameter.make<TH2F>("h2_n_dz_vs_eta","d_{z}/#sigma_{d_{z}} vs track #eta;track #eta;track d_{z}(PV)/#sigma_{d_{z}}",
1579  nBins_,-etaOfProbe_,etaOfProbe_,theDetails_.histobins,-dzmax_eta/100.,dzmax_eta/100.);
1580 
1581  MeanTrendsDir = fs->mkdir("MeanTrends");
1582  WidthTrendsDir = fs->mkdir("WidthTrends");
1583  MedianTrendsDir = fs->mkdir("MedianTrends");
1584  MADTrendsDir = fs->mkdir("MADTrends");
1585 
1586  Mean2DMapsDir = fs->mkdir("MeanMaps");
1587  Width2DMapsDir = fs->mkdir("WidthMaps");
1588 
1589  double highedge=nBins_-0.5;
1590  double lowedge=-0.5;
1591 
1592  // means and widths from the fit
1593 
1594  a_dxyPhiMeanTrend = MeanTrendsDir.make<TH1F> ("means_dxy_phi",
1595  "#LT d_{xy} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{xy} #GT [#mum]",
1596  nBins_,lowedge,highedge);
1597 
1598  a_dxyPhiWidthTrend = WidthTrendsDir.make<TH1F>("widths_dxy_phi",
1599  "#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees];#sigma_{d_{xy}} [#mum]",
1600  nBins_,lowedge,highedge);
1601 
1602  a_dzPhiMeanTrend = MeanTrendsDir.make<TH1F> ("means_dz_phi",
1603  "#LT d_{z} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{z} #GT [#mum]",
1604  nBins_,lowedge,highedge);
1605 
1606  a_dzPhiWidthTrend = WidthTrendsDir.make<TH1F>("widths_dz_phi","#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];#sigma_{d_{z}} [#mum]",
1607  nBins_,lowedge,highedge);
1608 
1609  a_dxyEtaMeanTrend = MeanTrendsDir.make<TH1F> ("means_dxy_eta",
1610  "#LT d_{xy} #GT vs #eta sector;#eta (sector);#LT d_{xy} #GT [#mum]",
1611  nBins_,lowedge,highedge);
1612 
1613  a_dxyEtaWidthTrend = WidthTrendsDir.make<TH1F>("widths_dxy_eta",
1614  "#sigma_{d_{xy}} vs #eta sector;#eta (sector);#sigma_{d_{xy}} [#mum]",
1615  nBins_,lowedge,highedge);
1616 
1617  a_dzEtaMeanTrend = MeanTrendsDir.make<TH1F> ("means_dz_eta",
1618  "#LT d_{z} #GT vs #eta sector;#eta (sector);#LT d_{z} #GT [#mum]"
1619  ,nBins_,lowedge,highedge);
1620 
1621  a_dzEtaWidthTrend = WidthTrendsDir.make<TH1F>("widths_dz_eta",
1622  "#sigma_{d_{xy}} vs #eta sector;#eta (sector);#sigma_{d_{z}} [#mum]",
1623  nBins_,lowedge,highedge);
1624 
1625  n_dxyPhiMeanTrend = MeanTrendsDir.make<TH1F> ("norm_means_dxy_phi",
1626  "#LT d_{xy}/#sigma_{d_{xy}} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{xy}/#sigma_{d_{xy}} #GT",
1627  nBins_,lowedge,highedge);
1628 
1629  n_dxyPhiWidthTrend = WidthTrendsDir.make<TH1F>("norm_widths_dxy_phi",
1630  "width(d_{xy}/#sigma_{d_{xy}}) vs #phi sector;#varphi (sector) [degrees]; width(d_{xy}/#sigma_{d_{xy}})",
1631  nBins_,lowedge,highedge);
1632 
1633  n_dzPhiMeanTrend = MeanTrendsDir.make<TH1F> ("norm_means_dz_phi",
1634  "#LT d_{z}/#sigma_{d_{z}} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{z}/#sigma_{d_{z}} #GT",
1635  nBins_,lowedge,highedge);
1636 
1637  n_dzPhiWidthTrend = WidthTrendsDir.make<TH1F>("norm_widths_dz_phi",
1638  "width(d_{z}/#sigma_{d_{z}}) vs #phi sector;#varphi (sector) [degrees];width(d_{z}/#sigma_{d_{z}})",
1639  nBins_,lowedge,highedge);
1640 
1641  n_dxyEtaMeanTrend = MeanTrendsDir.make<TH1F> ("norm_means_dxy_eta",
1642  "#LT d_{xy}/#sigma_{d_{xy}} #GT vs #eta sector;#eta (sector);#LT d_{xy}/#sigma_{d_{z}} #GT",
1643  nBins_,lowedge,highedge);
1644 
1645  n_dxyEtaWidthTrend = WidthTrendsDir.make<TH1F>("norm_widths_dxy_eta",
1646  "width(d_{xy}/#sigma_{d_{xy}}) vs #eta sector;#eta (sector);width(d_{xy}/#sigma_{d_{z}})",
1647  nBins_,lowedge,highedge);
1648 
1649  n_dzEtaMeanTrend = MeanTrendsDir.make<TH1F> ("norm_means_dz_eta",
1650  "#LT d_{z}/#sigma_{d_{z}} #GT vs #eta sector;#eta (sector);#LT d_{z}/#sigma_{d_{z}} #GT",
1651  nBins_,lowedge,highedge);
1652 
1653  n_dzEtaWidthTrend = WidthTrendsDir.make<TH1F>("norm_widths_dz_eta",
1654  "width(d_{z}/#sigma_{d_{z}}) vs #eta sector;#eta (sector);width(d_{z}/#sigma_{d_{z}})",
1655  nBins_,lowedge,highedge);
1656 
1657  // means and widhts vs pT and pTCentral
1658 
1659  a_dxypTMeanTrend = MeanTrendsDir.make<TH1F> ("means_dxy_pT",
1660  "#LT d_{xy} #GT vs pT;p_{T} [GeV];#LT d_{xy} #GT [#mum]",
1661  mypT_bins_.size()-1,mypT_bins_.data());
1662 
1663  a_dxypTWidthTrend = WidthTrendsDir.make<TH1F>("widths_dxy_pT",
1664  "#sigma_{d_{xy}} vs pT;p_{T} [GeV];#sigma_{d_{xy}} [#mum]",
1665  mypT_bins_.size()-1,mypT_bins_.data());
1666 
1667  a_dzpTMeanTrend = MeanTrendsDir.make<TH1F> ("means_dz_pT",
1668  "#LT d_{z} #GT vs pT;p_{T} [GeV];#LT d_{z} #GT [#mum]",
1669  mypT_bins_.size()-1,mypT_bins_.data());
1670 
1671  a_dzpTWidthTrend = WidthTrendsDir.make<TH1F>("widths_dz_pT","#sigma_{d_{z}} vs pT;p_{T} [GeV];#sigma_{d_{z}} [#mum]",
1672  mypT_bins_.size()-1,mypT_bins_.data());
1673 
1674 
1675  n_dxypTMeanTrend = MeanTrendsDir.make<TH1F> ("norm_means_dxy_pT",
1676  "#LT d_{xy}/#sigma_{d_{xy}} #GT vs pT;p_{T} [GeV];#LT d_{xy}/#sigma_{d_{xy}} #GT",
1677  mypT_bins_.size()-1,mypT_bins_.data());
1678 
1679  n_dxypTWidthTrend = WidthTrendsDir.make<TH1F>("norm_widths_dxy_pT",
1680  "width(d_{xy}/#sigma_{d_{xy}}) vs pT;p_{T} [GeV]; width(d_{xy}/#sigma_{d_{xy}})",
1681  mypT_bins_.size()-1,mypT_bins_.data());
1682 
1683  n_dzpTMeanTrend = MeanTrendsDir.make<TH1F> ("norm_means_dz_pT",
1684  "#LT d_{z}/#sigma_{d_{z}} #GT vs pT;p_{T} [GeV];#LT d_{z}/#sigma_{d_{z}} #GT",
1685  mypT_bins_.size()-1,mypT_bins_.data());
1686 
1687  n_dzpTWidthTrend = WidthTrendsDir.make<TH1F>("norm_widths_dz_pT",
1688  "width(d_{z}/#sigma_{d_{z}}) vs pT;p_{T} [GeV];width(d_{z}/#sigma_{d_{z}})",
1689  mypT_bins_.size()-1,mypT_bins_.data());
1690 
1691 
1692  a_dxypTCentralMeanTrend = MeanTrendsDir.make<TH1F> ("means_dxy_pTCentral",
1693  "#LT d_{xy} #GT vs p_{T};p_{T}(|#eta|<1.) [GeV];#LT d_{xy} #GT [#mum]",
1694  mypT_bins_.size()-1,mypT_bins_.data());
1695 
1696  a_dxypTCentralWidthTrend = WidthTrendsDir.make<TH1F>("widths_dxy_pTCentral",
1697  "#sigma_{d_{xy}} vs p_{T};p_{T}(|#eta|<1.) [GeV];#sigma_{d_{xy}} [#mum]",
1698  mypT_bins_.size()-1,mypT_bins_.data());
1699 
1700  a_dzpTCentralMeanTrend = MeanTrendsDir.make<TH1F> ("means_dz_pTCentral",
1701  "#LT d_{z} #GT vs p_{T};p_{T}(|#eta|<1.) [GeV];#LT d_{z} #GT [#mum]"
1702  ,mypT_bins_.size()-1,mypT_bins_.data());
1703 
1704  a_dzpTCentralWidthTrend = WidthTrendsDir.make<TH1F>("widths_dz_pTCentral",
1705  "#sigma_{d_{z}} vs p_{T};p_{T}(|#eta|<1.) [GeV];#sigma_{d_{z}} [#mum]",
1706  mypT_bins_.size()-1,mypT_bins_.data());
1707 
1708  n_dxypTCentralMeanTrend = MeanTrendsDir.make<TH1F> ("norm_means_dxy_pTCentral",
1709  "#LT d_{xy}/#sigma_{d_{xy}} #GT vs p_{T};p_{T}(|#eta|<1.) [GeV];#LT d_{xy}/#sigma_{d_{z}} #GT",
1710  mypT_bins_.size()-1,mypT_bins_.data());
1711 
1712  n_dxypTCentralWidthTrend = WidthTrendsDir.make<TH1F>("norm_widths_dxy_pTCentral",
1713  "width(d_{xy}/#sigma_{d_{xy}}) vs p_{T};p_{T}(|#eta|<1.) [GeV];width(d_{xy}/#sigma_{d_{z}})",
1714  mypT_bins_.size()-1,mypT_bins_.data());
1715 
1716  n_dzpTCentralMeanTrend = MeanTrendsDir.make<TH1F> ("norm_means_dz_pTCentral",
1717  "#LT d_{z}/#sigma_{d_{z}} #GT vs p_{T};p_{T}(|#eta|<1.) [GeV];#LT d_{z}/#sigma_{d_{z}} #GT",
1718  mypT_bins_.size()-1,mypT_bins_.data());
1719 
1720  n_dzpTCentralWidthTrend = WidthTrendsDir.make<TH1F>("norm_widths_dz_pTCentral",
1721  "width(d_{z}/#sigma_{d_{z}}) vs p_{T};p_{T}(|#eta|<1.) [GeV];width(d_{z}/#sigma_{d_{z}})",
1722  mypT_bins_.size()-1,mypT_bins_.data());
1723 
1724  // 2D maps
1725 
1726  a_dxyMeanMap = Mean2DMapsDir.make<TH2F> ("means_dxy_map",
1727  "#LT d_{xy} #GT map;#eta (sector);#varphi (sector) [degrees]",
1728  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1729 
1730  a_dzMeanMap = Mean2DMapsDir.make<TH2F> ("means_dz_map",
1731  "#LT d_{z} #GT map;#eta (sector);#varphi (sector) [degrees]",
1732  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1733 
1734  n_dxyMeanMap = Mean2DMapsDir.make<TH2F> ("norm_means_dxy_map",
1735  "#LT d_{xy}/#sigma_{d_{xy}} #GT map;#eta (sector);#varphi (sector) [degrees]",
1736  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1737 
1738  n_dzMeanMap = Mean2DMapsDir.make<TH2F> ("norm_means_dz_map",
1739  "#LT d_{z}/#sigma_{d_{z}} #GT map;#eta (sector);#varphi (sector) [degrees]",
1740  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1741 
1742  a_dxyWidthMap = Width2DMapsDir.make<TH2F> ("widths_dxy_map",
1743  "#sigma_{d_{xy}} map;#eta (sector);#varphi (sector) [degrees]",
1744  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1745 
1746  a_dzWidthMap = Width2DMapsDir.make<TH2F> ("widths_dz_map",
1747  "#sigma_{d_{z}} map;#eta (sector);#varphi (sector) [degrees]",
1748  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1749 
1750  n_dxyWidthMap = Width2DMapsDir.make<TH2F> ("norm_widths_dxy_map",
1751  "width(d_{xy}/#sigma_{d_{xy}}) map;#eta (sector);#varphi (sector) [degrees]",
1752  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1753 
1754  n_dzWidthMap = Width2DMapsDir.make<TH2F> ("norm_widths_dz_map",
1755  "width(d_{z}/#sigma_{d_{z}}) map;#eta (sector);#varphi (sector) [degrees]",
1756  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1757 
1758  // medians and MADs
1759 
1760  a_dxyPhiMedianTrend = MedianTrendsDir.make<TH1F>("medians_dxy_phi",
1761  "Median of d_{xy} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{xy}) [#mum]",
1762  nBins_,lowedge,highedge);
1763 
1764  a_dxyPhiMADTrend = MADTrendsDir.make<TH1F> ("MADs_dxy_phi",
1765  "Median absolute deviation of d_{xy} vs #phi sector;#varphi (sector) [degrees];MAD(d_{xy}) [#mum]",
1766  nBins_,lowedge,highedge);
1767 
1768  a_dzPhiMedianTrend = MedianTrendsDir.make<TH1F>("medians_dz_phi",
1769  "Median of d_{z} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{z}) [#mum]",
1770  nBins_,lowedge,highedge);
1771 
1772  a_dzPhiMADTrend = MADTrendsDir.make<TH1F> ("MADs_dz_phi",
1773  "Median absolute deviation of d_{z} vs #phi sector;#varphi (sector) [degrees];MAD(d_{z}) [#mum]",
1774  nBins_,lowedge,highedge);
1775 
1776  a_dxyEtaMedianTrend = MedianTrendsDir.make<TH1F>("medians_dxy_eta",
1777  "Median of d_{xy} vs #eta sector;#eta (sector);#mu_{1/2}(d_{xy}) [#mum]",
1778  nBins_,lowedge,highedge);
1779 
1780  a_dxyEtaMADTrend = MADTrendsDir.make<TH1F> ("MADs_dxy_eta",
1781  "Median absolute deviation of d_{xy} vs #eta sector;#eta (sector);MAD(d_{xy}) [#mum]",
1782  nBins_,lowedge,highedge);
1783 
1784  a_dzEtaMedianTrend = MedianTrendsDir.make<TH1F>("medians_dz_eta",
1785  "Median of d_{z} vs #eta sector;#eta (sector);#mu_{1/2}(d_{z}) [#mum]",
1786  nBins_,lowedge,highedge);
1787 
1788  a_dzEtaMADTrend = MADTrendsDir.make<TH1F> ("MADs_dz_eta",
1789  "Median absolute deviation of d_{z} vs #eta sector;#eta (sector);MAD(d_{z}) [#mum]",
1790  nBins_,lowedge,highedge);
1791 
1792  n_dxyPhiMedianTrend = MedianTrendsDir.make<TH1F>("norm_medians_dxy_phi",
1793  "Median of d_{xy}/#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{xy}/#sigma_{d_{xy}})",
1794  nBins_,lowedge,highedge);
1795 
1796  n_dxyPhiMADTrend = MADTrendsDir.make<TH1F> ("norm_MADs_dxy_phi",
1797  "Median absolute deviation of d_{xy}/#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees]; MAD(d_{xy}/#sigma_{d_{xy}})",
1798  nBins_,lowedge,highedge);
1799 
1800  n_dzPhiMedianTrend = MedianTrendsDir.make<TH1F>("norm_medians_dz_phi",
1801  "Median of d_{z}/#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{z}/#sigma_{d_{z}})",
1802  nBins_,lowedge,highedge);
1803 
1804  n_dzPhiMADTrend = MADTrendsDir.make<TH1F> ("norm_MADs_dz_phi",
1805  "Median absolute deviation of d_{z}/#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];MAD(d_{z}/#sigma_{d_{z}})",
1806  nBins_,lowedge,highedge);
1807 
1808  n_dxyEtaMedianTrend = MedianTrendsDir.make<TH1F>("norm_medians_dxy_eta",
1809  "Median of d_{xy}/#sigma_{d_{xy}} vs #eta sector;#eta (sector);#mu_{1/2}(d_{xy}/#sigma_{d_{z}})",
1810  nBins_,lowedge,highedge);
1811 
1812  n_dxyEtaMADTrend = MADTrendsDir.make<TH1F> ("norm_MADs_dxy_eta",
1813  "Median absolute deviation of d_{xy}/#sigma_{d_{xy}} vs #eta sector;#eta (sector);MAD(d_{xy}/#sigma_{d_{z}})",
1814  nBins_,lowedge,highedge);
1815 
1816  n_dzEtaMedianTrend = MedianTrendsDir.make<TH1F>("norm_medians_dz_eta",
1817  "Median of d_{z}/#sigma_{d_{z}} vs #eta sector;#eta (sector);#mu_{1/2}(d_{z}/#sigma_{d_{z}})",
1818  nBins_,lowedge,highedge);
1819 
1820  n_dzEtaMADTrend = MADTrendsDir.make<TH1F> ("norm_MADs_dz_eta",
1821  "Median absolute deviation of d_{z}/#sigma_{d_{z}} vs #eta sector;#eta (sector);MAD(d_{z}/#sigma_{d_{z}})",
1822  nBins_,lowedge,highedge);
1823 
1824 
1826  //
1827  // plots of biased residuals
1828  // The vertex includes the probe track
1829  //
1831 
1832  if (useTracksFromRecoVtx_){
1833 
1834  TFileDirectory AbsTransPhiBiasRes = fs->mkdir("Abs_Transv_Phi_BiasResiduals");
1835  a_dxyPhiBiasResiduals = bookResidualsHistogram(AbsTransPhiBiasRes,nBins_,PVValHelper::dxy,PVValHelper::phi);
1836 
1837  TFileDirectory AbsTransEtaBiasRes = fs->mkdir("Abs_Transv_Eta_BiasResiduals");
1838  a_dxyEtaBiasResiduals = bookResidualsHistogram(AbsTransEtaBiasRes,nBins_,PVValHelper::dxy,PVValHelper::eta);
1839 
1840  TFileDirectory AbsLongPhiBiasRes = fs->mkdir("Abs_Long_Phi_BiasResiduals");
1841  a_dzPhiBiasResiduals = bookResidualsHistogram(AbsLongPhiBiasRes,nBins_,PVValHelper::dz,PVValHelper::phi);
1842 
1843  TFileDirectory AbsLongEtaBiasRes = fs->mkdir("Abs_Long_Eta_BiasResiduals");
1844  a_dzEtaBiasResiduals = bookResidualsHistogram(AbsLongEtaBiasRes,nBins_,PVValHelper::dz,PVValHelper::eta);
1845 
1846  TFileDirectory NormTransPhiBiasRes = fs->mkdir("Norm_Transv_Phi_BiasResiduals");
1847  n_dxyPhiBiasResiduals = bookResidualsHistogram(NormTransPhiBiasRes,nBins_,PVValHelper::dxy,PVValHelper::phi);
1848 
1849  TFileDirectory NormTransEtaBiasRes = fs->mkdir("Norm_Transv_Eta_BiasResiduals");
1850  n_dxyEtaBiasResiduals = bookResidualsHistogram(NormTransEtaBiasRes,nBins_,PVValHelper::dxy,PVValHelper::eta);
1851 
1852  TFileDirectory NormLongPhiBiasRes = fs->mkdir("Norm_Long_Phi_BiasResiduals");
1853  n_dzPhiBiasResiduals = bookResidualsHistogram(NormLongPhiBiasRes,nBins_,PVValHelper::dz,PVValHelper::phi);
1854 
1855  TFileDirectory NormLongEtaBiasRes = fs->mkdir("Norm_Long_Eta_BiasResiduals");
1856  n_dzEtaBiasResiduals = bookResidualsHistogram(NormLongEtaBiasRes,nBins_,PVValHelper::dz,PVValHelper::eta);
1857 
1858  TFileDirectory AbsDoubleDiffBiasRes = fs->mkdir("Abs_DoubleDiffBiasResiduals");
1859  TFileDirectory NormDoubleDiffBiasRes = fs->mkdir("Norm_DoubleDiffBiasResiduals");
1860 
1861  for ( int i=0; i<nBins_; ++i ) {
1862 
1863  float phiF = theDetails_.trendbins[PVValHelper::phi][i];
1864  float phiL = theDetails_.trendbins[PVValHelper::phi][i+1];
1865 
1866  for ( int j=0; j<nBins_; ++j ) {
1867 
1868  float etaF = theDetails_.trendbins[PVValHelper::eta][j];
1869  float etaL = theDetails_.trendbins[PVValHelper::eta][j+1];
1870 
1871  a_dxyBiasResidualsMap[i][j] = AbsDoubleDiffBiasRes.make<TH1F>(Form("histo_dxy_eta_plot%i_phi_plot%i",i,j),
1872  Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{xy} [#mum];tracks",etaF,etaL,phiF,phiL),
1873  theDetails_.histobins,-dzmax_eta,dzmax_eta);
1874 
1875  a_dzBiasResidualsMap[i][j] = AbsDoubleDiffBiasRes.make<TH1F>(Form("histo_dxy_eta_plot%i_phi_plot%i",i,j),
1876  Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{z} [#mum];tracks",etaF,etaL,phiF,phiL),
1877  theDetails_.histobins,-dzmax_eta,dzmax_eta);
1878 
1879  n_dxyBiasResidualsMap[i][j] = NormDoubleDiffBiasRes.make<TH1F>(Form("histo_norm_dxy_eta_plot%i_phi_plot%i",i,j),
1880  Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{xy}/#sigma_{d_{xy}};tracks",etaF,etaL,phiF,phiL),
1881  theDetails_.histobins,-dzmax_eta/100,dzmax_eta/100);
1882 
1883  n_dzBiasResidualsMap[i][j] = NormDoubleDiffBiasRes.make<TH1F>(Form("histo_norm_dxy_eta_plot%i_phi_plot%i",i,j),
1884  Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{z}/#sigma_{d_{z}};tracks",etaF,etaL,phiF,phiL),
1885  theDetails_.histobins,-dzmax_eta/100,dzmax_eta/100);
1886  }
1887  }
1888 
1889  // declaration of the directories
1890 
1891  TFileDirectory MeanBiasTrendsDir = fs->mkdir("MeanBiasTrends");
1892  TFileDirectory WidthBiasTrendsDir = fs->mkdir("WidthBiasTrends");
1893  TFileDirectory MedianBiasTrendsDir = fs->mkdir("MedianBiasTrends");
1894  TFileDirectory MADBiasTrendsDir = fs->mkdir("MADBiasTrends");
1895 
1896  TFileDirectory Mean2DBiasMapsDir = fs->mkdir("MeanBiasMaps");
1897  TFileDirectory Width2DBiasMapsDir = fs->mkdir("WidthBiasMaps");
1898 
1899  // means and widths from the fit
1900 
1901  a_dxyPhiMeanBiasTrend = MeanBiasTrendsDir.make<TH1F> ("means_dxy_phi",
1902  "#LT d_{xy} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{xy} #GT [#mum]",
1903  nBins_,lowedge,highedge);
1904 
1905  a_dxyPhiWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>("widths_dxy_phi",
1906  "#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees];#sigma_{d_{xy}} [#mum]",
1907  nBins_,lowedge,highedge);
1908 
1909  a_dzPhiMeanBiasTrend = MeanBiasTrendsDir.make<TH1F> ("means_dz_phi",
1910  "#LT d_{z} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{z} #GT [#mum]",
1911  nBins_,lowedge,highedge);
1912 
1913  a_dzPhiWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>("widths_dz_phi",
1914  "#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];#sigma_{d_{z}} [#mum]",
1915  nBins_,lowedge,highedge);
1916 
1917  a_dxyEtaMeanBiasTrend = MeanBiasTrendsDir.make<TH1F> ("means_dxy_eta",
1918  "#LT d_{xy} #GT vs #eta sector;#eta (sector);#LT d_{xy} #GT [#mum]",
1919  nBins_,lowedge,highedge);
1920 
1921  a_dxyEtaWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>("widths_dxy_eta",
1922  "#sigma_{d_{xy}} vs #eta sector;#eta (sector);#sigma_{d_{xy}} [#mum]",
1923  nBins_,lowedge,highedge);
1924 
1925  a_dzEtaMeanBiasTrend = MeanBiasTrendsDir.make<TH1F> ("means_dz_eta",
1926  "#LT d_{z} #GT vs #eta sector;#eta (sector);#LT d_{z} #GT [#mum]",
1927  nBins_,lowedge,highedge);
1928 
1929  a_dzEtaWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>("widths_dz_eta",
1930  "#sigma_{d_{xy}} vs #eta sector;#eta (sector);#sigma_{d_{z}} [#mum]",
1931  nBins_,lowedge,highedge);
1932 
1933  n_dxyPhiMeanBiasTrend = MeanBiasTrendsDir.make<TH1F> ("norm_means_dxy_phi",
1934  "#LT d_{xy}/#sigma_{d_{xy}} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{xy}/#sigma_{d_{xy}} #GT",
1935  nBins_,lowedge,highedge);
1936 
1937  n_dxyPhiWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>("norm_widths_dxy_phi",
1938  "width(d_{xy}/#sigma_{d_{xy}}) vs #phi sector;#varphi (sector) [degrees]; width(d_{xy}/#sigma_{d_{xy}})",
1939  nBins_,lowedge,highedge);
1940 
1941  n_dzPhiMeanBiasTrend = MeanBiasTrendsDir.make<TH1F> ("norm_means_dz_phi",
1942  "#LT d_{z}/#sigma_{d_{z}} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{z}/#sigma_{d_{z}} #GT",
1943  nBins_,lowedge,highedge);
1944 
1945  n_dzPhiWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>("norm_widths_dz_phi",
1946  "width(d_{z}/#sigma_{d_{z}}) vs #phi sector;#varphi (sector) [degrees];width(d_{z}/#sigma_{d_{z}})",
1947  nBins_,lowedge,highedge);
1948 
1949  n_dxyEtaMeanBiasTrend = MeanBiasTrendsDir.make<TH1F> ("norm_means_dxy_eta",
1950  "#LT d_{xy}/#sigma_{d_{xy}} #GT vs #eta sector;#eta (sector);#LT d_{xy}/#sigma_{d_{z}} #GT",
1951  nBins_,lowedge,highedge);
1952 
1953  n_dxyEtaWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>("norm_widths_dxy_eta",
1954  "width(d_{xy}/#sigma_{d_{xy}}) vs #eta sector;#eta (sector);width(d_{xy}/#sigma_{d_{z}})",
1955  nBins_,lowedge,highedge);
1956 
1957  n_dzEtaMeanBiasTrend = MeanBiasTrendsDir.make<TH1F> ("norm_means_dz_eta",
1958  "#LT d_{z}/#sigma_{d_{z}} #GT vs #eta sector;#eta (sector);#LT d_{z}/#sigma_{d_{z}} #GT",
1959  nBins_,lowedge,highedge);
1960 
1961  n_dzEtaWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>("norm_widths_dz_eta",
1962  "width(d_{z}/#sigma_{d_{z}}) vs #eta sector;#eta (sector);width(d_{z}/#sigma_{d_{z}})",
1963  nBins_,lowedge,highedge);
1964 
1965  // 2D maps
1966 
1967  a_dxyMeanBiasMap = Mean2DBiasMapsDir.make<TH2F> ("means_dxy_map",
1968  "#LT d_{xy} #GT map;#eta (sector);#varphi (sector) [degrees]",
1969  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1970 
1971  a_dzMeanBiasMap = Mean2DBiasMapsDir.make<TH2F> ("means_dz_map",
1972  "#LT d_{z} #GT map;#eta (sector);#varphi (sector) [degrees]",
1973  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1974 
1975  n_dxyMeanBiasMap = Mean2DBiasMapsDir.make<TH2F> ("norm_means_dxy_map",
1976  "#LT d_{xy}/#sigma_{d_{xy}} #GT map;#eta (sector);#varphi (sector) [degrees]",
1977  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1978 
1979  n_dzMeanBiasMap = Mean2DBiasMapsDir.make<TH2F> ("norm_means_dz_map",
1980  "#LT d_{z}/#sigma_{d_{z}} #GT map;#eta (sector);#varphi (sector) [degrees]",
1981  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1982 
1983  a_dxyWidthBiasMap = Width2DBiasMapsDir.make<TH2F> ("widths_dxy_map",
1984  "#sigma_{d_{xy}} map;#eta (sector);#varphi (sector) [degrees]",
1985  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1986 
1987  a_dzWidthBiasMap = Width2DBiasMapsDir.make<TH2F> ("widths_dz_map",
1988  "#sigma_{d_{z}} map;#eta (sector);#varphi (sector) [degrees]",
1989  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1990 
1991  n_dxyWidthBiasMap = Width2DBiasMapsDir.make<TH2F> ("norm_widths_dxy_map",
1992  "width(d_{xy}/#sigma_{d_{xy}}) map;#eta (sector);#varphi (sector) [degrees]",
1993  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1994 
1995  n_dzWidthBiasMap = Width2DBiasMapsDir.make<TH2F> ("norm_widths_dz_map",
1996  "width(d_{z}/#sigma_{d_{z}}) map;#eta (sector);#varphi (sector) [degrees]",
1997  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1998 
1999  // medians and MADs
2000 
2001  a_dxyPhiMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>("medians_dxy_phi",
2002  "Median of d_{xy} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{xy}) [#mum]",
2003  nBins_,lowedge,highedge);
2004 
2005  a_dxyPhiMADBiasTrend = MADBiasTrendsDir.make<TH1F> ("MADs_dxy_phi",
2006  "Median absolute deviation of d_{xy} vs #phi sector;#varphi (sector) [degrees];MAD(d_{xy}) [#mum]",
2007  nBins_,lowedge,highedge);
2008 
2009  a_dzPhiMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>("medians_dz_phi",
2010  "Median of d_{z} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{z}) [#mum]",
2011  nBins_,lowedge,highedge);
2012 
2013  a_dzPhiMADBiasTrend = MADBiasTrendsDir.make<TH1F> ("MADs_dz_phi",
2014  "Median absolute deviation of d_{z} vs #phi sector;#varphi (sector) [degrees];MAD(d_{z}) [#mum]",
2015  nBins_,lowedge,highedge);
2016 
2017  a_dxyEtaMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>("medians_dxy_eta",
2018  "Median of d_{xy} vs #eta sector;#eta (sector);#mu_{1/2}(d_{xy}) [#mum]",
2019  nBins_,lowedge,highedge);
2020 
2021  a_dxyEtaMADBiasTrend = MADBiasTrendsDir.make<TH1F> ("MADs_dxy_eta",
2022  "Median absolute deviation of d_{xy} vs #eta sector;#eta (sector);MAD(d_{xy}) [#mum]",
2023  nBins_,lowedge,highedge);
2024 
2025  a_dzEtaMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>("medians_dz_eta",
2026  "Median of d_{z} vs #eta sector;#eta (sector);#mu_{1/2}(d_{z}) [#mum]",
2027  nBins_,lowedge,highedge);
2028 
2029  a_dzEtaMADBiasTrend = MADBiasTrendsDir.make<TH1F> ("MADs_dz_eta",
2030  "Median absolute deviation of d_{z} vs #eta sector;#eta (sector);MAD(d_{z}) [#mum]",
2031  nBins_,lowedge,highedge);
2032 
2033  n_dxyPhiMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>("norm_medians_dxy_phi",
2034  "Median of d_{xy}/#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{xy}/#sigma_{d_{xy}})",
2035  nBins_,lowedge,highedge);
2036 
2037  n_dxyPhiMADBiasTrend = MADBiasTrendsDir.make<TH1F> ("norm_MADs_dxy_phi",
2038  "Median absolute deviation of d_{xy}/#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees]; MAD(d_{xy}/#sigma_{d_{xy}})",
2039  nBins_,lowedge,highedge);
2040 
2041  n_dzPhiMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>("norm_medians_dz_phi",
2042  "Median of d_{z}/#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{z}/#sigma_{d_{z}})",
2043  nBins_,lowedge,highedge);
2044 
2045  n_dzPhiMADBiasTrend = MADBiasTrendsDir.make<TH1F> ("norm_MADs_dz_phi",
2046  "Median absolute deviation of d_{z}/#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];MAD(d_{z}/#sigma_{d_{z}})",
2047  nBins_,lowedge,highedge);
2048 
2049  n_dxyEtaMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>("norm_medians_dxy_eta",
2050  "Median of d_{xy}/#sigma_{d_{xy}} vs #eta sector;#eta (sector);#mu_{1/2}(d_{xy}/#sigma_{d_{z}})",
2051  nBins_,lowedge,highedge);
2052 
2053  n_dxyEtaMADBiasTrend = MADBiasTrendsDir.make<TH1F> ("norm_MADs_dxy_eta",
2054  "Median absolute deviation of d_{xy}/#sigma_{d_{xy}} vs #eta sector;#eta (sector);MAD(d_{xy}/#sigma_{d_{z}})",
2055  nBins_,lowedge,highedge);
2056 
2057  n_dzEtaMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>("norm_medians_dz_eta",
2058  "Median of d_{z}/#sigma_{d_{z}} vs #eta sector;#eta (sector);#mu_{1/2}(d_{z}/#sigma_{d_{z}})",
2059  nBins_,lowedge,highedge);
2060 
2061  n_dzEtaMADBiasTrend = MADBiasTrendsDir.make<TH1F> ("norm_MADs_dz_eta",
2062  "Median absolute deviation of d_{z}/#sigma_{d_{z}} vs #eta sector;#eta (sector);MAD(d_{z}/#sigma_{d_{z}})",
2063  nBins_,lowedge,highedge);
2064 
2065  }
2066 }
2067 // ------------ method called once each job just after ending the event loop ------------
2068 void PrimaryVertexValidation::endJob()
2069 {
2070 
2071  edm::LogInfo("PrimaryVertexValidation")
2072  <<"######################################\n"
2073  <<"# PrimaryVertexValidation::endJob()\n"
2074  <<"# Number of analyzed events: "<<Nevt_<<"\n"
2075  <<"######################################";
2076 
2077  // means and widhts vs ladder and module number
2078 
2079  a_dxymodZMeanTrend = MeanTrendsDir.make<TH1F> ("means_dxy_modZ",
2080  "#LT d_{xy} #GT vs modZ;module number (Z);#LT d_{xy} #GT [#mum]",
2081  nModZ_,0.,nModZ_);
2082 
2083  a_dxymodZWidthTrend = WidthTrendsDir.make<TH1F>("widths_dxy_modZ",
2084  "#sigma_{d_{xy}} vs modZ;module number (Z);#sigma_{d_{xy}} [#mum]",
2085  nModZ_,0.,nModZ_);
2086 
2087  a_dzmodZMeanTrend = MeanTrendsDir.make<TH1F> ("means_dz_modZ",
2088  "#LT d_{z} #GT vs modZ;module number (Z);#LT d_{z} #GT [#mum]",
2089  nModZ_,0.,nModZ_);
2090 
2091  a_dzmodZWidthTrend = WidthTrendsDir.make<TH1F>("widths_dz_modZ",
2092  "#sigma_{d_{z}} vs modZ;module number (Z);#sigma_{d_{z}} [#mum]",
2093  nModZ_,0.,nModZ_);
2094 
2095  a_dxyladderMeanTrend = MeanTrendsDir.make<TH1F> ("means_dxy_ladder",
2096  "#LT d_{xy} #GT vs ladder;ladder number (#phi);#LT d_{xy} #GT [#mum]",
2097  nLadders_,0.,nLadders_);
2098 
2099  a_dxyladderWidthTrend = WidthTrendsDir.make<TH1F>("widths_dxy_ladder",
2100  "#sigma_{d_{xy}} vs ladder;ladder number (#phi);#sigma_{d_{xy}} [#mum]",
2101  nLadders_,0.,nLadders_);
2102 
2103  a_dzladderMeanTrend = MeanTrendsDir.make<TH1F> ("means_dz_ladder",
2104  "#LT d_{z} #GT vs ladder;ladder number (#phi);#LT d_{z} #GT [#mum]"
2105  ,nLadders_,0.,nLadders_);
2106 
2107  a_dzladderWidthTrend = WidthTrendsDir.make<TH1F>("widths_dz_ladder",
2108  "#sigma_{d_{z}} vs ladder;ladder number (#phi);#sigma_{d_{z}} [#mum]",
2109  nLadders_,0.,nLadders_);
2110 
2111  n_dxymodZMeanTrend = MeanTrendsDir.make<TH1F> ("norm_means_dxy_modZ",
2112  "#LT d_{xy}/#sigma_{d_{xy}} #GT vs modZ;module number (Z);#LT d_{xy}/#sigma_{d_{xy}} #GT",
2113  nModZ_,0.,nModZ_);
2114 
2115  n_dxymodZWidthTrend = WidthTrendsDir.make<TH1F>("norm_widths_dxy_modZ",
2116  "width(d_{xy}/#sigma_{d_{xy}}) vs modZ;module number (Z); width(d_{xy}/#sigma_{d_{xy}})",
2117  nModZ_,0.,nModZ_);
2118 
2119  n_dzmodZMeanTrend = MeanTrendsDir.make<TH1F> ("norm_means_dz_modZ",
2120  "#LT d_{z}/#sigma_{d_{z}} #GT vs modZ;module number (Z);#LT d_{z}/#sigma_{d_{z}} #GT",
2121  nModZ_,0.,nModZ_);
2122 
2123  n_dzmodZWidthTrend = WidthTrendsDir.make<TH1F>("norm_widths_dz_modZ",
2124  "width(d_{z}/#sigma_{d_{z}}) vs pT;module number (Z);width(d_{z}/#sigma_{d_{z}})",
2125  nModZ_,0.,nModZ_);
2126 
2127  n_dxyladderMeanTrend = MeanTrendsDir.make<TH1F> ("norm_means_dxy_ladder",
2128  "#LT d_{xy}/#sigma_{d_{xy}} #GT vs ladder;ladder number (#phi);#LT d_{xy}/#sigma_{d_{z}} #GT",
2129  nLadders_,0.,nLadders_);
2130 
2131  n_dxyladderWidthTrend = WidthTrendsDir.make<TH1F>("norm_widths_dxy_ladder",
2132  "width(d_{xy}/#sigma_{d_{xy}}) vs ladder;ladder number (#phi);width(d_{xy}/#sigma_{d_{z}})",
2133  nLadders_,0.,nLadders_);
2134 
2135  n_dzladderMeanTrend = MeanTrendsDir.make<TH1F> ("norm_means_dz_ladder",
2136  "#LT d_{z}/#sigma_{d_{z}} #GT vs ladder;ladder number (#phi);#LT d_{z}/#sigma_{d_{z}} #GT",
2137  nLadders_,0.,nLadders_);
2138 
2139  n_dzladderWidthTrend = WidthTrendsDir.make<TH1F>("norm_widths_dz_ladder",
2140  "width(d_{z}/#sigma_{d_{z}}) vs ladder;ladder number (#phi);width(d_{z}/#sigma_{d_{z}})",
2141  nLadders_,0.,nLadders_);
2142 
2143  if(useTracksFromRecoVtx_){
2144 
2145  fillTrendPlotByIndex(a_dxyPhiMeanBiasTrend ,a_dxyPhiBiasResiduals,PVValHelper::MEAN,PVValHelper::phi);
2146  fillTrendPlotByIndex(a_dxyPhiWidthBiasTrend,a_dxyPhiBiasResiduals,PVValHelper::WIDTH,PVValHelper::phi);
2147  fillTrendPlotByIndex(a_dzPhiMeanBiasTrend ,a_dzPhiBiasResiduals ,PVValHelper::MEAN,PVValHelper::phi);
2148  fillTrendPlotByIndex(a_dzPhiWidthBiasTrend ,a_dzPhiBiasResiduals ,PVValHelper::WIDTH,PVValHelper::phi);
2149 
2150  fillTrendPlotByIndex(a_dxyEtaMeanBiasTrend ,a_dxyEtaBiasResiduals,PVValHelper::MEAN,PVValHelper::eta);
2151  fillTrendPlotByIndex(a_dxyEtaWidthBiasTrend,a_dxyEtaBiasResiduals,PVValHelper::WIDTH,PVValHelper::eta);
2152  fillTrendPlotByIndex(a_dzEtaMeanBiasTrend ,a_dzEtaBiasResiduals ,PVValHelper::MEAN,PVValHelper::eta);
2153  fillTrendPlotByIndex(a_dzEtaWidthBiasTrend ,a_dzEtaBiasResiduals ,PVValHelper::WIDTH,PVValHelper::eta);
2154 
2155  fillTrendPlotByIndex(n_dxyPhiMeanBiasTrend ,n_dxyPhiBiasResiduals,PVValHelper::MEAN,PVValHelper::phi);
2156  fillTrendPlotByIndex(n_dxyPhiWidthBiasTrend,n_dxyPhiBiasResiduals,PVValHelper::WIDTH,PVValHelper::phi);
2157  fillTrendPlotByIndex(n_dzPhiMeanBiasTrend ,n_dzPhiBiasResiduals ,PVValHelper::MEAN,PVValHelper::phi);
2158  fillTrendPlotByIndex(n_dzPhiWidthBiasTrend ,n_dzPhiBiasResiduals ,PVValHelper::WIDTH,PVValHelper::phi);
2159 
2160  fillTrendPlotByIndex(n_dxyEtaMeanBiasTrend ,n_dxyEtaBiasResiduals,PVValHelper::MEAN,PVValHelper::eta);
2161  fillTrendPlotByIndex(n_dxyEtaWidthBiasTrend,n_dxyEtaBiasResiduals,PVValHelper::WIDTH,PVValHelper::eta);
2162  fillTrendPlotByIndex(n_dzEtaMeanBiasTrend ,n_dzEtaBiasResiduals ,PVValHelper::MEAN,PVValHelper::eta);
2163  fillTrendPlotByIndex(n_dzEtaWidthBiasTrend ,n_dzEtaBiasResiduals ,PVValHelper::WIDTH,PVValHelper::eta);
2164 
2165  // medians and MADs
2166 
2167  fillTrendPlotByIndex(a_dxyPhiMedianBiasTrend,a_dxyPhiBiasResiduals,PVValHelper::MEDIAN,PVValHelper::phi);
2168  fillTrendPlotByIndex(a_dxyPhiMADBiasTrend ,a_dxyPhiBiasResiduals,PVValHelper::MAD,PVValHelper::phi);
2169  fillTrendPlotByIndex(a_dzPhiMedianBiasTrend ,a_dzPhiBiasResiduals ,PVValHelper::MEDIAN,PVValHelper::phi);
2170  fillTrendPlotByIndex(a_dzPhiMADBiasTrend ,a_dzPhiBiasResiduals ,PVValHelper::MAD,PVValHelper::phi);
2171 
2172  fillTrendPlotByIndex(a_dxyEtaMedianBiasTrend,a_dxyEtaBiasResiduals,PVValHelper::MEDIAN,PVValHelper::eta);
2173  fillTrendPlotByIndex(a_dxyEtaMADBiasTrend ,a_dxyEtaBiasResiduals,PVValHelper::MAD,PVValHelper::eta);
2174  fillTrendPlotByIndex(a_dzEtaMedianBiasTrend ,a_dzEtaBiasResiduals ,PVValHelper::MEDIAN,PVValHelper::eta);
2175  fillTrendPlotByIndex(a_dzEtaMADBiasTrend ,a_dzEtaBiasResiduals ,PVValHelper::MAD,PVValHelper::eta);
2176 
2177  fillTrendPlotByIndex(n_dxyPhiMedianBiasTrend,n_dxyPhiBiasResiduals,PVValHelper::MEDIAN,PVValHelper::phi);
2178  fillTrendPlotByIndex(n_dxyPhiMADBiasTrend ,n_dxyPhiBiasResiduals,PVValHelper::MAD,PVValHelper::phi);
2179  fillTrendPlotByIndex(n_dzPhiMedianBiasTrend ,n_dzPhiBiasResiduals ,PVValHelper::MEDIAN,PVValHelper::phi);
2180  fillTrendPlotByIndex(n_dzPhiMADBiasTrend ,n_dzPhiBiasResiduals ,PVValHelper::MAD,PVValHelper::phi);
2181 
2182  fillTrendPlotByIndex(n_dxyEtaMedianBiasTrend,n_dxyEtaBiasResiduals,PVValHelper::MEDIAN,PVValHelper::eta);
2183  fillTrendPlotByIndex(n_dxyEtaMADBiasTrend ,n_dxyEtaBiasResiduals,PVValHelper::MAD,PVValHelper::eta);
2184  fillTrendPlotByIndex(n_dzEtaMedianBiasTrend ,n_dzEtaBiasResiduals ,PVValHelper::MEDIAN,PVValHelper::eta);
2185  fillTrendPlotByIndex(n_dzEtaMADBiasTrend ,n_dzEtaBiasResiduals ,PVValHelper::MAD,PVValHelper::eta);
2186 
2187  // 2d Maps
2188 
2189  fillMap(a_dxyMeanBiasMap ,a_dxyBiasResidualsMap,PVValHelper::MEAN);
2190  fillMap(a_dxyWidthBiasMap,a_dxyBiasResidualsMap,PVValHelper::WIDTH);
2191  fillMap(a_dzMeanBiasMap ,a_dzBiasResidualsMap,PVValHelper::MEAN);
2192  fillMap(a_dzWidthBiasMap ,a_dzBiasResidualsMap,PVValHelper::WIDTH);
2193 
2194  fillMap(n_dxyMeanBiasMap ,n_dxyBiasResidualsMap,PVValHelper::MEAN);
2195  fillMap(n_dxyWidthBiasMap,n_dxyBiasResidualsMap,PVValHelper::WIDTH);
2196  fillMap(n_dzMeanBiasMap ,n_dzBiasResidualsMap,PVValHelper::MEAN);
2197  fillMap(n_dzWidthBiasMap ,n_dzBiasResidualsMap,PVValHelper::WIDTH);
2198 
2199  }
2200 
2201  // do profiles
2202 
2203  fillTrendPlotByIndex(a_dxyPhiMeanTrend, a_dxyPhiResiduals,PVValHelper::MEAN,PVValHelper::phi);
2204  fillTrendPlotByIndex(a_dxyPhiWidthTrend,a_dxyPhiResiduals,PVValHelper::WIDTH,PVValHelper::phi);
2205  fillTrendPlotByIndex(a_dzPhiMeanTrend ,a_dzPhiResiduals ,PVValHelper::MEAN,PVValHelper::phi);
2206  fillTrendPlotByIndex(a_dzPhiWidthTrend ,a_dzPhiResiduals ,PVValHelper::WIDTH,PVValHelper::phi);
2207 
2208  fillTrendPlotByIndex(a_dxyEtaMeanTrend ,a_dxyEtaResiduals,PVValHelper::MEAN,PVValHelper::eta);
2209  fillTrendPlotByIndex(a_dxyEtaWidthTrend,a_dxyEtaResiduals,PVValHelper::WIDTH,PVValHelper::eta);
2210  fillTrendPlotByIndex(a_dzEtaMeanTrend ,a_dzEtaResiduals ,PVValHelper::MEAN,PVValHelper::eta);
2211  fillTrendPlotByIndex(a_dzEtaWidthTrend ,a_dzEtaResiduals ,PVValHelper::WIDTH,PVValHelper::eta);
2212 
2213  fillTrendPlotByIndex(n_dxyPhiMeanTrend ,n_dxyPhiResiduals,PVValHelper::MEAN,PVValHelper::phi);
2214  fillTrendPlotByIndex(n_dxyPhiWidthTrend,n_dxyPhiResiduals,PVValHelper::WIDTH,PVValHelper::phi);
2215  fillTrendPlotByIndex(n_dzPhiMeanTrend ,n_dzPhiResiduals ,PVValHelper::MEAN,PVValHelper::phi);
2216  fillTrendPlotByIndex(n_dzPhiWidthTrend ,n_dzPhiResiduals ,PVValHelper::WIDTH);
2217 
2218  fillTrendPlotByIndex(n_dxyEtaMeanTrend ,n_dxyEtaResiduals,PVValHelper::MEAN,PVValHelper::eta);
2219  fillTrendPlotByIndex(n_dxyEtaWidthTrend,n_dxyEtaResiduals,PVValHelper::WIDTH,PVValHelper::eta);
2220  fillTrendPlotByIndex(n_dzEtaMeanTrend ,n_dzEtaResiduals ,PVValHelper::MEAN,PVValHelper::eta);
2221  fillTrendPlotByIndex(n_dzEtaWidthTrend ,n_dzEtaResiduals ,PVValHelper::WIDTH,PVValHelper::eta);
2222 
2223  // vs transverse momentum
2224 
2225  fillTrendPlotByIndex(a_dxypTMeanTrend ,h_dxy_pT_,PVValHelper::MEAN );
2226  fillTrendPlotByIndex(a_dxypTWidthTrend,h_dxy_pT_,PVValHelper::WIDTH);
2227  fillTrendPlotByIndex(a_dzpTMeanTrend ,h_dz_pT_ ,PVValHelper::MEAN );
2228  fillTrendPlotByIndex(a_dzpTWidthTrend ,h_dz_pT_ ,PVValHelper::WIDTH);
2229 
2230  fillTrendPlotByIndex(a_dxypTCentralMeanTrend ,h_dxy_Central_pT_,PVValHelper::MEAN );
2231  fillTrendPlotByIndex(a_dxypTCentralWidthTrend,h_dxy_Central_pT_,PVValHelper::WIDTH);
2232  fillTrendPlotByIndex(a_dzpTCentralMeanTrend ,h_dz_Central_pT_ ,PVValHelper::MEAN );
2233  fillTrendPlotByIndex(a_dzpTCentralWidthTrend ,h_dz_Central_pT_ ,PVValHelper::WIDTH);
2234 
2235  fillTrendPlotByIndex(n_dxypTMeanTrend ,h_norm_dxy_pT_,PVValHelper::MEAN );
2236  fillTrendPlotByIndex(n_dxypTWidthTrend,h_norm_dxy_pT_,PVValHelper::WIDTH);
2237  fillTrendPlotByIndex(n_dzpTMeanTrend ,h_norm_dz_pT_ ,PVValHelper::MEAN );
2238  fillTrendPlotByIndex(n_dzpTWidthTrend ,h_norm_dz_pT_ ,PVValHelper::WIDTH);
2239 
2240  fillTrendPlotByIndex(n_dxypTCentralMeanTrend ,h_norm_dxy_Central_pT_,PVValHelper::MEAN );
2241  fillTrendPlotByIndex(n_dxypTCentralWidthTrend,h_norm_dxy_Central_pT_,PVValHelper::WIDTH);
2242  fillTrendPlotByIndex(n_dzpTCentralMeanTrend ,h_norm_dz_Central_pT_ ,PVValHelper::MEAN );
2243  fillTrendPlotByIndex(n_dzpTCentralWidthTrend ,h_norm_dz_Central_pT_ ,PVValHelper::WIDTH);
2244 
2245  // vs ladder and module number
2246 
2247  fillTrendPlotByIndex(a_dxymodZMeanTrend ,h_dxy_modZ_,PVValHelper::MEAN);
2248  fillTrendPlotByIndex(a_dxymodZWidthTrend ,h_dxy_modZ_,PVValHelper::WIDTH);
2249  fillTrendPlotByIndex(a_dzmodZMeanTrend ,h_dz_modZ_,PVValHelper::MEAN);
2250  fillTrendPlotByIndex(a_dzmodZWidthTrend ,h_dz_modZ_,PVValHelper::WIDTH);
2251 
2252  fillTrendPlotByIndex(a_dxyladderMeanTrend ,h_dxy_ladder_,PVValHelper::MEAN);
2253  fillTrendPlotByIndex(a_dxyladderWidthTrend,h_dxy_ladder_,PVValHelper::WIDTH);
2254  fillTrendPlotByIndex(a_dzladderMeanTrend ,h_dz_ladder_,PVValHelper::MEAN);
2255  fillTrendPlotByIndex(a_dzladderWidthTrend ,h_dz_ladder_,PVValHelper::WIDTH);
2256 
2257  fillTrendPlotByIndex(n_dxymodZMeanTrend ,h_norm_dxy_modZ_,PVValHelper::MEAN);
2258  fillTrendPlotByIndex(n_dxymodZWidthTrend ,h_norm_dxy_modZ_,PVValHelper::WIDTH);
2259  fillTrendPlotByIndex(n_dzmodZMeanTrend ,h_norm_dz_modZ_,PVValHelper::MEAN);
2260  fillTrendPlotByIndex(n_dzmodZWidthTrend ,h_norm_dz_modZ_,PVValHelper::WIDTH);
2261 
2262  fillTrendPlotByIndex(n_dxyladderMeanTrend ,h_norm_dxy_ladder_,PVValHelper::MEAN);
2263  fillTrendPlotByIndex(n_dxyladderWidthTrend,h_norm_dxy_ladder_,PVValHelper::WIDTH);
2264  fillTrendPlotByIndex(n_dzladderMeanTrend ,h_norm_dz_ladder_,PVValHelper::MEAN);
2265  fillTrendPlotByIndex(n_dzladderWidthTrend ,h_norm_dz_ladder_,PVValHelper::WIDTH);
2266 
2267  // medians and MADs
2268 
2269  fillTrendPlotByIndex(a_dxyPhiMedianTrend,a_dxyPhiResiduals,PVValHelper::MEDIAN,PVValHelper::phi);
2270  fillTrendPlotByIndex(a_dxyPhiMADTrend ,a_dxyPhiResiduals,PVValHelper::MAD,PVValHelper::phi);
2271 
2272  fillTrendPlotByIndex(a_dzPhiMedianTrend ,a_dzPhiResiduals ,PVValHelper::MEDIAN,PVValHelper::phi);
2273  fillTrendPlotByIndex(a_dzPhiMADTrend ,a_dzPhiResiduals ,PVValHelper::MAD,PVValHelper::phi);
2274 
2275  fillTrendPlotByIndex(a_dxyEtaMedianTrend,a_dxyEtaResiduals,PVValHelper::MEDIAN,PVValHelper::eta);
2276  fillTrendPlotByIndex(a_dxyEtaMADTrend ,a_dxyEtaResiduals,PVValHelper::MAD,PVValHelper::eta);
2277  fillTrendPlotByIndex(a_dzEtaMedianTrend ,a_dzEtaResiduals ,PVValHelper::MEDIAN,PVValHelper::eta);
2278  fillTrendPlotByIndex(a_dzEtaMADTrend ,a_dzEtaResiduals ,PVValHelper::MAD,PVValHelper::eta);
2279 
2280  fillTrendPlotByIndex(n_dxyPhiMedianTrend,n_dxyPhiResiduals,PVValHelper::MEDIAN,PVValHelper::phi);
2281  fillTrendPlotByIndex(n_dxyPhiMADTrend ,n_dxyPhiResiduals,PVValHelper::MAD,PVValHelper::phi);
2282  fillTrendPlotByIndex(n_dzPhiMedianTrend ,n_dzPhiResiduals ,PVValHelper::MEDIAN,PVValHelper::phi);
2283  fillTrendPlotByIndex(n_dzPhiMADTrend ,n_dzPhiResiduals ,PVValHelper::MAD,PVValHelper::phi);
2284 
2285  fillTrendPlotByIndex(n_dxyEtaMedianTrend,n_dxyEtaResiduals,PVValHelper::MEDIAN,PVValHelper::eta);
2286  fillTrendPlotByIndex(n_dxyEtaMADTrend ,n_dxyEtaResiduals,PVValHelper::MAD,PVValHelper::eta);
2287  fillTrendPlotByIndex(n_dzEtaMedianTrend ,n_dzEtaResiduals ,PVValHelper::MEDIAN,PVValHelper::eta);
2288  fillTrendPlotByIndex(n_dzEtaMADTrend ,n_dzEtaResiduals ,PVValHelper::MAD,PVValHelper::eta);
2289 
2290  // 2D Maps
2291 
2292  fillMap(a_dxyMeanMap ,a_dxyResidualsMap,PVValHelper::MEAN);
2293  fillMap(a_dxyWidthMap,a_dxyResidualsMap,PVValHelper::WIDTH);
2294  fillMap(a_dzMeanMap ,a_dzResidualsMap,PVValHelper::MEAN);
2295  fillMap(a_dzWidthMap ,a_dzResidualsMap,PVValHelper::WIDTH);
2296 
2297  fillMap(n_dxyMeanMap ,n_dxyResidualsMap,PVValHelper::MEAN);
2298  fillMap(n_dxyWidthMap,n_dxyResidualsMap,PVValHelper::WIDTH);
2299  fillMap(n_dzMeanMap ,n_dzResidualsMap,PVValHelper::MEAN);
2300  fillMap(n_dzWidthMap ,n_dzResidualsMap,PVValHelper::WIDTH);
2301 
2302 }
2303 
2304 //*************************************************************
2305 bool PrimaryVertexValidation::isBFieldConsistentWithMode(const edm::EventSetup& iSetup) const
2306 //*************************************************************
2307 {
2308  edm::ESHandle<RunInfo> runInfo;
2309  iSetup.get<RunInfoRcd>().get(runInfo);
2310 
2311  double average_current = runInfo.product()->m_avg_current;
2312  bool isOn = (average_current > 2000.);
2313  bool is0T = (ptOfProbe_==0.);
2314 
2315  return ( (isOn && !is0T) || (!isOn && is0T) );
2316 }
2317 
2318 //*************************************************************
2319 void PrimaryVertexValidation::SetVarToZero()
2320 //*************************************************************
2321 {
2322  nTracks_ = 0;
2323  nClus_ = 0;
2324  nOfflineVertices_=0;
2325  RunNumber_ =0;
2326  LuminosityBlockNumber_=0;
2327  xOfflineVertex_ =-999.;
2328  yOfflineVertex_ =-999.;
2329  zOfflineVertex_ =-999.;
2330  xErrOfflineVertex_=0.;
2331  yErrOfflineVertex_=0.;
2332  zErrOfflineVertex_=0.;
2333  BSx0_ = -999.;
2334  BSy0_ = -999.;
2335  BSz0_ = -999.;
2336  Beamsigmaz_=-999.;
2337  Beamdxdz_=-999.;
2338  BeamWidthX_=-999.;
2339  BeamWidthY_=-999.;
2340  wxy2_=-999.;
2341 
2342  for ( int i=0; i<nMaxtracks_; ++i ) {
2343 
2344  pt_[i] = 0;
2345  p_[i] = 0;
2346  nhits_[i] = 0;
2347  nhits1D_[i] = 0;
2348  nhits2D_[i] = 0;
2349  nhitsBPIX_[i] = 0;
2350  nhitsFPIX_[i] = 0;
2351  nhitsTIB_[i] = 0;
2352  nhitsTID_[i] = 0;
2353  nhitsTOB_[i] = 0;
2354  nhitsTEC_[i] = 0;
2355  isHighPurity_[i] = 0;
2356  eta_[i] = 0;
2357  theta_[i] = 0;
2358  phi_[i] = 0;
2359  chi2_[i] = 0;
2360  chi2ndof_[i] = 0;
2361  charge_[i] = 0;
2362  qoverp_[i] = 0;
2363  dz_[i] = 0;
2364  dxy_[i] = 0;
2365  dzBs_[i] = 0;
2366  dxyBs_[i] = 0;
2367  xPCA_[i] = 0;
2368  yPCA_[i] = 0;
2369  zPCA_[i] = 0;
2370  xUnbiasedVertex_[i] =0;
2371  yUnbiasedVertex_[i] =0;
2372  zUnbiasedVertex_[i] =0;
2373  chi2normUnbiasedVertex_[i]=0;
2374  chi2UnbiasedVertex_[i]=0;
2375  chi2ProbUnbiasedVertex_[i]=0;
2376  DOFUnbiasedVertex_[i]=0;
2377  sumOfWeightsUnbiasedVertex_[i]=0;
2378  tracksUsedForVertexing_[i]=0;
2379  dxyFromMyVertex_[i]=0;
2380  dzFromMyVertex_[i]=0;
2381  d3DFromMyVertex_[i]=0;
2382  dxyErrorFromMyVertex_[i]=0;
2383  dzErrorFromMyVertex_[i]=0;
2384  d3DErrorFromMyVertex_[i]=0;
2385  IPTsigFromMyVertex_[i]=0;
2386  IPLsigFromMyVertex_[i]=0;
2387  IP3DsigFromMyVertex_[i]=0;
2388  hasRecVertex_[i] = 0;
2389  isGoodTrack_[i] = 0;
2390  }
2391 }
2392 
2393 //*************************************************************
2394 void PrimaryVertexValidation::fillTrendPlot(TH1F* trendPlot, TH1F* residualsPlot[100], PVValHelper::estimator fitPar_, const std::string& var_)
2395 //*************************************************************
2396 {
2397 
2398  for ( int i=0; i<nBins_; ++i ) {
2399 
2400  char phibincenter[129];
2401  auto phiBins = theDetails_.trendbins[PVValHelper::phi];
2402  sprintf(phibincenter,"%.f",(phiBins[i]+phiBins[i+1])/2.);
2403 
2404  char etabincenter[129];
2405  auto etaBins = theDetails_.trendbins[PVValHelper::eta];
2406  sprintf(etabincenter,"%.1f",(etaBins[i]+etaBins[i+1])/2.);
2407 
2408  switch(fitPar_)
2409  {
2410  case PVValHelper::MEAN:
2411  {
2412  float mean_ = PVValHelper::fitResiduals(residualsPlot[i]).first.value();
2413  float meanErr_ = PVValHelper::fitResiduals(residualsPlot[i]).first.error();
2414  trendPlot->SetBinContent(i+1,mean_);
2415  trendPlot->SetBinError(i+1,meanErr_);
2416  break;
2417  }
2418  case PVValHelper::WIDTH:
2419  {
2420  float width_ = PVValHelper::fitResiduals(residualsPlot[i]).second.value();
2421  float widthErr_ = PVValHelper::fitResiduals(residualsPlot[i]).second.error();
2422  trendPlot->SetBinContent(i+1,width_);
2423  trendPlot->SetBinError(i+1,widthErr_);
2424  break;
2425  }
2426  case PVValHelper::MEDIAN:
2427  {
2428  float median_ = PVValHelper::getMedian(residualsPlot[i]).value();
2429  float medianErr_ = PVValHelper::getMedian(residualsPlot[i]).error();
2430  trendPlot->SetBinContent(i+1,median_);
2431  trendPlot->SetBinError(i+1,medianErr_);
2432  break;
2433  }
2434  case PVValHelper::MAD:
2435  {
2436  float mad_ = PVValHelper::getMAD(residualsPlot[i]).value();
2437  float madErr_ = PVValHelper::getMAD(residualsPlot[i]).error();
2438  trendPlot->SetBinContent(i+1,mad_);
2439  trendPlot->SetBinError(i+1,madErr_);
2440  break;
2441  }
2442  default:
2443  edm::LogWarning("PrimaryVertexValidation")<<"fillTrendPlot() "<<fitPar_<<" unknown estimator!"<<std::endl;
2444  break;
2445  }
2446 
2447  if(var_.find("eta") != std::string::npos){
2448  trendPlot->GetXaxis()->SetBinLabel(i+1,etabincenter);
2449  } else if(var_.find("phi") != std::string::npos){
2450  trendPlot->GetXaxis()->SetBinLabel(i+1,phibincenter);
2451  } else {
2452  edm::LogWarning("PrimaryVertexValidation")<<"fillTrendPlot() "<<var_<<" unknown track parameter!"<<std::endl;
2453  }
2454  }
2455 }
2456 
2457 //*************************************************************
2458 void PrimaryVertexValidation::fillTrendPlotByIndex(TH1F* trendPlot,std::vector<TH1F*>& h, PVValHelper::estimator fitPar_, PVValHelper::plotVariable plotVar)
2459 //*************************************************************
2460 {
2461 
2462  for(auto iterator = h.begin(); iterator != h.end(); iterator++) {
2463 
2464  unsigned int bin = std::distance(h.begin(),iterator)+1;
2465  std::pair<Measurement1D, Measurement1D> myFit = PVValHelper::fitResiduals((*iterator));
2466 
2467  switch(fitPar_)
2468  {
2469  case PVValHelper::MEAN:
2470  {
2471  float mean_ = myFit.first.value();
2472  float meanErr_ = myFit.first.error();
2473  trendPlot->SetBinContent(bin,mean_);
2474  trendPlot->SetBinError(bin,meanErr_);
2475  break;
2476  }
2477  case PVValHelper::WIDTH:
2478  {
2479  float width_ = myFit.second.value();
2480  float widthErr_ = myFit.second.error();
2481  trendPlot->SetBinContent(bin,width_);
2482  trendPlot->SetBinError(bin,widthErr_);
2483  break;
2484  }
2485  case PVValHelper::MEDIAN:
2486  {
2487  float median_ = PVValHelper::getMedian(*iterator).value();
2488  float medianErr_ = PVValHelper::getMedian(*iterator).error();
2489  trendPlot->SetBinContent(bin,median_);
2490  trendPlot->SetBinError(bin,medianErr_);
2491  break;
2492  }
2493  case PVValHelper::MAD:
2494  {
2495  float mad_ = PVValHelper::getMAD(*iterator).value();
2496  float madErr_ = PVValHelper::getMAD(*iterator).error();
2497  trendPlot->SetBinContent(bin,mad_);
2498  trendPlot->SetBinError(bin,madErr_);
2499  break;
2500  }
2501  default:
2502  edm::LogWarning("PrimaryVertexValidation")<<"fillTrendPlotByIndex() "<<fitPar_<<" unknown estimator!"<<std::endl;
2503  break;
2504  }
2505 
2506  char bincenter[129];
2507  if(plotVar == PVValHelper::eta){
2508  auto etaBins = theDetails_.trendbins[PVValHelper::eta];
2509  sprintf(bincenter,"%.1f",(etaBins[bin-1]+etaBins[bin])/2.);
2510  trendPlot->GetXaxis()->SetBinLabel(bin,bincenter);
2511  } else if(plotVar == PVValHelper::phi){
2512  auto phiBins = theDetails_.trendbins[PVValHelper::phi];
2513  sprintf(bincenter,"%.f",(phiBins[bin-1]+phiBins[bin])/2.);
2514  trendPlot->GetXaxis()->SetBinLabel(bin,bincenter);
2515  } else {
2517  //edm::LogWarning("PrimaryVertexValidation")<<"fillTrendPlotByIndex() "<< plotVar <<" unknown track parameter!"<<std::endl;
2518  }
2519 
2520  }
2521 }
2522 
2523 //*************************************************************
2524 void PrimaryVertexValidation::fillMap(TH2F* trendMap, TH1F* residualsMapPlot[100][100], PVValHelper::estimator fitPar_)
2525 //*************************************************************
2526 {
2527 
2528  for ( int i=0; i<nBins_; ++i ) {
2529 
2530  char phibincenter[129];
2531  auto phiBins = theDetails_.trendbins[PVValHelper::phi];
2532  sprintf(phibincenter,"%.f",(phiBins[i]+phiBins[i+1])/2.);
2533 
2534  trendMap->GetYaxis()->SetBinLabel(i+1,phibincenter);
2535 
2536  for ( int j=0; j<nBins_; ++j ) {
2537 
2538  char etabincenter[129];
2539  auto etaBins = theDetails_.trendbins[PVValHelper::eta];
2540  sprintf(etabincenter,"%.1f",(etaBins[j]+etaBins[j+1])/2.);
2541 
2542  if(i==0) { trendMap->GetXaxis()->SetBinLabel(j+1,etabincenter); }
2543 
2544  switch (fitPar_)
2545  {
2546  case PVValHelper::MEAN:
2547  {
2548  float mean_ = PVValHelper::fitResiduals(residualsMapPlot[i][j]).first.value();
2549  float meanErr_ = PVValHelper::fitResiduals(residualsMapPlot[i][j]).first.error();
2550  trendMap->SetBinContent(j+1,i+1,mean_);
2551  trendMap->SetBinError(j+1,i+1,meanErr_);
2552  break;
2553  }
2554  case PVValHelper::WIDTH:
2555  {
2556  float width_ = PVValHelper::fitResiduals(residualsMapPlot[i][j]).second.value();
2557  float widthErr_ = PVValHelper::fitResiduals(residualsMapPlot[i][j]).second.error();
2558  trendMap->SetBinContent(j+1,i+1,width_);
2559  trendMap->SetBinError(j+1,i+1,widthErr_);
2560  break;
2561  }
2562  case PVValHelper::MEDIAN:
2563  {
2564  float median_ = PVValHelper::getMedian(residualsMapPlot[i][j]).value();
2565  float medianErr_ = PVValHelper::getMedian(residualsMapPlot[i][j]).error();
2566  trendMap->SetBinContent(j+1,i+1,median_);
2567  trendMap->SetBinError(j+1,i+1,medianErr_);
2568  break;
2569  }
2570  case PVValHelper::MAD:
2571  {
2572  float mad_ = PVValHelper::getMAD(residualsMapPlot[i][j]).value();
2573  float madErr_ = PVValHelper::getMAD(residualsMapPlot[i][j]).error();
2574  trendMap->SetBinContent(j+1,i+1,mad_);
2575  trendMap->SetBinError(j+1,i+1,madErr_);
2576  break;
2577  }
2578  default:
2579  edm::LogWarning("PrimaryVertexValidation:") <<" fillMap() "<<fitPar_<<" unknown estimator!"<<std::endl;
2580  }
2581  } // closes loop on eta bins
2582  } // cloeses loop on phi bins
2583 }
2584 
2585 //*************************************************************
2586 bool PrimaryVertexValidation::vtxSort( const reco::Vertex & a, const reco::Vertex & b )
2587 //*************************************************************
2588 {
2589  if( a.tracksSize() != b.tracksSize() )
2590  return a.tracksSize() > b.tracksSize() ? true : false ;
2591  else
2592  return a.chi2() < b.chi2() ? true : false ;
2593 }
2594 
2595 //*************************************************************
2596 bool PrimaryVertexValidation::passesTrackCuts(const reco::Track & track, const reco::Vertex & vertex,const std::string& qualityString_, double dxyErrMax_,double dzErrMax_, double ptErrMax_)
2597 //*************************************************************
2598 {
2599 
2600  math::XYZPoint vtxPoint(0.0,0.0,0.0);
2601  double vzErr =0.0, vxErr=0.0, vyErr=0.0;
2602  vtxPoint=vertex.position();
2603  vzErr=vertex.zError();
2604  vxErr=vertex.xError();
2605  vyErr=vertex.yError();
2606 
2607  double dxy=0.0, dz=0.0, dxysigma=0.0, dzsigma=0.0;
2608  dxy = track.dxy(vtxPoint);
2609  dz = track.dz(vtxPoint);
2610  dxysigma = sqrt(track.d0Error()*track.d0Error()+vxErr*vyErr);
2611  dzsigma = sqrt(track.dzError()*track.dzError()+vzErr*vzErr);
2612 
2613  if(track.quality(reco::TrackBase::qualityByName(qualityString_)) != 1)return false;
2614  if(std::abs(dxy/dxysigma) > dxyErrMax_) return false;
2615  if(std::abs(dz/dzsigma) > dzErrMax_) return false;
2616  if(track.ptError() / track.pt() > ptErrMax_) return false;
2617 
2618  return true;
2619 }
2620 
2621 
2622 //*************************************************************
2623 std::map<std::string, TH1*> PrimaryVertexValidation::bookVertexHistograms(const TFileDirectory& dir)
2624 //*************************************************************
2625 {
2626 
2627  TH1F::SetDefaultSumw2(kTRUE);
2628 
2629  std::map<std::string, TH1*> h;
2630 
2631  // histograms of track quality (Data and MC)
2632  std::string types[] = {"all","sel"};
2633  for(const auto & type : types){
2634  h["pseudorapidity_"+type] =dir.make <TH1F>(("rapidity_"+type).c_str(),"track pseudorapidity; track #eta; tracks",100,-3., 3.);
2635  h["z0_"+type] = dir.make<TH1F>(("z0_"+type).c_str(),"track z_{0};track z_{0} (cm);tracks",80,-40., 40.);
2636  h["phi_"+type] = dir.make<TH1F>(("phi_"+type).c_str(),"track #phi; track #phi;tracks",80,-M_PI, M_PI);
2637  h["eta_"+type] = dir.make<TH1F>(("eta_"+type).c_str(),"track #eta; track #eta;tracks",80,-4., 4.);
2638  h["pt_"+type] = dir.make<TH1F>(("pt_"+type).c_str(),"track p_{T}; track p_{T} [GeV];tracks",100,0., 20.);
2639  h["p_"+type] = dir.make<TH1F>(("p_"+type).c_str(),"track p; track p [GeV];tracks",100,0., 20.);
2640  h["found_"+type] = dir.make<TH1F>(("found_"+type).c_str(),"n. found hits;n^{found}_{hits};tracks",30, 0., 30.);
2641  h["lost_"+type] = dir.make<TH1F>(("lost_"+type).c_str(),"n. lost hits;n^{lost}_{hits};tracks",20, 0., 20.);
2642  h["nchi2_"+type] = dir.make<TH1F>(("nchi2_"+type).c_str(),"normalized track #chi^{2};track #chi^{2}/ndf;tracks",100, 0., 20.);
2643  h["rstart_"+type] = dir.make<TH1F>(("rstart_"+type).c_str(),"track start radius; track innermost radius r (cm);tracks",100, 0., 20.);
2644  h["expectedInner_"+type] = dir.make<TH1F>(("expectedInner_"+type).c_str(),"n. expected inner hits;n^{expected}_{inner};tracks",10, 0., 10.);
2645  h["expectedOuter_"+type] = dir.make<TH1F>(("expectedOuter_"+type).c_str(),"n. expected outer hits;n^{expected}_{outer};tracks ",10, 0., 10.);
2646  h["logtresxy_"+type] = dir.make<TH1F>(("logtresxy_"+type).c_str(),"log10(track r-#phi resolution/#mum);log10(track r-#phi resolution/#mum);tracks",100, 0., 5.);
2647  h["logtresz_"+type] = dir.make<TH1F>(("logtresz_"+type).c_str(),"log10(track z resolution/#mum);log10(track z resolution/#mum);tracks",100, 0., 5.);
2648  h["tpullxy_"+type] = dir.make<TH1F>(("tpullxy_"+type).c_str(),"track r-#phi pull;pull_{r-#phi};tracks",100, -10., 10.);
2649  h["tpullz_"+type] = dir.make<TH1F>(("tpullz_"+type).c_str(),"track r-z pull;pull_{r-z};tracks",100, -50., 50.);
2650  h["tlogDCAxy_"+type] = dir.make<TH1F>(("tlogDCAxy_"+type).c_str(),"track log_{10}(DCA_{r-#phi});track log_{10}(DCA_{r-#phi});tracks",200, -5., 3.);
2651  h["tlogDCAz_"+type] = dir.make<TH1F>(("tlogDCAz_"+type).c_str(),"track log_{10}(DCA_{r-z});track log_{10}(DCA_{r-z});tracks",200, -5., 5.);
2652  h["lvseta_"+type] = dir.make<TH2F>(("lvseta_"+type).c_str(),"cluster length vs #eta;track #eta;cluster length",60,-3., 3., 20, 0., 20);
2653  h["lvstanlambda_"+type] = dir.make<TH2F>(("lvstanlambda_"+type).c_str(),"cluster length vs tan #lambda; tan#lambda;cluster length",60,-6., 6., 20, 0., 20);
2654  h["restrkz_"+type] = dir.make<TH1F>(("restrkz_"+type).c_str(),"z-residuals (track vs vertex);res_{z} (cm);tracks", 200, -5., 5.);
2655  h["restrkzvsphi_"+type] = dir.make<TH2F>(("restrkzvsphi_"+type).c_str(),"z-residuals (track - vertex) vs track #phi;track #phi;res_{z} (cm)", 12,-M_PI,M_PI,100, -0.5,0.5);
2656  h["restrkzvseta_"+type] = dir.make<TH2F>(("restrkzvseta_"+type).c_str(),"z-residuals (track - vertex) vs track #eta;track #eta;res_{z} (cm)", 12,-3.,3.,200, -0.5,0.5);
2657  h["pulltrkzvsphi_"+type] = dir.make<TH2F>(("pulltrkzvsphi_"+type).c_str(),"normalized z-residuals (track - vertex) vs track #phi;track #phi;res_{z}/#sigma_{res_{z}}", 12,-M_PI,M_PI,100, -5., 5.);
2658  h["pulltrkzvseta_"+type] = dir.make<TH2F>(("pulltrkzvseta_"+type).c_str(),"normalized z-residuals (track - vertex) vs track #eta;track #eta;res_{z}/#sigma_{res_{z}}", 12,-3.,3.,100, -5., 5.);
2659  h["pulltrkz_"+type] = dir.make<TH1F>(("pulltrkz_"+type).c_str(),"normalized z-residuals (track vs vertex);res_{z}/#sigma_{res_{z}};tracks", 100, -5., 5.);
2660  h["sigmatrkz0_"+type] = dir.make<TH1F>(("sigmatrkz0_"+type).c_str(),"z-resolution (excluding beam);#sigma^{trk}_{z_{0}} (cm);tracks", 100, 0., 5.);
2661  h["sigmatrkz_"+type] = dir.make<TH1F>(("sigmatrkz_"+type).c_str(),"z-resolution (including beam);#sigma^{trk}_{z} (cm);tracks", 100,0., 5.);
2662  h["nbarrelhits_"+type] = dir.make<TH1F>(("nbarrelhits_"+type).c_str(),"number of pixel barrel hits;n. hits Barrel Pixel;tracks", 10, 0., 10.);
2663  h["nbarrelLayers_"+type] = dir.make<TH1F>(("nbarrelLayers_"+type).c_str(),"number of pixel barrel layers;n. layers Barrel Pixel;tracks", 10, 0., 10.);
2664  h["nPxLayers_"+type] = dir.make<TH1F>(("nPxLayers_"+type).c_str(),"number of pixel layers (barrel+endcap);n. Pixel layers;tracks", 10, 0., 10.);
2665  h["nSiLayers_"+type] = dir.make<TH1F>(("nSiLayers_"+type).c_str(),"number of Tracker layers;n. Tracker layers;tracks", 20, 0., 20.);
2666  h["trackAlgo_"+type] = dir.make<TH1F>(("trackAlgo_"+type).c_str(),"track algorithm;track algo;tracks", 30, 0., 30.);
2667  h["trackQuality_"+type] = dir.make<TH1F>(("trackQuality_"+type).c_str(),"track quality;track quality;tracks", 7, -1., 6.);
2668  }
2669 
2670  return h;
2671 
2672 }
2673 
2674 //*************************************************************
2675 // Generic booker function
2676 //*************************************************************
2677 std::vector<TH1F*> PrimaryVertexValidation::bookResidualsHistogram(const TFileDirectory& dir,
2678  unsigned int theNOfBins,
2679  PVValHelper::residualType resType,
2680  PVValHelper::plotVariable varType,
2681  bool isNormalized){
2682 
2683  TH1F::SetDefaultSumw2(kTRUE);
2684 
2685  auto hash = std::make_pair(resType,varType);
2686 
2687  double down = theDetails_.range[hash].first;
2688  double up = theDetails_.range[hash].second;
2689 
2690  if(isNormalized){
2691  up = up/100.;
2692  down = down/100.;
2693  }
2694 
2695  std::vector<TH1F*> h;
2696  h.reserve(theNOfBins);
2697 
2698  if (theNOfBins==0) {
2699  edm::LogError("PrimaryVertexValidation") <<"bookResidualsHistogram() The number of bins cannot be identically 0" << std::endl;
2700  assert(false);
2701  }
2702 
2703  std::string s_resType = std::get<0>(PVValHelper::getTypeString(resType));
2704  std::string s_varType = std::get<0>(PVValHelper::getVarString(varType));
2705 
2706  std::string t_resType = std::get<1>(PVValHelper::getTypeString(resType));
2707  std::string t_varType = std::get<1>(PVValHelper::getVarString(varType));
2708  std::string units = std::get<2>(PVValHelper::getTypeString(resType));
2709 
2710  for(unsigned int i=0; i<theNOfBins;i++){
2711 
2712  TString title = (varType == PVValHelper::phi || varType == PVValHelper::eta) ?
2713  Form("%s vs %s - bin %i (%f < %s < %f);%s %s;tracks",t_resType.c_str(),t_varType.c_str(),i, theDetails_.trendbins[varType][i],t_varType.c_str(),theDetails_.trendbins[varType][i+1],t_resType.c_str(),units.c_str()) : 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());
2714 
2715  TH1F* htemp = dir.make<TH1F>(Form("histo_%s_%s_plot%i",s_resType.c_str(),s_varType.c_str(),i),
2716  //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()),
2717  title.Data(),
2718  theDetails_.histobins,down,up);
2719  h.push_back(htemp);
2720  }
2721 
2722  return h;
2723 
2724 }
2725 
2726 //*************************************************************
2727 void PrimaryVertexValidation::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_)
2728 //*************************************************************
2729 {
2730 
2731  using namespace reco;
2732 
2733  PVValHelper::fill(h,"pseudorapidity_"+ttype,tt->track().eta());
2734  PVValHelper::fill(h,"z0_"+ttype,tt->track().vz());
2735  PVValHelper::fill(h,"phi_"+ttype,tt->track().phi());
2736  PVValHelper::fill(h,"eta_"+ttype,tt->track().eta());
2737  PVValHelper::fill(h,"pt_"+ttype,tt->track().pt());
2738  PVValHelper::fill(h,"p_"+ttype,tt->track().p());
2739  PVValHelper::fill(h,"found_"+ttype,tt->track().found());
2740  PVValHelper::fill(h,"lost_"+ttype,tt->track().lost());
2741  PVValHelper::fill(h,"nchi2_"+ttype,tt->track().normalizedChi2());
2742  PVValHelper::fill(h,"rstart_"+ttype,(tt->track().innerPosition()).Rho());
2743 
2744  double d0Error=tt->track().d0Error();
2745  double d0=tt->track().dxy(beamSpot.position());
2746  double dz=tt->track().dz(beamSpot.position());
2747  if (d0Error>0){
2748  PVValHelper::fill(h,"logtresxy_"+ttype,log(d0Error/0.0001)/log(10.));
2749  PVValHelper::fill(h,"tpullxy_"+ttype,d0/d0Error);
2750  PVValHelper::fill(h,"tlogDCAxy_"+ttype,log(std::abs(d0/d0Error)));
2751 
2752  }
2753  //double z0=tt->track().vz();
2754  double dzError=tt->track().dzError();
2755  if(dzError>0){
2756  PVValHelper::fill(h,"logtresz_"+ttype,log(dzError/0.0001)/log(10.));
2757  PVValHelper::fill(h,"tpullz_"+ttype,dz/dzError);
2758  PVValHelper::fill(h,"tlogDCAz_"+ttype,log(std::abs(dz/dzError)));
2759  }
2760 
2761  //
2762  double wxy2_=pow(beamSpot.BeamWidthX(),2)+pow(beamSpot.BeamWidthY(),2);
2763 
2764  PVValHelper::fill(h,"sigmatrkz_"+ttype,sqrt(pow(tt->track().dzError(),2)+wxy2_/pow(tan(tt->track().theta()),2)));
2765  PVValHelper::fill(h,"sigmatrkz0_"+ttype,tt->track().dzError());
2766 
2767  // track vs vertex
2768  if( v.isValid()){ // && (v.ndof()<10.)) {
2769  // emulate clusterizer input
2770  //const TransientTrack & tt = theB_->build(&t); wrong !!!!
2771  //reco::TransientTrack tt = theB_->build(&t);
2772  //ttt->track().setBeamSpot(beamSpot); // need the setBeamSpot !
2773  double z=(tt->stateAtBeamLine().trackStateAtPCA()).position().z();
2774  double tantheta=tan((tt->stateAtBeamLine().trackStateAtPCA()).momentum().theta());
2775  double dz2= pow(tt->track().dzError(),2)+wxy2_/pow(tantheta,2);
2776 
2777  PVValHelper::fill(h,"restrkz_"+ttype,z-v.position().z());
2778  PVValHelper::fill(h,"restrkzvsphi_"+ttype,tt->track().phi(), z-v.position().z());
2779  PVValHelper::fill(h,"restrkzvseta_"+ttype,tt->track().eta(), z-v.position().z());
2780  PVValHelper::fill(h,"pulltrkzvsphi_"+ttype,tt->track().phi(), (z-v.position().z())/sqrt(dz2));
2781  PVValHelper::fill(h,"pulltrkzvseta_"+ttype,tt->track().eta(), (z-v.position().z())/sqrt(dz2));
2782 
2783  PVValHelper::fill(h,"pulltrkz_"+ttype,(z-v.position().z())/sqrt(dz2));
2784 
2785  double x1=tt->track().vx()-beamSpot.x0(); double y1=tt->track().vy()-beamSpot.y0();
2786 
2787  double kappa=-0.002998*fBfield_*tt->track().qoverp()/cos(tt->track().theta());
2788  double D0=x1*sin(tt->track().phi())-y1*cos(tt->track().phi())-0.5*kappa*(x1*x1+y1*y1);
2789  double q=sqrt(1.-2.*kappa*D0);
2790  double s0=(x1*cos(tt->track().phi())+y1*sin(tt->track().phi()))/q;
2791  // double s1;
2792  if (std::abs(kappa*s0)>0.001){
2793  //s1=asin(kappa*s0)/kappa;
2794  }else{
2795  //double ks02=(kappa*s0)*(kappa*s0);
2796  //s1=s0*(1.+ks02/6.+3./40.*ks02*ks02+5./112.*pow(ks02,3));
2797  }
2798  // sp.ddcap=-2.*D0/(1.+q);
2799  //double zdcap=tt->track().vz()-s1/tan(tt->track().theta());
2800 
2801  }
2802  //
2803 
2804  // collect some info on hits and clusters
2805  PVValHelper::fill(h,"nbarrelLayers_"+ttype,tt->track().hitPattern().pixelBarrelLayersWithMeasurement());
2806  PVValHelper::fill(h,"nPxLayers_"+ttype,tt->track().hitPattern().pixelLayersWithMeasurement());
2807  PVValHelper::fill(h,"nSiLayers_"+ttype,tt->track().hitPattern().trackerLayersWithMeasurement());
2808  PVValHelper::fill(h,"expectedInner_"+ttype,tt->track().hitPattern().numberOfLostHits(HitPattern::MISSING_INNER_HITS));
2809  PVValHelper::fill(h,"expectedOuter_"+ttype,tt->track().hitPattern().numberOfLostHits(HitPattern::MISSING_OUTER_HITS));
2810  PVValHelper::fill(h,"trackAlgo_"+ttype,tt->track().algo());
2811  PVValHelper::fill(h,"trackQuality_"+ttype,tt->track().qualityMask());
2812 
2813  //
2814  int longesthit=0, nbarrel=0;
2815  for(trackingRecHit_iterator hit=tt->track().recHitsBegin(); hit!=tt->track().recHitsEnd(); hit++){
2816  if ((**hit).isValid() && (**hit).geographicalId().det() == DetId::Tracker ){
2817  bool barrel = DetId((**hit).geographicalId()).subdetId() == static_cast<int>(PixelSubdetector::PixelBarrel);
2818  //bool endcap = DetId::DetId((**hit).geographicalId()).subdetId() == static_cast<int>(PixelSubdetector::PixelEndcap);
2819  if (barrel){
2820  const SiPixelRecHit *pixhit = dynamic_cast<const SiPixelRecHit*>( &(**hit));
2821  edm::Ref<edmNew::DetSetVector<SiPixelCluster>, SiPixelCluster> const& clust = (*pixhit).cluster();
2822  if (clust.isNonnull()) {
2823  nbarrel++;
2824  if (clust->sizeY()-longesthit>0) longesthit=clust->sizeY();
2825  if (clust->sizeY()>20.){
2826  PVValHelper::fill(h,"lvseta_"+ttype,tt->track().eta(), 19.9);
2827  PVValHelper::fill(h,"lvstanlambda_"+ttype,tan(tt->track().lambda()), 19.9);
2828  }else{
2829  PVValHelper::fill(h,"lvseta_"+ttype,tt->track().eta(), float(clust->sizeY()));
2830  PVValHelper::fill(h,"lvstanlambda_"+ttype,tan(tt->track().lambda()), float(clust->sizeY()));
2831  }
2832  }
2833  }
2834  }
2835  }
2836  PVValHelper::fill(h,"nbarrelhits_"+ttype,float(nbarrel));
2837  //-------------------------------------------------------------------
2838 }
2839 
2840 //define this as a plug-in
2841 DEFINE_FWK_MODULE(PrimaryVertexValidation);
2842 
2843 
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:503
PrimaryVertexValidation::h_norm_dxy_modZ_
std::vector< TH1F * > h_norm_dxy_modZ_
Definition: PrimaryVertexValidation.h:638
PerigeeConversions::trajectoryStateClosestToPoint
TrajectoryStateClosestToPoint trajectoryStateClosestToPoint(const AlgebraicVector3 &momentum, const GlobalPoint &referencePoint, const TrackCharge &charge, const AlgebraicSymMatrix66 &theCovarianceMatrix, const MagneticField *field)
Definition: PerigeeConversions.cc:107
reco::TrackBase::qoverp
double qoverp() const
q / p
Definition: TrackBase.h:573
TFileService::kSharedResource
static const std::string kSharedResource
Definition: TFileService.h:76
up
Definition: BitonicSort.h:8
reco::TrackBase::p
double p() const
momentum vector magnitude
Definition: TrackBase.h:615
type
type
Definition: HCALResponse.h:21
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
edm::EventID::event
EventNumber_t event() const
Definition: EventID.h:41
reco::BeamSpot::z0
double z0() const
z coordinate
Definition: BeamSpot.h:68
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
PrimaryVertexValidation::n_IP3DPhiResiduals
std::vector< TH1F * > n_IP3DPhiResiduals
Definition: PrimaryVertexValidation.h:313
PVValHelper::getMedian
Measurement1D getMedian(TH1F *histo)
Definition: PVValidationHelpers.cc:183
PrimaryVertexValidation::a_dxyEtaResiduals
std::vector< TH1F * > a_dxyEtaResiduals
Definition: PrimaryVertexValidation.h:279
PrimaryVertexValidation::n_dzResidualsMap
TH1F * n_dzResidualsMap[nMaxBins_][nMaxBins_]
Definition: PrimaryVertexValidation.h:329
reco::HitPattern::getLayer
static uint32_t getLayer(uint16_t pattern)
Definition: HitPattern.h:700
PrimaryVertexValidation::a_d3DEtaResiduals
std::vector< TH1F * > a_d3DEtaResiduals
Definition: PrimaryVertexValidation.h:300
reco::TrackBase::d0Error
double d0Error() const
error on d0
Definition: TrackBase.h:802
PerigeeTrajectoryError::longitudinalImpactParameterError
double longitudinalImpactParameterError() const
Definition: PerigeeTrajectoryError.h:75
PrimaryVertexValidation::n_reszPhiResiduals
std::vector< TH1F * > n_reszPhiResiduals
Definition: PrimaryVertexValidation.h:316
PrimaryVertexValidation::runControlNumbers_
std::vector< unsigned int > runControlNumbers_
Definition: PrimaryVertexValidation.h:152
edm::Event::eventAuxiliary
EventAuxiliary const & eventAuxiliary() const override
Definition: Event.h:92
edm::Ref::isNonnull
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:253
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:269
reco::Track::lost
unsigned short lost() const
Number of lost (=invalid) hits on track.
Definition: Track.h:199
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
h
FWCore Framework interface EventSetupRecordImplementation h
Helper function to determine trigger accepts.
Definition: L1TUtmAlgorithmRcd.h:4
PrimaryVertexValidation::vtxSort
static bool vtxSort(const reco::Vertex &a, const reco::Vertex &b)
Definition: PrimaryVertexValidation.cc:2586
Reference_intrackfit_cff.barrel
barrel
Definition: Reference_intrackfit_cff.py:37
TrajectoryStateClosestToPoint::perigeeError
const PerigeeTrajectoryError & perigeeError() const
Definition: TrajectoryStateClosestToPoint.h:83
PVValHelper::fitResiduals
std::pair< Measurement1D, Measurement1D > fitResiduals(TH1 *hist)
Definition: PVValidationHelpers.cc:230
reco::Vertex::tracks_end
trackRef_iterator tracks_end() const
last iterator over tracks
Definition: Vertex.cc:81
reco::TransientTrack::setBeamSpot
void setBeamSpot(const reco::BeamSpot &beamSpot)
Definition: TransientTrack.h:76
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:561
PrimaryVertexValidation::fillTrendPlot
void fillTrendPlot(TH1F *trendPlot, TH1F *residualsPlot[100], PVValHelper::estimator fitPar_, const std::string &var_)
Definition: PrimaryVertexValidation.cc:2394
edm::Event::getByToken
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:519
PrimaryVertexValidation::h_dxy_pT_
std::vector< TH1F * > h_dxy_pT_
Definition: PrimaryVertexValidation.h:624
reco::Vertex::zError
double zError() const
error on z
Definition: Vertex.h:123
PrimaryVertexValidation::theTrackFilter_
TrackFilterForPVFindingBase * theTrackFilter_
Definition: PrimaryVertexValidation.h:115
PrimaryVertexValidation::a_dxyResidualsMap
TH1F * a_dxyResidualsMap[nMaxBins_][nMaxBins_]
Definition: PrimaryVertexValidation.h:324
PrimaryVertexValidation::hDA
std::map< std::string, TH1 * > hDA
Definition: PrimaryVertexValidation.h:620
reco::TrackBase::TrackQuality
TrackQuality
track quality
Definition: TrackBase.h:151
PrimaryVertexValidation::a_IP3DEtaResiduals
std::vector< TH1F * > a_IP3DEtaResiduals
Definition: PrimaryVertexValidation.h:294
VertexException
Common base class.
Definition: VertexException.h:12
TrajectoryStateClosestToPoint::theState
const FreeTrajectoryState & theState() const
Definition: TrajectoryStateClosestToPoint.h:106
reco::TrackBase::theta
double theta() const
polar angle
Definition: TrackBase.h:579
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
PVValHelper::generateBins
std::vector< float > generateBins(int n, float start, float range)
Definition: PVValidationHelpers.cc:164
PrimaryVertexValidation::n_dzPhiResiduals
std::vector< TH1F * > n_dzPhiResiduals
Definition: PrimaryVertexValidation.h:307
reco::Vertex::isValid
bool isValid() const
Tells whether the vertex is valid.
Definition: Vertex.h:68
funct::D0
Divides< arg, void > D0
Definition: Factorize.h:144
PrimaryVertexValidation::a_d3DPhiResiduals
std::vector< TH1F * > a_d3DPhiResiduals
Definition: PrimaryVertexValidation.h:299
PrimaryVertexValidation::n_d3DPhiResiduals
std::vector< TH1F * > n_d3DPhiResiduals
Definition: PrimaryVertexValidation.h:319
PrimaryVertexValidation::PrimaryVertexValidation
PrimaryVertexValidation(const edm::ParameterSet &)
Definition: PrimaryVertexValidation.cc:69
PrimaryVertexValidation::n_dxyPhiBiasResiduals
std::vector< TH1F * > n_dxyPhiBiasResiduals
Definition: PrimaryVertexValidation.h:462
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:50
hfClusterShapes_cfi.hits
hits
Definition: hfClusterShapes_cfi.py:5
PrimaryVertexValidation::n_dxyBiasResidualsMap
TH1F * n_dxyBiasResidualsMap[nMaxBins_][nMaxBins_]
Definition: PrimaryVertexValidation.h:473
PrimaryVertexValidation::n_d3DEtaResiduals
std::vector< TH1F * > n_d3DEtaResiduals
Definition: PrimaryVertexValidation.h:320
reco::TransientTrack::hitPattern
const HitPattern & hitPattern() const
Definition: TransientTrack.h:138
TrackerTopology::pxbLadder
unsigned int pxbLadder(const DetId &id) const
Definition: TrackerTopology.h:155
GlobalPoint
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
TransientTrackBuilder::build
reco::TransientTrack build(const reco::Track *p) const
Definition: TransientTrackBuilder.cc:20
PrimaryVertexValidation::h_norm_dz_Central_pT_
std::vector< TH1F * > h_norm_dz_Central_pT_
Definition: PrimaryVertexValidation.h:632
IPTools::absoluteImpactParameter3D
std::pair< bool, Measurement1D > absoluteImpactParameter3D(const reco::TransientTrack &transientTrack, const reco::Vertex &vertex)
Definition: IPTools.cc:37
PV3DBase::y
T y() const
Definition: PV3DBase.h:63
Measurement1D::error
double error() const
Definition: Measurement1D.h:30
edm::EventAuxiliary::run
RunNumber_t run() const
Definition: EventAuxiliary.h:67
IPTools::signedImpactParameter3D
std::pair< bool, Measurement1D > signedImpactParameter3D(const reco::TransientTrack &track, const GlobalVector &direction, const reco::Vertex &vertex)
Definition: IPTools.cc:71
electrons_cff.vertices
vertices
Definition: electrons_cff.py:11
reco::TrackBase::phi
double phi() const
azimuthal angle of momentum vector
Definition: TrackBase.h:645
PrimaryVertexValidation::DOFUnbiasedVertex_
float DOFUnbiasedVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:239
TrackerTopology::pxbModule
unsigned int pxbModule(const DetId &id) const
Definition: TrackerTopology.h:160
TFileService::make
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
PrimaryVertexValidation::d3DFromMyVertex_
double d3DFromMyVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:245
cms::Exception::what
char const * what() const override
Definition: Exception.cc:141
PrimaryVertexValidation::fillTrendPlotByIndex
void fillTrendPlotByIndex(TH1F *trendPlot, std::vector< TH1F * > &h, PVValHelper::estimator fitPar_, PVValHelper::plotVariable plotVar=PVValHelper::END_OF_PLOTS)
Definition: PrimaryVertexValidation.cc:2458
PrimaryVertexValidation::a_reszPhiResiduals
std::vector< TH1F * > a_reszPhiResiduals
Definition: PrimaryVertexValidation.h:296
DAClusterizerInZ_vect
Definition: DAClusterizerInZ_vect.h:21
PrimaryVertexValidation::a_IP2DPhiResiduals
std::vector< TH1F * > a_IP2DPhiResiduals
Definition: PrimaryVertexValidation.h:290
PrimaryVertexValidation::chi2ProbUnbiasedVertex_
float chi2ProbUnbiasedVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:238
reco::TrackBase::px
double px() const
x coordinate of momentum vector
Definition: TrackBase.h:627
reco::HitPattern::pixelLayersWithMeasurement
int pixelLayersWithMeasurement() const
Definition: HitPattern.cc:499
PrimaryVertexValidation::n_reszEtaResiduals
std::vector< TH1F * > n_reszEtaResiduals
Definition: PrimaryVertexValidation.h:317
PrimaryVertexValidation::dxyFromMyVertex_
double dxyFromMyVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:243
PVValHelper::getMAD
Measurement1D getMAD(TH1F *histo)
Definition: PVValidationHelpers.cc:200
PrimaryVertexValidation::bookVertexHistograms
std::map< std::string, TH1 * > bookVertexHistograms(const TFileDirectory &dir)
Definition: PrimaryVertexValidation.cc:2623
reco::HitPattern::trackerLayersWithMeasurement
int trackerLayersWithMeasurement() const
Definition: HitPattern.cc:518
reco::Vertex::position
const Point & position() const
position
Definition: Vertex.h:109
reco::HitPattern::pixelBarrelHitFilter
static bool pixelBarrelHitFilter(uint16_t pattern)
Definition: HitPattern.h:575
PrimaryVertexValidation::h_dz_modZ_
std::vector< TH1F * > h_dz_modZ_
Definition: PrimaryVertexValidation.h:637
PrimaryVertexValidation::h_dxy_modZ_
std::vector< TH1F * > h_dxy_modZ_
Definition: PrimaryVertexValidation.h:636
PrimaryVertexValidation::a_dxEtaResiduals
std::vector< TH1F * > a_dxEtaResiduals
Definition: PrimaryVertexValidation.h:282
PrimaryVertexValidation::isBFieldConsistentWithMode
bool isBFieldConsistentWithMode(const edm::EventSetup &iSetup) const
Definition: PrimaryVertexValidation.cc:2305
PrimaryVertexValidation::pixelHitsCheck
std::pair< bool, bool > pixelHitsCheck(const reco::TransientTrack &track)
Definition: PrimaryVertexValidation.cc:1114
reco::TransientTrack::stateAtBeamLine
TrajectoryStateClosestToBeamLine stateAtBeamLine() const
Definition: TransientTrack.h:124
reco::HitPattern::numberOfValidStripTOBHits
int numberOfValidStripTOBHits() const
Definition: HitPattern.h:868
ecalDrivenElectronSeedsParameters_cff.beamSpot
beamSpot
Definition: ecalDrivenElectronSeedsParameters_cff.py:11
PVValHelper::getVarString
plotLabels getVarString(plotVariable var)
Definition: PVValidationHelpers.cc:131
PrimaryVertexValidation::dzErrorFromMyVertex_
double dzErrorFromMyVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:248
types
types
Definition: AlignPCLThresholds_PayloadInspector.cc:30
SoftLeptonByDistance_cfi.distance
distance
Definition: SoftLeptonByDistance_cfi.py:6
PrimaryVertexValidation::n_dxyEtaBiasResiduals
std::vector< TH1F * > n_dxyEtaBiasResiduals
Definition: PrimaryVertexValidation.h:463
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:2727
edm::EventAuxiliary::luminosityBlock
LuminosityBlockNumber_t luminosityBlock() const
Definition: EventAuxiliary.h:63
reco::HitPattern::pixelEndcapHitFilter
static bool pixelEndcapHitFilter(uint16_t pattern)
Definition: HitPattern.h:585
PrimaryVertexValidation::n_dxyEtaResiduals
std::vector< TH1F * > n_dxyEtaResiduals
Definition: PrimaryVertexValidation.h:305
PrimaryVertexValidation::IPTsigFromMyVertex_
double IPTsigFromMyVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:251
reco::HitPattern::validHitFilter
static bool validHitFilter(uint16_t pattern)
Definition: HitPattern.h:787
TrackerGeometry::isThere
bool isThere(GeomDetEnumerators::SubDetector subdet) const
Definition: TrackerGeometry.cc:244
reco::Track::innerPosition
const math::XYZPoint & innerPosition() const
position of the innermost hit
Definition: Track.h:55
reco::TrackBase::algo
TrackAlgorithm algo() const
Definition: TrackBase.h:497
PrimaryVertexValidation::a_dzPhiBiasResiduals
std::vector< TH1F * > a_dzPhiBiasResiduals
Definition: PrimaryVertexValidation.h:457
PrimaryVertexValidation::h_norm_dxy_pT_
std::vector< TH1F * > h_norm_dxy_pT_
Definition: PrimaryVertexValidation.h:626
TrackerGeometry::dets
const DetContainer & dets() const override
Returm a vector of all GeomDet (including all GeomDetUnits)
Definition: TrackerGeometry.h:61
PrimaryVertexValidation::a_dzEtaResiduals
std::vector< TH1F * > a_dzEtaResiduals
Definition: PrimaryVertexValidation.h:288
reco::TrackBase::vertex
const Point & vertex() const
reference point on the track. This method is DEPRECATED, please use referencePoint() instead ...
Definition: TrackBase.h:687
PVValHelper::histodetails::getHigh
float getHigh(residualType type, plotVariable plot)
Definition: PVValidationHelpers.h:102
PrimaryVertexValidation::tracksUsedForVertexing_
int tracksUsedForVertexing_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:241
PVValHelper::histodetails::getLow
float getLow(residualType type, plotVariable plot)
Definition: PVValidationHelpers.h:93
iEvent
int iEvent
Definition: GenABIO.cc:230
reco::TrackBase::eta
double eta() const
pseudorapidity of momentum vector
Definition: TrackBase.h:651
PrimaryVertexValidation::a_dxyEtaBiasResiduals
std::vector< TH1F * > a_dxyEtaBiasResiduals
Definition: PrimaryVertexValidation.h:455
reco::HitPattern::numberOfValidPixelBarrelHits
int numberOfValidPixelBarrelHits() const
Definition: HitPattern.h:843
PrimaryVertexValidation::h_norm_dxy_ladder_
std::vector< TH1F * > h_norm_dxy_ladder_
Definition: PrimaryVertexValidation.h:646
TrajectoryStateClosestToPoint::perigeeParameters
const PerigeeTrajectoryParameters & perigeeParameters() const
Definition: TrajectoryStateClosestToPoint.h:68
PrimaryVertexValidation::a_dzPhiResiduals
std::vector< TH1F * > a_dzPhiResiduals
Definition: PrimaryVertexValidation.h:287
reco::TrackBase::chi2
double chi2() const
chi-squared of the fit
Definition: TrackBase.h:549
PrimaryVertexValidation::yUnbiasedVertex_
double yUnbiasedVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:234
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
TrackingRecHit::dimension
virtual int dimension() const =0
reco::TrackBase::pt
double pt() const
track transverse momentum
Definition: TrackBase.h:621
PrimaryVertexValidation::n_dxyPhiResiduals
std::vector< TH1F * > n_dxyPhiResiduals
Definition: PrimaryVertexValidation.h:304
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:763
reco::TrackBase::qualityMask
int qualityMask() const
Definition: TrackBase.h:862
reco::HitPattern::numberOfAllHits
int numberOfAllHits(HitCategory category) const
Definition: HitPattern.h:807
PVValHelper::histodetails::trendbins
std::map< plotVariable, std::vector< float > > trendbins
Definition: PVValidationHelpers.h:77
funct::tan
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
MetAnalyzer.pv
def pv(vc)
Definition: MetAnalyzer.py:6
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
reco::Vertex::chi2
double chi2() const
chi-squares
Definition: Vertex.h:98
PrimaryVertexValidation::zUnbiasedVertex_
double zUnbiasedVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:235
DAClusterizerInZ_vect.h
reco::HitPattern::numberOfValidStripTIDHits
int numberOfValidStripTIDHits() const
Definition: HitPattern.h:863
reco::TrackBase::lambda
double lambda() const
Lambda angle.
Definition: TrackBase.h:585
reco::BeamSpot::BeamWidthX
double BeamWidthX() const
beam width X
Definition: BeamSpot.h:86
PrimaryVertexValidation::h_dxy_ladder_
std::vector< TH1F * > h_dxy_ladder_
Definition: PrimaryVertexValidation.h:644
reco::TrackBase::numberOfValidHits
unsigned short numberOfValidHits() const
number of valid hits found
Definition: TrackBase.h:820
reco::HitPattern::numberOfValidStripTECHits
int numberOfValidStripTECHits() const
Definition: HitPattern.h:873
PrimaryVertexValidation::chi2normUnbiasedVertex_
float chi2normUnbiasedVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:236
PrimaryVertexValidation::n_d3DResidualsMap
TH1F * n_d3DResidualsMap[nMaxBins_][nMaxBins_]
Definition: PrimaryVertexValidation.h:330
TFileDirectory::make
T * make(const Args &...args) const
make new ROOT object
Definition: TFileDirectory.h:59
min
T min(T a, T b)
Definition: MathUtil.h:58
PrimaryVertexValidation::a_dxyPhiResiduals
std::vector< TH1F * > a_dxyPhiResiduals
Definition: PrimaryVertexValidation.h:278
reco::Track::recHitsBegin
trackingRecHit_iterator recHitsBegin() const
Iterator to first hit on the track.
Definition: Track.h:104
PrimaryVertexValidation::a_dxyPhiBiasResiduals
std::vector< TH1F * > a_dxyPhiBiasResiduals
Definition: PrimaryVertexValidation.h:454
PrimaryVertexValidation::theVertexCollectionToken
edm::EDGetTokenT< reco::VertexCollection > theVertexCollectionToken
Definition: PrimaryVertexValidation.h:158
PrimaryVertexValidation::analyze
void analyze(const edm::Event &, const edm::EventSetup &) override
Definition: PrimaryVertexValidation.cc:200
PrimaryVertexValidation::a_d3DResidualsMap
TH1F * a_d3DResidualsMap[nMaxBins_][nMaxBins_]
Definition: PrimaryVertexValidation.h:326
PrimaryVertexValidation::h_dxy_ladderNoOverlap_
std::vector< TH1F * > h_dxy_ladderNoOverlap_
Definition: PrimaryVertexValidation.h:642
edm::HandleBase::isValid
bool isValid() const
Definition: HandleBase.h:74
DetId::subdetId
int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:38
PrimaryVertexValidation::IPLsigFromMyVertex_
double IPLsigFromMyVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:252
LogTrace
#define LogTrace(id)
Definition: MessageLogger.h:504
globals_cff.x1
x1
Definition: globals_cff.py:24
PrimaryVertexValidation::hasFirstLayerPixelHits
bool hasFirstLayerPixelHits(const reco::TransientTrack &track)
Definition: PrimaryVertexValidation.cc:1132
PrimaryVertexValidation::n_dzEtaResiduals
std::vector< TH1F * > n_dzEtaResiduals
Definition: PrimaryVertexValidation.h:308
reco::BeamSpot::dxdz
double dxdz() const
dxdz slope
Definition: BeamSpot.h:82
reco::Vertex::ndof
double ndof() const
Definition: Vertex.h:105
PVValHelper::histodetails::setMap
void setMap(residualType type, plotVariable plot, float low, float high)
Definition: PVValidationHelpers.h:79
M_PI
#define M_PI
Definition: BXVectorInputProducer.cc:51
reco::TrackBase::pz
double pz() const
z coordinate of momentum vector
Definition: TrackBase.h:639
TrackClusterizerInZ::clusterize
virtual std::vector< std::vector< reco::TransientTrack > > clusterize(const std::vector< reco::TransientTrack > &tracks) const =0
stringResolutionProvider_cfi.bin
bin
set the eta bin as selection string.
Definition: stringResolutionProvider_cfi.py:11
PrimaryVertexValidation::n_IP3DEtaResiduals
std::vector< TH1F * > n_IP3DEtaResiduals
Definition: PrimaryVertexValidation.h:314
TrackerTopology::pxbLayer
unsigned int pxbLayer(const DetId &id) const
Definition: TrackerTopology.h:144
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:609
reco::TrackBase::dzError
double dzError() const
error on dz
Definition: TrackBase.h:814
electrons_cff.bool
bool
Definition: electrons_cff.py:231
PVValHelper::histodetails::range
std::map< std::pair< residualType, plotVariable >, std::pair< float, float > > range
Definition: PVValidationHelpers.h:76
reco::TrackBase::vz
double vz() const
z coordinate of the reference point on track
Definition: TrackBase.h:669
PrimaryVertexValidation::h_norm_dz_pT_
std::vector< TH1F * > h_norm_dz_pT_
Definition: PrimaryVertexValidation.h:627
DetId
Definition: DetId.h:18
FreeTrajectoryState::position
GlobalPoint position() const
Definition: FreeTrajectoryState.h:84
PrimaryVertexValidation::fillMap
void fillMap(TH2F *trendMap, TH1F *residualsMapPlot[100][100], PVValHelper::estimator fitPar_)
Definition: PrimaryVertexValidation.cc:2524
reco::TrackBase::qualityByName
static TrackQuality qualityByName(const std::string &name)
Definition: TrackBase.cc:125
reco::Vertex::xError
double xError() const
error on x
Definition: Vertex.h:119
PrimaryVertexValidation::h_dz_pT_
std::vector< TH1F * > h_dz_pT_
Definition: PrimaryVertexValidation.h:625
PVValHelper::getTypeString
plotLabels getTypeString(residualType type)
Definition: PVValidationHelpers.cc:62
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:55
fastPrimaryVertexProducer_cfi.clusters
clusters
Definition: fastPrimaryVertexProducer_cfi.py:5
PrimaryVertexValidation::n_dzEtaBiasResiduals
std::vector< TH1F * > n_dzEtaBiasResiduals
Definition: PrimaryVertexValidation.h:466
PrimaryVertexValidation::fs
edm::Service< TFileService > fs
Definition: PrimaryVertexValidation.h:258
PrimaryVertexValidation::h_dz_Central_pT_
std::vector< TH1F * > h_dz_Central_pT_
Definition: PrimaryVertexValidation.h:630
reco::TransientTrack::track
const Track & track() const
Definition: TransientTrack.h:120
reco::TrackBase::hitPattern
const HitPattern & hitPattern() const
Access the hit pattern, indicating in which Tracker layers the track has hits.
Definition: TrackBase.h:446
HiIsolationCommonParameters_cff.track
track
Definition: HiIsolationCommonParameters_cff.py:9
PrimaryVertexValidation::isHit2D
bool isHit2D(const TrackingRecHit &hit) const
Definition: PrimaryVertexValidation.cc:1083
PrimaryVertexValidation::h_dz_ladder_
std::vector< TH1F * > h_dz_ladder_
Definition: PrimaryVertexValidation.h:645
PrimaryVertexValidation::a_IP2DEtaResiduals
std::vector< TH1F * > a_IP2DEtaResiduals
Definition: PrimaryVertexValidation.h:291
math::XYZPoint
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
edm::EventSetup::get
const T & get() const
Definition: EventSetup.h:59
reco::BeamSpot::sigmaZ
double sigmaZ() const
sigma z
Definition: BeamSpot.h:80
reco::HitPattern::pixelBarrelLayersWithMeasurement
int pixelBarrelLayersWithMeasurement() const
Definition: HitPattern.cc:558
reco::BeamSpot::BeamWidthY
double BeamWidthY() const
beam width Y
Definition: BeamSpot.h:88
PrimaryVertexValidation::theTrackClusterizer_
TrackClusterizerInZ * theTrackClusterizer_
Definition: PrimaryVertexValidation.h:116
b
double b
Definition: hdecay.h:120
PrimaryVertexValidation::n_IP2DPhiResiduals
std::vector< TH1F * > n_IP2DPhiResiduals
Definition: PrimaryVertexValidation.h:310
PrimaryVertexValidation::xUnbiasedVertex_
double xUnbiasedVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:233
PrimaryVertexValidation::h_dxy_Central_pT_
std::vector< TH1F * > h_dxy_Central_pT_
Definition: PrimaryVertexValidation.h:629
Measurement1D::value
double value() const
Definition: Measurement1D.h:28
TrackFilterForPVFindingBase::select
virtual std::vector< reco::TransientTrack > select(const std::vector< reco::TransientTrack > &tracks) const =0
reco::HitPattern::numberOfValidStripTIBHits
int numberOfValidStripTIBHits() const
Definition: HitPattern.h:858
reco::TrackBase::quality
bool quality(const TrackQuality) const
Track quality.
Definition: TrackBase.h:510
reco::HitPattern::numberOfLostHits
int numberOfLostHits(HitCategory category) const
Definition: HitPattern.h:902
reco::HitPattern::numberOfValidPixelEndcapHits
int numberOfValidPixelEndcapHits() const
Definition: HitPattern.h:848
RunInfo::m_avg_current
float m_avg_current
Definition: RunInfo.h:29
PrimaryVertexValidation::theBeamspotToken
edm::EDGetTokenT< reco::BeamSpot > theBeamspotToken
Definition: PrimaryVertexValidation.h:159
PrimaryVertexValidation::theTrackCollectionToken
edm::EDGetTokenT< reco::TrackCollection > theTrackCollectionToken
Definition: PrimaryVertexValidation.h:157
fftjetcommon_cfi.title
title
Definition: fftjetcommon_cfi.py:32
edm::EventAuxiliary::id
EventID const & id() const
Definition: EventAuxiliary.h:59
SiPixelCluster
Pixel cluster – collection of neighboring pixels above threshold.
Definition: SiPixelCluster.h:27
reco::TrackBase::vy
double vy() const
y coordinate of the reference point on track
Definition: TrackBase.h:663
units
TString units(TString variable, Char_t axis)
PrimaryVertexValidation::dzFromMyVertex_
double dzFromMyVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:244
reco
fixed size matrix
Definition: AlignmentAlgorithmBase.h:43
PrimaryVertexValidation::mypT_bins_
std::array< float, nPtBins_+1 > mypT_bins_
Definition: PrimaryVertexValidation.h:179
a
double a
Definition: hdecay.h:121
allConversions_cfi.d0
d0
Definition: allConversions_cfi.py:30
position
static int position[264][3]
Definition: ReadPGInfo.cc:509
reco::Track::found
unsigned short found() const
Number of valid hits on track.
Definition: Track.h:194
edm::OwnVector::const_iterator
Definition: OwnVector.h:41
PrimaryVertexValidation::a_reszEtaResiduals
std::vector< TH1F * > a_reszEtaResiduals
Definition: PrimaryVertexValidation.h:297
PrimaryVertexValidation::bookResidualsHistogram
std::vector< TH1F * > bookResidualsHistogram(const TFileDirectory &dir, unsigned int theNOfBins, PVValHelper::residualType resType, PVValHelper::plotVariable varType, bool isNormalized=false)
Definition: PrimaryVertexValidation.cc:2677
reco::Vertex::trackRef_iterator
std::vector< TrackBaseRef >::const_iterator trackRef_iterator
The iteratator for the vector<TrackRef>
Definition: Vertex.h:37
reco::BeamSpot::y0
double y0() const
y coordinate
Definition: BeamSpot.h:66
PrimaryVertexValidation::a_dxyBiasResidualsMap
TH1F * a_dxyBiasResidualsMap[nMaxBins_][nMaxBins_]
Definition: PrimaryVertexValidation.h:470
PrimaryVertexValidation::a_dxPhiResiduals
std::vector< TH1F * > a_dxPhiResiduals
Definition: PrimaryVertexValidation.h:281
PrimaryVertexValidation::h_norm_dz_modZ_
std::vector< TH1F * > h_norm_dz_modZ_
Definition: PrimaryVertexValidation.h:639
reco::TrackBase::charge
int charge() const
track electric charge
Definition: TrackBase.h:567
reco::BeamSpot::position
const Point & position() const
position
Definition: BeamSpot.h:62
PrimaryVertexValidation::n_dxyResidualsMap
TH1F * n_dxyResidualsMap[nMaxBins_][nMaxBins_]
Definition: PrimaryVertexValidation.h:328
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:2596
kappa
static const G4double kappa
Definition: UrbanMscModel93.cc:34
reco::Vertex::normalizedChi2
double normalizedChi2() const
chi-squared divided by n.d.o.f.
Definition: Vertex.h:107
PrimaryVertexValidation::dxyErrorFromMyVertex_
double dxyErrorFromMyVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:247
TrackingRecHit::geographicalId
DetId geographicalId() const
Definition: TrackingRecHit.h:116
dir
dbl *** dir
Definition: mlp_gen.cc:35
PrimaryVertexValidation::a_dyEtaResiduals
std::vector< TH1F * > a_dyEtaResiduals
Definition: PrimaryVertexValidation.h:285
reco::Vertex::tracks_begin
trackRef_iterator tracks_begin() const
first iterator over tracks
Definition: Vertex.cc:76
PrimaryVertexValidation::n_IP2DEtaResiduals
std::vector< TH1F * > n_IP2DEtaResiduals
Definition: PrimaryVertexValidation.h:311
PrimaryVertexValidation::h_norm_dz_ladder_
std::vector< TH1F * > h_norm_dz_ladder_
Definition: PrimaryVertexValidation.h:647
reco::HitPattern::getHitPattern
uint16_t getHitPattern(HitCategory category, int position) const
Definition: HitPattern.h:515
PrimaryVertexValidation::sumOfWeightsUnbiasedVertex_
float sumOfWeightsUnbiasedVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:240
PrimaryVertexValidation::n_dzBiasResidualsMap
TH1F * n_dzBiasResidualsMap[nMaxBins_][nMaxBins_]
Definition: PrimaryVertexValidation.h:474
PrimaryVertexValidation::d3DErrorFromMyVertex_
double d3DErrorFromMyVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:249
PrimaryVertexValidation::a_dyPhiResiduals
std::vector< TH1F * > a_dyPhiResiduals
Definition: PrimaryVertexValidation.h:284
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:591
PrimaryVertexValidation::h_dxy_ladderOverlap_
std::vector< TH1F * > h_dxy_ladderOverlap_
Definition: PrimaryVertexValidation.h:641
PrimaryVertexValidation::a_dzBiasResidualsMap
TH1F * a_dzBiasResidualsMap[nMaxBins_][nMaxBins_]
Definition: PrimaryVertexValidation.h:471
PV3DBase::x
T x() const
Definition: PV3DBase.h:62
PrimaryVertexValidation::theDetails_
PVValHelper::histodetails theDetails_
Definition: PrimaryVertexValidation.h:127
edm::ESHandle::product
T const * product() const
Definition: ESHandle.h:86
CombinatorialSeedGeneratorForCosmicsRegionalReconstruction_cfi.zVertex
zVertex
Definition: CombinatorialSeedGeneratorForCosmicsRegionalReconstruction_cfi.py:44
PrimaryVertexValidation::a_dzResidualsMap
TH1F * a_dzResidualsMap[nMaxBins_][nMaxBins_]
Definition: PrimaryVertexValidation.h:325
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
reco::Vertex::yError
double yError() const
error on y
Definition: Vertex.h:121
PrimaryVertexValidation::chi2UnbiasedVertex_
float chi2UnbiasedVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:237
reco::Vertex::tracksSize
size_t tracksSize() const
number of tracks
Definition: Vertex.cc:71
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:631
reco::TrackBase::py
double py() const
y coordinate of momentum vector
Definition: TrackBase.h:633
reco::TrackBase::vx
double vx() const
x coordinate of the reference point on track
Definition: TrackBase.h:657
GlobalVector
Global3DVector GlobalVector
Definition: GlobalVector.h:10
PrimaryVertexValidation::a_dzEtaBiasResiduals
std::vector< TH1F * > a_dzEtaBiasResiduals
Definition: PrimaryVertexValidation.h:458
SiPixelRecHit
Our base class.
Definition: SiPixelRecHit.h:23
PrimaryVertexValidation::a_IP3DPhiResiduals
std::vector< TH1F * > a_IP3DPhiResiduals
Definition: PrimaryVertexValidation.h:293
PrimaryVertexValidation::n_dzPhiBiasResiduals
std::vector< TH1F * > n_dzPhiBiasResiduals
Definition: PrimaryVertexValidation.h:465
reco::BeamSpot::x0
double x0() const
x coordinate
Definition: BeamSpot.h:64
reco::Track::recHitsEnd
trackingRecHit_iterator recHitsEnd() const
Iterator to last hit on the track.
Definition: Track.h:109
PrimaryVertexValidation::IP3DsigFromMyVertex_
double IP3DsigFromMyVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:253
Generated for CMSSW Reference Manual by doxygen 1.8.11