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