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(
"");
67 inRange(
const T& value,
const T& minVal,
const T& maxVal) {
68 return not ((value < minVal)or(value > maxVal));
73 inRange(
const T& value,
const T& absoluteMax) {
74 return not (std::abs(value) > absoluteMax);
82 binIndex(
const T& val,
const std::vector<T>& partitions) {
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);
101 2.7, 3.5, std::numeric_limits<float>::infinity()
108 const IInterface* parent) :
126 {not m_vertexContainerName.key().empty() &&
127 not m_hardScatterSelectionTool.empty()}));
137 (
nullptr,
static_cast<std::string
>(
m_dirName) +
static_cast<std::string
>(
m_folder),
148 std::vector<std::string> trk_decorations;
164 std::vector<std::string> required_float_track_decorations {
"d0",
"hitResiduals_residualLocX",
"d0err"};
165 std::vector<std::string> required_int_track_decorations {};
166 std::vector<std::string> required_float_truth_decorations {
"d0"};
167 std::vector<std::string> required_int_truth_decorations {};
168 std::vector<std::string> required_int_jet_decorations {
"HadronConeExclTruthLabelID"};
173 m_datfile <<
"EvtNumber,PdgID,Px,Py,Pz,E,Pt,Eta,Phi,Mass,numPU,numvtx,MatchProb"<<std::endl;
175 std::string empty_prefix;
189 return StatusCode::SUCCESS;
318 if (not trackHandle.
isValid()) {
320 return StatusCode::FAILURE;
331 ATH_MSG_WARNING(
"Shouldn't happen - EventInfo is buggy, setting mu to 0");
339 return StatusCode::SUCCESS;
343 std::vector<const xAOD::TruthParticle*> truthParticlesVec =
getTruthParticles(ctx);
352 unsigned int truthMu = 0;
355 truthMu =
static_cast<int>( truthPileupEventContainer->size() );
357 if(pie.
isValid()) actualMu = pie->actualInteractionsPerCrossing();
360 float puEvents = truthMu>0 ? truthMu : actualMu;
363 unsigned int nVertices = 0;
364 float beamSpotWeight = 1;
367 ATH_MSG_DEBUG(
"Getting number of pu interactings per event");
373 nVertices = not vertices->empty() ? vertices->size() : 0;
374 beamSpotWeight = pie->beamSpotWeight();
375 ATH_MSG_DEBUG(
"beamSpotWeight is equal to " << beamSpotWeight);
379 if(readDecorHandle.
isAvailable()) prwWeight = readDecorHandle(*pie);
380 ATH_MSG_DEBUG(
"Applying pileup weight equal to " << prwWeight);
381 beamSpotWeight *= prwWeight;
384 if (vertices.
isValid() and not vertices->empty()) {
385 ATH_MSG_DEBUG(
"Number of vertices retrieved for this event " << vertices->size());
390 primaryvertex = vertex;
401 ATH_MSG_DEBUG(
"Failed to find a hard scatter vertex in this event.");
407 std::pair<std::vector<const xAOD::TruthVertex*>, std::vector<const xAOD::TruthVertex*>> truthVertices =
getTruthVertices(ctx);
408 std::vector<const xAOD::TruthVertex*> truthHSVertices = truthVertices.first;
409 std::vector<const xAOD::TruthVertex*> truthPUVertices = truthVertices.second;
416 m_monPlots->fill(*vertices, primaryvertex, truthHSVertices, truthPUVertices, actualMu, beamSpotWeight);
420 m_monPlots->fill(*vertices, truthMu, actualMu, beamSpotWeight);
437 const auto& stdVertexContainer = truthVrt->stdcont();
439 auto findVtx = std::find_if(stdVertexContainer.rbegin(), stdVertexContainer.rend(), acceptTruthVertex);
440 truthVertex = (findVtx == stdVertexContainer.rend()) ? nullptr : *findVtx;
444 if (not truthVertex)
ATH_MSG_INFO (
"Truth vertex did not pass cuts");
449 unsigned int nSelectedTruthTracks(0), nSelectedRecoTracks(0), nSelectedMatchedTracks(0), nAssociatedTruth(0), nMissingAssociatedTruth(0), nTruths(0);
458 std::map<const xAOD::TruthParticle*, std::vector<const xAOD::TrackParticle*>> cacheTruthMatching {};
460 std::vector<const xAOD::TrackParticle*> selectedTracks {};
461 selectedTracks.reserve(tracks->
size());
462 unsigned int nTrackTOT = 0;
463 unsigned int nTrackCentral = 0;
464 unsigned int nTrackPt1GeV = 0;
465 for (
const auto *
const thisTrack: *tracks) {
471 selectedTracks.push_back(thisTrack);
473 nSelectedRecoTracks++;
477 if (thisTrack->pt() >= (1 * Gaudi::Units::GeV))
479 if (std::abs(thisTrack->eta()) < 2.5)
482 m_monPlots->fill(*thisTrack, puEvents, nVertices, beamSpotWeight);
484 float prob = getMatchingProbability(*thisTrack);
487 if (associatedTruth) {
498 nSelectedMatchedTracks++;
499 bool truthIsFromB =
false;
503 m_monPlots->fill(*thisTrack, *associatedTruth, truthIsFromB, puEvents, beamSpotWeight);
507 const bool isAssociatedTruth = associatedTruth !=
nullptr;
508 const bool isFake = not std::isnan(prob) ? (prob <
m_lowProb) :
true;
510 if(!isAssociatedTruth) nMissingAssociatedTruth++;
511 m_monPlots->fillFakeRate(*thisTrack, isFake, puEvents, beamSpotWeight);
517 if (isAssociatedTruth) {
522 auto cachedAssoc = cacheTruthMatching.find(associatedTruth);
524 if (cachedAssoc == cacheTruthMatching.end()) {
526 cacheTruthMatching[associatedTruth] = {thisTrack};
530 cachedAssoc->second.push_back(thisTrack);
535 m_monPlots->fillNtuple(*thisTrack, primaryvertex);
543 for (
auto& cachedAssoc: cacheTruthMatching) {
550 std::sort(cachedAssoc.second.begin(), cachedAssoc.second.end(),
556 for (
int itrack = 0; itrack < (int) cachedAssoc.second.size(); itrack++) {
560 m_monPlots->fillNtuple(*thisTrack, *thisTruth, primaryvertex, itrack);
566 m_monPlots->fill(nTrackTOT, nTrackCentral, nTrackPt1GeV, truthMu, actualMu, nVertices, beamSpotWeight);
572 m_truthCutFlow.merge(std::move(tmp_truth_cutflow));
578 for (
int itruth = 0; itruth < (int) truthParticlesVec.size(); itruth++) {
585 ++nSelectedTruthTracks;
586 bool isEfficient(
false);
587 float matchingProbability{};
591 auto cachedAssoc = cacheTruthMatching.find(thisTruth);
593 if (cachedAssoc == cacheTruthMatching.end()) {
595 cacheTruthMatching[thisTruth] = {};
597 m_monPlots->fillDuplicate(*thisTruth, cacheTruthMatching[thisTruth], beamSpotWeight);
604 for (
const auto& thisTrack: selectedTracks) {
606 if (associatedTruth && associatedTruth == thisTruth) {
607 float prob = getMatchingProbability(*thisTrack);
608 if (not std::isnan(prob) && prob >
m_lowProb) {
609 matchingProbability = prob;
611 matchedTrack = thisTrack;
618 <<thisTruth->
py()/ Gaudi::Units::GeV<<
","<<thisTruth->
pz()/ Gaudi::Units::GeV<<
","
619 <<thisTruth->
e()/ Gaudi::Units::GeV<<
","<<thisTruth->
pt()/ Gaudi::Units::GeV<<
","
620 <<thisTruth->
eta()<<
","<<thisTruth->
phi()<<
","<<thisTruth->
m()/ Gaudi::Units::GeV<<
","
621 <<puEvents<<
","<<nVertices<<
","<<matchingProbability<<std::endl;
624 ATH_MSG_ERROR(
"An error occurred: Truth particle for tracking efficiency calculation is a nullptr");
626 else if (isEfficient && !matchedTrack){
627 ATH_MSG_ERROR(
"Something went wrong - we log a truth particle as reconstructed, but the reco track is a nullptr! Bailing out... ");
630 ATH_MSG_DEBUG(
"Filling efficiency plots info monitoring plots");
631 m_monPlots->fillEfficiency(*thisTruth, matchedTrack, isEfficient, truthMu, actualMu, beamSpotWeight);
633 ATH_MSG_DEBUG(
"Filling technical efficiency plots info monitoring plots");
637 m_monPlots->fillTechnicalEfficiency(*thisTruth, isEfficient,
638 truthMu, actualMu, beamSpotWeight);
641 ATH_MSG_DEBUG(
"Cannot fill technical efficiency. Missing si hit information for truth particle.");
659 if (nSelectedRecoTracks == nMissingAssociatedTruth) {
664 ATH_MSG_DEBUG(nAssociatedTruth <<
" tracks out of " << tracks->
size() <<
" had associated truth.");
683 return StatusCode::SUCCESS;
709 ATH_MSG_DEBUG(
"Selected by pileup switch decoration requested from a truth particle but not available");
716 for (
const auto& thisTruth: truthParticles) {
725 std::vector<HistData> hists =
m_monPlots->retrieveBookedHistograms();
726 for (
const auto& hist : hists) {
730 std::vector<EfficiencyData> effs =
m_monPlots->retrieveBookedEfficiencies();
731 for (
auto& eff : effs) {
738 std::vector<TreeData> trees =
m_monPlots->retrieveBookedTrees();
739 for (
auto& t : trees) {
744 return StatusCode::SUCCESS;
769 return StatusCode::SUCCESS;
772const std::vector<const xAOD::TruthParticle*>
775 std::vector<const xAOD::TruthParticle*> tempVec {};
782 if (not truthParticleContainer.
isValid()) {
785 tempVec.insert(tempVec.begin(), truthParticleContainer->begin(), truthParticleContainer->end());
795 const auto& links =
event->truthParticleLinks();
796 tempVec.reserve(event->nTruthParticles());
797 for (
const auto& link : links) {
799 tempVec.push_back(*link);
808 if (truthPileupEventContainer.
isValid()) {
809 const unsigned int nPileup = truthPileupEventContainer->size();
810 tempVec.reserve(nPileup * 200);
811 for (
unsigned int i(0); i != nPileup; ++i) {
812 const auto *eventPileup = truthPileupEventContainer->at(i);
814 int ntruth = eventPileup->nTruthParticles();
815 ATH_MSG_VERBOSE(
"Adding " << ntruth <<
" truth particles from TruthPileupEvents container");
816 const auto& links = eventPileup->truthParticleLinks();
817 for (
const auto& link : links) {
819 tempVec.push_back(*link);
834std::pair<std::vector<const xAOD::TruthVertex*>, std::vector<const xAOD::TruthVertex*>>
837 std::vector<const xAOD::TruthVertex*> truthHSVertices = {};
838 truthHSVertices.reserve(5);
839 std::vector<const xAOD::TruthVertex*> truthPUVertices = {};
840 truthPUVertices.reserve(100);
863 if (truthEventContainer.
isValid()) {
864 for (
const auto *
const evt : *truthEventContainer) {
865 truthVtx = evt->signalProcessVertex();
867 truthHSVertices.push_back(truthVtx);
881 if (truthPileupEventContainer.
isValid()) {
882 for (
const auto *
const evt : *truthPileupEventContainer) {
887 size_t i_vtx = 0;
size_t n_vtx = evt->nTruthVertices();
888 while(!truthVtx && i_vtx<n_vtx){
889 truthVtx = evt->truthVertex(i_vtx);
894 truthPUVertices.push_back(truthVtx);
904 return {std::move(truthHSVertices), std::move(truthPUVertices)};
914 std::vector<int>& cutFlow) {
916 if (cutFlow.empty()) {
917 names.emplace_back(
"preCut");
918 cutFlow.push_back(0);
919 for (
unsigned int i = 0; i != accept.getNCuts(); ++i) {
920 cutFlow.push_back(0);
921 names.push_back((std::string) accept.getCutName(i));
926 bool cutPositive =
true;
927 for (
unsigned int i = 0; i != (accept.getNCuts() + 1); ++i) {
931 if (accept.getCutResult(i)) {
940 double absEta = std::abs(truth.
eta());
943 ATH_MSG_INFO(
"Requesting cut value outside of configured eta range: clamping eta = "
944 << std::abs(truth.
eta()) <<
" to eta= " << absEta);
947 const auto pVal = std::lower_bound(
m_etaBins.value().begin(),
m_etaBins.value().end(), absEta);
948 const int bin = std::distance(
m_etaBins.value().begin(), pVal) - 1;
949 ATH_MSG_DEBUG(
"Checking (abs(eta)/bin) = (" << absEta <<
"," <<
bin <<
")");
955 const std::vector<const xAOD::TruthParticle*>& truthParticles,
958 float beamSpotWeight) {
963 if (truthParticles.empty()) {
964 ATH_MSG_WARNING(
"No entries in TruthParticles truth particle container. Skipping jet plots.");
965 return StatusCode::SUCCESS;
971 return StatusCode::SUCCESS;
976 for (
const xAOD::Jet *
const thisJet: *jets) {
984 isBjet = (btagLabel(*thisJet) == 5);
992 if (not el.isValid())
continue;
998 if (std::find(truthParticles.begin(), truthParticles.end(), truth) == truthParticles.end()) {
1009 if(!accept)
continue;
1011 bool isEfficient(
false);
1013 for (
const auto *thisTrack: tracks) {
1019 if (associatedTruth and associatedTruth == truth) {
1020 float prob = getMatchingProbability(*thisTrack);
1021 if (not std::isnan(prob) && prob >
m_lowProb) {
1028 bool truthIsFromB =
false;
1030 truthIsFromB =
true;
1032 m_monPlots->fillEfficiency(*truth, *thisJet, isEfficient, isBjet, truthIsFromB, beamSpotWeight);
1042 if (thisJet->p4().DeltaR(thisTrack->p4()) >
m_maxTrkJetDR) {
1046 float prob = getMatchingProbability(*thisTrack);
1047 if(std::isnan(prob)) prob = 0.0;
1050 const bool isFake = (associatedTruth && prob <
m_lowProb);
1051 bool truthIsFromB =
false;
1053 truthIsFromB =
true;
1055 m_monPlots->fill(*thisTrack, *thisJet, isBjet, isFake, truthIsFromB, beamSpotWeight);
1056 if (associatedTruth){
1057 m_monPlots->fillFakeRate(*thisTrack, *thisJet, isFake, isBjet, truthIsFromB, beamSpotWeight);
1063 return StatusCode::SUCCESS;
1068 const float jetPt =
jet.pt();
1069 const float jetEta = std::abs(
jet.eta());
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_DEBUG(x,...)
#define ATH_MSG_ERROR(x,...)
#define ATH_MSG_WARNING(x,...)
#define ATH_MSG_VERBOSE(x,...)
#define ATH_MSG_INFO(x,...)
Helper class to provide constant type-safe access to aux data.
virtual void lock()=0
Interface to allow an object to lock itself when made const in SG.
header file for class of same name
bool inRange(const double *boundaries, const double value, const double tolerance=0.02)
ServiceHandle< StoreGateSvc > & evtStore()
size_type size() const noexcept
Returns the number of elements in the collection.
ElementLink implementation for ROOT usage.
static void neededTrackParticleDecorations(std::vector< std::string > &decorations)
const xAOD::TruthParticle * getTruth(const xAOD::TrackParticle *const trackParticle)
Helper class to provide constant type-safe access to aux data.
bool isAvailable(const ELT &e) const
Test to see if this variable exists in the store.
Handle class for reading a decoration on an object.
bool isAvailable()
Test to see if this variable exists in the store, for the referenced object.
virtual bool isValid() override final
Can the handle be successfully dereferenced?
const_pointer_type cptr()
Dereference the pointer.
const_pointer_type get() const
Dereference the pointer, but don't cache anything.
@ IS_SIMULATION
true: simulation, false: data
uint64_t eventNumber() const
The current event's event number.
virtual double m() const override final
The mass of the particle.
int pdgId() const
PDG ID code.
float px() const
The x component of the particle's momentum.
virtual double e() const override final
The total energy of the particle.
virtual double pt() const override final
The transverse momentum ( ) of the particle.
float py() const
The y component of the particle's momentum.
virtual double eta() const override final
The pseudorapidity ( ) of the particle.
virtual double phi() const override final
The azimuthal angle ( ) of the particle.
virtual FourMom_t p4() const override final
The full 4-momentum of the particle.
float pz() const
The z component of the particle's momentum.
float z() const
Vertex longitudinal distance along the beam line form the origin.
float y() const
Vertex y displacement.
float x() const
Vertex x displacement.
void addReadDecoratorHandleKeys(T_Parent &parent, const SG::ReadHandleKey< T_Cont > &container_key, const std::string &prefix, const std::vector< std::string > &decor_names, std::vector< SG::ReadDecorHandleKey< T_Cont > > &decor_out)
float safelyGetEta(const T &pTrk, const float safePtThreshold=0.1)
Safely get eta.
unsigned int binIndex(const T &val, const std::vector< T > &partitions)
general utility function to return bin index given a value and the upper endpoints of each bin
HardScatterType classifyHardScatter(const xAOD::VertexContainer &vxContainer)
void sort(typename DataModel_detail::iterator< DVL > beg, typename DataModel_detail::iterator< DVL > end)
Specialization of sort for DataVector/List.
Jet_v1 Jet
Definition of the current "jet version".
EventInfo_v1 EventInfo
Definition of the latest event info version.
TruthVertex_v1 TruthVertex
Typedef to implementation.
TrackParticle_v1 TrackParticle
Reference the current persistent version:
Vertex_v1 Vertex
Define the latest version of the vertex class.
TruthEvent_v1 TruthEvent
Typedef to implementation.
TruthParticle_v1 TruthParticle
Typedef to implementation.
TrackParticleContainer_v1 TrackParticleContainer
Definition of the current "TrackParticle container version".
JetContainer_v1 JetContainer
Definition of the current "jet container version".
implementation file for function of same name
helper struct - steer the configuration from the parent tool's side
bool doTrkInJetPlots_matched_bjets
bool doHitsFakeTracksPlots
int detailLevel
detail level (kept for compatibility)
bool doNtupleTruthToReco
Ntuple functionality.
bool doEfficienciesPerAuthor
per author plots
bool doHitsRecoTracksPlots
bool doHitsRecoTracksPlotsPerAuthor
bool doEffPlots
Efficiency and duplicate plots - require truth, optionally matching reco.
bool doTrkInJetPlots_bjets
bool doTrkInJetPlots_matched
bool doFakePlots
Fake plots.
bool doResolutionPlotPrim
Resolution and "matched track" plots - filled if both reco and truth exist.
bool doTrkInJetPlots_fake_bjets
bool doVertexTruthMatchingPlots
Vertexing plots - truth requirement.
bool doVertexPlots
Vertexing plots - no truth requirement.
bool doHardScatterVertexTruthMatchingPlots
bool doTrkInJetPlots_truthFromB
bool doTrkInJetPlots_fake
bool doTrackParametersPerAuthor
bool doTrackParameters
Plots for (selected) tracks, not necessarily truth matched.
bool doTrkInJetPlots
Plots for tracks in jets.
bool doResolutionPlotPrim_truthFromB
bool doHardScatterVertexPlots
bool doResolutionPlotSecd
bool doResolutionsPerAuthor
bool doHitsMatchedTracksPlots