10 #include "GaudiKernel/SystemOfUnits.h"
49 float result(std::numeric_limits<float>::quiet_NaN());
52 if (truthMatchProbabilityAcc.isAvailable(trackParticle)) {
53 result = truthMatchProbabilityAcc(trackParticle);
60 safelyGetEta(
const T& pTrk,
const float safePtThreshold = 0.1) {
61 return (pTrk->pt() > safePtThreshold) ? (pTrk->eta()) : std::nan(
"");
68 return not ((
value < minVal)or(
value > maxVal));
74 return not (std::abs(
value) > absoluteMax);
89 return nf ?
i :
i - 1;
94 const float x(vtx->
x()),
y(vtx->
y()),
z(vtx->
z());
95 const float vrtR2 = (
x *
x +
y *
y);
97 return inRange(
z, 500.
f) and not (vrtR2 > 1);
100 const std::vector<float> ETA_PARTITIONS = {
101 2.7, 3.5, std::numeric_limits<float>::infinity()
108 const IInterface*
parent) :
133 (
nullptr,
static_cast<std::string
>(
m_dirName) +
static_cast<std::string
>(
m_folder),
144 std::vector<std::string> trk_decorations;
160 std::vector<std::string> required_float_track_decorations {
"d0",
"hitResiduals_residualLocX",
"d0err"};
161 std::vector<std::string> required_int_track_decorations {};
162 std::vector<std::string> required_float_truth_decorations {
"d0"};
163 std::vector<std::string> required_int_truth_decorations {};
164 std::vector<std::string> required_int_jet_decorations {
"HadronConeExclTruthLabelID"};
169 m_datfile <<
"EvtNumber,PdgID,Px,Py,Pz,E,Pt,Eta,Phi,Mass,numPU,numvtx,MatchProb"<<std::endl;
171 std::string empty_prefix;
185 return StatusCode::SUCCESS;
325 const EventContext& ctx = Gaudi::Hive::currentContext();
329 if (not trackHandle.
isValid()) {
331 return StatusCode::FAILURE;
342 ATH_MSG_WARNING(
"Shouldn't happen - EventInfo is buggy, setting mu to 0");
350 return StatusCode::SUCCESS;
354 std::vector<const xAOD::TruthParticle*> truthParticlesVec =
getTruthParticles(ctx);
363 unsigned int truthMu = 0;
366 truthMu =
static_cast<int>( truthPileupEventContainer->
size() );
371 float puEvents = truthMu>0 ? truthMu : actualMu;
374 unsigned int nVertices = 0;
375 float beamSpotWeight = 1;
378 ATH_MSG_DEBUG(
"Getting number of pu interactings per event");
384 nVertices = not vertices->
empty() ? vertices->
size() : 0;
386 ATH_MSG_DEBUG(
"beamSpotWeight is equal to " << beamSpotWeight);
390 if(readDecorHandle.
isAvailable()) prwWeight = readDecorHandle(*pie);
391 ATH_MSG_DEBUG(
"Applying pileup weight equal to " << prwWeight);
392 beamSpotWeight *= prwWeight;
396 ATH_MSG_DEBUG(
"Number of vertices retrieved for this event " << vertices->
size());
402 ATH_MSG_DEBUG(
"Failed to find a hard scatter vertex in this event.");
408 std::pair<std::vector<const xAOD::TruthVertex*>, std::vector<const xAOD::TruthVertex*>> truthVertices =
getTruthVertices(ctx);
409 std::vector<const xAOD::TruthVertex*> truthHSVertices = truthVertices.first;
410 std::vector<const xAOD::TruthVertex*> truthPUVertices = truthVertices.second;
417 m_monPlots->fill(*vertices, primaryvertex, truthHSVertices, truthPUVertices, actualMu, beamSpotWeight);
421 m_monPlots->fill(*vertices, truthMu, actualMu, beamSpotWeight);
438 const auto& stdVertexContainer = truthVrt->
stdcont();
440 auto findVtx = std::find_if(stdVertexContainer.rbegin(), stdVertexContainer.rend(), acceptTruthVertex);
441 truthVertex = (findVtx == stdVertexContainer.rend()) ?
nullptr : *findVtx;
445 if (not truthVertex)
ATH_MSG_INFO (
"Truth vertex did not pass cuts");
450 unsigned int nSelectedTruthTracks(0), nSelectedRecoTracks(0), nSelectedMatchedTracks(0), nAssociatedTruth(0), nMissingAssociatedTruth(0), nTruths(0);
459 std::map<const xAOD::TruthParticle*, std::vector<const xAOD::TrackParticle*>> cacheTruthMatching {};
461 std::vector<const xAOD::TrackParticle*> selectedTracks {};
462 selectedTracks.reserve(tracks->
size());
463 unsigned int nTrackTOT = 0;
464 unsigned int nTrackCentral = 0;
465 unsigned int nTrackPt1GeV = 0;
466 for (
const auto *
const thisTrack: *tracks) {
472 selectedTracks.push_back(thisTrack);
474 nSelectedRecoTracks++;
480 if (std::abs(thisTrack->eta()) < 2.5)
483 m_monPlots->fill(*thisTrack, puEvents, nVertices, beamSpotWeight);
485 float prob = getMatchingProbability(*thisTrack);
488 if (associatedTruth) {
499 nSelectedMatchedTracks++;
500 bool truthIsFromB =
false;
504 m_monPlots->fill(*thisTrack, *associatedTruth, truthIsFromB, puEvents, beamSpotWeight);
508 const bool isAssociatedTruth = associatedTruth !=
nullptr;
511 if(!isAssociatedTruth) nMissingAssociatedTruth++;
512 m_monPlots->fillFakeRate(*thisTrack,
isFake, isAssociatedTruth, puEvents, beamSpotWeight);
518 if (isAssociatedTruth) {
523 auto cachedAssoc = cacheTruthMatching.find(associatedTruth);
525 if (cachedAssoc == cacheTruthMatching.end()) {
527 cacheTruthMatching[associatedTruth] = {thisTrack};
531 cachedAssoc->second.push_back(thisTrack);
536 m_monPlots->fillNtuple(*thisTrack, primaryvertex);
544 for (
auto& cachedAssoc: cacheTruthMatching) {
551 std::sort(cachedAssoc.second.begin(), cachedAssoc.second.end(),
557 for (
int itrack = 0; itrack < (
int) cachedAssoc.second.size(); itrack++) {
561 m_monPlots->fillNtuple(*thisTrack, *thisTruth, primaryvertex, itrack);
567 m_monPlots->fill(nTrackTOT, nTrackCentral, nTrackPt1GeV, truthMu, actualMu, nVertices, beamSpotWeight);
573 m_truthCutFlow.merge(std::move(tmp_truth_cutflow));
579 for (
int itruth = 0; itruth < (
int) truthParticlesVec.size(); itruth++) {
586 ++nSelectedTruthTracks;
587 bool isEfficient(
false);
588 float matchingProbability{};
592 auto cachedAssoc = cacheTruthMatching.find(thisTruth);
594 if (cachedAssoc == cacheTruthMatching.end()) {
596 cacheTruthMatching[thisTruth] = {};
598 m_monPlots->fillDuplicate(*thisTruth, cacheTruthMatching[thisTruth], beamSpotWeight);
605 for (
const auto& thisTrack: selectedTracks) {
607 if (associatedTruth && associatedTruth == thisTruth) {
608 float prob = getMatchingProbability(*thisTrack);
610 matchingProbability =
prob;
612 matchedTrack = thisTrack;
622 <<puEvents<<
","<<nVertices<<
","<<matchingProbability<<std::endl;
625 ATH_MSG_ERROR(
"An error occurred: Truth particle for tracking efficiency calculation is a nullptr");
627 else if (isEfficient && !matchedTrack){
628 ATH_MSG_ERROR(
"Something went wrong - we log a truth particle as reconstructed, but the reco track is a nullptr! Bailing out... ");
631 ATH_MSG_DEBUG(
"Filling efficiency plots info monitoring plots");
632 m_monPlots->fillEfficiency(*thisTruth, matchedTrack, isEfficient, truthMu, actualMu, beamSpotWeight);
634 ATH_MSG_DEBUG(
"Filling technical efficiency plots info monitoring plots");
638 m_monPlots->fillTechnicalEfficiency(*thisTruth, isEfficient,
639 truthMu, actualMu, beamSpotWeight);
642 ATH_MSG_DEBUG(
"Cannot fill technical efficiency. Missing si hit information for truth particle.");
660 if (nSelectedRecoTracks == nMissingAssociatedTruth) {
665 ATH_MSG_DEBUG(nAssociatedTruth <<
" tracks out of " << tracks->size() <<
" had associated truth.");
684 return StatusCode::SUCCESS;
710 ATH_MSG_DEBUG(
"Selected by pileup switch decoration requested from a truth particle but not available");
717 for (
const auto& thisTruth: truthParticles) {
731 std::vector<EfficiencyData> effs =
m_monPlots->retrieveBookedEfficiencies();
732 for (
auto&
eff : effs) {
745 return StatusCode::SUCCESS;
770 return StatusCode::SUCCESS;
773 const std::vector<const xAOD::TruthParticle*>
776 std::vector<const xAOD::TruthParticle*> tempVec {};
783 if (not truthParticleContainer.
isValid()) {
786 tempVec.insert(tempVec.begin(), truthParticleContainer->
begin(), truthParticleContainer->
end());
796 const auto&
links =
event->truthParticleLinks();
797 tempVec.reserve(
event->nTruthParticles());
798 for (
const auto& link :
links) {
800 tempVec.push_back(*link);
809 if (truthPileupEventContainer.
isValid()) {
810 const unsigned int nPileup = truthPileupEventContainer->
size();
811 tempVec.reserve(nPileup * 200);
812 for (
unsigned int i(0);
i != nPileup; ++
i) {
813 const auto *eventPileup = truthPileupEventContainer->
at(
i);
815 int ntruth = eventPileup->nTruthParticles();
816 ATH_MSG_VERBOSE(
"Adding " << ntruth <<
" truth particles from TruthPileupEvents container");
817 const auto&
links = eventPileup->truthParticleLinks();
818 for (
const auto& link :
links) {
820 tempVec.push_back(*link);
835 std::pair<const std::vector<const xAOD::TruthVertex*>,
const std::vector<const xAOD::TruthVertex*>>
838 std::vector<const xAOD::TruthVertex*> truthHSVertices = {};
839 truthHSVertices.reserve(5);
840 std::vector<const xAOD::TruthVertex*> truthPUVertices = {};
841 truthPUVertices.reserve(100);
864 if (truthEventContainer.
isValid()) {
865 for (
const auto *
const evt : *truthEventContainer) {
866 truthVtx =
evt->signalProcessVertex();
868 truthHSVertices.push_back(truthVtx);
882 if (truthPileupEventContainer.
isValid()) {
883 for (
const auto *
const evt : *truthPileupEventContainer) {
888 size_t i_vtx = 0;
size_t n_vtx =
evt->nTruthVertices();
889 while(!truthVtx && i_vtx<n_vtx){
890 truthVtx =
evt->truthVertex(i_vtx);
895 truthPUVertices.push_back(truthVtx);
905 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);
916 std::vector<int>& cutFlow) {
918 if (cutFlow.empty()) {
919 names.emplace_back(
"preCut");
920 cutFlow.push_back(0);
921 for (
unsigned int i = 0;
i !=
accept.getNCuts(); ++
i) {
922 cutFlow.push_back(0);
928 bool cutPositive =
true;
929 for (
unsigned int i = 0;
i != (
accept.getNCuts() + 1); ++
i) {
945 ATH_MSG_INFO(
"Requesting cut value outside of configured eta range: clamping eta = "
946 << std::abs(truth.
eta()) <<
" to eta= " <<
absEta);
957 const std::vector<const xAOD::TruthParticle*>& truthParticles,
960 float beamSpotWeight) {
965 if (truthParticles.empty()) {
966 ATH_MSG_WARNING(
"No entries in TruthParticles truth particle container. Skipping jet plots.");
967 return StatusCode::SUCCESS;
973 return StatusCode::SUCCESS;
986 isBjet = (btagLabel(*thisJet) == 5);
994 if (not
el.isValid())
continue;
1005 bool isEfficient(
false);
1007 for (
const auto *thisTrack: tracks) {
1013 if (associatedTruth and associatedTruth == truth) {
1014 float prob = getMatchingProbability(*thisTrack);
1022 bool truthIsFromB =
false;
1024 truthIsFromB =
true;
1026 m_monPlots->fillEfficiency(*truth, *thisJet, isEfficient, isBjet, truthIsFromB, beamSpotWeight);
1036 if (thisJet->p4().DeltaR(thisTrack->p4()) >
m_maxTrkJetDR) {
1040 float prob = getMatchingProbability(*thisTrack);
1044 const bool unlinked = (associatedTruth==
nullptr);
1046 bool truthIsFromB =
false;
1048 truthIsFromB =
true;
1050 m_monPlots->fill(*thisTrack, *thisJet, isBjet,
isFake, unlinked, truthIsFromB, beamSpotWeight);
1051 if (associatedTruth){
1052 m_monPlots->fillFakeRate(*thisTrack, *thisJet,
isFake, isBjet, truthIsFromB, beamSpotWeight);
1058 return StatusCode::SUCCESS;
1063 const float jetPt =
jet.pt();
1064 const float jetEta = std::abs(
jet.eta());