10 #include "GaudiKernel/SystemOfUnits.h"
50 float result(std::numeric_limits<float>::quiet_NaN());
53 if (truthMatchProbabilityAcc.isAvailable(trackParticle)) {
54 result = truthMatchProbabilityAcc(trackParticle);
61 safelyGetEta(
const T& pTrk,
const float safePtThreshold = 0.1) {
62 return (pTrk->pt() > safePtThreshold) ? (pTrk->eta()) : std::nan(
"");
69 return not ((
value < minVal)or(
value > maxVal));
75 return not (std::abs(
value) > absoluteMax);
90 return nf ?
i :
i - 1;
95 const float x(vtx->
x()),
y(vtx->
y()),
z(vtx->
z());
96 const float vrtR2 = (
x *
x +
y *
y);
98 return inRange(
z, 500.
f) and not (vrtR2 > 1);
101 const std::vector<float> ETA_PARTITIONS = {
102 2.7, 3.5, std::numeric_limits<float>::infinity()
109 const IInterface*
parent) :
134 (
nullptr,
static_cast<std::string
>(
m_dirName) +
static_cast<std::string
>(
m_folder),
147 std::vector<std::string> required_float_track_decorations {
"d0",
"hitResiduals_residualLocX",
"d0err"};
148 std::vector<std::string> required_int_track_decorations {};
149 std::vector<std::string> required_float_truth_decorations {
"d0"};
150 std::vector<std::string> required_int_truth_decorations {};
151 std::vector<std::string> required_int_jet_decorations {
"HadronConeExclTruthLabelID"};
156 m_datfile <<
"EvtNumber,PdgID,Px,Py,Pz,E,Pt,Eta,Phi,Mass,numPU,numvtx,MatchProb"<<std::endl;
158 std::string empty_prefix;
172 return StatusCode::SUCCESS;
311 const EventContext& ctx = Gaudi::Hive::currentContext();
315 if (not trackHandle.
isValid()) {
317 return StatusCode::FAILURE;
328 ATH_MSG_WARNING(
"Shouldn't happen - EventInfo is buggy, setting mu to 0");
336 return StatusCode::SUCCESS;
340 std::vector<const xAOD::TruthParticle*> truthParticlesVec =
getTruthParticles(ctx);
349 unsigned int truthMu = 0;
352 truthMu =
static_cast<int>( truthPileupEventContainer->
size() );
357 float puEvents = truthMu>0 ? truthMu : actualMu;
361 float beamSpotWeight = 1;
364 ATH_MSG_DEBUG(
"Getting number of pu interactings per event");
372 ATH_MSG_DEBUG(
"beamSpotWeight is equal to " << beamSpotWeight);
376 if(readDecorHandle.
isAvailable()) prwWeight = readDecorHandle(*pie);
377 ATH_MSG_DEBUG(
"Applying pileup weight equal to " << prwWeight);
378 beamSpotWeight *= prwWeight;
382 ATH_MSG_DEBUG(
"Number of vertices retrieved for this event " << vertices->
size());
388 ATH_MSG_DEBUG(
"Failed to find a hard scatter vertex in this event.");
394 std::pair<std::vector<const xAOD::TruthVertex*>, std::vector<const xAOD::TruthVertex*>> truthVertices =
getTruthVertices(ctx);
395 std::vector<const xAOD::TruthVertex*> truthHSVertices = truthVertices.first;
396 std::vector<const xAOD::TruthVertex*> truthPUVertices = truthVertices.second;
403 m_monPlots->fill(*vertices, primaryvertex, truthHSVertices, truthPUVertices, actualMu, beamSpotWeight);
407 m_monPlots->fill(*vertices, truthMu, actualMu, beamSpotWeight);
424 const auto& stdVertexContainer = truthVrt->
stdcont();
426 auto findVtx = std::find_if(stdVertexContainer.rbegin(), stdVertexContainer.rend(), acceptTruthVertex);
427 truthVertex = (findVtx == stdVertexContainer.rend()) ?
nullptr : *findVtx;
431 if (not truthVertex)
ATH_MSG_INFO (
"Truth vertex did not pass cuts");
436 unsigned int nSelectedTruthTracks(0), nSelectedRecoTracks(0), nSelectedMatchedTracks(0), nAssociatedTruth(0), nMissingAssociatedTruth(0), nTruths(0);
445 std::map<const xAOD::TruthParticle*, std::vector<const xAOD::TrackParticle*>> cacheTruthMatching {};
447 std::vector<const xAOD::TrackParticle*> selectedTracks {};
448 selectedTracks.reserve(tracks->
size());
449 unsigned int nTrackTOT = 0;
450 unsigned int nTrackCentral = 0;
451 unsigned int nTrackPt1GeV = 0;
452 for (
const auto *
const thisTrack: *tracks) {
458 selectedTracks.push_back(thisTrack);
460 nSelectedRecoTracks++;
466 if (std::abs(thisTrack->eta()) < 2.5)
471 float prob = getMatchingProbability(*thisTrack);
474 if (associatedTruth) {
485 nSelectedMatchedTracks++;
486 bool truthIsFromB =
false;
490 m_monPlots->fill(*thisTrack, *associatedTruth, truthIsFromB, puEvents, beamSpotWeight);
494 const bool isAssociatedTruth = associatedTruth !=
nullptr;
497 if(!isAssociatedTruth) nMissingAssociatedTruth++;
498 m_monPlots->fillFakeRate(*thisTrack,
isFake, isAssociatedTruth, puEvents, beamSpotWeight);
504 if (isAssociatedTruth) {
509 auto cachedAssoc = cacheTruthMatching.find(associatedTruth);
511 if (cachedAssoc == cacheTruthMatching.end()) {
513 cacheTruthMatching[associatedTruth] = {thisTrack};
517 cachedAssoc->second.push_back(thisTrack);
522 m_monPlots->fillNtuple(*thisTrack, primaryvertex);
530 for (
auto& cachedAssoc: cacheTruthMatching) {
537 std::sort(cachedAssoc.second.begin(), cachedAssoc.second.end(),
543 for (
int itrack = 0; itrack < (
int) cachedAssoc.second.size(); itrack++) {
547 m_monPlots->fillNtuple(*thisTrack, *thisTruth, primaryvertex, itrack);
553 m_monPlots->fill(nTrackTOT, nTrackCentral, nTrackPt1GeV, truthMu, actualMu,
nVertices, beamSpotWeight);
558 std::lock_guard<std::mutex> lock(
m_mutex);
559 m_truthCutFlow.merge(std::move(tmp_truth_cutflow));
565 for (
int itruth = 0; itruth < (
int) truthParticlesVec.size(); itruth++) {
572 ++nSelectedTruthTracks;
573 bool isEfficient(
false);
574 float matchingProbability{};
578 auto cachedAssoc = cacheTruthMatching.find(thisTruth);
580 if (cachedAssoc == cacheTruthMatching.end()) {
582 cacheTruthMatching[thisTruth] = {};
584 m_monPlots->fillDuplicate(*thisTruth, cacheTruthMatching[thisTruth], beamSpotWeight);
591 for (
const auto& thisTrack: selectedTracks) {
593 if (associatedTruth && associatedTruth == thisTruth) {
594 float prob = getMatchingProbability(*thisTrack);
596 matchingProbability =
prob;
598 matchedTrack = thisTrack;
608 <<puEvents<<
","<<
nVertices<<
","<<matchingProbability<<std::endl;
611 ATH_MSG_ERROR(
"An error occurred: Truth particle for tracking efficiency calculation is a nullptr");
613 else if (isEfficient && !matchedTrack){
614 ATH_MSG_ERROR(
"Something went wrong - we log a truth particle as reconstructed, but the reco track is a nullptr! Bailing out... ");
617 ATH_MSG_DEBUG(
"Filling efficiency plots info monitoring plots");
618 m_monPlots->fillEfficiency(*thisTruth, matchedTrack, isEfficient, truthMu, actualMu, beamSpotWeight);
620 ATH_MSG_DEBUG(
"Filling technical efficiency plots info monitoring plots");
624 m_monPlots->fillTechnicalEfficiency(*thisTruth, isEfficient,
625 truthMu, actualMu, beamSpotWeight);
628 ATH_MSG_DEBUG(
"Cannot fill technical efficiency. Missing si hit information for truth particle.");
646 if (nSelectedRecoTracks == nMissingAssociatedTruth) {
651 ATH_MSG_DEBUG(nAssociatedTruth <<
" tracks out of " << tracks->size() <<
" had associated truth.");
670 return StatusCode::SUCCESS;
696 ATH_MSG_DEBUG(
"Selected by pileup switch decoration requested from a truth particle but not available");
703 for (
const auto& thisTruth: truthParticles) {
717 std::vector<EfficiencyData> effs =
m_monPlots->retrieveBookedEfficiencies();
718 for (
auto&
eff : effs) {
731 return StatusCode::SUCCESS;
756 return StatusCode::SUCCESS;
759 const std::vector<const xAOD::TruthParticle*>
762 std::vector<const xAOD::TruthParticle*> tempVec {};
769 if (not truthParticleContainer.
isValid()) {
772 tempVec.insert(tempVec.begin(), truthParticleContainer->
begin(), truthParticleContainer->
end());
782 const auto&
links =
event->truthParticleLinks();
783 tempVec.reserve(
event->nTruthParticles());
784 for (
const auto& link :
links) {
786 tempVec.push_back(*link);
795 if (truthPileupEventContainer.
isValid()) {
796 const unsigned int nPileup = truthPileupEventContainer->
size();
797 tempVec.reserve(nPileup * 200);
798 for (
unsigned int i(0);
i != nPileup; ++
i) {
799 const auto *eventPileup = truthPileupEventContainer->
at(
i);
801 int ntruth = eventPileup->nTruthParticles();
802 ATH_MSG_VERBOSE(
"Adding " << ntruth <<
" truth particles from TruthPileupEvents container");
803 const auto&
links = eventPileup->truthParticleLinks();
804 for (
const auto& link :
links) {
806 tempVec.push_back(*link);
821 std::pair<const std::vector<const xAOD::TruthVertex*>,
const std::vector<const xAOD::TruthVertex*>>
824 std::vector<const xAOD::TruthVertex*> truthHSVertices = {};
825 truthHSVertices.reserve(5);
826 std::vector<const xAOD::TruthVertex*> truthPUVertices = {};
827 truthPUVertices.reserve(100);
850 if (truthEventContainer.
isValid()) {
851 for (
const auto *
const evt : *truthEventContainer) {
852 truthVtx =
evt->signalProcessVertex();
854 truthHSVertices.push_back(truthVtx);
868 if (truthPileupEventContainer.
isValid()) {
869 for (
const auto *
const evt : *truthPileupEventContainer) {
874 size_t i_vtx = 0;
size_t n_vtx =
evt->nTruthVertices();
875 while(!truthVtx && i_vtx<n_vtx){
876 truthVtx =
evt->truthVertex(i_vtx);
881 truthPUVertices.push_back(truthVtx);
891 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);
902 std::vector<int>& cutFlow) {
904 if (cutFlow.empty()) {
905 names.emplace_back(
"preCut");
906 cutFlow.push_back(0);
907 for (
unsigned int i = 0;
i !=
accept.getNCuts(); ++
i) {
908 cutFlow.push_back(0);
914 bool cutPositive =
true;
915 for (
unsigned int i = 0;
i != (
accept.getNCuts() + 1); ++
i) {
931 ATH_MSG_INFO(
"Requesting cut value outside of configured eta range: clamping eta = "
932 << std::abs(truth.
eta()) <<
" to eta= " <<
absEta);
943 const std::vector<const xAOD::TruthParticle*>& truthParticles,
946 float beamSpotWeight) {
951 if (truthParticles.empty()) {
952 ATH_MSG_WARNING(
"No entries in TruthParticles truth particle container. Skipping jet plots.");
953 return StatusCode::SUCCESS;
959 return StatusCode::SUCCESS;
972 isBjet = (btagLabel(*thisJet) == 5);
980 if (not
el.isValid())
continue;
991 bool isEfficient(
false);
993 for (
const auto *thisTrack: tracks) {
999 if (associatedTruth and associatedTruth == truth) {
1000 float prob = getMatchingProbability(*thisTrack);
1008 bool truthIsFromB =
false;
1010 truthIsFromB =
true;
1012 m_monPlots->fillEfficiency(*truth, *thisJet, isEfficient, isBjet, truthIsFromB, beamSpotWeight);
1022 if (thisJet->p4().DeltaR(thisTrack->p4()) >
m_maxTrkJetDR) {
1026 float prob = getMatchingProbability(*thisTrack);
1030 const bool unlinked = (associatedTruth==
nullptr);
1032 bool truthIsFromB =
false;
1034 truthIsFromB =
true;
1036 m_monPlots->fill(*thisTrack, *thisJet, isBjet,
isFake, unlinked, truthIsFromB, beamSpotWeight);
1037 if (associatedTruth){
1038 m_monPlots->fillFakeRate(*thisTrack, *thisJet,
isFake, isBjet, truthIsFromB, beamSpotWeight);
1044 return StatusCode::SUCCESS;
1049 const float jetPt =
jet.pt();
1050 const float jetEta = std::abs(
jet.eta());