Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-11-07 23:57:44

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