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