OfflineValidation/plugins/PrimaryVertexValidation.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    Alignment/OfflineValidation
0004 // Class:      PrimaryVertexValidation
0005 //
0006 /**\class PrimaryVertexValidation PrimaryVertexValidation.cc Alignment/OfflineValidation/plugins/PrimaryVertexValidation.cc
0007
0008  Description: Validate alignment constants using unbiased vertex residuals
0009
0010  Implementation:
0011  <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Marco Musich
0015 //         Created:  Tue Mar 02 10:39:34 CDT 2010
0016 //
0017
0018 // system include files
0019 #include <memory>
0020 #include <vector>
0021 #include <regex>
0022 #include <cassert>
0023 #include <chrono>
0024 #include <iomanip>
0025 #include <boost/range/adaptor/indexed.hpp>
0026
0027 // user include files
0028 #include "Alignment/OfflineValidation/plugins/PrimaryVertexValidation.h"
0029
0030 // ROOT includes
0031 #include "TH1F.h"
0032 #include "TH2F.h"
0033 #include "TF1.h"
0034 #include "TVector3.h"
0035 #include "TFile.h"
0036 #include "TMath.h"
0037 #include "TROOT.h"
0038 #include "TChain.h"
0039 #include "TNtuple.h"
0040 #include "TMatrixD.h"
0041 #include "TVectorD.h"
0042
0043 // CMSSW includes
0044 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0045 #include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
0046 #include "DataFormats/TrackReco/interface/Track.h"
0047 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0048 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0049 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
0050 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0051 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0052 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0053 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0054 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0055 #include "Geometry/TrackerGeometryBuilder/interface/PixelTopologyMap.h"
0056 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0057 #include "RecoVertex/PrimaryVertexProducer/interface/TrackFilterForPVFinding.h"
0058 #include "RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZ_vect.h"
0059 #include "RecoVertex/PrimaryVertexProducer/interface/GapClusterizerInZ.h"
0060 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
0061 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0062 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateClosestToPoint.h"
0063
0064 const int PrimaryVertexValidation::nMaxtracks_;
0065
0066 // Constructor
0067 PrimaryVertexValidation::PrimaryVertexValidation(const edm::ParameterSet& iConfig)
0068     : magFieldToken_(esConsumes<MagneticField, IdealMagneticFieldRecord>()),
0069       trackingGeomToken_(esConsumes<GlobalTrackingGeometry, GlobalTrackingGeometryRecord>()),
0070       ttkToken_(esConsumes<TransientTrackBuilder, TransientTrackRecord>(edm::ESInputTag("", "TransientTrackBuilder"))),
0071       topoToken_(esConsumes<TrackerTopology, TrackerTopologyRcd>()),
0072       topoTokenBR_(esConsumes<TrackerTopology, TrackerTopologyRcd, edm::Transition::BeginRun>()),
0073       geomTokenBR_(esConsumes<TrackerGeometry, TrackerDigiGeometryRecord, edm::Transition::BeginRun>()),
0074       runInfoTokenBR_(esConsumes<RunInfo, RunInfoRcd, edm::Transition::BeginRun>()),
0075       compressionSettings_(iConfig.getUntrackedParameter<int>("compressionSettings", -1)),
0076       storeNtuple_(iConfig.getParameter<bool>("storeNtuple")),
0077       lightNtupleSwitch_(iConfig.getParameter<bool>("isLightNtuple")),
0078       useTracksFromRecoVtx_(iConfig.getParameter<bool>("useTracksFromRecoVtx")),
0079       vertexZMax_(iConfig.getUntrackedParameter<double>("vertexZMax", 99.)),
0080       intLumi_(iConfig.getUntrackedParameter<double>("intLumi", 0.)),
0081       askFirstLayerHit_(iConfig.getParameter<bool>("askFirstLayerHit")),
0082       doBPix_(iConfig.getUntrackedParameter<bool>("doBPix", true)),
0083       doFPix_(iConfig.getUntrackedParameter<bool>("doFPix", true)),
0084       ptOfProbe_(iConfig.getUntrackedParameter<double>("probePt", 0.)),
0085       pOfProbe_(iConfig.getUntrackedParameter<double>("probeP", 0.)),
0086       etaOfProbe_(iConfig.getUntrackedParameter<double>("probeEta", 2.4)),
0087       nHitsOfProbe_(iConfig.getUntrackedParameter<double>("probeNHits", 0.)),
0088       nBins_(iConfig.getUntrackedParameter<int>("numberOfBins", 24)),
0089       minPt_(iConfig.getUntrackedParameter<double>("minPt", 1.)),
0090       maxPt_(iConfig.getUntrackedParameter<double>("maxPt", 20.)),
0091       debug_(iConfig.getParameter<bool>("Debug")),
0092       runControl_(iConfig.getUntrackedParameter<bool>("runControl", false)),
0093       forceBeamSpotContraint_(iConfig.getUntrackedParameter<bool>("forceBeamSpot", false)) {
0094   // now do what ever initialization is needed
0095   // initialize phase space boundaries
0096
0097   usesResource(TFileService::kSharedResource);
0098
0099   std::vector<unsigned int> defaultRuns;
0100   defaultRuns.push_back(0);
0101   runControlNumbers_ = iConfig.getUntrackedParameter<std::vector<unsigned int>>("runControlNumber", defaultRuns);
0102
0103   edm::InputTag TrackCollectionTag_ = iConfig.getParameter<edm::InputTag>("TrackCollectionTag");
0104   theTrackCollectionToken_ = consumes<reco::TrackCollection>(TrackCollectionTag_);
0105
0106   edm::InputTag VertexCollectionTag_ = iConfig.getParameter<edm::InputTag>("VertexCollectionTag");
0107   theVertexCollectionToken_ = consumes<reco::VertexCollection>(VertexCollectionTag_);
0108
0109   edm::InputTag BeamspotTag_ = iConfig.getParameter<edm::InputTag>("BeamSpotTag");
0110   theBeamspotToken_ = consumes<reco::BeamSpot>(BeamspotTag_);
0111
0112   // select and configure the track filter
0113   theTrackFilter_ =
0114       std::make_unique<TrackFilterForPVFinding>(iConfig.getParameter<edm::ParameterSet>("TkFilterParameters"));
0115   // select and configure the track clusterizer
0116   std::string clusteringAlgorithm =
0117       iConfig.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<std::string>("algorithm");
0118   if (clusteringAlgorithm == "gap") {
0119     theTrackClusterizer_ =
0120         std::make_unique<GapClusterizerInZ>(iConfig.getParameter<edm::ParameterSet>("TkClusParameters")
0121                                                 .getParameter<edm::ParameterSet>("TkGapClusParameters"));
0122   } else if (clusteringAlgorithm == "DA_vect") {
0123     theTrackClusterizer_ =
0124         std::make_unique<DAClusterizerInZ_vect>(iConfig.getParameter<edm::ParameterSet>("TkClusParameters")
0125                                                     .getParameter<edm::ParameterSet>("TkDAClusParameters"));
0126   } else {
0127     throw VertexException("PrimaryVertexProducerAlgorithm: unknown clustering algorithm: " + clusteringAlgorithm);
0128   }
0129
0130   theDetails_.histobins = 500;
0131   theDetails_.setMap(PVValHelper::dxy, PVValHelper::phi, -2000., 2000.);
0132   theDetails_.setMap(PVValHelper::dxy, PVValHelper::eta, -3000., 3000.);
0133   theDetails_.setMap(PVValHelper::dxy, PVValHelper::pT, -1000., 1000.);
0134   theDetails_.setMap(PVValHelper::dxy, PVValHelper::pTCentral, -1000., 1000.);
0135   theDetails_.setMap(PVValHelper::dxy, PVValHelper::ladder, -1000., 1000.);
0136   theDetails_.setMap(PVValHelper::dxy, PVValHelper::modZ, -1000., 1000.);
0137
0138   for (int i = PVValHelper::phi; i < PVValHelper::END_OF_PLOTS; i++) {
0139     for (int j = PVValHelper::dx; j < PVValHelper::END_OF_TYPES; j++) {
0140       auto plot_index = static_cast<PVValHelper::plotVariable>(i);
0141       auto res_index = static_cast<PVValHelper::residualType>(j);
0142
0143       if (debug_) {
0144         edm::LogInfo("PrimaryVertexValidation")
0145             << "==> " << std::get<0>(PVValHelper::getTypeString(res_index)) << " " << std::setw(10)
0146             << std::get<0>(PVValHelper::getVarString(plot_index)) << std::endl;
0147       }
0148       if (res_index != PVValHelper::d3D && res_index != PVValHelper::norm_d3D)
0149         theDetails_.setMap(res_index,
0150                            plot_index,
0151                            theDetails_.getLow(PVValHelper::dxy, plot_index),
0152                            theDetails_.getHigh(PVValHelper::dxy, plot_index));
0153       else
0154         theDetails_.setMap(res_index, plot_index, 0., theDetails_.getHigh(PVValHelper::dxy, plot_index));
0155     }
0156   }
0157
0158   edm::LogVerbatim("PrimaryVertexValidation") << "######################################";
0159   for (const auto& it : theDetails_.range) {
0160     edm::LogVerbatim("PrimaryVertexValidation")
0161         << "|" << std::setw(10) << std::get<0>(PVValHelper::getTypeString(it.first.first)) << "|" << std::setw(10)
0162         << std::get<0>(PVValHelper::getVarString(it.first.second)) << "| (" << std::setw(5) << it.second.first << ";"
0163         << std::setw(5) << it.second.second << ") |" << std::endl;
0164   }
0165
0166   theDetails_.trendbins[PVValHelper::phi] = PVValHelper::generateBins(nBins_ + 1, -180., 360.);
0167   theDetails_.trendbins[PVValHelper::eta] = PVValHelper::generateBins(nBins_ + 1, -etaOfProbe_, 2 * etaOfProbe_);
0168
0169   if (debug_) {
0170     edm::LogVerbatim("PrimaryVertexValidation") << "etaBins: ";
0171     for (auto ieta : theDetails_.trendbins[PVValHelper::eta]) {
0172       edm::LogVerbatim("PrimaryVertexValidation") << ieta << " ";
0173     }
0174     edm::LogVerbatim("PrimaryVertexValidation") << "\n";
0175
0176     edm::LogVerbatim("PrimaryVertexValidation") << "phiBins: ";
0177     for (auto iphi : theDetails_.trendbins[PVValHelper::phi]) {
0178       edm::LogVerbatim("PrimaryVertexValidation") << iphi << " ";
0179     }
0180     edm::LogVerbatim("PrimaryVertexValidation") << "\n";
0181   }
0182
0183   // create the bins of the pT-binned distributions
0184
0185   mypT_bins_ = PVValHelper::makeLogBins<float, nPtBins_>(minPt_, maxPt_);
0186
0187   std::string toOutput = "";
0188   for (auto ptbin : mypT_bins_) {
0189     toOutput += " ";
0190     toOutput += std::to_string(ptbin);
0191     toOutput += ",";
0192   }
0193
0194   edm::LogVerbatim("PrimaryVertexValidation") << "######################################\n";
0195   edm::LogVerbatim("PrimaryVertexValidation") << "The pT binning is: [" << toOutput << "] \n";
0196 }
0197
0198 // Destructor
0199 PrimaryVertexValidation::~PrimaryVertexValidation() = default;
0200 //
0201 // member functions
0202 //
0203
0204 // ------------ method called to for each event  ------------
0205 void PrimaryVertexValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0206   using namespace std;
0207   using namespace reco;
0208   using namespace IPTools;
0209
0210   if (nBins_ != 24 && debug_) {
0211     edm::LogInfo("PrimaryVertexValidation") << "Using: " << nBins_ << " bins plots";
0212   }
0213
0214   bool passesRunControl = false;
0215
0216   if (runControl_) {
0217     for (const auto& runControlNumber : runControlNumbers_) {
0218       if (iEvent.eventAuxiliary().run() == runControlNumber) {
0219         if (debug_) {
0220           edm::LogInfo("PrimaryVertexValidation")
0221               << " run number: " << iEvent.eventAuxiliary().run() << " keeping run:" << runControlNumber;
0222         }
0223         passesRunControl = true;
0224         break;
0225       }
0226     }
0227     if (!passesRunControl)
0228       return;
0229   }
0230
0231   Nevt_++;
0232
0233   //=======================================================
0234   // Initialize Root-tuple variables
0235   //=======================================================
0236
0237   SetVarToZero();
0238
0239   //=======================================================
0240   // Retrieve tracker topology from geometry
0241   //=======================================================
0242
0243   const TrackerTopology* const tTopo = &iSetup.getData(topoToken_);
0244
0245   //=======================================================
0246   // Retrieve the Magnetic Field information
0247   //=======================================================
0248
0249   edm::ESHandle<MagneticField> theMGField = iSetup.getHandle(magFieldToken_);
0250
0251   //=======================================================
0252   // Retrieve the Tracking Geometry information
0253   //=======================================================
0254
0255   edm::ESHandle<GlobalTrackingGeometry> theTrackingGeometry = iSetup.getHandle(trackingGeomToken_);
0256
0257   //=======================================================
0258   // Retrieve the Transient Track Builder information
0259   //=======================================================
0260
0261   edm::ESHandle<TransientTrackBuilder> theB_ = iSetup.getHandle(ttkToken_);
0262   double fBfield_ = ((*theB_).field()->inTesla(GlobalPoint(0., 0., 0.))).z();
0263
0264   //=======================================================
0265   // Retrieve the Track information
0266   //=======================================================
0267
0268   edm::Handle<TrackCollection> trackCollectionHandle = iEvent.getHandle(theTrackCollectionToken_);
0269   if (!trackCollectionHandle.isValid())
0270     return;
0271   auto const& tracks = *trackCollectionHandle;
0272
0273   //=======================================================
0274   // Retrieve offline vartex information (only for reco)
0275   //=======================================================
0276
0277   //edm::Handle<VertexCollection> vertices;
0278   edm::Handle<std::vector<Vertex>> vertices;
0279   vertices = iEvent.getHandle(theVertexCollectionToken_);
0280   if (!vertices.isValid()) {
0281     edm::LogError("PrimaryVertexValidation") << "Vertex collection handle is not valid. Aborting!" << std::endl;
0282     return;
0283   }
0284
0285   std::vector<Vertex> vsorted = *(vertices);
0286   // sort the vertices by number of tracks in descending order
0287   // use chi2 as tiebreaker
0288   std::sort(vsorted.begin(), vsorted.end(), PrimaryVertexValidation::vtxSort);
0289
0290   // skip events with no PV, this should not happen
0291   if (vsorted.empty())
0292     return;
0293
0294   // skip events failing vertex cut
0295   if (std::abs(vsorted[0].z()) > vertexZMax_)
0296     return;
0297
0298   if (vsorted[0].isValid()) {
0299     xOfflineVertex_ = (vsorted)[0].x();
0300     yOfflineVertex_ = (vsorted)[0].y();
0301     zOfflineVertex_ = (vsorted)[0].z();
0302
0303     xErrOfflineVertex_ = (vsorted)[0].xError();
0304     yErrOfflineVertex_ = (vsorted)[0].yError();
0305     zErrOfflineVertex_ = (vsorted)[0].zError();
0306   }
0307
0308   h_xOfflineVertex->Fill(xOfflineVertex_);
0309   h_yOfflineVertex->Fill(yOfflineVertex_);
0310   h_zOfflineVertex->Fill(zOfflineVertex_);
0311   h_xErrOfflineVertex->Fill(xErrOfflineVertex_);
0312   h_yErrOfflineVertex->Fill(yErrOfflineVertex_);
0313   h_zErrOfflineVertex->Fill(zErrOfflineVertex_);
0314
0315   unsigned int vertexCollectionSize = vsorted.size();
0316   int nvvertex = 0;
0317
0318   for (unsigned int i = 0; i < vertexCollectionSize; i++) {
0319     const Vertex& vertex = vsorted.at(i);
0320     if (vertex.isValid())
0321       nvvertex++;
0322   }
0323
0324   nOfflineVertices_ = nvvertex;
0325   h_nOfflineVertices->Fill(nvvertex);
0326
0327   if (!vsorted.empty() && useTracksFromRecoVtx_) {
0328     double sumpt = 0;
0329     size_t ntracks = 0;
0330     double chi2ndf = 0.;
0331     double chi2prob = 0.;
0332
0333     if (!vsorted.at(0).isFake()) {
0334       Vertex pv = vsorted.at(0);
0335
0336       ntracks = pv.tracksSize();
0337       chi2ndf = pv.normalizedChi2();
0338       chi2prob = TMath::Prob(pv.chi2(), (int)pv.ndof());
0339
0340       h_recoVtxNtracks_->Fill(ntracks);
0341       h_recoVtxChi2ndf_->Fill(chi2ndf);
0342       h_recoVtxChi2Prob_->Fill(chi2prob);
0343
0344       for (Vertex::trackRef_iterator itrk = pv.tracks_begin(); itrk != pv.tracks_end(); ++itrk) {
0345         double pt = (**itrk).pt();
0346         sumpt += pt * pt;
0347
0348         const math::XYZPoint myVertex(pv.position().x(), pv.position().y(), pv.position().z());
0349
0350         double dxyRes = (**itrk).dxy(myVertex);
0351         double dzRes = (**itrk).dz(myVertex);
0352
0353         double dxy_err = (**itrk).dxyError();
0354         double dz_err = (**itrk).dzError();
0355
0356         float trackphi = ((**itrk).phi()) * (180 / M_PI);
0357         float tracketa = (**itrk).eta();
0358
0359         for (int i = 0; i < nBins_; i++) {
0360           float phiF = theDetails_.trendbins[PVValHelper::phi][i];
0361           float phiL = theDetails_.trendbins[PVValHelper::phi][i + 1];
0362
0363           float etaF = theDetails_.trendbins[PVValHelper::eta][i];
0364           float etaL = theDetails_.trendbins[PVValHelper::eta][i + 1];
0365
0366           if (tracketa >= etaF && tracketa < etaL) {
0367             PVValHelper::fillByIndex(a_dxyEtaBiasResiduals, i, dxyRes * cmToum, "1");
0368             PVValHelper::fillByIndex(a_dzEtaBiasResiduals, i, dzRes * cmToum, "2");
0369             PVValHelper::fillByIndex(n_dxyEtaBiasResiduals, i, (dxyRes) / dxy_err, "3");
0370             PVValHelper::fillByIndex(n_dzEtaBiasResiduals, i, (dzRes) / dz_err, "4");
0371           }
0372
0373           if (trackphi >= phiF && trackphi < phiL) {
0374             PVValHelper::fillByIndex(a_dxyPhiBiasResiduals, i, dxyRes * cmToum, "5");
0375             PVValHelper::fillByIndex(a_dzPhiBiasResiduals, i, dzRes * cmToum, "6");
0376             PVValHelper::fillByIndex(n_dxyPhiBiasResiduals, i, (dxyRes) / dxy_err, "7");
0377             PVValHelper::fillByIndex(n_dzPhiBiasResiduals, i, (dzRes) / dz_err, "8");
0378
0379             for (int j = 0; j < nBins_; j++) {
0380               float etaJ = theDetails_.trendbins[PVValHelper::eta][j];
0381               float etaK = theDetails_.trendbins[PVValHelper::eta][j + 1];
0382
0383               if (tracketa >= etaJ && tracketa < etaK) {
0384                 a_dxyBiasResidualsMap[i][j]->Fill(dxyRes * cmToum);
0385                 a_dzBiasResidualsMap[i][j]->Fill(dzRes * cmToum);
0386
0387                 n_dxyBiasResidualsMap[i][j]->Fill((dxyRes) / dxy_err);
0388                 n_dzBiasResidualsMap[i][j]->Fill((dzRes) / dz_err);
0389               }
0390             }
0391           }
0392         }
0393       }
0394
0395       h_recoVtxSumPt_->Fill(sumpt);
0396     }
0397   }
0398
0399   //=======================================================
0400   // Retrieve Beamspot information
0401   //=======================================================
0402
0403   BeamSpot beamSpot;
0404   edm::Handle<BeamSpot> beamSpotHandle = iEvent.getHandle(theBeamspotToken_);
0405
0406   if (beamSpotHandle.isValid()) {
0407     beamSpot = *beamSpotHandle;
0408     BSx0_ = beamSpot.x0();
0409     BSy0_ = beamSpot.y0();
0410     BSz0_ = beamSpot.z0();
0411     Beamsigmaz_ = beamSpot.sigmaZ();
0412     Beamdxdz_ = beamSpot.dxdz();
0413     BeamWidthX_ = beamSpot.BeamWidthX();
0414     BeamWidthY_ = beamSpot.BeamWidthY();
0415
0416     wxy2_ = TMath::Power(BeamWidthX_, 2) + TMath::Power(BeamWidthY_, 2);
0417
0418   } else {
0419     edm::LogWarning("PrimaryVertexValidation") << "No BeamSpot found!";
0420   }
0421
0422   h_BSx0->Fill(BSx0_);
0423   h_BSy0->Fill(BSy0_);
0424   h_BSz0->Fill(BSz0_);
0425   h_Beamsigmaz->Fill(Beamsigmaz_);
0426   h_BeamWidthX->Fill(BeamWidthX_);
0427   h_BeamWidthY->Fill(BeamWidthY_);
0428
0429   if (debug_)
0430     edm::LogInfo("PrimaryVertexValidation") << "Beamspot x:" << BSx0_ << " y:" << BSy0_ << " z:" << BSz0_;
0431
0432   //=======================================================
0433   // Starts here ananlysis
0434   //=======================================================
0435
0436   RunNumber_ = iEvent.eventAuxiliary().run();
0437   LuminosityBlockNumber_ = iEvent.eventAuxiliary().luminosityBlock();
0438   EventNumber_ = iEvent.eventAuxiliary().id().event();
0439
0440   if (debug_)
0441     edm::LogInfo("PrimaryVertexValidation") << " looping over " << trackCollectionHandle->size() << "tracks";
0442
0443   h_nTracks->Fill(trackCollectionHandle->size());
0444
0445   //======================================================
0446   // Interface RECO tracks to vertex reconstruction
0447   //======================================================
0448
0449   std::vector<TransientTrack> t_tks;
0450   for (const auto& track : tracks) {
0451     TransientTrack tt = theB_->build(&(track));
0452     tt.setBeamSpot(beamSpot);
0453     t_tks.push_back(tt);
0454   }
0455
0456   if (debug_) {
0457     edm::LogInfo("PrimaryVertexValidation") << "Found: " << t_tks.size() << " reconstructed tracks";
0458   }
0459
0460   //======================================================
0461   // select the tracks
0462   //======================================================
0463
0464   std::vector<TransientTrack> seltks = theTrackFilter_->select(t_tks);
0465
0466   //======================================================
0467   // clusterize tracks in Z
0468   //======================================================
0469
0470   vector<vector<TransientTrack>> clusters = theTrackClusterizer_->clusterize(seltks);
0471
0472   if (debug_) {
0473     edm::LogInfo("PrimaryVertexValidation")
0474         << " looping over: " << clusters.size() << " clusters  from " << t_tks.size() << " selected tracks";
0475   }
0476
0477   nClus_ = clusters.size();
0478   h_nClus->Fill(nClus_);
0479
0480   //======================================================
0481   // Starts loop on clusters
0482   //======================================================
0483   for (const auto& iclus : clusters) {
0484     nTracksPerClus_ = 0;
0485
0486     unsigned int i = 0;
0487     for (const auto& theTTrack : iclus) {
0488       i++;
0489
0490       if (nTracks_ >= nMaxtracks_) {
0491         edm::LogError("PrimaryVertexValidation")
0492             << " Warning - Number of tracks: " << nTracks_ << " , greater than " << nMaxtracks_;
0493         continue;
0494       }
0495
0496       const Track& theTrack = theTTrack.track();
0497
0498       pt_[nTracks_] = theTrack.pt();
0499       p_[nTracks_] = theTrack.p();
0500       nhits_[nTracks_] = theTrack.numberOfValidHits();
0501       eta_[nTracks_] = theTrack.eta();
0502       theta_[nTracks_] = theTrack.theta();
0503       phi_[nTracks_] = theTrack.phi();
0504       chi2_[nTracks_] = theTrack.chi2();
0505       chi2ndof_[nTracks_] = theTrack.normalizedChi2();
0506       charge_[nTracks_] = theTrack.charge();
0507       qoverp_[nTracks_] = theTrack.qoverp();
0508       dz_[nTracks_] = theTrack.dz();
0509       dxy_[nTracks_] = theTrack.dxy();
0510
0511       TrackBase::TrackQuality _trackQuality = TrackBase::qualityByName("highPurity");
0512       isHighPurity_[nTracks_] = theTrack.quality(_trackQuality);
0513
0514       math::XYZPoint point(BSx0_, BSy0_, BSz0_);
0515       dxyBs_[nTracks_] = theTrack.dxy(point);
0516       dzBs_[nTracks_] = theTrack.dz(point);
0517
0518       xPCA_[nTracks_] = theTrack.vertex().x();
0519       yPCA_[nTracks_] = theTrack.vertex().y();
0520       zPCA_[nTracks_] = theTrack.vertex().z();
0521
0522       //=======================================================
0523       // Retrieve rechit information
0524       //=======================================================
0525
0526       const reco::HitPattern& hits = theTrack.hitPattern();
0527
0528       int nRecHit1D = 0;
0529       int nRecHit2D = 0;
0530       int nhitinTIB = hits.numberOfValidStripTIBHits();
0531       int nhitinTOB = hits.numberOfValidStripTOBHits();
0532       int nhitinTID = hits.numberOfValidStripTIDHits();
0533       int nhitinTEC = hits.numberOfValidStripTECHits();
0534       int nhitinBPIX = hits.numberOfValidPixelBarrelHits();
0535       int nhitinFPIX = hits.numberOfValidPixelEndcapHits();
0536       for (trackingRecHit_iterator iHit = theTTrack.recHitsBegin(); iHit != theTTrack.recHitsEnd(); ++iHit) {
0537         if ((*iHit)->isValid()) {
0538           if (this->isHit2D(**iHit, phase_)) {
0539             ++nRecHit2D;
0540           } else {
0541             ++nRecHit1D;
0542           }
0543         }
0544       }
0545
0546       nhits1D_[nTracks_] = nRecHit1D;
0547       nhits2D_[nTracks_] = nRecHit2D;
0548       nhitsBPIX_[nTracks_] = nhitinBPIX;
0549       nhitsFPIX_[nTracks_] = nhitinFPIX;
0550       nhitsTIB_[nTracks_] = nhitinTIB;
0551       nhitsTID_[nTracks_] = nhitinTID;
0552       nhitsTOB_[nTracks_] = nhitinTOB;
0553       nhitsTEC_[nTracks_] = nhitinTEC;
0554
0555       //=======================================================
0556       // Good tracks for vertexing selection
0557       //=======================================================
0558
0559       bool pass = true;
0560       if (askFirstLayerHit_)
0561         pass = this->hasFirstLayerPixelHits(theTTrack);
0562       if (pass && (theTrack.pt() >= ptOfProbe_) && std::abs(theTrack.eta()) <= etaOfProbe_ &&
0563           (theTrack.numberOfValidHits()) >= nHitsOfProbe_ && (theTrack.p()) >= pOfProbe_) {
0564         isGoodTrack_[nTracks_] = 1;
0565       }
0566
0567       //=======================================================
0568       // Fit unbiased vertex
0569       //=======================================================
0570
0571       vector<TransientTrack> theFinalTracks;
0572       theFinalTracks.clear();
0573
0574       for (const auto& tk : iclus) {
0575         pass = this->hasFirstLayerPixelHits(tk);
0576         if (pass) {
0577           if (tk == theTTrack)
0578             continue;
0579           else {
0580             theFinalTracks.push_back(tk);
0581           }
0582         }
0583       }
0584
0585       if (theFinalTracks.size() > 1) {
0586         if (debug_)
0587           edm::LogInfo("PrimaryVertexValidation") << "Transient Track Collection size: " << theFinalTracks.size();
0588         try {
0589           //AdaptiveVertexFitter* theFitter = new AdaptiveVertexFitter;
0590           auto theFitter = std::unique_ptr<VertexFitter<5>>(new AdaptiveVertexFitter());
0591           TransientVertex theFittedVertex;
0592
0593           if (forceBeamSpotContraint_) {
0594             theFittedVertex = theFitter->vertex(theFinalTracks, beamSpot);  // if you want the beam constraint
0595           } else {
0596             theFittedVertex = theFitter->vertex(theFinalTracks);
0597           }
0598
0599           double totalTrackWeights = 0;
0600           if (theFittedVertex.isValid()) {
0601             if (theFittedVertex.hasTrackWeight()) {
0602               for (const auto& theFinalTrack : theFinalTracks) {
0603                 sumOfWeightsUnbiasedVertex_[nTracks_] += theFittedVertex.trackWeight(theFinalTrack);
0604                 totalTrackWeights += theFittedVertex.trackWeight(theFinalTrack);
0605                 h_fitVtxTrackWeights_->Fill(theFittedVertex.trackWeight(theFinalTrack));
0606               }
0607             }
0608
0609             h_fitVtxTrackAverageWeight_->Fill(totalTrackWeights / theFinalTracks.size());
0610
0611             const math::XYZPoint theRecoVertex(xOfflineVertex_, yOfflineVertex_, zOfflineVertex_);
0612             const math::XYZPoint myVertex(
0613                 theFittedVertex.position().x(), theFittedVertex.position().y(), theFittedVertex.position().z());
0614
0615             const Vertex vertex = theFittedVertex;
0616             fillTrackHistos(hDA, "all", &theTTrack, vertex, beamSpot, fBfield_);
0617
0618             hasRecVertex_[nTracks_] = 1;
0619             xUnbiasedVertex_[nTracks_] = theFittedVertex.position().x();
0620             yUnbiasedVertex_[nTracks_] = theFittedVertex.position().y();
0621             zUnbiasedVertex_[nTracks_] = theFittedVertex.position().z();
0622
0623             chi2normUnbiasedVertex_[nTracks_] = theFittedVertex.normalisedChiSquared();
0624             chi2UnbiasedVertex_[nTracks_] = theFittedVertex.totalChiSquared();
0625             DOFUnbiasedVertex_[nTracks_] = theFittedVertex.degreesOfFreedom();
0626             chi2ProbUnbiasedVertex_[nTracks_] =
0627                 TMath::Prob(theFittedVertex.totalChiSquared(), (int)theFittedVertex.degreesOfFreedom());
0628             tracksUsedForVertexing_[nTracks_] = theFinalTracks.size();
0629
0630             h_fitVtxNtracks_->Fill(theFinalTracks.size());
0631             h_fitVtxChi2_->Fill(theFittedVertex.totalChiSquared());
0632             h_fitVtxNdof_->Fill(theFittedVertex.degreesOfFreedom());
0633             h_fitVtxChi2ndf_->Fill(theFittedVertex.normalisedChiSquared());
0634             h_fitVtxChi2Prob_->Fill(
0635                 TMath::Prob(theFittedVertex.totalChiSquared(), (int)theFittedVertex.degreesOfFreedom()));
0636
0637             // from my Vertex
0638             double dxyFromMyVertex = theTrack.dxy(myVertex);
0639             double dzFromMyVertex = theTrack.dz(myVertex);
0640
0641             GlobalPoint vert(
0642                 theFittedVertex.position().x(), theFittedVertex.position().y(), theFittedVertex.position().z());
0643
0644             //FreeTrajectoryState theTrackNearVertex = theTTrack.trajectoryStateClosestToPoint(vert).theState();
0645             //double dz_err = sqrt(theFittedVertex.positionError().czz() + theTrackNearVertex.cartesianError().position().czz());
0646             //double dz_err = hypot(theTrack.dzError(),theFittedVertex.positionError().czz());
0647
0648             double dz_err = sqrt(std::pow(theTrack.dzError(), 2) + theFittedVertex.positionError().czz());
0649
0650             // PV2D
0651             std::pair<bool, Measurement1D> s_ip2dpv = signedTransverseImpactParameter(
0652                 theTTrack, GlobalVector(theTrack.px(), theTrack.py(), theTrack.pz()), theFittedVertex);
0653
0654             double s_ip2dpv_corr = s_ip2dpv.second.value();
0655             double s_ip2dpv_err = s_ip2dpv.second.error();
0656
0657             // PV3D
0658             std::pair<bool, Measurement1D> s_ip3dpv = signedImpactParameter3D(
0659                 theTTrack, GlobalVector(theTrack.px(), theTrack.py(), theTrack.pz()), theFittedVertex);
0660
0661             double s_ip3dpv_corr = s_ip3dpv.second.value();
0662             double s_ip3dpv_err = s_ip3dpv.second.error();
0663
0664             // PV3D absolute
0665             std::pair<bool, Measurement1D> ip3dpv = absoluteImpactParameter3D(theTTrack, theFittedVertex);
0666             double ip3d_corr = ip3dpv.second.value();
0667             double ip3d_err = ip3dpv.second.error();
0668
0669             // with respect to any specified vertex, such as primary vertex
0670             TrajectoryStateClosestToPoint traj = (theTTrack).trajectoryStateClosestToPoint(vert);
0671
0672             GlobalPoint refPoint = traj.position();
0673             GlobalPoint cPToVtx = traj.theState().position();
0674
0675             float my_dx = refPoint.x() - myVertex.x();
0676             float my_dy = refPoint.y() - myVertex.y();
0677
0678             float my_dx2 = cPToVtx.x() - myVertex.x();
0679             float my_dy2 = cPToVtx.y() - myVertex.y();
0680
0681             float my_dxy = std::sqrt(my_dx * my_dx + my_dy * my_dy);
0682
0683             double d0 = traj.perigeeParameters().transverseImpactParameter();
0684             //double d0_error = traj.perigeeError().transverseImpactParameterError();
0685             double z0 = traj.perigeeParameters().longitudinalImpactParameter();
0686             double z0_error = traj.perigeeError().longitudinalImpactParameterError();
0687
0688             if (debug_) {
0689               edm::LogInfo("PrimaryVertexValidation")
0690                   << "my_dx:" << my_dx << " my_dy:" << my_dy << " my_dxy:" << my_dxy << " my_dx2:" << my_dx2
0691                   << " my_dy2:" << my_dy2 << " d0: " << d0 << " dxyFromVtx:" << dxyFromMyVertex << "\n"
0692                   << " ============================== "
0693                   << "\n"
0694                   << "diff1:" << std::abs(d0) - std::abs(my_dxy) << "\n"
0695                   << "diff2:" << std::abs(d0) - std::abs(dxyFromMyVertex) << "\n"
0696                   << "diff3:" << (my_dx - my_dx2) << " " << (my_dy - my_dy2) << "\n"
0697                   << std::endl;
0698             }
0699
0700             // define IPs
0701
0702             dxyFromMyVertex_[nTracks_] = dxyFromMyVertex;
0703             dxyErrorFromMyVertex_[nTracks_] = s_ip2dpv_err;
0704             IPTsigFromMyVertex_[nTracks_] = dxyFromMyVertex / s_ip2dpv_err;
0705
0706             dzFromMyVertex_[nTracks_] = dzFromMyVertex;
0707             dzErrorFromMyVertex_[nTracks_] = dz_err;
0708             IPLsigFromMyVertex_[nTracks_] = dzFromMyVertex / dz_err;
0709
0710             d3DFromMyVertex_[nTracks_] = ip3d_corr;
0711             d3DErrorFromMyVertex_[nTracks_] = ip3d_err;
0712             IP3DsigFromMyVertex_[nTracks_] = (ip3d_corr / ip3d_err);
0713
0714             // fill directly the histograms of residuals
0715
0716             float trackphi = (theTrack.phi()) * (180. / M_PI);
0717             float tracketa = theTrack.eta();
0718             float trackpt = theTrack.pt();
0719             float trackp = theTrack.p();
0720             float tracknhits = theTrack.numberOfValidHits();
0721
0722             // determine the module number and ladder
0723
0724             int ladder_num = -1.;
0725             int module_num = -1.;
0726             int L1BPixHitCount = 0;
0727
0728             for (auto const& hit : theTrack.recHits()) {
0729               const DetId& detId = hit->geographicalId();
0730               unsigned int subid = detId.subdetId();
0731
0732               if (hit->isValid() && (subid == PixelSubdetector::PixelBarrel)) {
0733                 int layer = tTopo->pxbLayer(detId);
0734                 if (layer == 1) {
0735                   const SiPixelRecHit* prechit = dynamic_cast<const SiPixelRecHit*>(
0736                       hit);  //to be used to get the associated cluster and the cluster probability
0737                   double clusterProbability = prechit->clusterProbability(0);
0738                   if (clusterProbability > 0) {
0739                     h_probeL1ClusterProb_->Fill(log10(clusterProbability));
0740                   }
0741
0742                   L1BPixHitCount += 1;
0743                   ladder_num = tTopo->pxbLadder(detId);
0744                   module_num = tTopo->pxbModule(detId);
0745                 }
0746               }
0747             }
0748
0749             h_probeL1Ladder_->Fill(ladder_num);
0750             h_probeL1Module_->Fill(module_num);
0751             h2_probeLayer1Map_->Fill(module_num, ladder_num);
0752             h_probeHasBPixL1Overlap_->Fill(L1BPixHitCount);
0753
0754             // residuals vs ladder and module number for map
0755             if (module_num > 0 && ladder_num > 0) {  // only if we are on BPix Layer 1
0756               a_dxyL1ResidualsMap[ladder_num - 1][module_num - 1]->Fill(dxyFromMyVertex * cmToum);
0757               a_dzL1ResidualsMap[ladder_num - 1][module_num - 1]->Fill(dzFromMyVertex * cmToum);
0758               n_dxyL1ResidualsMap[ladder_num - 1][module_num - 1]->Fill(dxyFromMyVertex / s_ip2dpv_err);
0759               n_dzL1ResidualsMap[ladder_num - 1][module_num - 1]->Fill(dzFromMyVertex / dz_err);
0760             }
0761
0762             // filling the pT-binned distributions
0763
0764             for (int ipTBin = 0; ipTBin < nPtBins_; ipTBin++) {
0765               float pTF = mypT_bins_[ipTBin];
0766               float pTL = mypT_bins_[ipTBin + 1];
0767
0768               if (debug_)
0769                 edm::LogInfo("PrimaryVertexValidation") << "ipTBin:" << ipTBin << " " << mypT_bins_[ipTBin]
0770                                                         << " < pT < " << mypT_bins_[ipTBin + 1] << std::endl;
0771
0772               if (std::abs(tracketa) < 1.5 && (trackpt >= pTF && trackpt < pTL)) {
0773                 if (debug_)
0774                   edm::LogInfo("PrimaryVertexValidation") << "passes this cut: " << mypT_bins_[ipTBin] << std::endl;
0775                 PVValHelper::fillByIndex(h_dxy_pT_, ipTBin, dxyFromMyVertex * cmToum, "9");
0776                 PVValHelper::fillByIndex(h_dz_pT_, ipTBin, dzFromMyVertex * cmToum, "10");
0777                 PVValHelper::fillByIndex(h_norm_dxy_pT_, ipTBin, dxyFromMyVertex / s_ip2dpv_err, "11");
0778                 PVValHelper::fillByIndex(h_norm_dz_pT_, ipTBin, dzFromMyVertex / dz_err, "12");
0779
0780                 if (std::abs(tracketa) < 1.) {
0781                   if (debug_)
0782                     edm::LogInfo("PrimaryVertexValidation")
0783                         << "passes tight eta cut: " << mypT_bins_[ipTBin] << std::endl;
0784                   PVValHelper::fillByIndex(h_dxy_Central_pT_, ipTBin, dxyFromMyVertex * cmToum, "13");
0785                   PVValHelper::fillByIndex(h_dz_Central_pT_, ipTBin, dzFromMyVertex * cmToum, "14");
0786                   PVValHelper::fillByIndex(h_norm_dxy_Central_pT_, ipTBin, dxyFromMyVertex / s_ip2dpv_err, "15");
0787                   PVValHelper::fillByIndex(h_norm_dz_Central_pT_, ipTBin, dzFromMyVertex / dz_err, "16");
0788                 }
0789               }
0790             }
0791
0792             // checks on the probe track quality
0793             if (trackpt >= ptOfProbe_ && std::abs(tracketa) <= etaOfProbe_ && tracknhits >= nHitsOfProbe_ &&
0794                 trackp >= pOfProbe_) {
0795               std::pair<bool, bool> pixelOcc = pixelHitsCheck((theTTrack));
0796
0797               if (debug_) {
0798                 if (pixelOcc.first == true)
0799                   edm::LogInfo("PrimaryVertexValidation") << "has BPIx hits" << std::endl;
0800                 if (pixelOcc.second == true)
0801                   edm::LogInfo("PrimaryVertexValidation") << "has FPix hits" << std::endl;
0802               }
0803
0804               if (!doBPix_ && (pixelOcc.first == true))
0805                 continue;
0806               if (!doFPix_ && (pixelOcc.second == true))
0807                 continue;
0808
0809               fillTrackHistos(hDA, "sel", &(theTTrack), vertex, beamSpot, fBfield_);
0810
0811               // probe checks
0812               h_probePt_->Fill(theTrack.pt());
0813               h_probePtRebin_->Fill(theTrack.pt());
0814               h_probeP_->Fill(theTrack.p());
0815               h_probeEta_->Fill(theTrack.eta());
0816               h_probePhi_->Fill(theTrack.phi());
0817               h2_probeEtaPhi_->Fill(theTrack.eta(), theTrack.phi());
0818               h2_probeEtaPt_->Fill(theTrack.eta(), theTrack.pt());
0819
0820               h_probeChi2_->Fill(theTrack.chi2());
0821               h_probeNormChi2_->Fill(theTrack.normalizedChi2());
0822               h_probeCharge_->Fill(theTrack.charge());
0823               h_probeQoverP_->Fill(theTrack.qoverp());
0824               h_probeHits_->Fill(theTrack.numberOfValidHits());
0825               h_probeHits1D_->Fill(nRecHit1D);
0826               h_probeHits2D_->Fill(nRecHit2D);
0827               h_probeHitsInTIB_->Fill(nhitinTIB);
0828               h_probeHitsInTOB_->Fill(nhitinTOB);
0829               h_probeHitsInTID_->Fill(nhitinTID);
0830               h_probeHitsInTEC_->Fill(nhitinTEC);
0831               h_probeHitsInBPIX_->Fill(nhitinBPIX);
0832               h_probeHitsInFPIX_->Fill(nhitinFPIX);
0833
0834               float dxyRecoV = theTrack.dz(theRecoVertex);
0835               float dzRecoV = theTrack.dxy(theRecoVertex);
0836               float dxysigmaRecoV =
0837                   TMath::Sqrt(theTrack.d0Error() * theTrack.d0Error() + xErrOfflineVertex_ * yErrOfflineVertex_);
0838               float dzsigmaRecoV =
0839                   TMath::Sqrt(theTrack.dzError() * theTrack.dzError() + zErrOfflineVertex_ * zErrOfflineVertex_);
0840
0841               double zTrack = (theTTrack.stateAtBeamLine().trackStateAtPCA()).position().z();
0842               double zVertex = theFittedVertex.position().z();
0843               double tantheta = tan((theTTrack.stateAtBeamLine().trackStateAtPCA()).momentum().theta());
0844
0845               double dz2 = pow(theTrack.dzError(), 2) + wxy2_ / pow(tantheta, 2);
0846               double restrkz = zTrack - zVertex;
0847               double pulltrkz = (zTrack - zVertex) / TMath::Sqrt(dz2);
0848
0849               h_probedxyRecoV_->Fill(dxyRecoV);
0850               h_probedzRecoV_->Fill(dzRecoV);
0851
0852               h_probedzRefitV_->Fill(dxyFromMyVertex);
0853               h_probedxyRefitV_->Fill(dzFromMyVertex);
0854
0855               h_probed0RefitV_->Fill(d0);
0856               h_probez0RefitV_->Fill(z0);
0857
0858               h_probesignIP2DRefitV_->Fill(s_ip2dpv_corr);
0859               h_probed3DRefitV_->Fill(ip3d_corr);
0860               h_probereszRefitV_->Fill(restrkz);
0861
0862               h_probeRecoVSigZ_->Fill(dzRecoV / dzsigmaRecoV);
0863               h_probeRecoVSigXY_->Fill(dxyRecoV / dxysigmaRecoV);
0864               h_probeRefitVSigZ_->Fill(dzFromMyVertex / dz_err);
0865               h_probeRefitVSigXY_->Fill(dxyFromMyVertex / s_ip2dpv_err);
0866               h_probeRefitVSig3D_->Fill(ip3d_corr / ip3d_err);
0867               h_probeRefitVLogSig3D_->Fill(log10(ip3d_corr / ip3d_err));
0868               h_probeRefitVSigResZ_->Fill(pulltrkz);
0869
0870               a_dxyVsPhi->Fill(trackphi, dxyFromMyVertex * cmToum);
0871               a_dzVsPhi->Fill(trackphi, z0 * cmToum);
0872               n_dxyVsPhi->Fill(trackphi, dxyFromMyVertex / s_ip2dpv_err);
0873               n_dzVsPhi->Fill(trackphi, z0 / z0_error);
0874
0875               a_dxyVsEta->Fill(tracketa, dxyFromMyVertex * cmToum);
0876               a_dzVsEta->Fill(tracketa, z0 * cmToum);
0877               n_dxyVsEta->Fill(tracketa, dxyFromMyVertex / s_ip2dpv_err);
0878               n_dzVsEta->Fill(tracketa, z0 / z0_error);
0879
0880               if (ladder_num > 0 && module_num > 0) {
0881                 LogDebug("PrimaryVertexValidation")
0882                     << " ladder_num: " << ladder_num << " module_num: " << module_num << std::endl;
0883
0884                 PVValHelper::fillByIndex(h_dxy_modZ_, module_num - 1, dxyFromMyVertex * cmToum, "17");
0885                 PVValHelper::fillByIndex(h_dz_modZ_, module_num - 1, dzFromMyVertex * cmToum, "18");
0886                 PVValHelper::fillByIndex(h_norm_dxy_modZ_, module_num - 1, dxyFromMyVertex / s_ip2dpv_err, "19");
0887                 PVValHelper::fillByIndex(h_norm_dz_modZ_, module_num - 1, dzFromMyVertex / dz_err, "20");
0888
0889                 PVValHelper::fillByIndex(h_dxy_ladder_, ladder_num - 1, dxyFromMyVertex * cmToum, "21");
0890
0891                 LogDebug("PrimaryVertexValidation") << "h_dxy_ladder size:" << h_dxy_ladder_.size() << std::endl;
0892
0893                 if (L1BPixHitCount == 1) {
0894                   PVValHelper::fillByIndex(h_dxy_ladderNoOverlap_, ladder_num - 1, dxyFromMyVertex * cmToum);
0895                 } else {
0896                   PVValHelper::fillByIndex(h_dxy_ladderOverlap_, ladder_num - 1, dxyFromMyVertex * cmToum);
0897                 }
0898
0899                 h2_probePassingLayer1Map_->Fill(module_num, ladder_num);
0900
0901                 PVValHelper::fillByIndex(h_dz_ladder_, ladder_num - 1, dzFromMyVertex * cmToum, "22");
0902                 PVValHelper::fillByIndex(h_norm_dxy_ladder_, ladder_num - 1, dxyFromMyVertex / s_ip2dpv_err, "23");
0903                 PVValHelper::fillByIndex(h_norm_dz_ladder_, ladder_num - 1, dzFromMyVertex / dz_err, "24");
0904               }
0905
0906               // filling the binned distributions
0907               for (int i = 0; i < nBins_; i++) {
0908                 float phiF = theDetails_.trendbins[PVValHelper::phi][i];
0909                 float phiL = theDetails_.trendbins[PVValHelper::phi][i + 1];
0910
0911                 float etaF = theDetails_.trendbins[PVValHelper::eta][i];
0912                 float etaL = theDetails_.trendbins[PVValHelper::eta][i + 1];
0913
0914                 if (tracketa >= etaF && tracketa < etaL) {
0915                   PVValHelper::fillByIndex(a_dxyEtaResiduals, i, dxyFromMyVertex * cmToum, "25");
0916                   PVValHelper::fillByIndex(a_dxEtaResiduals, i, my_dx * cmToum, "26");
0917                   PVValHelper::fillByIndex(a_dyEtaResiduals, i, my_dy * cmToum, "27");
0918                   PVValHelper::fillByIndex(a_dzEtaResiduals, i, dzFromMyVertex * cmToum, "28");
0919                   PVValHelper::fillByIndex(n_dxyEtaResiduals, i, dxyFromMyVertex / s_ip2dpv_err, "29");
0920                   PVValHelper::fillByIndex(n_dzEtaResiduals, i, dzFromMyVertex / dz_err, "30");
0921                   PVValHelper::fillByIndex(a_IP2DEtaResiduals, i, s_ip2dpv_corr * cmToum, "31");
0922                   PVValHelper::fillByIndex(n_IP2DEtaResiduals, i, s_ip2dpv_corr / s_ip2dpv_err, "32");
0923                   PVValHelper::fillByIndex(a_reszEtaResiduals, i, restrkz * cmToum, "33");
0924                   PVValHelper::fillByIndex(n_reszEtaResiduals, i, pulltrkz, "34");
0925                   PVValHelper::fillByIndex(a_d3DEtaResiduals, i, ip3d_corr * cmToum, "35");
0926                   PVValHelper::fillByIndex(n_d3DEtaResiduals, i, ip3d_corr / ip3d_err, "36");
0927                   PVValHelper::fillByIndex(a_IP3DEtaResiduals, i, s_ip3dpv_corr * cmToum, "37");
0928                   PVValHelper::fillByIndex(n_IP3DEtaResiduals, i, s_ip3dpv_corr / s_ip3dpv_err, "38");
0929                 }
0930
0931                 if (trackphi >= phiF && trackphi < phiL) {
0932                   PVValHelper::fillByIndex(a_dxyPhiResiduals, i, dxyFromMyVertex * cmToum, "39");
0933                   PVValHelper::fillByIndex(a_dxPhiResiduals, i, my_dx * cmToum, "40");
0934                   PVValHelper::fillByIndex(a_dyPhiResiduals, i, my_dy * cmToum, "41");
0935                   PVValHelper::fillByIndex(a_dzPhiResiduals, i, dzFromMyVertex * cmToum, "42");
0936                   PVValHelper::fillByIndex(n_dxyPhiResiduals, i, dxyFromMyVertex / s_ip2dpv_err, "43");
0937                   PVValHelper::fillByIndex(n_dzPhiResiduals, i, dzFromMyVertex / dz_err, "44");
0938                   PVValHelper::fillByIndex(a_IP2DPhiResiduals, i, s_ip2dpv_corr * cmToum, "45");
0939                   PVValHelper::fillByIndex(n_IP2DPhiResiduals, i, s_ip2dpv_corr / s_ip2dpv_err, "46");
0940                   PVValHelper::fillByIndex(a_reszPhiResiduals, i, restrkz * cmToum, "47");
0941                   PVValHelper::fillByIndex(n_reszPhiResiduals, i, pulltrkz, "48");
0942                   PVValHelper::fillByIndex(a_d3DPhiResiduals, i, ip3d_corr * cmToum, "49");
0943                   PVValHelper::fillByIndex(n_d3DPhiResiduals, i, ip3d_corr / ip3d_err, "50");
0944                   PVValHelper::fillByIndex(a_IP3DPhiResiduals, i, s_ip3dpv_corr * cmToum, "51");
0945                   PVValHelper::fillByIndex(n_IP3DPhiResiduals, i, s_ip3dpv_corr / s_ip3dpv_err, "52");
0946
0947                   for (int j = 0; j < nBins_; j++) {
0948                     float etaJ = theDetails_.trendbins[PVValHelper::eta][j];
0949                     float etaK = theDetails_.trendbins[PVValHelper::eta][j + 1];
0950
0951                     if (tracketa >= etaJ && tracketa < etaK) {
0952                       a_dxyResidualsMap[i][j]->Fill(dxyFromMyVertex * cmToum);
0953                       a_dzResidualsMap[i][j]->Fill(dzFromMyVertex * cmToum);
0954                       n_dxyResidualsMap[i][j]->Fill(dxyFromMyVertex / s_ip2dpv_err);
0955                       n_dzResidualsMap[i][j]->Fill(dzFromMyVertex / dz_err);
0956                       a_d3DResidualsMap[i][j]->Fill(ip3d_corr * cmToum);
0957                       n_d3DResidualsMap[i][j]->Fill(ip3d_corr / ip3d_err);
0958                     }
0959                   }
0960                 }
0961               }
0962             }
0963
0964             if (debug_) {
0965               edm::LogInfo("PrimaryVertexValidation")
0966                   << " myVertex.x()= " << myVertex.x() << "\n"
0967                   << " myVertex.y()= " << myVertex.y() << " \n"
0968                   << " myVertex.z()= " << myVertex.z() << " \n"
0969                   << " theTrack.dz(myVertex)= " << theTrack.dz(myVertex) << " \n"
0970                   << " zPCA -myVertex.z() = " << (theTrack.vertex().z() - myVertex.z());
0971
0972             }  // ends if debug_
0973           }    // ends if the fitted vertex is Valid
0974
0975           //delete theFitter;
0976
0977         } catch (cms::Exception& er) {
0978           LogTrace("PrimaryVertexValidation") << "caught std::exception " << er.what() << std::endl;
0979         }
0980
0981       }  //ends if theFinalTracks.size() > 2
0982
0983       else {
0984         if (debug_)
0985           edm::LogInfo("PrimaryVertexValidation") << "Not enough tracks to make a vertex.  Returns no vertex info";
0986       }
0987
0988       ++nTracks_;
0989       ++nTracksPerClus_;
0990
0991       if (debug_)
0992         edm::LogInfo("PrimaryVertexValidation") << "Track " << i << " : pT = " << theTrack.pt();
0993
0994     }  // for loop on tracks
0995   }    // for loop on track clusters
0996
0997   // Fill the TTree if needed
0998
0999   if (storeNtuple_) {
1000     rootTree_->Fill();
1001   }
1002 }
1003
1004 // ------------ method called to discriminate 1D from 2D hits  ------------
1005 bool PrimaryVertexValidation::isHit2D(const TrackingRecHit& hit, const PVValHelper::detectorPhase& thePhase) const {
1006   if (hit.dimension() < 2) {
1007     return false;  // some (muon...) stuff really has RecHit1D
1008   } else {
1009     const DetId detId(hit.geographicalId());
1010     if (detId.det() == DetId::Tracker) {
1011       if (detId.subdetId() == PixelSubdetector::PixelBarrel || detId.subdetId() == PixelSubdetector::PixelEndcap) {
1012         return true;                                 // pixel is always 2D
1013       } else if (thePhase != PVValHelper::phase2) {  // should be SiStrip now
1014         if (dynamic_cast<const SiStripRecHit2D*>(&hit))
1015           return false;  // normal hit
1016         else if (dynamic_cast<const SiStripMatchedRecHit2D*>(&hit))
1017           return true;  // matched is 2D
1018         else if (dynamic_cast<const ProjectedSiStripRecHit2D*>(&hit))
1019           return false;  // crazy hit...
1020         else {
1021           edm::LogError("UnknownType") << "@SUB=PrimaryVertexValidation::isHit2D"
1022                                        << "Tracker hit not in pixel and neither SiStripRecHit2D nor "
1023                                        << "SiStripMatchedRecHit2D nor ProjectedSiStripRecHit2D.";
1024           return false;
1025         }
1026       } else {
1027         return false;
1028       }
1029     } else {  // not tracker??
1030       edm::LogWarning("DetectorMismatch") << "@SUB=PrimaryVertexValidation::isHit2D"
1031                                           << "Hit not in tracker with 'official' dimension >=2.";
1032       return true;  // dimension() >= 2 so accept that...
1033     }
1034   }
1035   // never reached...
1036 }
1037
1038 // ------------ method to check the presence of pixel hits  ------------
1039 std::pair<bool, bool> PrimaryVertexValidation::pixelHitsCheck(const reco::TransientTrack& track) {
1040   bool hasBPixHits = false;
1041   bool hasFPixHits = false;
1042
1043   const reco::HitPattern& p = track.hitPattern();
1044   if (p.numberOfValidPixelEndcapHits() != 0) {
1045     hasFPixHits = true;
1046   }
1047   if (p.numberOfValidPixelBarrelHits() != 0) {
1048     hasBPixHits = true;
1049   }
1050
1051   return std::make_pair(hasBPixHits, hasFPixHits);
1052 }
1053
1054 // ------------ method to check the presence of pixel hits  ------------
1055 bool PrimaryVertexValidation::hasFirstLayerPixelHits(const reco::TransientTrack& track) {
1056   using namespace reco;
1057   const HitPattern& p = track.hitPattern();
1058   for (int i = 0; i < p.numberOfAllHits(HitPattern::TRACK_HITS); i++) {
1059     uint32_t pattern = p.getHitPattern(HitPattern::TRACK_HITS, i);
1060     if (p.pixelBarrelHitFilter(pattern) || p.pixelEndcapHitFilter(pattern)) {
1061       if (p.getLayer(pattern) == 1) {
1062         if (p.validHitFilter(pattern)) {
1063           return true;
1064         }
1065       }
1066     }
1067   }
1068   return false;
1069 }
1070
1071 // ------------ method called once each job before begining the event loop  ------------
1072 void PrimaryVertexValidation::beginJob() {
1073   edm::LogInfo("PrimaryVertexValidation") << "######################################\n"
1074                                           << "Begin Job \n"
1075                                           << "######################################";
1076
1077   // Define TTree for output
1078   Nevt_ = 0;
1079   if (compressionSettings_ > 0) {
1080     fs->file().SetCompressionSettings(compressionSettings_);
1081   }
1082
1083   //  rootFile_ = new TFile(filename_.c_str(),"recreate");
1084   rootTree_ = fs->make<TTree>("tree", "PV Validation tree");
1085
1086   // Track Paramters
1087
1088   if (lightNtupleSwitch_) {
1089     rootTree_->Branch("EventNumber", &EventNumber_, "EventNumber/i");
1090     rootTree_->Branch("RunNumber", &RunNumber_, "RunNumber/i");
1091     rootTree_->Branch("LuminosityBlockNumber", &LuminosityBlockNumber_, "LuminosityBlockNumber/i");
1092     rootTree_->Branch("nOfflineVertices", &nOfflineVertices_, "nOfflineVertices/I");
1093     rootTree_->Branch("nTracks", &nTracks_, "nTracks/I");
1094     rootTree_->Branch("phi", &phi_, "phi[nTracks]/D");
1095     rootTree_->Branch("eta", &eta_, "eta[nTracks]/D");
1096     rootTree_->Branch("pt", &pt_, "pt[nTracks]/D");
1097     rootTree_->Branch("dxyFromMyVertex", &dxyFromMyVertex_, "dxyFromMyVertex[nTracks]/D");
1098     rootTree_->Branch("dzFromMyVertex", &dzFromMyVertex_, "dzFromMyVertex[nTracks]/D");
1099     rootTree_->Branch("d3DFromMyVertex", &d3DFromMyVertex_, "d3DFromMyVertex[nTracks]/D");
1100     rootTree_->Branch("IPTsigFromMyVertex", &IPTsigFromMyVertex_, "IPTsigFromMyVertex_[nTracks]/D");
1101     rootTree_->Branch("IPLsigFromMyVertex", &IPLsigFromMyVertex_, "IPLsigFromMyVertex_[nTracks]/D");
1102     rootTree_->Branch("IP3DsigFromMyVertex", &IP3DsigFromMyVertex_, "IP3DsigFromMyVertex_[nTracks]/D");
1103     rootTree_->Branch("hasRecVertex", &hasRecVertex_, "hasRecVertex[nTracks]/I");
1104     rootTree_->Branch("isGoodTrack", &isGoodTrack_, "isGoodTrack[nTracks]/I");
1105     rootTree_->Branch("isHighPurity", &isHighPurity_, "isHighPurity_[nTracks]/I");
1106
1107   } else {
1108     rootTree_->Branch("nTracks", &nTracks_, "nTracks/I");
1109     rootTree_->Branch("nTracksPerClus", &nTracksPerClus_, "nTracksPerClus/I");
1110     rootTree_->Branch("nClus", &nClus_, "nClus/I");
1111     rootTree_->Branch("xOfflineVertex", &xOfflineVertex_, "xOfflineVertex/D");
1112     rootTree_->Branch("yOfflineVertex", &yOfflineVertex_, "yOfflineVertex/D");
1113     rootTree_->Branch("zOfflineVertex", &zOfflineVertex_, "zOfflineVertex/D");
1114     rootTree_->Branch("BSx0", &BSx0_, "BSx0/D");
1115     rootTree_->Branch("BSy0", &BSy0_, "BSy0/D");
1116     rootTree_->Branch("BSz0", &BSz0_, "BSz0/D");
1117     rootTree_->Branch("Beamsigmaz", &Beamsigmaz_, "Beamsigmaz/D");
1118     rootTree_->Branch("Beamdxdz", &Beamdxdz_, "Beamdxdz/D");
1119     rootTree_->Branch("BeamWidthX", &BeamWidthX_, "BeamWidthX/D");
1120     rootTree_->Branch("BeamWidthY", &BeamWidthY_, "BeamWidthY/D");
1121     rootTree_->Branch("pt", &pt_, "pt[nTracks]/D");
1122     rootTree_->Branch("p", &p_, "p[nTracks]/D");
1123     rootTree_->Branch("nhits", &nhits_, "nhits[nTracks]/I");
1124     rootTree_->Branch("nhits1D", &nhits1D_, "nhits1D[nTracks]/I");
1125     rootTree_->Branch("nhits2D", &nhits2D_, "nhits2D[nTracks]/I");
1126     rootTree_->Branch("nhitsBPIX", &nhitsBPIX_, "nhitsBPIX[nTracks]/I");
1127     rootTree_->Branch("nhitsFPIX", &nhitsFPIX_, "nhitsFPIX[nTracks]/I");
1128     rootTree_->Branch("nhitsTIB", &nhitsTIB_, "nhitsTIB[nTracks]/I");
1129     rootTree_->Branch("nhitsTID", &nhitsTID_, "nhitsTID[nTracks]/I");
1130     rootTree_->Branch("nhitsTOB", &nhitsTOB_, "nhitsTOB[nTracks]/I");
1131     rootTree_->Branch("nhitsTEC", &nhitsTEC_, "nhitsTEC[nTracks]/I");
1132     rootTree_->Branch("eta", &eta_, "eta[nTracks]/D");
1133     rootTree_->Branch("theta", &theta_, "theta[nTracks]/D");
1134     rootTree_->Branch("phi", &phi_, "phi[nTracks]/D");
1135     rootTree_->Branch("chi2", &chi2_, "chi2[nTracks]/D");
1136     rootTree_->Branch("chi2ndof", &chi2ndof_, "chi2ndof[nTracks]/D");
1137     rootTree_->Branch("charge", &charge_, "charge[nTracks]/I");
1138     rootTree_->Branch("qoverp", &qoverp_, "qoverp[nTracks]/D");
1139     rootTree_->Branch("dz", &dz_, "dz[nTracks]/D");
1140     rootTree_->Branch("dxy", &dxy_, "dxy[nTracks]/D");
1141     rootTree_->Branch("dzBs", &dzBs_, "dzBs[nTracks]/D");
1142     rootTree_->Branch("dxyBs", &dxyBs_, "dxyBs[nTracks]/D");
1143     rootTree_->Branch("xPCA", &xPCA_, "xPCA[nTracks]/D");
1144     rootTree_->Branch("yPCA", &yPCA_, "yPCA[nTracks]/D");
1145     rootTree_->Branch("zPCA", &zPCA_, "zPCA[nTracks]/D");
1146     rootTree_->Branch("xUnbiasedVertex", &xUnbiasedVertex_, "xUnbiasedVertex[nTracks]/D");
1147     rootTree_->Branch("yUnbiasedVertex", &yUnbiasedVertex_, "yUnbiasedVertex[nTracks]/D");
1148     rootTree_->Branch("zUnbiasedVertex", &zUnbiasedVertex_, "zUnbiasedVertex[nTracks]/D");
1149     rootTree_->Branch("chi2normUnbiasedVertex", &chi2normUnbiasedVertex_, "chi2normUnbiasedVertex[nTracks]/F");
1150     rootTree_->Branch("chi2UnbiasedVertex", &chi2UnbiasedVertex_, "chi2UnbiasedVertex[nTracks]/F");
1151     rootTree_->Branch("DOFUnbiasedVertex", &DOFUnbiasedVertex_, " DOFUnbiasedVertex[nTracks]/F");
1152     rootTree_->Branch("chi2ProbUnbiasedVertex", &chi2ProbUnbiasedVertex_, "chi2ProbUnbiasedVertex[nTracks]/F");
1153     rootTree_->Branch(
1154         "sumOfWeightsUnbiasedVertex", &sumOfWeightsUnbiasedVertex_, "sumOfWeightsUnbiasedVertex[nTracks]/F");
1155     rootTree_->Branch("tracksUsedForVertexing", &tracksUsedForVertexing_, "tracksUsedForVertexing[nTracks]/I");
1156     rootTree_->Branch("dxyFromMyVertex", &dxyFromMyVertex_, "dxyFromMyVertex[nTracks]/D");
1157     rootTree_->Branch("dzFromMyVertex", &dzFromMyVertex_, "dzFromMyVertex[nTracks]/D");
1158     rootTree_->Branch("dxyErrorFromMyVertex", &dxyErrorFromMyVertex_, "dxyErrorFromMyVertex_[nTracks]/D");
1159     rootTree_->Branch("dzErrorFromMyVertex", &dzErrorFromMyVertex_, "dzErrorFromMyVertex_[nTracks]/D");
1160     rootTree_->Branch("IPTsigFromMyVertex", &IPTsigFromMyVertex_, "IPTsigFromMyVertex_[nTracks]/D");
1161     rootTree_->Branch("IPLsigFromMyVertex", &IPLsigFromMyVertex_, "IPLsigFromMyVertex_[nTracks]/D");
1162     rootTree_->Branch("hasRecVertex", &hasRecVertex_, "hasRecVertex[nTracks]/I");
1163     rootTree_->Branch("isGoodTrack", &isGoodTrack_, "isGoodTrack[nTracks]/I");
1164   }
1165
1166   // event histograms
1167   TFileDirectory EventFeatures = fs->mkdir("EventFeatures");
1168
1169   TH1F::SetDefaultSumw2(kTRUE);
1170
1171   h_lumiFromConfig =
1172       EventFeatures.make<TH1F>("h_lumiFromConfig", "luminosity from config;;luminosity of present run", 1, -0.5, 0.5);
1173   h_lumiFromConfig->SetBinContent(1, intLumi_);
1174
1175   h_runFromConfig = EventFeatures.make<TH1I>("h_runFromConfig",
1176                                              "run number from config;;run number (from configuration)",
1177                                              runControlNumbers_.size(),
1178                                              0.,
1179                                              runControlNumbers_.size());
1180
1181   for (const auto& run : runControlNumbers_ | boost::adaptors::indexed(1)) {
1182     h_runFromConfig->SetBinContent(run.index(), run.value());
1183   }
1184
1185   h_runFromEvent =
1186       EventFeatures.make<TH1I>("h_runFromEvent", "run number from event;;run number (from event)", 1, -0.5, 0.5);
1187   h_nTracks =
1188       EventFeatures.make<TH1F>("h_nTracks", "number of tracks per event;n_{tracks}/event;n_{events}", 300, -0.5, 299.5);
1189   h_nClus =
1190       EventFeatures.make<TH1F>("h_nClus", "number of track clusters;n_{clusters}/event;n_{events}", 50, -0.5, 49.5);
1191   h_nOfflineVertices = EventFeatures.make<TH1F>(
1192       "h_nOfflineVertices", "number of offline reconstructed vertices;n_{vertices}/event;n_{events}", 50, -0.5, 49.5);
1193   h_runNumber = EventFeatures.make<TH1F>("h_runNumber", "run number;run number;n_{events}", 100000, 250000., 350000.);
1194   h_xOfflineVertex = EventFeatures.make<TH1F>(
1195       "h_xOfflineVertex", "x-coordinate of offline vertex;x_{vertex};n_{events}", 100, -0.1, 0.1);
1196   h_yOfflineVertex = EventFeatures.make<TH1F>(
1197       "h_yOfflineVertex", "y-coordinate of offline vertex;y_{vertex};n_{events}", 100, -0.1, 0.1);
1198   h_zOfflineVertex = EventFeatures.make<TH1F>(
1199       "h_zOfflineVertex", "z-coordinate of offline vertex;z_{vertex};n_{events}", 100, -30., 30.);
1200   h_xErrOfflineVertex = EventFeatures.make<TH1F>(
1201       "h_xErrOfflineVertex", "x-coordinate error of offline vertex;err_{x}^{vtx};n_{events}", 100, 0., 0.01);
1202   h_yErrOfflineVertex = EventFeatures.make<TH1F>(
1203       "h_yErrOfflineVertex", "y-coordinate error of offline vertex;err_{y}^{vtx};n_{events}", 100, 0., 0.01);
1204   h_zErrOfflineVertex = EventFeatures.make<TH1F>(
1205       "h_zErrOfflineVertex", "z-coordinate error of offline vertex;err_{z}^{vtx};n_{events}", 100, 0., 10.);
1206   h_BSx0 = EventFeatures.make<TH1F>("h_BSx0", "x-coordinate of reco beamspot;x^{BS}_{0};n_{events}", 100, -0.1, 0.1);
1207   h_BSy0 = EventFeatures.make<TH1F>("h_BSy0", "y-coordinate of reco beamspot;y^{BS}_{0};n_{events}", 100, -0.1, 0.1);
1208   h_BSz0 = EventFeatures.make<TH1F>("h_BSz0", "z-coordinate of reco beamspot;z^{BS}_{0};n_{events}", 100, -1., 1.);
1209   h_Beamsigmaz =
1210       EventFeatures.make<TH1F>("h_Beamsigmaz", "z-coordinate beam width;#sigma_{Z}^{beam};n_{events}", 100, 0., 1.);
1211   h_BeamWidthX =
1212       EventFeatures.make<TH1F>("h_BeamWidthX", "x-coordinate beam width;#sigma_{X}^{beam};n_{events}", 100, 0., 0.01);
1213   h_BeamWidthY =
1214       EventFeatures.make<TH1F>("h_BeamWidthY", "y-coordinate beam width;#sigma_{Y}^{beam};n_{events}", 100, 0., 0.01);
1215
1216   h_etaMax = EventFeatures.make<TH1F>("etaMax", "etaMax", 1, -0.5, 0.5);
1217   h_pTinfo = EventFeatures.make<TH1F>("pTinfo", "pTinfo", 3, -1.5, 1.5);
1218   h_pTinfo->GetXaxis()->SetBinLabel(1, "n. bins");
1219   h_pTinfo->GetXaxis()->SetBinLabel(2, "pT min");
1220   h_pTinfo->GetXaxis()->SetBinLabel(3, "pT max");
1221
1222   h_nbins = EventFeatures.make<TH1F>("nbins", "nbins", 1, -0.5, 0.5);
1223   h_nLadders = EventFeatures.make<TH1F>("nladders", "n. ladders", 1, -0.5, 0.5);
1224   h_nModZ = EventFeatures.make<TH1F>("nModZ", "n. modules along z", 1, -0.5, 0.5);
1225
1226   // probe track histograms
1227   TFileDirectory ProbeFeatures = fs->mkdir("ProbeTrackFeatures");
1228
1229   h_probePt_ = ProbeFeatures.make<TH1F>("h_probePt", "p_{T} of probe track;track p_{T} (GeV); tracks", 100, 0., 50.);
1230   h_probePtRebin_ = ProbeFeatures.make<TH1F>(
1231       "h_probePtRebin", "p_{T} of probe track;track p_{T} (GeV); tracks", mypT_bins_.size() - 1, mypT_bins_.data());
1232   h_probeP_ = ProbeFeatures.make<TH1F>("h_probeP", "momentum of probe track;track p (GeV); tracks", 100, 0., 100.);
1233   h_probeEta_ = ProbeFeatures.make<TH1F>(
1234       "h_probeEta", "#eta of the probe track;track #eta;tracks", 54, -etaOfProbe_, etaOfProbe_);
1235   h_probePhi_ = ProbeFeatures.make<TH1F>("h_probePhi", "#phi of probe track;track #phi (rad);tracks", 100, -3.15, 3.15);
1236
1237   h2_probeEtaPhi_ =
1238       ProbeFeatures.make<TH2F>("h2_probeEtaPhi",
1239                                "probe track #phi vs #eta;#eta of probe track;track #phi of probe track (rad); tracks",
1240                                54,
1241                                -etaOfProbe_,
1242                                etaOfProbe_,
1243                                100,
1244                                -M_PI,
1245                                M_PI);
1246   h2_probeEtaPt_ = ProbeFeatures.make<TH2F>("h2_probeEtaPt",
1247                                             "probe track p_{T} vs #eta;#eta of probe track;track p_{T} (GeV); tracks",
1248                                             54,
1249                                             -etaOfProbe_,
1250                                             etaOfProbe_,
1251                                             100,
1252                                             0.,
1253                                             50.);
1254
1255   h_probeChi2_ =
1256       ProbeFeatures.make<TH1F>("h_probeChi2", "#chi^{2} of probe track;track #chi^{2}; tracks", 100, 0., 100.);
1257   h_probeNormChi2_ = ProbeFeatures.make<TH1F>(
1258       "h_probeNormChi2", " normalized #chi^{2} of probe track;track #chi^{2}/ndof; tracks", 100, 0., 10.);
1259   h_probeCharge_ =
1260       ProbeFeatures.make<TH1F>("h_probeCharge", "charge of probe track;track charge Q;tracks", 3, -1.5, 1.5);
1261   h_probeQoverP_ =
1262       ProbeFeatures.make<TH1F>("h_probeQoverP", "q/p of probe track; track Q/p (GeV^{-1});tracks", 200, -1., 1.);
1263   h_probedzRecoV_ = ProbeFeatures.make<TH1F>(
1264       "h_probedzRecoV", "d_{z}(V_{offline}) of probe track;track d_{z}(V_{off}) (cm);tracks", 200, -1., 1.);
1265   h_probedxyRecoV_ = ProbeFeatures.make<TH1F>(
1266       "h_probedxyRecoV", "d_{xy}(V_{offline}) of probe track;track d_{xy}(V_{off}) (cm);tracks", 200, -1., 1.);
1267   h_probedzRefitV_ = ProbeFeatures.make<TH1F>(
1268       "h_probedzRefitV", "d_{z}(V_{refit}) of probe track;track d_{z}(V_{fit}) (cm);tracks", 200, -0.5, 0.5);
1269   h_probesignIP2DRefitV_ = ProbeFeatures.make<TH1F>(
1270       "h_probesignIPRefitV", "ip_{2D}(V_{refit}) of probe track;track ip_{2D}(V_{fit}) (cm);tracks", 200, -1., 1.);
1271   h_probedxyRefitV_ = ProbeFeatures.make<TH1F>(
1272       "h_probedxyRefitV", "d_{xy}(V_{refit}) of probe track;track d_{xy}(V_{fit}) (cm);tracks", 200, -0.5, 0.5);
1273
1274   h_probez0RefitV_ = ProbeFeatures.make<TH1F>(
1275       "h_probez0RefitV", "z_{0}(V_{refit}) of probe track;track z_{0}(V_{fit}) (cm);tracks", 200, -1., 1.);
1276   h_probed0RefitV_ = ProbeFeatures.make<TH1F>(
1277       "h_probed0RefitV", "d_{0}(V_{refit}) of probe track;track d_{0}(V_{fit}) (cm);tracks", 200, -1., 1.);
1278
1279   h_probed3DRefitV_ = ProbeFeatures.make<TH1F>(
1280       "h_probed3DRefitV", "d_{3D}(V_{refit}) of probe track;track d_{3D}(V_{fit}) (cm);tracks", 200, 0., 1.);
1281   h_probereszRefitV_ = ProbeFeatures.make<TH1F>(
1282       "h_probeReszRefitV", "z_{track} -z_{V_{refit}};track res_{z}(V_{refit}) (cm);tracks", 200, -1., 1.);
1283
1284   h_probeRecoVSigZ_ = ProbeFeatures.make<TH1F>(
1285       "h_probeRecoVSigZ", "Longitudinal DCA Significance (reco);d_{z}(V_{off})/#sigma_{dz};tracks", 100, -8, 8);
1286   h_probeRecoVSigXY_ = ProbeFeatures.make<TH1F>(
1287       "h_probeRecoVSigXY", "Transverse DCA Significance (reco);d_{xy}(V_{off})/#sigma_{dxy};tracks", 100, -8, 8);
1288   h_probeRefitVSigZ_ = ProbeFeatures.make<TH1F>(
1289       "h_probeRefitVSigZ", "Longitudinal DCA Significance (refit);d_{z}(V_{fit})/#sigma_{dz};tracks", 100, -8, 8);
1290   h_probeRefitVSigXY_ = ProbeFeatures.make<TH1F>(
1291       "h_probeRefitVSigXY", "Transverse DCA Significance (refit);d_{xy}(V_{fit})/#sigma_{dxy};tracks", 100, -8, 8);
1292   h_probeRefitVSig3D_ = ProbeFeatures.make<TH1F>(
1293       "h_probeRefitVSig3D", "3D DCA Significance (refit);d_{3D}/#sigma_{3D};tracks", 100, 0., 20.);
1294   h_probeRefitVLogSig3D_ =
1295       ProbeFeatures.make<TH1F>("h_probeRefitVLogSig3D",
1296                                "log_{10}(3D DCA-Significance) (refit);log_{10}(d_{3D}/#sigma_{3D});tracks",
1297                                100,
1298                                -5.,
1299                                4.);
1300   h_probeRefitVSigResZ_ = ProbeFeatures.make<TH1F>(
1301       "h_probeRefitVSigResZ",
1302       "Longitudinal residual significance (refit);(z_{track} -z_{V_{fit}})/#sigma_{res_{z}};tracks",
1303       100,
1304       -8,
1305       8);
1306
1307   h_probeHits_ = ProbeFeatures.make<TH1F>("h_probeNRechits", "N_{hits}     ;N_{hits}    ;tracks", 40, -0.5, 39.5);
1308   h_probeHits1D_ = ProbeFeatures.make<TH1F>("h_probeNRechits1D", "N_{hits} 1D  ;N_{hits} 1D ;tracks", 40, -0.5, 39.5);
1309   h_probeHits2D_ = ProbeFeatures.make<TH1F>("h_probeNRechits2D", "N_{hits} 2D  ;N_{hits} 2D ;tracks", 40, -0.5, 39.5);
1310   h_probeHitsInTIB_ =
1311       ProbeFeatures.make<TH1F>("h_probeNRechitsTIB", "N_{hits} TIB ;N_{hits} TIB;tracks", 40, -0.5, 39.5);
1312   h_probeHitsInTOB_ =
1313       ProbeFeatures.make<TH1F>("h_probeNRechitsTOB", "N_{hits} TOB ;N_{hits} TOB;tracks", 40, -0.5, 39.5);
1314   h_probeHitsInTID_ =
1315       ProbeFeatures.make<TH1F>("h_probeNRechitsTID", "N_{hits} TID ;N_{hits} TID;tracks", 40, -0.5, 39.5);
1316   h_probeHitsInTEC_ =
1317       ProbeFeatures.make<TH1F>("h_probeNRechitsTEC", "N_{hits} TEC ;N_{hits} TEC;tracks", 40, -0.5, 39.5);
1318   h_probeHitsInBPIX_ =
1319       ProbeFeatures.make<TH1F>("h_probeNRechitsBPIX", "N_{hits} BPIX;N_{hits} BPIX;tracks", 40, -0.5, 39.5);
1320   h_probeHitsInFPIX_ =
1321       ProbeFeatures.make<TH1F>("h_probeNRechitsFPIX", "N_{hits} FPIX;N_{hits} FPIX;tracks", 40, -0.5, 39.5);
1322
1323   h_probeL1Ladder_ = ProbeFeatures.make<TH1F>(
1324       "h_probeL1Ladder", "Ladder number (L1 hit); ladder number", nLadders_ + 2, -1.5, nLadders_ + 0.5);
1325   h_probeL1Module_ = ProbeFeatures.make<TH1F>(
1326       "h_probeL1Module", "Module number (L1 hit); module number", nModZ_ + 2, -1.5, nModZ_ + 0.5);
1327
1328   h2_probeLayer1Map_ = ProbeFeatures.make<TH2F>("h2_probeLayer1Map",
1329                                                 "Position in Layer 1 of first hit;module number;ladder number",
1330                                                 nModZ_,
1331                                                 0.5,
1332                                                 nModZ_ + 0.5,
1333                                                 nLadders_,
1334                                                 0.5,
1335                                                 nLadders_ + 0.5);
1336
1337   h2_probePassingLayer1Map_ = ProbeFeatures.make<TH2F>("h2_probePassingLayer1Map",
1338                                                        "Position in Layer 1 of first hit;module number;ladder number",
1339                                                        nModZ_,
1340                                                        0.5,
1341                                                        nModZ_ + 0.5,
1342                                                        nLadders_,
1343                                                        0.5,
1344                                                        nLadders_ + 0.5);
1345   h_probeHasBPixL1Overlap_ =
1346       ProbeFeatures.make<TH1I>("h_probeHasBPixL1Overlap", "n. hits in L1;n. L1-BPix hits;tracks", 5, -0.5, 4.5);
1347   h_probeL1ClusterProb_ = ProbeFeatures.make<TH1F>(
1348       "h_probeL1ClusterProb",
1349       "log_{10}(Cluster Probability) for Layer1 hits;log_{10}(cluster probability); n. Layer1 hits",
1350       100,
1351       -10.,
1352       0.);
1353
1354   // refit vertex features
1355   TFileDirectory RefitVertexFeatures = fs->mkdir("RefitVertexFeatures");
1356   h_fitVtxNtracks_ = RefitVertexFeatures.make<TH1F>(
1357       "h_fitVtxNtracks", "N_{trks} used in vertex fit;N^{fit}_{tracks};vertices", 100, -0.5, 99.5);
1358   h_fitVtxNdof_ = RefitVertexFeatures.make<TH1F>(
1359       "h_fitVtxNdof", "N_{DOF} of vertex fit;N_{DOF} of refit vertex;vertices", 100, -0.5, 99.5);
1360   h_fitVtxChi2_ = RefitVertexFeatures.make<TH1F>(
1361       "h_fitVtxChi2", "#chi^{2} of vertex fit;vertex #chi^{2};vertices", 100, -0.5, 99.5);
1362   h_fitVtxChi2ndf_ = RefitVertexFeatures.make<TH1F>(
1363       "h_fitVtxChi2ndf", "#chi^{2}/ndf of vertex fit;vertex #chi^{2}/ndf;vertices", 100, -0.5, 9.5);
1364   h_fitVtxChi2Prob_ = RefitVertexFeatures.make<TH1F>(
1365       "h_fitVtxChi2Prob", "Prob(#chi^{2},ndf) of vertex fit;Prob(#chi^{2},ndf);vertices", 40, 0., 1.);
1366   h_fitVtxTrackWeights_ = RefitVertexFeatures.make<TH1F>(
1367       "h_fitVtxTrackWeights", "track weights associated to track;track weights;tracks", 40, 0., 1.);
1368   h_fitVtxTrackAverageWeight_ = RefitVertexFeatures.make<TH1F>(
1369       "h_fitVtxTrackAverageWeight_", "average track weight per vertex;#LT track weight #GT;vertices", 40, 0., 1.);
1370
1371   if (useTracksFromRecoVtx_) {
1372     TFileDirectory RecoVertexFeatures = fs->mkdir("RecoVertexFeatures");
1373     h_recoVtxNtracks_ =
1374         RecoVertexFeatures.make<TH1F>("h_recoVtxNtracks", "N^{vtx}_{trks};N^{vtx}_{trks};vertices", 100, -0.5, 99.5);
1375     h_recoVtxChi2ndf_ =
1376         RecoVertexFeatures.make<TH1F>("h_recoVtxChi2ndf", "#chi^{2}/ndf vtx;#chi^{2}/ndf vtx;vertices", 10, -0.5, 9.5);
1377     h_recoVtxChi2Prob_ = RecoVertexFeatures.make<TH1F>(
1378         "h_recoVtxChi2Prob", "Prob(#chi^{2},ndf);Prob(#chi^{2},ndf);vertices", 40, 0., 1.);
1379     h_recoVtxSumPt_ =
1380         RecoVertexFeatures.make<TH1F>("h_recoVtxSumPt", "Sum(p^{trks}_{T});Sum(p^{trks}_{T});vertices", 100, 0., 200.);
1381   }
1382
1383   TFileDirectory DA = fs->mkdir("DA");
1384   //DA.cd();
1385   hDA = bookVertexHistograms(DA);
1386   //for(std::map<std::string,TH1*>::const_iterator hist=hDA.begin(); hist!=hDA.end(); hist++){
1387   //hist->second->SetDirectory(DA);
1388   // DA.make<TH1F>(hist->second);
1389   // }
1390
1391   // initialize the residuals histograms
1392
1393   const float dxymax_phi = theDetails_.getHigh(PVValHelper::dxy, PVValHelper::phi);
1394   const float dzmax_phi = theDetails_.getHigh(PVValHelper::dz, PVValHelper::eta);
1395   const float dxymax_eta = theDetails_.getHigh(PVValHelper::dxy, PVValHelper::phi);
1396   const float dzmax_eta = theDetails_.getHigh(PVValHelper::dz, PVValHelper::eta);
1397   //const float d3Dmax_phi = theDetails_.getHigh(PVValHelper::d3D,PVValHelper::phi);
1398   const float d3Dmax_eta = theDetails_.getHigh(PVValHelper::d3D, PVValHelper::eta);
1399
1400   ///////////////////////////////////////////////////////////////////
1401   //
1402   // Unbiased track-to-vertex residuals
1403   // The vertex is refit without the probe track
1404   //
1405   ///////////////////////////////////////////////////////////////////
1406
1407   //    _   _            _      _         ___        _    _           _
1408   //   /_\ | |__ ___ ___| |_  _| |_ ___  | _ \___ __(_)__| |_  _ __ _| |___
1409   //  / _ \| '_ (_-</ _ \ | || |  _/ -_) |   / -_|_-< / _` | || / _` | (_-<
1410   // /_/ \_\_.__/__/\___/_|\_,_|\__\___| |_|_\___/__/_\__,_|\_,_\__,_|_/__/
1411   //
1412
1413   TFileDirectory AbsTransPhiRes = fs->mkdir("Abs_Transv_Phi_Residuals");
1414   a_dxyPhiResiduals = bookResidualsHistogram(AbsTransPhiRes, nBins_, PVValHelper::dxy, PVValHelper::phi);
1415   a_dxPhiResiduals = bookResidualsHistogram(AbsTransPhiRes, nBins_, PVValHelper::dx, PVValHelper::phi);
1416   a_dyPhiResiduals = bookResidualsHistogram(AbsTransPhiRes, nBins_, PVValHelper::dy, PVValHelper::phi);
1417   a_IP2DPhiResiduals = bookResidualsHistogram(AbsTransPhiRes, nBins_, PVValHelper::IP2D, PVValHelper::phi);
1418
1419   TFileDirectory AbsTransEtaRes = fs->mkdir("Abs_Transv_Eta_Residuals");
1420   a_dxyEtaResiduals = bookResidualsHistogram(AbsTransEtaRes, nBins_, PVValHelper::dxy, PVValHelper::eta);
1421   a_dxEtaResiduals = bookResidualsHistogram(AbsTransEtaRes, nBins_, PVValHelper::dx, PVValHelper::eta);
1422   a_dyEtaResiduals = bookResidualsHistogram(AbsTransEtaRes, nBins_, PVValHelper::dy, PVValHelper::eta);
1423   a_IP2DEtaResiduals = bookResidualsHistogram(AbsTransEtaRes, nBins_, PVValHelper::IP2D, PVValHelper::eta);
1424
1425   TFileDirectory AbsLongPhiRes = fs->mkdir("Abs_Long_Phi_Residuals");
1426   a_dzPhiResiduals = bookResidualsHistogram(AbsLongPhiRes, nBins_, PVValHelper::dz, PVValHelper::phi);
1427   a_reszPhiResiduals = bookResidualsHistogram(AbsLongPhiRes, nBins_, PVValHelper::resz, PVValHelper::phi);
1428
1429   TFileDirectory AbsLongEtaRes = fs->mkdir("Abs_Long_Eta_Residuals");
1430   a_dzEtaResiduals = bookResidualsHistogram(AbsLongEtaRes, nBins_, PVValHelper::dz, PVValHelper::eta);
1431   a_reszEtaResiduals = bookResidualsHistogram(AbsLongEtaRes, nBins_, PVValHelper::resz, PVValHelper::eta);
1432
1433   TFileDirectory Abs3DPhiRes = fs->mkdir("Abs_3D_Phi_Residuals");
1434   a_d3DPhiResiduals = bookResidualsHistogram(Abs3DPhiRes, nBins_, PVValHelper::d3D, PVValHelper::phi);
1435   a_IP3DPhiResiduals = bookResidualsHistogram(Abs3DPhiRes, nBins_, PVValHelper::IP3D, PVValHelper::phi);
1436
1437   TFileDirectory Abs3DEtaRes = fs->mkdir("Abs_3D_Eta_Residuals");
1438   a_d3DEtaResiduals = bookResidualsHistogram(Abs3DEtaRes, nBins_, PVValHelper::d3D, PVValHelper::eta);
1439   a_IP3DEtaResiduals = bookResidualsHistogram(Abs3DEtaRes, nBins_, PVValHelper::IP3D, PVValHelper::eta);
1440
1441   TFileDirectory NormTransPhiRes = fs->mkdir("Norm_Transv_Phi_Residuals");
1442   n_dxyPhiResiduals = bookResidualsHistogram(NormTransPhiRes, nBins_, PVValHelper::norm_dxy, PVValHelper::phi, true);
1443   n_IP2DPhiResiduals = bookResidualsHistogram(NormTransPhiRes, nBins_, PVValHelper::norm_IP2D, PVValHelper::phi, true);
1444
1445   TFileDirectory NormTransEtaRes = fs->mkdir("Norm_Transv_Eta_Residuals");
1446   n_dxyEtaResiduals = bookResidualsHistogram(NormTransEtaRes, nBins_, PVValHelper::norm_dxy, PVValHelper::eta, true);
1447   n_IP2DEtaResiduals = bookResidualsHistogram(NormTransEtaRes, nBins_, PVValHelper::norm_IP2D, PVValHelper::eta, true);
1448
1449   TFileDirectory NormLongPhiRes = fs->mkdir("Norm_Long_Phi_Residuals");
1450   n_dzPhiResiduals = bookResidualsHistogram(NormLongPhiRes, nBins_, PVValHelper::norm_dz, PVValHelper::phi, true);
1451   n_reszPhiResiduals = bookResidualsHistogram(NormLongPhiRes, nBins_, PVValHelper::norm_resz, PVValHelper::phi, true);
1452
1453   TFileDirectory NormLongEtaRes = fs->mkdir("Norm_Long_Eta_Residuals");
1454   n_dzEtaResiduals = bookResidualsHistogram(NormLongEtaRes, nBins_, PVValHelper::norm_dz, PVValHelper::eta, true);
1455   n_reszEtaResiduals = bookResidualsHistogram(NormLongEtaRes, nBins_, PVValHelper::norm_resz, PVValHelper::eta, true);
1456
1457   TFileDirectory Norm3DPhiRes = fs->mkdir("Norm_3D_Phi_Residuals");
1458   n_d3DPhiResiduals = bookResidualsHistogram(Norm3DPhiRes, nBins_, PVValHelper::norm_d3D, PVValHelper::phi, true);
1459   n_IP3DPhiResiduals = bookResidualsHistogram(Norm3DPhiRes, nBins_, PVValHelper::norm_IP3D, PVValHelper::phi, true);
1460
1461   TFileDirectory Norm3DEtaRes = fs->mkdir("Norm_3D_Eta_Residuals");
1462   n_d3DEtaResiduals = bookResidualsHistogram(Norm3DEtaRes, nBins_, PVValHelper::norm_d3D, PVValHelper::eta, true);
1463   n_IP3DEtaResiduals = bookResidualsHistogram(Norm3DEtaRes, nBins_, PVValHelper::norm_IP3D, PVValHelper::eta, true);
1464
1465   TFileDirectory AbsDoubleDiffRes = fs->mkdir("Abs_DoubleDiffResiduals");
1466   TFileDirectory NormDoubleDiffRes = fs->mkdir("Norm_DoubleDiffResiduals");
1467
1468   TFileDirectory AbsL1Map = fs->mkdir("Abs_L1Residuals");
1469   TFileDirectory NormL1Map = fs->mkdir("Norm_L1Residuals");
1470
1471   // book residuals vs pT histograms
1472
1473   TFileDirectory AbsTranspTRes = fs->mkdir("Abs_Transv_pT_Residuals");
1474   h_dxy_pT_ = bookResidualsHistogram(AbsTranspTRes, nPtBins_, PVValHelper::dxy, PVValHelper::pT);
1475
1476   TFileDirectory AbsLongpTRes = fs->mkdir("Abs_Long_pT_Residuals");
1477   h_dz_pT_ = bookResidualsHistogram(AbsLongpTRes, nPtBins_, PVValHelper::dz, PVValHelper::pT);
1478
1479   TFileDirectory NormTranspTRes = fs->mkdir("Norm_Transv_pT_Residuals");
1480   h_norm_dxy_pT_ = bookResidualsHistogram(NormTranspTRes, nPtBins_, PVValHelper::norm_dxy, PVValHelper::pT, true);
1481
1482   TFileDirectory NormLongpTRes = fs->mkdir("Norm_Long_pT_Residuals");
1483   h_norm_dz_pT_ = bookResidualsHistogram(NormLongpTRes, nPtBins_, PVValHelper::norm_dz, PVValHelper::pT, true);
1484
1485   // book residuals vs pT histograms in central region (|eta|<1.0)
1486
1487   TFileDirectory AbsTranspTCentralRes = fs->mkdir("Abs_Transv_pTCentral_Residuals");
1488   h_dxy_Central_pT_ = bookResidualsHistogram(AbsTranspTCentralRes, nPtBins_, PVValHelper::dxy, PVValHelper::pTCentral);
1489
1490   TFileDirectory AbsLongpTCentralRes = fs->mkdir("Abs_Long_pTCentral_Residuals");
1491   h_dz_Central_pT_ = bookResidualsHistogram(AbsLongpTCentralRes, nPtBins_, PVValHelper::dz, PVValHelper::pTCentral);
1492
1493   TFileDirectory NormTranspTCentralRes = fs->mkdir("Norm_Transv_pTCentral_Residuals");
1494   h_norm_dxy_Central_pT_ =
1495       bookResidualsHistogram(NormTranspTCentralRes, nPtBins_, PVValHelper::norm_dxy, PVValHelper::pTCentral, true);
1496
1497   TFileDirectory NormLongpTCentralRes = fs->mkdir("Norm_Long_pTCentral_Residuals");
1498   h_norm_dz_Central_pT_ =
1499       bookResidualsHistogram(NormLongpTCentralRes, nPtBins_, PVValHelper::norm_dz, PVValHelper::pTCentral, true);
1500
1501   // book residuals vs module number
1502
1503   TFileDirectory AbsTransModZRes = fs->mkdir("Abs_Transv_modZ_Residuals");
1504   h_dxy_modZ_ = bookResidualsHistogram(AbsTransModZRes, nModZ_, PVValHelper::dxy, PVValHelper::modZ);
1505
1506   TFileDirectory AbsLongModZRes = fs->mkdir("Abs_Long_modZ_Residuals");
1507   h_dz_modZ_ = bookResidualsHistogram(AbsLongModZRes, nModZ_, PVValHelper::dz, PVValHelper::modZ);
1508
1509   //  _  _                    _ _           _   ___        _    _           _
1510   // | \| |___ _ _ _ __  __ _| (_)______ __| | | _ \___ __(_)__| |_  _ __ _| |___
1511   // | .` / _ \ '_| '  \/ _` | | |_ / -_) _` | |   / -_|_-< / _` | || / _` | (_-<
1512   // |_|\_\___/_| |_|_|_\__,_|_|_/__\___\__,_| |_|_\___/__/_\__,_|\_,_\__,_|_/__/
1513   //
1514
1515   TFileDirectory NormTransModZRes = fs->mkdir("Norm_Transv_modZ_Residuals");
1516   h_norm_dxy_modZ_ = bookResidualsHistogram(NormTransModZRes, nModZ_, PVValHelper::norm_dxy, PVValHelper::modZ, true);
1517
1518   TFileDirectory NormLongModZRes = fs->mkdir("Norm_Long_modZ_Residuals");
1519   h_norm_dz_modZ_ = bookResidualsHistogram(NormLongModZRes, nModZ_, PVValHelper::norm_dz, PVValHelper::modZ, true);
1520
1521   TFileDirectory AbsTransLadderRes = fs->mkdir("Abs_Transv_ladder_Residuals");
1522   h_dxy_ladder_ = bookResidualsHistogram(AbsTransLadderRes, nLadders_, PVValHelper::dxy, PVValHelper::ladder);
1523
1524   TFileDirectory AbsTransLadderResOverlap = fs->mkdir("Abs_Transv_ladderOverlap_Residuals");
1525   h_dxy_ladderOverlap_ =
1526       bookResidualsHistogram(AbsTransLadderResOverlap, nLadders_, PVValHelper::dxy, PVValHelper::ladder);
1527
1528   TFileDirectory AbsTransLadderResNoOverlap = fs->mkdir("Abs_Transv_ladderNoOverlap_Residuals");
1529   h_dxy_ladderNoOverlap_ =
1530       bookResidualsHistogram(AbsTransLadderResNoOverlap, nLadders_, PVValHelper::dxy, PVValHelper::ladder);
1531
1532   TFileDirectory AbsLongLadderRes = fs->mkdir("Abs_Long_ladder_Residuals");
1533   h_dz_ladder_ = bookResidualsHistogram(AbsLongLadderRes, nLadders_, PVValHelper::dz, PVValHelper::ladder);
1534
1535   TFileDirectory NormTransLadderRes = fs->mkdir("Norm_Transv_ladder_Residuals");
1536   h_norm_dxy_ladder_ =
1537       bookResidualsHistogram(NormTransLadderRes, nLadders_, PVValHelper::norm_dxy, PVValHelper::ladder, true);
1538
1539   TFileDirectory NormLongLadderRes = fs->mkdir("Norm_Long_ladder_Residuals");
1540   h_norm_dz_ladder_ =
1541       bookResidualsHistogram(NormLongLadderRes, nLadders_, PVValHelper::norm_dz, PVValHelper::ladder, true);
1542
1543   // book residuals as function of nLadders and nModules
1544
1545   for (unsigned int iLadder = 0; iLadder < nLadders_; iLadder++) {
1546     for (unsigned int iModule = 0; iModule < nModZ_; iModule++) {
1547       a_dxyL1ResidualsMap[iLadder][iModule] =
1548           AbsL1Map.make<TH1F>(Form("histo_dxy_ladder%i_module%i", iLadder, iModule),
1549                               Form("d_{xy} ladder=%i module=%i;d_{xy} [#mum];tracks", iLadder, iModule),
1550                               theDetails_.histobins,
1551                               -dzmax_eta,
1552                               dzmax_eta);
1553
1554       a_dzL1ResidualsMap[iLadder][iModule] =
1555           AbsL1Map.make<TH1F>(Form("histo_dz_ladder%i_module%i", iLadder, iModule),
1556                               Form("d_{z} ladder=%i module=%i;d_{z} [#mum];tracks", iLadder, iModule),
1557                               theDetails_.histobins,
1558                               -dzmax_eta,
1559                               dzmax_eta);
1560
1561       n_dxyL1ResidualsMap[iLadder][iModule] =
1562           NormL1Map.make<TH1F>(Form("histo_norm_dxy_ladder%i_module%i", iLadder, iModule),
1563                                Form("d_{xy} ladder=%i module=%i;d_{xy}/#sigma_{d_{xy}};tracks", iLadder, iModule),
1564                                theDetails_.histobins,
1565                                -dzmax_eta / 100,
1566                                dzmax_eta / 100);
1567
1568       n_dzL1ResidualsMap[iLadder][iModule] =
1569           NormL1Map.make<TH1F>(Form("histo_norm_dz_ladder%i_module%i", iLadder, iModule),
1570                                Form("d_{z} ladder=%i module=%i;d_{z}/#sigma_{d_{z}};tracks", iLadder, iModule),
1571                                theDetails_.histobins,
1572                                -dzmax_eta / 100,
1573                                dzmax_eta / 100);
1574     }
1575   }
1576
1577   // book residuals as function of phi and eta
1578
1579   for (int i = 0; i < nBins_; ++i) {
1580     float phiF = theDetails_.trendbins[PVValHelper::phi][i];
1581     float phiL = theDetails_.trendbins[PVValHelper::phi][i + 1];
1582
1583     //  ___           _    _     ___  _  __  __   ___        _    _           _
1584     // |   \ ___ _  _| |__| |___|   \(_)/ _|/ _| | _ \___ __(_)__| |_  _ __ _| |___
1585     // | |) / _ \ || | '_ \ / -_) |) | |  _|  _| |   / -_|_-< / _` | || / _` | (_-<
1586     // |___/\___/\_,_|_.__/_\___|___/|_|_| |_|   |_|_\___/__/_\__,_|\_,_\__,_|_/__/
1587
1588     for (int j = 0; j < nBins_; ++j) {
1589       float etaF = theDetails_.trendbins[PVValHelper::eta][j];
1590       float etaL = theDetails_.trendbins[PVValHelper::eta][j + 1];
1591
1592       a_dxyResidualsMap[i][j] = AbsDoubleDiffRes.make<TH1F>(
1593           Form("histo_dxy_eta_plot%i_phi_plot%i", i, j),
1594           Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{xy};tracks", etaF, etaL, phiF, phiL),
1595           theDetails_.histobins,
1596           -dzmax_eta,
1597           dzmax_eta);
1598
1599       a_dzResidualsMap[i][j] = AbsDoubleDiffRes.make<TH1F>(
1600           Form("histo_dz_eta_plot%i_phi_plot%i", i, j),
1601           Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{z};tracks", etaF, etaL, phiF, phiL),
1602           theDetails_.histobins,
1603           -dzmax_eta,
1604           dzmax_eta);
1605
1606       a_d3DResidualsMap[i][j] = AbsDoubleDiffRes.make<TH1F>(
1607           Form("histo_d3D_eta_plot%i_phi_plot%i", i, j),
1608           Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{3D};tracks", etaF, etaL, phiF, phiL),
1609           theDetails_.histobins,
1610           0.,
1611           d3Dmax_eta);
1612
1613       n_dxyResidualsMap[i][j] = NormDoubleDiffRes.make<TH1F>(
1614           Form("histo_norm_dxy_eta_plot%i_phi_plot%i", i, j),
1615           Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{xy}/#sigma_{d_{xy}};tracks",
1616                etaF,
1617                etaL,
1618                phiF,
1619                phiL),
1620           theDetails_.histobins,
1621           -dzmax_eta / 100,
1622           dzmax_eta / 100);
1623
1624       n_dzResidualsMap[i][j] = NormDoubleDiffRes.make<TH1F>(
1625           Form("histo_norm_dz_eta_plot%i_phi_plot%i", i, j),
1626           Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{z}/#sigma_{d_{z}};tracks",
1627                etaF,
1628                etaL,
1629                phiF,
1630                phiL),
1631           theDetails_.histobins,
1632           -dzmax_eta / 100,
1633           dzmax_eta / 100);
1634
1635       n_d3DResidualsMap[i][j] = NormDoubleDiffRes.make<TH1F>(
1636           Form("histo_norm_d3D_eta_plot%i_phi_plot%i", i, j),
1637           Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{3D}/#sigma_{d_{3D}};tracks",
1638                etaF,
1639                etaL,
1640                phiF,
1641                phiL),
1642           theDetails_.histobins,
1643           0.,
1644           d3Dmax_eta);
1645     }
1646   }
1647
1648   // declaration of the directories
1649
1650   TFileDirectory BiasVsParameter = fs->mkdir("BiasVsParameter");
1651
1652   a_dxyVsPhi = BiasVsParameter.make<TH2F>("h2_dxy_vs_phi",
1653                                           "d_{xy} vs track #phi;track #phi [rad];track d_{xy}(PV) [#mum]",
1654                                           nBins_,
1655                                           -M_PI,
1656                                           M_PI,
1657                                           theDetails_.histobins,
1658                                           -dxymax_phi,
1659                                           dxymax_phi);
1660
1661   a_dzVsPhi = BiasVsParameter.make<TH2F>("h2_dz_vs_phi",
1662                                          "d_{z} vs track #phi;track #phi [rad];track d_{z}(PV) [#mum]",
1663                                          nBins_,
1664                                          -M_PI,
1665                                          M_PI,
1666                                          theDetails_.histobins,
1667                                          -dzmax_phi,
1668                                          dzmax_phi);
1669
1670   n_dxyVsPhi = BiasVsParameter.make<TH2F>(
1671       "h2_n_dxy_vs_phi",
1672       "d_{xy}/#sigma_{d_{xy}} vs track #phi;track #phi [rad];track d_{xy}(PV)/#sigma_{d_{xy}}",
1673       nBins_,
1674       -M_PI,
1675       M_PI,
1676       theDetails_.histobins,
1677       -dxymax_phi / 100.,
1678       dxymax_phi / 100.);
1679
1680   n_dzVsPhi =
1681       BiasVsParameter.make<TH2F>("h2_n_dz_vs_phi",
1682                                  "d_{z}/#sigma_{d_{z}} vs track #phi;track #phi [rad];track d_{z}(PV)/#sigma_{d_{z}}",
1683                                  nBins_,
1684                                  -M_PI,
1685                                  M_PI,
1686                                  theDetails_.histobins,
1687                                  -dzmax_phi / 100.,
1688                                  dzmax_phi / 100.);
1689
1690   a_dxyVsEta = BiasVsParameter.make<TH2F>("h2_dxy_vs_eta",
1691                                           "d_{xy} vs track #eta;track #eta;track d_{xy}(PV) [#mum]",
1692                                           nBins_,
1693                                           -etaOfProbe_,
1694                                           etaOfProbe_,
1695                                           theDetails_.histobins,
1696                                           -dxymax_eta,
1697                                           dzmax_eta);
1698
1699   a_dzVsEta = BiasVsParameter.make<TH2F>("h2_dz_vs_eta",
1700                                          "d_{z} vs track #eta;track #eta;track d_{z}(PV) [#mum]",
1701                                          nBins_,
1702                                          -etaOfProbe_,
1703                                          etaOfProbe_,
1704                                          theDetails_.histobins,
1705                                          -dzmax_eta,
1706                                          dzmax_eta);
1707
1708   n_dxyVsEta =
1709       BiasVsParameter.make<TH2F>("h2_n_dxy_vs_eta",
1710                                  "d_{xy}/#sigma_{d_{xy}} vs track #eta;track #eta;track d_{xy}(PV)/#sigma_{d_{xy}}",
1711                                  nBins_,
1712                                  -etaOfProbe_,
1713                                  etaOfProbe_,
1714                                  theDetails_.histobins,
1715                                  -dxymax_eta / 100.,
1716                                  dxymax_eta / 100.);
1717
1718   n_dzVsEta = BiasVsParameter.make<TH2F>("h2_n_dz_vs_eta",
1719                                          "d_{z}/#sigma_{d_{z}} vs track #eta;track #eta;track d_{z}(PV)/#sigma_{d_{z}}",
1720                                          nBins_,
1721                                          -etaOfProbe_,
1722                                          etaOfProbe_,
1723                                          theDetails_.histobins,
1724                                          -dzmax_eta / 100.,
1725                                          dzmax_eta / 100.);
1726
1727   MeanTrendsDir = fs->mkdir("MeanTrends");
1728   WidthTrendsDir = fs->mkdir("WidthTrends");
1729   MedianTrendsDir = fs->mkdir("MedianTrends");
1730   MADTrendsDir = fs->mkdir("MADTrends");
1731
1732   Mean2DMapsDir = fs->mkdir("MeanMaps");
1733   Width2DMapsDir = fs->mkdir("WidthMaps");
1734
1735   double highedge = nBins_ - 0.5;
1736   double lowedge = -0.5;
1737
1738   // means and widths from the fit
1739
1740   a_dxyPhiMeanTrend =
1741       MeanTrendsDir.make<TH1F>("means_dxy_phi",
1742                                "#LT d_{xy} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{xy} #GT [#mum]",
1743                                nBins_,
1744                                lowedge,
1745                                highedge);
1746
1747   a_dxyPhiWidthTrend =
1748       WidthTrendsDir.make<TH1F>("widths_dxy_phi",
1749                                 "#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees];#sigma_{d_{xy}} [#mum]",
1750                                 nBins_,
1751                                 lowedge,
1752                                 highedge);
1753
1754   a_dzPhiMeanTrend =
1755       MeanTrendsDir.make<TH1F>("means_dz_phi",
1756                                "#LT d_{z} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{z} #GT [#mum]",
1757                                nBins_,
1758                                lowedge,
1759                                highedge);
1760
1761   a_dzPhiWidthTrend =
1762       WidthTrendsDir.make<TH1F>("widths_dz_phi",
1763                                 "#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];#sigma_{d_{z}} [#mum]",
1764                                 nBins_,
1765                                 lowedge,
1766                                 highedge);
1767
1768   a_dxyEtaMeanTrend = MeanTrendsDir.make<TH1F>(
1769       "means_dxy_eta", "#LT d_{xy} #GT vs #eta sector;#eta (sector);#LT d_{xy} #GT [#mum]", nBins_, lowedge, highedge);
1770
1771   a_dxyEtaWidthTrend = WidthTrendsDir.make<TH1F>("widths_dxy_eta",
1772                                                  "#sigma_{d_{xy}} vs #eta sector;#eta (sector);#sigma_{d_{xy}} [#mum]",
1773                                                  nBins_,
1774                                                  lowedge,
1775                                                  highedge);
1776
1777   a_dzEtaMeanTrend = MeanTrendsDir.make<TH1F>(
1778       "means_dz_eta", "#LT d_{z} #GT vs #eta sector;#eta (sector);#LT d_{z} #GT [#mum]", nBins_, lowedge, highedge);
1779
1780   a_dzEtaWidthTrend = WidthTrendsDir.make<TH1F>(
1781       "widths_dz_eta", "#sigma_{d_{xy}} vs #eta sector;#eta (sector);#sigma_{d_{z}} [#mum]", nBins_, lowedge, highedge);
1782
1783   n_dxyPhiMeanTrend = MeanTrendsDir.make<TH1F>(
1784       "norm_means_dxy_phi",
1785       "#LT d_{xy}/#sigma_{d_{xy}} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{xy}/#sigma_{d_{xy}} #GT",
1786       nBins_,
1787       lowedge,
1788       highedge);
1789
1790   n_dxyPhiWidthTrend = WidthTrendsDir.make<TH1F>(
1791       "norm_widths_dxy_phi",
1792       "width(d_{xy}/#sigma_{d_{xy}}) vs #phi sector;#varphi (sector) [degrees]; width(d_{xy}/#sigma_{d_{xy}})",
1793       nBins_,
1794       lowedge,
1795       highedge);
1796
1797   n_dzPhiMeanTrend = MeanTrendsDir.make<TH1F>(
1798       "norm_means_dz_phi",
1799       "#LT d_{z}/#sigma_{d_{z}} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{z}/#sigma_{d_{z}} #GT",
1800       nBins_,
1801       lowedge,
1802       highedge);
1803
1804   n_dzPhiWidthTrend = WidthTrendsDir.make<TH1F>(
1805       "norm_widths_dz_phi",
1806       "width(d_{z}/#sigma_{d_{z}}) vs #phi sector;#varphi (sector) [degrees];width(d_{z}/#sigma_{d_{z}})",
1807       nBins_,
1808       lowedge,
1809       highedge);
1810
1811   n_dxyEtaMeanTrend = MeanTrendsDir.make<TH1F>(
1812       "norm_means_dxy_eta",
1813       "#LT d_{xy}/#sigma_{d_{xy}} #GT vs #eta sector;#eta (sector);#LT d_{xy}/#sigma_{d_{z}} #GT",
1814       nBins_,
1815       lowedge,
1816       highedge);
1817
1818   n_dxyEtaWidthTrend = WidthTrendsDir.make<TH1F>(
1819       "norm_widths_dxy_eta",
1820       "width(d_{xy}/#sigma_{d_{xy}}) vs #eta sector;#eta (sector);width(d_{xy}/#sigma_{d_{z}})",
1821       nBins_,
1822       lowedge,
1823       highedge);
1824
1825   n_dzEtaMeanTrend =
1826       MeanTrendsDir.make<TH1F>("norm_means_dz_eta",
1827                                "#LT d_{z}/#sigma_{d_{z}} #GT vs #eta sector;#eta (sector);#LT d_{z}/#sigma_{d_{z}} #GT",
1828                                nBins_,
1829                                lowedge,
1830                                highedge);
1831
1832   n_dzEtaWidthTrend =
1833       WidthTrendsDir.make<TH1F>("norm_widths_dz_eta",
1834                                 "width(d_{z}/#sigma_{d_{z}}) vs #eta sector;#eta (sector);width(d_{z}/#sigma_{d_{z}})",
1835                                 nBins_,
1836                                 lowedge,
1837                                 highedge);
1838
1839   // means and widhts vs pT and pTCentral
1840
1841   a_dxypTMeanTrend = MeanTrendsDir.make<TH1F>("means_dxy_pT",
1842                                               "#LT d_{xy} #GT vs pT;p_{T} [GeV];#LT d_{xy} #GT [#mum]",
1843                                               mypT_bins_.size() - 1,
1844                                               mypT_bins_.data());
1845
1846   a_dxypTWidthTrend = WidthTrendsDir.make<TH1F>("widths_dxy_pT",
1847                                                 "#sigma_{d_{xy}} vs pT;p_{T} [GeV];#sigma_{d_{xy}} [#mum]",
1848                                                 mypT_bins_.size() - 1,
1849                                                 mypT_bins_.data());
1850
1851   a_dzpTMeanTrend = MeanTrendsDir.make<TH1F>(
1852       "means_dz_pT", "#LT d_{z} #GT vs pT;p_{T} [GeV];#LT d_{z} #GT [#mum]", mypT_bins_.size() - 1, mypT_bins_.data());
1853
1854   a_dzpTWidthTrend = WidthTrendsDir.make<TH1F>("widths_dz_pT",
1855                                                "#sigma_{d_{z}} vs pT;p_{T} [GeV];#sigma_{d_{z}} [#mum]",
1856                                                mypT_bins_.size() - 1,
1857                                                mypT_bins_.data());
1858
1859   n_dxypTMeanTrend =
1860       MeanTrendsDir.make<TH1F>("norm_means_dxy_pT",
1861                                "#LT d_{xy}/#sigma_{d_{xy}} #GT vs pT;p_{T} [GeV];#LT d_{xy}/#sigma_{d_{xy}} #GT",
1862                                mypT_bins_.size() - 1,
1863                                mypT_bins_.data());
1864
1865   n_dxypTWidthTrend =
1866       WidthTrendsDir.make<TH1F>("norm_widths_dxy_pT",
1867                                 "width(d_{xy}/#sigma_{d_{xy}}) vs pT;p_{T} [GeV]; width(d_{xy}/#sigma_{d_{xy}})",
1868                                 mypT_bins_.size() - 1,
1869                                 mypT_bins_.data());
1870
1871   n_dzpTMeanTrend =
1872       MeanTrendsDir.make<TH1F>("norm_means_dz_pT",
1873                                "#LT d_{z}/#sigma_{d_{z}} #GT vs pT;p_{T} [GeV];#LT d_{z}/#sigma_{d_{z}} #GT",
1874                                mypT_bins_.size() - 1,
1875                                mypT_bins_.data());
1876
1877   n_dzpTWidthTrend =
1878       WidthTrendsDir.make<TH1F>("norm_widths_dz_pT",
1879                                 "width(d_{z}/#sigma_{d_{z}}) vs pT;p_{T} [GeV];width(d_{z}/#sigma_{d_{z}})",
1880                                 mypT_bins_.size() - 1,
1881                                 mypT_bins_.data());
1882
1883   a_dxypTCentralMeanTrend =
1884       MeanTrendsDir.make<TH1F>("means_dxy_pTCentral",
1885                                "#LT d_{xy} #GT vs p_{T};p_{T}(|#eta|<1.) [GeV];#LT d_{xy} #GT [#mum]",
1886                                mypT_bins_.size() - 1,
1887                                mypT_bins_.data());
1888
1889   a_dxypTCentralWidthTrend =
1890       WidthTrendsDir.make<TH1F>("widths_dxy_pTCentral",
1891                                 "#sigma_{d_{xy}} vs p_{T};p_{T}(|#eta|<1.) [GeV];#sigma_{d_{xy}} [#mum]",
1892                                 mypT_bins_.size() - 1,
1893                                 mypT_bins_.data());
1894
1895   a_dzpTCentralMeanTrend =
1896       MeanTrendsDir.make<TH1F>("means_dz_pTCentral",
1897                                "#LT d_{z} #GT vs p_{T};p_{T}(|#eta|<1.) [GeV];#LT d_{z} #GT [#mum]",
1898                                mypT_bins_.size() - 1,
1899                                mypT_bins_.data());
1900
1901   a_dzpTCentralWidthTrend =
1902       WidthTrendsDir.make<TH1F>("widths_dz_pTCentral",
1903                                 "#sigma_{d_{z}} vs p_{T};p_{T}(|#eta|<1.) [GeV];#sigma_{d_{z}} [#mum]",
1904                                 mypT_bins_.size() - 1,
1905                                 mypT_bins_.data());
1906
1907   n_dxypTCentralMeanTrend = MeanTrendsDir.make<TH1F>(
1908       "norm_means_dxy_pTCentral",
1909       "#LT d_{xy}/#sigma_{d_{xy}} #GT vs p_{T};p_{T}(|#eta|<1.) [GeV];#LT d_{xy}/#sigma_{d_{z}} #GT",
1910       mypT_bins_.size() - 1,
1911       mypT_bins_.data());
1912
1913   n_dxypTCentralWidthTrend = WidthTrendsDir.make<TH1F>(
1914       "norm_widths_dxy_pTCentral",
1915       "width(d_{xy}/#sigma_{d_{xy}}) vs p_{T};p_{T}(|#eta|<1.) [GeV];width(d_{xy}/#sigma_{d_{z}})",
1916       mypT_bins_.size() - 1,
1917       mypT_bins_.data());
1918
1919   n_dzpTCentralMeanTrend = MeanTrendsDir.make<TH1F>(
1920       "norm_means_dz_pTCentral",
1921       "#LT d_{z}/#sigma_{d_{z}} #GT vs p_{T};p_{T}(|#eta|<1.) [GeV];#LT d_{z}/#sigma_{d_{z}} #GT",
1922       mypT_bins_.size() - 1,
1923       mypT_bins_.data());
1924
1925   n_dzpTCentralWidthTrend = WidthTrendsDir.make<TH1F>(
1926       "norm_widths_dz_pTCentral",
1927       "width(d_{z}/#sigma_{d_{z}}) vs p_{T};p_{T}(|#eta|<1.) [GeV];width(d_{z}/#sigma_{d_{z}})",
1928       mypT_bins_.size() - 1,
1929       mypT_bins_.data());
1930
1931   // 2D maps
1932
1933   a_dxyMeanMap = Mean2DMapsDir.make<TH2F>("means_dxy_map",
1934                                           "#LT d_{xy} #GT map;#eta (sector);#varphi (sector) [degrees]",
1935                                           nBins_,
1936                                           lowedge,
1937                                           highedge,
1938                                           nBins_,
1939                                           lowedge,
1940                                           highedge);
1941
1942   a_dzMeanMap = Mean2DMapsDir.make<TH2F>("means_dz_map",
1943                                          "#LT d_{z} #GT map;#eta (sector);#varphi (sector) [degrees]",
1944                                          nBins_,
1945                                          lowedge,
1946                                          highedge,
1947                                          nBins_,
1948                                          lowedge,
1949                                          highedge);
1950
1951   n_dxyMeanMap = Mean2DMapsDir.make<TH2F>("norm_means_dxy_map",
1952                                           "#LT d_{xy}/#sigma_{d_{xy}} #GT map;#eta (sector);#varphi (sector) [degrees]",
1953                                           nBins_,
1954                                           lowedge,
1955                                           highedge,
1956                                           nBins_,
1957                                           lowedge,
1958                                           highedge);
1959
1960   n_dzMeanMap = Mean2DMapsDir.make<TH2F>("norm_means_dz_map",
1961                                          "#LT d_{z}/#sigma_{d_{z}} #GT map;#eta (sector);#varphi (sector) [degrees]",
1962                                          nBins_,
1963                                          lowedge,
1964                                          highedge,
1965                                          nBins_,
1966                                          lowedge,
1967                                          highedge);
1968
1969   a_dxyWidthMap = Width2DMapsDir.make<TH2F>("widths_dxy_map",
1970                                             "#sigma_{d_{xy}} map;#eta (sector);#varphi (sector) [degrees]",
1971                                             nBins_,
1972                                             lowedge,
1973                                             highedge,
1974                                             nBins_,
1975                                             lowedge,
1976                                             highedge);
1977
1978   a_dzWidthMap = Width2DMapsDir.make<TH2F>("widths_dz_map",
1979                                            "#sigma_{d_{z}} map;#eta (sector);#varphi (sector) [degrees]",
1980                                            nBins_,
1981                                            lowedge,
1982                                            highedge,
1983                                            nBins_,
1984                                            lowedge,
1985                                            highedge);
1986
1987   n_dxyWidthMap =
1988       Width2DMapsDir.make<TH2F>("norm_widths_dxy_map",
1989                                 "width(d_{xy}/#sigma_{d_{xy}}) map;#eta (sector);#varphi (sector) [degrees]",
1990                                 nBins_,
1991                                 lowedge,
1992                                 highedge,
1993                                 nBins_,
1994                                 lowedge,
1995                                 highedge);
1996
1997   n_dzWidthMap = Width2DMapsDir.make<TH2F>("norm_widths_dz_map",
1998                                            "width(d_{z}/#sigma_{d_{z}}) map;#eta (sector);#varphi (sector) [degrees]",
1999                                            nBins_,
2000                                            lowedge,
2001                                            highedge,
2002                                            nBins_,
2003                                            lowedge,
2004                                            highedge);
2005
2006   // medians and MADs
2007
2008   a_dxyPhiMedianTrend =
2009       MedianTrendsDir.make<TH1F>("medians_dxy_phi",
2010                                  "Median of d_{xy} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{xy}) [#mum]",
2011                                  nBins_,
2012                                  lowedge,
2013                                  highedge);
2014
2015   a_dxyPhiMADTrend = MADTrendsDir.make<TH1F>(
2016       "MADs_dxy_phi",
2017       "Median absolute deviation of d_{xy} vs #phi sector;#varphi (sector) [degrees];MAD(d_{xy}) [#mum]",
2018       nBins_,
2019       lowedge,
2020       highedge);
2021
2022   a_dzPhiMedianTrend =
2023       MedianTrendsDir.make<TH1F>("medians_dz_phi",
2024                                  "Median of d_{z} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{z}) [#mum]",
2025                                  nBins_,
2026                                  lowedge,
2027                                  highedge);
2028
2029   a_dzPhiMADTrend = MADTrendsDir.make<TH1F>(
2030       "MADs_dz_phi",
2031       "Median absolute deviation of d_{z} vs #phi sector;#varphi (sector) [degrees];MAD(d_{z}) [#mum]",
2032       nBins_,
2033       lowedge,
2034       highedge);
2035
2036   a_dxyEtaMedianTrend =
2037       MedianTrendsDir.make<TH1F>("medians_dxy_eta",
2038                                  "Median of d_{xy} vs #eta sector;#eta (sector);#mu_{1/2}(d_{xy}) [#mum]",
2039                                  nBins_,
2040                                  lowedge,
2041                                  highedge);
2042
2043   a_dxyEtaMADTrend =
2044       MADTrendsDir.make<TH1F>("MADs_dxy_eta",
2045                               "Median absolute deviation of d_{xy} vs #eta sector;#eta (sector);MAD(d_{xy}) [#mum]",
2046                               nBins_,
2047                               lowedge,
2048                               highedge);
2049
2050   a_dzEtaMedianTrend =
2051       MedianTrendsDir.make<TH1F>("medians_dz_eta",
2052                                  "Median of d_{z} vs #eta sector;#eta (sector);#mu_{1/2}(d_{z}) [#mum]",
2053                                  nBins_,
2054                                  lowedge,
2055                                  highedge);
2056
2057   a_dzEtaMADTrend =
2058       MADTrendsDir.make<TH1F>("MADs_dz_eta",
2059                               "Median absolute deviation of d_{z} vs #eta sector;#eta (sector);MAD(d_{z}) [#mum]",
2060                               nBins_,
2061                               lowedge,
2062                               highedge);
2063
2064   n_dxyPhiMedianTrend = MedianTrendsDir.make<TH1F>(
2065       "norm_medians_dxy_phi",
2066       "Median of d_{xy}/#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{xy}/#sigma_{d_{xy}})",
2067       nBins_,
2068       lowedge,
2069       highedge);
2070
2071   n_dxyPhiMADTrend = MADTrendsDir.make<TH1F>("norm_MADs_dxy_phi",
2072                                              "Median absolute deviation of d_{xy}/#sigma_{d_{xy}} vs #phi "
2073                                              "sector;#varphi (sector) [degrees]; MAD(d_{xy}/#sigma_{d_{xy}})",
2074                                              nBins_,
2075                                              lowedge,
2076                                              highedge);
2077
2078   n_dzPhiMedianTrend = MedianTrendsDir.make<TH1F>(
2079       "norm_medians_dz_phi",
2080       "Median of d_{z}/#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{z}/#sigma_{d_{z}})",
2081       nBins_,
2082       lowedge,
2083       highedge);
2084
2085   n_dzPhiMADTrend = MADTrendsDir.make<TH1F>("norm_MADs_dz_phi",
2086                                             "Median absolute deviation of d_{z}/#sigma_{d_{z}} vs #phi sector;#varphi "
2087                                             "(sector) [degrees];MAD(d_{z}/#sigma_{d_{z}})",
2088                                             nBins_,
2089                                             lowedge,
2090                                             highedge);
2091
2092   n_dxyEtaMedianTrend = MedianTrendsDir.make<TH1F>(
2093       "norm_medians_dxy_eta",
2094       "Median of d_{xy}/#sigma_{d_{xy}} vs #eta sector;#eta (sector);#mu_{1/2}(d_{xy}/#sigma_{d_{z}})",
2095       nBins_,
2096       lowedge,
2097       highedge);
2098
2099   n_dxyEtaMADTrend = MADTrendsDir.make<TH1F>(
2100       "norm_MADs_dxy_eta",
2101       "Median absolute deviation of d_{xy}/#sigma_{d_{xy}} vs #eta sector;#eta (sector);MAD(d_{xy}/#sigma_{d_{z}})",
2102       nBins_,
2103       lowedge,
2104       highedge);
2105
2106   n_dzEtaMedianTrend = MedianTrendsDir.make<TH1F>(
2107       "norm_medians_dz_eta",
2108       "Median of d_{z}/#sigma_{d_{z}} vs #eta sector;#eta (sector);#mu_{1/2}(d_{z}/#sigma_{d_{z}})",
2109       nBins_,
2110       lowedge,
2111       highedge);
2112
2113   n_dzEtaMADTrend = MADTrendsDir.make<TH1F>(
2114       "norm_MADs_dz_eta",
2115       "Median absolute deviation of d_{z}/#sigma_{d_{z}} vs #eta sector;#eta (sector);MAD(d_{z}/#sigma_{d_{z}})",
2116       nBins_,
2117       lowedge,
2118       highedge);
2119
2120   ///////////////////////////////////////////////////////////////////
2121   //
2122   // plots of biased residuals
2123   // The vertex includes the probe track
2124   //
2125   ///////////////////////////////////////////////////////////////////
2126
2127   if (useTracksFromRecoVtx_) {
2128     TFileDirectory AbsTransPhiBiasRes = fs->mkdir("Abs_Transv_Phi_BiasResiduals");
2129     a_dxyPhiBiasResiduals = bookResidualsHistogram(AbsTransPhiBiasRes, nBins_, PVValHelper::dxy, PVValHelper::phi);
2130
2131     TFileDirectory AbsTransEtaBiasRes = fs->mkdir("Abs_Transv_Eta_BiasResiduals");
2132     a_dxyEtaBiasResiduals = bookResidualsHistogram(AbsTransEtaBiasRes, nBins_, PVValHelper::dxy, PVValHelper::eta);
2133
2134     TFileDirectory AbsLongPhiBiasRes = fs->mkdir("Abs_Long_Phi_BiasResiduals");
2135     a_dzPhiBiasResiduals = bookResidualsHistogram(AbsLongPhiBiasRes, nBins_, PVValHelper::dz, PVValHelper::phi);
2136
2137     TFileDirectory AbsLongEtaBiasRes = fs->mkdir("Abs_Long_Eta_BiasResiduals");
2138     a_dzEtaBiasResiduals = bookResidualsHistogram(AbsLongEtaBiasRes, nBins_, PVValHelper::dz, PVValHelper::eta);
2139
2140     TFileDirectory NormTransPhiBiasRes = fs->mkdir("Norm_Transv_Phi_BiasResiduals");
2141     n_dxyPhiBiasResiduals = bookResidualsHistogram(NormTransPhiBiasRes, nBins_, PVValHelper::dxy, PVValHelper::phi);
2142
2143     TFileDirectory NormTransEtaBiasRes = fs->mkdir("Norm_Transv_Eta_BiasResiduals");
2144     n_dxyEtaBiasResiduals = bookResidualsHistogram(NormTransEtaBiasRes, nBins_, PVValHelper::dxy, PVValHelper::eta);
2145
2146     TFileDirectory NormLongPhiBiasRes = fs->mkdir("Norm_Long_Phi_BiasResiduals");
2147     n_dzPhiBiasResiduals = bookResidualsHistogram(NormLongPhiBiasRes, nBins_, PVValHelper::dz, PVValHelper::phi);
2148
2149     TFileDirectory NormLongEtaBiasRes = fs->mkdir("Norm_Long_Eta_BiasResiduals");
2150     n_dzEtaBiasResiduals = bookResidualsHistogram(NormLongEtaBiasRes, nBins_, PVValHelper::dz, PVValHelper::eta);
2151
2152     TFileDirectory AbsDoubleDiffBiasRes = fs->mkdir("Abs_DoubleDiffBiasResiduals");
2153     TFileDirectory NormDoubleDiffBiasRes = fs->mkdir("Norm_DoubleDiffBiasResiduals");
2154
2155     for (int i = 0; i < nBins_; ++i) {
2156       float phiF = theDetails_.trendbins[PVValHelper::phi][i];
2157       float phiL = theDetails_.trendbins[PVValHelper::phi][i + 1];
2158
2159       for (int j = 0; j < nBins_; ++j) {
2160         float etaF = theDetails_.trendbins[PVValHelper::eta][j];
2161         float etaL = theDetails_.trendbins[PVValHelper::eta][j + 1];
2162
2163         a_dxyBiasResidualsMap[i][j] = AbsDoubleDiffBiasRes.make<TH1F>(
2164             Form("histo_dxy_eta_plot%i_phi_plot%i", i, j),
2165             Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{xy} [#mum];tracks", etaF, etaL, phiF, phiL),
2166             theDetails_.histobins,
2167             -dzmax_eta,
2168             dzmax_eta);
2169
2170         a_dzBiasResidualsMap[i][j] = AbsDoubleDiffBiasRes.make<TH1F>(
2171             Form("histo_dxy_eta_plot%i_phi_plot%i", i, j),
2172             Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{z} [#mum];tracks", etaF, etaL, phiF, phiL),
2173             theDetails_.histobins,
2174             -dzmax_eta,
2175             dzmax_eta);
2176
2177         n_dxyBiasResidualsMap[i][j] = NormDoubleDiffBiasRes.make<TH1F>(
2178             Form("histo_norm_dxy_eta_plot%i_phi_plot%i", i, j),
2179             Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{xy}/#sigma_{d_{xy}};tracks",
2180                  etaF,
2181                  etaL,
2182                  phiF,
2183                  phiL),
2184             theDetails_.histobins,
2185             -dzmax_eta / 100,
2186             dzmax_eta / 100);
2187
2188         n_dzBiasResidualsMap[i][j] = NormDoubleDiffBiasRes.make<TH1F>(
2189             Form("histo_norm_dxy_eta_plot%i_phi_plot%i", i, j),
2190             Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{z}/#sigma_{d_{z}};tracks",
2191                  etaF,
2192                  etaL,
2193                  phiF,
2194                  phiL),
2195             theDetails_.histobins,
2196             -dzmax_eta / 100,
2197             dzmax_eta / 100);
2198       }
2199     }
2200
2201     // declaration of the directories
2202
2203     TFileDirectory MeanBiasTrendsDir = fs->mkdir("MeanBiasTrends");
2204     TFileDirectory WidthBiasTrendsDir = fs->mkdir("WidthBiasTrends");
2205     TFileDirectory MedianBiasTrendsDir = fs->mkdir("MedianBiasTrends");
2206     TFileDirectory MADBiasTrendsDir = fs->mkdir("MADBiasTrends");
2207
2208     TFileDirectory Mean2DBiasMapsDir = fs->mkdir("MeanBiasMaps");
2209     TFileDirectory Width2DBiasMapsDir = fs->mkdir("WidthBiasMaps");
2210
2211     // means and widths from the fit
2212
2213     a_dxyPhiMeanBiasTrend =
2214         MeanBiasTrendsDir.make<TH1F>("means_dxy_phi",
2215                                      "#LT d_{xy} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{xy} #GT [#mum]",
2216                                      nBins_,
2217                                      lowedge,
2218                                      highedge);
2219
2220     a_dxyPhiWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>(
2221         "widths_dxy_phi",
2222         "#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees];#sigma_{d_{xy}} [#mum]",
2223         nBins_,
2224         lowedge,
2225         highedge);
2226
2227     a_dzPhiMeanBiasTrend =
2228         MeanBiasTrendsDir.make<TH1F>("means_dz_phi",
2229                                      "#LT d_{z} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{z} #GT [#mum]",
2230                                      nBins_,
2231                                      lowedge,
2232                                      highedge);
2233
2234     a_dzPhiWidthBiasTrend =
2235         WidthBiasTrendsDir.make<TH1F>("widths_dz_phi",
2236                                       "#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];#sigma_{d_{z}} [#mum]",
2237                                       nBins_,
2238                                       lowedge,
2239                                       highedge);
2240
2241     a_dxyEtaMeanBiasTrend = MeanBiasTrendsDir.make<TH1F>(
2242         "means_dxy_eta", "#LT d_{xy} #GT vs #eta sector;#eta (sector);#LT d_{xy} #GT [#mum]", nBins_, lowedge, highedge);
2243
2244     a_dxyEtaWidthBiasTrend =
2245         WidthBiasTrendsDir.make<TH1F>("widths_dxy_eta",
2246                                       "#sigma_{d_{xy}} vs #eta sector;#eta (sector);#sigma_{d_{xy}} [#mum]",
2247                                       nBins_,
2248                                       lowedge,
2249                                       highedge);
2250
2251     a_dzEtaMeanBiasTrend = MeanBiasTrendsDir.make<TH1F>(
2252         "means_dz_eta", "#LT d_{z} #GT vs #eta sector;#eta (sector);#LT d_{z} #GT [#mum]", nBins_, lowedge, highedge);
2253
2254     a_dzEtaWidthBiasTrend =
2255         WidthBiasTrendsDir.make<TH1F>("widths_dz_eta",
2256                                       "#sigma_{d_{xy}} vs #eta sector;#eta (sector);#sigma_{d_{z}} [#mum]",
2257                                       nBins_,
2258                                       lowedge,
2259                                       highedge);
2260
2261     n_dxyPhiMeanBiasTrend = MeanBiasTrendsDir.make<TH1F>(
2262         "norm_means_dxy_phi",
2263         "#LT d_{xy}/#sigma_{d_{xy}} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{xy}/#sigma_{d_{xy}} #GT",
2264         nBins_,
2265         lowedge,
2266         highedge);
2267
2268     n_dxyPhiWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>(
2269         "norm_widths_dxy_phi",
2270         "width(d_{xy}/#sigma_{d_{xy}}) vs #phi sector;#varphi (sector) [degrees]; width(d_{xy}/#sigma_{d_{xy}})",
2271         nBins_,
2272         lowedge,
2273         highedge);
2274
2275     n_dzPhiMeanBiasTrend = MeanBiasTrendsDir.make<TH1F>(
2276         "norm_means_dz_phi",
2277         "#LT d_{z}/#sigma_{d_{z}} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{z}/#sigma_{d_{z}} #GT",
2278         nBins_,
2279         lowedge,
2280         highedge);
2281
2282     n_dzPhiWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>(
2283         "norm_widths_dz_phi",
2284         "width(d_{z}/#sigma_{d_{z}}) vs #phi sector;#varphi (sector) [degrees];width(d_{z}/#sigma_{d_{z}})",
2285         nBins_,
2286         lowedge,
2287         highedge);
2288
2289     n_dxyEtaMeanBiasTrend = MeanBiasTrendsDir.make<TH1F>(
2290         "norm_means_dxy_eta",
2291         "#LT d_{xy}/#sigma_{d_{xy}} #GT vs #eta sector;#eta (sector);#LT d_{xy}/#sigma_{d_{z}} #GT",
2292         nBins_,
2293         lowedge,
2294         highedge);
2295
2296     n_dxyEtaWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>(
2297         "norm_widths_dxy_eta",
2298         "width(d_{xy}/#sigma_{d_{xy}}) vs #eta sector;#eta (sector);width(d_{xy}/#sigma_{d_{z}})",
2299         nBins_,
2300         lowedge,
2301         highedge);
2302
2303     n_dzEtaMeanBiasTrend = MeanBiasTrendsDir.make<TH1F>(
2304         "norm_means_dz_eta",
2305         "#LT d_{z}/#sigma_{d_{z}} #GT vs #eta sector;#eta (sector);#LT d_{z}/#sigma_{d_{z}} #GT",
2306         nBins_,
2307         lowedge,
2308         highedge);
2309
2310     n_dzEtaWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>(
2311         "norm_widths_dz_eta",
2312         "width(d_{z}/#sigma_{d_{z}}) vs #eta sector;#eta (sector);width(d_{z}/#sigma_{d_{z}})",
2313         nBins_,
2314         lowedge,
2315         highedge);
2316
2317     // 2D maps
2318
2319     a_dxyMeanBiasMap = Mean2DBiasMapsDir.make<TH2F>("means_dxy_map",
2320                                                     "#LT d_{xy} #GT map;#eta (sector);#varphi (sector) [degrees]",
2321                                                     nBins_,
2322                                                     lowedge,
2323                                                     highedge,
2324                                                     nBins_,
2325                                                     lowedge,
2326                                                     highedge);
2327
2328     a_dzMeanBiasMap = Mean2DBiasMapsDir.make<TH2F>("means_dz_map",
2329                                                    "#LT d_{z} #GT map;#eta (sector);#varphi (sector) [degrees]",
2330                                                    nBins_,
2331                                                    lowedge,
2332                                                    highedge,
2333                                                    nBins_,
2334                                                    lowedge,
2335                                                    highedge);
2336
2337     n_dxyMeanBiasMap =
2338         Mean2DBiasMapsDir.make<TH2F>("norm_means_dxy_map",
2339                                      "#LT d_{xy}/#sigma_{d_{xy}} #GT map;#eta (sector);#varphi (sector) [degrees]",
2340                                      nBins_,
2341                                      lowedge,
2342                                      highedge,
2343                                      nBins_,
2344                                      lowedge,
2345                                      highedge);
2346
2347     n_dzMeanBiasMap =
2348         Mean2DBiasMapsDir.make<TH2F>("norm_means_dz_map",
2349                                      "#LT d_{z}/#sigma_{d_{z}} #GT map;#eta (sector);#varphi (sector) [degrees]",
2350                                      nBins_,
2351                                      lowedge,
2352                                      highedge,
2353                                      nBins_,
2354                                      lowedge,
2355                                      highedge);
2356
2357     a_dxyWidthBiasMap = Width2DBiasMapsDir.make<TH2F>("widths_dxy_map",
2358                                                       "#sigma_{d_{xy}} map;#eta (sector);#varphi (sector) [degrees]",
2359                                                       nBins_,
2360                                                       lowedge,
2361                                                       highedge,
2362                                                       nBins_,
2363                                                       lowedge,
2364                                                       highedge);
2365
2366     a_dzWidthBiasMap = Width2DBiasMapsDir.make<TH2F>("widths_dz_map",
2367                                                      "#sigma_{d_{z}} map;#eta (sector);#varphi (sector) [degrees]",
2368                                                      nBins_,
2369                                                      lowedge,
2370                                                      highedge,
2371                                                      nBins_,
2372                                                      lowedge,
2373                                                      highedge);
2374
2375     n_dxyWidthBiasMap =
2376         Width2DBiasMapsDir.make<TH2F>("norm_widths_dxy_map",
2377                                       "width(d_{xy}/#sigma_{d_{xy}}) map;#eta (sector);#varphi (sector) [degrees]",
2378                                       nBins_,
2379                                       lowedge,
2380                                       highedge,
2381                                       nBins_,
2382                                       lowedge,
2383                                       highedge);
2384
2385     n_dzWidthBiasMap =
2386         Width2DBiasMapsDir.make<TH2F>("norm_widths_dz_map",
2387                                       "width(d_{z}/#sigma_{d_{z}}) map;#eta (sector);#varphi (sector) [degrees]",
2388                                       nBins_,
2389                                       lowedge,
2390                                       highedge,
2391                                       nBins_,
2392                                       lowedge,
2393                                       highedge);
2394
2395     // medians and MADs
2396
2397     a_dxyPhiMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>(
2398         "medians_dxy_phi",
2399         "Median of d_{xy} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{xy}) [#mum]",
2400         nBins_,
2401         lowedge,
2402         highedge);
2403
2404     a_dxyPhiMADBiasTrend = MADBiasTrendsDir.make<TH1F>(
2405         "MADs_dxy_phi",
2406         "Median absolute deviation of d_{xy} vs #phi sector;#varphi (sector) [degrees];MAD(d_{xy}) [#mum]",
2407         nBins_,
2408         lowedge,
2409         highedge);
2410
2411     a_dzPhiMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>(
2412         "medians_dz_phi",
2413         "Median of d_{z} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{z}) [#mum]",
2414         nBins_,
2415         lowedge,
2416         highedge);
2417
2418     a_dzPhiMADBiasTrend = MADBiasTrendsDir.make<TH1F>(
2419         "MADs_dz_phi",
2420         "Median absolute deviation of d_{z} vs #phi sector;#varphi (sector) [degrees];MAD(d_{z}) [#mum]",
2421         nBins_,
2422         lowedge,
2423         highedge);
2424
2425     a_dxyEtaMedianBiasTrend =
2426         MedianBiasTrendsDir.make<TH1F>("medians_dxy_eta",
2427                                        "Median of d_{xy} vs #eta sector;#eta (sector);#mu_{1/2}(d_{xy}) [#mum]",
2428                                        nBins_,
2429                                        lowedge,
2430                                        highedge);
2431
2432     a_dxyEtaMADBiasTrend = MADBiasTrendsDir.make<TH1F>(
2433         "MADs_dxy_eta",
2434         "Median absolute deviation of d_{xy} vs #eta sector;#eta (sector);MAD(d_{xy}) [#mum]",
2435         nBins_,
2436         lowedge,
2437         highedge);
2438
2439     a_dzEtaMedianBiasTrend =
2440         MedianBiasTrendsDir.make<TH1F>("medians_dz_eta",
2441                                        "Median of d_{z} vs #eta sector;#eta (sector);#mu_{1/2}(d_{z}) [#mum]",
2442                                        nBins_,
2443                                        lowedge,
2444                                        highedge);
2445
2446     a_dzEtaMADBiasTrend =
2447         MADBiasTrendsDir.make<TH1F>("MADs_dz_eta",
2448                                     "Median absolute deviation of d_{z} vs #eta sector;#eta (sector);MAD(d_{z}) [#mum]",
2449                                     nBins_,
2450                                     lowedge,
2451                                     highedge);
2452
2453     n_dxyPhiMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>(
2454         "norm_medians_dxy_phi",
2455         "Median of d_{xy}/#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{xy}/#sigma_{d_{xy}})",
2456         nBins_,
2457         lowedge,
2458         highedge);
2459
2460     n_dxyPhiMADBiasTrend = MADBiasTrendsDir.make<TH1F>("norm_MADs_dxy_phi",
2461                                                        "Median absolute deviation of d_{xy}/#sigma_{d_{xy}} vs #phi "
2462                                                        "sector;#varphi (sector) [degrees]; MAD(d_{xy}/#sigma_{d_{xy}})",
2463                                                        nBins_,
2464                                                        lowedge,
2465                                                        highedge);
2466
2467     n_dzPhiMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>(
2468         "norm_medians_dz_phi",
2469         "Median of d_{z}/#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{z}/#sigma_{d_{z}})",
2470         nBins_,
2471         lowedge,
2472         highedge);
2473
2474     n_dzPhiMADBiasTrend = MADBiasTrendsDir.make<TH1F>("norm_MADs_dz_phi",
2475                                                       "Median absolute deviation of d_{z}/#sigma_{d_{z}} vs #phi "
2476                                                       "sector;#varphi (sector) [degrees];MAD(d_{z}/#sigma_{d_{z}})",
2477                                                       nBins_,
2478                                                       lowedge,
2479                                                       highedge);
2480
2481     n_dxyEtaMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>(
2482         "norm_medians_dxy_eta",
2483         "Median of d_{xy}/#sigma_{d_{xy}} vs #eta sector;#eta (sector);#mu_{1/2}(d_{xy}/#sigma_{d_{z}})",
2484         nBins_,
2485         lowedge,
2486         highedge);
2487
2488     n_dxyEtaMADBiasTrend = MADBiasTrendsDir.make<TH1F>(
2489         "norm_MADs_dxy_eta",
2490         "Median absolute deviation of d_{xy}/#sigma_{d_{xy}} vs #eta sector;#eta (sector);MAD(d_{xy}/#sigma_{d_{z}})",
2491         nBins_,
2492         lowedge,
2493         highedge);
2494
2495     n_dzEtaMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>(
2496         "norm_medians_dz_eta",
2497         "Median of d_{z}/#sigma_{d_{z}} vs #eta sector;#eta (sector);#mu_{1/2}(d_{z}/#sigma_{d_{z}})",
2498         nBins_,
2499         lowedge,
2500         highedge);
2501
2502     n_dzEtaMADBiasTrend = MADBiasTrendsDir.make<TH1F>(
2503         "norm_MADs_dz_eta",
2504         "Median absolute deviation of d_{z}/#sigma_{d_{z}} vs #eta sector;#eta (sector);MAD(d_{z}/#sigma_{d_{z}})",
2505         nBins_,
2506         lowedge,
2507         highedge);
2508   }
2509 }
2510
2511 // ------------ method called once every run before doing the event loop ----------------
2512 void PrimaryVertexValidation::beginRun(edm::Run const& iRun, edm::EventSetup const& iSetup) {
2513   //=======================================================
2514   // check on the b-field consistency
2515   //=======================================================
2516
2517   if (!isBFieldConsistentWithMode(iSetup)) {
2518     edm::LogWarning("PrimaryVertexValidation")
2519         << "*********************************************************************************\n"
2520         << "* The configuration (ptOfProbe > " << ptOfProbe_
2521         << "GeV) is not correctly set for current value of magnetic field \n"
2522         << "* Switching it to 0. !!! \n"
2523         << "*********************************************************************************" << std::endl;
2524     ptOfProbe_ = 0.;
2525   }
2526
2527   //=======================================================
2528   // Run numbers analysis
2529   //=======================================================
2530
2531   const auto runNumber_ = iRun.run();
2532   h_runNumber->Fill(runNumber_);
2533
2534   if (!runNumbersTimesLog_.count(runNumber_)) {
2535     auto times = getRunTime(iSetup);
2536
2537     if (debug_) {
2538       const time_t start_time = times.first / 1000000;
2539       edm::LogInfo("PrimaryVertexValidation")
2540           << runNumber_ << " has start time: " << times.first << " - " << times.second << std::endl;
2541       edm::LogInfo("PrimaryVertexValidation")
2542           << "human readable time: " << std::asctime(std::gmtime(&start_time)) << std::endl;
2543     }
2544
2545     runNumbersTimesLog_[runNumber_] = times;
2546   }
2547
2548   if (h_runFromEvent->GetEntries() == 0) {
2549     h_runFromEvent->SetBinContent(1, RunNumber_);
2550   }
2551
2552   //=======================================================
2553   // Retrieve tracker topology from geometry
2554   //=======================================================
2555
2556   const TrackerTopology* const tTopo = &iSetup.getData(topoTokenBR_);
2557
2558   //=======================================================
2559   // Retrieve geometry information
2560   //=======================================================
2561
2562   edm::LogInfo("read tracker geometry...");
2563   edm::ESHandle<TrackerGeometry> pDD = iSetup.getHandle(geomTokenBR_);
2564   edm::LogInfo("tracker geometry read") << "There are: " << pDD->dets().size() << " detectors";
2565
2566   PixelTopologyMap PTMap = PixelTopologyMap(pDD.product(), tTopo);
2567   nLadders_ = PTMap.getPXBLadders(1);
2568   nModZ_ = PTMap.getPXBModules(1);
2569
2570   //=======================================================
2571   // shrink to fit
2572   //=======================================================
2573
2574   if (h_dxy_ladderOverlap_.size() != nLadders_) {
2575     PVValHelper::shrinkHistVectorToFit(h_dxy_ladderOverlap_, nLadders_);
2576     PVValHelper::shrinkHistVectorToFit(h_dxy_ladderNoOverlap_, nLadders_);
2577     PVValHelper::shrinkHistVectorToFit(h_dxy_ladder_, nLadders_);
2578     PVValHelper::shrinkHistVectorToFit(h_dz_ladder_, nLadders_);
2579     PVValHelper::shrinkHistVectorToFit(h_norm_dxy_ladder_, nLadders_);
2580     PVValHelper::shrinkHistVectorToFit(h_norm_dz_ladder_, nLadders_);
2581
2582     if (debug_) {
2583       edm::LogInfo("PrimaryVertexValidation") << "checking size:" << h_dxy_ladder_.size() << std::endl;
2584     }
2585   }
2586
2587   if (h_dxy_modZ_.size() != nModZ_) {
2588     PVValHelper::shrinkHistVectorToFit(h_dxy_modZ_, nModZ_);
2589     PVValHelper::shrinkHistVectorToFit(h_dz_modZ_, nModZ_);
2590     PVValHelper::shrinkHistVectorToFit(h_norm_dxy_modZ_, nModZ_);
2591     PVValHelper::shrinkHistVectorToFit(h_norm_dxy_modZ_, nModZ_);
2592
2593     if (debug_) {
2594       edm::LogInfo("PrimaryVertexValidation") << "checking size:" << h_dxy_modZ_.size() << std::endl;
2595     }
2596   }
2597
2598   if ((pDD->isThere(GeomDetEnumerators::P2PXB)) || (pDD->isThere(GeomDetEnumerators::P2PXEC))) {
2599     // switch on the phase-2
2600     phase_ = PVValHelper::phase2;
2601     if (debug_) {
2602       edm::LogInfo("PrimaryVertexValidation")
2603           << " pixel phase2 setup, nLadders: " << nLadders_ << " nModules:" << nModZ_;
2604     }
2605
2606   } else if ((pDD->isThere(GeomDetEnumerators::P1PXB)) || (pDD->isThere(GeomDetEnumerators::P1PXEC))) {
2607     // switch on the phase-1
2608     phase_ = PVValHelper::phase1;
2609     if (debug_) {
2610       edm::LogInfo("PrimaryVertexValidation")
2611           << " pixel phase1 setup, nLadders: " << nLadders_ << " nModules:" << nModZ_;
2612     }
2613   } else {
2614     // switch on the phase-0
2615     phase_ = PVValHelper::phase0;
2616     if (debug_) {
2617       edm::LogInfo("PrimaryVertexValidation")
2618           << " pixel phase0 setup, nLadders: " << nLadders_ << " nModules:" << nModZ_;
2619     }
2620   }
2621
2622   switch (phase_) {
2623     case PVValHelper::phase0:
2624       etaOfProbe_ = std::min(etaOfProbe_, PVValHelper::max_eta_phase0);
2625       break;
2626     case PVValHelper::phase1:
2627       etaOfProbe_ = std::min(etaOfProbe_, PVValHelper::max_eta_phase1);
2628       break;
2629     case PVValHelper::phase2:
2630       etaOfProbe_ = std::min(etaOfProbe_, PVValHelper::max_eta_phase2);
2631       break;
2632     default:
2633       throw cms::Exception("LogicError") << "Unknown detector phase: " << phase_;
2634   }
2635
2636   if (h_etaMax->GetEntries() == 0.) {
2637     h_etaMax->SetBinContent(1., etaOfProbe_);
2638     h_nbins->SetBinContent(1., nBins_);
2639     h_nLadders->SetBinContent(1., nLadders_);
2640     h_nModZ->SetBinContent(1., nModZ_);
2641     h_pTinfo->SetBinContent(1., mypT_bins_.size());
2642     h_pTinfo->SetBinContent(2., minPt_);
2643     h_pTinfo->SetBinContent(3., maxPt_);
2644   }
2645 }
2646
2647 // ------------ method called once each job just after ending the event loop  ------------
2648 void PrimaryVertexValidation::endJob() {
2649   // shring the histograms to fit
2650   h_probeL1Ladder_->GetXaxis()->SetRangeUser(-1.5, nLadders_ + 0.5);
2651   h_probeL1Module_->GetXaxis()->SetRangeUser(-1.5, nModZ_ + 0.5);
2652   h2_probeLayer1Map_->GetXaxis()->SetRangeUser(0.5, nModZ_ + 0.5);
2653   h2_probeLayer1Map_->GetYaxis()->SetRangeUser(0.5, nLadders_ + 0.5);
2654   h2_probePassingLayer1Map_->GetXaxis()->SetRangeUser(0.5, nModZ_ + 0.5);
2655   h2_probePassingLayer1Map_->GetYaxis()->SetRangeUser(0.5, nLadders_ + 0.5);
2656
2657   TFileDirectory RunFeatures = fs->mkdir("RunFeatures");
2658   h_runStartTimes = RunFeatures.make<TH1I>(
2659       "runStartTimes", "run start times", runNumbersTimesLog_.size(), 0, runNumbersTimesLog_.size());
2660   h_runEndTimes =
2661       RunFeatures.make<TH1I>("runEndTimes", "run end times", runNumbersTimesLog_.size(), 0, runNumbersTimesLog_.size());
2662
2663   unsigned int count = 1;
2664   for (const auto& run : runNumbersTimesLog_) {
2665     // strip down the microseconds
2666     h_runStartTimes->SetBinContent(count, run.second.first / 10e6);
2667     h_runStartTimes->GetXaxis()->SetBinLabel(count, (std::to_string(run.first)).c_str());
2668
2669     h_runEndTimes->SetBinContent(count, run.second.second / 10e6);
2670     h_runEndTimes->GetXaxis()->SetBinLabel(count, (std::to_string(run.first)).c_str());
2671
2672     count++;
2673   }
2674
2675   edm::LogInfo("PrimaryVertexValidation") << "######################################\n"
2676                                           << "# PrimaryVertexValidation::endJob()\n"
2677                                           << "# Number of analyzed events: " << Nevt_ << "\n"
2678                                           << "######################################";
2679
2680   // means and widhts vs ladder and module number
2681
2682   a_dxymodZMeanTrend = MeanTrendsDir.make<TH1F>(
2683       "means_dxy_modZ", "#LT d_{xy} #GT vs modZ;module number (Z);#LT d_{xy} #GT [#mum]", nModZ_, 0., nModZ_);
2684
2685   a_dxymodZWidthTrend = WidthTrendsDir.make<TH1F>(
2686       "widths_dxy_modZ", "#sigma_{d_{xy}} vs modZ;module number (Z);#sigma_{d_{xy}} [#mum]", nModZ_, 0., nModZ_);
2687
2688   a_dzmodZMeanTrend = MeanTrendsDir.make<TH1F>(
2689       "means_dz_modZ", "#LT d_{z} #GT vs modZ;module number (Z);#LT d_{z} #GT [#mum]", nModZ_, 0., nModZ_);
2690
2691   a_dzmodZWidthTrend = WidthTrendsDir.make<TH1F>(
2692       "widths_dz_modZ", "#sigma_{d_{z}} vs modZ;module number (Z);#sigma_{d_{z}} [#mum]", nModZ_, 0., nModZ_);
2693
2694   a_dxyladderMeanTrend = MeanTrendsDir.make<TH1F>("means_dxy_ladder",
2695                                                   "#LT d_{xy} #GT vs ladder;ladder number (#phi);#LT d_{xy} #GT [#mum]",
2696                                                   nLadders_,
2697                                                   0.,
2698                                                   nLadders_);
2699
2700   a_dxyladderWidthTrend =
2701       WidthTrendsDir.make<TH1F>("widths_dxy_ladder",
2702                                 "#sigma_{d_{xy}} vs ladder;ladder number (#phi);#sigma_{d_{xy}} [#mum]",
2703                                 nLadders_,
2704                                 0.,
2705                                 nLadders_);
2706
2707   a_dzladderMeanTrend = MeanTrendsDir.make<TH1F>(
2708       "means_dz_ladder", "#LT d_{z} #GT vs ladder;ladder number (#phi);#LT d_{z} #GT [#mum]", nLadders_, 0., nLadders_);
2709
2710   a_dzladderWidthTrend =
2711       WidthTrendsDir.make<TH1F>("widths_dz_ladder",
2712                                 "#sigma_{d_{z}} vs ladder;ladder number (#phi);#sigma_{d_{z}} [#mum]",
2713                                 nLadders_,
2714                                 0.,
2715                                 nLadders_);
2716
2717   n_dxymodZMeanTrend = MeanTrendsDir.make<TH1F>(
2718       "norm_means_dxy_modZ",
2719       "#LT d_{xy}/#sigma_{d_{xy}} #GT vs modZ;module number (Z);#LT d_{xy}/#sigma_{d_{xy}} #GT",
2720       nModZ_,
2721       0.,
2722       nModZ_);
2723
2724   n_dxymodZWidthTrend = WidthTrendsDir.make<TH1F>(
2725       "norm_widths_dxy_modZ",
2726       "width(d_{xy}/#sigma_{d_{xy}}) vs modZ;module number (Z); width(d_{xy}/#sigma_{d_{xy}})",
2727       nModZ_,
2728       0.,
2729       nModZ_);
2730
2731   n_dzmodZMeanTrend =
2732       MeanTrendsDir.make<TH1F>("norm_means_dz_modZ",
2733                                "#LT d_{z}/#sigma_{d_{z}} #GT vs modZ;module number (Z);#LT d_{z}/#sigma_{d_{z}} #GT",
2734                                nModZ_,
2735                                0.,
2736                                nModZ_);
2737
2738   n_dzmodZWidthTrend =
2739       WidthTrendsDir.make<TH1F>("norm_widths_dz_modZ",
2740                                 "width(d_{z}/#sigma_{d_{z}}) vs pT;module number (Z);width(d_{z}/#sigma_{d_{z}})",
2741                                 nModZ_,
2742                                 0.,
2743                                 nModZ_);
2744
2745   n_dxyladderMeanTrend = MeanTrendsDir.make<TH1F>(
2746       "norm_means_dxy_ladder",
2747       "#LT d_{xy}/#sigma_{d_{xy}} #GT vs ladder;ladder number (#phi);#LT d_{xy}/#sigma_{d_{z}} #GT",
2748       nLadders_,
2749       0.,
2750       nLadders_);
2751
2752   n_dxyladderWidthTrend = WidthTrendsDir.make<TH1F>(
2753       "norm_widths_dxy_ladder",
2754       "width(d_{xy}/#sigma_{d_{xy}}) vs ladder;ladder number (#phi);width(d_{xy}/#sigma_{d_{z}})",
2755       nLadders_,
2756       0.,
2757       nLadders_);
2758
2759   n_dzladderMeanTrend = MeanTrendsDir.make<TH1F>(
2760       "norm_means_dz_ladder",
2761       "#LT d_{z}/#sigma_{d_{z}} #GT vs ladder;ladder number (#phi);#LT d_{z}/#sigma_{d_{z}} #GT",
2762       nLadders_,
2763       0.,
2764       nLadders_);
2765
2766   n_dzladderWidthTrend = WidthTrendsDir.make<TH1F>(
2767       "norm_widths_dz_ladder",
2768       "width(d_{z}/#sigma_{d_{z}}) vs ladder;ladder number (#phi);width(d_{z}/#sigma_{d_{z}})",
2769       nLadders_,
2770       0.,
2771       nLadders_);
2772
2773   // 2D maps of residuals in bins of L1 modules
2774
2775   a_dxyL1MeanMap = Mean2DMapsDir.make<TH2F>("means_dxy_l1map",
2776                                             "#LT d_{xy} #GT map;module number [z];ladder number [#varphi]",
2777                                             nModZ_,
2778                                             0.,
2779                                             nModZ_,
2780                                             nLadders_,
2781                                             0.,
2782                                             nLadders_);
2783
2784   a_dzL1MeanMap = Mean2DMapsDir.make<TH2F>("means_dz_l1map",
2785                                            "#LT d_{z} #GT map;module number [z];ladder number [#varphi]",
2786                                            nModZ_,
2787                                            0.,
2788                                            nModZ_,
2789                                            nLadders_,
2790                                            0.,
2791                                            nLadders_);
2792
2793   n_dxyL1MeanMap =
2794       Mean2DMapsDir.make<TH2F>("norm_means_dxy_l1map",
2795                                "#LT d_{xy}/#sigma_{d_{xy}} #GT map;module number [z];ladder number [#varphi]",
2796                                nModZ_,
2797                                0.,
2798                                nModZ_,
2799                                nLadders_,
2800                                0.,
2801                                nLadders_);
2802
2803   n_dzL1MeanMap = Mean2DMapsDir.make<TH2F>("norm_means_dz_l1map",
2804                                            "#LT d_{z}/#sigma_{d_{z}} #GT map;module number [z];ladder number [#varphi]",
2805                                            nModZ_,
2806                                            0.,
2807                                            nModZ_,
2808                                            nLadders_,
2809                                            0.,
2810                                            nLadders_);
2811
2812   a_dxyL1WidthMap = Width2DMapsDir.make<TH2F>("widths_dxy_l1map",
2813                                               "#sigma_{d_{xy}} map;module number [z];ladder number [#varphi]",
2814                                               nModZ_,
2815                                               0.,
2816                                               nModZ_,
2817                                               nLadders_,
2818                                               0.,
2819                                               nLadders_);
2820
2821   a_dzL1WidthMap = Width2DMapsDir.make<TH2F>("widths_dz_l1map",
2822                                              "#sigma_{d_{z}} map;module number [z];ladder number [#varphi]",
2823                                              nModZ_,
2824                                              0.,
2825                                              nModZ_,
2826                                              nLadders_,
2827                                              0.,
2828                                              nLadders_);
2829
2830   n_dxyL1WidthMap =
2831       Width2DMapsDir.make<TH2F>("norm_widths_dxy_l1map",
2832                                 "width(d_{xy}/#sigma_{d_{xy}}) map;module number [z];ladder number [#varphi]",
2833                                 nModZ_,
2834                                 0.,
2835                                 nModZ_,
2836                                 nLadders_,
2837                                 0.,
2838                                 nLadders_);
2839
2840   n_dzL1WidthMap =
2841       Width2DMapsDir.make<TH2F>("norm_widths_dz_l1map",
2842                                 "width(d_{z}/#sigma_{d_{z}}) map;module number [z];ladder number [#varphi]",
2843                                 nModZ_,
2844                                 0.,
2845                                 nModZ_,
2846                                 nLadders_,
2847                                 0.,
2848                                 nLadders_);
2849
2850   if (useTracksFromRecoVtx_) {
2851     fillTrendPlotByIndex(a_dxyPhiMeanBiasTrend, a_dxyPhiBiasResiduals, PVValHelper::MEAN, PVValHelper::phi);
2852     fillTrendPlotByIndex(a_dxyPhiWidthBiasTrend, a_dxyPhiBiasResiduals, PVValHelper::WIDTH, PVValHelper::phi);
2853     fillTrendPlotByIndex(a_dzPhiMeanBiasTrend, a_dzPhiBiasResiduals, PVValHelper::MEAN, PVValHelper::phi);
2854     fillTrendPlotByIndex(a_dzPhiWidthBiasTrend, a_dzPhiBiasResiduals, PVValHelper::WIDTH, PVValHelper::phi);
2855
2856     fillTrendPlotByIndex(a_dxyEtaMeanBiasTrend, a_dxyEtaBiasResiduals, PVValHelper::MEAN, PVValHelper::eta);
2857     fillTrendPlotByIndex(a_dxyEtaWidthBiasTrend, a_dxyEtaBiasResiduals, PVValHelper::WIDTH, PVValHelper::eta);
2858     fillTrendPlotByIndex(a_dzEtaMeanBiasTrend, a_dzEtaBiasResiduals, PVValHelper::MEAN, PVValHelper::eta);
2859     fillTrendPlotByIndex(a_dzEtaWidthBiasTrend, a_dzEtaBiasResiduals, PVValHelper::WIDTH, PVValHelper::eta);
2860
2861     fillTrendPlotByIndex(n_dxyPhiMeanBiasTrend, n_dxyPhiBiasResiduals, PVValHelper::MEAN, PVValHelper::phi);
2862     fillTrendPlotByIndex(n_dxyPhiWidthBiasTrend, n_dxyPhiBiasResiduals, PVValHelper::WIDTH, PVValHelper::phi);
2863     fillTrendPlotByIndex(n_dzPhiMeanBiasTrend, n_dzPhiBiasResiduals, PVValHelper::MEAN, PVValHelper::phi);
2864     fillTrendPlotByIndex(n_dzPhiWidthBiasTrend, n_dzPhiBiasResiduals, PVValHelper::WIDTH, PVValHelper::phi);
2865
2866     fillTrendPlotByIndex(n_dxyEtaMeanBiasTrend, n_dxyEtaBiasResiduals, PVValHelper::MEAN, PVValHelper::eta);
2867     fillTrendPlotByIndex(n_dxyEtaWidthBiasTrend, n_dxyEtaBiasResiduals, PVValHelper::WIDTH, PVValHelper::eta);
2868     fillTrendPlotByIndex(n_dzEtaMeanBiasTrend, n_dzEtaBiasResiduals, PVValHelper::MEAN, PVValHelper::eta);
2869     fillTrendPlotByIndex(n_dzEtaWidthBiasTrend, n_dzEtaBiasResiduals, PVValHelper::WIDTH, PVValHelper::eta);
2870
2871     // medians and MADs
2872
2873     fillTrendPlotByIndex(a_dxyPhiMedianBiasTrend, a_dxyPhiBiasResiduals, PVValHelper::MEDIAN, PVValHelper::phi);
2874     fillTrendPlotByIndex(a_dxyPhiMADBiasTrend, a_dxyPhiBiasResiduals, PVValHelper::MAD, PVValHelper::phi);
2875     fillTrendPlotByIndex(a_dzPhiMedianBiasTrend, a_dzPhiBiasResiduals, PVValHelper::MEDIAN, PVValHelper::phi);
2876     fillTrendPlotByIndex(a_dzPhiMADBiasTrend, a_dzPhiBiasResiduals, PVValHelper::MAD, PVValHelper::phi);
2877
2878     fillTrendPlotByIndex(a_dxyEtaMedianBiasTrend, a_dxyEtaBiasResiduals, PVValHelper::MEDIAN, PVValHelper::eta);
2879     fillTrendPlotByIndex(a_dxyEtaMADBiasTrend, a_dxyEtaBiasResiduals, PVValHelper::MAD, PVValHelper::eta);
2880     fillTrendPlotByIndex(a_dzEtaMedianBiasTrend, a_dzEtaBiasResiduals, PVValHelper::MEDIAN, PVValHelper::eta);
2881     fillTrendPlotByIndex(a_dzEtaMADBiasTrend, a_dzEtaBiasResiduals, PVValHelper::MAD, PVValHelper::eta);
2882
2883     fillTrendPlotByIndex(n_dxyPhiMedianBiasTrend, n_dxyPhiBiasResiduals, PVValHelper::MEDIAN, PVValHelper::phi);
2884     fillTrendPlotByIndex(n_dxyPhiMADBiasTrend, n_dxyPhiBiasResiduals, PVValHelper::MAD, PVValHelper::phi);
2885     fillTrendPlotByIndex(n_dzPhiMedianBiasTrend, n_dzPhiBiasResiduals, PVValHelper::MEDIAN, PVValHelper::phi);
2886     fillTrendPlotByIndex(n_dzPhiMADBiasTrend, n_dzPhiBiasResiduals, PVValHelper::MAD, PVValHelper::phi);
2887
2888     fillTrendPlotByIndex(n_dxyEtaMedianBiasTrend, n_dxyEtaBiasResiduals, PVValHelper::MEDIAN, PVValHelper::eta);
2889     fillTrendPlotByIndex(n_dxyEtaMADBiasTrend, n_dxyEtaBiasResiduals, PVValHelper::MAD, PVValHelper::eta);
2890     fillTrendPlotByIndex(n_dzEtaMedianBiasTrend, n_dzEtaBiasResiduals, PVValHelper::MEDIAN, PVValHelper::eta);
2891     fillTrendPlotByIndex(n_dzEtaMADBiasTrend, n_dzEtaBiasResiduals, PVValHelper::MAD, PVValHelper::eta);
2892
2893     // 2d Maps
2894
2895     fillMap(a_dxyMeanBiasMap, a_dxyBiasResidualsMap, PVValHelper::MEAN, nBins_, nBins_);
2896     fillMap(a_dxyWidthBiasMap, a_dxyBiasResidualsMap, PVValHelper::WIDTH, nBins_, nBins_);
2897     fillMap(a_dzMeanBiasMap, a_dzBiasResidualsMap, PVValHelper::MEAN, nBins_, nBins_);
2898     fillMap(a_dzWidthBiasMap, a_dzBiasResidualsMap, PVValHelper::WIDTH, nBins_, nBins_);
2899
2900     fillMap(n_dxyMeanBiasMap, n_dxyBiasResidualsMap, PVValHelper::MEAN, nBins_, nBins_);
2901     fillMap(n_dxyWidthBiasMap, n_dxyBiasResidualsMap, PVValHelper::WIDTH, nBins_, nBins_);
2902     fillMap(n_dzMeanBiasMap, n_dzBiasResidualsMap, PVValHelper::MEAN, nBins_, nBins_);
2903     fillMap(n_dzWidthBiasMap, n_dzBiasResidualsMap, PVValHelper::WIDTH, nBins_, nBins_);
2904   }
2905
2906   // do profiles
2907
2908   fillTrendPlotByIndex(a_dxyPhiMeanTrend, a_dxyPhiResiduals, PVValHelper::MEAN, PVValHelper::phi);
2909   fillTrendPlotByIndex(a_dxyPhiWidthTrend, a_dxyPhiResiduals, PVValHelper::WIDTH, PVValHelper::phi);
2910   fillTrendPlotByIndex(a_dzPhiMeanTrend, a_dzPhiResiduals, PVValHelper::MEAN, PVValHelper::phi);
2911   fillTrendPlotByIndex(a_dzPhiWidthTrend, a_dzPhiResiduals, PVValHelper::WIDTH, PVValHelper::phi);
2912
2913   fillTrendPlotByIndex(a_dxyEtaMeanTrend, a_dxyEtaResiduals, PVValHelper::MEAN, PVValHelper::eta);
2914   fillTrendPlotByIndex(a_dxyEtaWidthTrend, a_dxyEtaResiduals, PVValHelper::WIDTH, PVValHelper::eta);
2915   fillTrendPlotByIndex(a_dzEtaMeanTrend, a_dzEtaResiduals, PVValHelper::MEAN, PVValHelper::eta);
2916   fillTrendPlotByIndex(a_dzEtaWidthTrend, a_dzEtaResiduals, PVValHelper::WIDTH, PVValHelper::eta);
2917
2918   fillTrendPlotByIndex(n_dxyPhiMeanTrend, n_dxyPhiResiduals, PVValHelper::MEAN, PVValHelper::phi);
2919   fillTrendPlotByIndex(n_dxyPhiWidthTrend, n_dxyPhiResiduals, PVValHelper::WIDTH, PVValHelper::phi);
2920   fillTrendPlotByIndex(n_dzPhiMeanTrend, n_dzPhiResiduals, PVValHelper::MEAN, PVValHelper::phi);
2921   fillTrendPlotByIndex(n_dzPhiWidthTrend, n_dzPhiResiduals, PVValHelper::WIDTH);
2922
2923   fillTrendPlotByIndex(n_dxyEtaMeanTrend, n_dxyEtaResiduals, PVValHelper::MEAN, PVValHelper::eta);
2924   fillTrendPlotByIndex(n_dxyEtaWidthTrend, n_dxyEtaResiduals, PVValHelper::WIDTH, PVValHelper::eta);
2925   fillTrendPlotByIndex(n_dzEtaMeanTrend, n_dzEtaResiduals, PVValHelper::MEAN, PVValHelper::eta);
2926   fillTrendPlotByIndex(n_dzEtaWidthTrend, n_dzEtaResiduals, PVValHelper::WIDTH, PVValHelper::eta);
2927
2928   // vs transverse momentum
2929
2930   fillTrendPlotByIndex(a_dxypTMeanTrend, h_dxy_pT_, PVValHelper::MEAN);
2931   fillTrendPlotByIndex(a_dxypTWidthTrend, h_dxy_pT_, PVValHelper::WIDTH);
2932   fillTrendPlotByIndex(a_dzpTMeanTrend, h_dz_pT_, PVValHelper::MEAN);
2933   fillTrendPlotByIndex(a_dzpTWidthTrend, h_dz_pT_, PVValHelper::WIDTH);
2934
2935   fillTrendPlotByIndex(a_dxypTCentralMeanTrend, h_dxy_Central_pT_, PVValHelper::MEAN);
2936   fillTrendPlotByIndex(a_dxypTCentralWidthTrend, h_dxy_Central_pT_, PVValHelper::WIDTH);
2937   fillTrendPlotByIndex(a_dzpTCentralMeanTrend, h_dz_Central_pT_, PVValHelper::MEAN);
2938   fillTrendPlotByIndex(a_dzpTCentralWidthTrend, h_dz_Central_pT_, PVValHelper::WIDTH);
2939
2940   fillTrendPlotByIndex(n_dxypTMeanTrend, h_norm_dxy_pT_, PVValHelper::MEAN);
2941   fillTrendPlotByIndex(n_dxypTWidthTrend, h_norm_dxy_pT_, PVValHelper::WIDTH);
2942   fillTrendPlotByIndex(n_dzpTMeanTrend, h_norm_dz_pT_, PVValHelper::MEAN);
2943   fillTrendPlotByIndex(n_dzpTWidthTrend, h_norm_dz_pT_, PVValHelper::WIDTH);
2944
2945   fillTrendPlotByIndex(n_dxypTCentralMeanTrend, h_norm_dxy_Central_pT_, PVValHelper::MEAN);
2946   fillTrendPlotByIndex(n_dxypTCentralWidthTrend, h_norm_dxy_Central_pT_, PVValHelper::WIDTH);
2947   fillTrendPlotByIndex(n_dzpTCentralMeanTrend, h_norm_dz_Central_pT_, PVValHelper::MEAN);
2948   fillTrendPlotByIndex(n_dzpTCentralWidthTrend, h_norm_dz_Central_pT_, PVValHelper::WIDTH);
2949
2950   // vs ladder and module number
2951
2952   fillTrendPlotByIndex(a_dxymodZMeanTrend, h_dxy_modZ_, PVValHelper::MEAN);
2953   fillTrendPlotByIndex(a_dxymodZWidthTrend, h_dxy_modZ_, PVValHelper::WIDTH);
2954   fillTrendPlotByIndex(a_dzmodZMeanTrend, h_dz_modZ_, PVValHelper::MEAN);
2955   fillTrendPlotByIndex(a_dzmodZWidthTrend, h_dz_modZ_, PVValHelper::WIDTH);
2956
2957   fillTrendPlotByIndex(a_dxyladderMeanTrend, h_dxy_ladder_, PVValHelper::MEAN);
2958   fillTrendPlotByIndex(a_dxyladderWidthTrend, h_dxy_ladder_, PVValHelper::WIDTH);
2959   fillTrendPlotByIndex(a_dzladderMeanTrend, h_dz_ladder_, PVValHelper::MEAN);
2960   fillTrendPlotByIndex(a_dzladderWidthTrend, h_dz_ladder_, PVValHelper::WIDTH);
2961
2962   fillTrendPlotByIndex(n_dxymodZMeanTrend, h_norm_dxy_modZ_, PVValHelper::MEAN);
2963   fillTrendPlotByIndex(n_dxymodZWidthTrend, h_norm_dxy_modZ_, PVValHelper::WIDTH);
2964   fillTrendPlotByIndex(n_dzmodZMeanTrend, h_norm_dz_modZ_, PVValHelper::MEAN);
2965   fillTrendPlotByIndex(n_dzmodZWidthTrend, h_norm_dz_modZ_, PVValHelper::WIDTH);
2966
2967   fillTrendPlotByIndex(n_dxyladderMeanTrend, h_norm_dxy_ladder_, PVValHelper::MEAN);
2968   fillTrendPlotByIndex(n_dxyladderWidthTrend, h_norm_dxy_ladder_, PVValHelper::WIDTH);
2969   fillTrendPlotByIndex(n_dzladderMeanTrend, h_norm_dz_ladder_, PVValHelper::MEAN);
2970   fillTrendPlotByIndex(n_dzladderWidthTrend, h_norm_dz_ladder_, PVValHelper::WIDTH);
2971
2972   // medians and MADs
2973
2974   fillTrendPlotByIndex(a_dxyPhiMedianTrend, a_dxyPhiResiduals, PVValHelper::MEDIAN, PVValHelper::phi);
2975   fillTrendPlotByIndex(a_dxyPhiMADTrend, a_dxyPhiResiduals, PVValHelper::MAD, PVValHelper::phi);
2976
2977   fillTrendPlotByIndex(a_dzPhiMedianTrend, a_dzPhiResiduals, PVValHelper::MEDIAN, PVValHelper::phi);
2978   fillTrendPlotByIndex(a_dzPhiMADTrend, a_dzPhiResiduals, PVValHelper::MAD, PVValHelper::phi);
2979
2980   fillTrendPlotByIndex(a_dxyEtaMedianTrend, a_dxyEtaResiduals, PVValHelper::MEDIAN, PVValHelper::eta);
2981   fillTrendPlotByIndex(a_dxyEtaMADTrend, a_dxyEtaResiduals, PVValHelper::MAD, PVValHelper::eta);
2982   fillTrendPlotByIndex(a_dzEtaMedianTrend, a_dzEtaResiduals, PVValHelper::MEDIAN, PVValHelper::eta);
2983   fillTrendPlotByIndex(a_dzEtaMADTrend, a_dzEtaResiduals, PVValHelper::MAD, PVValHelper::eta);
2984
2985   fillTrendPlotByIndex(n_dxyPhiMedianTrend, n_dxyPhiResiduals, PVValHelper::MEDIAN, PVValHelper::phi);
2986   fillTrendPlotByIndex(n_dxyPhiMADTrend, n_dxyPhiResiduals, PVValHelper::MAD, PVValHelper::phi);
2987   fillTrendPlotByIndex(n_dzPhiMedianTrend, n_dzPhiResiduals, PVValHelper::MEDIAN, PVValHelper::phi);
2988   fillTrendPlotByIndex(n_dzPhiMADTrend, n_dzPhiResiduals, PVValHelper::MAD, PVValHelper::phi);
2989
2990   fillTrendPlotByIndex(n_dxyEtaMedianTrend, n_dxyEtaResiduals, PVValHelper::MEDIAN, PVValHelper::eta);
2991   fillTrendPlotByIndex(n_dxyEtaMADTrend, n_dxyEtaResiduals, PVValHelper::MAD, PVValHelper::eta);
2992   fillTrendPlotByIndex(n_dzEtaMedianTrend, n_dzEtaResiduals, PVValHelper::MEDIAN, PVValHelper::eta);
2993   fillTrendPlotByIndex(n_dzEtaMADTrend, n_dzEtaResiduals, PVValHelper::MAD, PVValHelper::eta);
2994
2995   // 2D Maps
2996
2997   fillMap(a_dxyMeanMap, a_dxyResidualsMap, PVValHelper::MEAN, nBins_, nBins_);
2998   fillMap(a_dxyWidthMap, a_dxyResidualsMap, PVValHelper::WIDTH, nBins_, nBins_);
2999   fillMap(a_dzMeanMap, a_dzResidualsMap, PVValHelper::MEAN, nBins_, nBins_);
3000   fillMap(a_dzWidthMap, a_dzResidualsMap, PVValHelper::WIDTH, nBins_, nBins_);
3001
3002   fillMap(n_dxyMeanMap, n_dxyResidualsMap, PVValHelper::MEAN, nBins_, nBins_);
3003   fillMap(n_dxyWidthMap, n_dxyResidualsMap, PVValHelper::WIDTH, nBins_, nBins_);
3004   fillMap(n_dzMeanMap, n_dzResidualsMap, PVValHelper::MEAN, nBins_, nBins_);
3005   fillMap(n_dzWidthMap, n_dzResidualsMap, PVValHelper::WIDTH, nBins_, nBins_);
3006
3007   // 2D Maps of residuals in bins of L1 modules
3008
3009   fillMap(a_dxyL1MeanMap, a_dxyL1ResidualsMap, PVValHelper::MEAN, nModZ_, nLadders_);
3010   fillMap(a_dxyL1WidthMap, a_dxyL1ResidualsMap, PVValHelper::WIDTH, nModZ_, nLadders_);
3011   fillMap(a_dzL1MeanMap, a_dzL1ResidualsMap, PVValHelper::MEAN, nModZ_, nLadders_);
3012   fillMap(a_dzL1WidthMap, a_dzL1ResidualsMap, PVValHelper::WIDTH, nModZ_, nLadders_);
3013
3014   fillMap(n_dxyL1MeanMap, n_dxyL1ResidualsMap, PVValHelper::MEAN, nModZ_, nLadders_);
3015   fillMap(n_dxyL1WidthMap, n_dxyL1ResidualsMap, PVValHelper::WIDTH, nModZ_, nLadders_);
3016   fillMap(n_dzL1MeanMap, n_dzL1ResidualsMap, PVValHelper::MEAN, nModZ_, nLadders_);
3017   fillMap(n_dzL1WidthMap, n_dzL1ResidualsMap, PVValHelper::WIDTH, nModZ_, nLadders_);
3018 }
3019
3020 //*************************************************************
3021 std::pair<long long, long long> PrimaryVertexValidation::getRunTime(const edm::EventSetup& iSetup) const
3022 //*************************************************************
3023 {
3024   edm::ESHandle<RunInfo> runInfo = iSetup.getHandle(runInfoTokenBR_);
3025   if (debug_) {
3026     edm::LogInfo("PrimaryVertexValidation")
3027         << runInfo.product()->m_start_time_str << " " << runInfo.product()->m_stop_time_str << std::endl;
3028   }
3029   return std::make_pair(runInfo.product()->m_start_time_ll, runInfo.product()->m_stop_time_ll);
3030 }
3031
3032 //*************************************************************
3033 bool PrimaryVertexValidation::isBFieldConsistentWithMode(const edm::EventSetup& iSetup) const
3034 //*************************************************************
3035 {
3036   edm::ESHandle<RunInfo> runInfo = iSetup.getHandle(runInfoTokenBR_);
3037   double average_current = runInfo.product()->m_avg_current;
3038   bool isOn = (average_current > 2000.);
3039   bool is0T = (ptOfProbe_ == 0.);
3040
3041   return ((isOn && !is0T) || (!isOn && is0T));
3042 }
3043
3044 //*************************************************************
3045 void PrimaryVertexValidation::SetVarToZero()
3046 //*************************************************************
3047 {
3048   nTracks_ = 0;
3049   nClus_ = 0;
3050   nOfflineVertices_ = 0;
3051   RunNumber_ = 0;
3052   LuminosityBlockNumber_ = 0;
3053   xOfflineVertex_ = -999.;
3054   yOfflineVertex_ = -999.;
3055   zOfflineVertex_ = -999.;
3056   xErrOfflineVertex_ = 0.;
3057   yErrOfflineVertex_ = 0.;
3058   zErrOfflineVertex_ = 0.;
3059   BSx0_ = -999.;
3060   BSy0_ = -999.;
3061   BSz0_ = -999.;
3062   Beamsigmaz_ = -999.;
3063   Beamdxdz_ = -999.;
3064   BeamWidthX_ = -999.;
3065   BeamWidthY_ = -999.;
3066   wxy2_ = -999.;
3067
3068   for (int i = 0; i < nMaxtracks_; ++i) {
3069     pt_[i] = 0;
3070     p_[i] = 0;
3071     nhits_[i] = 0;
3072     nhits1D_[i] = 0;
3073     nhits2D_[i] = 0;
3074     nhitsBPIX_[i] = 0;
3075     nhitsFPIX_[i] = 0;
3076     nhitsTIB_[i] = 0;
3077     nhitsTID_[i] = 0;
3078     nhitsTOB_[i] = 0;
3079     nhitsTEC_[i] = 0;
3080     isHighPurity_[i] = 0;
3081     eta_[i] = 0;
3082     theta_[i] = 0;
3083     phi_[i] = 0;
3084     chi2_[i] = 0;
3085     chi2ndof_[i] = 0;
3086     charge_[i] = 0;
3087     qoverp_[i] = 0;
3088     dz_[i] = 0;
3089     dxy_[i] = 0;
3090     dzBs_[i] = 0;
3091     dxyBs_[i] = 0;
3092     xPCA_[i] = 0;
3093     yPCA_[i] = 0;
3094     zPCA_[i] = 0;
3095     xUnbiasedVertex_[i] = 0;
3096     yUnbiasedVertex_[i] = 0;
3097     zUnbiasedVertex_[i] = 0;
3098     chi2normUnbiasedVertex_[i] = 0;
3099     chi2UnbiasedVertex_[i] = 0;
3100     chi2ProbUnbiasedVertex_[i] = 0;
3101     DOFUnbiasedVertex_[i] = 0;
3102     sumOfWeightsUnbiasedVertex_[i] = 0;
3103     tracksUsedForVertexing_[i] = 0;
3104     dxyFromMyVertex_[i] = 0;
3105     dzFromMyVertex_[i] = 0;
3106     d3DFromMyVertex_[i] = 0;
3107     dxyErrorFromMyVertex_[i] = 0;
3108     dzErrorFromMyVertex_[i] = 0;
3109     d3DErrorFromMyVertex_[i] = 0;
3110     IPTsigFromMyVertex_[i] = 0;
3111     IPLsigFromMyVertex_[i] = 0;
3112     IP3DsigFromMyVertex_[i] = 0;
3113     hasRecVertex_[i] = 0;
3114     isGoodTrack_[i] = 0;
3115   }
3116 }
3117
3118 //*************************************************************
3119 void PrimaryVertexValidation::fillTrendPlot(TH1F* trendPlot,
3120                                             TH1F* residualsPlot[100],
3121                                             PVValHelper::estimator fitPar_,
3122                                             const std::string& var_)
3123 //*************************************************************
3124 {
3125   for (int i = 0; i < nBins_; ++i) {
3126     char phibincenter[129];
3127     auto phiBins = theDetails_.trendbins[PVValHelper::phi];
3128     sprintf(phibincenter, "%.f", (phiBins[i] + phiBins[i + 1]) / 2.);
3129
3130     char etabincenter[129];
3131     auto etaBins = theDetails_.trendbins[PVValHelper::eta];
3132     sprintf(etabincenter, "%.1f", (etaBins[i] + etaBins[i + 1]) / 2.);
3133
3134     switch (fitPar_) {
3135       case PVValHelper::MEAN: {
3136         float mean_ = PVValHelper::fitResiduals(residualsPlot[i]).first.value();
3137         float meanErr_ = PVValHelper::fitResiduals(residualsPlot[i]).first.error();
3138         trendPlot->SetBinContent(i + 1, mean_);
3139         trendPlot->SetBinError(i + 1, meanErr_);
3140         break;
3141       }
3142       case PVValHelper::WIDTH: {
3143         float width_ = PVValHelper::fitResiduals(residualsPlot[i]).second.value();
3144         float widthErr_ = PVValHelper::fitResiduals(residualsPlot[i]).second.error();
3145         trendPlot->SetBinContent(i + 1, width_);
3146         trendPlot->SetBinError(i + 1, widthErr_);
3147         break;
3148       }
3149       case PVValHelper::MEDIAN: {
3150         float median_ = PVValHelper::getMedian(residualsPlot[i]).value();
3151         float medianErr_ = PVValHelper::getMedian(residualsPlot[i]).error();
3152         trendPlot->SetBinContent(i + 1, median_);
3153         trendPlot->SetBinError(i + 1, medianErr_);
3154         break;
3155       }
3156       case PVValHelper::MAD: {
3157         float mad_ = PVValHelper::getMAD(residualsPlot[i]).value();
3158         float madErr_ = PVValHelper::getMAD(residualsPlot[i]).error();
3159         trendPlot->SetBinContent(i + 1, mad_);
3160         trendPlot->SetBinError(i + 1, madErr_);
3161         break;
3162       }
3163       default:
3164         edm::LogWarning("PrimaryVertexValidation")
3165             << "fillTrendPlot() " << fitPar_ << " unknown estimator!" << std::endl;
3166         break;
3167     }
3168
3169     if (var_.find("eta") != std::string::npos) {
3170       trendPlot->GetXaxis()->SetBinLabel(i + 1, etabincenter);
3171     } else if (var_.find("phi") != std::string::npos) {
3172       trendPlot->GetXaxis()->SetBinLabel(i + 1, phibincenter);
3173     } else {
3174       edm::LogWarning("PrimaryVertexValidation")
3175           << "fillTrendPlot() " << var_ << " unknown track parameter!" << std::endl;
3176     }
3177   }
3178 }
3179
3180 //*************************************************************
3181 void PrimaryVertexValidation::fillTrendPlotByIndex(TH1F* trendPlot,
3182                                                    std::vector<TH1F*>& h,
3183                                                    PVValHelper::estimator fitPar_,
3184                                                    PVValHelper::plotVariable plotVar)
3185 //*************************************************************
3186 {
3187   for (auto iterator = h.begin(); iterator != h.end(); iterator++) {
3188     unsigned int bin = std::distance(h.begin(), iterator) + 1;
3189     std::pair<Measurement1D, Measurement1D> myFit = PVValHelper::fitResiduals((*iterator));
3190
3191     switch (fitPar_) {
3192       case PVValHelper::MEAN: {
3193         float mean_ = myFit.first.value();
3194         float meanErr_ = myFit.first.error();
3195         trendPlot->SetBinContent(bin, mean_);
3196         trendPlot->SetBinError(bin, meanErr_);
3197         break;
3198       }
3199       case PVValHelper::WIDTH: {
3200         float width_ = myFit.second.value();
3201         float widthErr_ = myFit.second.error();
3202         trendPlot->SetBinContent(bin, width_);
3203         trendPlot->SetBinError(bin, widthErr_);
3204         break;
3205       }
3206       case PVValHelper::MEDIAN: {
3207         float median_ = PVValHelper::getMedian(*iterator).value();
3208         float medianErr_ = PVValHelper::getMedian(*iterator).error();
3209         trendPlot->SetBinContent(bin, median_);
3210         trendPlot->SetBinError(bin, medianErr_);
3211         break;
3212       }
3213       case PVValHelper::MAD: {
3214         float mad_ = PVValHelper::getMAD(*iterator).value();
3215         float madErr_ = PVValHelper::getMAD(*iterator).error();
3216         trendPlot->SetBinContent(bin, mad_);
3217         trendPlot->SetBinError(bin, madErr_);
3218         break;
3219       }
3220       default:
3221         edm::LogWarning("PrimaryVertexValidation")
3222             << "fillTrendPlotByIndex() " << fitPar_ << " unknown estimator!" << std::endl;
3223         break;
3224     }
3225
3226     char bincenter[129];
3227     if (plotVar == PVValHelper::eta) {
3228       auto etaBins = theDetails_.trendbins[PVValHelper::eta];
3229       sprintf(bincenter, "%.1f", (etaBins[bin - 1] + etaBins[bin]) / 2.);
3230       trendPlot->GetXaxis()->SetBinLabel(bin, bincenter);
3231     } else if (plotVar == PVValHelper::phi) {
3232       auto phiBins = theDetails_.trendbins[PVValHelper::phi];
3233       sprintf(bincenter, "%.f", (phiBins[bin - 1] + phiBins[bin]) / 2.);
3234       trendPlot->GetXaxis()->SetBinLabel(bin, bincenter);
3235     } else {
3236       /// FIXME DO SOMETHING HERE
3237       //edm::LogWarning("PrimaryVertexValidation")<<"fillTrendPlotByIndex() "<< plotVar <<" unknown track parameter!"<<std::endl;
3238     }
3239   }
3240 }
3241
3242 //*************************************************************
3243 void PrimaryVertexValidation::fillMap(TH2F* trendMap,
3244                                       TH1F* residualsMapPlot[100][100],
3245                                       PVValHelper::estimator fitPar_,
3246                                       const int nXBins_,
3247                                       const int nYBins_)
3248 //*************************************************************
3249 {
3250   for (int i = 0; i < nYBins_; ++i) {
3251     char phibincenter[129];
3252     auto phiBins = theDetails_.trendbins[PVValHelper::phi];
3253     sprintf(phibincenter, "%.f", (phiBins[i] + phiBins[i + 1]) / 2.);
3254
3255     if (nXBins_ == nYBins_) {
3256       trendMap->GetYaxis()->SetBinLabel(i + 1, phibincenter);
3257     }
3258
3259     for (int j = 0; j < nXBins_; ++j) {
3260       char etabincenter[129];
3261       auto etaBins = theDetails_.trendbins[PVValHelper::eta];
3262       sprintf(etabincenter, "%.1f", (etaBins[j] + etaBins[j + 1]) / 2.);
3263
3264       if (i == 0) {
3265         if (nXBins_ == nYBins_) {
3266           trendMap->GetXaxis()->SetBinLabel(j + 1, etabincenter);
3267         }
3268       }
3269
3270       switch (fitPar_) {
3271         case PVValHelper::MEAN: {
3272           float mean_ = PVValHelper::fitResiduals(residualsMapPlot[i][j]).first.value();
3273           float meanErr_ = PVValHelper::fitResiduals(residualsMapPlot[i][j]).first.error();
3274           trendMap->SetBinContent(j + 1, i + 1, mean_);
3275           trendMap->SetBinError(j + 1, i + 1, meanErr_);
3276           break;
3277         }
3278         case PVValHelper::WIDTH: {
3279           float width_ = PVValHelper::fitResiduals(residualsMapPlot[i][j]).second.value();
3280           float widthErr_ = PVValHelper::fitResiduals(residualsMapPlot[i][j]).second.error();
3281           trendMap->SetBinContent(j + 1, i + 1, width_);
3282           trendMap->SetBinError(j + 1, i + 1, widthErr_);
3283           break;
3284         }
3285         case PVValHelper::MEDIAN: {
3286           float median_ = PVValHelper::getMedian(residualsMapPlot[i][j]).value();
3287           float medianErr_ = PVValHelper::getMedian(residualsMapPlot[i][j]).error();
3288           trendMap->SetBinContent(j + 1, i + 1, median_);
3289           trendMap->SetBinError(j + 1, i + 1, medianErr_);
3290           break;
3291         }
3292         case PVValHelper::MAD: {
3293           float mad_ = PVValHelper::getMAD(residualsMapPlot[i][j]).value();
3294           float madErr_ = PVValHelper::getMAD(residualsMapPlot[i][j]).error();
3295           trendMap->SetBinContent(j + 1, i + 1, mad_);
3296           trendMap->SetBinError(j + 1, i + 1, madErr_);
3297           break;
3298         }
3299         default:
3300           edm::LogWarning("PrimaryVertexValidation:") << " fillMap() " << fitPar_ << " unknown estimator!" << std::endl;
3301       }
3302     }  // closes loop on eta bins
3303   }    // cloeses loop on phi bins
3304 }
3305
3306 //*************************************************************
3307 bool PrimaryVertexValidation::vtxSort(const reco::Vertex& a, const reco::Vertex& b)
3308 //*************************************************************
3309 {
3310   if (a.tracksSize() != b.tracksSize())
3311     return a.tracksSize() > b.tracksSize() ? true : false;
3312   else
3313     return a.chi2() < b.chi2() ? true : false;
3314 }
3315
3316 //*************************************************************
3317 bool PrimaryVertexValidation::passesTrackCuts(const reco::Track& track,
3318                                               const reco::Vertex& vertex,
3319                                               const std::string& qualityString_,
3320                                               double dxyErrMax_,
3321                                               double dzErrMax_,
3322                                               double ptErrMax_)
3323 //*************************************************************
3324 {
3325   math::XYZPoint vtxPoint(0.0, 0.0, 0.0);
3326   double vzErr = 0.0, vxErr = 0.0, vyErr = 0.0;
3327   vtxPoint = vertex.position();
3328   vzErr = vertex.zError();
3329   vxErr = vertex.xError();
3330   vyErr = vertex.yError();
3331
3332   double dxy = 0.0, dz = 0.0, dxysigma = 0.0, dzsigma = 0.0;
3333   dxy = track.dxy(vtxPoint);
3334   dz = track.dz(vtxPoint);
3335   dxysigma = sqrt(track.d0Error() * track.d0Error() + vxErr * vyErr);
3336   dzsigma = sqrt(track.dzError() * track.dzError() + vzErr * vzErr);
3337
3338   if (track.quality(reco::TrackBase::qualityByName(qualityString_)) != 1)
3339     return false;
3340   if (std::abs(dxy / dxysigma) > dxyErrMax_)
3341     return false;
3342   if (std::abs(dz / dzsigma) > dzErrMax_)
3343     return false;
3344   if (track.ptError() / track.pt() > ptErrMax_)
3345     return false;
3346
3347   return true;
3348 }
3349
3350 //*************************************************************
3351 std::map<std::string, TH1*> PrimaryVertexValidation::bookVertexHistograms(const TFileDirectory& dir)
3352 //*************************************************************
3353 {
3354   TH1F::SetDefaultSumw2(kTRUE);
3355
3356   std::map<std::string, TH1*> h;
3357
3358   // histograms of track quality (Data and MC)
3359   std::string types[] = {"all", "sel"};
3360   for (const auto& type : types) {
3361     h["pseudorapidity_" + type] =
3362         dir.make<TH1F>(("rapidity_" + type).c_str(), "track pseudorapidity; track #eta; tracks", 100, -3., 3.);
3363     h["z0_" + type] = dir.make<TH1F>(("z0_" + type).c_str(), "track z_{0};track z_{0} (cm);tracks", 80, -40., 40.);
3364     h["phi_" + type] = dir.make<TH1F>(("phi_" + type).c_str(), "track #phi; track #phi;tracks", 80, -M_PI, M_PI);
3365     h["eta_" + type] = dir.make<TH1F>(("eta_" + type).c_str(), "track #eta; track #eta;tracks", 80, -4., 4.);
3366     h["pt_" + type] = dir.make<TH1F>(("pt_" + type).c_str(), "track p_{T}; track p_{T} [GeV];tracks", 100, 0., 20.);
3367     h["p_" + type] = dir.make<TH1F>(("p_" + type).c_str(), "track p; track p [GeV];tracks", 100, 0., 20.);
3368     h["found_" + type] =
3369         dir.make<TH1F>(("found_" + type).c_str(), "n. found hits;n^{found}_{hits};tracks", 30, 0., 30.);
3370     h["lost_" + type] = dir.make<TH1F>(("lost_" + type).c_str(), "n. lost hits;n^{lost}_{hits};tracks", 20, 0., 20.);
3371     h["nchi2_" + type] =
3372         dir.make<TH1F>(("nchi2_" + type).c_str(), "normalized track #chi^{2};track #chi^{2}/ndf;tracks", 100, 0., 20.);
3373     h["rstart_" + type] = dir.make<TH1F>(
3374         ("rstart_" + type).c_str(), "track start radius; track innermost radius r (cm);tracks", 100, 0., 20.);
3375     h["expectedInner_" + type] = dir.make<TH1F>(
3376         ("expectedInner_" + type).c_str(), "n. expected inner hits;n^{expected}_{inner};tracks", 10, 0., 10.);
3377     h["expectedOuter_" + type] = dir.make<TH1F>(
3378         ("expectedOuter_" + type).c_str(), "n. expected outer hits;n^{expected}_{outer};tracks ", 10, 0., 10.);
3379     h["logtresxy_" + type] =
3380         dir.make<TH1F>(("logtresxy_" + type).c_str(),
3381                        "log10(track r-#phi resolution/#mum);log10(track r-#phi resolution/#mum);tracks",
3382                        100,
3383                        0.,
3384                        5.);
3385     h["logtresz_" + type] = dir.make<TH1F>(("logtresz_" + type).c_str(),
3386                                            "log10(track z resolution/#mum);log10(track z resolution/#mum);tracks",
3387                                            100,
3388                                            0.,
3389                                            5.);
3390     h["tpullxy_" + type] =
3391         dir.make<TH1F>(("tpullxy_" + type).c_str(), "track r-#phi pull;pull_{r-#phi};tracks", 100, -10., 10.);
3392     h["tpullz_" + type] =
3393         dir.make<TH1F>(("tpullz_" + type).c_str(), "track r-z pull;pull_{r-z};tracks", 100, -50., 50.);
3394     h["tlogDCAxy_" + type] = dir.make<TH1F>(
3395         ("tlogDCAxy_" + type).c_str(), "track log_{10}(DCA_{r-#phi});track log_{10}(DCA_{r-#phi});tracks", 200, -5., 3.);
3396     h["tlogDCAz_" + type] = dir.make<TH1F>(
3397         ("tlogDCAz_" + type).c_str(), "track log_{10}(DCA_{r-z});track log_{10}(DCA_{r-z});tracks", 200, -5., 5.);
3398     h["lvseta_" + type] = dir.make<TH2F>(
3399         ("lvseta_" + type).c_str(), "cluster length vs #eta;track #eta;cluster length", 60, -3., 3., 20, 0., 20);
3400     h["lvstanlambda_" + type] = dir.make<TH2F>(("lvstanlambda_" + type).c_str(),
3401                                                "cluster length vs tan #lambda; tan#lambda;cluster length",
3402                                                60,
3403                                                -6.,
3404                                                6.,
3405                                                20,
3406                                                0.,
3407                                                20);
3408     h["restrkz_" + type] =
3409         dir.make<TH1F>(("restrkz_" + type).c_str(), "z-residuals (track vs vertex);res_{z} (cm);tracks", 200, -5., 5.);
3410     h["restrkzvsphi_" + type] = dir.make<TH2F>(("restrkzvsphi_" + type).c_str(),
3411                                                "z-residuals (track - vertex) vs track #phi;track #phi;res_{z} (cm)",
3412                                                12,
3413                                                -M_PI,
3414                                                M_PI,
3415                                                100,
3416                                                -0.5,
3417                                                0.5);
3418     h["restrkzvseta_" + type] = dir.make<TH2F>(("restrkzvseta_" + type).c_str(),
3419                                                "z-residuals (track - vertex) vs track #eta;track #eta;res_{z} (cm)",
3420                                                12,
3421                                                -3.,
3422                                                3.,
3423                                                200,
3424                                                -0.5,
3425                                                0.5);
3426     h["pulltrkzvsphi_" + type] =
3427         dir.make<TH2F>(("pulltrkzvsphi_" + type).c_str(),
3428                        "normalized z-residuals (track - vertex) vs track #phi;track #phi;res_{z}/#sigma_{res_{z}}",
3429                        12,
3430                        -M_PI,
3431                        M_PI,
3432                        100,
3433                        -5.,
3434                        5.);
3435     h["pulltrkzvseta_" + type] =
3436         dir.make<TH2F>(("pulltrkzvseta_" + type).c_str(),
3437                        "normalized z-residuals (track - vertex) vs track #eta;track #eta;res_{z}/#sigma_{res_{z}}",
3438                        12,
3439                        -3.,
3440                        3.,
3441                        100,
3442                        -5.,
3443                        5.);
3444     h["pulltrkz_" + type] = dir.make<TH1F>(("pulltrkz_" + type).c_str(),
3445                                            "normalized z-residuals (track vs vertex);res_{z}/#sigma_{res_{z}};tracks",
3446                                            100,
3447                                            -5.,
3448                                            5.);
3449     h["sigmatrkz0_" + type] = dir.make<TH1F>(
3450         ("sigmatrkz0_" + type).c_str(), "z-resolution (excluding beam);#sigma^{trk}_{z_{0}} (cm);tracks", 100, 0., 5.);
3451     h["sigmatrkz_" + type] = dir.make<TH1F>(
3452         ("sigmatrkz_" + type).c_str(), "z-resolution (including beam);#sigma^{trk}_{z} (cm);tracks", 100, 0., 5.);
3453     h["nbarrelhits_" + type] = dir.make<TH1F>(
3454         ("nbarrelhits_" + type).c_str(), "number of pixel barrel hits;n. hits Barrel Pixel;tracks", 10, 0., 10.);
3455     h["nbarrelLayers_" + type] = dir.make<TH1F>(
3456         ("nbarrelLayers_" + type).c_str(), "number of pixel barrel layers;n. layers Barrel Pixel;tracks", 10, 0., 10.);
3457     h["nPxLayers_" + type] = dir.make<TH1F>(
3458         ("nPxLayers_" + type).c_str(), "number of pixel layers (barrel+endcap);n. Pixel layers;tracks", 10, 0., 10.);
3459     h["nSiLayers_" + type] =
3460         dir.make<TH1F>(("nSiLayers_" + type).c_str(), "number of Tracker layers;n. Tracker layers;tracks", 20, 0., 20.);
3461     h["trackAlgo_" + type] =
3462         dir.make<TH1F>(("trackAlgo_" + type).c_str(), "track algorithm;track algo;tracks", 30, 0., 30.);
3463     h["trackQuality_" + type] =
3464         dir.make<TH1F>(("trackQuality_" + type).c_str(), "track quality;track quality;tracks", 7, -1., 6.);
3465   }
3466
3467   return h;
3468 }
3469
3470 //*************************************************************
3471 // Generic booker function
3472 //*************************************************************
3473 std::vector<TH1F*> PrimaryVertexValidation::bookResidualsHistogram(const TFileDirectory& dir,
3474                                                                    unsigned int theNOfBins,
3475                                                                    PVValHelper::residualType resType,
3476                                                                    PVValHelper::plotVariable varType,
3477                                                                    bool isNormalized) {
3478   TH1F::SetDefaultSumw2(kTRUE);
3479
3480   auto hash = std::make_pair(resType, varType);
3481
3482   double down = theDetails_.range[hash].first;
3483   double up = theDetails_.range[hash].second;
3484
3485   if (isNormalized) {
3486     up = up / 100.;
3487     down = down / 100.;
3488   }
3489
3490   std::vector<TH1F*> h;
3491   h.reserve(theNOfBins);
3492
3493   if (theNOfBins == 0) {
3494     edm::LogError("PrimaryVertexValidation")
3495         << "bookResidualsHistogram() The number of bins cannot be identically 0" << std::endl;
3496     assert(false);
3497   }
3498
3499   std::string s_resType = std::get<0>(PVValHelper::getTypeString(resType));
3500   std::string s_varType = std::get<0>(PVValHelper::getVarString(varType));
3501
3502   std::string t_resType = std::get<1>(PVValHelper::getTypeString(resType));
3503   std::string t_varType = std::get<1>(PVValHelper::getVarString(varType));
3504   std::string units = std::get<2>(PVValHelper::getTypeString(resType));
3505
3506   for (unsigned int i = 0; i < theNOfBins; i++) {
3507     TString title = (varType == PVValHelper::phi || varType == PVValHelper::eta)
3508                         ? Form("%s vs %s - bin %i (%f < %s < %f);%s %s;tracks",
3509                                t_resType.c_str(),
3510                                t_varType.c_str(),
3511                                i,
3512                                theDetails_.trendbins[varType][i],
3513                                t_varType.c_str(),
3514                                theDetails_.trendbins[varType][i + 1],
3515                                t_resType.c_str(),
3516                                units.c_str())
3517                         : Form("%s vs %s - bin %i;%s %s;tracks",
3518                                t_resType.c_str(),
3519                                t_varType.c_str(),
3520                                i,
3521                                t_resType.c_str(),
3522                                units.c_str());
3523
3524     TH1F* htemp = dir.make<TH1F>(
3525         Form("histo_%s_%s_plot%i", s_resType.c_str(), s_varType.c_str(), i),
3526         //Form("%s vs %s - bin %i;%s %s;tracks",t_resType.c_str(),t_varType.c_str(),i,t_resType.c_str(),units.c_str()),
3527         title.Data(),
3528         theDetails_.histobins,
3529         down,
3530         up);
3531     h.push_back(htemp);
3532   }
3533
3534   return h;
3535 }
3536
3537 //*************************************************************
3538 void PrimaryVertexValidation::fillTrackHistos(std::map<std::string, TH1*>& h,
3539                                               const std::string& ttype,
3540                                               const reco::TransientTrack* tt,
3541                                               const reco::Vertex& v,
3542                                               const reco::BeamSpot& beamSpot,
3543                                               double fBfield_)
3544 //*************************************************************
3545 {
3546   using namespace reco;
3547
3548   PVValHelper::fill(h, "pseudorapidity_" + ttype, tt->track().eta());
3549   PVValHelper::fill(h, "z0_" + ttype, tt->track().vz());
3550   PVValHelper::fill(h, "phi_" + ttype, tt->track().phi());
3551   PVValHelper::fill(h, "eta_" + ttype, tt->track().eta());
3552   PVValHelper::fill(h, "pt_" + ttype, tt->track().pt());
3553   PVValHelper::fill(h, "p_" + ttype, tt->track().p());
3554   PVValHelper::fill(h, "found_" + ttype, tt->track().found());
3555   PVValHelper::fill(h, "lost_" + ttype, tt->track().lost());
3556   PVValHelper::fill(h, "nchi2_" + ttype, tt->track().normalizedChi2());
3557   PVValHelper::fill(h, "rstart_" + ttype, (tt->track().innerPosition()).Rho());
3558
3559   double d0Error = tt->track().d0Error();
3560   double d0 = tt->track().dxy(beamSpot.position());
3561   double dz = tt->track().dz(beamSpot.position());
3562   if (d0Error > 0) {
3563     PVValHelper::fill(h, "logtresxy_" + ttype, log(d0Error / 0.0001) / log(10.));
3564     PVValHelper::fill(h, "tpullxy_" + ttype, d0 / d0Error);
3565     PVValHelper::fill(h, "tlogDCAxy_" + ttype, log(std::abs(d0 / d0Error)));
3566   }
3567   //double z0=tt->track().vz();
3568   double dzError = tt->track().dzError();
3569   if (dzError > 0) {
3570     PVValHelper::fill(h, "logtresz_" + ttype, log(dzError / 0.0001) / log(10.));
3571     PVValHelper::fill(h, "tpullz_" + ttype, dz / dzError);
3572     PVValHelper::fill(h, "tlogDCAz_" + ttype, log(std::abs(dz / dzError)));
3573   }
3574
3575   //
3576   double wxy2_ = pow(beamSpot.BeamWidthX(), 2) + pow(beamSpot.BeamWidthY(), 2);
3577
3578   PVValHelper::fill(
3579       h, "sigmatrkz_" + ttype, sqrt(pow(tt->track().dzError(), 2) + wxy2_ / pow(tan(tt->track().theta()), 2)));
3580   PVValHelper::fill(h, "sigmatrkz0_" + ttype, tt->track().dzError());
3581
3582   // track vs vertex
3583   if (v.isValid()) {  // && (v.ndof()<10.)) {
3584     // emulate clusterizer input
3585     //const TransientTrack & tt = theB_->build(&t); wrong !!!!
3586     //reco::TransientTrack tt = theB_->build(&t);
3587     //ttt->track().setBeamSpot(beamSpot); // need the setBeamSpot !
3588     double z = (tt->stateAtBeamLine().trackStateAtPCA()).position().z();
3589     double tantheta = tan((tt->stateAtBeamLine().trackStateAtPCA()).momentum().theta());
3590     double dz2 = pow(tt->track().dzError(), 2) + wxy2_ / pow(tantheta, 2);
3591
3592     PVValHelper::fill(h, "restrkz_" + ttype, z - v.position().z());
3593     PVValHelper::fill(h, "restrkzvsphi_" + ttype, tt->track().phi(), z - v.position().z());
3594     PVValHelper::fill(h, "restrkzvseta_" + ttype, tt->track().eta(), z - v.position().z());
3595     PVValHelper::fill(h, "pulltrkzvsphi_" + ttype, tt->track().phi(), (z - v.position().z()) / sqrt(dz2));
3596     PVValHelper::fill(h, "pulltrkzvseta_" + ttype, tt->track().eta(), (z - v.position().z()) / sqrt(dz2));
3597
3598     PVValHelper::fill(h, "pulltrkz_" + ttype, (z - v.position().z()) / sqrt(dz2));
3599
3600     double x1 = tt->track().vx() - beamSpot.x0();
3601     double y1 = tt->track().vy() - beamSpot.y0();
3602
3603     double kappa = -0.002998 * fBfield_ * tt->track().qoverp() / cos(tt->track().theta());
3604     double D0 = x1 * sin(tt->track().phi()) - y1 * cos(tt->track().phi()) - 0.5 * kappa * (x1 * x1 + y1 * y1);
3605     double q = sqrt(1. - 2. * kappa * D0);
3606     double s0 = (x1 * cos(tt->track().phi()) + y1 * sin(tt->track().phi())) / q;
3607     // double s1;
3608     if (std::abs(kappa * s0) > 0.001) {
3609       //s1=asin(kappa*s0)/kappa;
3610     } else {
3611       //double ks02=(kappa*s0)*(kappa*s0);
3612       //s1=s0*(1.+ks02/6.+3./40.*ks02*ks02+5./112.*pow(ks02,3));
3613     }
3614     // sp.ddcap=-2.*D0/(1.+q);
3615     //double zdcap=tt->track().vz()-s1/tan(tt->track().theta());
3616   }
3617   //
3618
3619   // collect some info on hits and clusters
3620   PVValHelper::fill(h, "nbarrelLayers_" + ttype, tt->track().hitPattern().pixelBarrelLayersWithMeasurement());
3621   PVValHelper::fill(h, "nPxLayers_" + ttype, tt->track().hitPattern().pixelLayersWithMeasurement());
3622   PVValHelper::fill(h, "nSiLayers_" + ttype, tt->track().hitPattern().trackerLayersWithMeasurement());
3623   PVValHelper::fill(
3624       h, "expectedInner_" + ttype, tt->track().hitPattern().numberOfLostHits(HitPattern::MISSING_INNER_HITS));
3625   PVValHelper::fill(
3626       h, "expectedOuter_" + ttype, tt->track().hitPattern().numberOfLostHits(HitPattern::MISSING_OUTER_HITS));
3627   PVValHelper::fill(h, "trackAlgo_" + ttype, tt->track().algo());
3628   PVValHelper::fill(h, "trackQuality_" + ttype, tt->track().qualityMask());
3629
3630   //
3631   int longesthit = 0, nbarrel = 0;
3632   for (auto const& hit : tt->track().recHits()) {
3633     if (hit->isValid() && hit->geographicalId().det() == DetId::Tracker) {
3634       bool barrel = DetId(hit->geographicalId()).subdetId() == static_cast<int>(PixelSubdetector::PixelBarrel);
3635       //bool endcap = DetId::DetId(hit->geographicalId()).subdetId() == static_cast<int>(PixelSubdetector::PixelEndcap);
3636       if (barrel) {
3637         const SiPixelRecHit* pixhit = dynamic_cast<const SiPixelRecHit*>(&(*hit));
3638         edm::Ref<edmNew::DetSetVector<SiPixelCluster>, SiPixelCluster> const& clust = (*pixhit).cluster();
3639         if (clust.isNonnull()) {
3640           nbarrel++;
3641           if (clust->sizeY() - longesthit > 0)
3642             longesthit = clust->sizeY();
3643           if (clust->sizeY() > 20.) {
3644             PVValHelper::fill(h, "lvseta_" + ttype, tt->track().eta(), 19.9);
3645             PVValHelper::fill(h, "lvstanlambda_" + ttype, tan(tt->track().lambda()), 19.9);
3646           } else {
3647             PVValHelper::fill(h, "lvseta_" + ttype, tt->track().eta(), float(clust->sizeY()));
3648             PVValHelper::fill(h, "lvstanlambda_" + ttype, tan(tt->track().lambda()), float(clust->sizeY()));
3649           }
3650         }
3651       }
3652     }
3653   }
3654   PVValHelper::fill(h, "nbarrelhits_" + ttype, float(nbarrel));
3655   //-------------------------------------------------------------------
3656 }
3657
3658 void PrimaryVertexValidation::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
3659   edm::ParameterSetDescription desc;
3660   desc.setComment("Validates alignment payloads by evaluating unbiased track paramter resisuals to vertices");
3661
3662   // PV Validation specific
3663
3664   desc.addUntracked<int>("compressionSettings", -1);
3665   desc.add<bool>("storeNtuple", false);
3666   desc.add<bool>("isLightNtuple", true);
3667   desc.add<bool>("useTracksFromRecoVtx", false);
3668   desc.addUntracked<double>("vertexZMax", 99);
3669   desc.addUntracked<double>("intLumi", 0.);
3670   desc.add<bool>("askFirstLayerHit", false);
3671   desc.addUntracked<bool>("doBPix", true);
3672   desc.addUntracked<bool>("doFPix", true);
3673   desc.addUntracked<double>("probePt", 0.);
3674   desc.addUntracked<double>("probeP", 0.);
3675   desc.addUntracked<double>("probeEta", 2.4);
3676   desc.addUntracked<double>("probeNHits", 0.);
3677   desc.addUntracked<int>("numberOfBins", 24);
3678   desc.addUntracked<double>("minPt", 1.);
3679   desc.addUntracked<double>("maxPt", 20.);
3680   desc.add<bool>("Debug", false);
3681   desc.addUntracked<bool>("runControl", false);
3682   desc.addUntracked<bool>("forceBeamSpot", false);
3683
3684   std::vector<unsigned int> defaultRuns;
3685   defaultRuns.push_back(0);
3686   desc.addUntracked<std::vector<unsigned int>>("runControlNumber", defaultRuns);
3687
3688   // event sources
3689
3690   desc.add<edm::InputTag>("TrackCollectionTag", edm::InputTag("ALCARECOTkAlMinBias"));
3691   desc.add<edm::InputTag>("VertexCollectionTag", edm::InputTag("offlinePrimaryVertices"));
3692   desc.add<edm::InputTag>("BeamSpotTag", edm::InputTag("offlineBeamSpot"));
3693
3694   // track filtering
3695   edm::ParameterSetDescription psd0;
3696   TrackFilterForPVFinding::fillPSetDescription(psd0);
3697   psd0.add<int>("numTracksThreshold", 0);  // HI only
3698   desc.add<edm::ParameterSetDescription>("TkFilterParameters", psd0);
3699
3700   // PV Clusterization
3701   {
3702     edm::ParameterSetDescription psd0;
3703     {
3704       edm::ParameterSetDescription psd1;
3705       DAClusterizerInZ_vect::fillPSetDescription(psd1);
3706       psd0.add<edm::ParameterSetDescription>("TkDAClusParameters", psd1);
3707
3708       edm::ParameterSetDescription psd2;
3709       GapClusterizerInZ::fillPSetDescription(psd2);
3710       psd0.add<edm::ParameterSetDescription>("TkGapClusParameters", psd2);
3711     }
3712     psd0.add<std::string>("algorithm", "DA_vect");
3713     desc.add<edm::ParameterSetDescription>("TkClusParameters", psd0);
3714   }
3715
3716   descriptions.add("primaryVertexValidation", desc);
3717 }
3718
3719 //define this as a plug-in
3720 DEFINE_FWK_MODULE(PrimaryVertexValidation);