15 #include "TFitResult.h"
16 #include "TFitResultPtr.h"
17 #include "GaudiKernel/PhysicalConstants.h"
19 using namespace IDPVM;
23 m_detailLevel(detailLevel),
24 m_vx_type_truth(nullptr),
25 m_vx_hs_classification(nullptr),
26 m_vx_nReco_vs_nTruth_inclusive(nullptr),
27 m_vx_nReco_vs_nTruth_matched(nullptr),
28 m_vx_nReco_vs_nTruth_merged(nullptr),
29 m_vx_nReco_vs_nTruth_split(nullptr),
30 m_vx_nReco_vs_nTruth_fake(nullptr),
31 m_vx_nReco_vs_nTruth_dummy(nullptr),
32 m_vx_nReco_vs_nTruth_clean(nullptr),
33 m_vx_nReco_vs_nTruth_lowpu(nullptr),
34 m_vx_nReco_vs_nTruth_highpu(nullptr),
35 m_vx_nReco_vs_nTruth_hssplit(nullptr),
36 m_vx_nReco_vs_nTruth_none(nullptr),
37 m_vx_hs_reco_eff(nullptr),
38 m_vx_hs_sel_eff(nullptr),
39 m_vx_hs_sel_eff_vs_nReco(nullptr),
40 m_vx_hs_reco_sel_eff(nullptr),
41 m_vx_hs_sel_eff_dist(nullptr),
42 m_vx_hs_sel_eff_mu(nullptr),
43 m_vx_hs_sel_eff_dist_vs_nReco(nullptr),
44 m_vx_hs_reco_eff_vs_ntruth(nullptr),
45 m_vx_hs_sel_eff_vs_ntruth(nullptr),
46 m_vx_hs_reco_sel_eff_vs_ntruth(nullptr),
47 m_vx_hs_reco_long_reso(nullptr),
48 m_vx_hs_reco_trans_reso(nullptr),
51 m_resHelper_PUdensity_hsVxTruthLong(nullptr),
52 m_resolution_vs_PUdensity_hsVxTruthLong(nullptr),
53 m_resmean_vs_PUdensity_hsVxTruthLong(nullptr),
54 m_resHelper_PUdensity_hsVxTruthTransv(nullptr),
55 m_resolution_vs_PUdensity_hsVxTruthTransv(nullptr),
56 m_resmean_vs_PUdensity_hsVxTruthTransv(nullptr),
58 m_vx_hs_z_pull(nullptr),
59 m_vx_hs_y_pull(nullptr),
60 m_vx_hs_x_pull(nullptr),
61 m_vx_all_z_pull(nullptr),
62 m_vx_all_y_pull(nullptr),
63 m_vx_all_x_pull(nullptr),
64 m_vx_hs_z_res(nullptr),
65 m_vx_hs_y_res(nullptr),
66 m_vx_hs_x_res(nullptr),
67 m_vx_all_z_res(nullptr),
68 m_vx_all_y_res(nullptr),
69 m_vx_all_x_res(nullptr),
70 m_vx_all_truth_z_res_vs_PU(nullptr),
71 m_vx_all_truth_x_res_vs_PU(nullptr),
72 m_vx_all_truth_y_res_vs_PU(nullptr),
73 m_vx_all_truth_z_pull_vs_PU(nullptr),
74 m_vx_all_truth_x_pull_vs_PU(nullptr),
75 m_vx_all_truth_y_pull_vs_PU(nullptr),
76 m_vx_all_truth_z_res_vs_nTrk(nullptr),
77 m_vx_all_truth_x_res_vs_nTrk(nullptr),
78 m_vx_all_truth_y_res_vs_nTrk(nullptr),
79 m_vx_all_truth_z_pull_vs_nTrk(nullptr),
80 m_vx_all_truth_x_pull_vs_nTrk(nullptr),
81 m_vx_all_truth_y_pull_vs_nTrk(nullptr),
82 m_vx_hs_truth_z_res_vs_PU(nullptr),
83 m_vx_hs_truth_x_res_vs_PU(nullptr),
84 m_vx_hs_truth_y_res_vs_PU(nullptr),
85 m_vx_hs_truth_z_pull_vs_PU(nullptr),
86 m_vx_hs_truth_x_pull_vs_PU(nullptr),
87 m_vx_hs_truth_y_pull_vs_PU(nullptr),
88 m_vx_hs_truth_z_res_vs_nTrk(nullptr),
89 m_vx_hs_truth_x_res_vs_nTrk(nullptr),
90 m_vx_hs_truth_y_res_vs_nTrk(nullptr),
91 m_vx_hs_truth_z_pull_vs_nTrk(nullptr),
92 m_vx_hs_truth_x_pull_vs_nTrk(nullptr),
93 m_vx_hs_truth_y_pull_vs_nTrk(nullptr),
96 m_vx_ntracks_matched(nullptr),
97 m_vx_ntracks_merged(nullptr),
98 m_vx_ntracks_split(nullptr),
99 m_vx_ntracks_HS_matched(nullptr),
100 m_vx_ntracks_HS_merged(nullptr),
101 m_vx_ntracks_HS_split(nullptr),
102 m_vx_ntracks_ALL_matched(nullptr),
103 m_vx_ntracks_ALL_merged(nullptr),
104 m_vx_ntracks_ALL_split(nullptr),
105 m_vx_sumpT_matched(nullptr),
106 m_vx_sumpT_merged(nullptr),
107 m_vx_sumpT_split(nullptr),
108 m_vx_sumpT_HS_matched(nullptr),
109 m_vx_sumpT_HS_merged(nullptr),
110 m_vx_sumpT_HS_split(nullptr),
112 m_vx_z_asym_matched(nullptr),
113 m_vx_z_asym_merged(nullptr),
114 m_vx_z_asym_split(nullptr),
115 m_vx_z_asym_HS_matched(nullptr),
116 m_vx_z_asym_HS_merged(nullptr),
117 m_vx_z_asym_HS_split(nullptr),
118 m_vx_z_asym_weighted_matched(nullptr),
119 m_vx_z_asym_weighted_merged(nullptr),
120 m_vx_z_asym_weighted_split(nullptr),
121 m_vx_z_asym_weighted_HS_matched(nullptr),
122 m_vx_z_asym_weighted_HS_merged(nullptr),
123 m_vx_z_asym_weighted_HS_split(nullptr),
125 m_vx_track_weight_matched(nullptr),
126 m_vx_track_weight_merged(nullptr),
127 m_vx_track_weight_split(nullptr),
128 m_vx_track_weight_HS_matched(nullptr),
129 m_vx_track_weight_HS_merged(nullptr),
130 m_vx_track_weight_HS_split(nullptr),
132 m_vx_normalised_track_weight_matched(nullptr),
133 m_vx_normalised_track_weight_merged(nullptr),
134 m_vx_normalised_track_weight_split(nullptr),
135 m_vx_normalised_track_weight_HS_matched(nullptr),
136 m_vx_normalised_track_weight_HS_merged(nullptr),
137 m_vx_normalised_track_weight_HS_split(nullptr),
139 m_vx_chi2Over_ndf_matched(nullptr),
140 m_vx_chi2Over_ndf_merged(nullptr),
141 m_vx_chi2Over_ndf_split(nullptr),
142 m_vx_chi2Over_ndf_HS_matched(nullptr),
143 m_vx_chi2Over_ndf_HS_merged(nullptr),
144 m_vx_chi2Over_ndf_HS_split(nullptr),
146 m_vx_z0_skewness_matched(nullptr),
147 m_vx_z0_skewness_merged(nullptr),
148 m_vx_z0_skewness_split(nullptr),
149 m_vx_z0_skewness_HS_matched(nullptr),
150 m_vx_z0_skewness_HS_merged(nullptr),
151 m_vx_z0_skewness_HS_split(nullptr),
153 m_vx_z0_kurtosis_matched(nullptr),
154 m_vx_z0_kurtosis_merged(nullptr),
155 m_vx_z0_kurtosis_split(nullptr),
156 m_vx_z0_kurtosis_HS_matched(nullptr),
157 m_vx_z0_kurtosis_HS_merged(nullptr),
158 m_vx_z0_kurtosis_HS_split(nullptr),
160 m_vx_sumpT_ALL_matched(nullptr),
161 m_vx_sumpT_ALL_merged(nullptr),
162 m_vx_sumpT_ALL_split(nullptr),
163 m_vx_z_asym_ALL_matched(nullptr),
164 m_vx_z_asym_ALL_merged(nullptr),
165 m_vx_z_asym_ALL_split(nullptr),
166 m_vx_z_asym_weighted_ALL_matched(nullptr),
167 m_vx_z_asym_weighted_ALL_merged(nullptr),
168 m_vx_z_asym_weighted_ALL_split(nullptr),
170 m_vx_track_weight_ALL_matched(nullptr),
171 m_vx_track_weight_ALL_merged(nullptr),
172 m_vx_track_weight_ALL_split(nullptr),
174 m_vx_normalised_track_weight_ALL_matched(nullptr),
175 m_vx_normalised_track_weight_ALL_merged(nullptr),
176 m_vx_normalised_track_weight_ALL_split(nullptr),
178 m_vx_chi2Over_ndf_ALL_matched(nullptr),
179 m_vx_chi2Over_ndf_ALL_merged(nullptr),
180 m_vx_chi2Over_ndf_ALL_split(nullptr),
182 m_vx_z0_skewness_ALL_matched(nullptr),
183 m_vx_z0_skewness_ALL_merged(nullptr),
184 m_vx_z0_skewness_ALL_split(nullptr),
186 m_vx_z0_kurtosis_ALL_matched(nullptr),
187 m_vx_z0_kurtosis_ALL_merged(nullptr),
188 m_vx_z0_kurtosis_ALL_split(nullptr),
190 m_vx_nVertices_ALL_matched(nullptr),
191 m_vx_nVertices_ALL_merged(nullptr),
192 m_vx_nVertices_ALL_split(nullptr),
193 m_vx_nVertices_ALL_fake(nullptr),
194 m_vx_nVertices_HS_matched(nullptr),
195 m_vx_nVertices_HS_merged(nullptr),
196 m_vx_nVertices_HS_split(nullptr),
197 m_vx_nVertices_HS_fake(nullptr),
198 m_vx_nVertices_matched(nullptr),
199 m_vx_nVertices_merged(nullptr),
200 m_vx_nVertices_split(nullptr),
201 m_vx_nVertices_fake(nullptr),
203 m_vx_all_dz(nullptr),
204 m_vx_hs_mindz(nullptr),
206 m_vx_PUdensity(nullptr),
207 m_vx_nTruth(nullptr),
208 m_vx_nTruth_vs_PUdensity(nullptr)
413 float sumPtMax = -1.;
416 for (
const auto& vtx : recoVertices.
stdcont()) {
419 for (
size_t i = 0;
i < vtx->nTrackParticles();
i++) {
420 trackTmp = vtx->trackParticle(
i);
425 if (sumPtTmp > sumPtMax) {
435 template<
typename U,
typename V>
437 return (
std::pow((vtx1->x() - vtx2->x()), 2) +
std::pow((vtx1->y() - vtx2->y()), 2) +
std::pow((vtx1->z() - vtx2->z()), 2));
441 float radialWindow2 =
std::pow(radialWindow, 2);
442 int nTracksInWindow = 0;
443 float localPUDensity;
445 for (
const auto& vtx : truthHSVertices) {
446 if (vtx != vtxOfInterest) {
448 if (radialDiff2 < radialWindow2) {
449 nTracksInWindow += 1;
453 for (
const auto& vtx : truthPUVertices) {
454 if (vtx != vtxOfInterest) {
456 if (radialDiff2 < radialWindow2) {
457 nTracksInWindow += 1;
461 localPUDensity = (
float)(nTracksInWindow) / (2 * radialWindow);
462 return localPUDensity;
466 return std::sqrt(recoVtx->covariancePosition()(2, 2));
470 float x = recoVtx->
x();
471 float y = recoVtx->
y();
472 float xErr2 = recoVtx->covariancePosition()(0, 0);
473 float yErr2 = recoVtx->covariancePosition()(1, 1);
474 float xyCov = recoVtx->covariancePosition()(0, 1);
476 return std::sqrt(
std::pow(
x, 2) / r2 * xErr2 +
std::pow(
y, 2) / r2 * yErr2 + 2 *
x *
y / r2 * xyCov);
482 TH1* projHist =
nullptr;
484 TFitResultPtr fitResult;
487 double itr_rms = -1.;
490 for (
int i = 1;
i < resoHist2D->GetNbinsX() + 1;
i++) {
492 projHist = resoHist2D->ProjectionY(
"projectionY",
i,
i);
494 if (projHist->GetEntries() == 0.) {
495 resoHist->SetBinContent(
i, 0.);
496 resoHist->SetBinError(
i, 0.);
502 fitResult = projHist->Fit(
"gaus",
"QS0");
503 if (!fitResult.Get()) {
506 resoHist->SetBinContent(
i, 0.);
507 resoHist->SetBinError(
i, 0.);
510 mean = fitResult->Parameter(1);
511 rms = fitResult->Parameter(2);
517 fitResult = projHist->Fit(
"gaus",
"QS0");
518 if (!fitResult.Get()) {
523 itr_rms = fitResult->Parameter(2);
524 itr_rms_err = fitResult->ParError(2);
526 if ((fabs(itr_rms -
rms) < 0.0001) || (safety_counter == 5)) {
531 mean = fitResult->Parameter(1);
536 resoHist->SetBinContent(
i, itr_rms);
537 resoHist->SetBinError(
i, itr_rms_err);
548 ATH_MSG_WARNING(
"TruthEventMatchingInfos DECORATOR not available -- returning nullptr!");
551 const std::vector<InDetVertexTruthMatchUtils::VertexTruthMatchInfo>& truthInfos = truthMatchingInfos(*recoVtx);
552 if (!truthInfos.empty()) {
556 if (truthEventLink.
isValid()) {
561 while(!truthVtx && i_vtx<n_vtx){
569 ATH_MSG_WARNING(
"TruthEventMatchingInfos DECORATOR yields empty vector -- returning nullptr!");
573 ATH_MSG_WARNING(
"TruthEventMatchingInfos DECORATOR yields empty vector -- returning nullptr!");
582 float diff_z=
vertex.z()-tvrt->
z();
594 if (
vertex.hasValidTime()) {
596 float err_time =
vertex.timeResolution();
608 matchType = recoVtxMatchTypeInfo(
vertex);
609 ATH_MSG_DEBUG(
"VERTEX DECORATOR ======= " << matchType <<
", with nTRACKS === " <<
vertex.nTrackParticles() <<
", vertex index = " <<
vertex.index() <<
" AT (x, y, z) = (" <<
vertex.x() <<
", " <<
vertex.y() <<
", " <<
vertex.z() <<
")");
613 ATH_MSG_WARNING(
"VertexMatchType DECORATOR seems to be available, but may be broken ===========");
617 ATH_MSG_WARNING(
"VertexMatchType DECORATOR is NOT available ===========");
627 int nTruthVertices = (
int)(truthHSVertices.size() + truthPUVertices.size());
628 int nRecoVertices = (
int)vertexContainer.
size()-1;
634 std::map<InDetVertexTruthMatchUtils::VertexMatchType, int> breakdown = {};
641 if (!recoHardScatter){
642 ATH_MSG_INFO(
"No recoHardScatter vertex - not filling vertex truth matching.");
651 if (!truthHSVertices.empty()) {
652 if (truthHSVertices.size() != 1) {
653 ATH_MSG_WARNING(
"Size of truth HS vertex vector is >1 -- only using the first one in the vector.");
655 truthHSVtx = truthHSVertices.at(0);
660 ATH_MSG_WARNING(
"Size of truth HS vertex vector is 0 -- assuming truth HS vertex to NOT be reconstructed.");
664 float localPUDensity =
getLocalPUDensity(truthHSVtx, truthHSVertices, truthPUVertices);
670 if (!bestRecoHSVtx_truth){
671 ATH_MSG_INFO(
"No bestRecoHS vertex - not filling vertex truth matching.");
681 bool truthHSVtxRecoed =
false;
683 float truthRecoRadialDiff2 = -1.;
687 truthRecoRadialDiff2 =
getRadialDiff2(bestRecoHSVtx_truth, truthHSVtx);
688 if (truthRecoRadialDiff2 < minTruthRecoRadialDiff2) {
689 truthHSVtxRecoed =
true;
690 minTruthRecoRadialDiff2 = truthRecoRadialDiff2;
695 float number_matched = 0;
696 float number_merged = 0;
697 float number_split = 0;
698 float number_fake = 0;
699 float number_matched_HS = 0;
700 float number_merged_HS = 0;
701 float number_split_HS = 0;
702 float number_fake_HS = 0;
703 float number_matched_PU = 0;
704 float number_merged_PU = 0;
705 float number_split_PU = 0;
706 float number_fake_PU = 0;
709 float vx_hs_mindz=9999.;
710 float min_fabs_dz = 9999.;
723 matchType = recoVtxMatchTypeInfo(*
vertex);
724 breakdown[matchType] += 1;
728 if(!matchVertex)
continue;
729 float residual_z = matchVertex->
z() -
vertex->z();
730 float residual_x = matchVertex->
x() -
vertex->x();
731 float residual_y = matchVertex->
y() -
vertex->y();
736 localPUDensity =
getLocalPUDensity(matchVertex, truthHSVertices, truthPUVertices);
772 float minpt = 20000 ;
780 float weighted_sumDZ = 0;
781 float weighted_deltaZ = 0;
782 float weighted_modsumDZ = 0;
783 float weighted_z_asym =0;
786 std::vector<float> track_deltaZ;
787 std::vector<float> track_deltaPt;
788 std::vector<float> track_deltaZ_weighted;
790 for (
size_t i = 0;
i <
vertex->nTrackParticles();
i++) {
791 trackTmp =
vertex->trackParticle(
i);
798 track_deltaZ.push_back(
deltaZ);
800 float trk_weight =
vertex->trackWeight(
i);
801 weighted_deltaZ =
deltaZ*trk_weight;
804 modsumDZ = modsumDZ + std::abs(
deltaZ);
805 weighted_sumDZ = weighted_sumDZ + weighted_deltaZ;
806 weighted_modsumDZ = weighted_modsumDZ + std::abs(weighted_deltaZ);
811 z_asym = sumDZ/modsumDZ;
813 if (weighted_modsumDZ >0) {
814 weighted_z_asym = weighted_sumDZ/weighted_modsumDZ;
820 mean_Dz=sumDZ/track_deltaZ.size();
821 double number_tracks =0;
822 number_tracks = track_deltaZ.size();
830 for (
auto i : track_deltaZ) {
832 z_zbar = (
i - mean_Dz);
833 z_var =(z_var + z_zbar*z_zbar);
834 z_skew =(z_skew + z_zbar*z_zbar*z_zbar);
835 z_kurt =(z_kurt + z_zbar*z_zbar*z_zbar*z_zbar);
838 z_var = z_var/(number_tracks -1);
839 z_sd = std::sqrt(z_var);
840 z_skew = z_skew/((number_tracks -1)*z_sd*z_sd*z_sd);
841 z_kurt = z_kurt/((number_tracks -1)*z_sd*z_sd*z_sd*z_sd);
847 if (
vertex == bestRecoHSVtx_truth) {
858 for (
const float& trkWeight :
vertex->trackWeights()) {
875 for (
const float& trkWeight :
vertex->trackWeights()) {
889 for (
const float& trkWeight :
vertex->trackWeights()) {
898 if (
vertex == bestRecoHSVtx_truth) {
907 for (
const float& trkWeight :
vertex->trackWeights()) {
923 for (
const float& trkWeight :
vertex->trackWeights()) {
937 for (
const float& trkWeight :
vertex->trackWeights()) {
946 if (
vertex == bestRecoHSVtx_truth) {
954 for (
const float& trkWeight :
vertex->trackWeights()) {
969 for (
const float& trkWeight :
vertex->trackWeights()) {
985 for (
const float& trkWeight :
vertex->trackWeights()) {
997 if (
vertex == bestRecoHSVtx_truth) {
1001 number_matched_PU++;
1003 if (sumPt > minpt) {
1009 if (
vertex == bestRecoHSVtx_truth) {
1015 if (sumPt > minpt) {
1021 if (
vertex == bestRecoHSVtx_truth) {
1027 if (sumPt > minpt) {
1033 if (
vertex == bestRecoHSVtx_truth) {
1039 if (sumPt > minpt) {
1049 if (
vertex == bestRecoHSVtx_truth) {
1084 if (sumPt > minpt) {
1104 float absd_hs_dz = std::abs(bestRecoHSVtx_truth->
z() -
vertex->z());
1105 if(bestRecoHSVtx_truth !=
vertex && absd_hs_dz < min_fabs_dz) {
1106 min_fabs_dz = absd_hs_dz;
1107 vx_hs_mindz = bestRecoHSVtx_truth->
z() -
vertex->z();
1110 for (
const auto& vertex2 : vertexContainer.
stdcont()) {
1112 if(vertex2 ==
vertex)
continue;
1136 if (!truthHSVertices.empty()) {
1137 if (truthHSVtxRecoed) {
1138 float residual_z = truthHSVtx->
z() - bestRecoHSVtx_truth->
z();
1139 float residual_r = std::sqrt(
std::pow(truthHSVtx->
x() - bestRecoHSVtx_truth->
x(), 2) +
std::pow(truthHSVtx->
y() - bestRecoHSVtx_truth->
y(), 2));
1140 float residual_x = truthHSVtx->
x() - bestRecoHSVtx_truth->
x();
1141 float residual_y = truthHSVtx->
y() - bestRecoHSVtx_truth->
y();
1152 const AmgSymMatrix(3)& covariance = bestRecoHSVtx_truth->covariancePosition();