Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-11-25 03:29:11

0001 // -*- C++ -*-
0002 //
0003 // Package:    Alignment/OfflineValidation
0004 // Class:      SplitVertexResolution
0005 //
0006 /**\class SplitVertexResolution SplitVertexResolution.cc Alignment/OfflineValidation/plugins/SplitVertexResolution.cc
0007 
0008 */
0009 //
0010 // Original Author:  Marco Musich
0011 //         Created:  Mon, 13 Jun 2016 15:07:11 GMT
0012 //
0013 //
0014 
0015 // system include files
0016 #include <memory>
0017 #include <algorithm>  // std::sort
0018 #include <vector>     // std::vector
0019 #include <chrono>
0020 #include <iostream>
0021 #include <random>
0022 #include <fmt/printf.h>
0023 #include <boost/range/adaptor/indexed.hpp>
0024 
0025 // ROOT include files
0026 #include "TTree.h"
0027 #include "TProfile.h"
0028 #include "TF1.h"
0029 #include "TMath.h"
0030 
0031 // user include files
0032 #include "FWCore/Framework/interface/Frameworkfwd.h"
0033 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0034 #include "FWCore/Utilities/interface/isFinite.h"
0035 
0036 #include "FWCore/Framework/interface/Event.h"
0037 #include "FWCore/Framework/interface/EventSetup.h"
0038 #include "FWCore/Framework/interface/MakerMacros.h"
0039 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0040 #include "FWCore/ServiceRegistry/interface/Service.h"
0041 #include "FWCore/Utilities/interface/InputTag.h"
0042 #include "FWCore/Utilities/interface/ESInputTag.h"
0043 
0044 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0045 
0046 #include "DataFormats/TrackReco/interface/Track.h"
0047 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0048 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0049 #include "DataFormats/VertexReco/interface/Vertex.h"
0050 
0051 #include "CondFormats/RunInfo/interface/RunInfo.h"
0052 #include "DataFormats/Common/interface/TriggerResults.h"
0053 #include "FWCore/Common/interface/TriggerNames.h"
0054 
0055 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0056 #include "RecoVertex/AdaptiveVertexFit/interface/AdaptiveVertexFitter.h"
0057 #include "RecoVertex/VertexTools/interface/VertexDistance3D.h"
0058 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0059 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0060 
0061 #include "Alignment/OfflineValidation/interface/PVValidationHelpers.h"
0062 #include "Alignment/OfflineValidation/interface/pvTree.h"
0063 
0064 //
0065 // useful code
0066 //
0067 
0068 namespace statmode {
0069   using fitParams = std::pair<Measurement1D, Measurement1D>;
0070 }
0071 
0072 //
0073 // class declaration
0074 //
0075 
0076 class SplitVertexResolution : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0077 public:
0078   explicit SplitVertexResolution(const edm::ParameterSet&);
0079   ~SplitVertexResolution() override;
0080 
0081   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0082   static bool mysorter(reco::Track i, reco::Track j) { return (i.pt() > j.pt()); }
0083 
0084 private:
0085   void beginJob() override;
0086   void beginRun(edm::Run const& iEvent, edm::EventSetup const&) override;
0087   virtual void beginEvent() final;
0088   void analyze(const edm::Event&, const edm::EventSetup&) override;
0089   void endJob() override;
0090   void endRun(edm::Run const&, edm::EventSetup const&) override{};
0091 
0092   template <std::size_t SIZE>
0093   bool checkBinOrdering(std::array<float, SIZE>& bins);
0094   std::vector<TH1F*> bookResidualsHistogram(TFileDirectory dir,
0095                                             unsigned int theNOfBins,
0096                                             TString resType,
0097                                             TString varType);
0098 
0099   void fillTrendPlotByIndex(TH1F* trendPlot, std::vector<TH1F*>& h, PVValHelper::estimator fitPar_);
0100   statmode::fitParams fitResiduals(TH1* hist, bool singleTime = false);
0101   statmode::fitParams fitResiduals_v0(TH1* hist);
0102 
0103   std::pair<long long, long long> getRunTime(const edm::EventSetup& iSetup) const;
0104 
0105   // counters
0106   int ievt;
0107   int itrks;
0108 
0109   // compression settings
0110   const int compressionSettings_;
0111 
0112   // switch to keep the ntuple
0113   bool storeNtuple_;
0114 
0115   // storing integrated luminosity
0116   double intLumi_;
0117   bool debug_;
0118 
0119   edm::InputTag pvsTag_;
0120   edm::EDGetTokenT<reco::VertexCollection> pvsToken_;
0121 
0122   edm::InputTag tracksTag_;
0123   edm::EDGetTokenT<reco::TrackCollection> tracksToken_;
0124 
0125   edm::InputTag triggerResultsTag_ = edm::InputTag("TriggerResults", "", "HLT");  //InputTag tag("TriggerResults");
0126   edm::EDGetTokenT<edm::TriggerResults> triggerResultsToken_;
0127 
0128   // ES Tokens
0129   edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> transientTrackBuilderToken_;
0130   edm::ESGetToken<RunInfo, RunInfoRcd> runInfoToken_;
0131 
0132   double minVtxNdf_;
0133   double minVtxWgt_;
0134 
0135   // checks on the run to be processed
0136 
0137   bool runControl_;
0138   std::vector<unsigned int> runControlNumbers_;
0139 
0140   static constexpr double cmToUm = 10000.;
0141 
0142   edm::Service<TFileService> outfile_;
0143 
0144   TH1F* h_lumiFromConfig;
0145   TH1I* h_runFromConfig;
0146 
0147   std::map<unsigned int, std::pair<long long, long long>> runNumbersTimesLog_;
0148   TH1I* h_runStartTimes;
0149   TH1I* h_runEndTimes;
0150 
0151   TH1F* h_diffX;
0152   TH1F* h_diffY;
0153   TH1F* h_diffZ;
0154 
0155   TH1F* h_OrigVertexErrX;
0156   TH1F* h_OrigVertexErrY;
0157   TH1F* h_OrigVertexErrZ;
0158 
0159   TH1F* h_errX;
0160   TH1F* h_errY;
0161   TH1F* h_errZ;
0162 
0163   TH1F* h_pullX;
0164   TH1F* h_pullY;
0165   TH1F* h_pullZ;
0166 
0167   TH1F* h_ntrks;
0168   TH1F* h_sumPt;
0169   TH1F* h_avgSumPt;
0170 
0171   TH1F* h_sumPt1;
0172   TH1F* h_sumPt2;
0173 
0174   TH1F* h_wTrks1;
0175   TH1F* h_wTrks2;
0176 
0177   TH1F* h_minWTrks1;
0178   TH1F* h_minWTrks2;
0179 
0180   TH1F* h_PVCL_subVtx1;
0181   TH1F* h_PVCL_subVtx2;
0182 
0183   TH1F* h_runNumber;
0184 
0185   TH1I* h_nOfflineVertices;
0186   TH1I* h_nVertices;
0187   TH1I* h_nNonFakeVertices;
0188   TH1I* h_nFinalVertices;
0189 
0190   // trigger results
0191 
0192   TH1D* tksByTrigger_;
0193   TH1D* evtsByTrigger_;
0194 
0195   // resolutions
0196 
0197   std::vector<TH1F*> h_resolX_sumPt_;
0198   std::vector<TH1F*> h_resolY_sumPt_;
0199   std::vector<TH1F*> h_resolZ_sumPt_;
0200 
0201   std::vector<TH1F*> h_resolX_Ntracks_;
0202   std::vector<TH1F*> h_resolY_Ntracks_;
0203   std::vector<TH1F*> h_resolZ_Ntracks_;
0204 
0205   std::vector<TH1F*> h_resolX_Nvtx_;
0206   std::vector<TH1F*> h_resolY_Nvtx_;
0207   std::vector<TH1F*> h_resolZ_Nvtx_;
0208 
0209   TH1F* p_resolX_vsSumPt;
0210   TH1F* p_resolY_vsSumPt;
0211   TH1F* p_resolZ_vsSumPt;
0212 
0213   TH1F* p_resolX_vsNtracks;
0214   TH1F* p_resolY_vsNtracks;
0215   TH1F* p_resolZ_vsNtracks;
0216 
0217   TH1F* p_resolX_vsNvtx;
0218   TH1F* p_resolY_vsNvtx;
0219   TH1F* p_resolZ_vsNvtx;
0220 
0221   // pulls
0222   std::vector<TH1F*> h_pullX_sumPt_;
0223   std::vector<TH1F*> h_pullY_sumPt_;
0224   std::vector<TH1F*> h_pullZ_sumPt_;
0225 
0226   std::vector<TH1F*> h_pullX_Ntracks_;
0227   std::vector<TH1F*> h_pullY_Ntracks_;
0228   std::vector<TH1F*> h_pullZ_Ntracks_;
0229 
0230   std::vector<TH1F*> h_pullX_Nvtx_;
0231   std::vector<TH1F*> h_pullY_Nvtx_;
0232   std::vector<TH1F*> h_pullZ_Nvtx_;
0233 
0234   TH1F* p_pullX_vsSumPt;
0235   TH1F* p_pullY_vsSumPt;
0236   TH1F* p_pullZ_vsSumPt;
0237 
0238   TH1F* p_pullX_vsNtracks;
0239   TH1F* p_pullY_vsNtracks;
0240   TH1F* p_pullZ_vsNtracks;
0241 
0242   TH1F* p_pullX_vsNvtx;
0243   TH1F* p_pullY_vsNvtx;
0244   TH1F* p_pullZ_vsNvtx;
0245 
0246   TH1F* h_profileBinnings;
0247 
0248   std::mt19937 engine_;
0249 
0250   pvEvent event_;
0251   TTree* tree_;
0252 
0253   // ----------member data ---------------------------
0254   const float sumpTStartScale_;
0255   const float sumpTEndScale_;
0256   static const int nPtBins_ = 30;
0257   std::array<float, nPtBins_ + 1> mypT_bins_;
0258 
0259   static const int nTrackBins_ = 120;
0260   const double nVisTrackBins_;  // will be configured
0261   std::array<float, nTrackBins_ + 1> myNTrack_bins_;
0262 
0263   static const int nVtxBins_ = 60;
0264   const double nVisVtxBins_;  // will be configured
0265   std::array<float, nVtxBins_ + 1> myNVtx_bins_;
0266 
0267   std::map<std::string, std::pair<int, int>> triggerMap_;
0268 };
0269 
0270 SplitVertexResolution::SplitVertexResolution(const edm::ParameterSet& iConfig)
0271     : compressionSettings_(iConfig.getUntrackedParameter<int>("compressionSettings", -1)),
0272       storeNtuple_(iConfig.getParameter<bool>("storeNtuple")),
0273       intLumi_(iConfig.getUntrackedParameter<double>("intLumi", 0.)),
0274       debug_(iConfig.getUntrackedParameter<bool>("Debug", false)),
0275       pvsTag_(iConfig.getParameter<edm::InputTag>("vtxCollection")),
0276       pvsToken_(consumes<reco::VertexCollection>(pvsTag_)),
0277       tracksTag_(iConfig.getParameter<edm::InputTag>("trackCollection")),
0278       tracksToken_(consumes<reco::TrackCollection>(tracksTag_)),
0279       triggerResultsToken_(consumes<edm::TriggerResults>(triggerResultsTag_)),
0280       transientTrackBuilderToken_(
0281           esConsumes<TransientTrackBuilder, TransientTrackRecord>(edm::ESInputTag("", "TransientTrackBuilder"))),
0282       runInfoToken_(esConsumes<RunInfo, RunInfoRcd, edm::Transition::BeginRun>()),
0283       minVtxNdf_(iConfig.getUntrackedParameter<double>("minVertexNdf")),
0284       minVtxWgt_(iConfig.getUntrackedParameter<double>("minVertexMeanWeight")),
0285       runControl_(iConfig.getUntrackedParameter<bool>("runControl", false)),
0286       // binning
0287       sumpTStartScale_(iConfig.getUntrackedParameter<double>("sumpTStartScale", 1.)),
0288       sumpTEndScale_(iConfig.getUntrackedParameter<double>("sumpTEndScale", 1e3)),
0289       nVisTrackBins_(iConfig.getUntrackedParameter<double>("nTrackBins", 120.)),
0290       nVisVtxBins_(iConfig.getUntrackedParameter<double>("nVtxBins", 60.)) {
0291   usesResource(TFileService::kSharedResource);
0292 
0293   std::vector<unsigned int> defaultRuns;
0294   defaultRuns.push_back(0);
0295   runControlNumbers_ = iConfig.getUntrackedParameter<std::vector<unsigned int>>("runControlNumber", defaultRuns);
0296 
0297   mypT_bins_ = PVValHelper::makeLogBins<float, nPtBins_>(sumpTStartScale_, sumpTEndScale_);
0298 
0299   // IMPORTANT
0300   // first argument is start, second argument is the range so it's [-0.5;nTracksBins_-0.5]
0301   std::vector<float> vect = PVValHelper::generateBins(nTrackBins_ + 1, -0.5, nTrackBins_);
0302   std::copy(vect.begin(), vect.begin() + nTrackBins_ + 1, myNTrack_bins_.begin());
0303 
0304   vect.clear();
0305 
0306   // IMPORTANT
0307   // first argument is start, second argument is the range so it's [-0.5;nVtxBins_-0.5]
0308   vect = PVValHelper::generateBins(nVtxBins_ + 1, -0.5, nVtxBins_);
0309   std::copy(vect.begin(), vect.begin() + nVtxBins_ + 1, myNVtx_bins_.begin());
0310 
0311   if (debug_) {
0312     std::string toOutput = "";
0313     for (const auto& pTbin : mypT_bins_ | boost::adaptors::indexed(1)) {
0314       if (pTbin.index() != 1) {
0315         toOutput += " ";
0316       }
0317       toOutput += fmt::sprintf("%.2f", pTbin.value());
0318       if (pTbin.index() != nPtBins_ + 1) {
0319         toOutput += ",";
0320       }
0321     }
0322     edm::LogVerbatim("SplitVertexResolution") << "sum(pT) bins = [" << toOutput << "] \n";
0323 
0324     toOutput.clear();
0325     for (const auto& tkbin : myNTrack_bins_ | boost::adaptors::indexed(1)) {
0326       if (tkbin.index() != 1) {
0327         toOutput += " ";
0328       }
0329       toOutput += fmt::sprintf("%.1f", tkbin.value());
0330       if (tkbin.index() != nTrackBins_ + 1) {
0331         toOutput += ",";
0332       }
0333     }
0334     edm::LogVerbatim("SplitVertexResolution") << "n. track bins = [" << toOutput << "] \n";
0335 
0336     toOutput.clear();
0337     for (const auto& vtxbin : myNVtx_bins_ | boost::adaptors::indexed(1)) {
0338       if (vtxbin.index() != 1) {
0339         toOutput += " ";
0340       }
0341       toOutput += fmt::sprintf("%.1f", vtxbin.value());
0342       if (vtxbin.index() != nVtxBins_ + 1) {
0343         toOutput += ",";
0344       }
0345     }
0346     edm::LogVerbatim("SplitVertexResolution") << "n. vertices bins = [" << toOutput << "] \n";
0347   }
0348 }
0349 
0350 SplitVertexResolution::~SplitVertexResolution() = default;
0351 
0352 // ------------ method called for each event  ------------
0353 void SplitVertexResolution::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0354   using namespace edm;
0355 
0356   // deterministic seed from the event number
0357   // should not bias the result as the event number is already
0358   // assigned randomly-enough
0359   engine_.seed(iEvent.id().event() + (iEvent.id().luminosityBlock() << 10) + (iEvent.id().run() << 20));
0360 
0361   // first check if the event passes the run control
0362   bool passesRunControl = false;
0363 
0364   if (runControl_) {
0365     for (const auto& runControlNumber : runControlNumbers_) {
0366       if (iEvent.eventAuxiliary().run() == runControlNumber) {
0367         if (debug_) {
0368           edm::LogInfo("SplitVertexResolution")
0369               << " run number: " << iEvent.eventAuxiliary().run() << " keeping run:" << runControlNumber;
0370         }
0371         passesRunControl = true;
0372         break;
0373       }
0374     }
0375     if (!passesRunControl)
0376       return;
0377   }
0378 
0379   // Fill general info
0380   h_runNumber->Fill(iEvent.id().run());
0381 
0382   ievt++;
0383   edm::Handle<edm::TriggerResults> hltresults;
0384   iEvent.getByToken(triggerResultsToken_, hltresults);
0385 
0386   const edm::TriggerNames& triggerNames_ = iEvent.triggerNames(*hltresults);
0387   int ntrigs = hltresults->size();
0388   //const std::vector<std::string>& triggernames = triggerNames_.triggerNames();
0389 
0390   beginEvent();
0391 
0392   // Fill general info
0393   event_.runNumber = iEvent.id().run();
0394   event_.luminosityBlockNumber = iEvent.id().luminosityBlock();
0395   event_.eventNumber = iEvent.id().event();
0396 
0397   TransientTrackBuilder const& theB = iSetup.getData(transientTrackBuilderToken_);
0398 
0399   edm::Handle<reco::VertexCollection> vertices;
0400   iEvent.getByToken(pvsToken_, vertices);
0401   const reco::VertexCollection pvtx = *(vertices.product());
0402 
0403   event_.nVtx = pvtx.size();
0404   int nOfflineVtx = pvtx.size();
0405   h_nOfflineVertices->Fill(nOfflineVtx);
0406 
0407   edm::Handle<reco::TrackCollection> tracks;
0408   iEvent.getByToken(tracksToken_, tracks);
0409   itrks += tracks.product()->size();
0410 
0411   for (int itrig = 0; itrig != ntrigs; ++itrig) {
0412     const std::string& trigName = triggerNames_.triggerName(itrig);
0413     bool accept = hltresults->accept(itrig);
0414     if (accept == 1) {
0415       triggerMap_[trigName].first += 1;
0416       triggerMap_[trigName].second += tracks.product()->size();
0417       // triggerInfo.push_back(pair <string, int> (trigName, accept));
0418     }
0419   }
0420 
0421   int counter = 0;
0422   int noFakecounter = 0;
0423   int goodcounter = 0;
0424 
0425   for (auto pvIt = pvtx.cbegin(); pvIt != pvtx.cend(); ++pvIt) {
0426     reco::Vertex iPV = *pvIt;
0427     counter++;
0428     if (iPV.isFake())
0429       continue;
0430     noFakecounter++;
0431 
0432     // vertex selection as in bs code
0433     if (iPV.ndof() < minVtxNdf_ || (iPV.ndof() + 3.) / iPV.tracksSize() < 2 * minVtxWgt_)
0434       continue;
0435 
0436     goodcounter++;
0437     reco::TrackCollection allTracks;
0438     reco::TrackCollection groupOne, groupTwo;
0439     for (auto trki = iPV.tracks_begin(); trki != iPV.tracks_end(); ++trki) {
0440       if (trki->isNonnull()) {
0441         reco::TrackRef trk_now(tracks, (*trki).key());
0442         allTracks.push_back(*trk_now);
0443       }
0444     }
0445 
0446     if (goodcounter > 1)
0447       continue;
0448 
0449     // order with decreasing pt
0450     std::sort(allTracks.begin(), allTracks.end(), mysorter);
0451 
0452     int ntrks = allTracks.size();
0453     h_ntrks->Fill(ntrks);
0454 
0455     // discard lowest pt track
0456     uint even_ntrks;
0457     ntrks % 2 == 0 ? even_ntrks = ntrks : even_ntrks = ntrks - 1;
0458 
0459     // split into two sets equally populated
0460     for (uint tracksIt = 0; tracksIt < even_ntrks; tracksIt = tracksIt + 2) {
0461       reco::Track firstTrk = allTracks.at(tracksIt);
0462       reco::Track secondTrk = allTracks.at(tracksIt + 1);
0463       auto dis = std::uniform_int_distribution<>(0, 1);  // [0, 1]
0464 
0465       if (dis(engine_) > 0.5) {
0466         groupOne.push_back(firstTrk);
0467         groupTwo.push_back(secondTrk);
0468       } else {
0469         groupOne.push_back(secondTrk);
0470         groupTwo.push_back(firstTrk);
0471       }
0472     }
0473 
0474     if (!(groupOne.size() >= 2 && groupTwo.size() >= 2))
0475       continue;
0476 
0477     h_OrigVertexErrX->Fill(iPV.xError() * cmToUm);
0478     h_OrigVertexErrY->Fill(iPV.yError() * cmToUm);
0479     h_OrigVertexErrZ->Fill(iPV.zError() * cmToUm);
0480 
0481     float sumPt = 0, sumPt1 = 0, sumPt2 = 0, avgSumPt = 0;
0482 
0483     // refit the two sets of tracks
0484     std::vector<reco::TransientTrack> groupOne_ttks;
0485     groupOne_ttks.clear();
0486     for (auto itrk = groupOne.cbegin(); itrk != groupOne.cend(); itrk++) {
0487       reco::TransientTrack tmpTransientTrack = theB.build(*itrk);
0488       groupOne_ttks.push_back(tmpTransientTrack);
0489       sumPt1 += itrk->pt();
0490       sumPt += itrk->pt();
0491     }
0492 
0493     AdaptiveVertexFitter pvFitter;
0494     TransientVertex pvOne = pvFitter.vertex(groupOne_ttks);
0495     if (!pvOne.isValid())
0496       continue;
0497 
0498     reco::Vertex onePV = pvOne;
0499 
0500     std::vector<reco::TransientTrack> groupTwo_ttks;
0501     groupTwo_ttks.clear();
0502     for (auto itrk = groupTwo.cbegin(); itrk != groupTwo.cend(); itrk++) {
0503       reco::TransientTrack tmpTransientTrack = theB.build(*itrk);
0504       groupTwo_ttks.push_back(tmpTransientTrack);
0505       sumPt2 += itrk->pt();
0506       sumPt += itrk->pt();
0507     }
0508 
0509     // average sumPt
0510     avgSumPt = (sumPt1 + sumPt2) / 2.;
0511     h_avgSumPt->Fill(avgSumPt);
0512 
0513     TransientVertex pvTwo = pvFitter.vertex(groupTwo_ttks);
0514     if (!pvTwo.isValid())
0515       continue;
0516 
0517     reco::Vertex twoPV = pvTwo;
0518 
0519     float theminW1 = 1.;
0520     float theminW2 = 1.;
0521     for (auto otrk = pvOne.originalTracks().cbegin(); otrk != pvOne.originalTracks().cend(); ++otrk) {
0522       h_wTrks1->Fill(pvOne.trackWeight(*otrk));
0523       if (pvOne.trackWeight(*otrk) < theminW1) {
0524         theminW1 = pvOne.trackWeight(*otrk);
0525       }
0526     }
0527     for (auto otrk = pvTwo.originalTracks().cbegin(); otrk != pvTwo.originalTracks().end(); ++otrk) {
0528       h_wTrks2->Fill(pvTwo.trackWeight(*otrk));
0529       if (pvTwo.trackWeight(*otrk) < theminW2) {
0530         theminW2 = pvTwo.trackWeight(*otrk);
0531       }
0532     }
0533 
0534     h_sumPt->Fill(sumPt);
0535 
0536     int half_trks = twoPV.nTracks();
0537 
0538     const double invSqrt2 = 1. / std::sqrt(2.);
0539 
0540     double deltaX = (twoPV.x() - onePV.x());
0541     double deltaY = (twoPV.y() - onePV.y());
0542     double deltaZ = (twoPV.z() - onePV.z());
0543 
0544     double resX = deltaX * invSqrt2;
0545     double resY = deltaY * invSqrt2;
0546     double resZ = deltaZ * invSqrt2;
0547 
0548     h_diffX->Fill(resX * cmToUm);
0549     h_diffY->Fill(resY * cmToUm);
0550     h_diffZ->Fill(resZ * cmToUm);
0551 
0552     double errX = sqrt(pow(twoPV.xError(), 2) + pow(onePV.xError(), 2));
0553     double errY = sqrt(pow(twoPV.yError(), 2) + pow(onePV.yError(), 2));
0554     double errZ = sqrt(pow(twoPV.zError(), 2) + pow(onePV.zError(), 2));
0555 
0556     h_errX->Fill(errX * cmToUm);
0557     h_errY->Fill(errY * cmToUm);
0558     h_errZ->Fill(errZ * cmToUm);
0559 
0560     h_pullX->Fill(deltaX / errX);
0561     h_pullY->Fill(deltaY / errY);
0562     h_pullZ->Fill(deltaZ / errZ);
0563 
0564     // filling the pT-binned distributions
0565 
0566     for (int ipTBin = 0; ipTBin < nPtBins_; ipTBin++) {
0567       float pTF = mypT_bins_[ipTBin];
0568       float pTL = mypT_bins_[ipTBin + 1];
0569 
0570       if (avgSumPt >= pTF && avgSumPt < pTL) {
0571         PVValHelper::fillByIndex(h_resolX_sumPt_, ipTBin, resX * cmToUm, "1");
0572         PVValHelper::fillByIndex(h_resolY_sumPt_, ipTBin, resY * cmToUm, "2");
0573         PVValHelper::fillByIndex(h_resolZ_sumPt_, ipTBin, resZ * cmToUm, "3");
0574 
0575         PVValHelper::fillByIndex(h_pullX_sumPt_, ipTBin, deltaX / errX, "4");
0576         PVValHelper::fillByIndex(h_pullY_sumPt_, ipTBin, deltaY / errY, "5");
0577         PVValHelper::fillByIndex(h_pullZ_sumPt_, ipTBin, deltaZ / errZ, "6");
0578       }
0579     }
0580 
0581     // filling the track multeplicity binned distributions
0582 
0583     for (int inTrackBin = 0; inTrackBin < nTrackBins_; inTrackBin++) {
0584       float nTrackF = myNTrack_bins_[inTrackBin];
0585       float nTrackL = myNTrack_bins_[inTrackBin + 1];
0586       if (ntrks >= nTrackF && ntrks < nTrackL) {
0587         //if (ntrks == inTrackBin) {
0588         PVValHelper::fillByIndex(h_resolX_Ntracks_, inTrackBin, resX * cmToUm, "7");
0589         PVValHelper::fillByIndex(h_resolY_Ntracks_, inTrackBin, resY * cmToUm, "8");
0590         PVValHelper::fillByIndex(h_resolZ_Ntracks_, inTrackBin, resZ * cmToUm, "9");
0591 
0592         PVValHelper::fillByIndex(h_pullX_Ntracks_, inTrackBin, deltaX / errX, "10");
0593         PVValHelper::fillByIndex(h_pullY_Ntracks_, inTrackBin, deltaY / errY, "11");
0594         PVValHelper::fillByIndex(h_pullZ_Ntracks_, inTrackBin, deltaZ / errZ, "12");
0595       }
0596     }
0597 
0598     // filling the vertex multeplicity binned distributions
0599 
0600     for (int inVtxBin = 0; inVtxBin < nVtxBins_; inVtxBin++) {
0601       float nVtxF = myNVtx_bins_[inVtxBin];
0602       float nVtxL = myNVtx_bins_[inVtxBin + 1];
0603       if (nOfflineVtx >= nVtxF && nOfflineVtx < nVtxL) {
0604         //if (nOfflineVtx == inVtxBin) {
0605         PVValHelper::fillByIndex(h_resolX_Nvtx_, inVtxBin, deltaX * cmToUm, "7");
0606         PVValHelper::fillByIndex(h_resolY_Nvtx_, inVtxBin, deltaY * cmToUm, "8");
0607         PVValHelper::fillByIndex(h_resolZ_Nvtx_, inVtxBin, deltaZ * cmToUm, "9");
0608 
0609         PVValHelper::fillByIndex(h_pullX_Nvtx_, inVtxBin, deltaX / errX, "10");
0610         PVValHelper::fillByIndex(h_pullY_Nvtx_, inVtxBin, deltaY / errY, "11");
0611         PVValHelper::fillByIndex(h_pullZ_Nvtx_, inVtxBin, deltaZ / errZ, "12");
0612       }
0613     }
0614 
0615     h_sumPt1->Fill(sumPt1);
0616     h_sumPt2->Fill(sumPt2);
0617 
0618     h_minWTrks1->Fill(theminW1);
0619     h_minWTrks2->Fill(theminW2);
0620 
0621     h_PVCL_subVtx1->Fill(TMath::Prob(pvOne.totalChiSquared(), (int)(pvOne.degreesOfFreedom())));
0622     h_PVCL_subVtx2->Fill(TMath::Prob(pvTwo.totalChiSquared(), (int)(pvTwo.degreesOfFreedom())));
0623 
0624     // fill ntuples
0625     pvCand thePV;
0626     thePV.ipos = counter;
0627     thePV.nTrks = ntrks;
0628 
0629     thePV.x_origVtx = iPV.x();
0630     thePV.y_origVtx = iPV.y();
0631     thePV.z_origVtx = iPV.z();
0632 
0633     thePV.xErr_origVtx = iPV.xError();
0634     thePV.yErr_origVtx = iPV.yError();
0635     thePV.zErr_origVtx = iPV.zError();
0636 
0637     thePV.n_subVtx1 = half_trks;
0638     thePV.x_subVtx1 = onePV.x();
0639     thePV.y_subVtx1 = onePV.y();
0640     thePV.z_subVtx1 = onePV.z();
0641 
0642     thePV.xErr_subVtx1 = onePV.xError();
0643     thePV.yErr_subVtx1 = onePV.yError();
0644     thePV.zErr_subVtx1 = onePV.zError();
0645     thePV.sumPt_subVtx1 = sumPt1;
0646 
0647     thePV.n_subVtx2 = half_trks;
0648     thePV.x_subVtx2 = twoPV.x();
0649     thePV.y_subVtx2 = twoPV.y();
0650     thePV.z_subVtx2 = twoPV.z();
0651 
0652     thePV.xErr_subVtx2 = twoPV.xError();
0653     thePV.yErr_subVtx2 = twoPV.yError();
0654     thePV.zErr_subVtx2 = twoPV.zError();
0655     thePV.sumPt_subVtx2 = sumPt2;
0656 
0657     thePV.CL_subVtx1 = TMath::Prob(pvOne.totalChiSquared(), (int)(pvOne.degreesOfFreedom()));
0658     thePV.CL_subVtx2 = TMath::Prob(pvTwo.totalChiSquared(), (int)(pvTwo.degreesOfFreedom()));
0659 
0660     thePV.minW_subVtx1 = theminW1;
0661     thePV.minW_subVtx2 = theminW2;
0662 
0663     event_.pvs.push_back(thePV);
0664 
0665   }  // loop on the vertices
0666 
0667   // fill the histogram of vertices per event
0668   h_nVertices->Fill(counter);
0669   h_nNonFakeVertices->Fill(noFakecounter);
0670   h_nFinalVertices->Fill(goodcounter);
0671 
0672   if (storeNtuple_) {
0673     tree_->Fill();
0674   }
0675 }
0676 
0677 void SplitVertexResolution::beginEvent() {
0678   event_.pvs.clear();
0679   event_.nVtx = -1;
0680 }
0681 
0682 void SplitVertexResolution::beginRun(edm::Run const& run, edm::EventSetup const& iSetup) {
0683   unsigned int RunNumber_ = run.run();
0684 
0685   if (!runNumbersTimesLog_.count(RunNumber_)) {
0686     auto times = getRunTime(iSetup);
0687 
0688     if (debug_) {
0689       const time_t start_time = times.first / 1.0e+6;
0690       edm::LogInfo("SplitVertexResolution")
0691           << RunNumber_ << " has start time: " << times.first << " - " << times.second << std::endl;
0692       edm::LogInfo("SplitVertexResolution")
0693           << "human readable time: " << std::asctime(std::gmtime(&start_time)) << std::endl;
0694     }
0695     runNumbersTimesLog_[RunNumber_] = times;
0696   }
0697 }
0698 
0699 // ------------ method called once each job just before starting event loop  ------------
0700 void SplitVertexResolution::beginJob() {
0701   ievt = 0;
0702   itrks = 0;
0703 
0704   if (compressionSettings_ > 0) {
0705     outfile_->file().SetCompressionSettings(compressionSettings_);
0706   }
0707 
0708   // binning histogram
0709   TFileDirectory BinningFeatures = outfile_->mkdir("BinningFeatures");
0710   // sumpT
0711   h_profileBinnings = BinningFeatures.make<TH1F>("h_profileBinnings", "profile binnings", 4, -0.5, 3.5);
0712   h_profileBinnings->GetXaxis()->SetBinLabel(1, "#sum p_{T} min");
0713   h_profileBinnings->GetXaxis()->SetBinLabel(2, "#sum p_{T} max");
0714   // n. tracks
0715   h_profileBinnings->GetXaxis()->SetBinLabel(3, "n. track bins");
0716   // n. vertices
0717   h_profileBinnings->GetXaxis()->SetBinLabel(4, "n. vertices bins");
0718 
0719   // fill the binning book-keeping
0720   h_profileBinnings->SetBinContent(1, sumpTStartScale_);
0721   h_profileBinnings->SetBinContent(2, sumpTEndScale_);
0722   h_profileBinnings->SetBinContent(3, nVisTrackBins_);
0723   h_profileBinnings->SetBinContent(4, nVisVtxBins_);
0724 
0725   // luminosity histo
0726   TFileDirectory EventFeatures = outfile_->mkdir("EventFeatures");
0727   h_lumiFromConfig =
0728       EventFeatures.make<TH1F>("h_lumiFromConfig", "luminosity from config;;luminosity of present run", 1, -0.5, 0.5);
0729   h_lumiFromConfig->SetBinContent(1, intLumi_);
0730 
0731   h_runFromConfig = EventFeatures.make<TH1I>("h_runFromConfig",
0732                                              "run number from config;;run number (from configuration)",
0733                                              runControlNumbers_.size(),
0734                                              0.,
0735                                              runControlNumbers_.size());
0736 
0737   for (const auto& run : runControlNumbers_ | boost::adaptors::indexed(1)) {
0738     h_runFromConfig->SetBinContent(run.index(), run.value());
0739   }
0740 
0741   // resolutions
0742 
0743   if (!checkBinOrdering(mypT_bins_)) {
0744     edm::LogError("SplitVertexResolution") << " Warning - the vector of pT bins is not ordered " << std::endl;
0745   }
0746 
0747   if (!checkBinOrdering(myNTrack_bins_)) {
0748     edm::LogError("SplitVertexResolution") << " Warning -the vector of n. tracks bins is not ordered " << std::endl;
0749   }
0750 
0751   if (!checkBinOrdering(myNVtx_bins_)) {
0752     edm::LogError("SplitVertexResolution") << " Warning -the vector of n. vertices bins is not ordered " << std::endl;
0753   }
0754 
0755   TFileDirectory xResolSumPt = outfile_->mkdir("xResolSumPt");
0756   h_resolX_sumPt_ = bookResidualsHistogram(xResolSumPt, nPtBins_, "resolX", "sumPt");
0757 
0758   TFileDirectory yResolSumPt = outfile_->mkdir("yResolSumPt");
0759   h_resolY_sumPt_ = bookResidualsHistogram(yResolSumPt, nPtBins_, "resolY", "sumPt");
0760 
0761   TFileDirectory zResolSumPt = outfile_->mkdir("zResolSumPt");
0762   h_resolZ_sumPt_ = bookResidualsHistogram(zResolSumPt, nPtBins_, "resolZ", "sumPt");
0763 
0764   TFileDirectory xResolNtracks_ = outfile_->mkdir("xResolNtracks");
0765   h_resolX_Ntracks_ = bookResidualsHistogram(xResolNtracks_, nTrackBins_, "resolX", "Ntracks");
0766 
0767   TFileDirectory yResolNtracks_ = outfile_->mkdir("yResolNtracks");
0768   h_resolY_Ntracks_ = bookResidualsHistogram(yResolNtracks_, nTrackBins_, "resolY", "Ntracks");
0769 
0770   TFileDirectory zResolNtracks_ = outfile_->mkdir("zResolNtracks");
0771   h_resolZ_Ntracks_ = bookResidualsHistogram(zResolNtracks_, nTrackBins_, "resolZ", "Ntracks");
0772 
0773   TFileDirectory xResolNvtx_ = outfile_->mkdir("xResolNvtx");
0774   h_resolX_Nvtx_ = bookResidualsHistogram(xResolNvtx_, nVtxBins_, "resolX", "Nvtx");
0775 
0776   TFileDirectory yResolNvtx_ = outfile_->mkdir("yResolNvtx");
0777   h_resolY_Nvtx_ = bookResidualsHistogram(yResolNvtx_, nVtxBins_, "resolY", "Nvtx");
0778 
0779   TFileDirectory zResolNvtx_ = outfile_->mkdir("zResolNvtx");
0780   h_resolZ_Nvtx_ = bookResidualsHistogram(zResolNvtx_, nVtxBins_, "resolZ", "Nvtx");
0781 
0782   // pulls
0783 
0784   TFileDirectory xPullSumPt = outfile_->mkdir("xPullSumPt");
0785   h_pullX_sumPt_ = bookResidualsHistogram(xPullSumPt, nPtBins_, "pullX", "sumPt");
0786 
0787   TFileDirectory yPullSumPt = outfile_->mkdir("yPullSumPt");
0788   h_pullY_sumPt_ = bookResidualsHistogram(yPullSumPt, nPtBins_, "pullY", "sumPt");
0789 
0790   TFileDirectory zPullSumPt = outfile_->mkdir("zPullSumPt");
0791   h_pullZ_sumPt_ = bookResidualsHistogram(zPullSumPt, nPtBins_, "pullZ", "sumPt");
0792 
0793   TFileDirectory xPullNtracks_ = outfile_->mkdir("xPullNtracks");
0794   h_pullX_Ntracks_ = bookResidualsHistogram(xPullNtracks_, nTrackBins_, "pullX", "Ntracks");
0795 
0796   TFileDirectory yPullNtracks_ = outfile_->mkdir("yPullNtracks");
0797   h_pullY_Ntracks_ = bookResidualsHistogram(yPullNtracks_, nTrackBins_, "pullY", "Ntracks");
0798 
0799   TFileDirectory zPullNtracks_ = outfile_->mkdir("zPullNtracks");
0800   h_pullZ_Ntracks_ = bookResidualsHistogram(zPullNtracks_, nTrackBins_, "pullZ", "Ntracks");
0801 
0802   TFileDirectory xPullNvtx_ = outfile_->mkdir("xPullNvtx");
0803   h_pullX_Nvtx_ = bookResidualsHistogram(xPullNvtx_, nVtxBins_, "pullX", "Nvtx");
0804 
0805   TFileDirectory yPullNvtx_ = outfile_->mkdir("yPullNvtx");
0806   h_pullY_Nvtx_ = bookResidualsHistogram(yPullNvtx_, nVtxBins_, "pullY", "Nvtx");
0807 
0808   TFileDirectory zPullNvtx_ = outfile_->mkdir("zPullNvtx");
0809   h_pullZ_Nvtx_ = bookResidualsHistogram(zPullNvtx_, nVtxBins_, "pullZ", "Nvtx");
0810 
0811   // control plots
0812   h_runNumber = outfile_->make<TH1F>("h_runNumber", "run number;run number;n_{events}", 100000, 250000., 350000.);
0813 
0814   h_nOfflineVertices = outfile_->make<TH1I>("h_nOfflineVertices", "n. of vertices;n. vertices; events", 100, 0, 100);
0815   h_nVertices = outfile_->make<TH1I>("h_nVertices", "n. of vertices;n. vertices; events", 100, 0, 100);
0816   h_nNonFakeVertices =
0817       outfile_->make<TH1I>("h_nRealVertices", "n. of non-fake vertices;n. vertices; events", 100, 0, 100);
0818   h_nFinalVertices =
0819       outfile_->make<TH1I>("h_nSelectedVertices", "n. of selected vertices vertices;n. vertices; events", 100, 0, 100);
0820 
0821   h_diffX = outfile_->make<TH1F>(
0822       "h_diffX", "x-coordinate vertex resolution;vertex resolution (x) [#mum];vertices", 100, -300, 300.);
0823   h_diffY = outfile_->make<TH1F>(
0824       "h_diffY", "y-coordinate vertex resolution;vertex resolution (y) [#mum];vertices", 100, -300, 300.);
0825   h_diffZ = outfile_->make<TH1F>(
0826       "h_diffZ", "z-coordinate vertex resolution;vertex resolution (z) [#mum];vertices", 100, -500, 500.);
0827 
0828   h_OrigVertexErrX = outfile_->make<TH1F>(
0829       "h_OrigVertexErrX", "x-coordinate vertex error;vertex error (x) [#mum];vertices", 300, 0., 300.);
0830   h_OrigVertexErrY = outfile_->make<TH1F>(
0831       "h_OrigVertexErrY", "y-coordinate vertex error;vertex error (y) [#mum];vertices", 300, 0., 300.);
0832   h_OrigVertexErrZ = outfile_->make<TH1F>(
0833       "h_OrigVertexErrZ", "z-coordinate vertex error;vertex error (z) [#mum];vertices", 500, 0., 500.);
0834 
0835   h_errX = outfile_->make<TH1F>(
0836       "h_errX", "x-coordinate vertex resolution error;vertex resoltuion error (x) [#mum];vertices", 300, 0., 300.);
0837   h_errY = outfile_->make<TH1F>(
0838       "h_errY", "y-coordinate vertex resolution error;vertex resolution error (y) [#mum];vertices", 300, 0., 300.);
0839   h_errZ = outfile_->make<TH1F>(
0840       "h_errZ", "z-coordinate vertex resolution error;vertex resolution error (z) [#mum];vertices", 500, 0., 500.);
0841 
0842   h_pullX = outfile_->make<TH1F>("h_pullX", "x-coordinate vertex pull;vertex pull (x);vertices", 500, -10, 10.);
0843   h_pullY = outfile_->make<TH1F>("h_pullY", "y-coordinate vertex pull;vertex pull (y);vertices", 500, -10, 10.);
0844   h_pullZ = outfile_->make<TH1F>("h_pullZ", "z-coordinate vertex pull;vertex pull (z);vertices", 500, -10, 10.);
0845 
0846   h_ntrks = outfile_->make<TH1F>("h_ntrks",
0847                                  "number of tracks in vertex;vertex multeplicity;vertices",
0848                                  myNTrack_bins_.size() - 1,
0849                                  myNTrack_bins_.data());
0850 
0851   h_sumPt = outfile_->make<TH1F>(
0852       "h_sumPt", "#Sigma p_{T};#sum p_{T} [GeV];vertices", mypT_bins_.size() - 1, mypT_bins_.data());
0853 
0854   h_avgSumPt = outfile_->make<TH1F>(
0855       "h_avgSumPt", "#LT #Sigma p_{T} #GT;#LT #sum p_{T} #GT [GeV];vertices", mypT_bins_.size() - 1, mypT_bins_.data());
0856 
0857   h_sumPt1 = outfile_->make<TH1F>("h_sumPt1",
0858                                   "#Sigma p_{T} sub-vertex 1;#sum p_{T} sub-vertex 1 [GeV];subvertices",
0859                                   mypT_bins_.size() - 1,
0860                                   mypT_bins_.data());
0861   h_sumPt2 = outfile_->make<TH1F>("h_sumPt2",
0862                                   "#Sigma p_{T} sub-vertex 2;#sum p_{T} sub-vertex 2 [GeV];subvertices",
0863                                   mypT_bins_.size() - 1,
0864                                   mypT_bins_.data());
0865 
0866   h_wTrks1 = outfile_->make<TH1F>("h_wTrks1", "weight of track for sub-vertex 1;track weight;subvertices", 500, 0., 1.);
0867   h_wTrks2 = outfile_->make<TH1F>("h_wTrks2", "weithg of track for sub-vertex 2;track weight;subvertices", 500, 0., 1.);
0868 
0869   h_minWTrks1 = outfile_->make<TH1F>(
0870       "h_minWTrks1", "minimum weight of track for sub-vertex 1;minimum track weight;subvertices", 500, 0., 1.);
0871   h_minWTrks2 = outfile_->make<TH1F>(
0872       "h_minWTrks2", "minimum weithg of track for sub-vertex 2;minimum track weight;subvertices", 500, 0., 1.);
0873 
0874   h_PVCL_subVtx1 =
0875       outfile_->make<TH1F>("h_PVCL_subVtx1",
0876                            "#chi^{2} probability for sub-vertex 1;Prob(#chi^{2},ndof) for sub-vertex 1;subvertices",
0877                            100,
0878                            0.,
0879                            1);
0880   h_PVCL_subVtx2 =
0881       outfile_->make<TH1F>("h_PVCL_subVtx2",
0882                            "#chi^{2} probability for sub-vertex 2;Prob(#chi^{2},ndof) for sub-vertex 2;subvertices",
0883                            100,
0884                            0.,
0885                            1);
0886 
0887   // resolutions
0888 
0889   p_resolX_vsSumPt = outfile_->make<TH1F>("p_resolX_vsSumPt",
0890                                           "x-resolution vs #Sigma p_{T};#sum p_{T} [GeV]; x vertex resolution [#mum]",
0891                                           mypT_bins_.size() - 1,
0892                                           mypT_bins_.data());
0893   p_resolY_vsSumPt = outfile_->make<TH1F>("p_resolY_vsSumPt",
0894                                           "y-resolution vs #Sigma p_{T};#sum p_{T} [GeV]; y vertex resolution [#mum]",
0895                                           mypT_bins_.size() - 1,
0896                                           mypT_bins_.data());
0897   p_resolZ_vsSumPt = outfile_->make<TH1F>("p_resolZ_vsSumPt",
0898                                           "z-resolution vs #Sigma p_{T};#sum p_{T} [GeV]; z vertex resolution [#mum]",
0899                                           mypT_bins_.size() - 1,
0900                                           mypT_bins_.data());
0901 
0902   p_resolX_vsNtracks = outfile_->make<TH1F>("p_resolX_vsNtracks",
0903                                             "x-resolution vs n_{tracks};n_{tracks}; x vertex resolution [#mum]",
0904                                             myNTrack_bins_.size() - 1,
0905                                             myNTrack_bins_.data());
0906   p_resolY_vsNtracks = outfile_->make<TH1F>("p_resolY_vsNtracks",
0907                                             "y-resolution vs n_{tracks};n_{tracks}; y vertex resolution [#mum]",
0908                                             myNTrack_bins_.size() - 1,
0909                                             myNTrack_bins_.data());
0910   p_resolZ_vsNtracks = outfile_->make<TH1F>("p_resolZ_vsNtracks",
0911                                             "z-resolution vs n_{tracks};n_{tracks}; z vertex resolution [#mum]",
0912                                             myNTrack_bins_.size() - 1,
0913                                             myNTrack_bins_.data());
0914 
0915   p_resolX_vsNvtx = outfile_->make<TH1F>("p_resolX_vsNvtx",
0916                                          "x-resolution vs n_{vertices};n_{vertices}; x vertex resolution [#mum]",
0917                                          myNVtx_bins_.size() - 1,
0918                                          myNVtx_bins_.data());
0919   p_resolY_vsNvtx = outfile_->make<TH1F>("p_resolY_vsNvtx",
0920                                          "y-resolution vs n_{vertices};n_{vertices}; y vertex resolution [#mum]",
0921                                          myNVtx_bins_.size() - 1,
0922                                          myNVtx_bins_.data());
0923   p_resolZ_vsNvtx = outfile_->make<TH1F>("p_resolZ_vsNvtx",
0924                                          "z-resolution vs n_{vertices};n_{vertices}; z vertex resolution [#mum]",
0925                                          myNVtx_bins_.size() - 1,
0926                                          myNVtx_bins_.data());
0927 
0928   // pulls
0929 
0930   p_pullX_vsSumPt = outfile_->make<TH1F>("p_pullX_vsSumPt",
0931                                          "x-pull vs #Sigma p_{T};#sum p_{T} [GeV]; x vertex pull",
0932                                          mypT_bins_.size() - 1,
0933                                          mypT_bins_.data());
0934   p_pullY_vsSumPt = outfile_->make<TH1F>("p_pullY_vsSumPt",
0935                                          "y-pull vs #Sigma p_{T};#sum p_{T} [GeV]; y vertex pull",
0936                                          mypT_bins_.size() - 1,
0937                                          mypT_bins_.data());
0938   p_pullZ_vsSumPt = outfile_->make<TH1F>("p_pullZ_vsSumPt",
0939                                          "z-pull vs #Sigma p_{T};#sum p_{T} [GeV]; z vertex pull",
0940                                          mypT_bins_.size() - 1,
0941                                          mypT_bins_.data());
0942 
0943   p_pullX_vsNtracks = outfile_->make<TH1F>("p_pullX_vsNtracks",
0944                                            "x-pull vs n_{tracks};n_{tracks}; x vertex pull",
0945                                            myNTrack_bins_.size() - 1,
0946                                            myNTrack_bins_.data());
0947   p_pullY_vsNtracks = outfile_->make<TH1F>("p_pullY_vsNtracks",
0948                                            "y-pull vs n_{tracks};n_{tracks}; y vertex pull",
0949                                            myNTrack_bins_.size() - 1,
0950                                            myNTrack_bins_.data());
0951   p_pullZ_vsNtracks = outfile_->make<TH1F>("p_pullZ_vsNtracks",
0952                                            "z-pull vs n_{tracks};n_{tracks}; z vertex pull",
0953                                            myNTrack_bins_.size() - 1,
0954                                            myNTrack_bins_.data());
0955 
0956   p_pullX_vsNvtx = outfile_->make<TH1F>("p_pullX_vsNvtx",
0957                                         "x-pull vs n_{vertices};n_{vertices}; x vertex pull",
0958                                         myNVtx_bins_.size() - 1,
0959                                         myNVtx_bins_.data());
0960   p_pullY_vsNvtx = outfile_->make<TH1F>("p_pullY_vsNvtx",
0961                                         "y-pull vs n_{vertices};n_{vertices}; y vertex pull",
0962                                         myNVtx_bins_.size() - 1,
0963                                         myNVtx_bins_.data());
0964   p_pullZ_vsNvtx = outfile_->make<TH1F>("p_pullZ_vsNvtx",
0965                                         "z-pull vs n_{vertices};n_{vertices}; z vertex pull",
0966                                         myNVtx_bins_.size() - 1,
0967                                         myNVtx_bins_.data());
0968 
0969   tree_ = outfile_->make<TTree>("pvTree", "pvTree");
0970   tree_->Branch("event", &event_, 64000, 2);
0971 }
0972 
0973 //*************************************************************
0974 // Generic booker function
0975 //*************************************************************
0976 std::vector<TH1F*> SplitVertexResolution::bookResidualsHistogram(TFileDirectory dir,
0977                                                                  unsigned int theNOfBins,
0978                                                                  TString resType,
0979                                                                  TString varType) {
0980   TH1F::SetDefaultSumw2(kTRUE);
0981 
0982   double up = 500.;
0983   double down = -500.;
0984 
0985   if (resType.Contains("pull")) {
0986     up *= 0.01;
0987     down *= 0.01;
0988   }
0989 
0990   std::vector<TH1F*> h;
0991   h.reserve(theNOfBins);
0992 
0993   const char* auxResType = (resType.ReplaceAll("_", "")).Data();
0994 
0995   for (unsigned int i = 0; i < theNOfBins; i++) {
0996     TH1F* htemp = dir.make<TH1F>(Form("histo_%s_%s_plot%i", resType.Data(), varType.Data(), i),
0997                                  Form("%s vs %s - bin %i;%s;vertices", auxResType, varType.Data(), i, auxResType),
0998                                  250,
0999                                  down,
1000                                  up);
1001     h.push_back(htemp);
1002   }
1003 
1004   return h;
1005 }
1006 
1007 // ------------ method called once each job just after ending the event loop  ------------
1008 void SplitVertexResolution::endJob() {
1009   edm::LogVerbatim("SplitVertexResolution") << "*******************************" << std::endl;
1010   edm::LogVerbatim("SplitVertexResolution") << "Events run in total: " << ievt << std::endl;
1011   edm::LogVerbatim("SplitVertexResolution") << "n. tracks: " << itrks << std::endl;
1012   edm::LogVerbatim("SplitVertexResolution") << "*******************************" << std::endl;
1013 
1014   int nFiringTriggers = triggerMap_.size();
1015   edm::LogVerbatim("SplitVertexResolution") << "firing triggers: " << nFiringTriggers << std::endl;
1016   edm::LogVerbatim("SplitVertexResolution") << "*******************************" << std::endl;
1017 
1018   tksByTrigger_ = outfile_->make<TH1D>(
1019       "tksByTrigger", "tracks by HLT path;;% of # traks", nFiringTriggers, -0.5, nFiringTriggers - 0.5);
1020   evtsByTrigger_ = outfile_->make<TH1D>(
1021       "evtsByTrigger", "events by HLT path;;% of # events", nFiringTriggers, -0.5, nFiringTriggers - 0.5);
1022 
1023   int i = 0;
1024   for (std::map<std::string, std::pair<int, int>>::iterator it = triggerMap_.begin(); it != triggerMap_.end(); ++it) {
1025     i++;
1026 
1027     double trkpercent = ((it->second).second) * 100. / double(itrks);
1028     double evtpercent = ((it->second).first) * 100. / double(ievt);
1029 
1030     edm::LogVerbatim("SplitVertexResolution")
1031         << "HLT path: " << std::setw(60) << std::left << it->first << " | events firing: " << std::right << std::setw(8)
1032         << (it->second).first << " (" << std::setw(8) << std::fixed << std::setprecision(4) << evtpercent << "%)"
1033         << " | tracks collected: " << std::setw(8) << (it->second).second << " (" << std::setw(8) << std::fixed
1034         << std::setprecision(4) << trkpercent << "%)";
1035 
1036     tksByTrigger_->SetBinContent(i, trkpercent);
1037     tksByTrigger_->GetXaxis()->SetBinLabel(i, (it->first).c_str());
1038 
1039     evtsByTrigger_->SetBinContent(i, evtpercent);
1040     evtsByTrigger_->GetXaxis()->SetBinLabel(i, (it->first).c_str());
1041   }
1042 
1043   TFileDirectory RunFeatures = outfile_->mkdir("RunFeatures");
1044   h_runStartTimes = RunFeatures.make<TH1I>(
1045       "runStartTimes", "run start times", runNumbersTimesLog_.size(), 0, runNumbersTimesLog_.size());
1046   h_runEndTimes =
1047       RunFeatures.make<TH1I>("runEndTimes", "run end times", runNumbersTimesLog_.size(), 0, runNumbersTimesLog_.size());
1048 
1049   unsigned int count = 1;
1050   for (const auto& run : runNumbersTimesLog_) {
1051     // strip down the microseconds
1052     h_runStartTimes->SetBinContent(count, run.second.first / 10e6);
1053     h_runStartTimes->GetXaxis()->SetBinLabel(count, (std::to_string(run.first)).c_str());
1054 
1055     h_runEndTimes->SetBinContent(count, run.second.second / 10e6);
1056     h_runEndTimes->GetXaxis()->SetBinLabel(count, (std::to_string(run.first)).c_str());
1057 
1058     count++;
1059   }
1060 
1061   // resolutions
1062 
1063   fillTrendPlotByIndex(p_resolX_vsSumPt, h_resolX_sumPt_, PVValHelper::WIDTH);
1064   fillTrendPlotByIndex(p_resolY_vsSumPt, h_resolY_sumPt_, PVValHelper::WIDTH);
1065   fillTrendPlotByIndex(p_resolZ_vsSumPt, h_resolZ_sumPt_, PVValHelper::WIDTH);
1066 
1067   fillTrendPlotByIndex(p_resolX_vsNtracks, h_resolX_Ntracks_, PVValHelper::WIDTH);
1068   fillTrendPlotByIndex(p_resolY_vsNtracks, h_resolY_Ntracks_, PVValHelper::WIDTH);
1069   fillTrendPlotByIndex(p_resolZ_vsNtracks, h_resolZ_Ntracks_, PVValHelper::WIDTH);
1070 
1071   fillTrendPlotByIndex(p_resolX_vsNvtx, h_resolX_Nvtx_, PVValHelper::WIDTH);
1072   fillTrendPlotByIndex(p_resolY_vsNvtx, h_resolY_Nvtx_, PVValHelper::WIDTH);
1073   fillTrendPlotByIndex(p_resolZ_vsNvtx, h_resolZ_Nvtx_, PVValHelper::WIDTH);
1074 
1075   // pulls
1076 
1077   fillTrendPlotByIndex(p_pullX_vsSumPt, h_pullX_sumPt_, PVValHelper::WIDTH);
1078   fillTrendPlotByIndex(p_pullY_vsSumPt, h_pullY_sumPt_, PVValHelper::WIDTH);
1079   fillTrendPlotByIndex(p_pullZ_vsSumPt, h_pullZ_sumPt_, PVValHelper::WIDTH);
1080 
1081   fillTrendPlotByIndex(p_pullX_vsNtracks, h_pullX_Ntracks_, PVValHelper::WIDTH);
1082   fillTrendPlotByIndex(p_pullY_vsNtracks, h_pullY_Ntracks_, PVValHelper::WIDTH);
1083   fillTrendPlotByIndex(p_pullZ_vsNtracks, h_pullZ_Ntracks_, PVValHelper::WIDTH);
1084 
1085   fillTrendPlotByIndex(p_pullX_vsNvtx, h_pullX_Nvtx_, PVValHelper::WIDTH);
1086   fillTrendPlotByIndex(p_pullY_vsNvtx, h_pullY_Nvtx_, PVValHelper::WIDTH);
1087   fillTrendPlotByIndex(p_pullZ_vsNvtx, h_pullZ_Nvtx_, PVValHelper::WIDTH);
1088 }
1089 
1090 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
1091 void SplitVertexResolution::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
1092   edm::ParameterSetDescription desc;
1093   desc.setComment(
1094       "Validates alignment payloads by evaluating the resulting Primary Vertex Resolution via vertex splitting method");
1095 
1096   desc.addUntracked<int>("compressionSettings", -1);
1097   desc.add<bool>("storeNtuple", false);
1098   desc.addUntracked<double>("intLumi", 0.);
1099   desc.addUntracked<bool>("Debug", false);
1100   desc.add<edm::InputTag>("vtxCollection", edm::InputTag("offlinePrimaryVertices"));
1101   desc.add<edm::InputTag>("trackCollection", edm::InputTag("generalTracks"));
1102   desc.addUntracked<double>("minVertexNdf", 10);
1103   desc.addUntracked<double>("minVertexMeanWeight", 0.5);
1104   desc.addUntracked<bool>("runControl", false);
1105   desc.addUntracked<std::vector<unsigned int>>("runControlNumber", {});
1106   desc.addUntracked<double>("sumpTStartScale", 1.);
1107   desc.addUntracked<double>("sumpTEndScale", 1e3);
1108   desc.addUntracked<double>("nTrackBins", 120.);
1109   desc.addUntracked<double>("nVtxBins", 60.);
1110   descriptions.addWithDefaultLabel(desc);
1111 }
1112 
1113 //*************************************************************
1114 std::pair<long long, long long> SplitVertexResolution::getRunTime(const edm::EventSetup& iSetup) const
1115 //*************************************************************
1116 {
1117   const auto& runInfo = iSetup.getData(runInfoToken_);
1118   if (debug_) {
1119     edm::LogInfo("SplitVertexResolution")
1120         << "start time: " << runInfo.m_start_time_str << " - stop time: " << runInfo.m_stop_time_str << std::endl;
1121   }
1122   return std::make_pair(runInfo.m_start_time_ll, runInfo.m_stop_time_ll);
1123 }
1124 
1125 //*************************************************************
1126 void SplitVertexResolution::fillTrendPlotByIndex(TH1F* trendPlot, std::vector<TH1F*>& h, PVValHelper::estimator fitPar_)
1127 //*************************************************************
1128 {
1129   for (auto iterator = h.begin(); iterator != h.end(); iterator++) {
1130     unsigned int bin = std::distance(h.begin(), iterator) + 1;
1131     statmode::fitParams myFit = fitResiduals((*iterator));
1132 
1133     switch (fitPar_) {
1134       case PVValHelper::MEAN: {
1135         float mean_ = myFit.first.value();
1136         float meanErr_ = myFit.first.error();
1137         trendPlot->SetBinContent(bin, mean_);
1138         trendPlot->SetBinError(bin, meanErr_);
1139         break;
1140       }
1141       case PVValHelper::WIDTH: {
1142         float width_ = myFit.second.value();
1143         float widthErr_ = myFit.second.error();
1144         trendPlot->SetBinContent(bin, width_);
1145         trendPlot->SetBinError(bin, widthErr_);
1146         break;
1147       }
1148       case PVValHelper::MEDIAN: {
1149         float median_ = PVValHelper::getMedian((*iterator)).value();
1150         float medianErr_ = PVValHelper::getMedian((*iterator)).error();
1151         trendPlot->SetBinContent(bin, median_);
1152         trendPlot->SetBinError(bin, medianErr_);
1153         break;
1154       }
1155       case PVValHelper::MAD: {
1156         float mad_ = PVValHelper::getMAD((*iterator)).value();
1157         float madErr_ = PVValHelper::getMAD((*iterator)).error();
1158         trendPlot->SetBinContent(bin, mad_);
1159         trendPlot->SetBinError(bin, madErr_);
1160         break;
1161       }
1162       default:
1163         edm::LogWarning("SplitVertexResolution")
1164             << "fillTrendPlotByIndex() " << fitPar_ << " unknown estimator!" << std::endl;
1165         break;
1166     }
1167   }
1168 }
1169 
1170 //*************************************************************
1171 statmode::fitParams SplitVertexResolution::fitResiduals(TH1* hist, bool singleTime)
1172 //*************************************************************
1173 {
1174   if (hist->GetEntries() < 10) {
1175     LogDebug("SplitVertexResolution") << "hist name: " << hist->GetName() << " has less than 10 entries" << std::endl;
1176     return std::make_pair(Measurement1D(0., 0.), Measurement1D(0., 0.));
1177   }
1178 
1179   float maxHist = hist->GetXaxis()->GetXmax();
1180   float minHist = hist->GetXaxis()->GetXmin();
1181   float mean = hist->GetMean();
1182   float sigma = hist->GetRMS();
1183 
1184   if (edm::isNotFinite(mean) || edm::isNotFinite(sigma)) {
1185     mean = 0;
1186     //sigma= - hist->GetXaxis()->GetBinLowEdge(1) + hist->GetXaxis()->GetBinLowEdge(hist->GetNbinsX()+1);
1187     sigma = -minHist + maxHist;
1188     edm::LogWarning("SplitVertexResolution")
1189         << "FitPVResiduals::fitResiduals(): histogram" << hist->GetName() << " mean or sigma are NaN!!" << std::endl;
1190   }
1191 
1192   TF1 func("tmp", "gaus", mean - 2. * sigma, mean + 2. * sigma);
1193   if (0 == hist->Fit(&func, "QNR")) {  // N: do not blow up file by storing fit!
1194     mean = func.GetParameter(1);
1195     sigma = func.GetParameter(2);
1196 
1197     if (!singleTime) {
1198       // second fit: three sigma of first fit around mean of first fit
1199       func.SetRange(std::max(mean - 3 * sigma, minHist), std::min(mean + 3 * sigma, maxHist));
1200       // I: integral gives more correct results if binning is too wide
1201       // L: Likelihood can treat empty bins correctly (if hist not weighted...)
1202       if (0 == hist->Fit(&func, "Q0LR")) {
1203         if (hist->GetFunction(func.GetName())) {  // Take care that it is later on drawn:
1204           hist->GetFunction(func.GetName())->ResetBit(TF1::kNotDraw);
1205         }
1206       }
1207     }
1208   }
1209 
1210   float res_mean = func.GetParameter(1);
1211   float res_width = func.GetParameter(2);
1212 
1213   float res_mean_err = func.GetParError(1);
1214   float res_width_err = func.GetParError(2);
1215 
1216   Measurement1D resultM(res_mean, res_mean_err);
1217   Measurement1D resultW(res_width, res_width_err);
1218 
1219   statmode::fitParams result = std::make_pair(resultM, resultW);
1220   return result;
1221 }
1222 
1223 //*************************************************************
1224 statmode::fitParams SplitVertexResolution::fitResiduals_v0(TH1* hist)
1225 //*************************************************************
1226 {
1227   float mean = hist->GetMean();
1228   float sigma = hist->GetRMS();
1229 
1230   TF1 func("tmp", "gaus", mean - 1.5 * sigma, mean + 1.5 * sigma);
1231   if (0 == hist->Fit(&func, "QNR")) {  // N: do not blow up file by storing fit!
1232     mean = func.GetParameter(1);
1233     sigma = func.GetParameter(2);
1234     // second fit: three sigma of first fit around mean of first fit
1235     func.SetRange(mean - 2 * sigma, mean + 2 * sigma);
1236     // I: integral gives more correct results if binning is too wide
1237     // L: Likelihood can treat empty bins correctly (if hist not weighted...)
1238     if (0 == hist->Fit(&func, "Q0LR")) {
1239       if (hist->GetFunction(func.GetName())) {  // Take care that it is later on drawn:
1240         hist->GetFunction(func.GetName())->ResetBit(TF1::kNotDraw);
1241       }
1242     }
1243   }
1244 
1245   float res_mean = func.GetParameter(1);
1246   float res_width = func.GetParameter(2);
1247 
1248   float res_mean_err = func.GetParError(1);
1249   float res_width_err = func.GetParError(2);
1250 
1251   Measurement1D resultM(res_mean, res_mean_err);
1252   Measurement1D resultW(res_width, res_width_err);
1253 
1254   statmode::fitParams result = std::make_pair(resultM, resultW);
1255   return result;
1256 }
1257 
1258 //*************************************************************
1259 template <std::size_t SIZE>
1260 bool SplitVertexResolution::checkBinOrdering(std::array<float, SIZE>& bins)
1261 //*************************************************************
1262 {
1263   int i = 1;
1264 
1265   if (std::is_sorted(bins.begin(), bins.end())) {
1266     return true;
1267   } else {
1268     for (const auto& bin : bins) {
1269       edm::LogInfo("SplitVertexResolution") << "bin: " << i << " : " << bin << std::endl;
1270       i++;
1271     }
1272     edm::LogInfo("SplitVertexResolution") << "--------------------------------" << std::endl;
1273     return false;
1274   }
1275 }
1276 
1277 //define this as a plug-in
1278 DEFINE_FWK_MODULE(SplitVertexResolution);