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) :
136 (
nullptr,
static_cast<std::string
>(
m_dirName) +
static_cast<std::string
>(
m_folder),
147 std::vector<std::string> trk_decorations;
163 std::vector<std::string> required_float_track_decorations {
"d0",
"hitResiduals_residualLocX",
"d0err"};
164 std::vector<std::string> required_int_track_decorations {};
165 std::vector<std::string> required_float_truth_decorations {
"d0"};
166 std::vector<std::string> required_int_truth_decorations {};
167 std::vector<std::string> required_int_jet_decorations {
"HadronConeExclTruthLabelID"};
172 m_datfile <<
"EvtNumber,PdgID,Px,Py,Pz,E,Pt,Eta,Phi,Mass,numPU,numvtx,MatchProb"<<std::endl;
174 std::string empty_prefix;
188 return StatusCode::SUCCESS;
316 if (not trackHandle.
isValid()) {
318 return StatusCode::FAILURE;
329 ATH_MSG_WARNING(
"Shouldn't happen - EventInfo is buggy, setting mu to 0");
337 return StatusCode::SUCCESS;
341 std::vector<const xAOD::TruthParticle*> truthParticlesVec =
getTruthParticles(ctx);
350 unsigned int truthMu = 0;
353 truthMu =
static_cast<int>( truthPileupEventContainer->size() );
355 if(pie.
isValid()) actualMu = pie->actualInteractionsPerCrossing();
358 float puEvents = truthMu>0 ? truthMu : actualMu;
361 unsigned int nVertices = 0;
362 float beamSpotWeight = 1;
365 ATH_MSG_DEBUG(
"Getting number of pu interactings per event");
371 nVertices = not vertices->empty() ? vertices->size() : 0;
372 beamSpotWeight = pie->beamSpotWeight();
373 ATH_MSG_DEBUG(
"beamSpotWeight is equal to " << beamSpotWeight);
377 if(readDecorHandle.
isAvailable()) prwWeight = readDecorHandle(*pie);
378 ATH_MSG_DEBUG(
"Applying pileup weight equal to " << prwWeight);
379 beamSpotWeight *= prwWeight;
382 if (vertices.
isValid() and not vertices->empty()) {
383 ATH_MSG_DEBUG(
"Number of vertices retrieved for this event " << vertices->size());
389 ATH_MSG_DEBUG(
"Failed to find a hard scatter vertex in this event.");
395 std::pair<std::vector<const xAOD::TruthVertex*>, std::vector<const xAOD::TruthVertex*>> truthVertices =
getTruthVertices(ctx);
396 std::vector<const xAOD::TruthVertex*> truthHSVertices = truthVertices.first;
397 std::vector<const xAOD::TruthVertex*> truthPUVertices = truthVertices.second;
404 m_monPlots->fill(*vertices, primaryvertex, truthHSVertices, truthPUVertices, actualMu, beamSpotWeight);
408 m_monPlots->fill(*vertices, truthMu, actualMu, beamSpotWeight);
425 const auto& stdVertexContainer = truthVrt->stdcont();
427 auto findVtx = std::find_if(stdVertexContainer.rbegin(), stdVertexContainer.rend(), acceptTruthVertex);
428 truthVertex = (findVtx == stdVertexContainer.rend()) ? nullptr : *findVtx;
432 if (not truthVertex)
ATH_MSG_INFO (
"Truth vertex did not pass cuts");
437 unsigned int nSelectedTruthTracks(0), nSelectedRecoTracks(0), nSelectedMatchedTracks(0), nAssociatedTruth(0), nMissingAssociatedTruth(0), nTruths(0);
446 std::map<const xAOD::TruthParticle*, std::vector<const xAOD::TrackParticle*>> cacheTruthMatching {};
448 std::vector<const xAOD::TrackParticle*> selectedTracks {};
449 selectedTracks.reserve(tracks->
size());
450 unsigned int nTrackTOT = 0;
451 unsigned int nTrackCentral = 0;
452 unsigned int nTrackPt1GeV = 0;
453 for (
const auto *
const thisTrack: *tracks) {
459 selectedTracks.push_back(thisTrack);
461 nSelectedRecoTracks++;
465 if (thisTrack->pt() >= (1 * Gaudi::Units::GeV))
467 if (std::abs(thisTrack->eta()) < 2.5)
470 m_monPlots->fill(*thisTrack, puEvents, nVertices, beamSpotWeight);
472 float prob = getMatchingProbability(*thisTrack);
475 if (associatedTruth) {
486 nSelectedMatchedTracks++;
487 bool truthIsFromB =
false;
491 m_monPlots->fill(*thisTrack, *associatedTruth, truthIsFromB, puEvents, beamSpotWeight);
495 const bool isAssociatedTruth = associatedTruth !=
nullptr;
496 const bool isFake = not std::isnan(prob) ? (prob <
m_lowProb) :
true;
498 if(!isAssociatedTruth) nMissingAssociatedTruth++;
499 m_monPlots->fillFakeRate(*thisTrack, isFake, puEvents, beamSpotWeight);
505 if (isAssociatedTruth) {
510 auto cachedAssoc = cacheTruthMatching.find(associatedTruth);
512 if (cachedAssoc == cacheTruthMatching.end()) {
514 cacheTruthMatching[associatedTruth] = {thisTrack};
518 cachedAssoc->second.push_back(thisTrack);
523 m_monPlots->fillNtuple(*thisTrack, primaryvertex);
531 for (
auto& cachedAssoc: cacheTruthMatching) {
538 std::sort(cachedAssoc.second.begin(), cachedAssoc.second.end(),
544 for (
int itrack = 0; itrack < (int) cachedAssoc.second.size(); itrack++) {
548 m_monPlots->fillNtuple(*thisTrack, *thisTruth, primaryvertex, itrack);
554 m_monPlots->fill(nTrackTOT, nTrackCentral, nTrackPt1GeV, truthMu, actualMu, nVertices, beamSpotWeight);
560 m_truthCutFlow.merge(std::move(tmp_truth_cutflow));
566 for (
int itruth = 0; itruth < (int) truthParticlesVec.size(); itruth++) {
573 ++nSelectedTruthTracks;
574 bool isEfficient(
false);
575 float matchingProbability{};
579 auto cachedAssoc = cacheTruthMatching.find(thisTruth);
581 if (cachedAssoc == cacheTruthMatching.end()) {
583 cacheTruthMatching[thisTruth] = {};
585 m_monPlots->fillDuplicate(*thisTruth, cacheTruthMatching[thisTruth], beamSpotWeight);
592 for (
const auto& thisTrack: selectedTracks) {
594 if (associatedTruth && associatedTruth == thisTruth) {
595 float prob = getMatchingProbability(*thisTrack);
596 if (not std::isnan(prob) && prob >
m_lowProb) {
597 matchingProbability = prob;
599 matchedTrack = thisTrack;
606 <<thisTruth->
py()/ Gaudi::Units::GeV<<
","<<thisTruth->
pz()/ Gaudi::Units::GeV<<
","
607 <<thisTruth->
e()/ Gaudi::Units::GeV<<
","<<thisTruth->
pt()/ Gaudi::Units::GeV<<
","
608 <<thisTruth->
eta()<<
","<<thisTruth->
phi()<<
","<<thisTruth->
m()/ Gaudi::Units::GeV<<
","
609 <<puEvents<<
","<<nVertices<<
","<<matchingProbability<<std::endl;
612 ATH_MSG_ERROR(
"An error occurred: Truth particle for tracking efficiency calculation is a nullptr");
614 else if (isEfficient && !matchedTrack){
615 ATH_MSG_ERROR(
"Something went wrong - we log a truth particle as reconstructed, but the reco track is a nullptr! Bailing out... ");
618 ATH_MSG_DEBUG(
"Filling efficiency plots info monitoring plots");
619 m_monPlots->fillEfficiency(*thisTruth, matchedTrack, isEfficient, truthMu, actualMu, beamSpotWeight);
621 ATH_MSG_DEBUG(
"Filling technical efficiency plots info monitoring plots");
625 m_monPlots->fillTechnicalEfficiency(*thisTruth, isEfficient,
626 truthMu, actualMu, beamSpotWeight);
629 ATH_MSG_DEBUG(
"Cannot fill technical efficiency. Missing si hit information for truth particle.");
647 if (nSelectedRecoTracks == nMissingAssociatedTruth) {
652 ATH_MSG_DEBUG(nAssociatedTruth <<
" tracks out of " << tracks->
size() <<
" had associated truth.");
671 return StatusCode::SUCCESS;
697 ATH_MSG_DEBUG(
"Selected by pileup switch decoration requested from a truth particle but not available");
704 for (
const auto& thisTruth: truthParticles) {
713 std::vector<HistData> hists =
m_monPlots->retrieveBookedHistograms();
714 for (
const auto& hist : hists) {
718 std::vector<EfficiencyData> effs =
m_monPlots->retrieveBookedEfficiencies();
719 for (
auto& eff : effs) {
726 std::vector<TreeData> trees =
m_monPlots->retrieveBookedTrees();
727 for (
auto& t : trees) {
732 return StatusCode::SUCCESS;
757 return StatusCode::SUCCESS;
760const std::vector<const xAOD::TruthParticle*>
763 std::vector<const xAOD::TruthParticle*> tempVec {};
770 if (not truthParticleContainer.
isValid()) {
773 tempVec.insert(tempVec.begin(), truthParticleContainer->begin(), truthParticleContainer->end());
783 const auto& links =
event->truthParticleLinks();
784 tempVec.reserve(event->nTruthParticles());
785 for (
const auto& link : links) {
787 tempVec.push_back(*link);
796 if (truthPileupEventContainer.
isValid()) {
797 const unsigned int nPileup = truthPileupEventContainer->size();
798 tempVec.reserve(nPileup * 200);
799 for (
unsigned int i(0); i != nPileup; ++i) {
800 const auto *eventPileup = truthPileupEventContainer->at(i);
802 int ntruth = eventPileup->nTruthParticles();
803 ATH_MSG_VERBOSE(
"Adding " << ntruth <<
" truth particles from TruthPileupEvents container");
804 const auto& links = eventPileup->truthParticleLinks();
805 for (
const auto& link : links) {
807 tempVec.push_back(*link);
822std::pair<const std::vector<const xAOD::TruthVertex*>,
const std::vector<const xAOD::TruthVertex*>>
825 std::vector<const xAOD::TruthVertex*> truthHSVertices = {};
826 truthHSVertices.reserve(5);
827 std::vector<const xAOD::TruthVertex*> truthPUVertices = {};
828 truthPUVertices.reserve(100);
851 if (truthEventContainer.
isValid()) {
852 for (
const auto *
const evt : *truthEventContainer) {
853 truthVtx = evt->signalProcessVertex();
855 truthHSVertices.push_back(truthVtx);
869 if (truthPileupEventContainer.
isValid()) {
870 for (
const auto *
const evt : *truthPileupEventContainer) {
875 size_t i_vtx = 0;
size_t n_vtx = evt->nTruthVertices();
876 while(!truthVtx && i_vtx<n_vtx){
877 truthVtx = evt->truthVertex(i_vtx);
882 truthPUVertices.push_back(truthVtx);
892 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);
903 std::vector<int>& cutFlow) {
905 if (cutFlow.empty()) {
906 names.emplace_back(
"preCut");
907 cutFlow.push_back(0);
908 for (
unsigned int i = 0; i != accept.getNCuts(); ++i) {
909 cutFlow.push_back(0);
910 names.push_back((std::string) accept.getCutName(i));
915 bool cutPositive =
true;
916 for (
unsigned int i = 0; i != (accept.getNCuts() + 1); ++i) {
920 if (accept.getCutResult(i)) {
929 double absEta = std::abs(truth.
eta());
932 ATH_MSG_INFO(
"Requesting cut value outside of configured eta range: clamping eta = "
933 << std::abs(truth.
eta()) <<
" to eta= " << absEta);
936 const auto pVal = std::lower_bound(
m_etaBins.value().begin(),
m_etaBins.value().end(), absEta);
937 const int bin = std::distance(
m_etaBins.value().begin(), pVal) - 1;
938 ATH_MSG_DEBUG(
"Checking (abs(eta)/bin) = (" << absEta <<
"," <<
bin <<
")");
944 const std::vector<const xAOD::TruthParticle*>& truthParticles,
947 float beamSpotWeight) {
952 if (truthParticles.empty()) {
953 ATH_MSG_WARNING(
"No entries in TruthParticles truth particle container. Skipping jet plots.");
954 return StatusCode::SUCCESS;
960 return StatusCode::SUCCESS;
965 for (
const xAOD::Jet *
const thisJet: *jets) {
973 isBjet = (btagLabel(*thisJet) == 5);
981 if (not el.isValid())
continue;
990 if(!accept)
continue;
992 bool isEfficient(
false);
994 for (
const auto *thisTrack: tracks) {
1000 if (associatedTruth and associatedTruth == truth) {
1001 float prob = getMatchingProbability(*thisTrack);
1002 if (not std::isnan(prob) && prob >
m_lowProb) {
1009 bool truthIsFromB =
false;
1011 truthIsFromB =
true;
1013 m_monPlots->fillEfficiency(*truth, *thisJet, isEfficient, isBjet, truthIsFromB, beamSpotWeight);
1023 if (thisJet->p4().DeltaR(thisTrack->p4()) >
m_maxTrkJetDR) {
1027 float prob = getMatchingProbability(*thisTrack);
1028 if(std::isnan(prob)) prob = 0.0;
1031 const bool isFake = (associatedTruth && prob <
m_lowProb);
1032 bool truthIsFromB =
false;
1034 truthIsFromB =
true;
1036 m_monPlots->fill(*thisTrack, *thisJet, isBjet, isFake, 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());
#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.
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