test
CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
PrimaryVertexValidation.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: Alignment/OfflineValidation
4 // Class: PrimaryVertexValidation
5 //
13 //
14 // Original Author: Marco Musich
15 // Created: Tue Mar 02 10:39:34 CDT 2010
16 //
17 
18 // system include files
19 #include <memory>
20 #include <vector>
21 
22 // user include files
23 #include "Alignment/OfflineValidation/plugins/PrimaryVertexValidation.h"
24 
25 // ROOT includes
26 #include "TH1F.h"
27 #include "TH2F.h"
28 #include "TF1.h"
29 #include "TVector3.h"
30 #include "TFile.h"
31 #include "TMath.h"
32 #include "TROOT.h"
33 #include "TChain.h"
34 #include "TNtuple.h"
35 #include "TMatrixD.h"
36 #include "TVectorD.h"
37 
38 // CMSSW includes
39 #include "CommonTools/UtilAlgos/interface/TFileService.h"
40 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
41 #include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
42 #include "DataFormats/TrackReco/interface/Track.h"
43 #include "DataFormats/TrackReco/interface/TrackFwd.h"
44 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
45 #include "FWCore/MessageLogger/interface/MessageLogger.h"
46 #include "FWCore/ParameterSet/interface/ParameterSet.h"
47 #include "FWCore/ServiceRegistry/interface/Service.h"
48 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
49 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
50 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
51 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
52 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
53 #include "MagneticField/Engine/interface/MagneticField.h"
54 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
55 #include "RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZ.h"
56 #include "RecoVertex/PrimaryVertexProducer/interface/GapClusterizerInZ.h"
57 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
58 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
59 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
60 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateClosestToPoint.h"
61 
62 const int kBPIX = PixelSubdetector::PixelBarrel;
63 const int kFPIX = PixelSubdetector::PixelEndcap;
64 
65 // Constructor
66 PrimaryVertexValidation::PrimaryVertexValidation(const edm::ParameterSet& iConfig):
67  storeNtuple_(iConfig.getParameter<bool>("storeNtuple")),
68  lightNtupleSwitch_(iConfig.getParameter<bool>("isLightNtuple")),
69  useTracksFromRecoVtx_(iConfig.getParameter<bool>("useTracksFromRecoVtx")),
70  vertexZMax_(iConfig.getUntrackedParameter<double>("vertexZMax",99.)),
71  askFirstLayerHit_(iConfig.getParameter<bool>("askFirstLayerHit")),
72  ptOfProbe_(iConfig.getUntrackedParameter<double>("probePt",0.)),
73  etaOfProbe_(iConfig.getUntrackedParameter<double>("probeEta",2.4)),
74  nBins_(iConfig.getUntrackedParameter<int>("numberOfBins",24)),
75  debug_(iConfig.getParameter<bool>("Debug")),
76  runControl_(iConfig.getUntrackedParameter<bool>("runControl",false))
77 {
78 
79  // now do what ever initialization is needed
80  // initialize phase space boundaries
81 
82  usesResource(TFileService::kSharedResource);
83 
84  std::vector<unsigned int> defaultRuns;
85  defaultRuns.push_back(0);
86  runControlNumbers_ = iConfig.getUntrackedParameter<std::vector<unsigned int> >("runControlNumber",defaultRuns);
87 
88  phiSect_ = (2*TMath::Pi())/nBins_;
89  etaSect_ = 5./nBins_;
90 
91  edm::InputTag TrackCollectionTag_ = iConfig.getParameter<edm::InputTag>("TrackCollectionTag");
92  theTrackCollectionToken = consumes<reco::TrackCollection>(TrackCollectionTag_);
93 
94  edm::InputTag VertexCollectionTag_ = iConfig.getParameter<edm::InputTag>("VertexCollectionTag");
95  theVertexCollectionToken = consumes<reco::VertexCollection>(VertexCollectionTag_);
96 
97  edm::InputTag BeamspotTag_ = edm::InputTag("offlineBeamSpot");
98  theBeamspotToken = consumes<reco::BeamSpot>(BeamspotTag_);
99 
100  // select and configure the track filter
101  theTrackFilter_= new TrackFilterForPVFinding(iConfig.getParameter<edm::ParameterSet>("TkFilterParameters") );
102  // select and configure the track clusterizer
103  std::string clusteringAlgorithm=iConfig.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<std::string>("algorithm");
104  if (clusteringAlgorithm=="gap"){
105  theTrackClusterizer_ = new GapClusterizerInZ(iConfig.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<edm::ParameterSet>("TkGapClusParameters"));
106  }else if(clusteringAlgorithm=="DA"){
107  theTrackClusterizer_ = new DAClusterizerInZ(iConfig.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<edm::ParameterSet>("TkDAClusParameters"));
108  }else{
109  throw VertexException("PrimaryVertexProducerAlgorithm: unknown clustering algorithm: " + clusteringAlgorithm);
110  }
111 }
112 
113 // Destructor
114 PrimaryVertexValidation::~PrimaryVertexValidation()
115 {
116  // do anything here that needs to be done at desctruction time
117  // (e.g. close files, deallocate resources etc.)
118 }
119 
120 
121 //
122 // member functions
123 //
124 
125 // ------------ method called to for each event ------------
126 void
127 PrimaryVertexValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
128 {
129 
130  using namespace std;
131  using namespace reco;
132  using namespace IPTools;
133 
134  if(nBins_!=24){
135  edm::LogInfo("PrimaryVertexValidation")<<"Using: "<<nBins_<<" bins plots";
136  }
137 
138  bool passesRunControl = false;
139 
140  if(runControl_){
141  for(unsigned int j=0;j<runControlNumbers_.size();j++){
142  if(iEvent.eventAuxiliary().run() == runControlNumbers_[j]){
143  if (debug_){
144  edm::LogInfo("PrimaryVertexValidation")<<" run number: "<<iEvent.eventAuxiliary().run()<<" keeping run:"<<runControlNumbers_[j];
145  }
146  passesRunControl = true;
147  break;
148  }
149  }
150  if (!passesRunControl) return;
151  }
152 
153  Nevt_++;
154 
155  //=======================================================
156  // Initialize Root-tuple variables
157  //=======================================================
158 
159  SetVarToZero();
160 
161  //=======================================================
162  // Retrieve the Magnetic Field information
163  //=======================================================
164 
165  edm::ESHandle<MagneticField> theMGField;
166  iSetup.get<IdealMagneticFieldRecord>().get( theMGField );
167 
168  //=======================================================
169  // Retrieve the Tracking Geometry information
170  //=======================================================
171 
172  edm::ESHandle<GlobalTrackingGeometry> theTrackingGeometry;
173  iSetup.get<GlobalTrackingGeometryRecord>().get( theTrackingGeometry );
174 
175  //=======================================================
176  // Retrieve the Transient Track Builder information
177  //=======================================================
178 
179  edm::ESHandle<TransientTrackBuilder> theB_;
180  iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder",theB_);
181  double fBfield_=((*theB_).field()->inTesla(GlobalPoint(0.,0.,0.))).z();
182 
183  //=======================================================
184  // Retrieve the Track information
185  //=======================================================
186 
187  edm::Handle<TrackCollection> trackCollectionHandle;
188  iEvent.getByToken(theTrackCollectionToken, trackCollectionHandle);
189 
190  //=======================================================
191  // Retrieve offline vartex information (only for reco)
192  //=======================================================
193 
194  //edm::Handle<VertexCollection> vertices;
195  edm::Handle<std::vector<Vertex> > vertices;
196 
197  try {
198  iEvent.getByToken(theVertexCollectionToken, vertices);
199  } catch ( cms::Exception& er ) {
200  LogTrace("PrimaryVertexValidation")<<"caught std::exception "<<er.what()<<std::endl;
201  }
202 
203  std::vector<Vertex> vsorted = *(vertices);
204  // sort the vertices by number of tracks in descending order
205  // use chi2 as tiebreaker
206  std::sort( vsorted.begin(), vsorted.end(), PrimaryVertexValidation::vtxSort );
207 
208  // skip events with no PV, this should not happen
209  if( vsorted.size() == 0) return;
210  // skip events failing vertex cut
211  if( fabs(vsorted[0].z()) > vertexZMax_ ) return;
212 
213  if ( vsorted[0].isValid() ) {
214  xOfflineVertex_ = (vsorted)[0].x();
215  yOfflineVertex_ = (vsorted)[0].y();
216  zOfflineVertex_ = (vsorted)[0].z();
217 
218  xErrOfflineVertex_ = (vsorted)[0].xError();
219  yErrOfflineVertex_ = (vsorted)[0].yError();
220  zErrOfflineVertex_ = (vsorted)[0].zError();
221  }
222 
223  h_xOfflineVertex->Fill(xOfflineVertex_);
224  h_yOfflineVertex->Fill(yOfflineVertex_);
225  h_zOfflineVertex->Fill(zOfflineVertex_);
226  h_xErrOfflineVertex->Fill(xErrOfflineVertex_);
227  h_yErrOfflineVertex->Fill(yErrOfflineVertex_);
228  h_zErrOfflineVertex->Fill(zErrOfflineVertex_);
229 
230  unsigned int vertexCollectionSize = vsorted.size();
231  int nvvertex = 0;
232 
233  for (unsigned int i=0; i<vertexCollectionSize; i++) {
234  const Vertex& vertex = vsorted.at(i);
235  if (vertex.isValid()) nvvertex++;
236  }
237 
238  nOfflineVertices_ = nvvertex;
239  h_nOfflineVertices->Fill(nvvertex);
240 
241  if ( vsorted.size() && useTracksFromRecoVtx_ ) {
242 
243  double sumpt = 0;
244  size_t ntracks = 0;
245  double chi2ndf = 0.;
246  double chi2prob = 0.;
247 
248  if (!vsorted.at(0).isFake()) {
249 
250  Vertex pv = vsorted.at(0);
251 
252  ntracks = pv.tracksSize();
253  chi2ndf = pv.normalizedChi2();
254  chi2prob = TMath::Prob(pv.chi2(),(int)pv.ndof());
255 
256  h_recoVtxNtracks_->Fill(ntracks);
257  h_recoVtxChi2ndf_->Fill(chi2ndf);
258  h_recoVtxChi2Prob_->Fill(chi2prob);
259 
260  for (Vertex::trackRef_iterator itrk = pv.tracks_begin();itrk != pv.tracks_end(); ++itrk) {
261  double pt = (**itrk).pt();
262  sumpt += pt*pt;
263 
264  const math::XYZPoint myVertex(pv.position().x(),pv.position().y(),pv.position().z());
265 
266  double dxyRes = (**itrk).dxy(myVertex);
267  double dzRes = (**itrk).dz(myVertex);
268 
269  double dxy_err = (**itrk).dxyError();
270  double dz_err = (**itrk).dzError();
271 
272  float trackphi = ((**itrk).phi())*(180/TMath::Pi());
273  float tracketa = (**itrk).eta();
274 
275  for(int i=0; i<nBins_; i++){
276 
277  float phiF = (-TMath::Pi()+i*phiSect_)*(180/TMath::Pi());
278  float phiL = (-TMath::Pi()+(i+1)*phiSect_)*(180/TMath::Pi());
279 
280  float etaF=-2.5+i*etaSect_;
281  float etaL=-2.5+(i+1)*etaSect_;
282 
283  if(tracketa >= etaF && tracketa < etaL ){
284 
285  a_dxyEtaBiasResiduals[i]->Fill(dxyRes*cmToum);
286  a_dzEtaBiasResiduals[i]->Fill(dzRes*cmToum);
287  n_dxyEtaBiasResiduals[i]->Fill((dxyRes)/dxy_err);
288  n_dzEtaBiasResiduals[i]->Fill((dzRes)/dz_err);
289 
290  }
291 
292  if(trackphi >= phiF && trackphi < phiL ){
293 
294  a_dxyPhiBiasResiduals[i]->Fill(dxyRes*cmToum);
295  a_dzPhiBiasResiduals[i]->Fill(dzRes*cmToum);
296  n_dxyPhiBiasResiduals[i]->Fill((dxyRes)/dxy_err);
297  n_dzPhiBiasResiduals[i]->Fill((dzRes)/dz_err);
298 
299  for(int j=0; j<nBins_; j++){
300 
301  float etaJ=-2.5+j*etaSect_;
302  float etaK=-2.5+(j+1)*etaSect_;
303 
304  if(tracketa >= etaJ && tracketa < etaK ){
305 
306  a_dxyBiasResidualsMap[i][j]->Fill(dxyRes*cmToum);
307  a_dzBiasResidualsMap[i][j]->Fill(dzRes*cmToum);
308 
309  n_dxyBiasResidualsMap[i][j]->Fill((dxyRes)/dxy_err);
310  n_dzBiasResidualsMap[i][j]->Fill((dzRes)/dz_err);
311 
312  }
313  }
314  }
315  }
316  }
317 
318  h_recoVtxSumPt_->Fill(sumpt);
319 
320  }
321  }
322 
323  //=======================================================
324  // Retrieve Beamspot information
325  //=======================================================
326 
327  BeamSpot beamSpot;
328  edm::Handle<BeamSpot> beamSpotHandle;
329  iEvent.getByToken(theBeamspotToken, beamSpotHandle);
330 
331  if ( beamSpotHandle.isValid() ) {
332  beamSpot = *beamSpotHandle;
333  BSx0_ = beamSpot.x0();
334  BSy0_ = beamSpot.y0();
335  BSz0_ = beamSpot.z0();
336  Beamsigmaz_ = beamSpot.sigmaZ();
337  Beamdxdz_ = beamSpot.dxdz();
338  BeamWidthX_ = beamSpot.BeamWidthX();
339  BeamWidthY_ = beamSpot.BeamWidthY();
340 
341  wxy2_=TMath::Power(BeamWidthX_,2)+TMath::Power(BeamWidthY_,2);
342 
343  } else {
344  edm::LogWarning("PrimaryVertexValidation")<<"No BeamSpot found!";
345  }
346 
347  h_BSx0->Fill(BSx0_);
348  h_BSy0->Fill(BSy0_);
349  h_BSz0->Fill(BSz0_);
350  h_Beamsigmaz->Fill(Beamsigmaz_);
351  h_BeamWidthX->Fill(BeamWidthX_);
352  h_BeamWidthY->Fill(BeamWidthY_);
353 
354  if(debug_)
355  edm::LogInfo("PrimaryVertexValidation")<<"Beamspot x:" <<BSx0_<<" y:"<<BSy0_<<" z:"<<BSz0_;
356 
357  //=======================================================
358  // Starts here ananlysis
359  //=======================================================
360 
361  RunNumber_=iEvent.eventAuxiliary().run();
362  h_runNumber->Fill(RunNumber_);
363  LuminosityBlockNumber_=iEvent.eventAuxiliary().luminosityBlock();
364  EventNumber_=iEvent.eventAuxiliary().id().event();
365 
366  if(debug_)
367  edm::LogInfo("PrimaryVertexValidation")<<" looping over "<<trackCollectionHandle->size()<< "tracks";
368 
369  h_nTracks->Fill(trackCollectionHandle->size());
370 
371  //======================================================
372  // Interface RECO tracks to vertex reconstruction
373  //======================================================
374 
375  std::vector<TransientTrack> t_tks;
376  unsigned int k = 0;
377  for(TrackCollection::const_iterator track = trackCollectionHandle->begin(); track!= trackCollectionHandle->end(); ++track, ++k){
378 
379  TransientTrack tt = theB_->build(&(*track));
380  tt.setBeamSpot(beamSpot);
381  t_tks.push_back(tt);
382 
383  }
384 
385  if(debug_) {
386  edm::LogInfo("PrimaryVertexValidation") << "Found: " << t_tks.size() << " reconstructed tracks";
387  }
388 
389  //======================================================
390  // select the tracks
391  //======================================================
392 
393  std::vector<TransientTrack> seltks = theTrackFilter_->select(t_tks);
394 
395  //======================================================
396  // clusterize tracks in Z
397  //======================================================
398 
399  vector< vector<TransientTrack> > clusters = theTrackClusterizer_->clusterize(seltks);
400 
401  if (debug_){
402  edm::LogInfo("PrimaryVertexValidation")<<" looping over: "<< clusters.size() << " clusters from " << t_tks.size() << " selected tracks";
403  }
404 
405  nClus_=clusters.size();
406  h_nClus->Fill(nClus_);
407 
408  //======================================================
409  // Starts loop on clusters
410  //======================================================
411 
412  for (vector< vector<TransientTrack> >::const_iterator iclus = clusters.begin(); iclus != clusters.end(); iclus++) {
413 
414  nTracksPerClus_=0;
415 
416  unsigned int i = 0;
417  for(vector<TransientTrack>::const_iterator theTTrack = iclus->begin(); theTTrack!= iclus->end(); ++theTTrack, ++i)
418  {
419  if ( nTracks_ >= nMaxtracks_ ) {
420  edm::LogError("PrimaryVertexValidation")<<" Warning - Number of tracks: " << nTracks_ << " , greater than " << nMaxtracks_;
421  continue;
422  }
423 
424  const Track & theTrack = theTTrack->track();
425 
426  pt_[nTracks_] = theTrack.pt();
427  p_[nTracks_] = theTrack.p();
428  nhits_[nTracks_] = theTrack.numberOfValidHits();
429  eta_[nTracks_] = theTrack.eta();
430  theta_[nTracks_] = theTrack.theta();
431  phi_[nTracks_] = theTrack.phi();
432  chi2_[nTracks_] = theTrack.chi2();
433  chi2ndof_[nTracks_] = theTrack.normalizedChi2();
434  charge_[nTracks_] = theTrack.charge();
435  qoverp_[nTracks_] = theTrack.qoverp();
436  dz_[nTracks_] = theTrack.dz();
437  dxy_[nTracks_] = theTrack.dxy();
438 
439  TrackBase::TrackQuality _trackQuality = TrackBase::qualityByName("highPurity");
440  isHighPurity_[nTracks_] = theTrack.quality(_trackQuality);
441 
442  math::XYZPoint point(BSx0_,BSy0_,BSz0_);
443  dxyBs_[nTracks_] = theTrack.dxy(point);
444  dzBs_[nTracks_] = theTrack.dz(point);
445 
446  xPCA_[nTracks_] = theTrack.vertex().x();
447  yPCA_[nTracks_] = theTrack.vertex().y();
448  zPCA_[nTracks_] = theTrack.vertex().z();
449 
450  //=======================================================
451  // Retrieve rechit information
452  //=======================================================
453 
454  int nRecHit1D=0;
455  int nRecHit2D=0;
456  int nhitinTIB=0;
457  int nhitinTOB=0;
458  int nhitinTID=0;
459  int nhitinTEC=0;
460  int nhitinBPIX=0;
461  int nhitinFPIX=0;
462 
463  for (trackingRecHit_iterator iHit = theTTrack->recHitsBegin(); iHit != theTTrack->recHitsEnd(); ++iHit) {
464  if((*iHit)->isValid()) {
465 
466  if (this->isHit2D(**iHit)) {++nRecHit2D;}
467  else {++nRecHit1D; }
468 
469  int type =(*iHit)->geographicalId().subdetId();
470 
471  if(type==int(StripSubdetector::TIB)){++nhitinTIB;}
472  if(type==int(StripSubdetector::TOB)){++nhitinTOB;}
473  if(type==int(StripSubdetector::TID)){++nhitinTID;}
474  if(type==int(StripSubdetector::TEC)){++nhitinTEC;}
475  if(type==int( kBPIX)){++nhitinBPIX;}
476  if(type==int( kFPIX)){++nhitinFPIX;}
477  }
478  }
479 
480  nhits1D_[nTracks_] = nRecHit1D;
481  nhits2D_[nTracks_] = nRecHit2D;
482  nhitsBPIX_[nTracks_] = nhitinBPIX;
483  nhitsFPIX_[nTracks_] = nhitinFPIX;
484  nhitsTIB_[nTracks_] = nhitinTIB;
485  nhitsTID_[nTracks_] = nhitinTID;
486  nhitsTOB_[nTracks_] = nhitinTOB;
487  nhitsTEC_[nTracks_] = nhitinTEC;
488 
489  //=======================================================
490  // Good tracks for vertexing selection
491  //=======================================================
492 
493  bool pass = true;
494  if(askFirstLayerHit_) pass = this->hasFirstLayerPixelHits((*theTTrack));
495  if (pass && (theTrack.pt() >=ptOfProbe_) && fabs(theTrack.eta()) <= etaOfProbe_){
496  isGoodTrack_[nTracks_]=1;
497  }
498 
499  //=======================================================
500  // Fit unbiased vertex
501  //=======================================================
502 
503  vector<TransientTrack> theFinalTracks;
504  theFinalTracks.clear();
505 
506  for(vector<TransientTrack>::const_iterator tk = iclus->begin(); tk!= iclus->end(); ++tk){
507 
508  pass = this->hasFirstLayerPixelHits((*tk));
509  if (pass){
510  if( tk == theTTrack ) continue;
511  else {
512  theFinalTracks.push_back((*tk));
513  }
514  }
515  }
516 
517  if(theFinalTracks.size() > 1){
518 
519  if(debug_)
520  edm::LogInfo("PrimaryVertexValidation")<<"Transient Track Collection size: "<<theFinalTracks.size();
521  try{
522 
523  VertexFitter<5>* theFitter = new AdaptiveVertexFitter;
524  TransientVertex theFittedVertex = theFitter->vertex(theFinalTracks);
525 
526  //AdaptiveVertexFitter* theFitter = new AdaptiveVertexFitter;
527  //TransientVertex theFittedVertex = theFitter->vertex(theFinalTracks,beamSpot); // if you want the beam constraint
528 
529  double totalTrackWeights=0;
530  if(theFittedVertex.isValid ()){
531 
532 
533  if(theFittedVertex.hasTrackWeight()){
534  for(size_t rtracks= 0; rtracks < theFinalTracks.size(); rtracks++){
535  sumOfWeightsUnbiasedVertex_[nTracks_] += theFittedVertex.trackWeight(theFinalTracks[rtracks]);
536  totalTrackWeights+= theFittedVertex.trackWeight(theFinalTracks[rtracks]);
537  h_fitVtxTrackWeights_->Fill(theFittedVertex.trackWeight(theFinalTracks[rtracks]));
538  }
539  }
540 
541  h_fitVtxTrackAverageWeight_->Fill(totalTrackWeights/theFinalTracks.size());
542 
543  const math::XYZPoint theRecoVertex(xOfflineVertex_,yOfflineVertex_,zOfflineVertex_);
544  const math::XYZPoint myVertex(theFittedVertex.position().x(),theFittedVertex.position().y(),theFittedVertex.position().z());
545 
546  const Vertex vertex = theFittedVertex;
547  fillTrackHistos(hDA,"all",&(*theTTrack),vertex,beamSpot,fBfield_);
548 
549  hasRecVertex_[nTracks_] = 1;
550  xUnbiasedVertex_[nTracks_] = theFittedVertex.position().x();
551  yUnbiasedVertex_[nTracks_] = theFittedVertex.position().y();
552  zUnbiasedVertex_[nTracks_] = theFittedVertex.position().z();
553 
554  chi2normUnbiasedVertex_[nTracks_] = theFittedVertex.normalisedChiSquared();
555  chi2UnbiasedVertex_[nTracks_] = theFittedVertex.totalChiSquared();
556  DOFUnbiasedVertex_[nTracks_] = theFittedVertex.degreesOfFreedom();
557  chi2ProbUnbiasedVertex_[nTracks_] = TMath::Prob(theFittedVertex.totalChiSquared(),(int)theFittedVertex.degreesOfFreedom());
558  tracksUsedForVertexing_[nTracks_] = theFinalTracks.size();
559 
560  h_fitVtxNtracks_->Fill(theFinalTracks.size());
561  h_fitVtxChi2_->Fill(theFittedVertex.totalChiSquared());
562  h_fitVtxNdof_->Fill(theFittedVertex.degreesOfFreedom());
563  h_fitVtxChi2ndf_->Fill(theFittedVertex.normalisedChiSquared());
564  h_fitVtxChi2Prob_->Fill(TMath::Prob(theFittedVertex.totalChiSquared(),(int)theFittedVertex.degreesOfFreedom()));
565 
566  // from my Vertex
567  double dxyFromMyVertex = theTrack.dxy(myVertex);
568  double dzFromMyVertex = theTrack.dz(myVertex);
569 
570  double dz_err = hypot(theTrack.dzError(),theFittedVertex.positionError().czz());
571 
572  // PV2D
573  std::pair<bool,Measurement1D> s_ip2dpv = signedTransverseImpactParameter(*theTTrack,
574  GlobalVector(theTrack.px(),
575  theTrack.py(),
576  theTrack.pz()),
577  theFittedVertex);
578 
579  double s_ip2dpv_corr = s_ip2dpv.second.value();
580  double s_ip2dpv_err = s_ip2dpv.second.error();
581 
582  // PV3D
583  std::pair<bool, Measurement1D> s_ip3dpv = signedImpactParameter3D(*theTTrack,
584  GlobalVector(theTrack.px(),
585  theTrack.py(),
586  theTrack.pz()),
587  theFittedVertex);
588 
589  double s_ip3dpv_corr = s_ip3dpv.second.value();
590  double s_ip3dpv_err = s_ip3dpv.second.error();
591 
592  // PV3D absolute
593  std::pair<bool,Measurement1D> ip3dpv = absoluteImpactParameter3D(*theTTrack,theFittedVertex);
594  double ip3d_corr = ip3dpv.second.value();
595  double ip3d_err = ip3dpv.second.error();
596 
597  // with respect to any specified vertex, such as primary vertex
598  GlobalPoint vert(theFittedVertex.position().x(),theFittedVertex.position().y(),theFittedVertex.position().z());
599  TrajectoryStateClosestToPoint traj = (*theTTrack).trajectoryStateClosestToPoint(vert);
600 
601  double d0 = traj.perigeeParameters().transverseImpactParameter();
602  //double d0_error = traj.perigeeError().transverseImpactParameterError();
603  double z0 = traj.perigeeParameters().longitudinalImpactParameter();
604  double z0_error = traj.perigeeError().longitudinalImpactParameterError();
605 
606  // define IPs
607 
608  dxyFromMyVertex_[nTracks_] = dxyFromMyVertex;
609  dxyErrorFromMyVertex_[nTracks_] = s_ip2dpv_err;
610  IPTsigFromMyVertex_[nTracks_] = dxyFromMyVertex/s_ip2dpv_err;
611 
612  dzFromMyVertex_[nTracks_] = dzFromMyVertex;
613  dzErrorFromMyVertex_[nTracks_] = dz_err;
614  IPLsigFromMyVertex_[nTracks_] = dzFromMyVertex/dz_err;
615 
616  d3DFromMyVertex_[nTracks_] = ip3d_corr;
617  d3DErrorFromMyVertex_[nTracks_] = ip3d_err;
618  IP3DsigFromMyVertex_[nTracks_] = (ip3d_corr/ip3d_err);
619 
620  // fill directly the histograms of residuals
621 
622  float trackphi = (theTrack.phi())*(180/TMath::Pi());
623  float tracketa = theTrack.eta();
624  float trackpt = theTrack.pt();
625 
626  // checks on the probe track quality
627  if(trackpt >= ptOfProbe_ && fabs(tracketa)<= etaOfProbe_){
628 
629  fillTrackHistos(hDA,"sel",&(*theTTrack),vertex,beamSpot,fBfield_);
630 
631  // probe checks
632  h_probePt_->Fill(theTrack.pt());
633  h_probeP_->Fill(theTrack.p());
634  h_probeEta_->Fill(theTrack.eta());
635  h_probePhi_->Fill(theTrack.phi());
636  h2_probeEtaPhi_->Fill(theTrack.eta(),theTrack.phi());
637  h2_probeEtaPt_->Fill(theTrack.eta(),theTrack.pt());
638 
639  h_probeChi2_->Fill(theTrack.chi2());
640  h_probeNormChi2_->Fill(theTrack.normalizedChi2());
641  h_probeCharge_->Fill(theTrack.charge());
642  h_probeQoverP_->Fill(theTrack.qoverp());
643  h_probeHits_->Fill(theTrack.numberOfValidHits());
644  h_probeHits1D_->Fill(nRecHit1D);
645  h_probeHits2D_->Fill(nRecHit2D);
646  h_probeHitsInTIB_->Fill(nhitinBPIX);
647  h_probeHitsInTOB_->Fill(nhitinFPIX);
648  h_probeHitsInTID_->Fill(nhitinTIB);
649  h_probeHitsInTEC_->Fill(nhitinTID);
650  h_probeHitsInBPIX_->Fill(nhitinTOB);
651  h_probeHitsInFPIX_->Fill(nhitinTEC);
652 
653  float dxyRecoV = theTrack.dz(theRecoVertex);
654  float dzRecoV = theTrack.dxy(theRecoVertex);
655  float dxysigmaRecoV = TMath::Sqrt(theTrack.d0Error()*theTrack.d0Error()+xErrOfflineVertex_*yErrOfflineVertex_);
656  float dzsigmaRecoV = TMath::Sqrt(theTrack.dzError()*theTrack.dzError()+zErrOfflineVertex_*zErrOfflineVertex_);
657 
658  double zTrack=(theTTrack->stateAtBeamLine().trackStateAtPCA()).position().z();
659  double zVertex=theFittedVertex.position().z();
660  double tantheta=tan((theTTrack->stateAtBeamLine().trackStateAtPCA()).momentum().theta());
661 
662  double dz2= pow(theTrack.dzError(),2)+wxy2_/pow(tantheta,2);
663  double restrkz = zTrack-zVertex;
664  double pulltrkz = (zTrack-zVertex)/TMath::Sqrt(dz2);
665 
666  h_probedxyRecoV_->Fill(dxyRecoV);
667  h_probedzRecoV_->Fill(dzRecoV);
668 
669  h_probedzRefitV_->Fill(dxyFromMyVertex);
670  h_probedxyRefitV_->Fill(dzFromMyVertex);
671 
672  h_probed0RefitV_->Fill(d0);
673  h_probez0RefitV_->Fill(z0);
674 
675  h_probesignIP2DRefitV_->Fill(s_ip2dpv_corr);
676  h_probed3DRefitV_->Fill(ip3d_corr);
677  h_probereszRefitV_->Fill(restrkz);
678 
679  h_probeRecoVSigZ_->Fill(dzRecoV/dzsigmaRecoV);
680  h_probeRecoVSigXY_->Fill(dxyRecoV/dxysigmaRecoV);
681  h_probeRefitVSigZ_->Fill(dzFromMyVertex/dz_err);
682  h_probeRefitVSigXY_->Fill(dxyFromMyVertex/s_ip2dpv_err);
683  h_probeRefitVSig3D_->Fill(ip3d_corr/ip3d_err);
684  h_probeRefitVLogSig3D_->Fill(log10(ip3d_corr/ip3d_err));
685  h_probeRefitVSigResZ_->Fill(pulltrkz);
686 
687  a_dxyVsPhi->Fill(trackphi,dxyFromMyVertex*cmToum);
688  a_dzVsPhi->Fill(trackphi,z0*cmToum);
689  n_dxyVsPhi->Fill(trackphi,dxyFromMyVertex/s_ip2dpv_err);
690  n_dzVsPhi->Fill(trackphi,z0/z0_error);
691 
692  a_dxyVsEta->Fill(tracketa,dxyFromMyVertex*cmToum);
693  a_dzVsEta->Fill(tracketa,z0*cmToum);
694  n_dxyVsEta->Fill(tracketa,dxyFromMyVertex/s_ip2dpv_err);
695  n_dzVsEta->Fill(tracketa,z0/z0_error);
696 
697  // filling the binned distributions
698  for(int i=0; i<nBins_; i++){
699 
700  float phiF = (-TMath::Pi()+i*phiSect_)*(180/TMath::Pi());
701  float phiL = (-TMath::Pi()+(i+1)*phiSect_)*(180/TMath::Pi());
702 
703  float etaF=-2.5+i*etaSect_;
704  float etaL=-2.5+(i+1)*etaSect_;
705 
706  if(tracketa >= etaF && tracketa < etaL ){
707 
708  a_dxyEtaResiduals[i]->Fill(dxyFromMyVertex*cmToum);
709  a_dzEtaResiduals[i]->Fill(dzFromMyVertex*cmToum);
710  n_dxyEtaResiduals[i]->Fill(dxyFromMyVertex/s_ip2dpv_err);
711  n_dzEtaResiduals[i]->Fill(dzFromMyVertex/dz_err);
712  a_IP2DEtaResiduals[i]->Fill(s_ip2dpv_corr*cmToum);
713  n_IP2DEtaResiduals[i]->Fill(s_ip2dpv_corr/s_ip2dpv_err);
714  a_reszEtaResiduals[i]->Fill(restrkz*cmToum);
715  n_reszEtaResiduals[i]->Fill(pulltrkz);
716  a_d3DEtaResiduals[i]->Fill(ip3d_corr*cmToum);
717  n_d3DEtaResiduals[i]->Fill(ip3d_corr/ip3d_err);
718  a_IP3DEtaResiduals[i]->Fill(s_ip3dpv_corr*cmToum);
719  n_IP3DEtaResiduals[i]->Fill(s_ip3dpv_corr/s_ip3dpv_err);
720 
721  }
722 
723  if(trackphi >= phiF && trackphi < phiL ){
724  a_dxyPhiResiduals[i]->Fill(dxyFromMyVertex*cmToum);
725  a_dzPhiResiduals[i]->Fill(dzFromMyVertex*cmToum);
726  n_dxyPhiResiduals[i]->Fill(dxyFromMyVertex/s_ip2dpv_err);
727  n_dzPhiResiduals[i]->Fill(dzFromMyVertex/dz_err);
728  a_IP2DPhiResiduals[i]->Fill(s_ip2dpv_corr*cmToum);
729  n_IP2DPhiResiduals[i]->Fill(s_ip2dpv_corr/s_ip2dpv_err);
730  a_reszPhiResiduals[i]->Fill(restrkz*cmToum);
731  n_reszPhiResiduals[i]->Fill(pulltrkz);
732  a_d3DPhiResiduals[i]->Fill(ip3d_corr*cmToum);
733  n_d3DPhiResiduals[i]->Fill(ip3d_corr/ip3d_err);
734  a_IP3DPhiResiduals[i]->Fill(s_ip3dpv_corr*cmToum);
735  n_IP3DPhiResiduals[i]->Fill(s_ip3dpv_corr/s_ip3dpv_err);
736 
737  for(int j=0; j<nBins_; j++){
738 
739  float etaJ=-2.5+j*etaSect_;
740  float etaK=-2.5+(j+1)*etaSect_;
741 
742  if(tracketa >= etaJ && tracketa < etaK ){
743  a_dxyResidualsMap[i][j]->Fill(dxyFromMyVertex*cmToum);
744  a_dzResidualsMap[i][j]->Fill(dzFromMyVertex*cmToum);
745  n_dxyResidualsMap[i][j]->Fill(dxyFromMyVertex/s_ip2dpv_err);
746  n_dzResidualsMap[i][j]->Fill(dzFromMyVertex/dz_err);
747  a_d3DResidualsMap[i][j]->Fill(ip3d_corr*cmToum);
748  n_d3DResidualsMap[i][j]->Fill(ip3d_corr/ip3d_err);
749 
750  }
751  }
752  }
753  }
754  }
755 
756  if(debug_){
757  edm::LogInfo("PrimaryVertexValidation")<<" myVertex.x()= "<<myVertex.x()<<"\n"
758  <<" myVertex.y()= "<<myVertex.y()<<" \n"
759  <<" myVertex.z()= "<<myVertex.z()<<" \n"
760  <<" theTrack.dz(myVertex)= "<<theTrack.dz(myVertex)<<" \n"
761  <<" zPCA -myVertex.z() = "<<(theTrack.vertex().z() -myVertex.z());
762 
763  }// ends if debug_
764  } // ends if the fitted vertex is Valid
765 
766  delete theFitter;
767 
768  } catch ( cms::Exception& er ) {
769  LogTrace("PrimaryVertexValidation")<<"caught std::exception "<<er.what()<<std::endl;
770  }
771 
772  } //ends if theFinalTracks.size() > 2
773 
774  else {
775  if(debug_)
776  edm::LogInfo("PrimaryVertexValidation")<<"Not enough tracks to make a vertex. Returns no vertex info";
777  }
778 
779  ++nTracks_;
780  ++nTracksPerClus_;
781 
782  if(debug_)
783  edm::LogInfo("PrimaryVertexValidation")<<"Track "<<i<<" : pT = "<<theTrack.pt();
784 
785  }// for loop on tracks
786  } // for loop on track clusters
787 
788  // Fill the TTree if needed
789 
790  if(storeNtuple_){
791  rootTree_->Fill();
792  }
793 
794 
795 }
796 
797 // ------------ method called to discriminate 1D from 2D hits ------------
798 bool PrimaryVertexValidation::isHit2D(const TrackingRecHit &hit) const
799 {
800  if (hit.dimension() < 2) {
801  return false; // some (muon...) stuff really has RecHit1D
802  } else {
803  const DetId detId(hit.geographicalId());
804  if (detId.det() == DetId::Tracker) {
805  if (detId.subdetId() == kBPIX || detId.subdetId() == kFPIX) {
806  return true; // pixel is always 2D
807  } else { // should be SiStrip now
808  if (dynamic_cast<const SiStripRecHit2D*>(&hit)) return false; // normal hit
809  else if (dynamic_cast<const SiStripMatchedRecHit2D*>(&hit)) return true; // matched is 2D
810  else if (dynamic_cast<const ProjectedSiStripRecHit2D*>(&hit)) return false; // crazy hit...
811  else {
812  edm::LogError("UnkownType") << "@SUB=AlignmentTrackSelector::isHit2D"
813  << "Tracker hit not in pixel and neither SiStripRecHit2D nor "
814  << "SiStripMatchedRecHit2D nor ProjectedSiStripRecHit2D.";
815  return false;
816  }
817  }
818  } else { // not tracker??
819  edm::LogWarning("DetectorMismatch") << "@SUB=AlignmentTrackSelector::isHit2D"
820  << "Hit not in tracker with 'official' dimension >=2.";
821  return true; // dimension() >= 2 so accept that...
822  }
823  }
824  // never reached...
825 }
826 
827 // ------------ method to check the presence of pixel hits ------------
828 bool PrimaryVertexValidation::hasFirstLayerPixelHits(const reco::TransientTrack track)
829 {
830  using namespace reco;
831  const HitPattern& p = track.hitPattern();
832  for (int i=0; i<p.numberOfHits(HitPattern::TRACK_HITS); i++) {
833  uint32_t pattern = p.getHitPattern(HitPattern::TRACK_HITS, i);
834  if (p.pixelBarrelHitFilter(pattern) || p.pixelEndcapHitFilter(pattern) ) {
835  if (p.getLayer(pattern) == 1) {
836  if (p.validHitFilter(pattern)) {
837  return true;
838  }
839  }
840  }
841  }
842  return false;
843 }
844 
845 // ------------ method called once each job before begining the event loop ------------
846 void PrimaryVertexValidation::beginJob()
847 {
848  edm::LogInfo("PrimaryVertexValidation")
849  <<"######################################\n"
850  <<"Begin Job \n"
851  <<"######################################";
852 
853  // Define TTree for output
854  Nevt_ = 0;
855 
856  // rootFile_ = new TFile(filename_.c_str(),"recreate");
857  edm::Service<TFileService> fs;
858  rootTree_ = fs->make<TTree>("tree","PV Validation tree");
859 
860  // Track Paramters
861 
862  if(lightNtupleSwitch_){
863 
864  rootTree_->Branch("EventNumber",&EventNumber_,"EventNumber/i");
865  rootTree_->Branch("RunNumber",&RunNumber_,"RunNumber/i");
866  rootTree_->Branch("LuminosityBlockNumber",&LuminosityBlockNumber_,"LuminosityBlockNumber/i");
867  rootTree_->Branch("nOfflineVertices",&nOfflineVertices_,"nOfflineVertices/I");
868  rootTree_->Branch("nTracks",&nTracks_,"nTracks/I");
869  rootTree_->Branch("phi",&phi_,"phi[nTracks]/D");
870  rootTree_->Branch("eta",&eta_,"eta[nTracks]/D");
871  rootTree_->Branch("pt",&pt_,"pt[nTracks]/D");
872  rootTree_->Branch("dxyFromMyVertex",&dxyFromMyVertex_,"dxyFromMyVertex[nTracks]/D");
873  rootTree_->Branch("dzFromMyVertex",&dzFromMyVertex_,"dzFromMyVertex[nTracks]/D");
874  rootTree_->Branch("d3DFromMyVertex",&d3DFromMyVertex_,"d3DFromMyVertex[nTracks]/D");
875  rootTree_->Branch("IPTsigFromMyVertex",&IPTsigFromMyVertex_,"IPTsigFromMyVertex_[nTracks]/D");
876  rootTree_->Branch("IPLsigFromMyVertex",&IPLsigFromMyVertex_,"IPLsigFromMyVertex_[nTracks]/D");
877  rootTree_->Branch("IP3DsigFromMyVertex",&IP3DsigFromMyVertex_,"IP3DsigFromMyVertex_[nTracks]/D");
878  rootTree_->Branch("hasRecVertex",&hasRecVertex_,"hasRecVertex[nTracks]/I");
879  rootTree_->Branch("isGoodTrack",&isGoodTrack_,"isGoodTrack[nTracks]/I");
880  rootTree_->Branch("isHighPurity",&isHighPurity_,"isHighPurity_[nTracks]/I");
881 
882  } else {
883 
884  rootTree_->Branch("nTracks",&nTracks_,"nTracks/I");
885  rootTree_->Branch("nTracksPerClus",&nTracksPerClus_,"nTracksPerClus/I");
886  rootTree_->Branch("nClus",&nClus_,"nClus/I");
887  rootTree_->Branch("xOfflineVertex",&xOfflineVertex_,"xOfflineVertex/D");
888  rootTree_->Branch("yOfflineVertex",&yOfflineVertex_,"yOfflineVertex/D");
889  rootTree_->Branch("zOfflineVertex",&zOfflineVertex_,"zOfflineVertex/D");
890  rootTree_->Branch("BSx0",&BSx0_,"BSx0/D");
891  rootTree_->Branch("BSy0",&BSy0_,"BSy0/D");
892  rootTree_->Branch("BSz0",&BSz0_,"BSz0/D");
893  rootTree_->Branch("Beamsigmaz",&Beamsigmaz_,"Beamsigmaz/D");
894  rootTree_->Branch("Beamdxdz",&Beamdxdz_,"Beamdxdz/D");
895  rootTree_->Branch("BeamWidthX",&BeamWidthX_,"BeamWidthX/D");
896  rootTree_->Branch("BeamWidthY",&BeamWidthY_,"BeamWidthY/D");
897  rootTree_->Branch("pt",&pt_,"pt[nTracks]/D");
898  rootTree_->Branch("p",&p_,"p[nTracks]/D");
899  rootTree_->Branch("nhits",&nhits_,"nhits[nTracks]/I");
900  rootTree_->Branch("nhits1D",&nhits1D_,"nhits1D[nTracks]/I");
901  rootTree_->Branch("nhits2D",&nhits2D_,"nhits2D[nTracks]/I");
902  rootTree_->Branch("nhitsBPIX",&nhitsBPIX_,"nhitsBPIX[nTracks]/I");
903  rootTree_->Branch("nhitsFPIX",&nhitsFPIX_,"nhitsFPIX[nTracks]/I");
904  rootTree_->Branch("nhitsTIB",&nhitsTIB_,"nhitsTIB[nTracks]/I");
905  rootTree_->Branch("nhitsTID",&nhitsTID_,"nhitsTID[nTracks]/I");
906  rootTree_->Branch("nhitsTOB",&nhitsTOB_,"nhitsTOB[nTracks]/I");
907  rootTree_->Branch("nhitsTEC",&nhitsTEC_,"nhitsTEC[nTracks]/I");
908  rootTree_->Branch("eta",&eta_,"eta[nTracks]/D");
909  rootTree_->Branch("theta",&theta_,"theta[nTracks]/D");
910  rootTree_->Branch("phi",&phi_,"phi[nTracks]/D");
911  rootTree_->Branch("chi2",&chi2_,"chi2[nTracks]/D");
912  rootTree_->Branch("chi2ndof",&chi2ndof_,"chi2ndof[nTracks]/D");
913  rootTree_->Branch("charge",&charge_,"charge[nTracks]/I");
914  rootTree_->Branch("qoverp",&qoverp_,"qoverp[nTracks]/D");
915  rootTree_->Branch("dz",&dz_,"dz[nTracks]/D");
916  rootTree_->Branch("dxy",&dxy_,"dxy[nTracks]/D");
917  rootTree_->Branch("dzBs",&dzBs_,"dzBs[nTracks]/D");
918  rootTree_->Branch("dxyBs",&dxyBs_,"dxyBs[nTracks]/D");
919  rootTree_->Branch("xPCA",&xPCA_,"xPCA[nTracks]/D");
920  rootTree_->Branch("yPCA",&yPCA_,"yPCA[nTracks]/D");
921  rootTree_->Branch("zPCA",&zPCA_,"zPCA[nTracks]/D");
922  rootTree_->Branch("xUnbiasedVertex",&xUnbiasedVertex_,"xUnbiasedVertex[nTracks]/D");
923  rootTree_->Branch("yUnbiasedVertex",&yUnbiasedVertex_,"yUnbiasedVertex[nTracks]/D");
924  rootTree_->Branch("zUnbiasedVertex",&zUnbiasedVertex_,"zUnbiasedVertex[nTracks]/D");
925  rootTree_->Branch("chi2normUnbiasedVertex",&chi2normUnbiasedVertex_,"chi2normUnbiasedVertex[nTracks]/F");
926  rootTree_->Branch("chi2UnbiasedVertex",&chi2UnbiasedVertex_,"chi2UnbiasedVertex[nTracks]/F");
927  rootTree_->Branch("DOFUnbiasedVertex",&DOFUnbiasedVertex_," DOFUnbiasedVertex[nTracks]/F");
928  rootTree_->Branch("chi2ProbUnbiasedVertex",&chi2ProbUnbiasedVertex_,"chi2ProbUnbiasedVertex[nTracks]/F");
929  rootTree_->Branch("sumOfWeightsUnbiasedVertex",&sumOfWeightsUnbiasedVertex_,"sumOfWeightsUnbiasedVertex[nTracks]/F");
930  rootTree_->Branch("tracksUsedForVertexing",&tracksUsedForVertexing_,"tracksUsedForVertexing[nTracks]/I");
931  rootTree_->Branch("dxyFromMyVertex",&dxyFromMyVertex_,"dxyFromMyVertex[nTracks]/D");
932  rootTree_->Branch("dzFromMyVertex",&dzFromMyVertex_,"dzFromMyVertex[nTracks]/D");
933  rootTree_->Branch("dxyErrorFromMyVertex",&dxyErrorFromMyVertex_,"dxyErrorFromMyVertex_[nTracks]/D");
934  rootTree_->Branch("dzErrorFromMyVertex",&dzErrorFromMyVertex_,"dzErrorFromMyVertex_[nTracks]/D");
935  rootTree_->Branch("IPTsigFromMyVertex",&IPTsigFromMyVertex_,"IPTsigFromMyVertex_[nTracks]/D");
936  rootTree_->Branch("IPLsigFromMyVertex",&IPLsigFromMyVertex_,"IPLsigFromMyVertex_[nTracks]/D");
937  rootTree_->Branch("hasRecVertex",&hasRecVertex_,"hasRecVertex[nTracks]/I");
938  rootTree_->Branch("isGoodTrack",&isGoodTrack_,"isGoodTrack[nTracks]/I");
939  }
940 
941  // event histograms
942  TFileDirectory EventFeatures = fs->mkdir("EventFeatures");
943 
944  TH1F::SetDefaultSumw2(kTRUE);
945 
946  h_nTracks = EventFeatures.make<TH1F>("h_nTracks","number of tracks per event;n_{tracks}/event;n_{events}",300,-0.5,299.5);
947  h_nClus = EventFeatures.make<TH1F>("h_nClus","number of track clusters;n_{clusters}/event;n_{events}",50,-0.5,49.5);
948  h_nOfflineVertices = EventFeatures.make<TH1F>("h_nOfflineVertices","number of offline reconstructed vertices;n_{vertices}/event;n_{events}",50,-0.5,49.5);
949  h_runNumber = EventFeatures.make<TH1F>("h_runNumber","run number;run number;n_{events}",100000,150000.,250000.);
950  h_xOfflineVertex = EventFeatures.make<TH1F>("h_xOfflineVertex","x-coordinate of offline vertex;x_{vertex};n_{events}",100,-0.1,0.1);
951  h_yOfflineVertex = EventFeatures.make<TH1F>("h_yOfflineVertex","y-coordinate of offline vertex;y_{vertex};n_{events}",100,-0.1,0.1);
952  h_zOfflineVertex = EventFeatures.make<TH1F>("h_zOfflineVertex","z-coordinate of offline vertex;z_{vertex};n_{events}",100,-30.,30.);
953  h_xErrOfflineVertex = EventFeatures.make<TH1F>("h_xErrOfflineVertex","x-coordinate error of offline vertex;err_{x}^{vtx};n_{events}",100,0.,0.01);
954  h_yErrOfflineVertex = EventFeatures.make<TH1F>("h_yErrOfflineVertex","y-coordinate error of offline vertex;err_{y}^{vtx};n_{events}",100,0.,0.01);
955  h_zErrOfflineVertex = EventFeatures.make<TH1F>("h_zErrOfflineVertex","z-coordinate error of offline vertex;err_{z}^{vtx};n_{events}",100,0.,10.);
956  h_BSx0 = EventFeatures.make<TH1F>("h_BSx0","x-coordinate of reco beamspot;x^{BS}_{0};n_{events}",100,-0.1,0.1);
957  h_BSy0 = EventFeatures.make<TH1F>("h_BSy0","y-coordinate of reco beamspot;y^{BS}_{0};n_{events}",100,-0.1,0.1);
958  h_BSz0 = EventFeatures.make<TH1F>("h_BSz0","z-coordinate of reco beamspot;z^{BS}_{0};n_{events}",100,-1.,1.);
959  h_Beamsigmaz = EventFeatures.make<TH1F>("h_Beamsigmaz","z-coordinate beam width;#sigma_{Z}^{beam};n_{events}",100,0.,1.);
960  h_BeamWidthX = EventFeatures.make<TH1F>("h_BeamWidthX","x-coordinate beam width;#sigma_{X}^{beam};n_{events}",100,0.,0.01);
961  h_BeamWidthY = EventFeatures.make<TH1F>("h_BeamWidthY","y-coordinate beam width;#sigma_{Y}^{beam};n_{events}",100,0.,0.01);
962 
963  // probe track histograms
964  TFileDirectory ProbeFeatures = fs->mkdir("ProbeTrackFeatures");
965 
966  h_probePt_ = ProbeFeatures.make<TH1F>("h_probePt","p_{T} of probe track;track p_{T} (GeV); tracks",100,0.,50.);
967  h_probeP_ = ProbeFeatures.make<TH1F>("h_probeP","momentum of probe track;track p (GeV); tracks",100,0.,100.);
968  h_probeEta_ = ProbeFeatures.make<TH1F>("h_probeEta","#eta of the probe track;track #eta;tracks",54,-2.7,2.7);
969  h_probePhi_ = ProbeFeatures.make<TH1F>("h_probePhi","#phi of probe track;track #phi (rad);tracks",100,-3.15,3.15);
970 
971  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.7,2.7,100,-3.15,3.15);
972  h2_probeEtaPt_ = ProbeFeatures.make<TH2F>("h2_probeEtaPt","probe track p_{T} vs #eta;#eta of probe track;track p_{T} (GeV); tracks",54,-2.7,2.7,100,0.,50.);
973 
974  h_probeChi2_ = ProbeFeatures.make<TH1F>("h_probeChi2","#chi^{2} of probe track;track #chi^{2}; tracks",100,0.,100.);
975  h_probeNormChi2_ = ProbeFeatures.make<TH1F>("h_probeNormChi2"," normalized #chi^{2} of probe track;track #chi^{2}/ndof; tracks",100,0.,10.);
976  h_probeCharge_ = ProbeFeatures.make<TH1F>("h_probeCharge","charge of profe track;track charge Q;tracks",3,-1.5,1.5);
977  h_probeQoverP_ = ProbeFeatures.make<TH1F>("h_probeQoverP","q/p of probe track; track Q/p (GeV^{-1});tracks",200,-1.,1.);
978  h_probedzRecoV_ = ProbeFeatures.make<TH1F>("h_probedzRecoV","d_{z}(V_{offline}) of probe track;track d_{z}(V_{off}) (cm);tracks",200,-1.,1.);
979  h_probedxyRecoV_ = ProbeFeatures.make<TH1F>("h_probedxyRecoV","d_{xy}(V_{offline}) of probe track;track d_{xy}(V_{off}) (cm);tracks",200,-1.,1.);
980  h_probedzRefitV_ = ProbeFeatures.make<TH1F>("h_probedzRefitV","d_{z}(V_{refit}) of probe track;track d_{z}(V_{fit}) (cm);tracks",200,-1.,1.);
981  h_probesignIP2DRefitV_ = ProbeFeatures.make<TH1F>("h_probesignIPRefitV","ip_{2D}(V_{refit}) of probe track;track ip_{2D}(V_{fit}) (cm);tracks",200,-1.,1.);
982  h_probedxyRefitV_ = ProbeFeatures.make<TH1F>("h_probedxyRefitV","d_{xy}(V_{refit}) of probe track;track d_{xy}(V_{fit}) (cm);tracks",200,-1.,1.);
983 
984  h_probez0RefitV_ = ProbeFeatures.make<TH1F>("h_probez0RefitV","z_{0}(V_{refit}) of probe track;track z_{0}(V_{fit}) (cm);tracks",200,-1.,1.);
985  h_probed0RefitV_ = ProbeFeatures.make<TH1F>("h_probed0RefitV","d_{0}(V_{refit}) of probe track;track d_{0}(V_{fit}) (cm);tracks",200,-1.,1.);
986 
987  h_probed3DRefitV_ = ProbeFeatures.make<TH1F>("h_probed3DRefitV","d_{3D}(V_{refit}) of probe track;track d_{3D}(V_{fit}) (cm);tracks",200,0.,1.);
988  h_probereszRefitV_ = ProbeFeatures.make<TH1F>("h_probeReszRefitV","z_{track} -z_{V_{refit}};track res_{z}(V_{refit}) (cm);tracks",200,-1.,1.);
989 
990  h_probeRecoVSigZ_ = ProbeFeatures.make<TH1F>("h_probeRecoVSigZ" ,"Longitudinal DCA Significance (reco);d_{z}(V_{off})/#sigma_{dz};tracks",100,-8,8);
991  h_probeRecoVSigXY_ = ProbeFeatures.make<TH1F>("h_probeRecoVSigXY" ,"Transverse DCA Significance (reco);d_{xy}(V_{off})/#sigma_{dxy};tracks",100,-8,8);
992  h_probeRefitVSigZ_ = ProbeFeatures.make<TH1F>("h_probeRefitVSigZ" ,"Longitudinal DCA Significance (refit);d_{z}(V_{fit})/#sigma_{dz};tracks",100,-8,8);
993  h_probeRefitVSigXY_= ProbeFeatures.make<TH1F>("h_probeRefitVSigXY","Transverse DCA Significance (refit);d_{xy}(V_{fit})/#sigma_{dxy};tracks",100,-8,8);
994  h_probeRefitVSig3D_= ProbeFeatures.make<TH1F>("h_probeRefitVSig3D","3D DCA Significance (refit);d_{3D}/#sigma_{3D};tracks",100,0.,20.);
995  h_probeRefitVLogSig3D_ = ProbeFeatures.make<TH1F>("h_probeRefitVLogSig3D","log_{10}(3D DCA-Significance) (refit);log_{10}(d_{3D}/#sigma_{3D});tracks",100,-5.,4.);
996  h_probeRefitVSigResZ_ = ProbeFeatures.make<TH1F>("h_probeRefitVSigResZ" ,"Longitudinal residual significance (refit);(z_{track} -z_{V_{fit}})/#sigma_{res_{z}};tracks",100,-8,8);
997 
998  h_probeHits_ = ProbeFeatures.make<TH1F>("h_probeNRechits" ,"N_{hits} ;N_{hits} ;tracks",40,-0.5,39.5);
999  h_probeHits1D_ = ProbeFeatures.make<TH1F>("h_probeNRechits1D" ,"N_{hits} 1D ;N_{hits} 1D ;tracks",40,-0.5,39.5);
1000  h_probeHits2D_ = ProbeFeatures.make<TH1F>("h_probeNRechits2D" ,"N_{hits} 2D ;N_{hits} 2D ;tracks",40,-0.5,39.5);
1001  h_probeHitsInTIB_ = ProbeFeatures.make<TH1F>("h_probeNRechitsTIB" ,"N_{hits} TIB ;N_{hits} TIB;tracks",40,-0.5,39.5);
1002  h_probeHitsInTOB_ = ProbeFeatures.make<TH1F>("h_probeNRechitsTOB" ,"N_{hits} TOB ;N_{hits} TOB;tracks",40,-0.5,39.5);
1003  h_probeHitsInTID_ = ProbeFeatures.make<TH1F>("h_probeNRechitsTID" ,"N_{hits} TID ;N_{hits} TID;tracks",40,-0.5,39.5);
1004  h_probeHitsInTEC_ = ProbeFeatures.make<TH1F>("h_probeNRechitsTEC" ,"N_{hits} TEC ;N_{hits} TEC;tracks",40,-0.5,39.5);
1005  h_probeHitsInBPIX_ = ProbeFeatures.make<TH1F>("h_probeNRechitsBPIX","N_{hits} BPIX;N_{hits} BPIX;tracks",40,-0.5,39.5);
1006  h_probeHitsInFPIX_ = ProbeFeatures.make<TH1F>("h_probeNRechitsFPIX","N_{hits} FPIX;N_{hits} FPIX;tracks",40,-0.5,39.5);
1007 
1008  // refit vertex features
1009  TFileDirectory RefitVertexFeatures = fs->mkdir("RefitVertexFeatures");
1010  h_fitVtxNtracks_ = RefitVertexFeatures.make<TH1F>("h_fitVtxNtracks" ,"N_{trks} used in vertex fit;N^{fit}_{tracks};vertices" ,100,-0.5,99.5);
1011  h_fitVtxNdof_ = RefitVertexFeatures.make<TH1F>("h_fitVtxNdof" ,"N_{DOF} of vertex fit;N_{DOF} of refit vertex;vertices" ,100,-0.5,99.5);
1012  h_fitVtxChi2_ = RefitVertexFeatures.make<TH1F>("h_fitVtxChi2" ,"#chi^{2} of vertex fit;vertex #chi^{2};vertices" ,100,-0.5,99.5);
1013  h_fitVtxChi2ndf_ = RefitVertexFeatures.make<TH1F>("h_fitVtxChi2ndf" ,"#chi^{2}/ndf of vertex fit;vertex #chi^{2}/ndf;vertices" ,100,-0.5,9.5);
1014  h_fitVtxChi2Prob_ = RefitVertexFeatures.make<TH1F>("h_fitVtxChi2Prob" ,"Prob(#chi^{2},ndf) of vertex fit;Prob(#chi^{2},ndf);vertices",40,0.,1.);
1015  h_fitVtxTrackWeights_ = RefitVertexFeatures.make<TH1F>("h_fitVtxTrackWeights","track weights associated to track;track weights;tracks",40,0.,1.);
1016  h_fitVtxTrackAverageWeight_ = RefitVertexFeatures.make<TH1F>("h_fitVtxTrackAverageWeight_","average track weight per vertex;#LT track weight #GT;vertices",40,0.,1.);
1017 
1018  if(useTracksFromRecoVtx_) {
1019 
1020  TFileDirectory RecoVertexFeatures = fs->mkdir("RecoVertexFeatures");
1021  h_recoVtxNtracks_ = RecoVertexFeatures.make<TH1F>("h_recoVtxNtracks" ,"N^{vtx}_{trks};N^{vtx}_{trks};vertices" ,100,-0.5,99.5);
1022  h_recoVtxChi2ndf_ = RecoVertexFeatures.make<TH1F>("h_recoVtxChi2ndf" ,"#chi^{2}/ndf vtx;#chi^{2}/ndf vtx;vertices" ,10,-0.5,9.5);
1023  h_recoVtxChi2Prob_ = RecoVertexFeatures.make<TH1F>("h_recoVtxChi2Prob" ,"Prob(#chi^{2},ndf);Prob(#chi^{2},ndf);vertices",40,0.,1.);
1024  h_recoVtxSumPt_ = RecoVertexFeatures.make<TH1F>("h_recoVtxSumPt" ,"Sum(p^{trks}_{T});Sum(p^{trks}_{T});vertices" ,100,0.,200.);
1025 
1026  }
1027 
1028 
1029  TFileDirectory DA = fs->mkdir("DA");
1030  //DA.cd();
1031  hDA=bookVertexHistograms(DA);
1032  //for(std::map<std::string,TH1*>::const_iterator hist=hDA.begin(); hist!=hDA.end(); hist++){
1033  //hist->second->SetDirectory(DA);
1034  // DA.make<TH1F>(hist->second);
1035  // }
1036 
1037  // initialize the residuals histograms
1038 
1039  float dxymax_phi = 2000;
1040  float dzmax_phi = 2000;
1041  float dxymax_eta = 3000;
1042  float dzmax_eta = 3000;
1043 
1044  float d3Dmax_phi = hypot(dxymax_phi,dzmax_phi);
1045  float d3Dmax_eta = hypot(dxymax_eta,dzmax_eta);
1046 
1047  const int mybins_ = 500;
1048 
1050  //
1051  // Unbiased track-to-vertex residuals
1052  // The vertex is refit without the probe track
1053  //
1055 
1056  TFileDirectory AbsTransPhiRes = fs->mkdir("Abs_Transv_Phi_Residuals");
1057  TFileDirectory AbsTransEtaRes = fs->mkdir("Abs_Transv_Eta_Residuals");
1058 
1059  TFileDirectory AbsLongPhiRes = fs->mkdir("Abs_Long_Phi_Residuals");
1060  TFileDirectory AbsLongEtaRes = fs->mkdir("Abs_Long_Eta_Residuals");
1061 
1062  TFileDirectory Abs3DPhiRes = fs->mkdir("Abs_3D_Phi_Residuals");
1063  TFileDirectory Abs3DEtaRes = fs->mkdir("Abs_3D_Eta_Residuals");
1064 
1065  TFileDirectory NormTransPhiRes = fs->mkdir("Norm_Transv_Phi_Residuals");
1066  TFileDirectory NormTransEtaRes = fs->mkdir("Norm_Transv_Eta_Residuals");
1067 
1068  TFileDirectory NormLongPhiRes = fs->mkdir("Norm_Long_Phi_Residuals");
1069  TFileDirectory NormLongEtaRes = fs->mkdir("Norm_Long_Eta_Residuals");
1070 
1071  TFileDirectory Norm3DPhiRes = fs->mkdir("Norm_3D_Phi_Residuals");
1072  TFileDirectory Norm3DEtaRes = fs->mkdir("Norm_3D_Eta_Residuals");
1073 
1074  TFileDirectory AbsDoubleDiffRes = fs->mkdir("Abs_DoubleDiffResiduals");
1075  TFileDirectory NormDoubleDiffRes = fs->mkdir("Norm_DoubleDiffResiduals");
1076 
1077  for ( int i=0; i<nBins_; ++i ) {
1078 
1079  float phiF = (-TMath::Pi()+i*phiSect_)*(180/TMath::Pi());
1080  float phiL = (-TMath::Pi()+(i+1)*phiSect_)*(180/TMath::Pi());
1081 
1082  float etaF=-2.5+i*etaSect_;
1083  float etaL=-2.5+(i+1)*etaSect_;
1084 
1085  // dxy vs phi and eta
1086 
1087  a_dxyPhiResiduals[i] = AbsTransPhiRes.make<TH1F>(Form("histo_dxy_phi_plot%i",i),
1088  Form("%.2f#circ<#varphi^{probe}_{tk}<%.2f#circ;d_{xy} [#mum];tracks",phiF,phiL),
1089  mybins_,-dxymax_phi,dxymax_phi);
1090 
1091  a_dxyEtaResiduals[i] = AbsTransEtaRes.make<TH1F>(Form("histo_dxy_eta_plot%i",i),
1092  Form("%.2f<#eta^{probe}_{tk}<%.2f;d_{xy} [#mum];tracks",etaF,etaL),
1093  mybins_,-dxymax_eta,dxymax_eta);
1094 
1095  // IP2D vs phi and eta
1096 
1097  a_IP2DPhiResiduals[i] = AbsTransPhiRes.make<TH1F>(Form("histo_IP2D_phi_plot%i",i),
1098  Form("%.2f#circ<#varphi^{probe}_{tk}<%.2f#circ;IP_{2D} [#mum];tracks",phiF,phiL),
1099  mybins_,-dxymax_phi,dxymax_phi);
1100 
1101  a_IP2DEtaResiduals[i] = AbsTransEtaRes.make<TH1F>(Form("histo_IP2D_eta_plot%i",i),
1102  Form("%.2f<#eta^{probe}_{tk}<%.2f;IP_{2D} [#mum];tracks",etaF,etaL),
1103  mybins_,-dxymax_eta,dxymax_eta);
1104 
1105  // IP3D vs phi and eta
1106 
1107  a_IP3DPhiResiduals[i] = Abs3DPhiRes.make<TH1F>(Form("histo_IP3D_phi_plot%i",i),
1108  Form("%.2f#circ<#varphi^{probe}_{tk}<%.2f#circ;IP_{3D} [#mum];tracks",phiF,phiL),
1109  mybins_,-dxymax_phi,dxymax_phi);
1110 
1111  a_IP3DEtaResiduals[i] = Abs3DEtaRes.make<TH1F>(Form("histo_IP3D_eta_plot%i",i),
1112  Form("%.2f<#eta^{probe}_{tk}<%.2f;IP_{3D} [#mum];tracks",etaF,etaL),
1113  mybins_,-dxymax_eta,dxymax_eta);
1114 
1115  // dz vs phi and eta
1116 
1117  a_dzPhiResiduals[i] = AbsLongPhiRes.make<TH1F>(Form("histo_dz_phi_plot%i",i),
1118  Form("%.2f#circ<#varphi^{probe}_{tk}<%.2f#circ;d_{z} [#mum];tracks",phiF,phiL),
1119  mybins_,-dzmax_phi,dzmax_phi);
1120 
1121  a_dzEtaResiduals[i] = AbsLongEtaRes.make<TH1F>(Form("histo_dz_eta_plot%i",i),
1122  Form("%.2f<#eta^{probe}_{tk}<%.2f;d_{z} [#mum];tracks",etaF,etaL),
1123  mybins_,-dzmax_eta,dzmax_eta);
1124 
1125 
1126  // resz vs phi and eta
1127 
1128  a_reszPhiResiduals[i] = AbsLongPhiRes.make<TH1F>(Form("histo_resz_phi_plot%i",i),
1129  Form("%.2f#circ<#varphi^{probe}_{tk}<%.2f#circ;z_{trk} - z_{vtx} [#mum];tracks",phiF,phiL),
1130  mybins_,-dzmax_phi,dzmax_phi);
1131 
1132  a_reszEtaResiduals[i] = AbsLongEtaRes.make<TH1F>(Form("histo_resz_eta_plot%i",i),
1133  Form("%.2f<#eta^{probe}_{tk}<%.2f;z_{trk} - z_{vtx} [#mum];tracks",etaF,etaL),
1134  mybins_,-dzmax_eta,dzmax_eta);
1135 
1136  // d3D vs phi and eta
1137 
1138  a_d3DPhiResiduals[i] = Abs3DPhiRes.make<TH1F>(Form("histo_d3D_phi_plot%i",i),
1139  Form("%.2f#circ<#varphi^{probe}_{tk}<%.2f#circ;d_{3D} [#mum];tracks",phiF,phiL),
1140  mybins_,0.,d3Dmax_phi);
1141 
1142  a_d3DEtaResiduals[i] = Abs3DEtaRes.make<TH1F>(Form("histo_d3D_eta_plot%i",i),
1143  Form("%.2f<#eta^{probe}_{tk}<%.2f;d_{3D} [#mum];tracks",etaF,etaL),
1144  mybins_,0.,d3Dmax_eta);
1145 
1146  // normalized dxy vs eta and phi
1147 
1148  n_dxyPhiResiduals[i] = NormTransPhiRes.make<TH1F>(Form("histo_norm_dxy_phi_plot%i",i),
1149  Form("%.2f#circ<#varphi^{probe}_{tk}<%.2f#circ;d_{xy}/#sigma_{d_{xy}};tracks",phiF,phiL),
1150  mybins_,-dxymax_phi/100.,dxymax_phi/100.);
1151 
1152  n_dxyEtaResiduals[i] = NormTransEtaRes.make<TH1F>(Form("histo_norm_dxy_eta_plot%i",i),
1153  Form("%.2f<#eta^{probe}_{tk}<%.2f;d_{xy}/#sigma_{d_{xy}};tracks",etaF,etaL),
1154  mybins_,-dxymax_eta/100.,dxymax_eta/100.);
1155 
1156  // normalized IP2d vs eta and phi
1157 
1158  n_IP2DPhiResiduals[i] = NormTransPhiRes.make<TH1F>(Form("histo_norm_IP2D_phi_plot%i",i),
1159  Form("%.2f#circ<#varphi^{probe}_{tk}<%.2f#circ;IP_{2D}/#sigma_{IP_{2D}};tracks",phiF,phiL),
1160  mybins_,-dxymax_phi/100.,dxymax_phi/100.);
1161 
1162  n_IP2DEtaResiduals[i] = NormTransEtaRes.make<TH1F>(Form("histo_norm_IP2D_eta_plot%i",i),
1163  Form("%.2f<#eta^{probe}_{tk}<%.2f;IP_{2D}/#sigma_{IP_{2D}};tracks",etaF,etaL),
1164  mybins_,-dxymax_eta/100.,dxymax_eta/100.);
1165 
1166  // normalized IP3d vs eta and phi
1167 
1168  n_IP3DPhiResiduals[i] = Norm3DPhiRes.make<TH1F>(Form("histo_norm_IP3D_phi_plot%i",i),
1169  Form("%.2f#circ<#varphi^{probe}_{tk}<%.2f#circ;IP_{3D}/#sigma_{IP_{3D}};tracks",phiF,phiL),
1170  mybins_,-dxymax_phi/100.,dxymax_phi/100.);
1171 
1172  n_IP3DEtaResiduals[i] = Norm3DEtaRes.make<TH1F>(Form("histo_norm_IP3D_eta_plot%i",i),
1173  Form("%.2f<#eta^{probe}_{tk}<%.2f;IP_{3D}/#sigma_{IP_{3D}};tracks",etaF,etaL),
1174  mybins_,-dxymax_eta/100.,dxymax_eta/100.);
1175 
1176  // normalized dz vs phi and eta
1177 
1178  n_dzPhiResiduals[i] = NormLongPhiRes.make<TH1F>(Form("histo_norm_dz_phi_plot%i",i),
1179  Form("%.2f#circ<#varphi^{probe}_{tk}<%.2f#circ;d_{z}/#sigma_{d_{z}};tracks",phiF,phiL),
1180  mybins_,-dzmax_phi/100.,dzmax_phi/100.);
1181 
1182  n_dzEtaResiduals[i] = NormLongEtaRes.make<TH1F>(Form("histo_norm_dz_eta_plot%i",i),
1183  Form("%.2f<#eta^{probe}_{tk}<%.2f;d_{z}/#sigma_{d_{z}};tracks",etaF,etaL),
1184  mybins_,-dzmax_eta/100.,dzmax_eta/100.);
1185 
1186  // pull of resz
1187 
1188  n_reszPhiResiduals[i] = NormLongPhiRes.make<TH1F>(Form("histo_norm_resz_phi_plot%i",i),
1189  Form("%.2f#circ<#varphi^{probe}_{tk}<%.2f#circ;(z_{trk}-z_{vtx})/#sigma_{res_{z}};tracks",phiF,phiL),
1190  mybins_,-dzmax_phi/100.,dzmax_phi/100.);
1191 
1192  n_reszEtaResiduals[i] = NormLongEtaRes.make<TH1F>(Form("histo_norm_resz_eta_plot%i",i),
1193  Form("%.2f<#eta^{probe}_{tk}<%.2f;(z_{trk}-z_{vtx})/#sigma_{res_{z}};tracks",etaF,etaL),
1194  mybins_,-dzmax_eta/100.,dzmax_eta/100.);
1195 
1196  // normalized d3D vs phi and eta
1197 
1198  n_d3DPhiResiduals[i] = Norm3DPhiRes.make<TH1F>(Form("histo_norm_d3D_phi_plot%i",i),
1199  Form("%.2f#circ<#varphi^{probe}_{tk}<%.2f#circ;d_{3D}/#sigma_{d_{3D}};tracks",phiF,phiL),
1200  mybins_,0.,d3Dmax_phi/100.);
1201 
1202  n_d3DEtaResiduals[i] = Norm3DEtaRes.make<TH1F>(Form("histo_norm_d3D_eta_plot%i",i),
1203  Form("%.2f<#eta^{probe}_{tk}<%.2f;d_{3D}/#sigma_{d_{3D}};tracks",etaF,etaL),
1204  mybins_,0.,d3Dmax_eta/100.);
1205 
1206  for ( int j=0; j<nBins_; ++j ) {
1207 
1208  a_dxyResidualsMap[i][j] = AbsDoubleDiffRes.make<TH1F>(Form("histo_dxy_eta_plot%i_phi_plot%i",i,j),
1209  Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{xy};tracks",etaF,etaL,phiF,phiL),
1210  mybins_,-dzmax_eta,dzmax_eta);
1211 
1212  a_dzResidualsMap[i][j] = AbsDoubleDiffRes.make<TH1F>(Form("histo_dz_eta_plot%i_phi_plot%i",i,j),
1213  Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{z};tracks",etaF,etaL,phiF,phiL),
1214  mybins_,-dzmax_eta,dzmax_eta);
1215 
1216  a_d3DResidualsMap[i][j] = AbsDoubleDiffRes.make<TH1F>(Form("histo_d3D_eta_plot%i_phi_plot%i",i,j),
1217  Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{3D};tracks",etaF,etaL,phiF,phiL),
1218  mybins_,0.,d3Dmax_eta);
1219 
1220  n_dxyResidualsMap[i][j] = NormDoubleDiffRes.make<TH1F>(Form("histo_norm_dxy_eta_plot%i_phi_plot%i",i,j),
1221  Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{xy}/#sigma_{d_{xy}};tracks",etaF,etaL,phiF,phiL),
1222  mybins_,-dzmax_eta/100,dzmax_eta/100);
1223 
1224  n_dzResidualsMap[i][j] = NormDoubleDiffRes.make<TH1F>(Form("histo_norm_dz_eta_plot%i_phi_plot%i",i,j),
1225  Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{z}/#sigma_{d_{z}};tracks",etaF,etaL,phiF,phiL),
1226  mybins_,-dzmax_eta/100,dzmax_eta/100);
1227 
1228  n_d3DResidualsMap[i][j] = NormDoubleDiffRes.make<TH1F>(Form("histo_norm_d3D_eta_plot%i_phi_plot%i",i,j),
1229  Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{3D}/#sigma_{d_{3D}};tracks",etaF,etaL,phiF,phiL),
1230  mybins_,0.,d3Dmax_eta);
1231 
1232  }
1233  }
1234 
1235  // declaration of the directories
1236 
1237  TFileDirectory BiasVsParameter = fs->mkdir("BiasVsParameter");
1238 
1239  a_dxyVsPhi = BiasVsParameter.make<TH2F>("h2_dxy_vs_phi","d_{xy} vs track #phi;track #phi [rad];track d_{xy}(PV) [#mum]",
1240  48,-TMath::Pi(),TMath::Pi(),mybins_,-dxymax_phi,dxymax_phi);
1241 
1242  a_dzVsPhi = BiasVsParameter.make<TH2F>("h2_dz_vs_phi","d_{z} vs track #phi;track #phi [rad];track d_{z}(PV) [#mum]",
1243  48,-TMath::Pi(),TMath::Pi(),mybins_,-dzmax_phi,dzmax_phi);
1244 
1245  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}}",
1246  48,-TMath::Pi(),TMath::Pi(),mybins_,-dxymax_phi/100.,dxymax_phi/100.);
1247 
1248  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}}",
1249  48,-TMath::Pi(),TMath::Pi(),mybins_,-dzmax_phi/100.,dzmax_phi/100.);
1250 
1251  a_dxyVsEta = BiasVsParameter.make<TH2F>("h2_dxy_vs_eta","d_{xy} vs track #eta;track #eta;track d_{xy}(PV) [#mum]",
1252  48,-2.5,2.5,mybins_,-dxymax_eta,dzmax_eta);
1253 
1254  a_dzVsEta = BiasVsParameter.make<TH2F>("h2_dz_vs_eta","d_{z} vs track #eta;track #eta;track d_{z}(PV) [#mum]",
1255  48,-2.5,2.5,mybins_,-dzmax_eta,dzmax_eta);
1256 
1257  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}}",
1258  48,-2.5,2.5,mybins_,-dxymax_eta/100.,dxymax_eta/100.);
1259 
1260  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}}",
1261  48,-2.5,2.5,mybins_,-dzmax_eta/100.,dzmax_eta/100.);
1262 
1263  TFileDirectory MeanTrendsDir = fs->mkdir("MeanTrends");
1264  TFileDirectory WidthTrendsDir = fs->mkdir("WidthTrends");
1265  TFileDirectory MedianTrendsDir = fs->mkdir("MedianTrends");
1266  TFileDirectory MADTrendsDir = fs->mkdir("MADTrends");
1267 
1268  TFileDirectory Mean2DMapsDir = fs->mkdir("MeanMaps");
1269  TFileDirectory Width2DMapsDir = fs->mkdir("WidthMaps");
1270 
1271  Double_t highedge=nBins_-0.5;
1272  Double_t lowedge=-0.5;
1273 
1274  // means and widths from the fit
1275 
1276  a_dxyPhiMeanTrend = MeanTrendsDir.make<TH1F> ("means_dxy_phi",
1277  "#LT d_{xy} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{xy} #GT [#mum]",
1278  nBins_,lowedge,highedge);
1279 
1280  a_dxyPhiWidthTrend = WidthTrendsDir.make<TH1F>("widths_dxy_phi",
1281  "#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees];#sigma_{d_{xy}} [#mum]",
1282  nBins_,lowedge,highedge);
1283 
1284  a_dzPhiMeanTrend = MeanTrendsDir.make<TH1F> ("means_dz_phi",
1285  "#LT d_{z} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{z} #GT [#mum]",
1286  nBins_,lowedge,highedge);
1287 
1288  a_dzPhiWidthTrend = WidthTrendsDir.make<TH1F>("widths_dz_phi","#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];#sigma_{d_{z}} [#mum]",
1289  nBins_,lowedge,highedge);
1290 
1291  a_dxyEtaMeanTrend = MeanTrendsDir.make<TH1F> ("means_dxy_eta",
1292  "#LT d_{xy} #GT vs #eta sector;#eta (sector);#LT d_{xy} #GT [#mum]",
1293  nBins_,lowedge,highedge);
1294 
1295  a_dxyEtaWidthTrend = WidthTrendsDir.make<TH1F>("widths_dxy_eta",
1296  "#sigma_{d_{xy}} vs #eta sector;#eta (sector);#sigma_{d_{xy}} [#mum]",
1297  nBins_,lowedge,highedge);
1298 
1299  a_dzEtaMeanTrend = MeanTrendsDir.make<TH1F> ("means_dz_eta",
1300  "#LT d_{z} #GT vs #eta sector;#eta (sector);#LT d_{z} #GT [#mum]"
1301  ,nBins_,lowedge,highedge);
1302 
1303  a_dzEtaWidthTrend = WidthTrendsDir.make<TH1F>("widths_dz_eta",
1304  "#sigma_{d_{xy}} vs #eta sector;#eta (sector);#sigma_{d_{z}} [#mum]",
1305  nBins_,lowedge,highedge);
1306 
1307  n_dxyPhiMeanTrend = MeanTrendsDir.make<TH1F> ("norm_means_dxy_phi",
1308  "#LT d_{xy}/#sigma_{d_{xy}} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{xy}/#sigma_{d_{xy}} #GT",
1309  nBins_,lowedge,highedge);
1310 
1311  n_dxyPhiWidthTrend = WidthTrendsDir.make<TH1F>("norm_widths_dxy_phi",
1312  "width(d_{xy}/#sigma_{d_{xy}}) vs #phi sector;#varphi (sector) [degrees]; width(d_{xy}/#sigma_{d_{xy}})",
1313  nBins_,lowedge,highedge);
1314 
1315  n_dzPhiMeanTrend = MeanTrendsDir.make<TH1F> ("norm_means_dz_phi",
1316  "#LT d_{z}/#sigma_{d_{z}} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{z}/#sigma_{d_{z}} #GT",
1317  nBins_,lowedge,highedge);
1318 
1319  n_dzPhiWidthTrend = WidthTrendsDir.make<TH1F>("norm_widths_dz_phi",
1320  "width(d_{z}/#sigma_{d_{z}}) vs #phi sector;#varphi (sector) [degrees];width(d_{z}/#sigma_{d_{z}})",
1321  nBins_,lowedge,highedge);
1322 
1323  n_dxyEtaMeanTrend = MeanTrendsDir.make<TH1F> ("norm_means_dxy_eta",
1324  "#LT d_{xy}/#sigma_{d_{xy}} #GT vs #eta sector;#eta (sector);#LT d_{xy}/#sigma_{d_{z}} #GT",
1325  nBins_,lowedge,highedge);
1326 
1327  n_dxyEtaWidthTrend = WidthTrendsDir.make<TH1F>("norm_widths_dxy_eta",
1328  "width(d_{xy}/#sigma_{d_{xy}}) vs #eta sector;#eta (sector);width(d_{xy}/#sigma_{d_{z}})",
1329  nBins_,lowedge,highedge);
1330 
1331  n_dzEtaMeanTrend = MeanTrendsDir.make<TH1F> ("norm_means_dz_eta",
1332  "#LT d_{z}/#sigma_{d_{z}} #GT vs #eta sector;#eta (sector);#LT d_{z}/#sigma_{d_{z}} #GT",
1333  nBins_,lowedge,highedge);
1334 
1335  n_dzEtaWidthTrend = WidthTrendsDir.make<TH1F>("norm_widths_dz_eta",
1336  "width(d_{z}/#sigma_{d_{z}}) vs #eta sector;#eta (sector);width(d_{z}/#sigma_{d_{z}})",
1337  nBins_,lowedge,highedge);
1338 
1339  // 2D maps
1340 
1341  a_dxyMeanMap = Mean2DMapsDir.make<TH2F> ("means_dxy_map",
1342  "#LT d_{xy} #GT map;#eta (sector);#varphi (sector) [degrees]",
1343  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1344 
1345  a_dzMeanMap = Mean2DMapsDir.make<TH2F> ("means_dz_map",
1346  "#LT d_{z} #GT map;#eta (sector);#varphi (sector) [degrees]",
1347  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1348 
1349  n_dxyMeanMap = Mean2DMapsDir.make<TH2F> ("norm_means_dxy_map",
1350  "#LT d_{xy}/#sigma_{d_{xy}} #GT map;#eta (sector);#varphi (sector) [degrees]",
1351  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1352 
1353  n_dzMeanMap = Mean2DMapsDir.make<TH2F> ("norm_means_dz_map",
1354  "#LT d_{z}/#sigma_{d_{z}} #GT map;#eta (sector);#varphi (sector) [degrees]",
1355  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1356 
1357  a_dxyWidthMap = Width2DMapsDir.make<TH2F> ("widths_dxy_map",
1358  "#sigma_{d_{xy}} map;#eta (sector);#varphi (sector) [degrees]",
1359  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1360 
1361  a_dzWidthMap = Width2DMapsDir.make<TH2F> ("widths_dz_map",
1362  "#sigma_{d_{z}} map;#eta (sector);#varphi (sector) [degrees]",
1363  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1364 
1365  n_dxyWidthMap = Width2DMapsDir.make<TH2F> ("norm_widths_dxy_map",
1366  "width(d_{xy}/#sigma_{d_{xy}}) map;#eta (sector);#varphi (sector) [degrees]",
1367  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1368 
1369  n_dzWidthMap = Width2DMapsDir.make<TH2F> ("norm_widths_dz_map",
1370  "width(d_{z}/#sigma_{d_{z}}) map;#eta (sector);#varphi (sector) [degrees]",
1371  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1372 
1373  // medians and MADs
1374 
1375  a_dxyPhiMedianTrend = MedianTrendsDir.make<TH1F>("medians_dxy_phi",
1376  "Median of d_{xy} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{xy}) [#mum]",
1377  nBins_,lowedge,highedge);
1378 
1379  a_dxyPhiMADTrend = MADTrendsDir.make<TH1F> ("MADs_dxy_phi",
1380  "Median absolute deviation of d_{xy} vs #phi sector;#varphi (sector) [degrees];MAD(d_{xy}) [#mum]",
1381  nBins_,lowedge,highedge);
1382 
1383  a_dzPhiMedianTrend = MedianTrendsDir.make<TH1F>("medians_dz_phi",
1384  "Median of d_{z} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{z}) [#mum]",
1385  nBins_,lowedge,highedge);
1386 
1387  a_dzPhiMADTrend = MADTrendsDir.make<TH1F> ("MADs_dz_phi",
1388  "Median absolute deviation of d_{z} vs #phi sector;#varphi (sector) [degrees];MAD(d_{z}) [#mum]",
1389  nBins_,lowedge,highedge);
1390 
1391  a_dxyEtaMedianTrend = MedianTrendsDir.make<TH1F>("medians_dxy_eta",
1392  "Median of d_{xy} vs #eta sector;#eta (sector);#mu_{1/2}(d_{xy}) [#mum]",
1393  nBins_,lowedge,highedge);
1394 
1395  a_dxyEtaMADTrend = MADTrendsDir.make<TH1F> ("MADs_dxy_eta",
1396  "Median absolute deviation of d_{xy} vs #eta sector;#eta (sector);MAD(d_{xy}) [#mum]",
1397  nBins_,lowedge,highedge);
1398 
1399  a_dzEtaMedianTrend = MedianTrendsDir.make<TH1F>("medians_dz_eta",
1400  "Median of d_{z} vs #eta sector;#eta (sector);#mu_{1/2}(d_{z}) [#mum]",
1401  nBins_,lowedge,highedge);
1402 
1403  a_dzEtaMADTrend = MADTrendsDir.make<TH1F> ("MADs_dz_eta",
1404  "Median absolute deviation of d_{z} vs #eta sector;#eta (sector);MAD(d_{z}) [#mum]",
1405  nBins_,lowedge,highedge);
1406 
1407  n_dxyPhiMedianTrend = MedianTrendsDir.make<TH1F>("norm_medians_dxy_phi",
1408  "Median of d_{xy}/#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{xy}/#sigma_{d_{xy}})",
1409  nBins_,lowedge,highedge);
1410 
1411  n_dxyPhiMADTrend = MADTrendsDir.make<TH1F> ("norm_MADs_dxy_phi",
1412  "Median absolute deviation of d_{xy}/#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees]; MAD(d_{xy}/#sigma_{d_{xy}})",
1413  nBins_,lowedge,highedge);
1414 
1415  n_dzPhiMedianTrend = MedianTrendsDir.make<TH1F>("norm_medians_dz_phi",
1416  "Median of d_{z}/#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{z}/#sigma_{d_{z}})",
1417  nBins_,lowedge,highedge);
1418 
1419  n_dzPhiMADTrend = MADTrendsDir.make<TH1F> ("norm_MADs_dz_phi",
1420  "Median absolute deviation of d_{z}/#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];MAD(d_{z}/#sigma_{d_{z}})",
1421  nBins_,lowedge,highedge);
1422 
1423  n_dxyEtaMedianTrend = MedianTrendsDir.make<TH1F>("norm_medians_dxy_eta",
1424  "Median of d_{xy}/#sigma_{d_{xy}} vs #eta sector;#eta (sector);#mu_{1/2}(d_{xy}/#sigma_{d_{z}})",
1425  nBins_,lowedge,highedge);
1426 
1427  n_dxyEtaMADTrend = MADTrendsDir.make<TH1F> ("norm_MADs_dxy_eta",
1428  "Median absolute deviation of d_{xy}/#sigma_{d_{xy}} vs #eta sector;#eta (sector);MAD(d_{xy}/#sigma_{d_{z}})",
1429  nBins_,lowedge,highedge);
1430 
1431  n_dzEtaMedianTrend = MedianTrendsDir.make<TH1F>("norm_medians_dz_eta",
1432  "Median of d_{z}/#sigma_{d_{z}} vs #eta sector;#eta (sector);#mu_{1/2}(d_{z}/#sigma_{d_{z}})",
1433  nBins_,lowedge,highedge);
1434 
1435  n_dzEtaMADTrend = MADTrendsDir.make<TH1F> ("norm_MADs_dz_eta",
1436  "Median absolute deviation of d_{z}/#sigma_{d_{z}} vs #eta sector;#eta (sector);MAD(d_{z}/#sigma_{d_{z}})",
1437  nBins_,lowedge,highedge);
1438 
1439 
1441  //
1442  // plots of biased residuals
1443  // The vertex includes the probe track
1444  //
1446 
1447  if (useTracksFromRecoVtx_){
1448 
1449  TFileDirectory AbsTransPhiBiasRes = fs->mkdir("Abs_Transv_Phi_BiasResiduals");
1450  TFileDirectory AbsTransEtaBiasRes = fs->mkdir("Abs_Transv_Eta_BiasResiduals");
1451 
1452  TFileDirectory AbsLongPhiBiasRes = fs->mkdir("Abs_Long_Phi_BiasResiduals");
1453  TFileDirectory AbsLongEtaBiasRes = fs->mkdir("Abs_Long_Eta_BiasResiduals");
1454 
1455  TFileDirectory NormTransPhiBiasRes = fs->mkdir("Norm_Transv_Phi_BiasResiduals");
1456  TFileDirectory NormTransEtaBiasRes = fs->mkdir("Norm_Transv_Eta_BiasResiduals");
1457 
1458  TFileDirectory NormLongPhiBiasRes = fs->mkdir("Norm_Long_Phi_BiasResiduals");
1459  TFileDirectory NormLongEtaBiasRes = fs->mkdir("Norm_Long_Eta_BiasResiduals");
1460 
1461  TFileDirectory AbsDoubleDiffBiasRes = fs->mkdir("Abs_DoubleDiffBiasResiduals");
1462  TFileDirectory NormDoubleDiffBiasRes = fs->mkdir("Norm_DoubleDiffBiasResiduals");
1463 
1464  for ( int i=0; i<nBins_; ++i ) {
1465 
1466  float phiF = (-TMath::Pi()+i*phiSect_)*(180/TMath::Pi());
1467  float phiL = (-TMath::Pi()+(i+1)*phiSect_)*(180/TMath::Pi());
1468 
1469  float etaF=-2.5+i*etaSect_;
1470  float etaL=-2.5+(i+1)*etaSect_;
1471 
1472  a_dxyPhiBiasResiduals[i] = AbsTransPhiBiasRes.make<TH1F>(Form("histo_dxy_phi_plot%i",i),
1473  Form("%.2f#circ<#varphi^{probe}_{tk}<%.2f#circ;d_{xy} [#mum];tracks",phiF,phiL),
1474  mybins_,-dxymax_phi,dxymax_phi);
1475 
1476  a_dxyEtaBiasResiduals[i] = AbsTransEtaBiasRes.make<TH1F>(Form("histo_dxy_eta_plot%i",i),
1477  Form("%.2f<#eta^{probe}_{tk}<%.2f;d_{xy} [#mum];tracks",etaF,etaL),
1478  mybins_,-dxymax_eta,dxymax_eta);
1479 
1480  a_dzPhiBiasResiduals[i] = AbsLongPhiBiasRes.make<TH1F>(Form("histo_dz_phi_plot%i",i),
1481  Form("%.2f#circ<#varphi^{probe}_{tk}<%.2f #circ;d_{z} [#mum];tracks",phiF,phiL),
1482  mybins_,-dzmax_phi,dzmax_phi);
1483 
1484  a_dzEtaBiasResiduals[i] = AbsLongEtaBiasRes.make<TH1F>(Form("histo_dz_eta_plot%i",i),
1485  Form("%.2f<#eta^{probe}_{tk}<%.2f;d_{z} [#mum];tracks",etaF,etaL),
1486  mybins_,-dzmax_eta,dzmax_eta);
1487 
1488  n_dxyPhiBiasResiduals[i] = NormTransPhiBiasRes.make<TH1F>(Form("histo_norm_dxy_phi_plot%i",i),
1489  Form("%.2f#circ<#varphi^{probe}_{tk}<%.2f#circ;d_{xy}/#sigma_{d_{xy}};tracks",phiF,phiL),
1490  mybins_,-dxymax_phi/100.,dxymax_phi/100.);
1491 
1492  n_dxyEtaBiasResiduals[i] = NormTransEtaBiasRes.make<TH1F>(Form("histo_norm_dxy_eta_plot%i",i),
1493  Form("%.2f<#eta^{probe}_{tk}<%.2f;d_{xy}/#sigma_{d_{xy}};tracks",etaF,etaL),
1494  mybins_,-dxymax_eta/100.,dxymax_eta/100.);
1495 
1496  n_dzPhiBiasResiduals[i] = NormLongPhiBiasRes.make<TH1F>(Form("histo_norm_dz_phi_plot%i",i),
1497  Form("%.2f#circ<#varphi^{probe}_{tk}<%.2f#circ;d_{z}/#sigma_{d_{z}};tracks",phiF,phiL),
1498  mybins_,-dzmax_phi/100.,dzmax_phi/100.);
1499 
1500  n_dzEtaBiasResiduals[i] = NormLongEtaBiasRes.make<TH1F>(Form("histo_norm_dz_eta_plot%i",i),
1501  Form("%.2f<#eta^{probe}_{tk}<%.2f;d_{z}/#sigma_{d_{z}};tracks",etaF,etaL),
1502  mybins_,-dzmax_eta/100.,dzmax_eta/100.);
1503 
1504  for ( int j=0; j<nBins_; ++j ) {
1505 
1506  a_dxyBiasResidualsMap[i][j] = AbsDoubleDiffBiasRes.make<TH1F>(Form("histo_dxy_eta_plot%i_phi_plot%i",i,j),
1507  Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{xy} [#mum];tracks",etaF,etaL,phiF,phiL),
1508  mybins_,-dzmax_eta,dzmax_eta);
1509 
1510  a_dzBiasResidualsMap[i][j] = AbsDoubleDiffBiasRes.make<TH1F>(Form("histo_dxy_eta_plot%i_phi_plot%i",i,j),
1511  Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{z} [#mum];tracks",etaF,etaL,phiF,phiL),
1512  mybins_,-dzmax_eta,dzmax_eta);
1513 
1514  n_dxyBiasResidualsMap[i][j] = NormDoubleDiffBiasRes.make<TH1F>(Form("histo_norm_dxy_eta_plot%i_phi_plot%i",i,j),
1515  Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{xy}/#sigma_{d_{xy}};tracks",etaF,etaL,phiF,phiL),
1516  mybins_,-dzmax_eta/100,dzmax_eta/100);
1517 
1518  n_dzBiasResidualsMap[i][j] = NormDoubleDiffBiasRes.make<TH1F>(Form("histo_norm_dxy_eta_plot%i_phi_plot%i",i,j),
1519  Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{z}/#sigma_{d_{z}};tracks",etaF,etaL,phiF,phiL),
1520  mybins_,-dzmax_eta/100,dzmax_eta/100);
1521  }
1522  }
1523 
1524  // declaration of the directories
1525 
1526  TFileDirectory MeanBiasTrendsDir = fs->mkdir("MeanBiasTrends");
1527  TFileDirectory WidthBiasTrendsDir = fs->mkdir("WidthBiasTrends");
1528  TFileDirectory MedianBiasTrendsDir = fs->mkdir("MedianBiasTrends");
1529  TFileDirectory MADBiasTrendsDir = fs->mkdir("MADBiasTrends");
1530 
1531  TFileDirectory Mean2DBiasMapsDir = fs->mkdir("MeanBiasMaps");
1532  TFileDirectory Width2DBiasMapsDir = fs->mkdir("WidthBiasMaps");
1533 
1534  // means and widths from the fit
1535 
1536  a_dxyPhiMeanBiasTrend = MeanBiasTrendsDir.make<TH1F> ("means_dxy_phi",
1537  "#LT d_{xy} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{xy} #GT [#mum]",
1538  nBins_,lowedge,highedge);
1539 
1540  a_dxyPhiWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>("widths_dxy_phi",
1541  "#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees];#sigma_{d_{xy}} [#mum]",
1542  nBins_,lowedge,highedge);
1543 
1544  a_dzPhiMeanBiasTrend = MeanBiasTrendsDir.make<TH1F> ("means_dz_phi",
1545  "#LT d_{z} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{z} #GT [#mum]",
1546  nBins_,lowedge,highedge);
1547 
1548  a_dzPhiWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>("widths_dz_phi",
1549  "#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];#sigma_{d_{z}} [#mum]",
1550  nBins_,lowedge,highedge);
1551 
1552  a_dxyEtaMeanBiasTrend = MeanBiasTrendsDir.make<TH1F> ("means_dxy_eta",
1553  "#LT d_{xy} #GT vs #eta sector;#eta (sector);#LT d_{xy} #GT [#mum]",
1554  nBins_,lowedge,highedge);
1555 
1556  a_dxyEtaWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>("widths_dxy_eta",
1557  "#sigma_{d_{xy}} vs #eta sector;#eta (sector);#sigma_{d_{xy}} [#mum]",
1558  nBins_,lowedge,highedge);
1559 
1560  a_dzEtaMeanBiasTrend = MeanBiasTrendsDir.make<TH1F> ("means_dz_eta",
1561  "#LT d_{z} #GT vs #eta sector;#eta (sector);#LT d_{z} #GT [#mum]",
1562  nBins_,lowedge,highedge);
1563 
1564  a_dzEtaWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>("widths_dz_eta",
1565  "#sigma_{d_{xy}} vs #eta sector;#eta (sector);#sigma_{d_{z}} [#mum]",
1566  nBins_,lowedge,highedge);
1567 
1568  n_dxyPhiMeanBiasTrend = MeanBiasTrendsDir.make<TH1F> ("norm_means_dxy_phi",
1569  "#LT d_{xy}/#sigma_{d_{xy}} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{xy}/#sigma_{d_{xy}} #GT",
1570  nBins_,lowedge,highedge);
1571 
1572  n_dxyPhiWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>("norm_widths_dxy_phi",
1573  "width(d_{xy}/#sigma_{d_{xy}}) vs #phi sector;#varphi (sector) [degrees]; width(d_{xy}/#sigma_{d_{xy}})",
1574  nBins_,lowedge,highedge);
1575 
1576  n_dzPhiMeanBiasTrend = MeanBiasTrendsDir.make<TH1F> ("norm_means_dz_phi",
1577  "#LT d_{z}/#sigma_{d_{z}} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{z}/#sigma_{d_{z}} #GT",
1578  nBins_,lowedge,highedge);
1579 
1580  n_dzPhiWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>("norm_widths_dz_phi",
1581  "width(d_{z}/#sigma_{d_{z}}) vs #phi sector;#varphi (sector) [degrees];width(d_{z}/#sigma_{d_{z}})",
1582  nBins_,lowedge,highedge);
1583 
1584  n_dxyEtaMeanBiasTrend = MeanBiasTrendsDir.make<TH1F> ("norm_means_dxy_eta",
1585  "#LT d_{xy}/#sigma_{d_{xy}} #GT vs #eta sector;#eta (sector);#LT d_{xy}/#sigma_{d_{z}} #GT",
1586  nBins_,lowedge,highedge);
1587 
1588  n_dxyEtaWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>("norm_widths_dxy_eta",
1589  "width(d_{xy}/#sigma_{d_{xy}}) vs #eta sector;#eta (sector);width(d_{xy}/#sigma_{d_{z}})",
1590  nBins_,lowedge,highedge);
1591 
1592  n_dzEtaMeanBiasTrend = MeanBiasTrendsDir.make<TH1F> ("norm_means_dz_eta",
1593  "#LT d_{z}/#sigma_{d_{z}} #GT vs #eta sector;#eta (sector);#LT d_{z}/#sigma_{d_{z}} #GT",
1594  nBins_,lowedge,highedge);
1595 
1596  n_dzEtaWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>("norm_widths_dz_eta",
1597  "width(d_{z}/#sigma_{d_{z}}) vs #eta sector;#eta (sector);width(d_{z}/#sigma_{d_{z}})",
1598  nBins_,lowedge,highedge);
1599 
1600  // 2D maps
1601 
1602  a_dxyMeanBiasMap = Mean2DBiasMapsDir.make<TH2F> ("means_dxy_map",
1603  "#LT d_{xy} #GT map;#eta (sector);#varphi (sector) [degrees]",
1604  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1605 
1606  a_dzMeanBiasMap = Mean2DBiasMapsDir.make<TH2F> ("means_dz_map",
1607  "#LT d_{z} #GT map;#eta (sector);#varphi (sector) [degrees]",
1608  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1609 
1610  n_dxyMeanBiasMap = Mean2DBiasMapsDir.make<TH2F> ("norm_means_dxy_map",
1611  "#LT d_{xy}/#sigma_{d_{xy}} #GT map;#eta (sector);#varphi (sector) [degrees]",
1612  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1613 
1614  n_dzMeanBiasMap = Mean2DBiasMapsDir.make<TH2F> ("norm_means_dz_map",
1615  "#LT d_{z}/#sigma_{d_{z}} #GT map;#eta (sector);#varphi (sector) [degrees]",
1616  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1617 
1618  a_dxyWidthBiasMap = Width2DBiasMapsDir.make<TH2F> ("widths_dxy_map",
1619  "#sigma_{d_{xy}} map;#eta (sector);#varphi (sector) [degrees]",
1620  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1621 
1622  a_dzWidthBiasMap = Width2DBiasMapsDir.make<TH2F> ("widths_dz_map",
1623  "#sigma_{d_{z}} map;#eta (sector);#varphi (sector) [degrees]",
1624  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1625 
1626  n_dxyWidthBiasMap = Width2DBiasMapsDir.make<TH2F> ("norm_widths_dxy_map",
1627  "width(d_{xy}/#sigma_{d_{xy}}) map;#eta (sector);#varphi (sector) [degrees]",
1628  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1629 
1630  n_dzWidthBiasMap = Width2DBiasMapsDir.make<TH2F> ("norm_widths_dz_map",
1631  "width(d_{z}/#sigma_{d_{z}}) map;#eta (sector);#varphi (sector) [degrees]",
1632  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1633 
1634  // medians and MADs
1635 
1636  a_dxyPhiMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>("medians_dxy_phi",
1637  "Median of d_{xy} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{xy}) [#mum]",
1638  nBins_,lowedge,highedge);
1639 
1640  a_dxyPhiMADBiasTrend = MADBiasTrendsDir.make<TH1F> ("MADs_dxy_phi",
1641  "Median absolute deviation of d_{xy} vs #phi sector;#varphi (sector) [degrees];MAD(d_{xy}) [#mum]",
1642  nBins_,lowedge,highedge);
1643 
1644  a_dzPhiMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>("medians_dz_phi",
1645  "Median of d_{z} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{z}) [#mum]",
1646  nBins_,lowedge,highedge);
1647 
1648  a_dzPhiMADBiasTrend = MADBiasTrendsDir.make<TH1F> ("MADs_dz_phi",
1649  "Median absolute deviation of d_{z} vs #phi sector;#varphi (sector) [degrees];MAD(d_{z}) [#mum]",
1650  nBins_,lowedge,highedge);
1651 
1652  a_dxyEtaMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>("medians_dxy_eta",
1653  "Median of d_{xy} vs #eta sector;#eta (sector);#mu_{1/2}(d_{xy}) [#mum]",
1654  nBins_,lowedge,highedge);
1655 
1656  a_dxyEtaMADBiasTrend = MADBiasTrendsDir.make<TH1F> ("MADs_dxy_eta",
1657  "Median absolute deviation of d_{xy} vs #eta sector;#eta (sector);MAD(d_{xy}) [#mum]",
1658  nBins_,lowedge,highedge);
1659 
1660  a_dzEtaMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>("medians_dz_eta",
1661  "Median of d_{z} vs #eta sector;#eta (sector);#mu_{1/2}(d_{z}) [#mum]",
1662  nBins_,lowedge,highedge);
1663 
1664  a_dzEtaMADBiasTrend = MADBiasTrendsDir.make<TH1F> ("MADs_dz_eta",
1665  "Median absolute deviation of d_{z} vs #eta sector;#eta (sector);MAD(d_{z}) [#mum]",
1666  nBins_,lowedge,highedge);
1667 
1668  n_dxyPhiMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>("norm_medians_dxy_phi",
1669  "Median of d_{xy}/#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{xy}/#sigma_{d_{xy}})",
1670  nBins_,lowedge,highedge);
1671 
1672  n_dxyPhiMADBiasTrend = MADBiasTrendsDir.make<TH1F> ("norm_MADs_dxy_phi",
1673  "Median absolute deviation of d_{xy}/#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees]; MAD(d_{xy}/#sigma_{d_{xy}})",
1674  nBins_,lowedge,highedge);
1675 
1676  n_dzPhiMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>("norm_medians_dz_phi",
1677  "Median of d_{z}/#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{z}/#sigma_{d_{z}})",
1678  nBins_,lowedge,highedge);
1679 
1680  n_dzPhiMADBiasTrend = MADBiasTrendsDir.make<TH1F> ("norm_MADs_dz_phi",
1681  "Median absolute deviation of d_{z}/#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];MAD(d_{z}/#sigma_{d_{z}})",
1682  nBins_,lowedge,highedge);
1683 
1684  n_dxyEtaMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>("norm_medians_dxy_eta",
1685  "Median of d_{xy}/#sigma_{d_{xy}} vs #eta sector;#eta (sector);#mu_{1/2}(d_{xy}/#sigma_{d_{z}})",
1686  nBins_,lowedge,highedge);
1687 
1688  n_dxyEtaMADBiasTrend = MADBiasTrendsDir.make<TH1F> ("norm_MADs_dxy_eta",
1689  "Median absolute deviation of d_{xy}/#sigma_{d_{xy}} vs #eta sector;#eta (sector);MAD(d_{xy}/#sigma_{d_{z}})",
1690  nBins_,lowedge,highedge);
1691 
1692  n_dzEtaMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>("norm_medians_dz_eta",
1693  "Median of d_{z}/#sigma_{d_{z}} vs #eta sector;#eta (sector);#mu_{1/2}(d_{z}/#sigma_{d_{z}})",
1694  nBins_,lowedge,highedge);
1695 
1696  n_dzEtaMADBiasTrend = MADBiasTrendsDir.make<TH1F> ("norm_MADs_dz_eta",
1697  "Median absolute deviation of d_{z}/#sigma_{d_{z}} vs #eta sector;#eta (sector);MAD(d_{z}/#sigma_{d_{z}})",
1698  nBins_,lowedge,highedge);
1699 
1700  }
1701 }
1702 // ------------ method called once each job just after ending the event loop ------------
1703 void PrimaryVertexValidation::endJob()
1704 {
1705 
1706  edm::LogInfo("PrimaryVertexValidation")
1707  <<"######################################\n"
1708  <<"# PrimaryVertexValidation::endJob()\n"
1709  <<"# Number of analyzed events: "<<Nevt_<<"\n"
1710  <<"######################################";
1711 
1712  if(useTracksFromRecoVtx_){
1713 
1714  fillTrendPlot(a_dxyPhiMeanBiasTrend ,a_dxyPhiBiasResiduals,"mean","phi");
1715  fillTrendPlot(a_dxyPhiWidthBiasTrend,a_dxyPhiBiasResiduals,"width","phi");
1716  fillTrendPlot(a_dzPhiMeanBiasTrend ,a_dzPhiBiasResiduals ,"mean","phi");
1717  fillTrendPlot(a_dzPhiWidthBiasTrend ,a_dzPhiBiasResiduals ,"width","phi");
1718 
1719  fillTrendPlot(a_dxyEtaMeanBiasTrend ,a_dxyEtaBiasResiduals,"mean","eta");
1720  fillTrendPlot(a_dxyEtaWidthBiasTrend,a_dxyEtaBiasResiduals,"width","eta");
1721  fillTrendPlot(a_dzEtaMeanBiasTrend ,a_dzEtaBiasResiduals ,"mean","eta");
1722  fillTrendPlot(a_dzEtaWidthBiasTrend ,a_dzEtaBiasResiduals ,"width","eta");
1723 
1724  fillTrendPlot(n_dxyPhiMeanBiasTrend ,n_dxyPhiBiasResiduals,"mean","phi");
1725  fillTrendPlot(n_dxyPhiWidthBiasTrend,n_dxyPhiBiasResiduals,"width","phi");
1726  fillTrendPlot(n_dzPhiMeanBiasTrend ,n_dzPhiBiasResiduals ,"mean","phi");
1727  fillTrendPlot(n_dzPhiWidthBiasTrend ,n_dzPhiBiasResiduals ,"width","phi");
1728 
1729  fillTrendPlot(n_dxyEtaMeanBiasTrend ,n_dxyEtaBiasResiduals,"mean","eta");
1730  fillTrendPlot(n_dxyEtaWidthBiasTrend,n_dxyEtaBiasResiduals,"width","eta");
1731  fillTrendPlot(n_dzEtaMeanBiasTrend ,n_dzEtaBiasResiduals ,"mean","eta");
1732  fillTrendPlot(n_dzEtaWidthBiasTrend ,n_dzEtaBiasResiduals ,"width","eta");
1733 
1734  // medians and MADs
1735 
1736  fillTrendPlot(a_dxyPhiMedianBiasTrend,a_dxyPhiBiasResiduals,"median","phi");
1737  fillTrendPlot(a_dxyPhiMADBiasTrend ,a_dxyPhiBiasResiduals,"mad","phi");
1738  fillTrendPlot(a_dzPhiMedianBiasTrend ,a_dzPhiBiasResiduals ,"median","phi");
1739  fillTrendPlot(a_dzPhiMADBiasTrend ,a_dzPhiBiasResiduals ,"mad","phi");
1740 
1741  fillTrendPlot(a_dxyEtaMedianBiasTrend,a_dxyEtaBiasResiduals,"median","eta");
1742  fillTrendPlot(a_dxyEtaMADBiasTrend ,a_dxyEtaBiasResiduals,"mad","eta");
1743  fillTrendPlot(a_dzEtaMedianBiasTrend ,a_dzEtaBiasResiduals ,"median","eta");
1744  fillTrendPlot(a_dzEtaMADBiasTrend ,a_dzEtaBiasResiduals ,"mad","eta");
1745 
1746  fillTrendPlot(n_dxyPhiMedianBiasTrend,n_dxyPhiBiasResiduals,"median","phi");
1747  fillTrendPlot(n_dxyPhiMADBiasTrend ,n_dxyPhiBiasResiduals,"mad","phi");
1748  fillTrendPlot(n_dzPhiMedianBiasTrend ,n_dzPhiBiasResiduals ,"median","phi");
1749  fillTrendPlot(n_dzPhiMADBiasTrend ,n_dzPhiBiasResiduals ,"mad","phi");
1750 
1751  fillTrendPlot(n_dxyEtaMedianBiasTrend,n_dxyEtaBiasResiduals,"median","eta");
1752  fillTrendPlot(n_dxyEtaMADBiasTrend ,n_dxyEtaBiasResiduals,"mad","eta");
1753  fillTrendPlot(n_dzEtaMedianBiasTrend ,n_dzEtaBiasResiduals ,"median","eta");
1754  fillTrendPlot(n_dzEtaMADBiasTrend ,n_dzEtaBiasResiduals ,"mad","eta");
1755 
1756  // 2d Maps
1757 
1758  fillMap(a_dxyMeanBiasMap ,a_dxyBiasResidualsMap,"mean");
1759  fillMap(a_dxyWidthBiasMap,a_dxyBiasResidualsMap,"width");
1760  fillMap(a_dzMeanBiasMap ,a_dzBiasResidualsMap,"mean");
1761  fillMap(a_dzWidthBiasMap ,a_dzBiasResidualsMap,"width");
1762 
1763  fillMap(n_dxyMeanBiasMap ,n_dxyBiasResidualsMap,"mean");
1764  fillMap(n_dxyWidthBiasMap,n_dxyBiasResidualsMap,"width");
1765  fillMap(n_dzMeanBiasMap ,n_dzBiasResidualsMap,"mean");
1766  fillMap(n_dzWidthBiasMap ,n_dzBiasResidualsMap,"width");
1767 
1768  }
1769 
1770  fillTrendPlot(a_dxyPhiMeanTrend ,a_dxyPhiResiduals,"mean","phi");
1771  fillTrendPlot(a_dxyPhiWidthTrend,a_dxyPhiResiduals,"width","phi");
1772  fillTrendPlot(a_dzPhiMeanTrend ,a_dzPhiResiduals ,"mean","phi");
1773  fillTrendPlot(a_dzPhiWidthTrend ,a_dzPhiResiduals ,"width","phi");
1774 
1775  fillTrendPlot(a_dxyEtaMeanTrend ,a_dxyEtaResiduals,"mean","eta");
1776  fillTrendPlot(a_dxyEtaWidthTrend,a_dxyEtaResiduals,"width","eta");
1777  fillTrendPlot(a_dzEtaMeanTrend ,a_dzEtaResiduals ,"mean","eta");
1778  fillTrendPlot(a_dzEtaWidthTrend ,a_dzEtaResiduals ,"width","eta");
1779 
1780  fillTrendPlot(n_dxyPhiMeanTrend ,n_dxyPhiResiduals,"mean","phi");
1781  fillTrendPlot(n_dxyPhiWidthTrend,n_dxyPhiResiduals,"width","phi");
1782  fillTrendPlot(n_dzPhiMeanTrend ,n_dzPhiResiduals ,"mean","phi");
1783  fillTrendPlot(n_dzPhiWidthTrend ,n_dzPhiResiduals ,"width","phi");
1784 
1785  fillTrendPlot(n_dxyEtaMeanTrend ,n_dxyEtaResiduals,"mean","eta");
1786  fillTrendPlot(n_dxyEtaWidthTrend,n_dxyEtaResiduals,"width","eta");
1787  fillTrendPlot(n_dzEtaMeanTrend ,n_dzEtaResiduals ,"mean","eta");
1788  fillTrendPlot(n_dzEtaWidthTrend ,n_dzEtaResiduals ,"width","eta");
1789 
1790  // medians and MADs
1791 
1792  fillTrendPlot(a_dxyPhiMedianTrend,a_dxyPhiResiduals,"median","phi");
1793  fillTrendPlot(a_dxyPhiMADTrend ,a_dxyPhiResiduals,"mad","phi");
1794  fillTrendPlot(a_dzPhiMedianTrend ,a_dzPhiResiduals ,"median","phi");
1795  fillTrendPlot(a_dzPhiMADTrend ,a_dzPhiResiduals ,"mad","phi");
1796 
1797  fillTrendPlot(a_dxyEtaMedianTrend,a_dxyEtaResiduals,"median","eta");
1798  fillTrendPlot(a_dxyEtaMADTrend ,a_dxyEtaResiduals,"mad","eta");
1799  fillTrendPlot(a_dzEtaMedianTrend ,a_dzEtaResiduals ,"median","eta");
1800  fillTrendPlot(a_dzEtaMADTrend ,a_dzEtaResiduals ,"mad","eta");
1801 
1802  fillTrendPlot(n_dxyPhiMedianTrend,n_dxyPhiResiduals,"median","phi");
1803  fillTrendPlot(n_dxyPhiMADTrend ,n_dxyPhiResiduals,"mad","phi");
1804  fillTrendPlot(n_dzPhiMedianTrend ,n_dzPhiResiduals ,"median","phi");
1805  fillTrendPlot(n_dzPhiMADTrend ,n_dzPhiResiduals ,"mad","phi");
1806 
1807  fillTrendPlot(n_dxyEtaMedianTrend,n_dxyEtaResiduals,"median","eta");
1808  fillTrendPlot(n_dxyEtaMADTrend ,n_dxyEtaResiduals,"mad","eta");
1809  fillTrendPlot(n_dzEtaMedianTrend ,n_dzEtaResiduals ,"median","eta");
1810  fillTrendPlot(n_dzEtaMADTrend ,n_dzEtaResiduals ,"mad","eta");
1811 
1812  // 2D Maps
1813 
1814  fillMap(a_dxyMeanMap ,a_dxyResidualsMap,"mean");
1815  fillMap(a_dxyWidthMap,a_dxyResidualsMap,"width");
1816  fillMap(a_dzMeanMap ,a_dzResidualsMap,"mean");
1817  fillMap(a_dzWidthMap ,a_dzResidualsMap,"width");
1818 
1819  fillMap(n_dxyMeanMap ,n_dxyResidualsMap,"mean");
1820  fillMap(n_dxyWidthMap,n_dxyResidualsMap,"width");
1821  fillMap(n_dzMeanMap ,n_dzResidualsMap,"mean");
1822  fillMap(n_dzWidthMap ,n_dzResidualsMap,"width");
1823 
1824 }
1825 
1826 //*************************************************************
1827 void PrimaryVertexValidation::SetVarToZero()
1828 //*************************************************************
1829 {
1830  nTracks_ = 0;
1831  nClus_ = 0;
1832  nOfflineVertices_=0;
1833  RunNumber_ =0;
1834  LuminosityBlockNumber_=0;
1835  xOfflineVertex_ =-999.;
1836  yOfflineVertex_ =-999.;
1837  zOfflineVertex_ =-999.;
1838  xErrOfflineVertex_=0.;
1839  yErrOfflineVertex_=0.;
1840  zErrOfflineVertex_=0.;
1841  BSx0_ = -999.;
1842  BSy0_ = -999.;
1843  BSz0_ = -999.;
1844  Beamsigmaz_=-999.;
1845  Beamdxdz_=-999.;
1846  BeamWidthX_=-999.;
1847  BeamWidthY_=-999.;
1848  wxy2_=-999.;
1849 
1850  for ( int i=0; i<nMaxtracks_; ++i ) {
1851 
1852  pt_[i] = 0;
1853  p_[i] = 0;
1854  nhits_[i] = 0;
1855  nhits1D_[i] = 0;
1856  nhits2D_[i] = 0;
1857  nhitsBPIX_[i] = 0;
1858  nhitsFPIX_[i] = 0;
1859  nhitsTIB_[i] = 0;
1860  nhitsTID_[i] = 0;
1861  nhitsTOB_[i] = 0;
1862  nhitsTEC_[i] = 0;
1863  isHighPurity_[i] = 0;
1864  eta_[i] = 0;
1865  theta_[i] = 0;
1866  phi_[i] = 0;
1867  chi2_[i] = 0;
1868  chi2ndof_[i] = 0;
1869  charge_[i] = 0;
1870  qoverp_[i] = 0;
1871  dz_[i] = 0;
1872  dxy_[i] = 0;
1873  dzBs_[i] = 0;
1874  dxyBs_[i] = 0;
1875  xPCA_[i] = 0;
1876  yPCA_[i] = 0;
1877  zPCA_[i] = 0;
1878  xUnbiasedVertex_[i] =0;
1879  yUnbiasedVertex_[i] =0;
1880  zUnbiasedVertex_[i] =0;
1881  chi2normUnbiasedVertex_[i]=0;
1882  chi2UnbiasedVertex_[i]=0;
1883  chi2ProbUnbiasedVertex_[i]=0;
1884  DOFUnbiasedVertex_[i]=0;
1885  sumOfWeightsUnbiasedVertex_[i]=0;
1886  tracksUsedForVertexing_[i]=0;
1887  dxyFromMyVertex_[i]=0;
1888  dzFromMyVertex_[i]=0;
1889  d3DFromMyVertex_[i]=0;
1890  dxyErrorFromMyVertex_[i]=0;
1891  dzErrorFromMyVertex_[i]=0;
1892  d3DErrorFromMyVertex_[i]=0;
1893  IPTsigFromMyVertex_[i]=0;
1894  IPLsigFromMyVertex_[i]=0;
1895  IP3DsigFromMyVertex_[i]=0;
1896  hasRecVertex_[i] = 0;
1897  isGoodTrack_[i] = 0;
1898  }
1899 }
1900 
1901 //*************************************************************
1902 std::pair<Double_t,Double_t> PrimaryVertexValidation::getMedian(TH1F *histo)
1903 //*************************************************************
1904 {
1905  Double_t median = 999;
1906  int nbins = histo->GetNbinsX();
1907 
1908  //extract median from histogram
1909  double *x = new double[nbins];
1910  double *y = new double[nbins];
1911  for (int j = 0; j < nbins; j++) {
1912  x[j] = histo->GetBinCenter(j+1);
1913  y[j] = histo->GetBinContent(j+1);
1914  }
1915  median = TMath::Median(nbins, x, y);
1916 
1917  delete[] x; x = 0;
1918  delete[] y; y = 0;
1919 
1920  std::pair<Double_t,Double_t> result;
1921  result = std::make_pair(median,median/TMath::Sqrt(histo->GetEntries()));
1922 
1923  return result;
1924 
1925 }
1926 
1927 //*************************************************************
1928 std::pair<Double_t,Double_t> PrimaryVertexValidation::getMAD(TH1F *histo)
1929 //*************************************************************
1930 {
1931 
1932  int nbins = histo->GetNbinsX();
1933  Double_t median = getMedian(histo).first;
1934  Double_t x_lastBin = histo->GetBinLowEdge(nbins+1);
1935  const char *HistoName =histo->GetName();
1936  TString Finalname = Form("resMed%s",HistoName);
1937  TH1F *newHisto = new TH1F(Finalname,Finalname,nbins,0.,x_lastBin);
1938  Double_t *residuals = new Double_t[nbins];
1939  Double_t *weights = new Double_t[nbins];
1940 
1941  for (int j = 0; j < nbins; j++) {
1942  residuals[j] = TMath::Abs(median - histo->GetBinCenter(j+1));
1943  weights[j]=histo->GetBinContent(j+1);
1944  newHisto->Fill(residuals[j],weights[j]);
1945  }
1946 
1947  Double_t theMAD = (getMedian(newHisto).first)*1.4826;
1948 
1949  delete[] residuals; residuals=0;
1950  delete[] weights; weights=0;
1951  newHisto->Delete("");
1952 
1953  std::pair<Double_t,Double_t> result;
1954  result = std::make_pair(theMAD,theMAD/histo->GetEntries());
1955 
1956  return result;
1957 
1958 }
1959 
1960 //*************************************************************
1961 std::pair<std::pair<Double_t,Double_t>, std::pair<Double_t,Double_t> > PrimaryVertexValidation::fitResiduals(TH1 *hist)
1962 //*************************************************************
1963 {
1964  //float fitResult(9999);
1965  //if (hist->GetEntries() < 20) return ;
1966 
1967  float mean = hist->GetMean();
1968  float sigma = hist->GetRMS();
1969 
1970  TF1 func("tmp", "gaus", mean - 1.5*sigma, mean + 1.5*sigma);
1971  if (0 == hist->Fit(&func,"QNR")) { // N: do not blow up file by storing fit!
1972  mean = func.GetParameter(1);
1973  sigma = func.GetParameter(2);
1974  // second fit: three sigma of first fit around mean of first fit
1975  func.SetRange(mean - 2*sigma, mean + 2*sigma);
1976  // I: integral gives more correct results if binning is too wide
1977  // L: Likelihood can treat empty bins correctly (if hist not weighted...)
1978  if (0 == hist->Fit(&func, "Q0LR")) {
1979  if (hist->GetFunction(func.GetName())) { // Take care that it is later on drawn:
1980  hist->GetFunction(func.GetName())->ResetBit(TF1::kNotDraw);
1981  }
1982  }
1983  }
1984 
1985  float res_mean = func.GetParameter(1);
1986  float res_width = func.GetParameter(2);
1987 
1988  float res_mean_err = func.GetParError(1);
1989  float res_width_err = func.GetParError(2);
1990 
1991  std::pair<Double_t,Double_t> resultM;
1992  std::pair<Double_t,Double_t> resultW;
1993 
1994  resultM = std::make_pair(res_mean,res_mean_err);
1995  resultW = std::make_pair(res_width,res_width_err);
1996 
1997  std::pair<std::pair<Double_t,Double_t>, std::pair<Double_t,Double_t> > result;
1998 
1999  result = std::make_pair(resultM,resultW);
2000  return result;
2001 }
2002 
2003 //*************************************************************
2004 void PrimaryVertexValidation::fillTrendPlot(TH1F* trendPlot, TH1F* residualsPlot[100], TString fitPar_, TString var_)
2005 //*************************************************************
2006 {
2007 
2008  for ( int i=0; i<nBins_; ++i ) {
2009 
2010  char phipositionString[129];
2011  float phiInterval = phiSect_*(180/TMath::Pi());
2012  float phiposition = (-180+i*phiInterval)+(phiInterval/2);
2013  sprintf(phipositionString,"%.f",phiposition);
2014 
2015  char etapositionString[129];
2016  float etaposition = (-2.5+i*etaSect_)+(etaSect_/2);
2017  sprintf(etapositionString,"%.1f",etaposition);
2018 
2019  if(fitPar_=="mean"){
2020  float mean_ = fitResiduals(residualsPlot[i]).first.first;
2021  float meanErr_ = fitResiduals(residualsPlot[i]).first.second;
2022  trendPlot->SetBinContent(i+1,mean_);
2023  trendPlot->SetBinError(i+1,meanErr_);
2024  } else if (fitPar_=="width"){
2025  float width_ = fitResiduals(residualsPlot[i]).second.first;
2026  float widthErr_ = fitResiduals(residualsPlot[i]).second.second;
2027  trendPlot->SetBinContent(i+1,width_);
2028  trendPlot->SetBinError(i+1,widthErr_);
2029  } else if (fitPar_=="median"){
2030  float median_ = getMedian(residualsPlot[i]).first;
2031  float medianErr_ = getMedian(residualsPlot[i]).second;
2032  trendPlot->SetBinContent(i+1,median_);
2033  trendPlot->SetBinError(i+1,medianErr_);
2034  } else if (fitPar_=="mad"){
2035  float mad_ = getMAD(residualsPlot[i]).first;
2036  float madErr_ = getMAD(residualsPlot[i]).second;
2037  trendPlot->SetBinContent(i+1,mad_);
2038  trendPlot->SetBinError(i+1,madErr_);
2039  } else {
2040  std::cout<<"PrimaryVertexValidation::fillTrendPlot() "<<fitPar_<<" unknown estimator!"<<std::endl;
2041  }
2042 
2043  if(var_=="eta"){
2044  trendPlot->GetXaxis()->SetBinLabel(i+1,etapositionString);
2045  } else if(var_=="phi"){
2046  trendPlot->GetXaxis()->SetBinLabel(i+1,phipositionString);
2047  } else {
2048  std::cout<<"PrimaryVertexValidation::fillTrendPlot() "<<var_<<" unknown track parameter!"<<std::endl;
2049  }
2050  }
2051 }
2052 
2053 //*************************************************************
2054 void PrimaryVertexValidation::fillMap(TH2F* trendMap, TH1F* residualsMapPlot[100][100], TString fitPar_)
2055 //*************************************************************
2056 {
2057 
2058  for ( int i=0; i<nBins_; ++i ) {
2059 
2060  char phipositionString[129];
2061  float phiInterval = phiSect_*(180/TMath::Pi());
2062  float phiposition = (-180+i*phiInterval)+(phiInterval/2);
2063  sprintf(phipositionString,"%.f",phiposition);
2064 
2065  trendMap->GetYaxis()->SetBinLabel(i+1,phipositionString);
2066 
2067  for ( int j=0; j<nBins_; ++j ) {
2068 
2069  char etapositionString[129];
2070  float etaposition = (-2.5+j*etaSect_)+(etaSect_/2);
2071  sprintf(etapositionString,"%.1f",etaposition);
2072 
2073  if(i==0) { trendMap->GetXaxis()->SetBinLabel(j+1,etapositionString); }
2074 
2075  if(fitPar_=="mean"){
2076  float mean_ = fitResiduals(residualsMapPlot[i][j]).first.first;
2077  float meanErr_ = fitResiduals(residualsMapPlot[i][j]).first.second;
2078  trendMap->SetBinContent(j+1,i+1,mean_);
2079  trendMap->SetBinError(j+1,i+1,meanErr_);
2080  } else if (fitPar_=="width"){
2081  float width_ = fitResiduals(residualsMapPlot[i][j]).second.first;
2082  float widthErr_ = fitResiduals(residualsMapPlot[i][j]).second.second;
2083  trendMap->SetBinContent(j+1,i+1,width_);
2084  trendMap->SetBinError(j+1,i+1,widthErr_);
2085  } else if (fitPar_=="median"){
2086  float median_ = getMedian(residualsMapPlot[i][j]).first;
2087  float medianErr_ = getMedian(residualsMapPlot[i][j]).second;
2088  trendMap->SetBinContent(j+1,i+1,median_);
2089  trendMap->SetBinError(j+1,i+1,medianErr_);
2090  } else if (fitPar_=="mad"){
2091  float mad_ = getMAD(residualsMapPlot[i][j]).first;
2092  float madErr_ = getMAD(residualsMapPlot[i][j]).second;
2093  trendMap->SetBinContent(j+1,i+1,mad_);
2094  trendMap->SetBinError(j+1,i+1,madErr_);
2095  } else {
2096  std::cout<<"PrimaryVertexValidation::fillMap() "<<fitPar_<<" unknown estimator!"<<std::endl;
2097  }
2098  } // closes loop on eta bins
2099  } // cloeses loop on phi bins
2100 }
2101 
2102 //*************************************************************
2103 bool PrimaryVertexValidation::vtxSort( const reco::Vertex & a, const reco::Vertex & b )
2104 //*************************************************************
2105 {
2106  if( a.tracksSize() != b.tracksSize() )
2107  return a.tracksSize() > b.tracksSize() ? true : false ;
2108  else
2109  return a.chi2() < b.chi2() ? true : false ;
2110 }
2111 
2112 //*************************************************************
2113 bool PrimaryVertexValidation::passesTrackCuts(const reco::Track & track, const reco::Vertex & vertex,std::string qualityString_, double dxyErrMax_,double dzErrMax_, double ptErrMax_)
2114 //*************************************************************
2115 {
2116 
2117  math::XYZPoint vtxPoint(0.0,0.0,0.0);
2118  double vzErr =0.0, vxErr=0.0, vyErr=0.0;
2119  vtxPoint=vertex.position();
2120  vzErr=vertex.zError();
2121  vxErr=vertex.xError();
2122  vyErr=vertex.yError();
2123 
2124  double dxy=0.0, dz=0.0, dxysigma=0.0, dzsigma=0.0;
2125  dxy = track.dxy(vtxPoint);
2126  dz = track.dz(vtxPoint);
2127  dxysigma = sqrt(track.d0Error()*track.d0Error()+vxErr*vyErr);
2128  dzsigma = sqrt(track.dzError()*track.dzError()+vzErr*vzErr);
2129 
2130  if(track.quality(reco::TrackBase::qualityByName(qualityString_)) != 1)return false;
2131  if(fabs(dxy/dxysigma) > dxyErrMax_) return false;
2132  if(fabs(dz/dzsigma) > dzErrMax_) return false;
2133  if(track.ptError() / track.pt() > ptErrMax_) return false;
2134 
2135  return true;
2136 }
2137 
2138 
2139 //*************************************************************
2140 std::map<std::string, TH1*> PrimaryVertexValidation::bookVertexHistograms(TFileDirectory dir)
2141 //*************************************************************
2142 {
2143 
2144  TH1F::SetDefaultSumw2(kTRUE);
2145 
2146  std::map<std::string, TH1*> h;
2147 
2148  // histograms of track quality (Data and MC)
2149  std::string types[] = {"all","sel"};
2150  for(int t=0; t<2; t++){
2151  h["pseudorapidity_"+types[t]] =dir.make <TH1F>(("rapidity_"+types[t]).c_str(),"track pseudorapidity; track #eta; tracks",100,-3., 3.);
2152  h["z0_"+types[t]] = dir.make<TH1F>(("z0_"+types[t]).c_str(),"track z_{0};track z_{0} (cm);tracks",80,-40., 40.);
2153  h["phi_"+types[t]] = dir.make<TH1F>(("phi_"+types[t]).c_str(),"track #phi; track #phi;tracks",80,-TMath::Pi(), TMath::Pi());
2154  h["eta_"+types[t]] = dir.make<TH1F>(("eta_"+types[t]).c_str(),"track #eta; track #eta;tracks",80,-4., 4.);
2155  h["pt_"+types[t]] = dir.make<TH1F>(("pt_"+types[t]).c_str(),"track p_{T}; track p_{T} [GeV];tracks",100,0., 20.);
2156  h["p_"+types[t]] = dir.make<TH1F>(("p_"+types[t]).c_str(),"track p; track p [GeV];tracks",100,0., 20.);
2157  h["found_"+types[t]] = dir.make<TH1F>(("found_"+types[t]).c_str(),"n. found hits;n^{found}_{hits};tracks",30, 0., 30.);
2158  h["lost_"+types[t]] = dir.make<TH1F>(("lost_"+types[t]).c_str(),"n. lost hits;n^{lost}_{hits};tracks",20, 0., 20.);
2159  h["nchi2_"+types[t]] = dir.make<TH1F>(("nchi2_"+types[t]).c_str(),"normalized track #chi^{2};track #chi^{2}/ndf;tracks",100, 0., 20.);
2160  h["rstart_"+types[t]] = dir.make<TH1F>(("rstart_"+types[t]).c_str(),"track start radius; track innermost radius r (cm);tracks",100, 0., 20.);
2161  h["expectedInner_"+types[t]] = dir.make<TH1F>(("expectedInner_"+types[t]).c_str(),"n. expected inner hits;n^{expected}_{inner};tracks",10, 0., 10.);
2162  h["expectedOuter_"+types[t]] = dir.make<TH1F>(("expectedOuter_"+types[t]).c_str(),"n. expected outer hits;n^{expected}_{outer};tracks ",10, 0., 10.);
2163  h["logtresxy_"+types[t]] = dir.make<TH1F>(("logtresxy_"+types[t]).c_str(),"log10(track r-#phi resolution/#mum);log10(track r-#phi resolution/#mum);tracks",100, 0., 5.);
2164  h["logtresz_"+types[t]] = dir.make<TH1F>(("logtresz_"+types[t]).c_str(),"log10(track z resolution/#mum);log10(track z resolution/#mum);tracks",100, 0., 5.);
2165  h["tpullxy_"+types[t]] = dir.make<TH1F>(("tpullxy_"+types[t]).c_str(),"track r-#phi pull;pull_{r-#phi};tracks",100, -10., 10.);
2166  h["tpullz_"+types[t]] = dir.make<TH1F>(("tpullz_"+types[t]).c_str(),"track r-z pull;pull_{r-z};tracks",100, -50., 50.);
2167  h["tlogDCAxy_"+types[t]] = dir.make<TH1F>(("tlogDCAxy_"+types[t]).c_str(),"track log_{10}(DCA_{r-#phi});track log_{10}(DCA_{r-#phi});tracks",200, -5., 3.);
2168  h["tlogDCAz_"+types[t]] = dir.make<TH1F>(("tlogDCAz_"+types[t]).c_str(),"track log_{10}(DCA_{r-z});track log_{10}(DCA_{r-z});tracks",200, -5., 5.);
2169  h["lvseta_"+types[t]] = dir.make<TH2F>(("lvseta_"+types[t]).c_str(),"cluster length vs #eta;track #eta;cluster length",60,-3., 3., 20, 0., 20);
2170  h["lvstanlambda_"+types[t]] = dir.make<TH2F>(("lvstanlambda_"+types[t]).c_str(),"cluster length vs tan #lambda; tan#lambda;cluster length",60,-6., 6., 20, 0., 20);
2171  h["restrkz_"+types[t]] = dir.make<TH1F>(("restrkz_"+types[t]).c_str(),"z-residuals (track vs vertex);res_{z} (cm);tracks", 200, -5., 5.);
2172  h["restrkzvsphi_"+types[t]] = dir.make<TH2F>(("restrkzvsphi_"+types[t]).c_str(),"z-residuals (track - vertex) vs track #phi;track #phi;res_{z} (cm)", 12,-TMath::Pi(),TMath::Pi(),100, -0.5,0.5);
2173  h["restrkzvseta_"+types[t]] = dir.make<TH2F>(("restrkzvseta_"+types[t]).c_str(),"z-residuals (track - vertex) vs track #eta;track #eta;res_{z} (cm)", 12,-3.,3.,200, -0.5,0.5);
2174  h["pulltrkzvsphi_"+types[t]] = dir.make<TH2F>(("pulltrkzvsphi_"+types[t]).c_str(),"normalized z-residuals (track - vertex) vs track #phi;track #phi;res_{z}/#sigma_{res_{z}}", 12,-TMath::Pi(),TMath::Pi(),100, -5., 5.);
2175  h["pulltrkzvseta_"+types[t]] = dir.make<TH2F>(("pulltrkzvseta_"+types[t]).c_str(),"normalized z-residuals (track - vertex) vs track #eta;track #eta;res_{z}/#sigma_{res_{z}}", 12,-3.,3.,100, -5., 5.);
2176  h["pulltrkz_"+types[t]] = dir.make<TH1F>(("pulltrkz_"+types[t]).c_str(),"normalized z-residuals (track vs vertex);res_{z}/#sigma_{res_{z}};tracks", 100, -5., 5.);
2177  h["sigmatrkz0_"+types[t]] = dir.make<TH1F>(("sigmatrkz0_"+types[t]).c_str(),"z-resolution (excluding beam);#sigma^{trk}_{z_{0}} (cm);tracks", 100, 0., 5.);
2178  h["sigmatrkz_"+types[t]] = dir.make<TH1F>(("sigmatrkz_"+types[t]).c_str(),"z-resolution (including beam);#sigma^{trk}_{z} (cm);tracks", 100,0., 5.);
2179  h["nbarrelhits_"+types[t]] = dir.make<TH1F>(("nbarrelhits_"+types[t]).c_str(),"number of pixel barrel hits;n. hits Barrel Pixel;tracks", 10, 0., 10.);
2180  h["nbarrelLayers_"+types[t]] = dir.make<TH1F>(("nbarrelLayers_"+types[t]).c_str(),"number of pixel barrel layers;n. layers Barrel Pixel;tracks", 10, 0., 10.);
2181  h["nPxLayers_"+types[t]] = dir.make<TH1F>(("nPxLayers_"+types[t]).c_str(),"number of pixel layers (barrel+endcap);n. Pixel layers;tracks", 10, 0., 10.);
2182  h["nSiLayers_"+types[t]] = dir.make<TH1F>(("nSiLayers_"+types[t]).c_str(),"number of Tracker layers;n. Tracker layers;tracks", 20, 0., 20.);
2183  h["trackAlgo_"+types[t]] = dir.make<TH1F>(("trackAlgo_"+types[t]).c_str(),"track algorithm;track algo;tracks", 30, 0., 30.);
2184  h["trackQuality_"+types[t]] = dir.make<TH1F>(("trackQuality_"+types[t]).c_str(),"track quality;track quality;tracks", 7, -1., 6.);
2185  }
2186 
2187  return h;
2188 
2189 }
2190 
2191 //*************************************************************
2192 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_)
2193 //*************************************************************
2194 {
2195 
2196  using namespace reco;
2197 
2198  fill(h,"pseudorapidity_"+ttype,tt->track().eta());
2199  fill(h,"z0_"+ttype,tt->track().vz());
2200  fill(h,"phi_"+ttype,tt->track().phi());
2201  fill(h,"eta_"+ttype,tt->track().eta());
2202  fill(h,"pt_"+ttype,tt->track().pt());
2203  fill(h,"p_"+ttype,tt->track().p());
2204  fill(h,"found_"+ttype,tt->track().found());
2205  fill(h,"lost_"+ttype,tt->track().lost());
2206  fill(h,"nchi2_"+ttype,tt->track().normalizedChi2());
2207  fill(h,"rstart_"+ttype,(tt->track().innerPosition()).Rho());
2208 
2209  double d0Error=tt->track().d0Error();
2210  double d0=tt->track().dxy(beamSpot.position());
2211  double dz=tt->track().dz(beamSpot.position());
2212  if (d0Error>0){
2213  fill(h,"logtresxy_"+ttype,log(d0Error/0.0001)/log(10.));
2214  fill(h,"tpullxy_"+ttype,d0/d0Error);
2215  fill(h,"tlogDCAxy_"+ttype,log(fabs(d0/d0Error)));
2216 
2217  }
2218  //double z0=tt->track().vz();
2219  double dzError=tt->track().dzError();
2220  if(dzError>0){
2221  fill(h,"logtresz_"+ttype,log(dzError/0.0001)/log(10.));
2222  fill(h,"tpullz_"+ttype,dz/dzError);
2223  fill(h,"tlogDCAz_"+ttype,log(fabs(dz/dzError)));
2224  }
2225 
2226  //
2227  double wxy2_=pow(beamSpot.BeamWidthX(),2)+pow(beamSpot.BeamWidthY(),2);
2228 
2229  fill(h,"sigmatrkz_"+ttype,sqrt(pow(tt->track().dzError(),2)+wxy2_/pow(tan(tt->track().theta()),2)));
2230  fill(h,"sigmatrkz0_"+ttype,tt->track().dzError());
2231 
2232  // track vs vertex
2233  if( v.isValid()){ // && (v.ndof()<10.)) {
2234  // emulate clusterizer input
2235  //const TransientTrack & tt = theB_->build(&t); wrong !!!!
2236  //reco::TransientTrack tt = theB_->build(&t);
2237  //ttt->track().setBeamSpot(beamSpot); // need the setBeamSpot !
2238  double z=(tt->stateAtBeamLine().trackStateAtPCA()).position().z();
2239  double tantheta=tan((tt->stateAtBeamLine().trackStateAtPCA()).momentum().theta());
2240  double dz2= pow(tt->track().dzError(),2)+wxy2_/pow(tantheta,2);
2241 
2242  fill(h,"restrkz_"+ttype,z-v.position().z());
2243  fill(h,"restrkzvsphi_"+ttype,tt->track().phi(), z-v.position().z());
2244  fill(h,"restrkzvseta_"+ttype,tt->track().eta(), z-v.position().z());
2245  fill(h,"pulltrkzvsphi_"+ttype,tt->track().phi(), (z-v.position().z())/sqrt(dz2));
2246  fill(h,"pulltrkzvseta_"+ttype,tt->track().eta(), (z-v.position().z())/sqrt(dz2));
2247 
2248  fill(h,"pulltrkz_"+ttype,(z-v.position().z())/sqrt(dz2));
2249 
2250  double x1=tt->track().vx()-beamSpot.x0(); double y1=tt->track().vy()-beamSpot.y0();
2251 
2252  double kappa=-0.002998*fBfield_*tt->track().qoverp()/cos(tt->track().theta());
2253  double D0=x1*sin(tt->track().phi())-y1*cos(tt->track().phi())-0.5*kappa*(x1*x1+y1*y1);
2254  double q=sqrt(1.-2.*kappa*D0);
2255  double s0=(x1*cos(tt->track().phi())+y1*sin(tt->track().phi()))/q;
2256  // double s1;
2257  if (fabs(kappa*s0)>0.001){
2258  //s1=asin(kappa*s0)/kappa;
2259  }else{
2260  //double ks02=(kappa*s0)*(kappa*s0);
2261  //s1=s0*(1.+ks02/6.+3./40.*ks02*ks02+5./112.*pow(ks02,3));
2262  }
2263  // sp.ddcap=-2.*D0/(1.+q);
2264  //double zdcap=tt->track().vz()-s1/tan(tt->track().theta());
2265 
2266  }
2267  //
2268 
2269  // collect some info on hits and clusters
2270  fill(h,"nbarrelLayers_"+ttype,tt->track().hitPattern().pixelBarrelLayersWithMeasurement());
2271  fill(h,"nPxLayers_"+ttype,tt->track().hitPattern().pixelLayersWithMeasurement());
2272  fill(h,"nSiLayers_"+ttype,tt->track().hitPattern().trackerLayersWithMeasurement());
2273  fill(h,"expectedInner_"+ttype,tt->track().hitPattern().numberOfHits(HitPattern::MISSING_INNER_HITS));
2274  fill(h,"expectedOuter_"+ttype,tt->track().hitPattern().numberOfHits(HitPattern::MISSING_OUTER_HITS));
2275  fill(h,"trackAlgo_"+ttype,tt->track().algo());
2276  fill(h,"trackQuality_"+ttype,tt->track().qualityMask());
2277 
2278  //
2279  int longesthit=0, nbarrel=0;
2280  for(trackingRecHit_iterator hit=tt->track().recHitsBegin(); hit!=tt->track().recHitsEnd(); hit++){
2281  if ((**hit).isValid() && (**hit).geographicalId().det() == DetId::Tracker ){
2282  bool barrel = DetId((**hit).geographicalId()).subdetId() == static_cast<int>(PixelSubdetector::PixelBarrel);
2283  //bool endcap = DetId::DetId((**hit).geographicalId()).subdetId() == static_cast<int>(PixelSubdetector::PixelEndcap);
2284  if (barrel){
2285  const SiPixelRecHit *pixhit = dynamic_cast<const SiPixelRecHit*>( &(**hit));
2286  edm::Ref<edmNew::DetSetVector<SiPixelCluster>, SiPixelCluster> const& clust = (*pixhit).cluster();
2287  if (clust.isNonnull()) {
2288  nbarrel++;
2289  if (clust->sizeY()-longesthit>0) longesthit=clust->sizeY();
2290  if (clust->sizeY()>20.){
2291  fill(h,"lvseta_"+ttype,tt->track().eta(), 19.9);
2292  fill(h,"lvstanlambda_"+ttype,tan(tt->track().lambda()), 19.9);
2293  }else{
2294  fill(h,"lvseta_"+ttype,tt->track().eta(), float(clust->sizeY()));
2295  fill(h,"lvstanlambda_"+ttype,tan(tt->track().lambda()), float(clust->sizeY()));
2296  }
2297  }
2298  }
2299  }
2300  }
2301  fill(h,"nbarrelhits_"+ttype,float(nbarrel));
2302  //-------------------------------------------------------------------
2303 }
2304 
2305 //*************************************************************
2306 void PrimaryVertexValidation::add(std::map<std::string, TH1*>& h, TH1* hist)
2307 //*************************************************************
2308 {
2309  h[hist->GetName()]=hist;
2310  hist->StatOverflows(kTRUE);
2311 }
2312 
2313 //*************************************************************
2314 void PrimaryVertexValidation::fill(std::map<std::string, TH1*>& h, std::string s, double x)
2315 //*************************************************************
2316 {
2317  // cout << "fill1 " << s << endl;
2318  if(h.count(s)==0){
2319  std::cout << "Trying to fill non-exiting Histogram named " << s << std::endl;
2320  return;
2321  }
2322  h[s]->Fill(x);
2323 }
2324 
2325 //*************************************************************
2326 void PrimaryVertexValidation::fill(std::map<std::string, TH1*>& h, std::string s, double x, double y)
2327 //*************************************************************
2328 {
2329  // cout << "fill2 " << s << endl;
2330  if(h.count(s)==0){
2331  std::cout << "Trying to fill non-exiting Histogram named " << s << std::endl;
2332  return;
2333  }
2334  h[s]->Fill(x,y);
2335 }
2336 
2337 //define this as a plug-in
2338 DEFINE_FWK_MODULE(PrimaryVertexValidation);
reco::TrackBase::qoverp
double qoverp() const
q / p
Definition: TrackBase.h:568
cms::Exception::what
virtual char const * what() const
Definition: Exception.cc:141
TFileService::kSharedResource
static const std::string kSharedResource
Definition: TFileService.h:76
TransientVertex::positionError
GlobalError positionError() const
Definition: TransientVertex.h:115
reco::TrackBase::p
double p() const
momentum vector magnitude
Definition: TrackBase.h:610
Pi
const double Pi
Definition: CosmicMuonParameters.h:18
type
type
Definition: HCALResponse.h:21
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
edm::EventID::event
EventNumber_t event() const
Definition: EventID.h:41
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
i
int i
Definition: DBlmapReader.cc:9
TrackingRecHit::dimension
virtual int dimension() const =0
PrimaryVertexValidation::n_dzResidualsMap
TH1F * n_dzResidualsMap[nMaxBins_][nMaxBins_]
Definition: PrimaryVertexValidation.h:259
reco::TrackBase::d0Error
double d0Error() const
error on d0
Definition: TrackBase.h:797
PrimaryVertexValidation::a_IP3DEtaResiduals
TH1F * a_IP3DEtaResiduals[nMaxBins_]
Definition: PrimaryVertexValidation.h:224
PrimaryVertexValidation::runControlNumbers_
std::vector< unsigned int > runControlNumbers_
Definition: PrimaryVertexValidation.h:112
SiPixelRawToDigiRegional_cfi.beamSpot
tuple beamSpot
Definition: SiPixelRawToDigiRegional_cfi.py:11
edm::Ref::isNonnull
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:252
reco::Track::lost
unsigned short lost() const
Number of lost (=invalid) hits on track.
Definition: Track.h:199
h
FWCore Framework interface EventSetupRecordImplementation h
Helper function to determine trigger accepts.
Definition: L1TUtmAlgorithmRcd.h:4
PrimaryVertexValidation::vtxSort
static bool vtxSort(const reco::Vertex &a, const reco::Vertex &b)
Definition: PrimaryVertexValidation.cc:2103
Reference_intrackfit_cff.barrel
barrel
Definition: Reference_intrackfit_cff.py:37
reco::Vertex::tracks_end
trackRef_iterator tracks_end() const
last iterator over tracks
Definition: Vertex.cc:44
HLT_25ns10e33_v2_cff.vertices
tuple vertices
Definition: HLT_25ns10e33_v2_cff.py:6505
reco::TransientTrack::setBeamSpot
void setBeamSpot(const reco::BeamSpot &beamSpot)
Definition: TransientTrack.h:69
reco::TrackBase::normalizedChi2
double normalizedChi2() const
chi-squared divided by n.d.o.f. (or chi-squared * 1e6 if n.d.o.f. is zero)
Definition: TrackBase.h:556
edm::Event::getByToken
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:462
reco::Vertex::zError
double zError() const
error on z
Definition: Vertex.h:111
PrimaryVertexValidation::n_IP2DPhiResiduals
TH1F * n_IP2DPhiResiduals[nMaxBins_]
Definition: PrimaryVertexValidation.h:240
PrimaryVertexValidation::theTrackFilter_
TrackFilterForPVFindingBase * theTrackFilter_
Definition: PrimaryVertexValidation.h:93
PrimaryVertexValidation::a_dxyResidualsMap
TH1F * a_dxyResidualsMap[nMaxBins_][nMaxBins_]
Definition: PrimaryVertexValidation.h:254
PrimaryVertexValidation::hDA
std::map< std::string, TH1 * > hDA
Definition: PrimaryVertexValidation.h:496
reco::TrackBase::TrackQuality
TrackQuality
track quality
Definition: TrackBase.h:151
PrimaryVertexValidation::n_dxyEtaBiasResiduals
TH1F * n_dxyEtaBiasResiduals[nMaxBins_]
Definition: PrimaryVertexValidation.h:347
VertexException
Common base class.
Definition: VertexException.h:12
reco::TrackBase::theta
double theta() const
polar angle
Definition: TrackBase.h:574
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
reco::Vertex::isValid
bool isValid() const
Tells whether the vertex is valid.
Definition: Vertex.h:60
relval_2017.k
list k
Definition: relval_2017.py:23
funct::D0
Divides< arg, void > D0
Definition: Factorize.h:143
PrimaryVertexValidation::PrimaryVertexValidation
PrimaryVertexValidation(const edm::ParameterSet &)
Definition: PrimaryVertexValidation.cc:66
PrimaryVertexValidation::getMAD
std::pair< Double_t, Double_t > getMAD(TH1F *histo)
Definition: PrimaryVertexValidation.cc:1928
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
IPTools::signedTransverseImpactParameter
std::pair< bool, Measurement1D > signedTransverseImpactParameter(const reco::TransientTrack &track, const GlobalVector &direction, const reco::Vertex &vertex)
Definition: IPTools.cc:50
PrimaryVertexValidation::n_dxyBiasResidualsMap
TH1F * n_dxyBiasResidualsMap[nMaxBins_][nMaxBins_]
Definition: PrimaryVertexValidation.h:357
TransientVertex::totalChiSquared
float totalChiSquared() const
Definition: TransientVertex.h:125
reco::TransientTrack::hitPattern
const HitPattern & hitPattern() const
Definition: TransientTrack.h:128
GlobalPoint
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
IPTools::absoluteImpactParameter3D
std::pair< bool, Measurement1D > absoluteImpactParameter3D(const reco::TransientTrack &transientTrack, const reco::Vertex &vertex)
Definition: IPTools.cc:37
PV3DBase::y
T y() const
Definition: PV3DBase.h:63
edm::EventAuxiliary::run
RunNumber_t run() const
Definition: EventAuxiliary.h:67
IPTools::signedImpactParameter3D
std::pair< bool, Measurement1D > signedImpactParameter3D(const reco::TransientTrack &track, const GlobalVector &direction, const reco::Vertex &vertex)
Definition: IPTools.cc:71
reco::TrackBase::phi
double phi() const
azimuthal angle of momentum vector
Definition: TrackBase.h:640
PrimaryVertexValidation::DOFUnbiasedVertex_
float DOFUnbiasedVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:191
PrimaryVertexValidation::a_dzPhiBiasResiduals
TH1F * a_dzPhiBiasResiduals[nMaxBins_]
Definition: PrimaryVertexValidation.h:341
PrimaryVertexValidation::passesTrackCuts
bool passesTrackCuts(const reco::Track &track, const reco::Vertex &vertex, std::string qualityString_, double dxyErrMax_, double dzErrMax_, double ptErrMax_)
Definition: PrimaryVertexValidation.cc:2113
TFileService::make
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
PrimaryVertexValidation::d3DFromMyVertex_
double d3DFromMyVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:197
PrimaryVertexValidation::n_dzEtaBiasResiduals
TH1F * n_dzEtaBiasResiduals[nMaxBins_]
Definition: PrimaryVertexValidation.h:350
PrimaryVertexValidation::analyze
virtual void analyze(const edm::Event &, const edm::EventSetup &)
Definition: PrimaryVertexValidation.cc:127
PrimaryVertexValidation::chi2ProbUnbiasedVertex_
float chi2ProbUnbiasedVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:190
reco::TrackBase::px
double px() const
x coordinate of momentum vector
Definition: TrackBase.h:622
reco::HitPattern::pixelLayersWithMeasurement
int pixelLayersWithMeasurement() const
Definition: HitPattern.cc:508
PrimaryVertexValidation::dxyFromMyVertex_
double dxyFromMyVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:195
dt_dqm_sourceclient_common_cff.reco
tuple reco
Definition: dt_dqm_sourceclient_common_cff.py:107
RecoJets_EventContent_cff.func
tuple func
Definition: RecoJets_EventContent_cff.py:145
reco::HitPattern::trackerLayersWithMeasurement
int trackerLayersWithMeasurement() const
Definition: HitPattern.cc:527
reco::Vertex::position
const Point & position() const
position
Definition: Vertex.h:99
reco::TransientTrack::stateAtBeamLine
TrajectoryStateClosestToBeamLine stateAtBeamLine() const
Definition: TransientTrack.h:114
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
PrimaryVertexValidation::dzErrorFromMyVertex_
double dzErrorFromMyVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:200
PrimaryVertexValidation::fillTrackHistos
void fillTrackHistos(std::map< std::string, TH1 * > &h, const std::string &ttype, const reco::TransientTrack *tt, const reco::Vertex &v, const reco::BeamSpot &beamSpot, double fBfield)
Definition: PrimaryVertexValidation.cc:2192
edm::EventAuxiliary::luminosityBlock
LuminosityBlockNumber_t luminosityBlock() const
Definition: EventAuxiliary.h:63
PrimaryVertexValidation::a_dxyPhiBiasResiduals
TH1F * a_dxyPhiBiasResiduals[nMaxBins_]
Definition: PrimaryVertexValidation.h:338
mps_fire.result
tuple result
Definition: mps_fire.py:84
PrimaryVertexValidation::IPTsigFromMyVertex_
double IPTsigFromMyVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:203
reco::Track::innerPosition
const math::XYZPoint & innerPosition() const
position of the innermost hit
Definition: Track.h:55
reco::TrackBase::algo
TrackAlgorithm algo() const
Definition: TrackBase.h:492
TrackClusterizerInZ::clusterize
virtual std::vector< std::vector< reco::TransientTrack > > clusterize(const std::vector< reco::TransientTrack > &tracks) const =0
PrimaryVertexValidation::getMedian
std::pair< Double_t, Double_t > getMedian(TH1F *histo)
Definition: PrimaryVertexValidation.cc:1902
kFPIX
const int kFPIX
Definition: AlignmentTrackSelector.cc:29
HLT_25ns10e33_v2_cff.InputTag
tuple InputTag
Definition: HLT_25ns10e33_v2_cff.py:58467
PrimaryVertexValidation::fillMap
void fillMap(TH2F *trendMap, TH1F *residualsMapPlot[100][100], TString fitPar_)
Definition: PrimaryVertexValidation.cc:2054
PrimaryVertexValidation::a_IP2DEtaResiduals
TH1F * a_IP2DEtaResiduals[nMaxBins_]
Definition: PrimaryVertexValidation.h:221
TransientVertex::normalisedChiSquared
float normalisedChiSquared() const
Definition: TransientVertex.h:126
reco::TrackBase::vertex
const Point & vertex() const
reference point on the track. This method is DEPRECATED, please use referencePoint() instead ...
Definition: TrackBase.h:682
PrimaryVertexValidation::tracksUsedForVertexing_
int tracksUsedForVertexing_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:193
iEvent
int iEvent
Definition: GenABIO.cc:230
PrimaryVertexValidation::n_IP2DEtaResiduals
TH1F * n_IP2DEtaResiduals[nMaxBins_]
Definition: PrimaryVertexValidation.h:241
reco::TrackBase::eta
double eta() const
pseudorapidity of momentum vector
Definition: TrackBase.h:646
TrajectoryStateClosestToPoint::perigeeParameters
const PerigeeTrajectoryParameters & perigeeParameters() const
Definition: TrajectoryStateClosestToPoint.h:68
PrimaryVertexValidation::hasFirstLayerPixelHits
bool hasFirstLayerPixelHits(const reco::TransientTrack track)
Definition: PrimaryVertexValidation.cc:828
PrimaryVertexValidation::n_reszEtaResiduals
TH1F * n_reszEtaResiduals[nMaxBins_]
Definition: PrimaryVertexValidation.h:247
VertexFitter::vertex
virtual CachingVertex< N > vertex(const std::vector< reco::TransientTrack > &tracks) const =0
reco::TrackBase::chi2
double chi2() const
chi-squared of the fit
Definition: TrackBase.h:544
PrimaryVertexValidation::n_dzPhiBiasResiduals
TH1F * n_dzPhiBiasResiduals[nMaxBins_]
Definition: PrimaryVertexValidation.h:349
TransientVertex::degreesOfFreedom
float degreesOfFreedom() const
Definition: TransientVertex.h:129
PrimaryVertexValidation::yUnbiasedVertex_
double yUnbiasedVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:186
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
PrimaryVertexValidation::fitResiduals
std::pair< std::pair< Double_t, Double_t >, std::pair< Double_t, Double_t > > fitResiduals(TH1 *hist)
Definition: PrimaryVertexValidation.cc:1961
TransientVertex::position
GlobalPoint position() const
Definition: TransientVertex.h:114
reco::TrackBase::pt
double pt() const
track transverse momentum
Definition: TrackBase.h:616
PV3DBase::z
T z() const
Definition: PV3DBase.h:64
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
reco::TrackBase::ptError
double ptError() const
error on Pt (set to 1000 TeV if charge==0 for safety)
Definition: TrackBase.h:758
Abs
T Abs(T a)
Definition: MathUtil.h:49
PrimaryVertexValidation::n_IP3DEtaResiduals
TH1F * n_IP3DEtaResiduals[nMaxBins_]
Definition: PrimaryVertexValidation.h:244
reco::TrackBase::qualityMask
int qualityMask() const
Definition: TrackBase.h:857
funct::tan
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
reco::Vertex::chi2
double chi2() const
chi-squares
Definition: Vertex.h:88
j
int j
Definition: DBlmapReader.cc:9
PrimaryVertexValidation::zUnbiasedVertex_
double zUnbiasedVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:187
reco::TrackBase::lambda
double lambda() const
Lambda angle.
Definition: TrackBase.h:580
PrimaryVertexValidation::a_reszPhiResiduals
TH1F * a_reszPhiResiduals[nMaxBins_]
Definition: PrimaryVertexValidation.h:226
reco::BeamSpot::BeamWidthX
double BeamWidthX() const
beam width X
Definition: BeamSpot.h:86
reco::TrackBase::numberOfValidHits
unsigned short numberOfValidHits() const
number of valid hits found
Definition: TrackBase.h:815
PrimaryVertexValidation::chi2normUnbiasedVertex_
float chi2normUnbiasedVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:188
PrimaryVertexValidation::n_d3DResidualsMap
TH1F * n_d3DResidualsMap[nMaxBins_][nMaxBins_]
Definition: PrimaryVertexValidation.h:260
TFileDirectory::make
T * make(const Args &...args) const
make new ROOT object
Definition: TFileDirectory.h:59
TrackFilterForPVFindingBase::select
virtual std::vector< reco::TransientTrack > select(const std::vector< reco::TransientTrack > &tracks) const =0
reco::Track::recHitsBegin
trackingRecHit_iterator recHitsBegin() const
Iterator to first hit on the track.
Definition: Track.h:104
PrimaryVertexValidation::theVertexCollectionToken
edm::EDGetTokenT< reco::VertexCollection > theVertexCollectionToken
Definition: PrimaryVertexValidation.h:118
PrimaryVertexValidation::a_d3DResidualsMap
TH1F * a_d3DResidualsMap[nMaxBins_][nMaxBins_]
Definition: PrimaryVertexValidation.h:256
TransientVertex::trackWeight
float trackWeight(const reco::TransientTrack &track) const
Definition: TransientVertex.cc:138
PrimaryVertexValidation::IPLsigFromMyVertex_
double IPLsigFromMyVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:204
LogTrace
#define LogTrace(id)
Definition: MessageLogger.h:502
reco::Vertex::ndof
double ndof() const
Definition: Vertex.h:95
reco::TrackBase::pz
double pz() const
z coordinate of momentum vector
Definition: TrackBase.h:634
reco::TrackBase::dz
double dz() const
dz parameter (= dsz/cos(lambda)). This is the track z0 w.r.t (0,0,0) only if the refPoint is close to...
Definition: TrackBase.h:604
reco::TrackBase::dzError
double dzError() const
error on dz
Definition: TrackBase.h:809
reco::TrackBase::vz
double vz() const
z coordinate of the reference point on track
Definition: TrackBase.h:664
edm::Event::eventAuxiliary
EventAuxiliary const & eventAuxiliary() const
Definition: Event.h:77
DetId
Definition: DetId.h:18
PrimaryVertexValidation::a_IP2DPhiResiduals
TH1F * a_IP2DPhiResiduals[nMaxBins_]
Definition: PrimaryVertexValidation.h:220
reco::TrackBase::qualityByName
static TrackQuality qualityByName(const std::string &name)
Definition: TrackBase.cc:125
reco::Vertex::xError
double xError() const
error on x
Definition: Vertex.h:107
TFileService::mkdir
TFileDirectory mkdir(const std::string &dir, const std::string &descr="")
create a new subdirectory
Definition: TFileService.h:69
reco::TransientTrack::track
const Track & track() const
Definition: TransientTrack.h:110
TransientVertex::hasTrackWeight
bool hasTrackWeight() const
Definition: TransientVertex.h:181
reco::TrackBase::hitPattern
const HitPattern & hitPattern() const
Access the hit pattern, indicating in which Tracker layers the track has hits.
Definition: TrackBase.h:445
PrimaryVertexValidation::isHit2D
bool isHit2D(const TrackingRecHit &hit) const
Definition: PrimaryVertexValidation.cc:798
math::XYZPoint
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
edm::EventSetup::get
const T & get() const
Definition: EventSetup.h:56
PrimaryVertexValidation::n_dxyPhiBiasResiduals
TH1F * n_dxyPhiBiasResiduals[nMaxBins_]
Definition: PrimaryVertexValidation.h:346
reco::HitPattern::pixelBarrelLayersWithMeasurement
int pixelBarrelLayersWithMeasurement() const
Definition: HitPattern.cc:567
reco::BeamSpot::BeamWidthY
double BeamWidthY() const
beam width Y
Definition: BeamSpot.h:88
PrimaryVertexValidation::theTrackClusterizer_
TrackClusterizerInZ * theTrackClusterizer_
Definition: PrimaryVertexValidation.h:94
b
double b
Definition: hdecay.h:120
PrimaryVertexValidation::xUnbiasedVertex_
double xUnbiasedVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:185
reco::TrackBase::quality
bool quality(const TrackQuality) const
Track quality.
Definition: TrackBase.h:505
create_public_pileup_plots.weights
list weights
Definition: create_public_pileup_plots.py:205
cscdqm::HistoName
std::string HistoName
Definition: CSCDQM_HistoDef.h:32
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
PrimaryVertexValidation::a_dxyEtaBiasResiduals
TH1F * a_dxyEtaBiasResiduals[nMaxBins_]
Definition: PrimaryVertexValidation.h:339
PrimaryVertexValidation::theBeamspotToken
edm::EDGetTokenT< reco::BeamSpot > theBeamspotToken
Definition: PrimaryVertexValidation.h:119
PrimaryVertexValidation::theTrackCollectionToken
edm::EDGetTokenT< reco::TrackCollection > theTrackCollectionToken
Definition: PrimaryVertexValidation.h:117
edm::EventAuxiliary::id
EventID const & id() const
Definition: EventAuxiliary.h:59
SiPixelCluster
Pixel cluster – collection of neighboring pixels above threshold.
Definition: SiPixelCluster.h:33
reco::TrackBase::vy
double vy() const
y coordinate of the reference point on track
Definition: TrackBase.h:658
PrimaryVertexValidation::dzFromMyVertex_
double dzFromMyVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:196
PrimaryVertexValidation::n_IP3DPhiResiduals
TH1F * n_IP3DPhiResiduals[nMaxBins_]
Definition: PrimaryVertexValidation.h:243
a
double a
Definition: hdecay.h:121
position
static int position[264][3]
Definition: ReadPGInfo.cc:509
reco::Track::found
unsigned short found() const
Number of valid hits on track.
Definition: Track.h:194
PrimaryVertexValidation::add
void add(std::map< std::string, TH1 * > &h, TH1 *hist)
Definition: PrimaryVertexValidation.cc:2306
PrimaryVertexValidation::fill
void fill(std::map< std::string, TH1 * > &h, std::string s, double x)
Definition: PrimaryVertexValidation.cc:2314
reco::Vertex::trackRef_iterator
std::vector< TrackBaseRef >::const_iterator trackRef_iterator
The iteratator for the vector&lt;TrackRef&gt;
Definition: Vertex.h:37
reco::BeamSpot::y0
double y0() const
y coordinate
Definition: BeamSpot.h:66
PrimaryVertexValidation::a_reszEtaResiduals
TH1F * a_reszEtaResiduals[nMaxBins_]
Definition: PrimaryVertexValidation.h:227
PrimaryVertexValidation::a_dxyBiasResidualsMap
TH1F * a_dxyBiasResidualsMap[nMaxBins_][nMaxBins_]
Definition: PrimaryVertexValidation.h:354
PrimaryVertexValidation::a_IP3DPhiResiduals
TH1F * a_IP3DPhiResiduals[nMaxBins_]
Definition: PrimaryVertexValidation.h:223
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:145
reco::TrackBase::charge
int charge() const
track electric charge
Definition: TrackBase.h:562
reco::BeamSpot::position
const Point & position() const
position
Definition: BeamSpot.h:62
PrimaryVertexValidation::n_dxyResidualsMap
TH1F * n_dxyResidualsMap[nMaxBins_][nMaxBins_]
Definition: PrimaryVertexValidation.h:258
hlt_dqm_clientPB-live_cfg.nbins
tuple nbins
Definition: hlt_dqm_clientPB-live_cfg.py:35
kappa
static const G4double kappa
Definition: UrbanMscModel93.cc:34
reco::Vertex::normalizedChi2
double normalizedChi2() const
chi-squared divided by n.d.o.f.
Definition: Vertex.h:97
PrimaryVertexValidation::dxyErrorFromMyVertex_
double dxyErrorFromMyVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:199
TrackingRecHit::geographicalId
DetId geographicalId() const
Definition: TrackingRecHit.h:116
dir
dbl *** dir
Definition: mlp_gen.cc:35
edm::false
volatile std::atomic< bool > shutdown_flag false
Definition: UnixSignalHandlers.cc:22
HLT_25ns10e33_v2_cff.clusters
tuple clusters
Definition: HLT_25ns10e33_v2_cff.py:8795
reco::Vertex::tracks_begin
trackRef_iterator tracks_begin() const
first iterator over tracks
Definition: Vertex.cc:39
PrimaryVertexValidation::n_reszPhiResiduals
TH1F * n_reszPhiResiduals[nMaxBins_]
Definition: PrimaryVertexValidation.h:246
PrimaryVertexValidation::sumOfWeightsUnbiasedVertex_
float sumOfWeightsUnbiasedVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:192
PrimaryVertexValidation::n_dzBiasResidualsMap
TH1F * n_dzBiasResidualsMap[nMaxBins_][nMaxBins_]
Definition: PrimaryVertexValidation.h:358
PrimaryVertexValidation::d3DErrorFromMyVertex_
double d3DErrorFromMyVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:201
reco::TrackBase::dxy
double dxy() const
dxy parameter. (This is the transverse impact parameter w.r.t. to (0,0,0) ONLY if refPoint is close t...
Definition: TrackBase.h:586
PrimaryVertexValidation::a_dzBiasResidualsMap
TH1F * a_dzBiasResidualsMap[nMaxBins_][nMaxBins_]
Definition: PrimaryVertexValidation.h:355
PV3DBase::x
T x() const
Definition: PV3DBase.h:62
PrimaryVertexValidation::a_dzEtaBiasResiduals
TH1F * a_dzEtaBiasResiduals[nMaxBins_]
Definition: PrimaryVertexValidation.h:342
TransientVertex::isValid
bool isValid() const
Definition: TransientVertex.h:134
PrimaryVertexValidation::fillTrendPlot
void fillTrendPlot(TH1F *trendPlot, TH1F *residualsPlot[100], TString fitPar_, TString var_)
Definition: PrimaryVertexValidation.cc:2004
PrimaryVertexValidation::a_dzResidualsMap
TH1F * a_dzResidualsMap[nMaxBins_][nMaxBins_]
Definition: PrimaryVertexValidation.h:255
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
reco::Vertex::yError
double yError() const
error on y
Definition: Vertex.h:109
PrimaryVertexValidation::chi2UnbiasedVertex_
float chi2UnbiasedVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:189
reco::Vertex::tracksSize
size_t tracksSize() const
number of tracks
Definition: Vertex.cc:34
point
*vegas h *****************************************************used in the default bin number in original ***version of VEGAS is ***a higher bin number might help to derive a more precise ***grade subtle point
Definition: invegas.h:5
reco::TrackBase::py
double py() const
y coordinate of momentum vector
Definition: TrackBase.h:628
trackingRecHit_iterator
TrackingRecHitCollection::base::const_iterator trackingRecHit_iterator
iterator over a vector of reference to TrackingRecHit in the same collection
Definition: TrackingRecHitFwd.h:19
reco::TrackBase::vx
double vx() const
x coordinate of the reference point on track
Definition: TrackBase.h:652
GlobalVector
Global3DVector GlobalVector
Definition: GlobalVector.h:10
reco::HitPattern::numberOfHits
int numberOfHits(HitCategory category) const
Definition: HitPattern.h:807
PrimaryVertexValidation::bookVertexHistograms
std::map< std::string, TH1 * > bookVertexHistograms(TFileDirectory dir)
Definition: PrimaryVertexValidation.cc:2140
SiPixelRecHit
Our base class.
Definition: SiPixelRecHit.h:23
reco::BeamSpot::x0
double x0() const
x coordinate
Definition: BeamSpot.h:64
kBPIX
const int kBPIX
Definition: AlignmentTrackSelector.cc:28
reco::Track::recHitsEnd
trackingRecHit_iterator recHitsEnd() const
Iterator to last hit on the track.
Definition: Track.h:109
PrimaryVertexValidation::IP3DsigFromMyVertex_
double IP3DsigFromMyVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:205
Generated for CMSSW Reference Manual by doxygen 1.8.5