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) :
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() );
368 if(pie.
isValid()) actualMu = pie->actualInteractionsPerCrossing();
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;
385 beamSpotWeight = pie->beamSpotWeight();
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;
395 if (vertices.
isValid() and not vertices->empty()) {
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++;
478 if (thisTrack->pt() >= (1 * Gaudi::Units::GeV))
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;
509 const bool isFake = not std::isnan(prob) ? (prob <
m_lowProb) :
true;
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);
572 std::lock_guard<std::mutex> lock(
m_mutex);
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);
609 if (not std::isnan(prob) && prob >
m_lowProb) {
610 matchingProbability = prob;
612 matchedTrack = thisTrack;
619 <<thisTruth->
py()/ Gaudi::Units::GeV<<
","<<thisTruth->
pz()/ Gaudi::Units::GeV<<
","
620 <<thisTruth->
e()/ Gaudi::Units::GeV<<
","<<thisTruth->
pt()/ Gaudi::Units::GeV<<
","
621 <<thisTruth->
eta()<<
","<<thisTruth->
phi()<<
","<<thisTruth->
m()/ Gaudi::Units::GeV<<
","
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) {
726 std::vector<HistData> hists =
m_monPlots->retrieveBookedHistograms();
727 for (
const auto& hist : hists) {
731 std::vector<EfficiencyData> effs =
m_monPlots->retrieveBookedEfficiencies();
732 for (
auto& eff : effs) {
739 std::vector<TreeData> trees =
m_monPlots->retrieveBookedTrees();
740 for (
auto& t : trees) {
745 return StatusCode::SUCCESS;
770 return StatusCode::SUCCESS;
773const 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);
835std::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);
923 names.push_back((std::string) accept.getCutName(i));
928 bool cutPositive =
true;
929 for (
unsigned int i = 0; i != (accept.getNCuts() + 1); ++i) {
933 if (accept.getCutResult(i)) {
942 double absEta = std::abs(truth.
eta());
945 ATH_MSG_INFO(
"Requesting cut value outside of configured eta range: clamping eta = "
946 << std::abs(truth.
eta()) <<
" to eta= " << absEta);
949 const auto pVal = std::lower_bound(
m_etaBins.value().begin(),
m_etaBins.value().end(), absEta);
950 const int bin = std::distance(
m_etaBins.value().begin(), pVal) - 1;
951 ATH_MSG_DEBUG(
"Checking (abs(eta)/bin) = (" << absEta <<
"," <<
bin <<
")");
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;
978 for (
const xAOD::Jet *
const thisJet: *jets) {
986 isBjet = (btagLabel(*thisJet) == 5);
994 if (not el.isValid())
continue;
1003 if(!accept)
continue;
1005 bool isEfficient(
false);
1007 for (
const auto *thisTrack: tracks) {
1013 if (associatedTruth and associatedTruth == truth) {
1014 float prob = getMatchingProbability(*thisTrack);
1015 if (not std::isnan(prob) && prob >
m_lowProb) {
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);
1041 if(std::isnan(prob)) prob = 0.0;
1044 const bool unlinked = (associatedTruth==
nullptr);
1045 const bool isFake = (associatedTruth && prob <
m_lowProb);
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());
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_VERBOSE(x)
#define ATH_MSG_WARNING(x)
Helper class to provide constant type-safe access to aux data.
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 doTrkInJetPlots_unlinked_bjets
bool doHitsFakeTracksPlots
int detailLevel
detail level (kept for compatibility)
bool doNtupleTruthToReco
Ntuple functionality.
bool doEfficienciesPerAuthor
per author plots
bool doMissingTruthFakePlots
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 (and unlinked)
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 doHitsUnlinkedTracksPlots
bool doTrackParametersPerAuthor
bool doTrkInJetPlots_unlinked
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