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: PrimaryVertexValidation
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 
21 // user include files
22 #include "Alignment/OfflineValidation/plugins/PrimaryVertexValidation.h"
23 #include "FWCore/Framework/interface/Frameworkfwd.h"
24 #include "FWCore/Framework/interface/EDAnalyzer.h"
25 
26 #include "SimDataFormats/EncodedEventId/interface/EncodedEventId.h"
27 #include "FWCore/Framework/interface/Event.h"
28 #include <FWCore/Framework/interface/ESHandle.h>
29 #include "FWCore/Framework/interface/MakerMacros.h"
30 
31 #include "FWCore/ParameterSet/interface/ParameterSet.h"
32 
33 #include <SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h>
34 #include <SimDataFormats/TrackingHit/interface/PSimHit.h>
35 #include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
36 
37 #include "TH1F.h"
38 #include "TH2F.h"
39 #include "TF1.h"
40 #include "TVector3.h"
41 #include "TFile.h"
42 #include "TROOT.h"
43 #include "TChain.h"
44 #include "TNtuple.h"
45 #include <TMatrixD.h>
46 #include <TVectorD.h>
47 #include "FWCore/ServiceRegistry/interface/Service.h"
48 #include "CommonTools/UtilAlgos/interface/TFileService.h"
49 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
50 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
51 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
52 #include <Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h>
53 #include <DataFormats/GeometrySurface/interface/Surface.h>
54 #include <DataFormats/GeometrySurface/interface/GloballyPositioned.h>
55 #include <Geometry/CommonDetUnit/interface/GeomDet.h>
56 #include "MagneticField/Engine/interface/MagneticField.h"
57 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
58 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
59 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
60 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
61 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
62 
63 #include "DataFormats/TrackReco/interface/Track.h"
64 #include "DataFormats/TrackReco/interface/TrackFwd.h"
65 #include "DataFormats/VertexReco/interface/Vertex.h"
66 #include "DataFormats/VertexReco/interface/VertexFwd.h"
67 #include "RecoVertex/PrimaryVertexProducer/interface/GapClusterizerInZ.h"
68 #include "RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZ.h"
69 
70  const int kBPIX = PixelSubdetector::PixelBarrel;
71  const int kFPIX = PixelSubdetector::PixelEndcap;
72 
73 // Constructor
74 
75 PrimaryVertexValidation::PrimaryVertexValidation(const edm::ParameterSet& iConfig):
76  storeNtuple_(iConfig.getParameter<bool>("storeNtuple")),
77  lightNtupleSwitch_(iConfig.getParameter<bool>("isLightNtuple")),
78  useTracksFromRecoVtx_(iConfig.getParameter<bool>("useTracksFromRecoVtx")),
79  askFirstLayerHit_(iConfig.getParameter<bool>("askFirstLayerHit")),
80  ptOfProbe_(iConfig.getUntrackedParameter<double>("probePt",0.)),
81  etaOfProbe_(iConfig.getUntrackedParameter<double>("probeEta",2.4)),
82  nBins_(iConfig.getUntrackedParameter<int>("numberOfBins",24)),
83  debug_(iConfig.getParameter<bool>("Debug")),
84  TrackCollectionTag_(iConfig.getParameter<edm::InputTag>("TrackCollectionTag"))
85 {
86 
87  // now do what ever initialization is needed
88  // initialize phase space boundaries
89 
90  phipitch_ = (2*TMath::Pi())/nBins_;
91  etapitch_ = 5./nBins_;
92 
93  // old version
94  // theTrackClusterizer_ = new GapClusterizerInZ(iConfig.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<edm::ParameterSet>("TkGapClusParameters"));
95 
96  // select and configure the track filter
97  theTrackFilter_= new TrackFilterForPVFinding(iConfig.getParameter<edm::ParameterSet>("TkFilterParameters") );
98  // select and configure the track clusterizer
99  std::string clusteringAlgorithm=iConfig.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<std::string>("algorithm");
100  if (clusteringAlgorithm=="gap"){
101  theTrackClusterizer_ = new GapClusterizerInZ(iConfig.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<edm::ParameterSet>("TkGapClusParameters"));
102  }else if(clusteringAlgorithm=="DA"){
103  theTrackClusterizer_ = new DAClusterizerInZ(iConfig.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<edm::ParameterSet>("TkDAClusParameters"));
104  }else{
105  throw VertexException("PrimaryVertexProducerAlgorithm: unknown clustering algorithm: " + clusteringAlgorithm);
106  }
107 }
108 
109 // Destructor
110 PrimaryVertexValidation::~PrimaryVertexValidation()
111 {
112  // do anything here that needs to be done at desctruction time
113  // (e.g. close files, deallocate resources etc.)
114 }
115 
116 
117 //
118 // member functions
119 //
120 
121 // ------------ method called to for each event ------------
122 void
123 PrimaryVertexValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
124 {
125 
126  using namespace std;
127  using namespace IPTools;
128 
129  Nevt_++;
130 
131  //=======================================================
132  // Initialize Root-tuple variables
133  //=======================================================
134 
135  SetVarToZero();
136 
137  //=======================================================
138  // Retrieve the Magnetic Field information
139  //=======================================================
140 
141  edm::ESHandle<MagneticField> theMGField;
142  iSetup.get<IdealMagneticFieldRecord>().get( theMGField );
143 
144  //=======================================================
145  // Retrieve the Tracking Geometry information
146  //=======================================================
147 
148  edm::ESHandle<GlobalTrackingGeometry> theTrackingGeometry;
149  iSetup.get<GlobalTrackingGeometryRecord>().get( theTrackingGeometry );
150 
151  //=======================================================
152  // Retrieve the Track information
153  //=======================================================
154 
155  edm::Handle<reco::TrackCollection> trackCollectionHandle;
156  iEvent.getByLabel(TrackCollectionTag_, trackCollectionHandle);
157 
158  //=======================================================
159  // Retrieve offline vartex information (only for reco)
160  //=======================================================
161 
162  edm::Handle<reco::VertexCollection> vertices;
163  try {
164  iEvent.getByLabel("offlinePrimaryVertices", vertices);
165  } catch (...) {
166  if(debug_)
167  cout << "No offlinePrimaryVertices found!" << endl;
168  }
169  if ( vertices.isValid() ) {
170  xOfflineVertex_ = (*vertices)[0].x();
171  yOfflineVertex_ = (*vertices)[0].y();
172  zOfflineVertex_ = (*vertices)[0].z();
173  }
174 
175  unsigned int vertexCollectionSize = vertices.product()->size();
176  int nvvertex = 0;
177 
178  for (unsigned int i=0; i<vertexCollectionSize; i++) {
179  const reco::Vertex& vertex = vertices->at(i);
180  if (vertex.isValid()) nvvertex++;
181  }
182 
183  nOfflineVertices_ = nvvertex;
184 
185 
186  if ( vertices->size() && useTracksFromRecoVtx_ ) {
187 
188  double sumpt = 0;
189  size_t ntracks = 0;
190  double chi2ndf = 0.;
191  double chi2prob = 0.;
192 
193  if (!vertices->at(0).isFake()) {
194 
195  reco::Vertex pv = vertices->at(0);
196 
197  ntracks = pv.tracksSize();
198  chi2ndf = pv.normalizedChi2();
199  chi2prob = TMath::Prob(pv.chi2(),(int)pv.ndof());
200 
201  h_recoVtxNtracks_->Fill(ntracks);
202  h_recoVtxChi2ndf_->Fill(chi2ndf);
203  h_recoVtxChi2Prob_->Fill(chi2prob);
204 
205  for (reco::Vertex::trackRef_iterator itrk = pv.tracks_begin();itrk != pv.tracks_end(); ++itrk) {
206  double pt = (**itrk).pt();
207  sumpt += pt*pt;
208 
209  const math::XYZPoint myVertex(pv.position().x(),pv.position().y(),pv.position().z());
210 
211  double dxyRes = (**itrk).dxy(myVertex);
212  double dzRes = (**itrk).dz(myVertex);
213 
214  double dxy_err = (**itrk).dxyError();
215  double dz_err = (**itrk).dzError();
216 
217  float_t trackphi = ((**itrk).phi())*(180/TMath::Pi());
218  float_t tracketa = (**itrk).eta();
219 
220  for(Int_t i=0; i<nBins_; i++){
221 
222  float phiF = (-TMath::Pi()+i*phipitch_)*(180/TMath::Pi());
223  float phiL = (-TMath::Pi()+(i+1)*phipitch_)*(180/TMath::Pi());
224 
225  float etaF=-2.5+i*etapitch_;
226  float etaL=-2.5+(i+1)*etapitch_;
227 
228  if(trackphi >= phiF && trackphi < phiL ){
229 
230  a_dxyPhiBiasResiduals[i]->Fill(dxyRes*cmToum);
231  a_dzPhiBiasResiduals[i]->Fill(dzRes*cmToum);
232  n_dxyPhiBiasResiduals[i]->Fill((dxyRes)/dxy_err);
233  n_dzPhiBiasResiduals[i]->Fill((dzRes)/dz_err);
234 
235  for(Int_t j=0; j<nBins_; j++){
236 
237  float etaJ=-2.5+j*etapitch_;
238  float etaK=-2.5+(j+1)*etapitch_;
239 
240  if(tracketa >= etaJ && tracketa < etaK ){
241 
242  a_dxyBiasResidualsMap[i][j]->Fill(dxyRes*cmToum);
243  a_dzBiasResidualsMap[i][j]->Fill(dzRes*cmToum);
244 
245  n_dxyBiasResidualsMap[i][j]->Fill((dxyRes)/dxy_err);
246  n_dzBiasResidualsMap[i][j]->Fill((dzRes)/dz_err);
247 
248  }
249  }
250  }
251 
252  if(tracketa >= etaF && tracketa < etaL ){
253  a_dxyEtaBiasResiduals[i]->Fill(dxyRes*cmToum);
254  a_dzEtaBiasResiduals[i]->Fill(dzRes*cmToum);
255  n_dxyEtaBiasResiduals[i]->Fill((dxyRes)/dxy_err);
256  n_dzEtaBiasResiduals[i]->Fill((dzRes)/dz_err);
257  }
258  }
259  }
260 
261  h_recoVtxSumPt_->Fill(sumpt);
262 
263  }
264  }
265 
266  //=======================================================
267  // Retrieve Beamspot information
268  //=======================================================
269 
270  reco::BeamSpot beamSpot;
271  edm::Handle<reco::BeamSpot> beamSpotHandle;
272  iEvent.getByLabel("offlineBeamSpot", beamSpotHandle);
273 
274  if ( beamSpotHandle.isValid() )
275  {
276  beamSpot = *beamSpotHandle;
277  BSx0_ = beamSpot.x0();
278  BSy0_ = beamSpot.y0();
279  BSz0_ = beamSpot.z0();
280  Beamsigmaz_ = beamSpot.sigmaZ();
281  Beamdxdz_ = beamSpot.dxdz();
282  BeamWidthX_ = beamSpot.BeamWidthX();
283  BeamWidthY_ = beamSpot.BeamWidthY();
284  } else
285  {
286  if(debug_)
287  cout << "No BeamSpot found!" << endl;
288  }
289 
290  if(debug_)
291  std::cout<<"Beamspot x:"<<BSx0_<<" y:"<<BSy0_<<" z:"<<BSz0_<<std::endl;
292 
293  //double sigmaz = beamSpot.sigmaZ();
294  //double dxdz = beamSpot.dxdz();
295  //double BeamWidth = beamSpot.BeamWidth();
296 
297  //=======================================================
298  // Starts here ananlysis
299  //=======================================================
300 
301  RunNumber_=iEvent.eventAuxiliary().run();
302  LuminosityBlockNumber_=iEvent.eventAuxiliary().luminosityBlock();
303  EventNumber_=iEvent.eventAuxiliary().id().event();
304 
305  if(debug_)
306  std::cout<<"PrimaryVertexValidation::analyze() looping over "<<trackCollectionHandle->size()<< "tracks." <<std::endl;
307 
308  //======================================================
309  // Interface RECO tracks to vertex reconstruction
310  //======================================================
311 
312  std::vector<reco::TransientTrack> t_tks;
313  unsigned int k = 0;
314  for(reco::TrackCollection::const_iterator track = trackCollectionHandle->begin(); track!= trackCollectionHandle->end(); ++track, ++k){
315 
316  reco::TrackRef trackref(trackCollectionHandle,k);
317  reco::TransientTrack theTTRef = reco::TransientTrack(trackref, &*theMGField, theTrackingGeometry );
318  t_tks.push_back(theTTRef);
319 
320  }
321 
322  if(debug_) {cout << "PrimaryVertexValidation"
323  << "Found: " << t_tks.size() << " reconstructed tracks" << "\n";
324  }
325 
326  //======================================================
327  // select the tracks
328  //======================================================
329 
330  std::vector<reco::TransientTrack> seltks = theTrackFilter_->select(t_tks);
331 
332  //======================================================
333  // clusterize tracks in Z
334  //======================================================
335 
336  vector< vector<reco::TransientTrack> > clusters = theTrackClusterizer_->clusterize(seltks);
337 
338  if (debug_){
339  cout << " clustering returned "<< clusters.size() << " clusters from " << t_tks.size() << " selected tracks" <<endl;
340  }
341 
342  nClus_=clusters.size();
343 
344  //======================================================
345  // Starts loop on clusters
346  //======================================================
347 
348  for (vector< vector<reco::TransientTrack> >::const_iterator iclus = clusters.begin(); iclus != clusters.end(); iclus++) {
349 
350  nTracksPerClus_=0;
351 
352  unsigned int i = 0;
353  for(vector<reco::TransientTrack>::const_iterator theTrack = iclus->begin(); theTrack!= iclus->end(); ++theTrack, ++i)
354  {
355  if ( nTracks_ >= nMaxtracks_ ) {
356  std::cout << " PrimaryVertexValidation::analyze() : Warning - Number of tracks: " << nTracks_ << " , greater than " << nMaxtracks_ << std::endl;
357  continue;
358  }
359 
360  pt_[nTracks_] = theTrack->track().pt();
361  p_[nTracks_] = theTrack->track().p();
362  nhits_[nTracks_] = theTrack->track().numberOfValidHits();
363  eta_[nTracks_] = theTrack->track().eta();
364  theta_[nTracks_] = theTrack->track().theta();
365  phi_[nTracks_] = theTrack->track().phi();
366  chi2_[nTracks_] = theTrack->track().chi2();
367  chi2ndof_[nTracks_] = theTrack->track().normalizedChi2();
368  charge_[nTracks_] = theTrack->track().charge();
369  qoverp_[nTracks_] = theTrack->track().qoverp();
370  dz_[nTracks_] = theTrack->track().dz();
371  dxy_[nTracks_] = theTrack->track().dxy();
372 
373  reco::TrackBase::TrackQuality _trackQuality = reco::TrackBase::qualityByName("highPurity");
374  isHighPurity_[nTracks_] = theTrack->track().quality(_trackQuality);
375 
376  math::XYZPoint point(beamSpot.x0(),beamSpot.y0(), beamSpot.z0());
377  dxyBs_[nTracks_] = theTrack->track().dxy(point);
378  dzBs_[nTracks_] = theTrack->track().dz(point);
379 
380  xPCA_[nTracks_] = theTrack->track().vertex().x();
381  yPCA_[nTracks_] = theTrack->track().vertex().y();
382  zPCA_[nTracks_] = theTrack->track().vertex().z();
383 
384  //=======================================================
385  // Retrieve rechit information
386  //=======================================================
387 
388  int nRecHit1D =0;
389  int nRecHit2D =0;
390  int nhitinTIB =0;
391  int nhitinTOB =0;
392  int nhitinTID =0;
393  int nhitinTEC =0;
394  int nhitinBPIX=0;
395  int nhitinFPIX=0;
396 
397  for (trackingRecHit_iterator iHit = theTrack->recHitsBegin(); iHit != theTrack->recHitsEnd(); ++iHit) {
398  if((*iHit)->isValid()) {
399 
400  if (this->isHit2D(**iHit)) {++nRecHit2D;}
401  else {++nRecHit1D; }
402 
403  int type =(*iHit)->geographicalId().subdetId();
404 
405  if(type==int(StripSubdetector::TIB)){++nhitinTIB;}
406  if(type==int(StripSubdetector::TOB)){++nhitinTOB;}
407  if(type==int(StripSubdetector::TID)){++nhitinTID;}
408  if(type==int(StripSubdetector::TEC)){++nhitinTEC;}
409  if(type==int( kBPIX)){++nhitinBPIX;}
410  if(type==int( kFPIX)){++nhitinFPIX;}
411  }
412  }
413 
414  nhits1D_[nTracks_] = nRecHit1D;
415  nhits2D_[nTracks_] = nRecHit2D;
416  nhitsBPIX_[nTracks_] = nhitinBPIX;
417  nhitsFPIX_[nTracks_] = nhitinFPIX;
418  nhitsTIB_[nTracks_] = nhitinTIB;
419  nhitsTID_[nTracks_] = nhitinTID;
420  nhitsTOB_[nTracks_] = nhitinTOB;
421  nhitsTEC_[nTracks_] = nhitinTEC;
422 
423  //=======================================================
424  // Good tracks for vertexing selection
425  //=======================================================
426 
427  bool pass = true;
428  if(askFirstLayerHit_) pass = this->hasFirstLayerPixelHits((*theTrack));
429  if (pass){
430  isGoodTrack_[nTracks_]=1;
431  }
432 
433  //=======================================================
434  // Fit unbiased vertex
435  //=======================================================
436 
437  vector<reco::TransientTrack> theFinalTracks;
438  theFinalTracks.clear();
439 
440  for(vector<reco::TransientTrack>::const_iterator tk = iclus->begin(); tk!= iclus->end(); ++tk){
441 
442  pass = this->hasFirstLayerPixelHits((*tk));
443  if (pass){
444  if( tk == theTrack ) continue;
445  else {
446  theFinalTracks.push_back((*tk));
447  }
448  }
449  }
450 
451  if(theFinalTracks.size() > 2){
452 
453  if(debug_)
454  std::cout <<"PrimaryVertexValidation::analyze() :Transient Track Collection size: "<<theFinalTracks.size()<<std::endl;
455 
456  try{
457 
458  VertexFitter<5>* fitter = new AdaptiveVertexFitter;
459  TransientVertex theFittedVertex = fitter->vertex(theFinalTracks);
460 
461  if(theFittedVertex.isValid ()){
462 
463  if(theFittedVertex.hasTrackWeight()){
464  for(size_t rtracks= 0; rtracks < theFinalTracks.size(); rtracks++){
465  sumOfWeightsUnbiasedVertex_[nTracks_] += theFittedVertex.trackWeight(theFinalTracks[rtracks]);
466  }
467  }
468 
469  const math::XYZPoint theRecoVertex(xOfflineVertex_,yOfflineVertex_,zOfflineVertex_);
470 
471  const math::XYZPoint myVertex(theFittedVertex.position().x(),theFittedVertex.position().y(),theFittedVertex.position().z());
472  hasRecVertex_[nTracks_] = 1;
473  xUnbiasedVertex_[nTracks_] = theFittedVertex.position().x();
474  yUnbiasedVertex_[nTracks_] = theFittedVertex.position().y();
475  zUnbiasedVertex_[nTracks_] = theFittedVertex.position().z();
476 
477  chi2normUnbiasedVertex_[nTracks_] = theFittedVertex.normalisedChiSquared();
478  chi2UnbiasedVertex_[nTracks_] = theFittedVertex.totalChiSquared();
479  DOFUnbiasedVertex_[nTracks_] = theFittedVertex.degreesOfFreedom();
480  tracksUsedForVertexing_[nTracks_] = theFinalTracks.size();
481 
482  h_fitVtxNtracks_->Fill(theFinalTracks.size());
483  h_fitVtxChi2_->Fill(theFittedVertex.totalChiSquared());
484  h_fitVtxNdof_->Fill(theFittedVertex.degreesOfFreedom());
485  h_fitVtxChi2ndf_->Fill(theFittedVertex.normalisedChiSquared());
486  h_fitVtxChi2Prob_->Fill(TMath::Prob(theFittedVertex.totalChiSquared(),(int)theFittedVertex.degreesOfFreedom()));
487 
488  dszFromMyVertex_[nTracks_] = theTrack->track().dsz(myVertex);
489  dxyFromMyVertex_[nTracks_] = theTrack->track().dxy(myVertex);
490  dzFromMyVertex_[nTracks_] = theTrack->track().dz(myVertex);
491 
492  double dz_err = hypot(theTrack->track().dzError(),theFittedVertex.positionError().czz());
493 
494  // PV2D
495 
496  // std::pair<bool,Measurement1D> ip2dpv = absoluteTransverseImpactParameter(*theTrack,theFittedVertex);
497  // double ip2d_corr = ip2dpv.second.value();
498  // double ip2d_err = ip2dpv.second.error();
499 
500  std::pair<bool,Measurement1D> s_ip2dpv =
501  signedTransverseImpactParameter(*theTrack,GlobalVector(theTrack->track().px(),
502  theTrack->track().py(),
503  theTrack->track().pz()),
504  theFittedVertex);
505 
506  // double s_ip2dpv_corr = s_ip2dpv.second.value();
507  double s_ip2dpv_err = s_ip2dpv.second.error();
508 
509  // PV3D
510 
511  // std::pair<bool,Measurement1D> ip3dpv = absoluteImpactParameter3D(*theTrack,theFittedVertex);
512  // double ip3d_corr = ip3dpv.second.value();
513  // double ip3d_err = ip3dpv.second.error();
514 
515  // std::pair<bool,Measurement1D> s_ip3dpv =
516  // signedImpactParameter3D(*theTrack,GlobalVector(theTrack->track().px(),
517  // theTrack->track().py(),
518  // theTrack->track().pz()),
519  // theFittedVertex);
520 
521  // double s_ip3dpv_corr = s_ip3dpv.second.value();
522  // double s_ip3dpv_err = s_ip3dpv.second.error();
523 
524  dxyErrorFromMyVertex_[nTracks_] = s_ip2dpv_err;
525  dzErrorFromMyVertex_[nTracks_] = dz_err;
526 
527  IPTsigFromMyVertex_[nTracks_] = (theTrack->track().dxy(myVertex))/s_ip2dpv_err;
528  IPLsigFromMyVertex_[nTracks_] = (theTrack->track().dz(myVertex))/dz_err;
529 
530  // fill directly the histograms of residuals
531 
532  float_t trackphi = (theTrack->track().phi())*(180/TMath::Pi());
533  float_t tracketa = theTrack->track().eta();
534  float_t trackpt = theTrack->track().pt();
535 
536  // checks on the probe track quality
537  if(trackpt >= ptOfProbe_ && fabs(tracketa)<= etaOfProbe_){
538 
539  h_probePt_->Fill(theTrack->track().pt());
540  h_probeEta_->Fill(theTrack->track().eta());
541  h_probePhi_->Fill(theTrack->track().phi());
542  h_probeChi2_->Fill(theTrack->track().chi2());
543  h_probeNormChi2_->Fill(theTrack->track().normalizedChi2());
544  h_probeCharge_->Fill(theTrack->track().charge());
545  h_probeQoverP_->Fill(theTrack->track().qoverp());
546  h_probedz_->Fill(theTrack->track().dz(theRecoVertex));
547  h_probedxy_->Fill((theTrack->track().dxy(theRecoVertex)));
548 
549  h_probeHits_->Fill(theTrack->track().numberOfValidHits());
550  h_probeHits1D_->Fill(nRecHit1D);
551  h_probeHits2D_->Fill(nRecHit2D);
552  h_probeHitsInTIB_->Fill(nhitinBPIX);
553  h_probeHitsInTOB_->Fill(nhitinFPIX);
554  h_probeHitsInTID_->Fill(nhitinTIB);
555  h_probeHitsInTEC_->Fill(nhitinTID);
556  h_probeHitsInBPIX_->Fill(nhitinTOB);
557  h_probeHitsInFPIX_->Fill(nhitinTEC);
558 
559  for(Int_t i=0; i<nBins_; i++){
560 
561  float phiF = (-TMath::Pi()+i*phipitch_)*(180/TMath::Pi());
562  float phiL = (-TMath::Pi()+(i+1)*phipitch_)*(180/TMath::Pi());
563 
564  float etaF=-2.5+i*etapitch_;
565  float etaL=-2.5+(i+1)*etapitch_;
566 
567  if(trackphi >= phiF && trackphi < phiL ){
568  a_dxyPhiResiduals[i]->Fill(theTrack->track().dxy(myVertex)*cmToum);
569  a_dzPhiResiduals[i]->Fill(theTrack->track().dz(myVertex)*cmToum);
570  n_dxyPhiResiduals[i]->Fill((theTrack->track().dxy(myVertex))/s_ip2dpv_err);
571  n_dzPhiResiduals[i]->Fill((theTrack->track().dz(myVertex))/dz_err);
572 
573  for(Int_t j=0; j<nBins_; j++){
574 
575  float etaJ=-2.5+j*etapitch_;
576  float etaK=-2.5+(j+1)*etapitch_;
577 
578  if(tracketa >= etaJ && tracketa < etaK ){
579 
580  a_dxyResidualsMap[i][j]->Fill(theTrack->track().dxy(myVertex)*cmToum);
581  a_dzResidualsMap[i][j]->Fill(theTrack->track().dz(myVertex)*cmToum);
582 
583  n_dxyResidualsMap[i][j]->Fill((theTrack->track().dxy(myVertex))/s_ip2dpv_err);
584  n_dzResidualsMap[i][j]->Fill((theTrack->track().dz(myVertex))/dz_err);
585 
586  }
587  }
588  }
589 
590  if(tracketa >= etaF && tracketa < etaL ){
591  a_dxyEtaResiduals[i]->Fill(theTrack->track().dxy(myVertex)*cmToum);
592  a_dzEtaResiduals[i]->Fill(theTrack->track().dz(myVertex)*cmToum);
593  n_dxyEtaResiduals[i]->Fill((theTrack->track().dxy(myVertex))/s_ip2dpv_err);
594  n_dzEtaResiduals[i]->Fill((theTrack->track().dz(myVertex))/dz_err);
595  }
596  }
597  }
598 
599  if(debug_){
600  std::cout<<"PrimaryVertexValidation::analyze() : myVertex.x()= "<<myVertex.x()<<" myVertex.y()= "<<myVertex.y()<<" theFittedVertex.z()= "<<myVertex.z()<<std::endl;
601  std::cout<<"PrimaryVertexValidation::analyze() : theTrack->track().dz(myVertex)= "<<theTrack->track().dz(myVertex)<<std::endl;
602  std::cout<<"PrimaryVertexValidation::analyze() : zPCA -myVertex.z() = "<<(theTrack->track().vertex().z() -myVertex.z() )<<std::endl;
603  }// ends if debug_
604  } // ends if the fitted vertex is Valid
605 
606  delete fitter;
607 
608  } catch ( cms::Exception& er ) {
609  LogTrace("PrimaryVertexValidation::analyze RECO")<<"caught std::exception "<<er.what()<<std::endl;
610  }
611 
612  } //ends if theFinalTracks.size() > 2
613 
614  else {
615  if(debug_)
616  std::cout << "PrimaryVertexValidation::analyze() :Not enough tracks to make a vertex. Returns no vertex info" << std::endl;
617  }
618 
619  ++nTracks_;
620  ++nTracksPerClus_;
621 
622  if(debug_)
623  cout<< "Track "<<i<<" : pT = "<<theTrack->track().pt()<<endl;
624 
625  }// for loop on tracks
626 
627  } // for loop on track clusters
628 
629 
630  // Fill the TTree if needed
631 
632  if(storeNtuple_){
633  rootTree_->Fill();
634  }
635 
636 
637 }
638 
639 // ------------ method called to discriminate 1D from 2D hits ------------
640 bool PrimaryVertexValidation::isHit2D(const TrackingRecHit &hit) const
641 {
642  if (hit.dimension() < 2) {
643  return false; // some (muon...) stuff really has RecHit1D
644  } else {
645  const DetId detId(hit.geographicalId());
646  if (detId.det() == DetId::Tracker) {
647  if (detId.subdetId() == kBPIX || detId.subdetId() == kFPIX) {
648  return true; // pixel is always 2D
649  } else { // should be SiStrip now
650  if (dynamic_cast<const SiStripRecHit2D*>(&hit)) return false; // normal hit
651  else if (dynamic_cast<const SiStripMatchedRecHit2D*>(&hit)) return true; // matched is 2D
652  else if (dynamic_cast<const ProjectedSiStripRecHit2D*>(&hit)) return false; // crazy hit...
653  else {
654  edm::LogError("UnkownType") << "@SUB=AlignmentTrackSelector::isHit2D"
655  << "Tracker hit not in pixel and neither SiStripRecHit2D nor "
656  << "SiStripMatchedRecHit2D nor ProjectedSiStripRecHit2D.";
657  return false;
658  }
659  }
660  } else { // not tracker??
661  edm::LogWarning("DetectorMismatch") << "@SUB=AlignmentTrackSelector::isHit2D"
662  << "Hit not in tracker with 'official' dimension >=2.";
663  return true; // dimension() >= 2 so accept that...
664  }
665  }
666  // never reached...
667 }
668 
669 // ------------ method to check the presence of pixel hits ------------
670 bool PrimaryVertexValidation::hasFirstLayerPixelHits(const reco::TransientTrack track)
671 {
672  using namespace reco;
673  const HitPattern &p = track.hitPattern();
674  for (int i = 0; i < p.numberOfHits(HitPattern::TRACK_HITS); i++) {
675  uint32_t pattern = p.getHitPattern(HitPattern::TRACK_HITS, i);
676  if (p.pixelBarrelHitFilter(pattern) || p.pixelEndcapHitFilter(pattern) ) {
677  if (p.getLayer(pattern) == 1) {
678  if (p.validHitFilter(pattern)) {
679  return true;
680  }
681  }
682  }
683  }
684  return false;
685 }
686 
687 // ------------ method called once each job before begining the event loop ------------
688 void PrimaryVertexValidation::beginJob()
689 {
690  edm::LogInfo("beginJob") << "Begin Job" << std::endl;
691  // Define TTree for output
692  Nevt_ = 0;
693 
694  // rootFile_ = new TFile(filename_.c_str(),"recreate");
695  edm::Service<TFileService> fs;
696  rootTree_ = fs->make<TTree>("tree","PV Validation tree");
697 
698  // Track Paramters
699 
700  if(lightNtupleSwitch_){
701 
702  rootTree_->Branch("EventNumber",&EventNumber_,"EventNumber/i");
703  rootTree_->Branch("RunNumber",&RunNumber_,"RunNumber/i");
704  rootTree_->Branch("LuminosityBlockNumber",&LuminosityBlockNumber_,"LuminosityBlockNumber/i");
705  rootTree_->Branch("nOfflineVertices",&nOfflineVertices_,"nOfflineVertices/I");
706  rootTree_->Branch("nTracks",&nTracks_,"nTracks/I");
707  rootTree_->Branch("phi",&phi_,"phi[nTracks]/D");
708  rootTree_->Branch("eta",&eta_,"eta[nTracks]/D");
709  rootTree_->Branch("pt",&pt_,"pt[nTracks]/D");
710  rootTree_->Branch("dxyFromMyVertex",&dxyFromMyVertex_,"dxyFromMyVertex[nTracks]/D");
711  rootTree_->Branch("dzFromMyVertex",&dzFromMyVertex_,"dzFromMyVertex[nTracks]/D");
712  rootTree_->Branch("IPTsigFromMyVertex",&IPTsigFromMyVertex_,"IPTsigFromMyVertex_[nTracks]/D");
713  rootTree_->Branch("IPLsigFromMyVertex",&IPLsigFromMyVertex_,"IPLsigFromMyVertex_[nTracks]/D");
714  rootTree_->Branch("hasRecVertex",&hasRecVertex_,"hasRecVertex[nTracks]/I");
715  rootTree_->Branch("isGoodTrack",&isGoodTrack_,"isGoodTrack[nTracks]/I");
716  rootTree_->Branch("isHighPurity",&isHighPurity_,"isHighPurity_[nTracks]/I");
717 
718  } else {
719 
720  rootTree_->Branch("nTracks",&nTracks_,"nTracks/I");
721  rootTree_->Branch("nTracksPerClus",&nTracksPerClus_,"nTracksPerClus/I");
722  rootTree_->Branch("nClus",&nClus_,"nClus/I");
723  rootTree_->Branch("xOfflineVertex",&xOfflineVertex_,"xOfflineVertex/D");
724  rootTree_->Branch("yOfflineVertex",&yOfflineVertex_,"yOfflineVertex/D");
725  rootTree_->Branch("zOfflineVertex",&zOfflineVertex_,"zOfflineVertex/D");
726  rootTree_->Branch("BSx0",&BSx0_,"BSx0/D");
727  rootTree_->Branch("BSy0",&BSy0_,"BSy0/D");
728  rootTree_->Branch("BSz0",&BSz0_,"BSz0/D");
729  rootTree_->Branch("Beamsigmaz",&Beamsigmaz_,"Beamsigmaz/D");
730  rootTree_->Branch("Beamdxdz",&Beamdxdz_,"Beamdxdz/D");
731  rootTree_->Branch("BeamWidthX",&BeamWidthX_,"BeamWidthX/D");
732  rootTree_->Branch("BeamWidthY",&BeamWidthY_,"BeamWidthY/D");
733  rootTree_->Branch("pt",&pt_,"pt[nTracks]/D");
734  rootTree_->Branch("p",&p_,"p[nTracks]/D");
735  rootTree_->Branch("nhits",&nhits_,"nhits[nTracks]/I");
736  rootTree_->Branch("nhits1D",&nhits1D_,"nhits1D[nTracks]/I");
737  rootTree_->Branch("nhits2D",&nhits2D_,"nhits2D[nTracks]/I");
738  rootTree_->Branch("nhitsBPIX",&nhitsBPIX_,"nhitsBPIX[nTracks]/I");
739  rootTree_->Branch("nhitsFPIX",&nhitsFPIX_,"nhitsFPIX[nTracks]/I");
740  rootTree_->Branch("nhitsTIB",&nhitsTIB_,"nhitsTIB[nTracks]/I");
741  rootTree_->Branch("nhitsTID",&nhitsTID_,"nhitsTID[nTracks]/I");
742  rootTree_->Branch("nhitsTOB",&nhitsTOB_,"nhitsTOB[nTracks]/I");
743  rootTree_->Branch("nhitsTEC",&nhitsTEC_,"nhitsTEC[nTracks]/I");
744  rootTree_->Branch("eta",&eta_,"eta[nTracks]/D");
745  rootTree_->Branch("theta",&theta_,"theta[nTracks]/D");
746  rootTree_->Branch("phi",&phi_,"phi[nTracks]/D");
747  rootTree_->Branch("chi2",&chi2_,"chi2[nTracks]/D");
748  rootTree_->Branch("chi2ndof",&chi2ndof_,"chi2ndof[nTracks]/D");
749  rootTree_->Branch("charge",&charge_,"charge[nTracks]/I");
750  rootTree_->Branch("qoverp",&qoverp_,"qoverp[nTracks]/D");
751  rootTree_->Branch("dz",&dz_,"dz[nTracks]/D");
752  rootTree_->Branch("dxy",&dxy_,"dxy[nTracks]/D");
753  rootTree_->Branch("dzBs",&dzBs_,"dzBs[nTracks]/D");
754  rootTree_->Branch("dxyBs",&dxyBs_,"dxyBs[nTracks]/D");
755  rootTree_->Branch("xPCA",&xPCA_,"xPCA[nTracks]/D");
756  rootTree_->Branch("yPCA",&yPCA_,"yPCA[nTracks]/D");
757  rootTree_->Branch("zPCA",&zPCA_,"zPCA[nTracks]/D");
758  rootTree_->Branch("xUnbiasedVertex",&xUnbiasedVertex_,"xUnbiasedVertex[nTracks]/D");
759  rootTree_->Branch("yUnbiasedVertex",&yUnbiasedVertex_,"yUnbiasedVertex[nTracks]/D");
760  rootTree_->Branch("zUnbiasedVertex",&zUnbiasedVertex_,"zUnbiasedVertex[nTracks]/D");
761  rootTree_->Branch("chi2normUnbiasedVertex",&chi2normUnbiasedVertex_,"chi2normUnbiasedVertex[nTracks]/F");
762  rootTree_->Branch("chi2UnbiasedVertex",&chi2UnbiasedVertex_,"chi2UnbiasedVertex[nTracks]/F");
763  rootTree_->Branch("DOFUnbiasedVertex",&DOFUnbiasedVertex_," DOFUnbiasedVertex[nTracks]/F");
764  rootTree_->Branch("sumOfWeightsUnbiasedVertex",&sumOfWeightsUnbiasedVertex_,"sumOfWeightsUnbiasedVertex[nTracks]/F");
765  rootTree_->Branch("tracksUsedForVertexing",&tracksUsedForVertexing_,"tracksUsedForVertexing[nTracks]/I");
766  rootTree_->Branch("dxyFromMyVertex",&dxyFromMyVertex_,"dxyFromMyVertex[nTracks]/D");
767  rootTree_->Branch("dzFromMyVertex",&dzFromMyVertex_,"dzFromMyVertex[nTracks]/D");
768  rootTree_->Branch("dszFromMyVertex",&dszFromMyVertex_,"dszFromMyVertex[nTracks]/D");
769  rootTree_->Branch("dxyErrorFromMyVertex",&dxyErrorFromMyVertex_,"dxyErrorFromMyVertex_[nTracks]/D");
770  rootTree_->Branch("dzErrorFromMyVertex",&dzErrorFromMyVertex_,"dzErrorFromMyVertex_[nTracks]/D");
771  rootTree_->Branch("IPTsigFromMyVertex",&IPTsigFromMyVertex_,"IPTsigFromMyVertex_[nTracks]/D");
772  rootTree_->Branch("IPLsigFromMyVertex",&IPLsigFromMyVertex_,"IPLsigFromMyVertex_[nTracks]/D");
773  rootTree_->Branch("hasRecVertex",&hasRecVertex_,"hasRecVertex[nTracks]/I");
774  rootTree_->Branch("isGoodTrack",&isGoodTrack_,"isGoodTrack[nTracks]/I");
775  }
776 
777  // probe track histograms
778 
779  TFileDirectory ProbeFeatures = fs->mkdir("ProbeTrackFeatures");
780 
781  h_probePt_ = ProbeFeatures.make<TH1F>("h_probePt","p_{T} of probe track;track p_{T} (GeV); tracks",100,0.,50.);
782  h_probeEta_ = ProbeFeatures.make<TH1F>("h_probeEta","#eta of probe track;track #eta; tracks",54,-2.7,2.7);
783  h_probePhi_ = ProbeFeatures.make<TH1F>("h_probePhi","#phi of probe track;track #phi [rad]; tracks",100,-3.15,3.15);
784  h_probeChi2_ = ProbeFeatures.make<TH1F>("h_probeChi2","#chi^{2} of probe track;track #chi^{2}; tracks",100,0.,100.);
785  h_probeNormChi2_ = ProbeFeatures.make<TH1F>("h_probeNormChi2"," normalized #chi^{2} of probe track;track #chi^{2}/ndof; tracks",100,0.,10.);
786  h_probeCharge_ = ProbeFeatures.make<TH1F>("h_probeCharge","charge of profe track;track charge Q;tracks",3,-1.5,1.5);
787  h_probeQoverP_ = ProbeFeatures.make<TH1F>("h_probeQoverP","q/p of probe track; track Q/p (GeV^{-1})",200,-1.,1.);
788  h_probedz_ = ProbeFeatures.make<TH1F>("h_probedz","d_{z} of probe track;track d_{z} (cm);tracks",100,-5.,5.);
789  h_probedxy_ = ProbeFeatures.make<TH1F>("h_probedxy","d_{xy} of probe track;track d_{xy} (#mum);tracks",200,-1.,1.);
790 
791  h_probeHits_ = ProbeFeatures.make<TH1F>("h_probeNRechits" ,"N_{hits} ;N_{hits} ;tracks",40,-0.5,39.5);
792  h_probeHits1D_ = ProbeFeatures.make<TH1F>("h_probeNRechits1D" ,"N_{hits} 1D ;N_{hits} 1D ;tracks",40,-0.5,39.5);
793  h_probeHits2D_ = ProbeFeatures.make<TH1F>("h_probeNRechits2D" ,"N_{hits} 2D ;N_{hits} 2D ;tracks",40,-0.5,39.5);
794  h_probeHitsInTIB_ = ProbeFeatures.make<TH1F>("h_probeNRechitsTIB" ,"N_{hits} TIB ;N_{hits} TIB;tracks",40,-0.5,39.5);
795  h_probeHitsInTOB_ = ProbeFeatures.make<TH1F>("h_probeNRechitsTOB" ,"N_{hits} TOB ;N_{hits} TOB;tracks",40,-0.5,39.5);
796  h_probeHitsInTID_ = ProbeFeatures.make<TH1F>("h_probeNRechitsTID" ,"N_{hits} TID ;N_{hits} TID;tracks",40,-0.5,39.5);
797  h_probeHitsInTEC_ = ProbeFeatures.make<TH1F>("h_probeNRechitsTEC" ,"N_{hits} TEC ;N_{hits} TEC;tracks",40,-0.5,39.5);
798  h_probeHitsInBPIX_ = ProbeFeatures.make<TH1F>("h_probeNRechitsBPIX","N_{hits} BPIX;N_{hits} BPIX;tracks",40,-0.5,39.5);
799  h_probeHitsInFPIX_ = ProbeFeatures.make<TH1F>("h_probeNRechitsFPIX","N_{hits} FPIX;N_{hits} FPIX;tracks",40,-0.5,39.5);
800 
801  TFileDirectory RefitVertexFeatures = fs->mkdir("RefitVertexFeatures");
802  h_fitVtxNtracks_ = RefitVertexFeatures.make<TH1F>("h_fitVtxNtracks" ,"N^{vtx}_{trks};N^{vtx}_{trks};vertices" ,100,-0.5,99.5);
803  h_fitVtxNdof_ = RefitVertexFeatures.make<TH1F>("h_fitVtxNdof" ,"N^{vtx}_{DOF};N^{vtx}_{DOF};vertices" ,100,-0.5,99.5);
804  h_fitVtxChi2_ = RefitVertexFeatures.make<TH1F>("h_fitVtxChi2" ,"#chi^{2} vtx;#chi^{2} vtx;vertices" ,100,-0.5,99.5);
805  h_fitVtxChi2ndf_ = RefitVertexFeatures.make<TH1F>("h_fitVtxChi2ndf" ,"#chi^{2}/ndf vtx;#chi^{2}/ndf vtx;vertices" ,100,-0.5,9.5);
806  h_fitVtxChi2Prob_ = RefitVertexFeatures.make<TH1F>("h_fitVtxChi2Prob" ,"Prob(#chi^{2},ndf);Prob(#chi^{2},ndf);vertices",40,0.,1.);
807 
808  if(useTracksFromRecoVtx_) {
809 
810  TFileDirectory RecoVertexFeatures = fs->mkdir("RecoVertexFeatures");
811  h_recoVtxNtracks_ = RecoVertexFeatures.make<TH1F>("h_recoVtxNtracks" ,"N^{vtx}_{trks};N^{vtx}_{trks};vertices" ,100,-0.5,99.5);
812  h_recoVtxChi2ndf_ = RecoVertexFeatures.make<TH1F>("h_recoVtxChi2ndf" ,"#chi^{2}/ndf vtx;#chi^{2}/ndf vtx;vertices" ,10,-0.5,9.5);
813  h_recoVtxChi2Prob_ = RecoVertexFeatures.make<TH1F>("h_recoVtxChi2Prob" ,"Prob(#chi^{2},ndf);Prob(#chi^{2},ndf);vertices",40,0.,1.);
814  h_recoVtxSumPt_ = RecoVertexFeatures.make<TH1F>("h_recoVtxSumPt" ,"Sum(p^{trks}_{T});Sum(p^{trks}_{T});vertices" ,100,0.,200.);
815 
816  }
817 
818  // initialize the residuals histograms
819 
820  float dxymax_phi = 2000;
821  float dzmax_phi = 2000;
822  float dxymax_eta = 3000;
823  float dzmax_eta = 3000;
824 
825  const Int_t mybins_ = 500;
826 
828  //
829  // Usual plots from refitting the vertex
830  // The vertex is refit without the probe track
831  //
833 
834  TFileDirectory AbsTransPhiRes = fs->mkdir("Abs_Transv_Phi_Residuals");
835  TFileDirectory AbsTransEtaRes = fs->mkdir("Abs_Transv_Eta_Residuals");
836 
837  TFileDirectory AbsLongPhiRes = fs->mkdir("Abs_Long_Phi_Residuals");
838  TFileDirectory AbsLongEtaRes = fs->mkdir("Abs_Long_Eta_Residuals");
839 
840  TFileDirectory NormTransPhiRes = fs->mkdir("Norm_Transv_Phi_Residuals");
841  TFileDirectory NormTransEtaRes = fs->mkdir("Norm_Transv_Eta_Residuals");
842 
843  TFileDirectory NormLongPhiRes = fs->mkdir("Norm_Long_Phi_Residuals");
844  TFileDirectory NormLongEtaRes = fs->mkdir("Norm_Long_Eta_Residuals");
845 
846  TFileDirectory AbsDoubleDiffRes = fs->mkdir("Abs_DoubleDiffResiduals");
847  TFileDirectory NormDoubleDiffRes = fs->mkdir("Norm_DoubleDiffResiduals");
848 
849  for ( int i=0; i<nBins_; ++i ) {
850 
851  float phiF = (-TMath::Pi()+i*phipitch_)*(180/TMath::Pi());
852  float phiL = (-TMath::Pi()+(i+1)*phipitch_)*(180/TMath::Pi());
853 
854  float etaF=-2.5+i*etapitch_;
855  float etaL=-2.5+(i+1)*etapitch_;
856 
857  a_dxyPhiResiduals[i] = AbsTransPhiRes.make<TH1F>(Form("histo_dxy_phi_plot%i",i),Form("%.2f#circ<#varphi^{probe}_{tk}<%.2f#circ;d_{xy} [#mum];tracks",phiF,phiL),mybins_,-dxymax_phi,dxymax_phi);
858  a_dxyEtaResiduals[i] = AbsTransEtaRes.make<TH1F>(Form("histo_dxy_eta_plot%i",i),Form("%.2f<#eta^{probe}_{tk}<%.2f;d_{xy} [#mum];tracks",etaF,etaL),mybins_,-dxymax_eta,dxymax_eta);
859 
860  a_dzPhiResiduals[i] = AbsLongPhiRes.make<TH1F>(Form("histo_dz_phi_plot%i",i),Form("%.2f#circ<#varphi^{probe}_{tk}<%.2f #circ;d_{z} [#mum];tracks",phiF,phiL),mybins_,-dzmax_phi,dzmax_phi);
861  a_dzEtaResiduals[i] = AbsLongEtaRes.make<TH1F>(Form("histo_dz_eta_plot%i",i),Form("%.2f<#eta^{probe}_{tk}<%.2f;d_{z} [#mum];tracks",etaF,etaL),mybins_,-dzmax_eta,dzmax_eta);
862 
863  n_dxyPhiResiduals[i] = NormTransPhiRes.make<TH1F>(Form("histo_norm_dxy_phi_plot%i",i),Form("%.2f#circ<#varphi^{probe}_{tk}<%.2f#circ;d_{xy}/#sigma_{d_{xy}} [#mum];tracks",phiF,phiL),mybins_,-dxymax_phi/100.,dxymax_phi/100.);
864  n_dxyEtaResiduals[i] = NormTransEtaRes.make<TH1F>(Form("histo_norm_dxy_eta_plot%i",i),Form("%.2f<#eta^{probe}_{tk}<%.2f;d_{xy}/#sigma_{d_{xy}} [#mum];tracks",etaF,etaL),mybins_,-dxymax_eta/100.,dxymax_eta/100.);
865 
866  n_dzPhiResiduals[i] = NormLongPhiRes.make<TH1F>(Form("histo_norm_dz_phi_plot%i",i),Form("%.2f#circ<#varphi^{probe}_{tk}<%.2f#circ;d_{z}/#sigma_{d_{z}} [#mum];tracks",phiF,phiL),mybins_,-dzmax_phi/100.,dzmax_phi/100.);
867  n_dzEtaResiduals[i] = NormLongEtaRes.make<TH1F>(Form("histo_norm_dz_eta_plot%i",i),Form("%.2f<#eta^{probe}_{tk}<%.2f;d_{z}/#sigma_{d_{z}} [#mum];tracks",etaF,etaL),mybins_,-dzmax_eta/100.,dzmax_eta/100.);
868 
869  for ( int j=0; j<nBins_; ++j ) {
870 
871  a_dxyResidualsMap[i][j] = AbsDoubleDiffRes.make<TH1F>(Form("histo_dxy_eta_plot%i_phi_plot%i",i,j),Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{xy} [#mum];tracks",etaF,etaL,phiF,phiL),mybins_,-dzmax_eta,dzmax_eta);
872  a_dzResidualsMap[i][j] = AbsDoubleDiffRes.make<TH1F>(Form("histo_dxy_eta_plot%i_phi_plot%i",i,j),Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{z} [#mum];tracks",etaF,etaL,phiF,phiL),mybins_,-dzmax_eta,dzmax_eta);
873 
874  n_dxyResidualsMap[i][j] = NormDoubleDiffRes.make<TH1F>(Form("histo_norm_dxy_eta_plot%i_phi_plot%i",i,j),Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{xy}/#sigma_{d_{xy}} [#mum];tracks",etaF,etaL,phiF,phiL),mybins_,-dzmax_eta/100,dzmax_eta/100);
875  n_dzResidualsMap[i][j] = NormDoubleDiffRes.make<TH1F>(Form("histo_norm_dxy_eta_plot%i_phi_plot%i",i,j),Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{z}/#sigma_{d_{z}} [#mum];tracks",etaF,etaL,phiF,phiL),mybins_,-dzmax_eta/100,dzmax_eta/100);
876  }
877 
878  }
879 
880  // declaration of the directories
881 
882  TFileDirectory MeanTrendsDir = fs->mkdir("MeanTrends");
883  TFileDirectory WidthTrendsDir = fs->mkdir("WidthTrends");
884  TFileDirectory MedianTrendsDir = fs->mkdir("MedianTrends");
885  TFileDirectory MADTrendsDir = fs->mkdir("MADTrends");
886 
887  TFileDirectory Mean2DMapsDir = fs->mkdir("MeanMaps");
888  TFileDirectory Width2DMapsDir = fs->mkdir("WidthMaps");
889 
890  Double_t highedge=nBins_-0.5;
891  Double_t lowedge=-0.5;
892 
893  // means and widths from the fit
894 
895  a_dxyPhiMeanTrend = MeanTrendsDir.make<TH1F> ("means_dxy_phi","#LT d_{xy} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{xy} #GT [#mum]",nBins_,lowedge,highedge);
896  a_dxyPhiWidthTrend = WidthTrendsDir.make<TH1F>("widths_dxy_phi","#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees];#sigma_{d_{xy}} [#mum]",nBins_,lowedge,highedge);
897  a_dzPhiMeanTrend = MeanTrendsDir.make<TH1F> ("means_dz_phi","#LT d_{z} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{z} #GT [#mum]",nBins_,lowedge,highedge);
898  a_dzPhiWidthTrend = WidthTrendsDir.make<TH1F>("widths_dz_phi","#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];#sigma_{d_{z}} [#mum]",nBins_,lowedge,highedge);
899 
900  a_dxyEtaMeanTrend = MeanTrendsDir.make<TH1F> ("means_dxy_eta","#LT d_{xy} #GT vs #eta sector;#eta (sector);#LT d_{xy} #GT [#mum]",nBins_,lowedge,highedge);
901  a_dxyEtaWidthTrend = WidthTrendsDir.make<TH1F>("widths_dxy_eta","#sigma_{d_{xy}} vs #eta sector;#eta (sector);#sigma_{d_{xy}} [#mum]",nBins_,lowedge,highedge);
902  a_dzEtaMeanTrend = MeanTrendsDir.make<TH1F> ("means_dz_eta","#LT d_{z} #GT vs #eta sector;#eta (sector);#LT d_{z} #GT [#mum]",nBins_,lowedge,highedge);
903  a_dzEtaWidthTrend = WidthTrendsDir.make<TH1F>("widths_dz_eta","#sigma_{d_{xy}} vs #eta sector;#eta (sector);#sigma_{d_{z}} [#mum]",nBins_,lowedge,highedge);
904 
905  n_dxyPhiMeanTrend = MeanTrendsDir.make<TH1F> ("norm_means_dxy_phi","#LT d_{xy}/#sigma_{d_{xy}} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{xy}/#sigma_{d_{xy}} #GT",nBins_,lowedge,highedge);
906  n_dxyPhiWidthTrend = WidthTrendsDir.make<TH1F>("norm_widths_dxy_phi","width(d_{xy}/#sigma_{d_{xy}}) vs #phi sector;#varphi (sector) [degrees]; width(d_{xy}/#sigma_{d_{xy}})",nBins_,lowedge,highedge);
907  n_dzPhiMeanTrend = MeanTrendsDir.make<TH1F> ("norm_means_dz_phi","#LT d_{z}/#sigma_{d_{z}} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{z}/#sigma_{d_{z}} #GT",nBins_,lowedge,highedge);
908  n_dzPhiWidthTrend = WidthTrendsDir.make<TH1F>("norm_widths_dz_phi","width(d_{z}/#sigma_{d_{z}}) vs #phi sector;#varphi (sector) [degrees];width(d_{z}/#sigma_{d_{z}})",nBins_,lowedge,highedge);
909 
910  n_dxyEtaMeanTrend = MeanTrendsDir.make<TH1F> ("norm_means_dxy_eta","#LT d_{xy}/#sigma_{d_{xy}} #GT vs #eta sector;#eta (sector);#LT d_{xy}/#sigma_{d_{z}} #GT",nBins_,lowedge,highedge);
911  n_dxyEtaWidthTrend = WidthTrendsDir.make<TH1F>("norm_widths_dxy_eta","width(d_{xy}/#sigma_{d_{xy}}) vs #eta sector;#eta (sector);width(d_{xy}/#sigma_{d_{z}})",nBins_,lowedge,highedge);
912  n_dzEtaMeanTrend = MeanTrendsDir.make<TH1F> ("norm_means_dz_eta","#LT d_{z}/#sigma_{d_{z}} #GT vs #eta sector;#eta (sector);#LT d_{z}/#sigma_{d_{z}} #GT",nBins_,lowedge,highedge);
913  n_dzEtaWidthTrend = WidthTrendsDir.make<TH1F>("norm_widths_dz_eta","width(d_{z}/#sigma_{d_{z}}) vs #eta sector;#eta (sector);width(d_{z}/#sigma_{d_{z}})",nBins_,lowedge,highedge);
914 
915  // 2D maps
916 
917  a_dxyMeanMap = Mean2DMapsDir.make<TH2F> ("means_dxy_map","#LT d_{xy} #GT map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
918  a_dzMeanMap = Mean2DMapsDir.make<TH2F> ("means_dz_map","#LT d_{z} #GT map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
919 
920  n_dxyMeanMap = Mean2DMapsDir.make<TH2F> ("norm_means_dxy_map","#LT d_{xy}/#sigma_{d_{xy}} #GT map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
921  n_dzMeanMap = Mean2DMapsDir.make<TH2F> ("norm_means_dz_map","#LT d_{z}/#sigma_{d_{z}} #GT map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
922 
923  a_dxyWidthMap = Width2DMapsDir.make<TH2F> ("widths_dxy_map","#sigma_{d_{xy}} map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
924  a_dzWidthMap = Width2DMapsDir.make<TH2F> ("widths_dz_map","#sigma_{d_{z}} map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
925 
926  n_dxyWidthMap = Width2DMapsDir.make<TH2F> ("norm_widths_dxy_map","width(d_{xy}/#sigma_{d_{xy}}) map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
927  n_dzWidthMap = Width2DMapsDir.make<TH2F> ("norm_widths_dz_map","width(d_{z}/#sigma_{d_{z}}) map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
928 
929  // medians and MADs
930 
931  a_dxyPhiMedianTrend = MedianTrendsDir.make<TH1F>("medians_dxy_phi","Median of d_{xy} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{xy}) [#mum]",nBins_,lowedge,highedge);
932  a_dxyPhiMADTrend = MADTrendsDir.make<TH1F> ("MADs_dxy_phi","Median absolute deviation of d_{xy} vs #phi sector;#varphi (sector) [degrees];MAD(d_{xy}) [#mum]",nBins_,lowedge,highedge);
933  a_dzPhiMedianTrend = MedianTrendsDir.make<TH1F>("medians_dz_phi","Median of d_{z} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{z}) [#mum]",nBins_,lowedge,highedge);
934  a_dzPhiMADTrend = MADTrendsDir.make<TH1F> ("MADs_dz_phi","Median absolute deviation of d_{z} vs #phi sector;#varphi (sector) [degrees];MAD(d_{z}) [#mum]",nBins_,lowedge,highedge);
935 
936  a_dxyEtaMedianTrend = MedianTrendsDir.make<TH1F>("medians_dxy_eta","Median of d_{xy} vs #eta sector;#eta (sector);#mu_{1/2}(d_{xy}) [#mum]",nBins_,lowedge,highedge);
937  a_dxyEtaMADTrend = MADTrendsDir.make<TH1F> ("MADs_dxy_eta","Median absolute deviation of d_{xy} vs #eta sector;#eta (sector);MAD(d_{xy}) [#mum]",nBins_,lowedge,highedge);
938  a_dzEtaMedianTrend = MedianTrendsDir.make<TH1F>("medians_dz_eta","Median of d_{z} vs #eta sector;#eta (sector);#mu_{1/2}(d_{z}) [#mum]",nBins_,lowedge,highedge);
939  a_dzEtaMADTrend = MADTrendsDir.make<TH1F> ("MADs_dz_eta","Median absolute deviation of d_{z} vs #eta sector;#eta (sector);MAD(d_{z}) [#mum]",nBins_,lowedge,highedge);
940 
941  n_dxyPhiMedianTrend = MedianTrendsDir.make<TH1F>("norm_medians_dxy_phi","Median of d_{xy}/#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{xy}/#sigma_{d_{xy}})",nBins_,lowedge,highedge);
942  n_dxyPhiMADTrend = MADTrendsDir.make<TH1F> ("norm_MADs_dxy_phi","Median absolute deviation of d_{xy}/#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees]; MAD(d_{xy}/#sigma_{d_{xy}})",nBins_,lowedge,highedge);
943  n_dzPhiMedianTrend = MedianTrendsDir.make<TH1F>("norm_medians_dz_phi","Median of d_{z}/#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{z}/#sigma_{d_{z}})",nBins_,lowedge,highedge);
944  n_dzPhiMADTrend = MADTrendsDir.make<TH1F> ("norm_MADs_dz_phi","Median absolute deviation of d_{z}/#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];MAD(d_{z}/#sigma_{d_{z}})",nBins_,lowedge,highedge);
945 
946  n_dxyEtaMedianTrend = MedianTrendsDir.make<TH1F>("norm_medians_dxy_eta","Median of d_{xy}/#sigma_{d_{xy}} vs #eta sector;#eta (sector);#mu_{1/2}(d_{xy}/#sigma_{d_{z}})",nBins_,lowedge,highedge);
947  n_dxyEtaMADTrend = MADTrendsDir.make<TH1F> ("norm_MADs_dxy_eta","Median absolute deviation of d_{xy}/#sigma_{d_{xy}} vs #eta sector;#eta (sector);MAD(d_{xy}/#sigma_{d_{z}})",nBins_,lowedge,highedge);
948  n_dzEtaMedianTrend = MedianTrendsDir.make<TH1F>("norm_medians_dz_eta","Median of d_{z}/#sigma_{d_{z}} vs #eta sector;#eta (sector);#mu_{1/2}(d_{z}/#sigma_{d_{z}})",nBins_,lowedge,highedge);
949  n_dzEtaMADTrend = MADTrendsDir.make<TH1F> ("norm_MADs_dz_eta","Median absolute deviation of d_{z}/#sigma_{d_{z}} vs #eta sector;#eta (sector);MAD(d_{z}/#sigma_{d_{z}})",nBins_,lowedge,highedge);
950 
951 
953  //
954  // plots of biased residuals
955  // The vertex is refit without the probe track
956  //
958 
959  if (useTracksFromRecoVtx_){
960 
961  TFileDirectory AbsTransPhiBiasRes = fs->mkdir("Abs_Transv_Phi_BiasResiduals");
962  TFileDirectory AbsTransEtaBiasRes = fs->mkdir("Abs_Transv_Eta_BiasResiduals");
963 
964  TFileDirectory AbsLongPhiBiasRes = fs->mkdir("Abs_Long_Phi_BiasResiduals");
965  TFileDirectory AbsLongEtaBiasRes = fs->mkdir("Abs_Long_Eta_BiasResiduals");
966 
967  TFileDirectory NormTransPhiBiasRes = fs->mkdir("Norm_Transv_Phi_BiasResiduals");
968  TFileDirectory NormTransEtaBiasRes = fs->mkdir("Norm_Transv_Eta_BiasResiduals");
969 
970  TFileDirectory NormLongPhiBiasRes = fs->mkdir("Norm_Long_Phi_BiasResiduals");
971  TFileDirectory NormLongEtaBiasRes = fs->mkdir("Norm_Long_Eta_BiasResiduals");
972 
973  TFileDirectory AbsDoubleDiffBiasRes = fs->mkdir("Abs_DoubleDiffBiasResiduals");
974  TFileDirectory NormDoubleDiffBiasRes = fs->mkdir("Norm_DoubleDiffBiasResiduals");
975 
976  for ( int i=0; i<nBins_; ++i ) {
977 
978  float phiF = (-TMath::Pi()+i*phipitch_)*(180/TMath::Pi());
979  float phiL = (-TMath::Pi()+(i+1)*phipitch_)*(180/TMath::Pi());
980 
981  float etaF=-2.5+i*etapitch_;
982  float etaL=-2.5+(i+1)*etapitch_;
983 
984  a_dxyPhiBiasResiduals[i] = AbsTransPhiBiasRes.make<TH1F>(Form("histo_dxy_phi_plot%i",i),Form("%.2f#circ<#varphi^{probe}_{tk}<%.2f#circ;d_{xy} [#mum];tracks",phiF,phiL),mybins_,-dxymax_phi,dxymax_phi);
985  a_dxyEtaBiasResiduals[i] = AbsTransEtaBiasRes.make<TH1F>(Form("histo_dxy_eta_plot%i",i),Form("%.2f<#eta^{probe}_{tk}<%.2f;d_{xy} [#mum];tracks",etaF,etaL),mybins_,-dxymax_eta,dxymax_eta);
986 
987  a_dzPhiBiasResiduals[i] = AbsLongPhiBiasRes.make<TH1F>(Form("histo_dz_phi_plot%i",i),Form("%.2f#circ<#varphi^{probe}_{tk}<%.2f #circ;d_{z} [#mum];tracks",phiF,phiL),mybins_,-dzmax_phi,dzmax_phi);
988  a_dzEtaBiasResiduals[i] = AbsLongEtaBiasRes.make<TH1F>(Form("histo_dz_eta_plot%i",i),Form("%.2f<#eta^{probe}_{tk}<%.2f;d_{z} [#mum];tracks",etaF,etaL),mybins_,-dzmax_eta,dzmax_eta);
989 
990  n_dxyPhiBiasResiduals[i] = NormTransPhiBiasRes.make<TH1F>(Form("histo_norm_dxy_phi_plot%i",i),Form("%.2f#circ<#varphi^{probe}_{tk}<%.2f#circ;d_{xy}/#sigma_{d_{xy}} [#mum];tracks",phiF,phiL),mybins_,-dxymax_phi/100.,dxymax_phi/100.);
991  n_dxyEtaBiasResiduals[i] = NormTransEtaBiasRes.make<TH1F>(Form("histo_norm_dxy_eta_plot%i",i),Form("%.2f<#eta^{probe}_{tk}<%.2f;d_{xy}/#sigma_{d_{xy}} [#mum];tracks",etaF,etaL),mybins_,-dxymax_eta/100.,dxymax_eta/100.);
992 
993  n_dzPhiBiasResiduals[i] = NormLongPhiBiasRes.make<TH1F>(Form("histo_norm_dz_phi_plot%i",i),Form("%.2f#circ<#varphi^{probe}_{tk}<%.2f#circ;d_{z}/#sigma_{d_{z}} [#mum];tracks",phiF,phiL),mybins_,-dzmax_phi/100.,dzmax_phi/100.);
994  n_dzEtaBiasResiduals[i] = NormLongEtaBiasRes.make<TH1F>(Form("histo_norm_dz_eta_plot%i",i),Form("%.2f<#eta^{probe}_{tk}<%.2f;d_{z}/#sigma_{d_{z}} [#mum];tracks",etaF,etaL),mybins_,-dzmax_eta/100.,dzmax_eta/100.);
995 
996  for ( int j=0; j<nBins_; ++j ) {
997 
998  a_dxyBiasResidualsMap[i][j] = AbsDoubleDiffBiasRes.make<TH1F>(Form("histo_dxy_eta_plot%i_phi_plot%i",i,j),Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{xy} [#mum];tracks",etaF,etaL,phiF,phiL),mybins_,-dzmax_eta,dzmax_eta);
999  a_dzBiasResidualsMap[i][j] = AbsDoubleDiffBiasRes.make<TH1F>(Form("histo_dxy_eta_plot%i_phi_plot%i",i,j),Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{z} [#mum];tracks",etaF,etaL,phiF,phiL),mybins_,-dzmax_eta,dzmax_eta);
1000 
1001  n_dxyBiasResidualsMap[i][j] = NormDoubleDiffBiasRes.make<TH1F>(Form("histo_norm_dxy_eta_plot%i_phi_plot%i",i,j),Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{xy}/#sigma_{d_{xy}} [#mum];tracks",etaF,etaL,phiF,phiL),mybins_,-dzmax_eta/100,dzmax_eta/100);
1002  n_dzBiasResidualsMap[i][j] = NormDoubleDiffBiasRes.make<TH1F>(Form("histo_norm_dxy_eta_plot%i_phi_plot%i",i,j),Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{z}/#sigma_{d_{z}} [#mum];tracks",etaF,etaL,phiF,phiL),mybins_,-dzmax_eta/100,dzmax_eta/100);
1003  }
1004 
1005  }
1006 
1007  // declaration of the directories
1008 
1009  TFileDirectory MeanBiasTrendsDir = fs->mkdir("MeanBiasTrends");
1010  TFileDirectory WidthBiasTrendsDir = fs->mkdir("WidthBiasTrends");
1011  TFileDirectory MedianBiasTrendsDir = fs->mkdir("MedianBiasTrends");
1012  TFileDirectory MADBiasTrendsDir = fs->mkdir("MADBiasTrends");
1013 
1014  TFileDirectory Mean2DBiasMapsDir = fs->mkdir("MeanBiasMaps");
1015  TFileDirectory Width2DBiasMapsDir = fs->mkdir("WidthBiasMaps");
1016 
1017  // means and widths from the fit
1018 
1019  a_dxyPhiMeanBiasTrend = MeanBiasTrendsDir.make<TH1F> ("means_dxy_phi","#LT d_{xy} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{xy} #GT [#mum]",nBins_,lowedge,highedge);
1020  a_dxyPhiWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>("widths_dxy_phi","#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees];#sigma_{d_{xy}} [#mum]",nBins_,lowedge,highedge);
1021  a_dzPhiMeanBiasTrend = MeanBiasTrendsDir.make<TH1F> ("means_dz_phi","#LT d_{z} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{z} #GT [#mum]",nBins_,lowedge,highedge);
1022  a_dzPhiWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>("widths_dz_phi","#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];#sigma_{d_{z}} [#mum]",nBins_,lowedge,highedge);
1023 
1024  a_dxyEtaMeanBiasTrend = MeanBiasTrendsDir.make<TH1F> ("means_dxy_eta","#LT d_{xy} #GT vs #eta sector;#eta (sector);#LT d_{xy} #GT [#mum]",nBins_,lowedge,highedge);
1025  a_dxyEtaWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>("widths_dxy_eta","#sigma_{d_{xy}} vs #eta sector;#eta (sector);#sigma_{d_{xy}} [#mum]",nBins_,lowedge,highedge);
1026  a_dzEtaMeanBiasTrend = MeanBiasTrendsDir.make<TH1F> ("means_dz_eta","#LT d_{z} #GT vs #eta sector;#eta (sector);#LT d_{z} #GT [#mum]",nBins_,lowedge,highedge);
1027  a_dzEtaWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>("widths_dz_eta","#sigma_{d_{xy}} vs #eta sector;#eta (sector);#sigma_{d_{z}} [#mum]",nBins_,lowedge,highedge);
1028 
1029  n_dxyPhiMeanBiasTrend = MeanBiasTrendsDir.make<TH1F> ("norm_means_dxy_phi","#LT d_{xy}/#sigma_{d_{xy}} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{xy}/#sigma_{d_{xy}} #GT",nBins_,lowedge,highedge);
1030  n_dxyPhiWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>("norm_widths_dxy_phi","width(d_{xy}/#sigma_{d_{xy}}) vs #phi sector;#varphi (sector) [degrees]; width(d_{xy}/#sigma_{d_{xy}})",nBins_,lowedge,highedge);
1031  n_dzPhiMeanBiasTrend = MeanBiasTrendsDir.make<TH1F> ("norm_means_dz_phi","#LT d_{z}/#sigma_{d_{z}} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{z}/#sigma_{d_{z}} #GT",nBins_,lowedge,highedge);
1032  n_dzPhiWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>("norm_widths_dz_phi","width(d_{z}/#sigma_{d_{z}}) vs #phi sector;#varphi (sector) [degrees];width(d_{z}/#sigma_{d_{z}})",nBins_,lowedge,highedge);
1033 
1034  n_dxyEtaMeanBiasTrend = MeanBiasTrendsDir.make<TH1F> ("norm_means_dxy_eta","#LT d_{xy}/#sigma_{d_{xy}} #GT vs #eta sector;#eta (sector);#LT d_{xy}/#sigma_{d_{z}} #GT",nBins_,lowedge,highedge);
1035  n_dxyEtaWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>("norm_widths_dxy_eta","width(d_{xy}/#sigma_{d_{xy}}) vs #eta sector;#eta (sector);width(d_{xy}/#sigma_{d_{z}})",nBins_,lowedge,highedge);
1036  n_dzEtaMeanBiasTrend = MeanBiasTrendsDir.make<TH1F> ("norm_means_dz_eta","#LT d_{z}/#sigma_{d_{z}} #GT vs #eta sector;#eta (sector);#LT d_{z}/#sigma_{d_{z}} #GT",nBins_,lowedge,highedge);
1037  n_dzEtaWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>("norm_widths_dz_eta","width(d_{z}/#sigma_{d_{z}}) vs #eta sector;#eta (sector);width(d_{z}/#sigma_{d_{z}})",nBins_,lowedge,highedge);
1038 
1039  // 2D maps
1040 
1041  a_dxyMeanBiasMap = Mean2DBiasMapsDir.make<TH2F> ("means_dxy_map","#LT d_{xy} #GT map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1042  a_dzMeanBiasMap = Mean2DBiasMapsDir.make<TH2F> ("means_dz_map","#LT d_{z} #GT map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1043 
1044  n_dxyMeanBiasMap = Mean2DBiasMapsDir.make<TH2F> ("norm_means_dxy_map","#LT d_{xy}/#sigma_{d_{xy}} #GT map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1045  n_dzMeanBiasMap = Mean2DBiasMapsDir.make<TH2F> ("norm_means_dz_map","#LT d_{z}/#sigma_{d_{z}} #GT map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1046 
1047  a_dxyWidthBiasMap = Width2DBiasMapsDir.make<TH2F> ("widths_dxy_map","#sigma_{d_{xy}} map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1048  a_dzWidthBiasMap = Width2DBiasMapsDir.make<TH2F> ("widths_dz_map","#sigma_{d_{z}} map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1049 
1050  n_dxyWidthBiasMap = Width2DBiasMapsDir.make<TH2F> ("norm_widths_dxy_map","width(d_{xy}/#sigma_{d_{xy}}) map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1051  n_dzWidthBiasMap = Width2DBiasMapsDir.make<TH2F> ("norm_widths_dz_map","width(d_{z}/#sigma_{d_{z}}) map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1052 
1053  // medians and MADs
1054 
1055  a_dxyPhiMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>("medians_dxy_phi","Median of d_{xy} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{xy}) [#mum]",nBins_,lowedge,highedge);
1056  a_dxyPhiMADBiasTrend = MADBiasTrendsDir.make<TH1F> ("MADs_dxy_phi","Median absolute deviation of d_{xy} vs #phi sector;#varphi (sector) [degrees];MAD(d_{xy}) [#mum]",nBins_,lowedge,highedge);
1057  a_dzPhiMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>("medians_dz_phi","Median of d_{z} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{z}) [#mum]",nBins_,lowedge,highedge);
1058  a_dzPhiMADBiasTrend = MADBiasTrendsDir.make<TH1F> ("MADs_dz_phi","Median absolute deviation of d_{z} vs #phi sector;#varphi (sector) [degrees];MAD(d_{z}) [#mum]",nBins_,lowedge,highedge);
1059 
1060  a_dxyEtaMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>("medians_dxy_eta","Median of d_{xy} vs #eta sector;#eta (sector);#mu_{1/2}(d_{xy}) [#mum]",nBins_,lowedge,highedge);
1061  a_dxyEtaMADBiasTrend = MADBiasTrendsDir.make<TH1F> ("MADs_dxy_eta","Median absolute deviation of d_{xy} vs #eta sector;#eta (sector);MAD(d_{xy}) [#mum]",nBins_,lowedge,highedge);
1062  a_dzEtaMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>("medians_dz_eta","Median of d_{z} vs #eta sector;#eta (sector);#mu_{1/2}(d_{z}) [#mum]",nBins_,lowedge,highedge);
1063  a_dzEtaMADBiasTrend = MADBiasTrendsDir.make<TH1F> ("MADs_dz_eta","Median absolute deviation of d_{z} vs #eta sector;#eta (sector);MAD(d_{z}) [#mum]",nBins_,lowedge,highedge);
1064 
1065  n_dxyPhiMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>("norm_medians_dxy_phi","Median of d_{xy}/#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{xy}/#sigma_{d_{xy}})",nBins_,lowedge,highedge);
1066  n_dxyPhiMADBiasTrend = MADBiasTrendsDir.make<TH1F> ("norm_MADs_dxy_phi","Median absolute deviation of d_{xy}/#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees]; MAD(d_{xy}/#sigma_{d_{xy}})",nBins_,lowedge,highedge);
1067  n_dzPhiMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>("norm_medians_dz_phi","Median of d_{z}/#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{z}/#sigma_{d_{z}})",nBins_,lowedge,highedge);
1068  n_dzPhiMADBiasTrend = MADBiasTrendsDir.make<TH1F> ("norm_MADs_dz_phi","Median absolute deviation of d_{z}/#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];MAD(d_{z}/#sigma_{d_{z}})",nBins_,lowedge,highedge);
1069 
1070  n_dxyEtaMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>("norm_medians_dxy_eta","Median of d_{xy}/#sigma_{d_{xy}} vs #eta sector;#eta (sector);#mu_{1/2}(d_{xy}/#sigma_{d_{z}})",nBins_,lowedge,highedge);
1071  n_dxyEtaMADBiasTrend = MADBiasTrendsDir.make<TH1F> ("norm_MADs_dxy_eta","Median absolute deviation of d_{xy}/#sigma_{d_{xy}} vs #eta sector;#eta (sector);MAD(d_{xy}/#sigma_{d_{z}})",nBins_,lowedge,highedge);
1072  n_dzEtaMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>("norm_medians_dz_eta","Median of d_{z}/#sigma_{d_{z}} vs #eta sector;#eta (sector);#mu_{1/2}(d_{z}/#sigma_{d_{z}})",nBins_,lowedge,highedge);
1073  n_dzEtaMADBiasTrend = MADBiasTrendsDir.make<TH1F> ("norm_MADs_dz_eta","Median absolute deviation of d_{z}/#sigma_{d_{z}} vs #eta sector;#eta (sector);MAD(d_{z}/#sigma_{d_{z}})",nBins_,lowedge,highedge);
1074 
1075  }
1076 }
1077 // ------------ method called once each job just after ending the event loop ------------
1078 void PrimaryVertexValidation::endJob()
1079 {
1080 
1081  std::cout<<"######################################"<<std::endl;
1082  std::cout<<"Number of analyzed events: "<<Nevt_<<std::endl;
1083  std::cout<<"######################################"<<std::endl;
1084 
1085  if(useTracksFromRecoVtx_){
1086 
1087  FillTrendPlot(a_dxyPhiMeanBiasTrend ,a_dxyPhiBiasResiduals,"mean","phi");
1088  FillTrendPlot(a_dxyPhiWidthBiasTrend,a_dxyPhiBiasResiduals,"width","phi");
1089  FillTrendPlot(a_dzPhiMeanBiasTrend ,a_dzPhiBiasResiduals ,"mean","phi");
1090  FillTrendPlot(a_dzPhiWidthBiasTrend ,a_dzPhiBiasResiduals ,"width","phi");
1091 
1092  FillTrendPlot(a_dxyEtaMeanBiasTrend ,a_dxyEtaBiasResiduals,"mean","eta");
1093  FillTrendPlot(a_dxyEtaWidthBiasTrend,a_dxyEtaBiasResiduals,"width","eta");
1094  FillTrendPlot(a_dzEtaMeanBiasTrend ,a_dzEtaBiasResiduals ,"mean","eta");
1095  FillTrendPlot(a_dzEtaWidthBiasTrend ,a_dzEtaBiasResiduals ,"width","eta");
1096 
1097  FillTrendPlot(n_dxyPhiMeanBiasTrend ,n_dxyPhiBiasResiduals,"mean","phi");
1098  FillTrendPlot(n_dxyPhiWidthBiasTrend,n_dxyPhiBiasResiduals,"width","phi");
1099  FillTrendPlot(n_dzPhiMeanBiasTrend ,n_dzPhiBiasResiduals ,"mean","phi");
1100  FillTrendPlot(n_dzPhiWidthBiasTrend ,n_dzPhiBiasResiduals ,"width","phi");
1101 
1102  FillTrendPlot(n_dxyEtaMeanBiasTrend ,n_dxyEtaBiasResiduals,"mean","eta");
1103  FillTrendPlot(n_dxyEtaWidthBiasTrend,n_dxyEtaBiasResiduals,"width","eta");
1104  FillTrendPlot(n_dzEtaMeanBiasTrend ,n_dzEtaBiasResiduals ,"mean","eta");
1105  FillTrendPlot(n_dzEtaWidthBiasTrend ,n_dzEtaBiasResiduals ,"width","eta");
1106 
1107  // medians and MADs
1108 
1109  FillTrendPlot(a_dxyPhiMedianBiasTrend,a_dxyPhiBiasResiduals,"median","phi");
1110  FillTrendPlot(a_dxyPhiMADBiasTrend ,a_dxyPhiBiasResiduals,"mad","phi");
1111  FillTrendPlot(a_dzPhiMedianBiasTrend ,a_dzPhiBiasResiduals ,"median","phi");
1112  FillTrendPlot(a_dzPhiMADBiasTrend ,a_dzPhiBiasResiduals ,"mad","phi");
1113 
1114  FillTrendPlot(a_dxyEtaMedianBiasTrend,a_dxyEtaBiasResiduals,"median","eta");
1115  FillTrendPlot(a_dxyEtaMADBiasTrend ,a_dxyEtaBiasResiduals,"mad","eta");
1116  FillTrendPlot(a_dzEtaMedianBiasTrend ,a_dzEtaBiasResiduals ,"median","eta");
1117  FillTrendPlot(a_dzEtaMADBiasTrend ,a_dzEtaBiasResiduals ,"mad","eta");
1118 
1119  FillTrendPlot(n_dxyPhiMedianBiasTrend,n_dxyPhiBiasResiduals,"median","phi");
1120  FillTrendPlot(n_dxyPhiMADBiasTrend ,n_dxyPhiBiasResiduals,"mad","phi");
1121  FillTrendPlot(n_dzPhiMedianBiasTrend ,n_dzPhiBiasResiduals ,"median","phi");
1122  FillTrendPlot(n_dzPhiMADBiasTrend ,n_dzPhiBiasResiduals ,"mad","phi");
1123 
1124  FillTrendPlot(n_dxyEtaMedianBiasTrend,n_dxyEtaBiasResiduals,"median","eta");
1125  FillTrendPlot(n_dxyEtaMADBiasTrend ,n_dxyEtaBiasResiduals,"mad","eta");
1126  FillTrendPlot(n_dzEtaMedianBiasTrend ,n_dzEtaBiasResiduals ,"median","eta");
1127  FillTrendPlot(n_dzEtaMADBiasTrend ,n_dzEtaBiasResiduals ,"mad","eta");
1128 
1129  }
1130 
1131  FillTrendPlot(a_dxyPhiMeanTrend ,a_dxyPhiResiduals,"mean","phi");
1132  FillTrendPlot(a_dxyPhiWidthTrend,a_dxyPhiResiduals,"width","phi");
1133  FillTrendPlot(a_dzPhiMeanTrend ,a_dzPhiResiduals ,"mean","phi");
1134  FillTrendPlot(a_dzPhiWidthTrend ,a_dzPhiResiduals ,"width","phi");
1135 
1136  FillTrendPlot(a_dxyEtaMeanTrend ,a_dxyEtaResiduals,"mean","eta");
1137  FillTrendPlot(a_dxyEtaWidthTrend,a_dxyEtaResiduals,"width","eta");
1138  FillTrendPlot(a_dzEtaMeanTrend ,a_dzEtaResiduals ,"mean","eta");
1139  FillTrendPlot(a_dzEtaWidthTrend ,a_dzEtaResiduals ,"width","eta");
1140 
1141  FillTrendPlot(n_dxyPhiMeanTrend ,n_dxyPhiResiduals,"mean","phi");
1142  FillTrendPlot(n_dxyPhiWidthTrend,n_dxyPhiResiduals,"width","phi");
1143  FillTrendPlot(n_dzPhiMeanTrend ,n_dzPhiResiduals ,"mean","phi");
1144  FillTrendPlot(n_dzPhiWidthTrend ,n_dzPhiResiduals ,"width","phi");
1145 
1146  FillTrendPlot(n_dxyEtaMeanTrend ,n_dxyEtaResiduals,"mean","eta");
1147  FillTrendPlot(n_dxyEtaWidthTrend,n_dxyEtaResiduals,"width","eta");
1148  FillTrendPlot(n_dzEtaMeanTrend ,n_dzEtaResiduals ,"mean","eta");
1149  FillTrendPlot(n_dzEtaWidthTrend ,n_dzEtaResiduals ,"width","eta");
1150 
1151  // medians and MADs
1152 
1153  FillTrendPlot(a_dxyPhiMedianTrend,a_dxyPhiResiduals,"median","phi");
1154  FillTrendPlot(a_dxyPhiMADTrend ,a_dxyPhiResiduals,"mad","phi");
1155  FillTrendPlot(a_dzPhiMedianTrend ,a_dzPhiResiduals ,"median","phi");
1156  FillTrendPlot(a_dzPhiMADTrend ,a_dzPhiResiduals ,"mad","phi");
1157 
1158  FillTrendPlot(a_dxyEtaMedianTrend,a_dxyEtaResiduals,"median","eta");
1159  FillTrendPlot(a_dxyEtaMADTrend ,a_dxyEtaResiduals,"mad","eta");
1160  FillTrendPlot(a_dzEtaMedianTrend ,a_dzEtaResiduals ,"median","eta");
1161  FillTrendPlot(a_dzEtaMADTrend ,a_dzEtaResiduals ,"mad","eta");
1162 
1163  FillTrendPlot(n_dxyPhiMedianTrend,n_dxyPhiResiduals,"median","phi");
1164  FillTrendPlot(n_dxyPhiMADTrend ,n_dxyPhiResiduals,"mad","phi");
1165  FillTrendPlot(n_dzPhiMedianTrend ,n_dzPhiResiduals ,"median","phi");
1166  FillTrendPlot(n_dzPhiMADTrend ,n_dzPhiResiduals ,"mad","phi");
1167 
1168  FillTrendPlot(n_dxyEtaMedianTrend,n_dxyEtaResiduals,"median","eta");
1169  FillTrendPlot(n_dxyEtaMADTrend ,n_dxyEtaResiduals,"mad","eta");
1170  FillTrendPlot(n_dzEtaMedianTrend ,n_dzEtaResiduals ,"median","eta");
1171  FillTrendPlot(n_dzEtaMADTrend ,n_dzEtaResiduals ,"mad","eta");
1172 
1173 }
1174 
1175 //*************************************************************
1176 void PrimaryVertexValidation::SetVarToZero()
1177 //*************************************************************
1178 {
1179  nTracks_ = 0;
1180  nClus_ = 0;
1181  nOfflineVertices_=0;
1182  RunNumber_ =0;
1183  LuminosityBlockNumber_=0;
1184  xOfflineVertex_ =-999.;
1185  yOfflineVertex_ =-999.;
1186  zOfflineVertex_ =-999.;
1187  BSx0_ = -999.;
1188  BSy0_ = -999.;
1189  BSz0_ = -999.;
1190  Beamsigmaz_=-999.;
1191  Beamdxdz_=-999.;
1192  BeamWidthX_=-999.;
1193  BeamWidthY_=-999.;
1194 
1195  for ( int i=0; i<nMaxtracks_; ++i ) {
1196 
1197  pt_[i] = 0;
1198  p_[i] = 0;
1199  nhits_[i] = 0;
1200  nhits1D_[i] = 0;
1201  nhits2D_[i] = 0;
1202  nhitsBPIX_[i] = 0;
1203  nhitsFPIX_[i] = 0;
1204  nhitsTIB_[i] = 0;
1205  nhitsTID_[i] = 0;
1206  nhitsTOB_[i] = 0;
1207  nhitsTEC_[i] = 0;
1208  isHighPurity_[i] = 0;
1209  eta_[i] = 0;
1210  theta_[i] = 0;
1211  phi_[i] = 0;
1212  chi2_[i] = 0;
1213  chi2ndof_[i] = 0;
1214  charge_[i] = 0;
1215  qoverp_[i] = 0;
1216  dz_[i] = 0;
1217  dxy_[i] = 0;
1218  dzBs_[i] = 0;
1219  dxyBs_[i] = 0;
1220  xPCA_[i] = 0;
1221  yPCA_[i] = 0;
1222  zPCA_[i] = 0;
1223  xUnbiasedVertex_[i] =0;
1224  yUnbiasedVertex_[i] =0;
1225  zUnbiasedVertex_[i] =0;
1226  chi2normUnbiasedVertex_[i]=0;
1227  chi2UnbiasedVertex_[i]=0;
1228  DOFUnbiasedVertex_[i]=0;
1229  sumOfWeightsUnbiasedVertex_[i]=0;
1230  tracksUsedForVertexing_[i]=0;
1231  dxyFromMyVertex_[i]=0;
1232  dzFromMyVertex_[i]=0;
1233  dszFromMyVertex_[i]=0;
1234  dxyErrorFromMyVertex_[i]=0;
1235  dzErrorFromMyVertex_[i]=0;
1236  IPTsigFromMyVertex_[i]=0;
1237  IPLsigFromMyVertex_[i]=0;
1238  hasRecVertex_[i] = 0;
1239  isGoodTrack_[i] = 0;
1240  }
1241 }
1242 
1243 //*************************************************************
1244 std::pair<Double_t,Double_t> PrimaryVertexValidation::getMedian(TH1F *histo)
1245 //*************************************************************
1246 {
1247  Double_t median = 999;
1248  int nbins = histo->GetNbinsX();
1249 
1250  //extract median from histogram
1251  double *x = new double[nbins];
1252  double *y = new double[nbins];
1253  for (int j = 0; j < nbins; j++) {
1254  x[j] = histo->GetBinCenter(j+1);
1255  y[j] = histo->GetBinContent(j+1);
1256  }
1257  median = TMath::Median(nbins, x, y);
1258 
1259  delete[] x; x = 0;
1260  delete [] y; y = 0;
1261 
1262  std::pair<Double_t,Double_t> result;
1263  result = std::make_pair(median,median/TMath::Sqrt(histo->GetEntries()));
1264 
1265  return result;
1266 
1267 }
1268 
1269 //*************************************************************
1270 std::pair<Double_t,Double_t> PrimaryVertexValidation::getMAD(TH1F *histo)
1271 //*************************************************************
1272 {
1273 
1274  int nbins = histo->GetNbinsX();
1275  Double_t median = getMedian(histo).first;
1276  Double_t x_lastBin = histo->GetBinLowEdge(nbins+1);
1277  const char *HistoName =histo->GetName();
1278  TString Finalname = Form("resMed%s",HistoName);
1279  TH1F *newHisto = new TH1F(Finalname,Finalname,nbins,0.,x_lastBin);
1280  Double_t *residuals = new Double_t[nbins];
1281  Double_t *weights = new Double_t[nbins];
1282 
1283  for (int j = 0; j < nbins; j++) {
1284  residuals[j] = TMath::Abs(median - histo->GetBinCenter(j+1));
1285  weights[j]=histo->GetBinContent(j+1);
1286  newHisto->Fill(residuals[j],weights[j]);
1287  }
1288 
1289  Double_t theMAD = (getMedian(newHisto).first)*1.4826;
1290  newHisto->Delete("");
1291 
1292  std::pair<Double_t,Double_t> result;
1293  result = std::make_pair(theMAD,theMAD/histo->GetEntries());
1294 
1295  return result;
1296 
1297 }
1298 
1299 //*************************************************************
1300 std::pair<std::pair<Double_t,Double_t>, std::pair<Double_t,Double_t> > PrimaryVertexValidation::fitResiduals(TH1 *hist)
1301 //*************************************************************
1302 {
1303  //float fitResult(9999);
1304  //if (hist->GetEntries() < 20) return ;
1305 
1306  float mean = hist->GetMean();
1307  float sigma = hist->GetRMS();
1308 
1309  TF1 func("tmp", "gaus", mean - 1.5*sigma, mean + 1.5*sigma);
1310  if (0 == hist->Fit(&func,"QNR")) { // N: do not blow up file by storing fit!
1311  mean = func.GetParameter(1);
1312  sigma = func.GetParameter(2);
1313  // second fit: three sigma of first fit around mean of first fit
1314  func.SetRange(mean - 2*sigma, mean + 2*sigma);
1315  // I: integral gives more correct results if binning is too wide
1316  // L: Likelihood can treat empty bins correctly (if hist not weighted...)
1317  if (0 == hist->Fit(&func, "Q0LR")) {
1318  if (hist->GetFunction(func.GetName())) { // Take care that it is later on drawn:
1319  hist->GetFunction(func.GetName())->ResetBit(TF1::kNotDraw);
1320  }
1321  }
1322  }
1323 
1324  float res_mean = func.GetParameter(1);
1325  float res_width = func.GetParameter(2);
1326 
1327  float res_mean_err = func.GetParError(1);
1328  float res_width_err = func.GetParError(2);
1329 
1330  std::pair<Double_t,Double_t> resultM;
1331  std::pair<Double_t,Double_t> resultW;
1332 
1333  resultM = std::make_pair(res_mean,res_mean_err);
1334  resultW = std::make_pair(res_width,res_width_err);
1335 
1336  std::pair<std::pair<Double_t,Double_t>, std::pair<Double_t,Double_t> > result;
1337 
1338  result = std::make_pair(resultM,resultW);
1339  return result;
1340 }
1341 
1342 //*************************************************************
1343 void PrimaryVertexValidation::FillTrendPlot(TH1F* trendPlot, TH1F* residualsPlot[100], const TString& fitPar_, const TString& var_)
1344 //*************************************************************
1345 {
1346 
1347  float phiInterval = (360.)/nBins_;
1348  float etaInterval = 5./nBins_;
1349 
1350  for ( int i=0; i<nBins_; ++i ) {
1351 
1352  char phipositionString[129];
1353  float phiposition = (-180+i*phiInterval)+(phiInterval/2);
1354  sprintf(phipositionString,"%.f",phiposition);
1355 
1356  char etapositionString[129];
1357  float etaposition = (-2.5+i*etaInterval)+(etaInterval/2);
1358  sprintf(etapositionString,"%.1f",etaposition);
1359 
1360  if(fitPar_=="mean"){
1361  float mean_ = fitResiduals(residualsPlot[i]).first.first;
1362  float meanErr_ = fitResiduals(residualsPlot[i]).first.second;
1363  trendPlot->SetBinContent(i+1,mean_);
1364  trendPlot->SetBinError(i+1,meanErr_);
1365  } else if (fitPar_=="width"){
1366  float width_ = fitResiduals(residualsPlot[i]).second.first;
1367  float widthErr_ = fitResiduals(residualsPlot[i]).second.second;
1368  trendPlot->SetBinContent(i+1,width_);
1369  trendPlot->SetBinError(i+1,widthErr_);
1370  } else if (fitPar_=="median"){
1371  float median_ = getMedian(residualsPlot[i]).first;
1372  float medianErr_ = getMedian(residualsPlot[i]).second;
1373  trendPlot->SetBinContent(i+1,median_);
1374  trendPlot->SetBinError(i+1,medianErr_);
1375  } else if (fitPar_=="mad"){
1376  float mad_ = getMAD(residualsPlot[i]).first;
1377  float madErr_ = getMAD(residualsPlot[i]).second;
1378  trendPlot->SetBinContent(i+1,mad_);
1379  trendPlot->SetBinError(i+1,madErr_);
1380  } else {
1381  std::cout<<"PrimaryVertexValidation::FillTrendPlot() "<<fitPar_<<" unknown estimator!"<<std::endl;
1382  }
1383 
1384  if(var_=="eta"){
1385  trendPlot->GetXaxis()->SetBinLabel(i+1,phipositionString);
1386  } else if(var_=="phi"){
1387  trendPlot->GetXaxis()->SetBinLabel(i+1,etapositionString);
1388  } else {
1389  std::cout<<"PrimaryVertexValidation::FillTrendPlot() "<<var_<<" unknown track parameter!"<<std::endl;
1390  }
1391  }
1392 }
1393 
1394 
1395 
1396 //define this as a plug-in
1397 DEFINE_FWK_MODULE(PrimaryVertexValidation);
cms::Exception::what
virtual char const * what() const
Definition: Exception.cc:141
TransientVertex::positionError
GlobalError positionError() const
Definition: TransientVertex.h:115
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
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:226
SiPixelRawToDigiRegional_cfi.beamSpot
tuple beamSpot
Definition: SiPixelRawToDigiRegional_cfi.py:11
reco::Vertex::tracks_end
trackRef_iterator tracks_end() const
last iterator over tracks
Definition: Vertex.cc:44
PrimaryVertexValidation::theTrackFilter_
TrackFilterForPVFindingBase * theTrackFilter_
Definition: PrimaryVertexValidation.h:101
PrimaryVertexValidation::a_dxyResidualsMap
TH1F * a_dxyResidualsMap[nMaxBins_][nMaxBins_]
Definition: PrimaryVertexValidation.h:222
reco::TrackBase::TrackQuality
TrackQuality
track quality
Definition: TrackBase.h:151
PrimaryVertexValidation::n_dxyEtaBiasResiduals
TH1F * n_dxyEtaBiasResiduals[nMaxBins_]
Definition: PrimaryVertexValidation.h:299
VertexException
Common base class.
Definition: VertexException.h:12
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
PrimaryVertexValidation::PrimaryVertexValidation
PrimaryVertexValidation(const edm::ParameterSet &)
Definition: PrimaryVertexValidation.cc:75
PrimaryVertexValidation::getMAD
std::pair< Double_t, Double_t > getMAD(TH1F *histo)
Definition: PrimaryVertexValidation.cc:1270
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:309
TransientVertex::totalChiSquared
float totalChiSquared() const
Definition: TransientVertex.h:125
reco::TransientTrack::hitPattern
const HitPattern & hitPattern() const
Definition: TransientTrack.h:128
PV3DBase::y
T y() const
Definition: PV3DBase.h:63
edm::EventAuxiliary::run
RunNumber_t run() const
Definition: EventAuxiliary.h:67
PrimaryVertexValidation::DOFUnbiasedVertex_
float DOFUnbiasedVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:185
PrimaryVertexValidation::a_dzPhiBiasResiduals
TH1F * a_dzPhiBiasResiduals[nMaxBins_]
Definition: PrimaryVertexValidation.h:293
TFileService::make
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
PrimaryVertexValidation::n_dzEtaBiasResiduals
TH1F * n_dzEtaBiasResiduals[nMaxBins_]
Definition: PrimaryVertexValidation.h:302
PrimaryVertexValidation::analyze
virtual void analyze(const edm::Event &, const edm::EventSetup &)
Definition: PrimaryVertexValidation.cc:123
PrimaryVertexValidation::dszFromMyVertex_
double dszFromMyVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:198
PrimaryVertexValidation::dxyFromMyVertex_
double dxyFromMyVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:189
dt_dqm_sourceclient_common_cff.reco
tuple reco
Definition: dt_dqm_sourceclient_common_cff.py:107
reco::Vertex::position
const Point & position() const
position
Definition: Vertex.h:99
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
PrimaryVertexValidation::dzErrorFromMyVertex_
double dzErrorFromMyVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:193
edm::EventAuxiliary::luminosityBlock
LuminosityBlockNumber_t luminosityBlock() const
Definition: EventAuxiliary.h:63
PrimaryVertexValidation::a_dxyPhiBiasResiduals
TH1F * a_dxyPhiBiasResiduals[nMaxBins_]
Definition: PrimaryVertexValidation.h:290
mps_fire.result
tuple result
Definition: mps_fire.py:95
PrimaryVertexValidation::IPTsigFromMyVertex_
double IPTsigFromMyVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:195
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:1244
kFPIX
const int kFPIX
Definition: AlignmentTrackSelector.cc:29
TransientVertex::normalisedChiSquared
float normalisedChiSquared() const
Definition: TransientVertex.h:126
HLT_FULL_cff.clusters
tuple clusters
Definition: HLT_FULL_cff.py:9398
PrimaryVertexValidation::tracksUsedForVertexing_
int tracksUsedForVertexing_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:187
iEvent
int iEvent
Definition: GenABIO.cc:230
PrimaryVertexValidation::hasFirstLayerPixelHits
bool hasFirstLayerPixelHits(const reco::TransientTrack track)
Definition: PrimaryVertexValidation.cc:670
VertexFitter::vertex
virtual CachingVertex< N > vertex(const std::vector< reco::TransientTrack > &tracks) const =0
PrimaryVertexValidation::n_dzPhiBiasResiduals
TH1F * n_dzPhiBiasResiduals[nMaxBins_]
Definition: PrimaryVertexValidation.h:301
TransientVertex::degreesOfFreedom
float degreesOfFreedom() const
Definition: TransientVertex.h:129
PrimaryVertexValidation::yUnbiasedVertex_
double yUnbiasedVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:181
PrimaryVertexValidation::fitResiduals
std::pair< std::pair< Double_t, Double_t >, std::pair< Double_t, Double_t > > fitResiduals(TH1 *hist)
Definition: PrimaryVertexValidation.cc:1300
TransientVertex::position
GlobalPoint position() const
Definition: TransientVertex.h:114
PV3DBase::z
T z() const
Definition: PV3DBase.h:64
Abs
T Abs(T a)
Definition: MathUtil.h:49
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:182
PrimaryVertexValidation::chi2normUnbiasedVertex_
float chi2normUnbiasedVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:183
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
TransientVertex::trackWeight
float trackWeight(const reco::TransientTrack &track) const
Definition: TransientVertex.cc:138
edm::Event::getByLabel
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:418
LogTrace
#define LogTrace(id)
Definition: MessageLogger.h:502
PrimaryVertexValidation::IPLsigFromMyVertex_
double IPLsigFromMyVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:196
reco::Vertex::ndof
double ndof() const
Definition: Vertex.h:95
edm::Event::eventAuxiliary
EventAuxiliary const & eventAuxiliary() const
Definition: Event.h:77
DetId
Definition: DetId.h:18
reco::TrackBase::qualityByName
static TrackQuality qualityByName(const std::string &name)
Definition: TrackBase.cc:125
TFileService::mkdir
TFileDirectory mkdir(const std::string &dir, const std::string &descr="")
create a new subdirectory
Definition: TFileService.h:69
TransientVertex::hasTrackWeight
bool hasTrackWeight() const
Definition: TransientVertex.h:181
PrimaryVertexValidation::isHit2D
bool isHit2D(const TrackingRecHit &hit) const
Definition: PrimaryVertexValidation.cc:640
math::XYZPoint
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
PrimaryVertexValidation::FillTrendPlot
void FillTrendPlot(TH1F *trendPlot, TH1F *residualsPlot[100], const TString &fitPar_, const TString &var_)
Definition: PrimaryVertexValidation.cc:1343
edm::EventSetup::get
const T & get() const
Definition: EventSetup.h:56
PrimaryVertexValidation::n_dxyPhiBiasResiduals
TH1F * n_dxyPhiBiasResiduals[nMaxBins_]
Definition: PrimaryVertexValidation.h:298
PrimaryVertexValidation::theTrackClusterizer_
TrackClusterizerInZ * theTrackClusterizer_
Definition: PrimaryVertexValidation.h:102
PrimaryVertexValidation::xUnbiasedVertex_
double xUnbiasedVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:180
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
HLT_FULL_cff.InputTag
tuple InputTag
Definition: HLT_FULL_cff.py:72666
PrimaryVertexValidation::a_dxyEtaBiasResiduals
TH1F * a_dxyEtaBiasResiduals[nMaxBins_]
Definition: PrimaryVertexValidation.h:291
edm::EventAuxiliary::id
EventID const & id() const
Definition: EventAuxiliary.h:59
PrimaryVertexValidation::dzFromMyVertex_
double dzFromMyVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:190
reco::Vertex::trackRef_iterator
std::vector< TrackBaseRef >::const_iterator trackRef_iterator
The iteratator for the vector&lt;TrackRef&gt;
Definition: Vertex.h:37
PrimaryVertexValidation::a_dxyBiasResidualsMap
TH1F * a_dxyBiasResidualsMap[nMaxBins_][nMaxBins_]
Definition: PrimaryVertexValidation.h:306
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:145
PrimaryVertexValidation::n_dxyResidualsMap
TH1F * n_dxyResidualsMap[nMaxBins_][nMaxBins_]
Definition: PrimaryVertexValidation.h:225
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:192
TrackingRecHit::geographicalId
DetId geographicalId() const
Definition: TrackingRecHit.h:116
HLT_FULL_cff.vertices
tuple vertices
Definition: HLT_FULL_cff.py:7118
reco::Vertex::tracks_begin
trackRef_iterator tracks_begin() const
first iterator over tracks
Definition: Vertex.cc:39
PrimaryVertexValidation::sumOfWeightsUnbiasedVertex_
float sumOfWeightsUnbiasedVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:186
PrimaryVertexValidation::n_dzBiasResidualsMap
TH1F * n_dzBiasResidualsMap[nMaxBins_][nMaxBins_]
Definition: PrimaryVertexValidation.h:310
PrimaryVertexValidation::a_dzBiasResidualsMap
TH1F * a_dzBiasResidualsMap[nMaxBins_][nMaxBins_]
Definition: PrimaryVertexValidation.h:307
PV3DBase::x
T x() const
Definition: PV3DBase.h:62
PrimaryVertexValidation::a_dzEtaBiasResiduals
TH1F * a_dzEtaBiasResiduals[nMaxBins_]
Definition: PrimaryVertexValidation.h:294
TransientVertex::isValid
bool isValid() const
Definition: TransientVertex.h:134
PrimaryVertexValidation::a_dzResidualsMap
TH1F * a_dzResidualsMap[nMaxBins_][nMaxBins_]
Definition: PrimaryVertexValidation.h:223
PrimaryVertexValidation::chi2UnbiasedVertex_
float chi2UnbiasedVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:184
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
trackingRecHit_iterator
TrackingRecHitCollection::base::const_iterator trackingRecHit_iterator
iterator over a vector of reference to TrackingRecHit in the same collection
Definition: TrackingRecHitFwd.h:19
GlobalVector
Global3DVector GlobalVector
Definition: GlobalVector.h:10
kBPIX
const int kBPIX
Definition: AlignmentTrackSelector.cc:28
Generated for CMSSW Reference Manual by doxygen 1.8.5