File indexing completed on 2024-04-06 11:57:08
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
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
0026 #include "TTree.h"
0027 #include "TProfile.h"
0028 #include "TF1.h"
0029 #include "TMath.h"
0030
0031
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
0066
0067
0068 namespace statmode {
0069 using fitParams = std::pair<Measurement1D, Measurement1D>;
0070 }
0071
0072
0073
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
0106 int ievt;
0107 int itrks;
0108
0109
0110 const int compressionSettings_;
0111
0112
0113 bool storeNtuple_;
0114
0115
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");
0126 edm::EDGetTokenT<edm::TriggerResults> triggerResultsToken_;
0127
0128
0129 edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> transientTrackBuilderToken_;
0130 edm::ESGetToken<RunInfo, RunInfoRcd> runInfoToken_;
0131
0132 double minVtxNdf_;
0133 double minVtxWgt_;
0134
0135
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
0191
0192 TH1D* tksByTrigger_;
0193 TH1D* evtsByTrigger_;
0194
0195
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
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
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_;
0261 std::array<float, nTrackBins_ + 1> myNTrack_bins_;
0262
0263 static const int nVtxBins_ = 60;
0264 const double nVisVtxBins_;
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
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
0300
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
0307
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
0353 void SplitVertexResolution::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0354 using namespace edm;
0355
0356
0357
0358
0359 engine_.seed(iEvent.id().event() + (iEvent.id().luminosityBlock() << 10) + (iEvent.id().run() << 20));
0360
0361
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
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
0389
0390 beginEvent();
0391
0392
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
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 const reco::Vertex& iPV = *pvIt;
0427 counter++;
0428 if (iPV.isFake())
0429 continue;
0430 noFakecounter++;
0431
0432
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
0450 std::sort(allTracks.begin(), allTracks.end(), mysorter);
0451
0452 int ntrks = allTracks.size();
0453 h_ntrks->Fill(ntrks);
0454
0455
0456 uint even_ntrks;
0457 ntrks % 2 == 0 ? even_ntrks = ntrks : even_ntrks = ntrks - 1;
0458
0459
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);
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
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
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
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
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
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
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
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
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 }
0666
0667
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
0700 void SplitVertexResolution::beginJob() {
0701 ievt = 0;
0702 itrks = 0;
0703
0704 if (compressionSettings_ > 0) {
0705 outfile_->file().SetCompressionSettings(compressionSettings_);
0706 }
0707
0708
0709 TFileDirectory BinningFeatures = outfile_->mkdir("BinningFeatures");
0710
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
0715 h_profileBinnings->GetXaxis()->SetBinLabel(3, "n. track bins");
0716
0717 h_profileBinnings->GetXaxis()->SetBinLabel(4, "n. vertices bins");
0718
0719
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
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
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
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
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
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
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
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
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
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
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
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
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
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")) {
1194 mean = func.GetParameter(1);
1195 sigma = func.GetParameter(2);
1196
1197 if (!singleTime) {
1198
1199 func.SetRange(std::max(mean - 3 * sigma, minHist), std::min(mean + 3 * sigma, maxHist));
1200
1201
1202 if (0 == hist->Fit(&func, "Q0LR")) {
1203 if (hist->GetFunction(func.GetName())) {
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")) {
1232 mean = func.GetParameter(1);
1233 sigma = func.GetParameter(2);
1234
1235 func.SetRange(mean - 2 * sigma, mean + 2 * sigma);
1236
1237
1238 if (0 == hist->Fit(&func, "Q0LR")) {
1239 if (hist->GetFunction(func.GetName())) {
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
1278 DEFINE_FWK_MODULE(SplitVertexResolution);
