22 #include <boost/range/adaptor/indexed.hpp>
68 using fitParams = std::pair<Measurement1D, Measurement1D>;
91 template <std::
size_t SIZE>
94 unsigned int theNOfBins,
252 std::array<float, nPtBins_ + 1>
mypT_bins_ = PVValHelper::makeLogBins<float, nPtBins_>(1., 1
e3);
264 : compressionSettings_(iConfig.getUntrackedParameter<
int>(
"compressionSettings", -1)),
265 storeNtuple_(iConfig.getParameter<
bool>(
"storeNtuple")),
266 intLumi_(iConfig.getUntrackedParameter<double>(
"intLumi", 0.)),
267 debug_(iConfig.getUntrackedParameter<
bool>(
"Debug",
false)),
268 pvsTag_(iConfig.getParameter<
edm::
InputTag>(
"vtxCollection")),
270 tracksTag_(iConfig.getParameter<
edm::
InputTag>(
"trackCollection")),
273 transientTrackBuilderToken_(
276 minVtxNdf_(iConfig.getUntrackedParameter<double>(
"minVertexNdf")),
277 minVtxWgt_(iConfig.getUntrackedParameter<double>(
"minVertexMeanWeight")),
278 runControl_(iConfig.getUntrackedParameter<
bool>(
"runControl",
false)) {
281 std::vector<unsigned int> defaultRuns;
282 defaultRuns.push_back(0);
305 bool passesRunControl =
false;
309 if (
iEvent.eventAuxiliary().run() == runControlNumber) {
312 <<
" run number: " <<
iEvent.eventAuxiliary().run() <<
" keeping run:" << runControlNumber;
314 passesRunControl =
true;
318 if (!passesRunControl)
347 int nOfflineVtx = pvtx.size();
354 for (
int itrig = 0; itrig != ntrigs; ++itrig) {
365 int noFakecounter = 0;
368 for (
auto pvIt = pvtx.cbegin(); pvIt != pvtx.cend(); ++pvIt) {
383 if (trki->isNonnull()) {
400 ntrks % 2 == 0 ? even_ntrks = ntrks : even_ntrks = ntrks - 1;
403 for (
uint tracksIt = 0; tracksIt < even_ntrks; tracksIt = tracksIt + 2) {
406 auto dis = std::uniform_int_distribution<>(0, 1);
409 groupOne.push_back(firstTrk);
410 groupTwo.push_back(secondTrk);
412 groupOne.push_back(secondTrk);
413 groupTwo.push_back(firstTrk);
417 if (!(groupOne.size() >= 2 && groupTwo.size() >= 2))
424 float sumPt = 0, sumPt1 = 0, sumPt2 = 0, avgSumPt = 0;
427 std::vector<reco::TransientTrack> groupOne_ttks;
428 groupOne_ttks.clear();
429 for (
auto itrk = groupOne.cbegin(); itrk != groupOne.cend(); itrk++) {
431 groupOne_ttks.push_back(tmpTransientTrack);
432 sumPt1 += itrk->pt();
443 std::vector<reco::TransientTrack> groupTwo_ttks;
444 groupTwo_ttks.clear();
445 for (
auto itrk = groupTwo.cbegin(); itrk != groupTwo.cend(); itrk++) {
447 groupTwo_ttks.push_back(tmpTransientTrack);
448 sumPt2 += itrk->pt();
453 avgSumPt = (sumPt1 + sumPt2) / 2.;
479 int half_trks = twoPV.
nTracks();
481 const double invSqrt2 = 1. /
std::sqrt(2.);
483 double deltaX = (twoPV.
x() - onePV.
x());
484 double deltaY = (twoPV.
y() - onePV.
y());
485 double deltaZ = (twoPV.
z() - onePV.
z());
487 double resX = deltaX * invSqrt2;
488 double resY = deltaY * invSqrt2;
489 double resZ = deltaZ * invSqrt2;
509 for (
int ipTBin = 0; ipTBin <
nPtBins_; ipTBin++) {
513 if (avgSumPt >= pTF && avgSumPt < pTL) {
526 for (
int inTrackBin = 0; inTrackBin <
nTrackBins_; inTrackBin++) {
530 if (ntrks >= nTrackF && ntrks < nTrackL) {
543 for (
int inVtxBin = 0; inVtxBin <
nVtxBins_; inVtxBin++) {
550 if (nOfflineVtx == inVtxBin) {
567 h_PVCL_subVtx1->Fill(TMath::Prob(pvOne.totalChiSquared(), (
int)(pvOne.degreesOfFreedom())));
568 h_PVCL_subVtx2->Fill(TMath::Prob(pvTwo.totalChiSquared(), (
int)(pvTwo.degreesOfFreedom())));
603 thePV.
CL_subVtx1 = TMath::Prob(pvOne.totalChiSquared(), (
int)(pvOne.degreesOfFreedom()));
604 thePV.
CL_subVtx2 = TMath::Prob(pvTwo.totalChiSquared(), (
int)(pvTwo.degreesOfFreedom()));
629 unsigned int RunNumber_ =
run.run();
635 const time_t start_time = times.first / 1.0e+6;
637 << RunNumber_ <<
" has start time: " << times.first <<
" - " << times.second << std::endl;
639 <<
"human readable time: " << std::asctime(std::gmtime(&start_time)) << std::endl;
657 EventFeatures.
make<TH1F>(
"h_lumiFromConfig",
"luminosity from config;;luminosity of present run", 1, -0.5, 0.5);
661 "run number from config;;run number (from configuration)",
673 edm::LogError(
"SplitVertexResolution") <<
" Warning - the vector of pT bins is not ordered " << std::endl;
677 edm::LogError(
"SplitVertexResolution") <<
" Warning -the vector of n. tracks bins is not ordered " << std::endl;
681 edm::LogError(
"SplitVertexResolution") <<
" Warning -the vector of n. vertices bins is not ordered " << std::endl;
741 h_runNumber =
outfile_->
make<TH1F>(
"h_runNumber",
"run number;run number;n_{events}", 100000, 250000., 350000.);
746 outfile_->
make<TH1I>(
"h_nRealVertices",
"n. of non-fake vertices;n. vertices; events", 100, 0, 100);
748 outfile_->
make<TH1I>(
"h_nSelectedVertices",
"n. of selected vertices vertices;n. vertices; events", 100, 0, 100);
751 "h_diffX",
"x-coordinate vertex resolution;vertex resolution (x) [#mum];vertices", 100, -300, 300.);
753 "h_diffY",
"y-coordinate vertex resolution;vertex resolution (y) [#mum];vertices", 100, -300, 300.);
755 "h_diffZ",
"z-coordinate vertex resolution;vertex resolution (z) [#mum];vertices", 100, -500, 500.);
758 "h_OrigVertexErrX",
"x-coordinate vertex error;vertex error (x) [#mum];vertices", 300, 0., 300.);
760 "h_OrigVertexErrY",
"y-coordinate vertex error;vertex error (y) [#mum];vertices", 300, 0., 300.);
762 "h_OrigVertexErrZ",
"z-coordinate vertex error;vertex error (z) [#mum];vertices", 500, 0., 500.);
765 "h_errX",
"x-coordinate vertex resolution error;vertex resoltuion error (x) [#mum];vertices", 300, 0., 300.);
767 "h_errY",
"y-coordinate vertex resolution error;vertex resolution error (y) [#mum];vertices", 300, 0., 300.);
769 "h_errZ",
"z-coordinate vertex resolution error;vertex resolution error (z) [#mum];vertices", 500, 0., 500.);
771 h_pullX =
outfile_->
make<TH1F>(
"h_pullX",
"x-coordinate vertex pull;vertex pull (x);vertices", 500, -10, 10.);
772 h_pullY =
outfile_->
make<TH1F>(
"h_pullY",
"y-coordinate vertex pull;vertex pull (y);vertices", 500, -10, 10.);
773 h_pullZ =
outfile_->
make<TH1F>(
"h_pullZ",
"z-coordinate vertex pull;vertex pull (z);vertices", 500, -10, 10.);
776 "number of tracks in vertex;vertex multeplicity;vertices",
784 "h_avgSumPt",
"#LT #Sigma p_{T} #GT;#LT #sum p_{T} #GT [GeV];vertices",
mypT_bins_.size() - 1,
mypT_bins_.data());
787 "#Sigma p_{T} sub-vertex 1;#sum p_{T} sub-vertex 1 [GeV];subvertices",
791 "#Sigma p_{T} sub-vertex 2;#sum p_{T} sub-vertex 2 [GeV];subvertices",
795 h_wTrks1 =
outfile_->
make<TH1F>(
"h_wTrks1",
"weight of track for sub-vertex 1;track weight;subvertices", 500, 0., 1.);
796 h_wTrks2 =
outfile_->
make<TH1F>(
"h_wTrks2",
"weithg of track for sub-vertex 2;track weight;subvertices", 500, 0., 1.);
799 "h_minWTrks1",
"minimum weight of track for sub-vertex 1;minimum track weight;subvertices", 500, 0., 1.);
801 "h_minWTrks2",
"minimum weithg of track for sub-vertex 2;minimum track weight;subvertices", 500, 0., 1.);
805 "#chi^{2} probability for sub-vertex 1;Prob(#chi^{2},ndof) for sub-vertex 1;subvertices",
811 "#chi^{2} probability for sub-vertex 2;Prob(#chi^{2},ndof) for sub-vertex 2;subvertices",
819 "x-resolution vs #Sigma p_{T};#sum p_{T}; x vertex resolution [#mum]",
823 "y-resolution vs #Sigma p_{T};#sum p_{T}; y vertex resolution [#mum]",
827 "z-resolution vs #Sigma p_{T};#sum p_{T}; z vertex resolution [#mum]",
832 "x-resolution vs n_{tracks};n_{tracks}; x vertex resolution [#mum]",
836 "y-resolution vs n_{tracks};n_{tracks}; y vertex resolution [#mum]",
840 "z-resolution vs n_{tracks};n_{tracks}; z vertex resolution [#mum]",
845 "x-resolution vs n_{vertices};n_{vertices}; x vertex resolution [#mum]",
849 "y-resolution vs n_{vertices};n_{vertices}; y vertex resolution [#mum]",
853 "z-resolution vs n_{vertices};n_{vertices}; z vertex resolution [#mum]",
860 "p_pullX_vsSumPt",
"x-pull vs #Sigma p_{T};#sum p_{T}; x vertex pull",
mypT_bins_.size() - 1,
mypT_bins_.data());
862 "p_pullY_vsSumPt",
"y-pull vs #Sigma p_{T};#sum p_{T}; y vertex pull",
mypT_bins_.size() - 1,
mypT_bins_.data());
864 "p_pullZ_vsSumPt",
"z-pull vs #Sigma p_{T};#sum p_{T}; z vertex pull",
mypT_bins_.size() - 1,
mypT_bins_.data());
867 "x-pull vs n_{tracks};n_{tracks}; x vertex pull",
871 "y-pull vs n_{tracks};n_{tracks}; y vertex pull",
875 "z-pull vs n_{tracks};n_{tracks}; z vertex pull",
880 "x-pull vs n_{vertices};n_{vertices}; x vertex pull",
884 "y-pull vs n_{vertices};n_{vertices}; y vertex pull",
888 "z-pull vs n_{vertices};n_{vertices}; z vertex pull",
900 unsigned int theNOfBins,
903 TH1F::SetDefaultSumw2(kTRUE);
908 if (resType.Contains(
"pull")) {
913 std::vector<TH1F*>
h;
914 h.reserve(theNOfBins);
916 const char* auxResType = (resType.ReplaceAll(
"_",
"")).Data();
918 for (
unsigned int i = 0;
i < theNOfBins;
i++) {
919 TH1F* htemp =
dir.make<TH1F>(Form(
"histo_%s_%s_plot%i", resType.Data(),
varType.Data(),
i),
920 Form(
"%s vs %s - bin %i;%s;vertices", auxResType,
varType.Data(),
i, auxResType),
932 edm::LogVerbatim(
"SplitVertexResolution") <<
"*******************************" << std::endl;
935 edm::LogVerbatim(
"SplitVertexResolution") <<
"*******************************" << std::endl;
938 edm::LogVerbatim(
"SplitVertexResolution") <<
"firing triggers: " << nFiringTriggers << std::endl;
939 edm::LogVerbatim(
"SplitVertexResolution") <<
"*******************************" << std::endl;
942 "tksByTrigger",
"tracks by HLT path;;% of # traks", nFiringTriggers, -0.5, nFiringTriggers - 0.5);
944 "evtsByTrigger",
"events by HLT path;;% of # events", nFiringTriggers, -0.5, nFiringTriggers - 0.5);
950 double trkpercent = ((it->second).
second) * 100. / double(
itrks);
951 double evtpercent = ((it->second).
first) * 100. / double(
ievt);
954 <<
"HLT path: " << std::setw(60) << std::left << it->first <<
" | events firing: " << std::right << std::setw(8)
955 << (it->second).
first <<
" (" << std::setw(8) <<
std::fixed << std::setprecision(4) << evtpercent <<
"%)"
956 <<
" | tracks collected: " << std::setw(8) << (it->second).
second <<
" (" << std::setw(8) <<
std::fixed
957 << std::setprecision(4) << trkpercent <<
"%)";
972 unsigned int count = 1;
1030 return std::make_pair(
runInfo.m_start_time_ll,
runInfo.m_stop_time_ll);
1037 for (
auto iterator =
h.begin(); iterator !=
h.end(); iterator++) {
1043 float mean_ = myFit.first.value();
1044 float meanErr_ = myFit.first.error();
1045 trendPlot->SetBinContent(
bin, mean_);
1046 trendPlot->SetBinError(
bin, meanErr_);
1050 float width_ = myFit.second.value();
1051 float widthErr_ = myFit.second.error();
1052 trendPlot->SetBinContent(
bin, width_);
1053 trendPlot->SetBinError(
bin, widthErr_);
1059 trendPlot->SetBinContent(
bin, median_);
1060 trendPlot->SetBinError(
bin, medianErr_);
1066 trendPlot->SetBinContent(
bin, mad_);
1067 trendPlot->SetBinError(
bin, madErr_);
1072 <<
"fillTrendPlotByIndex() " << fitPar_ <<
" unknown estimator!" << std::endl;
1082 if (
hist->GetEntries() < 10) {
1083 LogDebug(
"SplitVertexResolution") <<
"hist name: " <<
hist->GetName() <<
" has less than 10 entries" << std::endl;
1087 float maxHist =
hist->GetXaxis()->GetXmax();
1088 float minHist =
hist->GetXaxis()->GetXmin();
1090 float sigma =
hist->GetRMS();
1095 sigma = -minHist + maxHist;
1097 <<
"FitPVResiduals::fitResiduals(): histogram" <<
hist->GetName() <<
" mean or sigma are NaN!!" << std::endl;
1100 TF1
func(
"tmp",
"gaus",
mean - 2. * sigma,
mean + 2. * sigma);
1101 if (0 ==
hist->Fit(&
func,
"QNR")) {
1103 sigma =
func.GetParameter(2);
1110 if (0 ==
hist->Fit(&
func,
"Q0LR")) {
1111 if (
hist->GetFunction(
func.GetName())) {
1112 hist->GetFunction(
func.GetName())->ResetBit(TF1::kNotDraw);
1118 float res_mean =
func.GetParameter(1);
1119 float res_width =
func.GetParameter(2);
1121 float res_mean_err =
func.GetParError(1);
1122 float res_width_err =
func.GetParError(2);
1136 float sigma =
hist->GetRMS();
1138 TF1
func(
"tmp",
"gaus",
mean - 1.5 * sigma,
mean + 1.5 * sigma);
1139 if (0 ==
hist->Fit(&
func,
"QNR")) {
1141 sigma =
func.GetParameter(2);
1146 if (0 ==
hist->Fit(&
func,
"Q0LR")) {
1147 if (
hist->GetFunction(
func.GetName())) {
1148 hist->GetFunction(
func.GetName())->ResetBit(TF1::kNotDraw);
1153 float res_mean =
func.GetParameter(1);
1154 float res_width =
func.GetParameter(2);
1156 float res_mean_err =
func.GetParError(1);
1157 float res_width_err =
func.GetParError(2);
1167 template <std::
size_t SIZE>
1173 if (std::is_sorted(
bins.begin(),
bins.end())) {
1176 for (
const auto&
bin :
bins) {
1177 edm::LogInfo(
"SplitVertexResolution") <<
"bin: " <<
i <<
" : " <<
bin << std::endl;
1180 edm::LogInfo(
"SplitVertexResolution") <<
"--------------------------------" << std::endl;