Skip to content
Snippets Groups Projects
Commit 9919e87f authored by Frank Winklmeier's avatar Frank Winklmeier
Browse files

Merge branch 'r22-idpvm-add-vtx-reso-and-effi-plots' into 'master'

Add HS Vertex Reco/Selection Efficiency and Reso Plots to IDPVM

See merge request !32794
parents 40e4dc24 d01304f6
No related branches found
No related tags found
6 merge requests!58791DataQualityConfigurations: Modify L1Calo config for web display,!46784MuonCondInterface: Enable thread-safety checking.,!46776Updated LArMonitoring config file for WD to match new files produced using MT,!45405updated ART test cron job,!42417Draft: DIRE and VINCIA Base Fragments for Pythia 8.3,!32794Add HS Vertex Reco/Selection Efficiency and Reso Plots to IDPVM
Showing with 457 additions and 126 deletions
......@@ -78,8 +78,8 @@ private:
void fillTrackCutFlow(const asg::AcceptData& accept);
void fillCutFlow(const asg::AcceptData& accept, std::vector<std::string> & names, std::vector<int> & cutFlow);
// Get truth particles into a vector, possibly using the pileup from the event
const std::vector<const xAOD::TruthParticle *> getTruthParticles();
const std::vector<const xAOD::TruthVertex*> getTruthVertices();
const std::vector<const xAOD::TruthParticle *> getTruthParticles() const;
std::pair<const std::vector<const xAOD::TruthVertex*>, const std::vector<const xAOD::TruthVertex*>> getTruthVertices() const;
//
const Trk::TrackParameters* getUnbiasedTrackParameters(const Trk::TrackParameters* trkParameters, const Trk::MeasurementBase* measurement );
......
......@@ -780,6 +780,38 @@
<x title="Number of interactions" n="81" lo="-0.5" hi="80.5"/>
<y title="Number of reco vertices" lo="0.0" hi="60.0"/>
</h>
<h id="vx_hs_reco_eff" type="TEfficiency" title="HS vertex reconstruction efficiency">
<x title="Local PU density [vertices/mm]" n="4" lo="0.0" hi="2.0"/>
<y title="HS vertex reconstruction efficiency" lo="0.8" hi="1.2"/>
</h>
<h id="vx_hs_sel_eff" type="TEfficiency" title="HS vertex selection efficiency (sum[pt,track^2])">
<x title="Local PU density [vertices/mm]" n="4" lo="0.0" hi="2.0"/>
<y title="HS vertex selection efficiency" lo="0.6" hi="1.2"/>
</h>
<h id="vx_hs_reco_long_reso" type="TProfile" title="Reco Longitudinal Resolution">
<x title="Local PU density [vertices/mm]" n="4" lo="0.0" hi="2.0"/>
<y title="Reco longitudinal resolution [mm]" lo="0.0" hi="1.0"/>
</h>
<h id="vx_hs_reco_trans_reso" type="TProfile" title="Reco Transverse Resolution">
<x title="Local PU density [vertices/mm]" n="4" lo="0.0" hi="2.0"/>
<y title="Reco transverse resolution [mm]" lo="0.0" hi="1.0"/>
</h>
<h id="vx_hs_truth_long_reso_vs_PU" type="TH2" title="Truth Longitudinal Resolution vs. PU (2D)">
<x title="Local PU density [vertices/mm]" n="4" lo="0.0" hi="2.0"/>
<y title="|Reco-truth| longitudinal resolution [mm]" n="20" lo="-0.08" hi="0.08"/>
</h>
<h id="vx_hs_truth_trans_reso_vs_PU" type="TH2" title="Truth Transverse Resolution vs. PU (2D)">
<x title="Local PU density [vertices/mm]" n="4" lo="0.0" hi="2.0"/>
<y title="|Reco-truth| transverse resolution [mm]" n="20" lo="0.0" hi="0.02"/>
</h>
<h id="vx_hs_truth_long_reso" type="TH1" title="Truth Longitudinal Resolution">
<x title="Local PU density [vertices/mm]" n="4" lo="0.0" hi="2.0"/>
<y title="|Reco-truth| longitudinal resolution [mm]" lo="0.0" hi="1.0"/>
</h>
<h id="vx_hs_truth_trans_reso" type="TH1" title="Truth Transverse Resolution">
<x title="Local PU density [vertices/mm]" n="4" lo="0.0" hi="2.0"/>
<y title="|Reco-truth| transverse resolution [mm]" lo="0.0" hi="1.0"/>
</h>
<h id="vx_nTracks" type="TH1F" title="Number of tracks at vertex">
<x title="Number of Tracks" n="150" lo="0" hi="150"/>
<y title="Entries"/>
......
......@@ -126,6 +126,14 @@ TProfile vx_nReco_vs_nTruth_lowpu "Number of reco vertices (LOWPU) vs Number of
TProfile vx_nReco_vs_nTruth_highpu "Number of reco vertices (HIGHPU) vs Number of interactions" 81 -0.5 80.5 0.0 60.0 "Number of interactions" "Number of reco vertices" default
TProfile vx_nReco_vs_nTruth_hssplit "Number of reco vertices (HSSPLIT) vs Number of interactions" 81 -0.5 80.5 0.0 60.0 "Number of interactions" "Number of reco vertices" default
TProfile vx_nReco_vs_nTruth_none "Number of reco vertices (NONE) vs Number of interactions" 81 -0.5 80.5 0.0 60.0 "Number of interactions" "Number of reco vertices" default
TEfficiency m_vx_hs_reco_eff "HS vertex reconstruction efficiency" 4 0.0 2.0 0.8 1.2 "Local PU density [vertices/mm]" "HS vertex reconstruction efficiency" default
TEfficiency m_vx_hs_sel_eff "HS vertex selection efficiency (sum[pt,track^2])" 4 0.0 2.0 0.6 1.2 "Local PU density [vertices/mm]" "HS vertex selection efficiency" default
TProfile m_vx_hs_reco_long_reso "Reco Longitudinal Resolution" 4 0.0 2.0 0.0 0.012 "Local PU density [vertices/mm]" "Reco longitudinal resolution [mm]" default
TProfile m_vx_hs_reco_trans_reso "Reco Transverse Resolution" 4 0.0 2.0 0.0 0.012 "Local PU density [vertices/mm]" "Reco transverse resolution [mm]" default
TH2F vx_hs_truth_long_reso_vs_PU "Truth Longitudinal Resolution vs. PU (2D)" 4 0.0 2.0 20 -0.08 0.08 "Local PU density [vertices/mm]" "|Reco-truth| longitudinal resolution [mm]" default
TH2F vx_hs_truth_trans_reso_vs_PU "Truth Transverse Resolution vs. PU (2D)" 4 0.0 2.0 20 0.0 0.02 "Local PU density [vertices/mm]" "|Reco-truth| transverse resolution [mm]" default
TH1F m_vx_hs_truth_long_reso "Truth Longitudinal Resolution" 4 0.0 2.0 "Local PU density [vertices/mm]" "(Reco-truth) longitudinal resolution [mm]" default
TH1F m_vx_hs_truth_trans_reso "Truth Transverse Resolution" 4 0.0 2.0 "Local PU density [vertices/mm]" "(Reco-truth) transverse resolution [mm]" default
TH1F vx_nTracks "Number of tracks at vertex" 150 0 150 "Number of Tracks" "Entries" default
TH1F vx_track_weights "Weights of tracks at vertex" 100 0. 10.0 "Weight" "Entries" default
TH1F vx_track_pt "Tracks at vertex p_{T}" 100 0 20. "p_{T} (GeV)" "Entries" default
......
......@@ -11,8 +11,9 @@
using namespace IDPVM;
InDetPerfPlot_VertexTruthMatching::InDetPerfPlot_VertexTruthMatching(InDetPlotBase* pParent, const std::string& sDir) :
InDetPerfPlot_VertexTruthMatching::InDetPerfPlot_VertexTruthMatching(InDetPlotBase* pParent, const std::string& sDir, const int iDetailLevel) :
InDetPlotBase(pParent, sDir),
m_iDetailLevel(iDetailLevel),
m_vx_type_truth(nullptr),
m_vx_hs_classification(nullptr),
m_vx_nReco_vs_nTruth_inclusive(nullptr),
......@@ -25,39 +26,211 @@ InDetPerfPlot_VertexTruthMatching::InDetPerfPlot_VertexTruthMatching(InDetPlotBa
m_vx_nReco_vs_nTruth_lowpu(nullptr),
m_vx_nReco_vs_nTruth_highpu(nullptr),
m_vx_nReco_vs_nTruth_hssplit(nullptr),
m_vx_nReco_vs_nTruth_none(nullptr)
m_vx_nReco_vs_nTruth_none(nullptr),
m_vx_hs_reco_eff(nullptr),
m_vx_hs_sel_eff(nullptr),
m_vx_hs_reco_long_reso(nullptr),
m_vx_hs_reco_trans_reso(nullptr),
m_vx_hs_truth_long_reso_vs_PU(nullptr),
m_vx_hs_truth_trans_reso_vs_PU(nullptr),
m_vx_hs_truth_long_reso(nullptr),
m_vx_hs_truth_trans_reso(nullptr)
{
// nop
}
void
InDetPerfPlot_VertexTruthMatching::initializePlots() {
void InDetPerfPlot_VertexTruthMatching::initializePlots() {
IDPVM_BOOK(m_vx_type_truth);
IDPVM_BOOK(m_vx_hs_classification);
IDPVM_BOOK(m_vx_nReco_vs_nTruth_inclusive);
IDPVM_BOOK(m_vx_nReco_vs_nTruth_matched);
IDPVM_BOOK(m_vx_nReco_vs_nTruth_merged);
IDPVM_BOOK(m_vx_nReco_vs_nTruth_split);
IDPVM_BOOK(m_vx_nReco_vs_nTruth_fake);
IDPVM_BOOK(m_vx_nReco_vs_nTruth_dummy);
IDPVM_BOOK(m_vx_nReco_vs_nTruth_clean);
IDPVM_BOOK(m_vx_nReco_vs_nTruth_lowpu);
IDPVM_BOOK(m_vx_nReco_vs_nTruth_highpu);
IDPVM_BOOK(m_vx_nReco_vs_nTruth_hssplit);
IDPVM_BOOK(m_vx_nReco_vs_nTruth_none);
if (m_iDetailLevel >= 200) {
IDPVM_BOOK(m_vx_hs_classification);
IDPVM_BOOK(m_vx_nReco_vs_nTruth_inclusive);
IDPVM_BOOK(m_vx_nReco_vs_nTruth_matched);
IDPVM_BOOK(m_vx_nReco_vs_nTruth_merged);
IDPVM_BOOK(m_vx_nReco_vs_nTruth_split);
IDPVM_BOOK(m_vx_nReco_vs_nTruth_fake);
IDPVM_BOOK(m_vx_nReco_vs_nTruth_dummy);
IDPVM_BOOK(m_vx_nReco_vs_nTruth_clean);
IDPVM_BOOK(m_vx_nReco_vs_nTruth_lowpu);
IDPVM_BOOK(m_vx_nReco_vs_nTruth_highpu);
IDPVM_BOOK(m_vx_nReco_vs_nTruth_hssplit);
IDPVM_BOOK(m_vx_nReco_vs_nTruth_none);
IDPVM_BOOK(m_vx_hs_reco_eff);
IDPVM_BOOK(m_vx_hs_sel_eff);
IDPVM_BOOK(m_vx_hs_reco_long_reso);
IDPVM_BOOK(m_vx_hs_reco_trans_reso);
IDPVM_BOOK(m_vx_hs_truth_long_reso);
IDPVM_BOOK(m_vx_hs_truth_trans_reso);
IDPVM_BOOK(m_vx_hs_truth_long_reso_vs_PU);
IDPVM_BOOK(m_vx_hs_truth_trans_reso_vs_PU);
}
}
const xAOD::Vertex* InDetPerfPlot_VertexTruthMatching::getHSRecoVertexSumPt2(const xAOD::VertexContainer& recoVertices) const {
const xAOD::Vertex* recoHSVertex = nullptr;
float sumPtMax = -1.;
const xAOD::TrackParticle* trackTmp = nullptr;
float sumPtTmp;
for (const auto& vtx : recoVertices.stdcont()) {
if (vtx) {
sumPtTmp = 0.;
for (size_t i = 0; i < vtx->nTrackParticles(); i++) {
trackTmp = vtx->trackParticle(i);
if (trackTmp) {
sumPtTmp += std::pow(trackTmp->pt(), 2);
}
}
if (sumPtTmp > sumPtMax) {
sumPtMax = sumPtTmp;
recoHSVertex = vtx;
}
}
}
return recoHSVertex;
}
template<typename U, typename V>
float InDetPerfPlot_VertexTruthMatching::getRadialDiff2(const U* vtx1, const V* vtx2) const {
return (std::pow((vtx1->x() - vtx2->x()), 2) + std::pow((vtx1->y() - vtx2->y()), 2) + std::pow((vtx1->z() - vtx2->z()), 2));
}
float InDetPerfPlot_VertexTruthMatching::getLocalPUDensity(const xAOD::TruthVertex* vtxOfInterest, const std::vector<const xAOD::TruthVertex*>& truthHSVertices, const std::vector<const xAOD::TruthVertex*>& truthPUVertices, const float radialWindow) const {
float radialWindow2 = std::pow(radialWindow, 2);
int nTracksInWindow = 0;
float localPUDensity;
float radialDiff2;
for (const auto& vtx : truthHSVertices) {
if (vtx != vtxOfInterest) {
radialDiff2 = getRadialDiff2(vtxOfInterest, vtx);
if (radialDiff2 < radialWindow2) {
nTracksInWindow += 1;
}
}
}
for (const auto& vtx : truthPUVertices) {
if (vtx != vtxOfInterest) {
radialDiff2 = getRadialDiff2(vtxOfInterest, vtx);
if (radialDiff2 < radialWindow2) {
nTracksInWindow += 1;
}
}
}
localPUDensity = (float)(nTracksInWindow) / (2 * radialWindow);
return localPUDensity;
}
float InDetPerfPlot_VertexTruthMatching::getRecoLongitudinalReso(const xAOD::Vertex* recoVtx) const {
return std::sqrt(recoVtx->covariancePosition()(2, 2));
}
float InDetPerfPlot_VertexTruthMatching::getRecoTransverseReso(const xAOD::Vertex* recoVtx) const {
float x = recoVtx->x();
float y = recoVtx->y();
float xErr2 = recoVtx->covariancePosition()(0, 0);
float yErr2 = recoVtx->covariancePosition()(1, 1);
float xyCov = recoVtx->covariancePosition()(0, 1);
float r2 = std::pow(x, 2) + std::pow(y, 2);
return std::sqrt(std::pow(x, 2) / r2 * xErr2 + std::pow(y, 2) / r2 * yErr2 + 2 * x * y / r2 * xyCov);
}
void
InDetPerfPlot_VertexTruthMatching::fill(const xAOD::Vertex& vertex) {
// Copied from Graham:
void InDetPerfPlot_VertexTruthMatching::fillResoHist(TH1* resoHist, const TH2* resoHist2D) {
TH1* projHist = nullptr;
int safety_counter;
TFitResultPtr fitResult;
double mean;
double rms;
double itr_rms = -1.;
double itr_rms_err;
for (int i = 1; i < resoHist2D->GetNbinsX() + 1; i++) {
projHist = resoHist2D->ProjectionY("projectionY", i, i);
if (projHist->GetEntries() == 0.) {
resoHist->SetBinContent(i, 0.);
resoHist->SetBinError(i, 0.);
continue;
}
safety_counter = 0;
fitResult = projHist->Fit("gaus", "QS0");
if (!fitResult.Get()) {
// Is it necessary to also require fitResult->Status() % 1000 == 0 for a successful fit?
// --> fitStatus = migradResult + 10 * minosResult + 100 * hesseResult + 1000 * improveResult
resoHist->SetBinContent(i, 0.);
resoHist->SetBinError(i, 0.);
continue;
}
mean = fitResult->Parameter(1);
rms = fitResult->Parameter(2);
while (true) {
projHist->SetAxisRange(mean - 3 * rms, mean + 3 * rms, "X");
fitResult = projHist->Fit("gaus", "QS0");
if (!fitResult.Get()) {
itr_rms = 0.;
itr_rms_err = 0.;
break;
}
itr_rms = fitResult->Parameter(2);
itr_rms_err = fitResult->ParError(2);
if ((fabs(itr_rms - rms) < 0.0001) || (safety_counter == 5)) {
break;
}
safety_counter++;
mean = fitResult->Parameter(1);
rms = itr_rms;
continue;
}
resoHist->SetBinContent(i, itr_rms);
resoHist->SetBinError(i, itr_rms_err);
// fill vertex truth match type
}
}
const xAOD::TruthVertex* InDetPerfPlot_VertexTruthMatching::getTruthVertex(const xAOD::Vertex* recoVtx) const {
const xAOD::TruthVertex* truthVtx = nullptr;
if (recoVtx) {
const static xAOD::Vertex::Decorator<std::vector<InDetVertexTruthMatchUtils::VertexTruthMatchInfo>> truthMatchingInfos("TruthEventMatchingInfos");
const std::vector<InDetVertexTruthMatchUtils::VertexTruthMatchInfo>& truthInfos = truthMatchingInfos(*recoVtx);
if (!truthInfos.empty()) {
const InDetVertexTruthMatchUtils::VertexTruthMatchInfo& truthInfo = truthInfos.at(0);
const ElementLink<xAOD::TruthEventBaseContainer> truthEventLink = std::get<0>(truthInfo);
const xAOD::TruthEvent* truthEvent = nullptr;
if (truthEventLink.isValid()) {
truthEvent = static_cast<const xAOD::TruthEvent*>(*truthEventLink);
if (truthEvent) {
truthVtx = truthEvent->truthVertex(0);
}
}
}
else {
ATH_MSG_WARNING("TruthEventMatchingInfos DECORATOR yields empty vector -- returning nullptr!");
}
}
return truthVtx;
}
void InDetPerfPlot_VertexTruthMatching::fill(const xAOD::Vertex& vertex) {
// Get the match type info for each vertex:
const static xAOD::Vertex::Decorator<InDetVertexTruthMatchUtils::VertexMatchType> recoVtxMatchTypeInfo("VertexMatchType");
InDetVertexTruthMatchUtils::VertexMatchType matchType;
if (recoVtxMatchTypeInfo.isAvailable(vertex)) {
try {
ATH_MSG_DEBUG("VERTEX DECORATOR ======= "<< recoVtxMatchTypeInfo(vertex)<< ", with nTRACKS === " << vertex.nTrackParticles() << ", vertex index = " << vertex.index()<< " AT (" << vertex.x() << ", " << vertex.y() << ", " << vertex.z() << ")");
fillHisto(m_vx_type_truth, recoVtxMatchTypeInfo(vertex));
matchType = recoVtxMatchTypeInfo(vertex);
ATH_MSG_DEBUG("VERTEX DECORATOR ======= " << matchType << ", with nTRACKS === " << vertex.nTrackParticles() << ", vertex index = " << vertex.index() << " AT (x, y, z) = (" << vertex.x() << ", " << vertex.y() << ", " << vertex.z() << ")");
fillHisto(m_vx_type_truth, matchType);
}
catch (SG::ExcBadAuxVar &) {
ATH_MSG_WARNING("VertexMatchType DECORATOR seems to be available, but may be broken ===========");
......@@ -65,64 +238,149 @@ InDetPerfPlot_VertexTruthMatching::fill(const xAOD::Vertex& vertex) {
}
else {
ATH_MSG_WARNING("VertexMatchType DECORATOR is NOT available ===========");
}
}
} //end InDetPerfPlot_VertexTruthMatching::fill(const xAOD::Vertex& vertex)
} // void InDetPerfPlot_VertexTruthMatching::fill(const xAOD::Vertex& vertex) {
void
InDetPerfPlot_VertexTruthMatching::fill(const xAOD::VertexContainer& vertexContainer, const std::vector<const xAOD::TruthVertex*>& truthVertices) {
void InDetPerfPlot_VertexTruthMatching::fill(const xAOD::VertexContainer& vertexContainer, const std::vector<const xAOD::TruthVertex*>& truthHSVertices, const std::vector<const xAOD::TruthVertex*>& truthPUVertices) {
// Fill our histogram
// We just need to add the size of the vertex container as an argument, are there other things we might need to fill in this kind of function?
// Inclusive:
fillHisto(m_vx_nReco_vs_nTruth_inclusive, (float)truthVertices.size(), (float)vertexContainer.size());
// Let's also plot the vertices by vertex match type:
const static xAOD::Vertex::Decorator<InDetVertexTruthMatchUtils::VertexMatchType> recoVtxMatchTypeInfo("VertexMatchType");
std::map<InDetVertexTruthMatchUtils::VertexMatchType, int> breakdown = {};
breakdown[InDetVertexTruthMatchUtils::VertexMatchType::MATCHED] = 0;
breakdown[InDetVertexTruthMatchUtils::VertexMatchType::MERGED] = 0;
breakdown[InDetVertexTruthMatchUtils::VertexMatchType::SPLIT] = 0;
breakdown[InDetVertexTruthMatchUtils::VertexMatchType::FAKE] = 0;
breakdown[InDetVertexTruthMatchUtils::VertexMatchType::DUMMY] = 0;
for (const auto& vertex : vertexContainer.stdcont()) {
if (recoVtxMatchTypeInfo.isAvailable(*vertex)) {
breakdown[recoVtxMatchTypeInfo(*vertex)] += 1;
}
}
fillHisto(m_vx_nReco_vs_nTruth_matched, (float)truthVertices.size(), (float)breakdown[InDetVertexTruthMatchUtils::VertexMatchType::MATCHED]);
fillHisto(m_vx_nReco_vs_nTruth_merged, (float)truthVertices.size(), (float)breakdown[InDetVertexTruthMatchUtils::VertexMatchType::MERGED]);
fillHisto(m_vx_nReco_vs_nTruth_split, (float)truthVertices.size(), (float)breakdown[InDetVertexTruthMatchUtils::VertexMatchType::SPLIT]);
fillHisto(m_vx_nReco_vs_nTruth_fake, (float)truthVertices.size(), (float)breakdown[InDetVertexTruthMatchUtils::VertexMatchType::FAKE]);
fillHisto(m_vx_nReco_vs_nTruth_dummy, (float)truthVertices.size(), (float)breakdown[InDetVertexTruthMatchUtils::VertexMatchType::DUMMY]);
// And by hardscatter type:
InDetVertexTruthMatchUtils::HardScatterType hsType = InDetVertexTruthMatchUtils::classifyHardScatter(vertexContainer);
fillHisto(m_vx_hs_classification, hsType);
switch(hsType) {
case InDetVertexTruthMatchUtils::HardScatterType::CLEAN: {
fillHisto(m_vx_nReco_vs_nTruth_clean, (float)truthVertices.size(), (float)vertexContainer.size());
break;
}
case InDetVertexTruthMatchUtils::HardScatterType::LOWPU: {
fillHisto(m_vx_nReco_vs_nTruth_lowpu, (float)truthVertices.size(), (float)vertexContainer.size());
break;
}
case InDetVertexTruthMatchUtils::HardScatterType::HIGHPU: {
fillHisto(m_vx_nReco_vs_nTruth_highpu, (float)truthVertices.size(), (float)vertexContainer.size());
break;
if (m_iDetailLevel >= 200) {
// Fill our histograms
// Inclusive:
int nTruthVertices = (int)(truthHSVertices.size() + truthPUVertices.size());
int nRecoVertices = (int)vertexContainer.size();
fillHisto(m_vx_nReco_vs_nTruth_inclusive, nTruthVertices, nRecoVertices);
// Let's also plot the vertices by vertex match type:
const static xAOD::Vertex::Decorator<InDetVertexTruthMatchUtils::VertexMatchType> recoVtxMatchTypeInfo("VertexMatchType");
std::map<InDetVertexTruthMatchUtils::VertexMatchType, int> breakdown = {};
breakdown[InDetVertexTruthMatchUtils::VertexMatchType::MATCHED] = 0;
breakdown[InDetVertexTruthMatchUtils::VertexMatchType::MERGED] = 0;
breakdown[InDetVertexTruthMatchUtils::VertexMatchType::SPLIT] = 0;
breakdown[InDetVertexTruthMatchUtils::VertexMatchType::FAKE] = 0;
breakdown[InDetVertexTruthMatchUtils::VertexMatchType::DUMMY] = 0;
const xAOD::TruthVertex* truthVtx = nullptr;
float localPUDensity;
// Best reco HS vertex identified via sumpt2
const xAOD::Vertex* bestRecoHSVtx_sumpt2 = getHSRecoVertexSumPt2(vertexContainer); // Could potentially use the first vertex in the container if sumpt2-ordered
// Best reco HS vertex identified via truth HS weights
const xAOD::Vertex* bestRecoHSVtx_truth = InDetVertexTruthMatchUtils::bestHardScatterMatch(vertexContainer);
// Did we correctly select the best reco HS vertex using sumpt2?
truthVtx = getTruthVertex(bestRecoHSVtx_sumpt2);
localPUDensity = getLocalPUDensity(truthVtx, truthHSVertices, truthPUVertices);
fillHisto(m_vx_hs_sel_eff, localPUDensity, (bestRecoHSVtx_sumpt2 == bestRecoHSVtx_truth));
// Did we successfully reconstruct our truth HS vertex?
bool truthHSVtxRecoed = false;
float minTruthRecoRadialDiff2 = std::pow(m_cutMinTruthRecoRadialDiff, 2);
float truthRecoRadialDiff2 = -1.;
const xAOD::TruthVertex* truthHSVtx = nullptr;
// Check that we have *exactly* 1 truth HS vertex
if (truthHSVertices.size() != 0) {
if (truthHSVertices.size() != 1) {
ATH_MSG_WARNING("Size of truth HS vertex vector is >1 -- only using the first one in the vector.");
}
truthHSVtx = truthHSVertices.at(0);
// If the radial difference between the truth-pkg-selected best reco HS vertex and the truth HS vertex is
// less than some cut (e.g., 0.1 mm), then we say the truth HS vertex is reconstructed
truthRecoRadialDiff2 = getRadialDiff2(bestRecoHSVtx_truth, truthHSVtx);
if (truthRecoRadialDiff2 < minTruthRecoRadialDiff2) {
truthHSVtxRecoed = true;
minTruthRecoRadialDiff2 = truthRecoRadialDiff2;
}
}
case InDetVertexTruthMatchUtils::HardScatterType::HSSPLIT: {
fillHisto(m_vx_nReco_vs_nTruth_hssplit, (float)truthVertices.size(), (float)vertexContainer.size());
break;
else {
ATH_MSG_WARNING("Size of truth HS vertex vector is 0 -- assuming truth HS vertex to NOT be reconstructed.");
}
case InDetVertexTruthMatchUtils::HardScatterType::NONE: {
fillHisto(m_vx_nReco_vs_nTruth_none, (float)truthVertices.size(), (float)vertexContainer.size());
break;
// Iterate over vertices:
InDetVertexTruthMatchUtils::VertexMatchType matchType;
for (const auto& vertex : vertexContainer.stdcont()) {
// Skip dummy vertex (last one in the container)
if (vertex->vertexType() == xAOD::VxType::NoVtx) {
continue;
}
fill(*vertex);
matchType = recoVtxMatchTypeInfo(*vertex);
breakdown[matchType] += 1;
// If we have reconstructed the truth HS vertex but we have a different reco vertex closer to the truth HS vertex
// than the best one identified by the truth pkg, we say we have NOT successfully reconstructed the truth HS vertex
if (truthHSVtxRecoed && (vertex != bestRecoHSVtx_truth)) {
truthRecoRadialDiff2 = getRadialDiff2(vertex, truthHSVtx);
if (truthRecoRadialDiff2 < minTruthRecoRadialDiff2) {
truthHSVtxRecoed = false;
minTruthRecoRadialDiff2 = truthRecoRadialDiff2;
}
}
} // end loop over vertices
// Now fill plots relating to the reconstruction of our truth HS vertex (efficiency and resolutions)
if (truthHSVertices.size() != 0) {
localPUDensity = getLocalPUDensity(truthHSVtx, truthHSVertices, truthPUVertices);
if (truthHSVtxRecoed) {
fillHisto(m_vx_hs_reco_eff, localPUDensity, 1);
fillHisto(m_vx_hs_reco_long_reso, localPUDensity, getRecoLongitudinalReso(bestRecoHSVtx_truth));
fillHisto(m_vx_hs_reco_trans_reso, localPUDensity, getRecoTransverseReso(bestRecoHSVtx_truth));
fillHisto(m_vx_hs_truth_long_reso_vs_PU, localPUDensity, truthHSVtx->z() - bestRecoHSVtx_truth->z());
fillHisto(m_vx_hs_truth_trans_reso_vs_PU, localPUDensity, std::sqrt(std::pow(truthHSVtx->x() - bestRecoHSVtx_truth->x(), 2) + std::pow(truthHSVtx->y() - bestRecoHSVtx_truth->y(), 2)));
}
else {
fillHisto(m_vx_hs_reco_eff, localPUDensity, 0);
}
}
default: {
break;
fillHisto(m_vx_nReco_vs_nTruth_matched, nTruthVertices, breakdown[InDetVertexTruthMatchUtils::VertexMatchType::MATCHED]);
fillHisto(m_vx_nReco_vs_nTruth_merged, nTruthVertices, breakdown[InDetVertexTruthMatchUtils::VertexMatchType::MERGED]);
fillHisto(m_vx_nReco_vs_nTruth_split, nTruthVertices, breakdown[InDetVertexTruthMatchUtils::VertexMatchType::SPLIT]);
fillHisto(m_vx_nReco_vs_nTruth_fake, nTruthVertices, breakdown[InDetVertexTruthMatchUtils::VertexMatchType::FAKE]);
fillHisto(m_vx_nReco_vs_nTruth_dummy, nTruthVertices, breakdown[InDetVertexTruthMatchUtils::VertexMatchType::DUMMY]);
// And by hardscatter type:
InDetVertexTruthMatchUtils::HardScatterType hsType = InDetVertexTruthMatchUtils::classifyHardScatter(vertexContainer);
fillHisto(m_vx_hs_classification, hsType);
switch (hsType) {
case InDetVertexTruthMatchUtils::HardScatterType::CLEAN: {
fillHisto(m_vx_nReco_vs_nTruth_clean, nTruthVertices, nRecoVertices);
break;
}
case InDetVertexTruthMatchUtils::HardScatterType::LOWPU: {
fillHisto(m_vx_nReco_vs_nTruth_lowpu, nTruthVertices, nRecoVertices);
break;
}
case InDetVertexTruthMatchUtils::HardScatterType::HIGHPU: {
fillHisto(m_vx_nReco_vs_nTruth_highpu, nTruthVertices, nRecoVertices);
break;
}
case InDetVertexTruthMatchUtils::HardScatterType::HSSPLIT: {
fillHisto(m_vx_nReco_vs_nTruth_hssplit, nTruthVertices, nRecoVertices);
break;
}
case InDetVertexTruthMatchUtils::HardScatterType::NONE: {
fillHisto(m_vx_nReco_vs_nTruth_none, nTruthVertices, nRecoVertices);
break;
}
default: {
break;
}
}
}
}//end InDetPerfPlot_VertexTruthMatching::fill(const xAOD::VertexContainer& vertexContainer, const xAOD::TruthEventContainer& truthHSEvents)
} // end InDetPerfPlot_VertexTruthMatching::fill(const xAOD::VertexContainer& vertexContainer, const std::vector<const xAOD::TruthVertex*>& truthHSVertices, const std::vector<const xAOD::TruthVertex*>& truthPUVertices)
void InDetPerfPlot_VertexTruthMatching::finalizePlots() {
if (m_iDetailLevel >= 200) {
fillResoHist(m_vx_hs_truth_long_reso, m_vx_hs_truth_long_reso_vs_PU);
fillResoHist(m_vx_hs_truth_trans_reso, m_vx_hs_truth_trans_reso_vs_PU);
}
} // end InDetPerfPlot_VertexTruthMatching::finalizePlots()
......@@ -14,6 +14,8 @@
#include "InDetPhysValMonitoringUtilities.h"
// Tracking includes:
#include "xAODTracking/TrackParticleContainer.h"
#include "xAODTracking/VertexContainer.h"
#include "xAODTracking/VertexFwd.h"
// xAOD truth object includes:
......@@ -28,35 +30,63 @@
// std includes
#include <string>
class TH1;
// root includes
#include "TF1.h"
#include "TH1.h"
#include "TH2.h"
#include "TFitResult.h"
#include "TFitResultPtr.h"
///class holding plots for truth matched vertices
class InDetPerfPlot_VertexTruthMatching: public InDetPlotBase {
public:
InDetPerfPlot_VertexTruthMatching(InDetPlotBase* pParent, const std::string& dirName);
void fill(const xAOD::Vertex& vertex);
void fill(const xAOD::VertexContainer& vertexContainer, const std::vector<const xAOD::TruthVertex*>& truthVertices);
InDetPerfPlot_VertexTruthMatching(InDetPlotBase* pParent, const std::string& dirName, const int iDetailLevel = 10);
void fill(const xAOD::Vertex& vertex);
void fill(const xAOD::VertexContainer& vertexContainer, const std::vector<const xAOD::TruthVertex*>& truthHSVertices, const std::vector<const xAOD::TruthVertex*>& truthPUVertices);
private:
///truth type
TH1* m_vx_type_truth;
///hardscatter classification
TH1* m_vx_hs_classification;
///vertex reco efficiency
TProfile* m_vx_nReco_vs_nTruth_inclusive;
TProfile* m_vx_nReco_vs_nTruth_matched;
TProfile* m_vx_nReco_vs_nTruth_merged;
TProfile* m_vx_nReco_vs_nTruth_split;
TProfile* m_vx_nReco_vs_nTruth_fake;
TProfile* m_vx_nReco_vs_nTruth_dummy;
TProfile* m_vx_nReco_vs_nTruth_clean;
TProfile* m_vx_nReco_vs_nTruth_lowpu;
TProfile* m_vx_nReco_vs_nTruth_highpu;
TProfile* m_vx_nReco_vs_nTruth_hssplit;
TProfile* m_vx_nReco_vs_nTruth_none;
///@}
// plot base has no default implementation of this; we use it to book the histos
void initializePlots();
int m_iDetailLevel;
float m_cutMinTruthRecoRadialDiff = 0.1;
///truth type
TH1* m_vx_type_truth;
///hardscatter classification
TH1* m_vx_hs_classification;
///vertex reco efficiency
TProfile* m_vx_nReco_vs_nTruth_inclusive;
TProfile* m_vx_nReco_vs_nTruth_matched;
TProfile* m_vx_nReco_vs_nTruth_merged;
TProfile* m_vx_nReco_vs_nTruth_split;
TProfile* m_vx_nReco_vs_nTruth_fake;
TProfile* m_vx_nReco_vs_nTruth_dummy;
TProfile* m_vx_nReco_vs_nTruth_clean;
TProfile* m_vx_nReco_vs_nTruth_lowpu;
TProfile* m_vx_nReco_vs_nTruth_highpu;
TProfile* m_vx_nReco_vs_nTruth_hssplit;
TProfile* m_vx_nReco_vs_nTruth_none;
// HS vertex reconstruction efficiency vs PU
TEfficiency* m_vx_hs_reco_eff;
// HS vertex selection efficiency vs PU
TEfficiency* m_vx_hs_sel_eff;
// For reco (covariance) resolutions:
TProfile* m_vx_hs_reco_long_reso;
TProfile* m_vx_hs_reco_trans_reso;
// For reco-truth resolutions:
TH2* m_vx_hs_truth_long_reso_vs_PU;
TH2* m_vx_hs_truth_trans_reso_vs_PU;
TH1* m_vx_hs_truth_long_reso;
TH1* m_vx_hs_truth_trans_reso;
///@}
private:
// plot base has no default implementation of this; we use it to book the histos
void initializePlots();
const xAOD::Vertex* getHSRecoVertexSumPt2(const xAOD::VertexContainer& recoVertices) const;
template<typename U, typename V>
float getRadialDiff2(const U* vtx1, const V* vtx2) const;
float getLocalPUDensity(const xAOD::TruthVertex* vtxOfInterest, const std::vector<const xAOD::TruthVertex*>& truthHSVertices, const std::vector<const xAOD::TruthVertex*>& truthPUVertices, const float radialWindow = 2.0) const;
float getRecoLongitudinalReso(const xAOD::Vertex* recoVtx) const;
float getRecoTransverseReso(const xAOD::Vertex* recoVtx) const;
const xAOD::TruthVertex* getTruthVertex(const xAOD::Vertex* recoVtx) const;
void fillResoHist(TH1* resoHist, const TH2* resoHist2D);
void finalizePlots();
};
#endif
......@@ -160,15 +160,15 @@ InDetPhysValMonitoringTool::initialize() {
if (m_truthSelectionTool.get() ) {
m_truthCutFlow = CutFlow(m_truthSelectionTool->nCuts());
}
m_monPlots = std::make_unique<InDetRttPlots> (nullptr, m_dirName + m_folder);
m_monPlots = std::make_unique<InDetRttPlots> (nullptr, m_dirName + m_folder, m_detailLevel); // m_detailLevel := DEBUG, enable expert histograms
ATH_CHECK( m_trkParticleName.initialize() );
ATH_CHECK( m_truthParticleName.initialize(m_pileupSwitch == "All" and not m_truthParticleName.key().empty() ) );
ATH_CHECK( m_truthParticleName.initialize( (m_pileupSwitch == "HardScatter" or m_pileupSwitch == "All") and not m_truthParticleName.key().empty() ) );
ATH_CHECK( m_vertexContainerName.initialize() );
ATH_CHECK( m_truthVertexContainerName.initialize( not m_truthVertexContainerName.key().empty() ) );
ATH_CHECK( m_eventInfoContainerName.initialize() );
ATH_CHECK( m_truthEventName.initialize( m_pileupSwitch == "HardScatter" and not m_truthEventName.key().empty()) );
ATH_CHECK( m_truthEventName.initialize( (m_pileupSwitch == "HardScatter" or m_pileupSwitch == "All") and not m_truthEventName.key().empty() ) );
ATH_CHECK( m_truthPileUpEventName.initialize( (m_pileupSwitch == "PileUp" or m_pileupSwitch == "All") and not m_truthPileUpEventName.key().empty() ) );
ATH_CHECK( m_jetContainerName.initialize( not m_jetContainerName.key().empty()) );
......@@ -202,9 +202,8 @@ InDetPhysValMonitoringTool::fillHistograms() {
SG::ReadHandle<xAOD::EventInfo> pie = SG::ReadHandle<xAOD::EventInfo>(m_eventInfoContainerName);
std::vector<const xAOD::TruthParticle*> truthParticlesVec = getTruthParticles();
std::vector<const xAOD::TruthVertex*> truthVertices = getTruthVertices();
IDPVM::CachedGetAssocTruth getAsTruth; // only cache one way, track->truth, not truth->tracks
if (not tracks.isValid()) {
return StatusCode::FAILURE;
}
......@@ -229,12 +228,17 @@ InDetPhysValMonitoringTool::fillHistograms() {
//Filling plots for all reconstructed vertices and the hard-scatter
ATH_MSG_DEBUG("Filling vertices info monitoring plots");
//Decorate vertices
// Fill vectors of truth HS and PU vertices
std::pair<std::vector<const xAOD::TruthVertex*>, std::vector<const xAOD::TruthVertex*>> truthVertices = getTruthVertices();
std::vector<const xAOD::TruthVertex*> truthHSVertices = truthVertices.first;
std::vector<const xAOD::TruthVertex*> truthPUVertices = truthVertices.second;
// Decorate vertices
if (m_useVertexTruthMatchTool && m_vtxValidTool) {
ATH_CHECK(m_vtxValidTool->matchVertices(*vertices));
ATH_MSG_DEBUG("Hard scatter classification type: " << InDetVertexTruthMatchUtils::classifyHardScatter(*vertices) << ", vertex container size = " << vertices->size());
}
m_monPlots->fill(*vertices, truthVertices);
m_monPlots->fill(*vertices, truthHSVertices, truthPUVertices);
ATH_MSG_DEBUG("Filling vertex/event info monitoring plots");
//Filling vertexing plots for the reconstructed hard-scatter as a function of mu
......@@ -452,7 +456,6 @@ InDetPhysValMonitoringTool::fillHistograms() {
StatusCode
InDetPhysValMonitoringTool::bookHistograms() {
ATH_MSG_INFO("Booking hists " << name() << "with detailed level: " << m_detailLevel);
m_monPlots->setDetailLevel(m_detailLevel); // DEBUG, enable expert histograms
m_monPlots->initialize();
std::vector<HistData> hists = m_monPlots->retrieveBookedHistograms();
for (auto hist : hists) {
......@@ -494,7 +497,7 @@ InDetPhysValMonitoringTool::procHistograms() {
}
const std::vector<const xAOD::TruthParticle*>
InDetPhysValMonitoringTool::getTruthParticles() {
InDetPhysValMonitoringTool::getTruthParticles() const {
// truthParticles.clear();
std::vector<const xAOD::TruthParticle*> tempVec {};
if (m_pileupSwitch == "All") {
......@@ -553,11 +556,13 @@ InDetPhysValMonitoringTool::getTruthParticles() {
return tempVec;
}
const std::vector<const xAOD::TruthVertex*>
InDetPhysValMonitoringTool::getTruthVertices() {
std::pair<const std::vector<const xAOD::TruthVertex*>, const std::vector<const xAOD::TruthVertex*>>
InDetPhysValMonitoringTool::getTruthVertices() const {
std::vector<const xAOD::TruthVertex*> truthVertices = {};
truthVertices.reserve(100);
std::vector<const xAOD::TruthVertex*> truthHSVertices = {};
truthHSVertices.reserve(5);
std::vector<const xAOD::TruthVertex*> truthPUVertices = {};
truthPUVertices.reserve(100);
const xAOD::TruthVertex* truthVtx = nullptr;
bool doHS = false;
......@@ -577,14 +582,14 @@ InDetPhysValMonitoringTool::getTruthVertices() {
}
if (doHS) {
if (!m_truthEventName.key().empty()) {
if (not m_truthEventName.key().empty()) {
ATH_MSG_VERBOSE("Getting TruthEvents container.");
SG::ReadHandle<xAOD::TruthEventContainer> truthEventContainer(m_truthEventName);
if (truthEventContainer.isValid()) {
for (const auto& evt : *truthEventContainer) {
truthVtx = evt->truthVertex(0);
if (truthVtx) {
truthVertices.push_back(truthVtx);
truthHSVertices.push_back(truthVtx);
}
}
}
......@@ -595,14 +600,14 @@ InDetPhysValMonitoringTool::getTruthVertices() {
}
if (doPU) {
if (!m_truthPileUpEventName.key().empty()) {
if (not m_truthPileUpEventName.key().empty()) {
ATH_MSG_VERBOSE("Getting TruthEvents container.");
SG::ReadHandle<xAOD::TruthPileupEventContainer> truthPileupEventContainer(m_truthPileUpEventName);
if (truthPileupEventContainer.isValid()) {
for (const auto& evt : *truthPileupEventContainer) {
truthVtx = evt->truthVertex(0);
if (truthVtx) {
truthVertices.push_back(truthVtx);
truthPUVertices.push_back(truthVtx);
}
}
}
......@@ -612,7 +617,7 @@ InDetPhysValMonitoringTool::getTruthVertices() {
}
}
return truthVertices;
return std::make_pair<const std::vector<const xAOD::TruthVertex*>, const std::vector<const xAOD::TruthVertex*>>((const std::vector<const xAOD::TruthVertex*>)truthHSVertices, (const std::vector<const xAOD::TruthVertex*>)truthPUVertices);
}
......
......@@ -15,7 +15,7 @@
#include <cmath> // std::isnan()
#include <limits>
InDetRttPlots::InDetRttPlots(InDetPlotBase* pParent, const std::string& sDir) : InDetPlotBase(pParent, sDir),
InDetRttPlots::InDetRttPlots(InDetPlotBase* pParent, const std::string& sDir, const int iDetailLevel) : InDetPlotBase(pParent, sDir),
m_trackParameters(this, "Tracks/Selected/Parameters"),
m_nTracks(this, "Tracks/Tracks"),
m_hitResidualPlot(this, "Tracks/Hits/Residuals"),
......@@ -32,13 +32,14 @@ InDetRttPlots::InDetRttPlots(InDetPlotBase* pParent, const std::string& sDir) :
m_resolutionPlotSecd(nullptr),
m_doTrackInJetPlots(true) //FIX CONFIGURATION
{
this->m_iDetailLevel = iDetailLevel;
m_trackParticleTruthProbKey = "truthMatchProbability";
m_truthProbLowThreshold = 0.5;
if(m_iDetailLevel >= 200){
m_resolutionPlotSecd = std::unique_ptr<InDetPerfPlot_Resolution>(new InDetPerfPlot_Resolution(this, "Tracks/Matched/Resolutions/Secondary"));
m_hitsMatchedTracksPlots = std::unique_ptr<InDetPerfPlot_Hits>(new InDetPerfPlot_Hits(this, "Tracks/Matched/HitsOnTracks"));
m_vertexTruthMatchingPlots = std::unique_ptr<InDetPerfPlot_VertexTruthMatching>(new InDetPerfPlot_VertexTruthMatching(this, "Vertices/AllPrimaryVertices"));
m_vertexTruthMatchingPlots = std::unique_ptr<InDetPerfPlot_VertexTruthMatching>(new InDetPerfPlot_VertexTruthMatching(this, "Vertices/AllPrimaryVertices", m_iDetailLevel));
}
//A lot of Jets... do we need these at all???
......@@ -125,7 +126,7 @@ InDetRttPlots::fillFakeRate(const xAOD::TrackParticle& track, const bool isFake,
//Fill Vertexing Plots
//
void
InDetRttPlots::fill(const xAOD::VertexContainer& vertexContainer, const std::vector<const xAOD::TruthVertex*>& truthVertices) {
InDetRttPlots::fill(const xAOD::VertexContainer& vertexContainer, const std::vector<const xAOD::TruthVertex*>& truthHSVertices, const std::vector<const xAOD::TruthVertex*>& truthPUVertices) {
// fill vertex container general properties
// m_verticesVsMuPlots.fill(vertexContainer); //if ever needed
// fill vertex-specific properties, for all vertices and for hard-scattering vertex
......@@ -135,9 +136,6 @@ InDetRttPlots::fill(const xAOD::VertexContainer& vertexContainer, const std::vec
continue; // skip dummy vertex
}
m_vertexPlots.fill(*vtx);
if(m_iDetailLevel >= 200){
m_vertexTruthMatchingPlots->fill(*vtx);
}
ATH_MSG_DEBUG("IN InDetRttPlots::fill, filling for all vertices");
if (vtx->vertexType() == xAOD::VxType::PriVtx) {
m_hardScatterVertexPlots.fill(*vtx);
......@@ -146,7 +144,7 @@ InDetRttPlots::fill(const xAOD::VertexContainer& vertexContainer, const std::vec
}
}
if(m_iDetailLevel >= 200){
m_vertexTruthMatchingPlots->fill(vertexContainer, truthVertices);
m_vertexTruthMatchingPlots->fill(vertexContainer, truthHSVertices, truthPUVertices);
}
}
......
......@@ -47,7 +47,7 @@
///class holding all plots for Inner Detector RTT Validation and implementing fill methods
class InDetRttPlots: public InDetPlotBase {
public:
InDetRttPlots(InDetPlotBase* pParent, const std::string& dirName);
InDetRttPlots(InDetPlotBase* pParent, const std::string& dirName, const int iDetailLevel = 10);
///fill for things needing truth and track only
void fill(const xAOD::TrackParticle& particle, const xAOD::TruthParticle& truthParticle);
......@@ -61,7 +61,7 @@ public:
///fill for things needing all truth - not just the ones from the reco tracks
///fill reco-vertex related plots
void fill(const xAOD::VertexContainer& vertexContainer, const std::vector<const xAOD::TruthVertex*>& truthVertices);
void fill(const xAOD::VertexContainer& vertexContainer, const std::vector<const xAOD::TruthVertex*>& truthHSVertices, const std::vector<const xAOD::TruthVertex*>& truthPUVertices);
///fill reco-vertex related plots that need EventInfo
void fill(const xAOD::VertexContainer& vertexContainer, unsigned int nPU);
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment