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  bool accepted = false;
673  // hit pattern of the track
674  const reco::HitPattern& p = track.hitPattern();
675  for (int i=0; i<p.numberOfHits(); i++) {
676  uint32_t pattern = p.getHitPattern(i);
677  if (p.pixelBarrelHitFilter(pattern) || p.pixelEndcapHitFilter(pattern) ) {
678  if (p.getLayer(pattern) == 1) {
679  if (p.validHitFilter(pattern)) {
680  accepted = true;
681  }
682  }
683  }
684  }
685  return accepted;
686 }
687 
688 // ------------ method called once each job before begining the event loop ------------
689 void PrimaryVertexValidation::beginJob()
690 {
691  edm::LogInfo("beginJob") << "Begin Job" << std::endl;
692  // Define TTree for output
693  Nevt_ = 0;
694 
695  // rootFile_ = new TFile(filename_.c_str(),"recreate");
696  edm::Service<TFileService> fs;
697  rootTree_ = fs->make<TTree>("tree","PV Validation tree");
698 
699  // Track Paramters
700 
701  if(lightNtupleSwitch_){
702 
703  rootTree_->Branch("EventNumber",&EventNumber_,"EventNumber/i");
704  rootTree_->Branch("RunNumber",&RunNumber_,"RunNumber/i");
705  rootTree_->Branch("LuminosityBlockNumber",&LuminosityBlockNumber_,"LuminosityBlockNumber/i");
706  rootTree_->Branch("nOfflineVertices",&nOfflineVertices_,"nOfflineVertices/I");
707  rootTree_->Branch("nTracks",&nTracks_,"nTracks/I");
708  rootTree_->Branch("phi",&phi_,"phi[nTracks]/D");
709  rootTree_->Branch("eta",&eta_,"eta[nTracks]/D");
710  rootTree_->Branch("pt",&pt_,"pt[nTracks]/D");
711  rootTree_->Branch("dxyFromMyVertex",&dxyFromMyVertex_,"dxyFromMyVertex[nTracks]/D");
712  rootTree_->Branch("dzFromMyVertex",&dzFromMyVertex_,"dzFromMyVertex[nTracks]/D");
713  rootTree_->Branch("IPTsigFromMyVertex",&IPTsigFromMyVertex_,"IPTsigFromMyVertex_[nTracks]/D");
714  rootTree_->Branch("IPLsigFromMyVertex",&IPLsigFromMyVertex_,"IPLsigFromMyVertex_[nTracks]/D");
715  rootTree_->Branch("hasRecVertex",&hasRecVertex_,"hasRecVertex[nTracks]/I");
716  rootTree_->Branch("isGoodTrack",&isGoodTrack_,"isGoodTrack[nTracks]/I");
717  rootTree_->Branch("isHighPurity",&isHighPurity_,"isHighPurity_[nTracks]/I");
718 
719  } else {
720 
721  rootTree_->Branch("nTracks",&nTracks_,"nTracks/I");
722  rootTree_->Branch("nTracksPerClus",&nTracksPerClus_,"nTracksPerClus/I");
723  rootTree_->Branch("nClus",&nClus_,"nClus/I");
724  rootTree_->Branch("xOfflineVertex",&xOfflineVertex_,"xOfflineVertex/D");
725  rootTree_->Branch("yOfflineVertex",&yOfflineVertex_,"yOfflineVertex/D");
726  rootTree_->Branch("zOfflineVertex",&zOfflineVertex_,"zOfflineVertex/D");
727  rootTree_->Branch("BSx0",&BSx0_,"BSx0/D");
728  rootTree_->Branch("BSy0",&BSy0_,"BSy0/D");
729  rootTree_->Branch("BSz0",&BSz0_,"BSz0/D");
730  rootTree_->Branch("Beamsigmaz",&Beamsigmaz_,"Beamsigmaz/D");
731  rootTree_->Branch("Beamdxdz",&Beamdxdz_,"Beamdxdz/D");
732  rootTree_->Branch("BeamWidthX",&BeamWidthX_,"BeamWidthX/D");
733  rootTree_->Branch("BeamWidthY",&BeamWidthY_,"BeamWidthY/D");
734  rootTree_->Branch("pt",&pt_,"pt[nTracks]/D");
735  rootTree_->Branch("p",&p_,"p[nTracks]/D");
736  rootTree_->Branch("nhits",&nhits_,"nhits[nTracks]/I");
737  rootTree_->Branch("nhits1D",&nhits1D_,"nhits1D[nTracks]/I");
738  rootTree_->Branch("nhits2D",&nhits2D_,"nhits2D[nTracks]/I");
739  rootTree_->Branch("nhitsBPIX",&nhitsBPIX_,"nhitsBPIX[nTracks]/I");
740  rootTree_->Branch("nhitsFPIX",&nhitsFPIX_,"nhitsFPIX[nTracks]/I");
741  rootTree_->Branch("nhitsTIB",&nhitsTIB_,"nhitsTIB[nTracks]/I");
742  rootTree_->Branch("nhitsTID",&nhitsTID_,"nhitsTID[nTracks]/I");
743  rootTree_->Branch("nhitsTOB",&nhitsTOB_,"nhitsTOB[nTracks]/I");
744  rootTree_->Branch("nhitsTEC",&nhitsTEC_,"nhitsTEC[nTracks]/I");
745  rootTree_->Branch("eta",&eta_,"eta[nTracks]/D");
746  rootTree_->Branch("theta",&theta_,"theta[nTracks]/D");
747  rootTree_->Branch("phi",&phi_,"phi[nTracks]/D");
748  rootTree_->Branch("chi2",&chi2_,"chi2[nTracks]/D");
749  rootTree_->Branch("chi2ndof",&chi2ndof_,"chi2ndof[nTracks]/D");
750  rootTree_->Branch("charge",&charge_,"charge[nTracks]/I");
751  rootTree_->Branch("qoverp",&qoverp_,"qoverp[nTracks]/D");
752  rootTree_->Branch("dz",&dz_,"dz[nTracks]/D");
753  rootTree_->Branch("dxy",&dxy_,"dxy[nTracks]/D");
754  rootTree_->Branch("dzBs",&dzBs_,"dzBs[nTracks]/D");
755  rootTree_->Branch("dxyBs",&dxyBs_,"dxyBs[nTracks]/D");
756  rootTree_->Branch("xPCA",&xPCA_,"xPCA[nTracks]/D");
757  rootTree_->Branch("yPCA",&yPCA_,"yPCA[nTracks]/D");
758  rootTree_->Branch("zPCA",&zPCA_,"zPCA[nTracks]/D");
759  rootTree_->Branch("xUnbiasedVertex",&xUnbiasedVertex_,"xUnbiasedVertex[nTracks]/D");
760  rootTree_->Branch("yUnbiasedVertex",&yUnbiasedVertex_,"yUnbiasedVertex[nTracks]/D");
761  rootTree_->Branch("zUnbiasedVertex",&zUnbiasedVertex_,"zUnbiasedVertex[nTracks]/D");
762  rootTree_->Branch("chi2normUnbiasedVertex",&chi2normUnbiasedVertex_,"chi2normUnbiasedVertex[nTracks]/F");
763  rootTree_->Branch("chi2UnbiasedVertex",&chi2UnbiasedVertex_,"chi2UnbiasedVertex[nTracks]/F");
764  rootTree_->Branch("DOFUnbiasedVertex",&DOFUnbiasedVertex_," DOFUnbiasedVertex[nTracks]/F");
765  rootTree_->Branch("sumOfWeightsUnbiasedVertex",&sumOfWeightsUnbiasedVertex_,"sumOfWeightsUnbiasedVertex[nTracks]/F");
766  rootTree_->Branch("tracksUsedForVertexing",&tracksUsedForVertexing_,"tracksUsedForVertexing[nTracks]/I");
767  rootTree_->Branch("dxyFromMyVertex",&dxyFromMyVertex_,"dxyFromMyVertex[nTracks]/D");
768  rootTree_->Branch("dzFromMyVertex",&dzFromMyVertex_,"dzFromMyVertex[nTracks]/D");
769  rootTree_->Branch("dszFromMyVertex",&dszFromMyVertex_,"dszFromMyVertex[nTracks]/D");
770  rootTree_->Branch("dxyErrorFromMyVertex",&dxyErrorFromMyVertex_,"dxyErrorFromMyVertex_[nTracks]/D");
771  rootTree_->Branch("dzErrorFromMyVertex",&dzErrorFromMyVertex_,"dzErrorFromMyVertex_[nTracks]/D");
772  rootTree_->Branch("IPTsigFromMyVertex",&IPTsigFromMyVertex_,"IPTsigFromMyVertex_[nTracks]/D");
773  rootTree_->Branch("IPLsigFromMyVertex",&IPLsigFromMyVertex_,"IPLsigFromMyVertex_[nTracks]/D");
774  rootTree_->Branch("hasRecVertex",&hasRecVertex_,"hasRecVertex[nTracks]/I");
775  rootTree_->Branch("isGoodTrack",&isGoodTrack_,"isGoodTrack[nTracks]/I");
776  }
777 
778  // probe track histograms
779 
780  TFileDirectory ProbeFeatures = fs->mkdir("ProbeTrackFeatures");
781 
782  h_probePt_ = ProbeFeatures.make<TH1F>("h_probePt","p_{T} of probe track;track p_{T} (GeV); tracks",100,0.,50.);
783  h_probeEta_ = ProbeFeatures.make<TH1F>("h_probeEta","#eta of probe track;track #eta; tracks",54,-2.7,2.7);
784  h_probePhi_ = ProbeFeatures.make<TH1F>("h_probePhi","#phi of probe track;track #phi [rad]; tracks",100,-3.15,3.15);
785  h_probeChi2_ = ProbeFeatures.make<TH1F>("h_probeChi2","#chi^{2} of probe track;track #chi^{2}; tracks",100,0.,100.);
786  h_probeNormChi2_ = ProbeFeatures.make<TH1F>("h_probeNormChi2"," normalized #chi^{2} of probe track;track #chi^{2}/ndof; tracks",100,0.,10.);
787  h_probeCharge_ = ProbeFeatures.make<TH1F>("h_probeCharge","charge of profe track;track charge Q;tracks",3,-1.5,1.5);
788  h_probeQoverP_ = ProbeFeatures.make<TH1F>("h_probeQoverP","q/p of probe track; track Q/p (GeV^{-1})",200,-1.,1.);
789  h_probedz_ = ProbeFeatures.make<TH1F>("h_probedz","d_{z} of probe track;track d_{z} (cm);tracks",100,-5.,5.);
790  h_probedxy_ = ProbeFeatures.make<TH1F>("h_probedxy","d_{xy} of probe track;track d_{xy} (#mum);tracks",200,-1.,1.);
791 
792  h_probeHits_ = ProbeFeatures.make<TH1F>("h_probeNRechits" ,"N_{hits} ;N_{hits} ;tracks",40,-0.5,39.5);
793  h_probeHits1D_ = ProbeFeatures.make<TH1F>("h_probeNRechits1D" ,"N_{hits} 1D ;N_{hits} 1D ;tracks",40,-0.5,39.5);
794  h_probeHits2D_ = ProbeFeatures.make<TH1F>("h_probeNRechits2D" ,"N_{hits} 2D ;N_{hits} 2D ;tracks",40,-0.5,39.5);
795  h_probeHitsInTIB_ = ProbeFeatures.make<TH1F>("h_probeNRechitsTIB" ,"N_{hits} TIB ;N_{hits} TIB;tracks",40,-0.5,39.5);
796  h_probeHitsInTOB_ = ProbeFeatures.make<TH1F>("h_probeNRechitsTOB" ,"N_{hits} TOB ;N_{hits} TOB;tracks",40,-0.5,39.5);
797  h_probeHitsInTID_ = ProbeFeatures.make<TH1F>("h_probeNRechitsTID" ,"N_{hits} TID ;N_{hits} TID;tracks",40,-0.5,39.5);
798  h_probeHitsInTEC_ = ProbeFeatures.make<TH1F>("h_probeNRechitsTEC" ,"N_{hits} TEC ;N_{hits} TEC;tracks",40,-0.5,39.5);
799  h_probeHitsInBPIX_ = ProbeFeatures.make<TH1F>("h_probeNRechitsBPIX","N_{hits} BPIX;N_{hits} BPIX;tracks",40,-0.5,39.5);
800  h_probeHitsInFPIX_ = ProbeFeatures.make<TH1F>("h_probeNRechitsFPIX","N_{hits} FPIX;N_{hits} FPIX;tracks",40,-0.5,39.5);
801 
802  TFileDirectory RefitVertexFeatures = fs->mkdir("RefitVertexFeatures");
803  h_fitVtxNtracks_ = RefitVertexFeatures.make<TH1F>("h_fitVtxNtracks" ,"N^{vtx}_{trks};N^{vtx}_{trks};vertices" ,100,-0.5,99.5);
804  h_fitVtxNdof_ = RefitVertexFeatures.make<TH1F>("h_fitVtxNdof" ,"N^{vtx}_{DOF};N^{vtx}_{DOF};vertices" ,100,-0.5,99.5);
805  h_fitVtxChi2_ = RefitVertexFeatures.make<TH1F>("h_fitVtxChi2" ,"#chi^{2} vtx;#chi^{2} vtx;vertices" ,100,-0.5,99.5);
806  h_fitVtxChi2ndf_ = RefitVertexFeatures.make<TH1F>("h_fitVtxChi2ndf" ,"#chi^{2}/ndf vtx;#chi^{2}/ndf vtx;vertices" ,100,-0.5,9.5);
807  h_fitVtxChi2Prob_ = RefitVertexFeatures.make<TH1F>("h_fitVtxChi2Prob" ,"Prob(#chi^{2},ndf);Prob(#chi^{2},ndf);vertices",40,0.,1.);
808 
809  if(useTracksFromRecoVtx_) {
810 
811  TFileDirectory RecoVertexFeatures = fs->mkdir("RecoVertexFeatures");
812  h_recoVtxNtracks_ = RecoVertexFeatures.make<TH1F>("h_recoVtxNtracks" ,"N^{vtx}_{trks};N^{vtx}_{trks};vertices" ,100,-0.5,99.5);
813  h_recoVtxChi2ndf_ = RecoVertexFeatures.make<TH1F>("h_recoVtxChi2ndf" ,"#chi^{2}/ndf vtx;#chi^{2}/ndf vtx;vertices" ,10,-0.5,9.5);
814  h_recoVtxChi2Prob_ = RecoVertexFeatures.make<TH1F>("h_recoVtxChi2Prob" ,"Prob(#chi^{2},ndf);Prob(#chi^{2},ndf);vertices",40,0.,1.);
815  h_recoVtxSumPt_ = RecoVertexFeatures.make<TH1F>("h_recoVtxSumPt" ,"Sum(p^{trks}_{T});Sum(p^{trks}_{T});vertices" ,100,0.,200.);
816 
817  }
818 
819  // initialize the residuals histograms
820 
821  float dxymax_phi = 2000;
822  float dzmax_phi = 2000;
823  float dxymax_eta = 3000;
824  float dzmax_eta = 3000;
825 
826  const Int_t mybins_ = 500;
827 
829  //
830  // Usual plots from refitting the vertex
831  // The vertex is refit without the probe track
832  //
834 
835  TFileDirectory AbsTransPhiRes = fs->mkdir("Abs_Transv_Phi_Residuals");
836  TFileDirectory AbsTransEtaRes = fs->mkdir("Abs_Transv_Eta_Residuals");
837 
838  TFileDirectory AbsLongPhiRes = fs->mkdir("Abs_Long_Phi_Residuals");
839  TFileDirectory AbsLongEtaRes = fs->mkdir("Abs_Long_Eta_Residuals");
840 
841  TFileDirectory NormTransPhiRes = fs->mkdir("Norm_Transv_Phi_Residuals");
842  TFileDirectory NormTransEtaRes = fs->mkdir("Norm_Transv_Eta_Residuals");
843 
844  TFileDirectory NormLongPhiRes = fs->mkdir("Norm_Long_Phi_Residuals");
845  TFileDirectory NormLongEtaRes = fs->mkdir("Norm_Long_Eta_Residuals");
846 
847  TFileDirectory AbsDoubleDiffRes = fs->mkdir("Abs_DoubleDiffResiduals");
848  TFileDirectory NormDoubleDiffRes = fs->mkdir("Norm_DoubleDiffResiduals");
849 
850  for ( int i=0; i<nBins_; ++i ) {
851 
852  float phiF = (-TMath::Pi()+i*phipitch_)*(180/TMath::Pi());
853  float phiL = (-TMath::Pi()+(i+1)*phipitch_)*(180/TMath::Pi());
854 
855  float etaF=-2.5+i*etapitch_;
856  float etaL=-2.5+(i+1)*etapitch_;
857 
858  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);
859  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);
860 
861  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);
862  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);
863 
864  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.);
865  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.);
866 
867  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.);
868  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.);
869 
870  for ( int j=0; j<nBins_; ++j ) {
871 
872  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);
873  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);
874 
875  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);
876  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);
877  }
878 
879  }
880 
881  // declaration of the directories
882 
883  TFileDirectory MeanTrendsDir = fs->mkdir("MeanTrends");
884  TFileDirectory WidthTrendsDir = fs->mkdir("WidthTrends");
885  TFileDirectory MedianTrendsDir = fs->mkdir("MedianTrends");
886  TFileDirectory MADTrendsDir = fs->mkdir("MADTrends");
887 
888  TFileDirectory Mean2DMapsDir = fs->mkdir("MeanMaps");
889  TFileDirectory Width2DMapsDir = fs->mkdir("WidthMaps");
890 
891  Double_t highedge=nBins_-0.5;
892  Double_t lowedge=-0.5;
893 
894  // means and widths from the fit
895 
896  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);
897  a_dxyPhiWidthTrend = WidthTrendsDir.make<TH1F>("widths_dxy_phi","#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees];#sigma_{d_{xy}} [#mum]",nBins_,lowedge,highedge);
898  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);
899  a_dzPhiWidthTrend = WidthTrendsDir.make<TH1F>("widths_dz_phi","#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];#sigma_{d_{z}} [#mum]",nBins_,lowedge,highedge);
900 
901  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);
902  a_dxyEtaWidthTrend = WidthTrendsDir.make<TH1F>("widths_dxy_eta","#sigma_{d_{xy}} vs #eta sector;#eta (sector);#sigma_{d_{xy}} [#mum]",nBins_,lowedge,highedge);
903  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);
904  a_dzEtaWidthTrend = WidthTrendsDir.make<TH1F>("widths_dz_eta","#sigma_{d_{xy}} vs #eta sector;#eta (sector);#sigma_{d_{z}} [#mum]",nBins_,lowedge,highedge);
905 
906  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);
907  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);
908  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);
909  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);
910 
911  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);
912  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);
913  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);
914  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);
915 
916  // 2D maps
917 
918  a_dxyMeanMap = Mean2DMapsDir.make<TH2F> ("means_dxy_map","#LT d_{xy} #GT map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
919  a_dzMeanMap = Mean2DMapsDir.make<TH2F> ("means_dz_map","#LT d_{z} #GT map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
920 
921  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);
922  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);
923 
924  a_dxyWidthMap = Width2DMapsDir.make<TH2F> ("widths_dxy_map","#sigma_{d_{xy}} map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
925  a_dzWidthMap = Width2DMapsDir.make<TH2F> ("widths_dz_map","#sigma_{d_{z}} map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
926 
927  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);
928  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);
929 
930  // medians and MADs
931 
932  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);
933  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);
934  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);
935  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);
936 
937  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);
938  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);
939  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);
940  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);
941 
942  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);
943  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);
944  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);
945  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);
946 
947  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);
948  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);
949  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);
950  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);
951 
952 
954  //
955  // plots of biased residuals
956  // The vertex is refit without the probe track
957  //
959 
960  if (useTracksFromRecoVtx_){
961 
962  TFileDirectory AbsTransPhiBiasRes = fs->mkdir("Abs_Transv_Phi_BiasResiduals");
963  TFileDirectory AbsTransEtaBiasRes = fs->mkdir("Abs_Transv_Eta_BiasResiduals");
964 
965  TFileDirectory AbsLongPhiBiasRes = fs->mkdir("Abs_Long_Phi_BiasResiduals");
966  TFileDirectory AbsLongEtaBiasRes = fs->mkdir("Abs_Long_Eta_BiasResiduals");
967 
968  TFileDirectory NormTransPhiBiasRes = fs->mkdir("Norm_Transv_Phi_BiasResiduals");
969  TFileDirectory NormTransEtaBiasRes = fs->mkdir("Norm_Transv_Eta_BiasResiduals");
970 
971  TFileDirectory NormLongPhiBiasRes = fs->mkdir("Norm_Long_Phi_BiasResiduals");
972  TFileDirectory NormLongEtaBiasRes = fs->mkdir("Norm_Long_Eta_BiasResiduals");
973 
974  TFileDirectory AbsDoubleDiffBiasRes = fs->mkdir("Abs_DoubleDiffBiasResiduals");
975  TFileDirectory NormDoubleDiffBiasRes = fs->mkdir("Norm_DoubleDiffBiasResiduals");
976 
977  for ( int i=0; i<nBins_; ++i ) {
978 
979  float phiF = (-TMath::Pi()+i*phipitch_)*(180/TMath::Pi());
980  float phiL = (-TMath::Pi()+(i+1)*phipitch_)*(180/TMath::Pi());
981 
982  float etaF=-2.5+i*etapitch_;
983  float etaL=-2.5+(i+1)*etapitch_;
984 
985  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);
986  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);
987 
988  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);
989  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);
990 
991  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.);
992  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.);
993 
994  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.);
995  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.);
996 
997  for ( int j=0; j<nBins_; ++j ) {
998 
999  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);
1000  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);
1001 
1002  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);
1003  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);
1004  }
1005 
1006  }
1007 
1008  // declaration of the directories
1009 
1010  TFileDirectory MeanBiasTrendsDir = fs->mkdir("MeanBiasTrends");
1011  TFileDirectory WidthBiasTrendsDir = fs->mkdir("WidthBiasTrends");
1012  TFileDirectory MedianBiasTrendsDir = fs->mkdir("MedianBiasTrends");
1013  TFileDirectory MADBiasTrendsDir = fs->mkdir("MADBiasTrends");
1014 
1015  TFileDirectory Mean2DBiasMapsDir = fs->mkdir("MeanBiasMaps");
1016  TFileDirectory Width2DBiasMapsDir = fs->mkdir("WidthBiasMaps");
1017 
1018  // means and widths from the fit
1019 
1020  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);
1021  a_dxyPhiWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>("widths_dxy_phi","#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees];#sigma_{d_{xy}} [#mum]",nBins_,lowedge,highedge);
1022  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);
1023  a_dzPhiWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>("widths_dz_phi","#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];#sigma_{d_{z}} [#mum]",nBins_,lowedge,highedge);
1024 
1025  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);
1026  a_dxyEtaWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>("widths_dxy_eta","#sigma_{d_{xy}} vs #eta sector;#eta (sector);#sigma_{d_{xy}} [#mum]",nBins_,lowedge,highedge);
1027  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);
1028  a_dzEtaWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>("widths_dz_eta","#sigma_{d_{xy}} vs #eta sector;#eta (sector);#sigma_{d_{z}} [#mum]",nBins_,lowedge,highedge);
1029 
1030  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);
1031  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);
1032  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);
1033  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);
1034 
1035  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);
1036  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);
1037  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);
1038  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);
1039 
1040  // 2D maps
1041 
1042  a_dxyMeanBiasMap = Mean2DBiasMapsDir.make<TH2F> ("means_dxy_map","#LT d_{xy} #GT map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1043  a_dzMeanBiasMap = Mean2DBiasMapsDir.make<TH2F> ("means_dz_map","#LT d_{z} #GT map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1044 
1045  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);
1046  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);
1047 
1048  a_dxyWidthBiasMap = Width2DBiasMapsDir.make<TH2F> ("widths_dxy_map","#sigma_{d_{xy}} map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1049  a_dzWidthBiasMap = Width2DBiasMapsDir.make<TH2F> ("widths_dz_map","#sigma_{d_{z}} map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
1050 
1051  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);
1052  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);
1053 
1054  // medians and MADs
1055 
1056  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);
1057  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);
1058  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);
1059  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);
1060 
1061  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);
1062  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);
1063  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);
1064  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);
1065 
1066  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);
1067  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);
1068  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);
1069  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);
1070 
1071  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);
1072  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);
1073  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);
1074  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);
1075 
1076  }
1077 }
1078 // ------------ method called once each job just after ending the event loop ------------
1079 void PrimaryVertexValidation::endJob()
1080 {
1081 
1082  std::cout<<"######################################"<<std::endl;
1083  std::cout<<"Number of analyzed events: "<<Nevt_<<std::endl;
1084  std::cout<<"######################################"<<std::endl;
1085 
1086  if(useTracksFromRecoVtx_){
1087 
1088  FillTrendPlot(a_dxyPhiMeanBiasTrend ,a_dxyPhiBiasResiduals,"mean","phi");
1089  FillTrendPlot(a_dxyPhiWidthBiasTrend,a_dxyPhiBiasResiduals,"width","phi");
1090  FillTrendPlot(a_dzPhiMeanBiasTrend ,a_dzPhiBiasResiduals ,"mean","phi");
1091  FillTrendPlot(a_dzPhiWidthBiasTrend ,a_dzPhiBiasResiduals ,"width","phi");
1092 
1093  FillTrendPlot(a_dxyEtaMeanBiasTrend ,a_dxyEtaBiasResiduals,"mean","eta");
1094  FillTrendPlot(a_dxyEtaWidthBiasTrend,a_dxyEtaBiasResiduals,"width","eta");
1095  FillTrendPlot(a_dzEtaMeanBiasTrend ,a_dzEtaBiasResiduals ,"mean","eta");
1096  FillTrendPlot(a_dzEtaWidthBiasTrend ,a_dzEtaBiasResiduals ,"width","eta");
1097 
1098  FillTrendPlot(n_dxyPhiMeanBiasTrend ,n_dxyPhiBiasResiduals,"mean","phi");
1099  FillTrendPlot(n_dxyPhiWidthBiasTrend,n_dxyPhiBiasResiduals,"width","phi");
1100  FillTrendPlot(n_dzPhiMeanBiasTrend ,n_dzPhiBiasResiduals ,"mean","phi");
1101  FillTrendPlot(n_dzPhiWidthBiasTrend ,n_dzPhiBiasResiduals ,"width","phi");
1102 
1103  FillTrendPlot(n_dxyEtaMeanBiasTrend ,n_dxyEtaBiasResiduals,"mean","eta");
1104  FillTrendPlot(n_dxyEtaWidthBiasTrend,n_dxyEtaBiasResiduals,"width","eta");
1105  FillTrendPlot(n_dzEtaMeanBiasTrend ,n_dzEtaBiasResiduals ,"mean","eta");
1106  FillTrendPlot(n_dzEtaWidthBiasTrend ,n_dzEtaBiasResiduals ,"width","eta");
1107 
1108  // medians and MADs
1109 
1110  FillTrendPlot(a_dxyPhiMedianBiasTrend,a_dxyPhiBiasResiduals,"median","phi");
1111  FillTrendPlot(a_dxyPhiMADBiasTrend ,a_dxyPhiBiasResiduals,"mad","phi");
1112  FillTrendPlot(a_dzPhiMedianBiasTrend ,a_dzPhiBiasResiduals ,"median","phi");
1113  FillTrendPlot(a_dzPhiMADBiasTrend ,a_dzPhiBiasResiduals ,"mad","phi");
1114 
1115  FillTrendPlot(a_dxyEtaMedianBiasTrend,a_dxyEtaBiasResiduals,"median","eta");
1116  FillTrendPlot(a_dxyEtaMADBiasTrend ,a_dxyEtaBiasResiduals,"mad","eta");
1117  FillTrendPlot(a_dzEtaMedianBiasTrend ,a_dzEtaBiasResiduals ,"median","eta");
1118  FillTrendPlot(a_dzEtaMADBiasTrend ,a_dzEtaBiasResiduals ,"mad","eta");
1119 
1120  FillTrendPlot(n_dxyPhiMedianBiasTrend,n_dxyPhiBiasResiduals,"median","phi");
1121  FillTrendPlot(n_dxyPhiMADBiasTrend ,n_dxyPhiBiasResiduals,"mad","phi");
1122  FillTrendPlot(n_dzPhiMedianBiasTrend ,n_dzPhiBiasResiduals ,"median","phi");
1123  FillTrendPlot(n_dzPhiMADBiasTrend ,n_dzPhiBiasResiduals ,"mad","phi");
1124 
1125  FillTrendPlot(n_dxyEtaMedianBiasTrend,n_dxyEtaBiasResiduals,"median","eta");
1126  FillTrendPlot(n_dxyEtaMADBiasTrend ,n_dxyEtaBiasResiduals,"mad","eta");
1127  FillTrendPlot(n_dzEtaMedianBiasTrend ,n_dzEtaBiasResiduals ,"median","eta");
1128  FillTrendPlot(n_dzEtaMADBiasTrend ,n_dzEtaBiasResiduals ,"mad","eta");
1129 
1130  }
1131 
1132  FillTrendPlot(a_dxyPhiMeanTrend ,a_dxyPhiResiduals,"mean","phi");
1133  FillTrendPlot(a_dxyPhiWidthTrend,a_dxyPhiResiduals,"width","phi");
1134  FillTrendPlot(a_dzPhiMeanTrend ,a_dzPhiResiduals ,"mean","phi");
1135  FillTrendPlot(a_dzPhiWidthTrend ,a_dzPhiResiduals ,"width","phi");
1136 
1137  FillTrendPlot(a_dxyEtaMeanTrend ,a_dxyEtaResiduals,"mean","eta");
1138  FillTrendPlot(a_dxyEtaWidthTrend,a_dxyEtaResiduals,"width","eta");
1139  FillTrendPlot(a_dzEtaMeanTrend ,a_dzEtaResiduals ,"mean","eta");
1140  FillTrendPlot(a_dzEtaWidthTrend ,a_dzEtaResiduals ,"width","eta");
1141 
1142  FillTrendPlot(n_dxyPhiMeanTrend ,n_dxyPhiResiduals,"mean","phi");
1143  FillTrendPlot(n_dxyPhiWidthTrend,n_dxyPhiResiduals,"width","phi");
1144  FillTrendPlot(n_dzPhiMeanTrend ,n_dzPhiResiduals ,"mean","phi");
1145  FillTrendPlot(n_dzPhiWidthTrend ,n_dzPhiResiduals ,"width","phi");
1146 
1147  FillTrendPlot(n_dxyEtaMeanTrend ,n_dxyEtaResiduals,"mean","eta");
1148  FillTrendPlot(n_dxyEtaWidthTrend,n_dxyEtaResiduals,"width","eta");
1149  FillTrendPlot(n_dzEtaMeanTrend ,n_dzEtaResiduals ,"mean","eta");
1150  FillTrendPlot(n_dzEtaWidthTrend ,n_dzEtaResiduals ,"width","eta");
1151 
1152  // medians and MADs
1153 
1154  FillTrendPlot(a_dxyPhiMedianTrend,a_dxyPhiResiduals,"median","phi");
1155  FillTrendPlot(a_dxyPhiMADTrend ,a_dxyPhiResiduals,"mad","phi");
1156  FillTrendPlot(a_dzPhiMedianTrend ,a_dzPhiResiduals ,"median","phi");
1157  FillTrendPlot(a_dzPhiMADTrend ,a_dzPhiResiduals ,"mad","phi");
1158 
1159  FillTrendPlot(a_dxyEtaMedianTrend,a_dxyEtaResiduals,"median","eta");
1160  FillTrendPlot(a_dxyEtaMADTrend ,a_dxyEtaResiduals,"mad","eta");
1161  FillTrendPlot(a_dzEtaMedianTrend ,a_dzEtaResiduals ,"median","eta");
1162  FillTrendPlot(a_dzEtaMADTrend ,a_dzEtaResiduals ,"mad","eta");
1163 
1164  FillTrendPlot(n_dxyPhiMedianTrend,n_dxyPhiResiduals,"median","phi");
1165  FillTrendPlot(n_dxyPhiMADTrend ,n_dxyPhiResiduals,"mad","phi");
1166  FillTrendPlot(n_dzPhiMedianTrend ,n_dzPhiResiduals ,"median","phi");
1167  FillTrendPlot(n_dzPhiMADTrend ,n_dzPhiResiduals ,"mad","phi");
1168 
1169  FillTrendPlot(n_dxyEtaMedianTrend,n_dxyEtaResiduals,"median","eta");
1170  FillTrendPlot(n_dxyEtaMADTrend ,n_dxyEtaResiduals,"mad","eta");
1171  FillTrendPlot(n_dzEtaMedianTrend ,n_dzEtaResiduals ,"median","eta");
1172  FillTrendPlot(n_dzEtaMADTrend ,n_dzEtaResiduals ,"mad","eta");
1173 
1174 }
1175 
1176 //*************************************************************
1177 void PrimaryVertexValidation::SetVarToZero()
1178 //*************************************************************
1179 {
1180  nTracks_ = 0;
1181  nClus_ = 0;
1182  nOfflineVertices_=0;
1183  RunNumber_ =0;
1184  LuminosityBlockNumber_=0;
1185  xOfflineVertex_ =-999.;
1186  yOfflineVertex_ =-999.;
1187  zOfflineVertex_ =-999.;
1188  BSx0_ = -999.;
1189  BSy0_ = -999.;
1190  BSz0_ = -999.;
1191  Beamsigmaz_=-999.;
1192  Beamdxdz_=-999.;
1193  BeamWidthX_=-999.;
1194  BeamWidthY_=-999.;
1195 
1196  for ( int i=0; i<nMaxtracks_; ++i ) {
1197 
1198  pt_[i] = 0;
1199  p_[i] = 0;
1200  nhits_[i] = 0;
1201  nhits1D_[i] = 0;
1202  nhits2D_[i] = 0;
1203  nhitsBPIX_[i] = 0;
1204  nhitsFPIX_[i] = 0;
1205  nhitsTIB_[i] = 0;
1206  nhitsTID_[i] = 0;
1207  nhitsTOB_[i] = 0;
1208  nhitsTEC_[i] = 0;
1209  isHighPurity_[i] = 0;
1210  eta_[i] = 0;
1211  theta_[i] = 0;
1212  phi_[i] = 0;
1213  chi2_[i] = 0;
1214  chi2ndof_[i] = 0;
1215  charge_[i] = 0;
1216  qoverp_[i] = 0;
1217  dz_[i] = 0;
1218  dxy_[i] = 0;
1219  dzBs_[i] = 0;
1220  dxyBs_[i] = 0;
1221  xPCA_[i] = 0;
1222  yPCA_[i] = 0;
1223  zPCA_[i] = 0;
1224  xUnbiasedVertex_[i] =0;
1225  yUnbiasedVertex_[i] =0;
1226  zUnbiasedVertex_[i] =0;
1227  chi2normUnbiasedVertex_[i]=0;
1228  chi2UnbiasedVertex_[i]=0;
1229  DOFUnbiasedVertex_[i]=0;
1230  sumOfWeightsUnbiasedVertex_[i]=0;
1231  tracksUsedForVertexing_[i]=0;
1232  dxyFromMyVertex_[i]=0;
1233  dzFromMyVertex_[i]=0;
1234  dszFromMyVertex_[i]=0;
1235  dxyErrorFromMyVertex_[i]=0;
1236  dzErrorFromMyVertex_[i]=0;
1237  IPTsigFromMyVertex_[i]=0;
1238  IPLsigFromMyVertex_[i]=0;
1239  hasRecVertex_[i] = 0;
1240  isGoodTrack_[i] = 0;
1241  }
1242 }
1243 
1244 //*************************************************************
1245 std::pair<Double_t,Double_t> PrimaryVertexValidation::getMedian(TH1F *histo)
1246 //*************************************************************
1247 {
1248  Double_t median = 999;
1249  int nbins = histo->GetNbinsX();
1250 
1251  //extract median from histogram
1252  double *x = new double[nbins];
1253  double *y = new double[nbins];
1254  for (int j = 0; j < nbins; j++) {
1255  x[j] = histo->GetBinCenter(j+1);
1256  y[j] = histo->GetBinContent(j+1);
1257  }
1258  median = TMath::Median(nbins, x, y);
1259 
1260  delete[] x; x = 0;
1261  delete [] y; y = 0;
1262 
1263  std::pair<Double_t,Double_t> result;
1264  result = std::make_pair(median,median/TMath::Sqrt(histo->GetEntries()));
1265 
1266  return result;
1267 
1268 }
1269 
1270 //*************************************************************
1271 std::pair<Double_t,Double_t> PrimaryVertexValidation::getMAD(TH1F *histo)
1272 //*************************************************************
1273 {
1274 
1275  int nbins = histo->GetNbinsX();
1276  Double_t median = getMedian(histo).first;
1277  Double_t x_lastBin = histo->GetBinLowEdge(nbins+1);
1278  const char *HistoName =histo->GetName();
1279  TString Finalname = Form("resMed%s",HistoName);
1280  TH1F *newHisto = new TH1F(Finalname,Finalname,nbins,0.,x_lastBin);
1281  Double_t *residuals = new Double_t[nbins];
1282  Double_t *weights = new Double_t[nbins];
1283 
1284  for (int j = 0; j < nbins; j++) {
1285  residuals[j] = TMath::Abs(median - histo->GetBinCenter(j+1));
1286  weights[j]=histo->GetBinContent(j+1);
1287  newHisto->Fill(residuals[j],weights[j]);
1288  }
1289 
1290  Double_t theMAD = (getMedian(newHisto).first)*1.4826;
1291  newHisto->Delete("");
1292 
1293  std::pair<Double_t,Double_t> result;
1294  result = std::make_pair(theMAD,theMAD/histo->GetEntries());
1295 
1296  return result;
1297 
1298 }
1299 
1300 //*************************************************************
1301 std::pair<std::pair<Double_t,Double_t>, std::pair<Double_t,Double_t> > PrimaryVertexValidation::fitResiduals(TH1 *hist)
1302 //*************************************************************
1303 {
1304  //float fitResult(9999);
1305  //if (hist->GetEntries() < 20) return ;
1306 
1307  float mean = hist->GetMean();
1308  float sigma = hist->GetRMS();
1309 
1310  TF1 func("tmp", "gaus", mean - 1.5*sigma, mean + 1.5*sigma);
1311  if (0 == hist->Fit(&func,"QNR")) { // N: do not blow up file by storing fit!
1312  mean = func.GetParameter(1);
1313  sigma = func.GetParameter(2);
1314  // second fit: three sigma of first fit around mean of first fit
1315  func.SetRange(mean - 2*sigma, mean + 2*sigma);
1316  // I: integral gives more correct results if binning is too wide
1317  // L: Likelihood can treat empty bins correctly (if hist not weighted...)
1318  if (0 == hist->Fit(&func, "Q0LR")) {
1319  if (hist->GetFunction(func.GetName())) { // Take care that it is later on drawn:
1320  hist->GetFunction(func.GetName())->ResetBit(TF1::kNotDraw);
1321  }
1322  }
1323  }
1324 
1325  float res_mean = func.GetParameter(1);
1326  float res_width = func.GetParameter(2);
1327 
1328  float res_mean_err = func.GetParError(1);
1329  float res_width_err = func.GetParError(2);
1330 
1331  std::pair<Double_t,Double_t> resultM;
1332  std::pair<Double_t,Double_t> resultW;
1333 
1334  resultM = std::make_pair(res_mean,res_mean_err);
1335  resultW = std::make_pair(res_width,res_width_err);
1336 
1337  std::pair<std::pair<Double_t,Double_t>, std::pair<Double_t,Double_t> > result;
1338 
1339  result = std::make_pair(resultM,resultW);
1340  return result;
1341 }
1342 
1343 //*************************************************************
1344 void PrimaryVertexValidation::FillTrendPlot(TH1F* trendPlot, TH1F* residualsPlot[100], const TString& fitPar_, const TString& var_)
1345 //*************************************************************
1346 {
1347 
1348  float phiInterval = (360.)/nBins_;
1349  float etaInterval = 5./nBins_;
1350 
1351  for ( int i=0; i<nBins_; ++i ) {
1352 
1353  char phipositionString[129];
1354  float phiposition = (-180+i*phiInterval)+(phiInterval/2);
1355  sprintf(phipositionString,"%.f",phiposition);
1356 
1357  char etapositionString[129];
1358  float etaposition = (-2.5+i*etaInterval)+(etaInterval/2);
1359  sprintf(etapositionString,"%.1f",etaposition);
1360 
1361  if(fitPar_=="mean"){
1362  float mean_ = fitResiduals(residualsPlot[i]).first.first;
1363  float meanErr_ = fitResiduals(residualsPlot[i]).first.second;
1364  trendPlot->SetBinContent(i+1,mean_);
1365  trendPlot->SetBinError(i+1,meanErr_);
1366  } else if (fitPar_=="width"){
1367  float width_ = fitResiduals(residualsPlot[i]).second.first;
1368  float widthErr_ = fitResiduals(residualsPlot[i]).second.second;
1369  trendPlot->SetBinContent(i+1,width_);
1370  trendPlot->SetBinError(i+1,widthErr_);
1371  } else if (fitPar_=="median"){
1372  float median_ = getMedian(residualsPlot[i]).first;
1373  float medianErr_ = getMedian(residualsPlot[i]).second;
1374  trendPlot->SetBinContent(i+1,median_);
1375  trendPlot->SetBinError(i+1,medianErr_);
1376  } else if (fitPar_=="mad"){
1377  float mad_ = getMAD(residualsPlot[i]).first;
1378  float madErr_ = getMAD(residualsPlot[i]).second;
1379  trendPlot->SetBinContent(i+1,mad_);
1380  trendPlot->SetBinError(i+1,madErr_);
1381  } else {
1382  std::cout<<"PrimaryVertexValidation::FillTrendPlot() "<<fitPar_<<" unknown estimator!"<<std::endl;
1383  }
1384 
1385  if(var_=="eta"){
1386  trendPlot->GetXaxis()->SetBinLabel(i+1,phipositionString);
1387  } else if(var_=="phi"){
1388  trendPlot->GetXaxis()->SetBinLabel(i+1,etapositionString);
1389  } else {
1390  std::cout<<"PrimaryVertexValidation::FillTrendPlot() "<<var_<<" unknown track parameter!"<<std::endl;
1391  }
1392  }
1393 }
1394 
1395 
1396 
1397 //define this as a plug-in
1398 DEFINE_FWK_MODULE(PrimaryVertexValidation);
cms::Exception::what
virtual char const * what() const
Definition: Exception.cc:141
TransientVertex::positionError
GlobalError positionError() const
Definition: TransientVertex.h:114
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:44
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:232
SiPixelRawToDigiRegional_cfi.beamSpot
tuple beamSpot
Definition: SiPixelRawToDigiRegional_cfi.py:11
reco::HitPattern::getLayer
static uint32_t getLayer(uint32_t pattern)
Definition: HitPattern.h:520
reco::Vertex::tracks_end
trackRef_iterator tracks_end() const
last iterator over tracks
Definition: Vertex.cc:44
PrimaryVertexValidation::theTrackFilter_
TrackFilterForPVFindingBase * theTrackFilter_
Definition: PrimaryVertexValidation.h:107
PrimaryVertexValidation::a_dxyResidualsMap
TH1F * a_dxyResidualsMap[nMaxBins_][nMaxBins_]
Definition: PrimaryVertexValidation.h:228
reco::TrackBase::TrackQuality
TrackQuality
track quality
Definition: TrackBase.h:93
PrimaryVertexValidation::n_dxyEtaBiasResiduals
TH1F * n_dxyEtaBiasResiduals[nMaxBins_]
Definition: PrimaryVertexValidation.h:305
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
PrimaryVertexValidation::PrimaryVertexValidation
PrimaryVertexValidation(const edm::ParameterSet &)
Definition: PrimaryVertexValidation.cc:75
PrimaryVertexValidation::getMAD
std::pair< Double_t, Double_t > getMAD(TH1F *histo)
Definition: PrimaryVertexValidation.cc:1271
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:315
TransientVertex::totalChiSquared
float totalChiSquared() const
Definition: TransientVertex.h:124
reco::TransientTrack::hitPattern
const HitPattern & hitPattern() const
Definition: TransientTrack.h:132
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:191
PrimaryVertexValidation::a_dzPhiBiasResiduals
TH1F * a_dzPhiBiasResiduals[nMaxBins_]
Definition: PrimaryVertexValidation.h:299
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:308
PrimaryVertexValidation::analyze
virtual void analyze(const edm::Event &, const edm::EventSetup &)
Definition: PrimaryVertexValidation.cc:123
PrimaryVertexValidation::dszFromMyVertex_
double dszFromMyVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:204
PrimaryVertexValidation::dxyFromMyVertex_
double dxyFromMyVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:195
reco::Vertex::position
const Point & position() const
position
Definition: Vertex.h:106
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
PrimaryVertexValidation::dzErrorFromMyVertex_
double dzErrorFromMyVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:199
edm::EventAuxiliary::luminosityBlock
LuminosityBlockNumber_t luminosityBlock() const
Definition: EventAuxiliary.h:63
PrimaryVertexValidation::a_dxyPhiBiasResiduals
TH1F * a_dxyPhiBiasResiduals[nMaxBins_]
Definition: PrimaryVertexValidation.h:296
PrimaryVertexValidation::IPTsigFromMyVertex_
double IPTsigFromMyVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:201
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:1245
kFPIX
const int kFPIX
Definition: AlignmentTrackSelector.cc:31
TransientVertex::normalisedChiSquared
float normalisedChiSquared() const
Definition: TransientVertex.h:125
PrimaryVertexValidation::tracksUsedForVertexing_
int tracksUsedForVertexing_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:193
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:307
TransientVertex::degreesOfFreedom
float degreesOfFreedom() const
Definition: TransientVertex.h:128
PrimaryVertexValidation::yUnbiasedVertex_
double yUnbiasedVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:187
PrimaryVertexValidation::fitResiduals
std::pair< std::pair< Double_t, Double_t >, std::pair< Double_t, Double_t > > fitResiduals(TH1 *hist)
Definition: PrimaryVertexValidation.cc:1301
TransientVertex::position
GlobalPoint position() const
Definition: TransientVertex.h:113
PV3DBase::z
T z() const
Definition: PV3DBase.h:64
query.result
tuple result
Definition: query.py:137
Abs
T Abs(T a)
Definition: MathUtil.h:49
reco::HitPattern::numberOfHits
int numberOfHits() const
Definition: HitPattern.cc:211
reco::Vertex::chi2
double chi2() const
chi-squares
Definition: Vertex.h:95
j
int j
Definition: DBlmapReader.cc:9
PrimaryVertexValidation::zUnbiasedVertex_
double zUnbiasedVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:188
reco::HitPattern::validHitFilter
static bool validHitFilter(uint32_t pattern)
Definition: HitPattern.h:564
PrimaryVertexValidation::chi2normUnbiasedVertex_
float chi2normUnbiasedVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:189
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:136
edm::Event::getByLabel
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:390
LogTrace
#define LogTrace(id)
Definition: MessageLogger.h:502
PrimaryVertexValidation::IPLsigFromMyVertex_
double IPLsigFromMyVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:202
reco::Vertex::ndof
double ndof() const
Definition: Vertex.h:102
gen::k
int k[5][pyjets_maxn]
Definition: Cascade2Hadronizer.cc:79
edm::Event::eventAuxiliary
EventAuxiliary const & eventAuxiliary() const
Definition: Event.h:72
DetId
Definition: DetId.h:18
reco::TrackBase::qualityByName
static TrackQuality qualityByName(const std::string &name)
Definition: TrackBase.cc:46
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:180
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:1344
edm::EventSetup::get
const T & get() const
Definition: EventSetup.h:55
PrimaryVertexValidation::n_dxyPhiBiasResiduals
TH1F * n_dxyPhiBiasResiduals[nMaxBins_]
Definition: PrimaryVertexValidation.h:304
PrimaryVertexValidation::theTrackClusterizer_
TrackClusterizerInZ * theTrackClusterizer_
Definition: PrimaryVertexValidation.h:108
PrimaryVertexValidation::xUnbiasedVertex_
double xUnbiasedVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:186
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:297
edm::EventAuxiliary::id
EventID const & id() const
Definition: EventAuxiliary.h:59
PrimaryVertexValidation::dzFromMyVertex_
double dzFromMyVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:196
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:312
reco::HitPattern::pixelEndcapHitFilter
static bool pixelEndcapHitFilter(uint32_t pattern)
Definition: HitPattern.h:434
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
PrimaryVertexValidation::n_dxyResidualsMap
TH1F * n_dxyResidualsMap[nMaxBins_][nMaxBins_]
Definition: PrimaryVertexValidation.h:231
reco::Vertex::normalizedChi2
double normalizedChi2() const
chi-squared divided by n.d.o.f.
Definition: Vertex.h:104
PrimaryVertexValidation::dxyErrorFromMyVertex_
double dxyErrorFromMyVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:198
TrackingRecHit::geographicalId
DetId geographicalId() const
Definition: TrackingRecHit.h:120
x
Definition: DDAxes.h:10
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:192
PrimaryVertexValidation::n_dzBiasResidualsMap
TH1F * n_dzBiasResidualsMap[nMaxBins_][nMaxBins_]
Definition: PrimaryVertexValidation.h:316
reco::HitPattern::getHitPattern
uint32_t getHitPattern(int position) const
Definition: HitPattern.cc:142
PrimaryVertexValidation::a_dzBiasResidualsMap
TH1F * a_dzBiasResidualsMap[nMaxBins_][nMaxBins_]
Definition: PrimaryVertexValidation.h:313
PV3DBase::x
T x() const
Definition: PV3DBase.h:62
PrimaryVertexValidation::a_dzEtaBiasResiduals
TH1F * a_dzEtaBiasResiduals[nMaxBins_]
Definition: PrimaryVertexValidation.h:300
TransientVertex::isValid
bool isValid() const
Definition: TransientVertex.h:133
PrimaryVertexValidation::a_dzResidualsMap
TH1F * a_dzResidualsMap[nMaxBins_][nMaxBins_]
Definition: PrimaryVertexValidation.h:229
reco::HitPattern::pixelBarrelHitFilter
static bool pixelBarrelHitFilter(uint32_t pattern)
Definition: HitPattern.h:427
PrimaryVertexValidation::chi2UnbiasedVertex_
float chi2UnbiasedVertex_[nMaxtracks_]
Definition: PrimaryVertexValidation.h:190
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
GlobalVector
Global3DVector GlobalVector
Definition: GlobalVector.h:10
kBPIX
const int kBPIX
Definition: AlignmentTrackSelector.cc:30
Generated for CMSSW Reference Manual by doxygen 1.8.5