diff --git a/README.md b/README.md index 103171fab8d..fd0d2a89b2e 100644 --- a/README.md +++ b/README.md @@ -77,3 +77,9 @@ To apply formatting to a file: ``` clang-format -i -style=file /path/to/file.cpp ``` +Note that this will reformat the entire file. +If you are preparing a PR, you can instead use +``` + git clang-format --style=file $(git merge-base upstream/master HEAD) +``` +to only format the lines you changed (avoids bloating your PR - unless you wish to clean up the whole file 😉). diff --git a/analyzers/dataframe/FCCAnalyses/MCParticle.h b/analyzers/dataframe/FCCAnalyses/MCParticle.h index 93dcdf2e296..c3e627f301e 100644 --- a/analyzers/dataframe/FCCAnalyses/MCParticle.h +++ b/analyzers/dataframe/FCCAnalyses/MCParticle.h @@ -3,8 +3,10 @@ #define MCPARTICLE_ANALYZERS_H #include +#include #include +#include "FCCAnalyses/Utils.h" #include "ROOT/RVec.hxx" #include "TLorentzVector.h" #include "edm4hep/MCParticleData.h" @@ -30,26 +32,32 @@ namespace MCParticle{ bool operator() (ROOT::VecOps::RVec in); }; + /// Template implementation for basic selection function on MC particle data. + using MCParticleSelection = Utils::selByPredicate; + /// select MCParticles with transverse momentum greater than a minimum value [GeV] - struct sel_pt { + struct sel_pt : MCParticleSelection { sel_pt(float arg_min_pt); - float m_min_pt = 20; //> transverse momentum threshold [GeV] - ROOT::VecOps::RVec operator() (ROOT::VecOps::RVec in); + }; + + /// select MCParticles with absolute pseudorapidity less than a max value + struct sel_eta : MCParticleSelection { + sel_eta(float arg_max_eta); }; /// select MCParticles with their status - struct sel_genStatus { + struct sel_genStatus : MCParticleSelection { sel_genStatus(int arg_status); - int m_status = 1; //> Generator status - ROOT::VecOps::RVec operator() (ROOT::VecOps::RVec in); }; /// select MCParticles with their PDG id - struct sel_pdgID { + struct sel_pdgID : MCParticleSelection { sel_pdgID(int arg_pdg, bool arg_chargeconjugate); - int m_pdg = 13; - bool m_chargeconjugate = true; - ROOT::VecOps::RVec operator() (ROOT::VecOps::RVec in); + }; + + /// select MCParticles with a non-zero charge + struct sel_charged : MCParticleSelection { + sel_charged(); }; /// get MC history tree for a given MCParticle index @@ -116,7 +124,6 @@ namespace MCParticle{ ROOT::VecOps::RVec in , ROOT::VecOps::RVec ind); - /// return the parent index of a given list of MC particles ROOT::VecOps::RVec get_parentid(ROOT::VecOps::RVec mcind, ROOT::VecOps::RVec mc, ROOT::VecOps::RVec parents); @@ -210,6 +217,16 @@ namespace MCParticle{ /// return the list of stable particles from the decay of a mother particle, looking at the full decay chain recursively. i is the mother index in the Particle block std::vector get_list_of_stable_particles_from_decay( int i, ROOT::VecOps::RVec in, ROOT::VecOps::RVec ind) ; + /// return the list of stable particles from the decays of a mother particle, + /// looking at the full decay chain recursively. i is the list of mother + /// indices to process in the Particle block. Will return a vector of vectors + /// - each vector is the set of children for one of the mothers in the input + /// vector + ROOT::VecOps::RVec> + get_lists_of_stable_particles_from_decays( + ROOT::VecOps::RVec i, ROOT::VecOps::RVec in, + ROOT::VecOps::RVec ind); + /// return the list of particles from the decay of a mother particle. i is the mother index in the Particle block. std::vector get_list_of_particles_from_decay( int i, ROOT::VecOps::RVec in, diff --git a/analyzers/dataframe/FCCAnalyses/ReconstructedParticle2MC.h b/analyzers/dataframe/FCCAnalyses/ReconstructedParticle2MC.h index 61845d91999..d1c5f1a20fc 100644 --- a/analyzers/dataframe/FCCAnalyses/ReconstructedParticle2MC.h +++ b/analyzers/dataframe/FCCAnalyses/ReconstructedParticle2MC.h @@ -125,12 +125,44 @@ namespace ReconstructedParticle2MC{ ROOT::VecOps::RVec reco, ROOT::VecOps::RVec mc) ; - /// select ReconstructedParticles matched to the (stable) MC particles whose indices are passed in a list - ROOT::VecOps::RVec selRP_matched_to_list( ROOT::VecOps::RVec mcParticles_indices, - ROOT::VecOps::RVec recind, - ROOT::VecOps::RVec mcind, - ROOT::VecOps::RVec reco, - ROOT::VecOps::RVec mc) ; + /// select ReconstructedParticles matched to the MC particles whose indices + /// are passed in a list + /// @param mcParticles_indices indices of the MC particles to look up + /// @param recind: reco index component of the MCRecoAssociations + /// @param mcind: mc index component of the MCRecoAssociations + /// @param reco: full reco particle list (ReconstructedParticles) + /// @param mc: full mc particle list (Particles) + /// @param require_stable: if set to true, will only match stable particles. + /// @return List of ReconstructedParticle candidates with length corresponding + /// to the number of *stable* MC particles in the mcParticles_indices vector. + /// In presence of unstable particles, no 1:1 correspondence. For + /// non-reconstructed stable MC particles, a dummy particle will be inserted. + /// If 1:1 length correspondence is required, set require_stable to false. + ROOT::VecOps::RVec selRP_matched_to_list( + const ROOT::VecOps::RVec &mcParticles_indices, + const ROOT::VecOps::RVec &recind, + const ROOT::VecOps::RVec &mcind, + const ROOT::VecOps::RVec &reco, + const ROOT::VecOps::RVec &mc, + bool require_stable = true); + /// select indices of ReconstructedParticles matched to the MC particles whose + /// indices are passed in a list + /// @param mcParticles_indices indices of the MC particles to look up + /// @param recind: reco index component of the MCRecoAssociations + /// @param mcind: mc index component of the MCRecoAssociations + /// @param mc: full mc particle list (Particles) + /// @param require_stable: if set to true, will only match stable particles. + /// @return List of ReconstructedParticle candidates with length corresponding + /// to the number of *stable* MC particles in the mcParticles_indices vector. + /// In presence of unstable particles, no 1:1 correspondence. For + /// non-reconstructed stable MC particles, a "-1" entry will be inserted. If + /// 1:1 length correspondence is required, set require_stable to false. + ROOT::VecOps::RVec selRP_indices_matched_to_list( + const ROOT::VecOps::RVec &mcParticles_indices, + const ROOT::VecOps::RVec &recind, + const ROOT::VecOps::RVec &mcind, + const ROOT::VecOps::RVec &mc, + bool require_stable = true); /// return the index of the MC particle that is associated to a given track (via the track-reco association) int getTrack2MC_index ( int track_index, diff --git a/analyzers/dataframe/FCCAnalyses/ReconstructedParticle2Track.h b/analyzers/dataframe/FCCAnalyses/ReconstructedParticle2Track.h index 44f898c220f..bd862994810 100644 --- a/analyzers/dataframe/FCCAnalyses/ReconstructedParticle2Track.h +++ b/analyzers/dataframe/FCCAnalyses/ReconstructedParticle2Track.h @@ -26,151 +26,190 @@ namespace edm4hep { namespace FCCAnalyses{ -namespace ReconstructedParticle2Track{ - - /// Return the momentum of a track to a reconstructed particle - ROOT::VecOps::RVec getRP2TRK_mom (ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks); - - /// Return the charge of a track to a reconstructed particle - ROOT::VecOps::RVec getRP2TRK_charge(ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks); - - //compute the magnetic field Bz - ROOT::VecOps::RVec getRP2TRK_Bz(const ROOT::VecOps::RVec& rps, - const ROOT::VecOps::RVec& tracks); //here computed for all particles passed - - float Bz(const ROOT::VecOps::RVec& rps, - const ROOT::VecOps::RVec& tracks); //here only computed for the first charged particle encountered - - ROOT::VecOps::RVec XPtoPar_dxy(const ROOT::VecOps::RVec& in, - const ROOT::VecOps::RVec& tracks, - const TLorentzVector& V, // primary vertex - const float& Bz); - - ROOT::VecOps::RVec XPtoPar_dz(const ROOT::VecOps::RVec& in, - const ROOT::VecOps::RVec& tracks, - const TLorentzVector& V, // primary vertex - const float& Bz); - - ROOT::VecOps::RVec XPtoPar_phi(const ROOT::VecOps::RVec& in, - const ROOT::VecOps::RVec& tracks, - const TLorentzVector& V, // primary vertex - const float& Bz); - - ROOT::VecOps::RVec XPtoPar_C(const ROOT::VecOps::RVec& in, - const ROOT::VecOps::RVec& tracks, - const float& Bz); - - ROOT::VecOps::RVec XPtoPar_ct(const ROOT::VecOps::RVec& in, - const ROOT::VecOps::RVec& tracks, - const float& Bz); - - /// Return the D0 of a track to a reconstructed particle - ROOT::VecOps::RVec getRP2TRK_D0 (ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks); - - /// Return the Z0 of a track to a reconstructed particle - ROOT::VecOps::RVec getRP2TRK_Z0 (ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks); - - /// Return the Phi of a track to a reconstructed particle - ROOT::VecOps::RVec getRP2TRK_phi (ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks); - - /// Return the omega of a track to a reconstructed particle - ROOT::VecOps::RVec getRP2TRK_omega (ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks); - - /// Return the tanLambda of a track to a reconstructed particle - ROOT::VecOps::RVec getRP2TRK_tanLambda (ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks); - - /// Return the D0 significance of a track to a reconstructed particle - ROOT::VecOps::RVec getRP2TRK_D0_sig (ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks); - - /// Return the Z0 significance of a track to a reconstructed particle - ROOT::VecOps::RVec getRP2TRK_Z0_sig (ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks); - - - /// Return the variance (not the sigma) of the the D0 of a track to a reconstructed particle - ROOT::VecOps::RVec getRP2TRK_D0_cov (ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks); - - /// Return the variance (not the sigma) of the the Z0 of a track to a reconstructed particle - ROOT::VecOps::RVec getRP2TRK_Z0_cov (ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks); - - /// Return the variance (not the sigma) of the the Phi of a track to a reconstructed particle - ROOT::VecOps::RVec getRP2TRK_phi_cov (ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks); - - /// Return the variance (not the sigma) of the omega of a track to a reconstructed particle - ROOT::VecOps::RVec getRP2TRK_omega_cov (ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks); - - /// Return the variance (not the sigma) of the tanLambda of a track to a reconstructed particle - ROOT::VecOps::RVec getRP2TRK_tanLambda_cov (ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks); - - /// Return the off-diag term (d0, phi0) of the covariance matrix - ROOT::VecOps::RVec getRP2TRK_d0_phi0_cov(ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks); - - /// Return the off-diag term (d0, omega) of the covariance matrix - ROOT::VecOps::RVec getRP2TRK_d0_omega_cov(ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks); - - /// Return the off-diag term (d0,z0) of the covariance matrix - ROOT::VecOps::RVec getRP2TRK_d0_z0_cov(ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks); - - /// Return the off-diag term (d0,tanlambda) of the covariance matrix - ROOT::VecOps::RVec getRP2TRK_d0_tanlambda_cov(ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks); - - /// Return the off-diag term (phi0,omega) of the covariance matrix - ROOT::VecOps::RVec getRP2TRK_phi0_omega_cov(ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks); - - /// Return the off-diag term (phi0,z0) of the covariance matrix - ROOT::VecOps::RVec getRP2TRK_phi0_z0_cov(ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks); - - /// Return the off-diag term (phi0,tanlambda) of the covariance matrix - ROOT::VecOps::RVec getRP2TRK_phi0_tanlambda_cov(ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks) ; - - /// Return the off-diag term (omega,z0) of the covariance matrix - ROOT::VecOps::RVec getRP2TRK_omega_z0_cov(ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks) ; - - /// Return the off-diag term (omega,tanlambda) of the covariance matrix - ROOT::VecOps::RVec getRP2TRK_omega_tanlambda_cov(ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks) ; - - /// Return the off-diag term (z0,tanlambda) of the covariance matrix - ROOT::VecOps::RVec getRP2TRK_z0_tanlambda_cov(ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks); - - - /// Return the tracks associated to reco'ed particles - ROOT::VecOps::RVec getRP2TRK( ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks ) ; - - /// Return the reco indices of particles that have tracks - ROOT::VecOps::RVec get_recoindTRK( ROOT::VecOps::RVec in, - ROOT::VecOps::RVec tracks ) ; - - /// Return the size of a collection of TrackStates - int getTK_n(ROOT::VecOps::RVec x) ; - - /// Return if a Reco particle have an associated track - ROOT::VecOps::RVec hasTRK( ROOT::VecOps::RVec in ) ; - -}//end NS ReconstructedParticle2Track +namespace ReconstructedParticle2Track { + +/// Return the momentum of a track to a reconstructed particle +ROOT::VecOps::RVec +getRP2TRK_mom(ROOT::VecOps::RVec in, + ROOT::VecOps::RVec tracks); + +/// Return the charge of a track to a reconstructed particle +ROOT::VecOps::RVec +getRP2TRK_charge(ROOT::VecOps::RVec in, + ROOT::VecOps::RVec tracks); + +// compute the magnetic field Bz +ROOT::VecOps::RVec +getRP2TRK_Bz(const ROOT::VecOps::RVec &rps, + const ROOT::VecOps::RVec + &tracks); // here computed for all particles passed + +float Bz(const ROOT::VecOps::RVec &rps, + const ROOT::VecOps::RVec + &tracks); // here only computed for the first charged particle + // encountered + +ROOT::VecOps::RVec +XPtoPar_dxy(const ROOT::VecOps::RVec &in, + const ROOT::VecOps::RVec &tracks, + const TLorentzVector &V, // primary vertex + const float &Bz); + +ROOT::VecOps::RVec +XPtoPar_dz(const ROOT::VecOps::RVec &in, + const ROOT::VecOps::RVec &tracks, + const TLorentzVector &V, // primary vertex + const float &Bz); + +ROOT::VecOps::RVec +XPtoPar_phi(const ROOT::VecOps::RVec &in, + const ROOT::VecOps::RVec &tracks, + const TLorentzVector &V, // primary vertex + const float &Bz); + +ROOT::VecOps::RVec +XPtoPar_C(const ROOT::VecOps::RVec &in, + const ROOT::VecOps::RVec &tracks, + const float &Bz); + +ROOT::VecOps::RVec +XPtoPar_ct(const ROOT::VecOps::RVec &in, + const ROOT::VecOps::RVec &tracks, + const float &Bz); + +/// Return the D0 of a track to a reconstructed particle +ROOT::VecOps::RVec +getRP2TRK_D0(ROOT::VecOps::RVec in, + ROOT::VecOps::RVec tracks); + +/// Return the Z0 of a track to a reconstructed particle +ROOT::VecOps::RVec +getRP2TRK_Z0(ROOT::VecOps::RVec in, + ROOT::VecOps::RVec tracks); + +/// Return the Phi of a track to a reconstructed particle +ROOT::VecOps::RVec +getRP2TRK_phi(ROOT::VecOps::RVec in, + ROOT::VecOps::RVec tracks); + +/// Return the omega of a track to a reconstructed particle +ROOT::VecOps::RVec +getRP2TRK_omega(ROOT::VecOps::RVec in, + ROOT::VecOps::RVec tracks); + +/// Return the tanLambda of a track to a reconstructed particle +ROOT::VecOps::RVec +getRP2TRK_tanLambda(ROOT::VecOps::RVec in, + ROOT::VecOps::RVec tracks); + +/// Return the D0 significance of a track to a reconstructed particle +ROOT::VecOps::RVec +getRP2TRK_D0_sig(ROOT::VecOps::RVec in, + ROOT::VecOps::RVec tracks); + +/// Return the Z0 significance of a track to a reconstructed particle +ROOT::VecOps::RVec +getRP2TRK_Z0_sig(ROOT::VecOps::RVec in, + ROOT::VecOps::RVec tracks); + +/// Return the variance (not the sigma) of the the D0 of a track to a +/// reconstructed particle +ROOT::VecOps::RVec +getRP2TRK_D0_cov(ROOT::VecOps::RVec in, + ROOT::VecOps::RVec tracks); + +/// Return the variance (not the sigma) of the the Z0 of a track to a +/// reconstructed particle +ROOT::VecOps::RVec +getRP2TRK_Z0_cov(ROOT::VecOps::RVec in, + ROOT::VecOps::RVec tracks); + +/// Return the variance (not the sigma) of the the Phi of a track to a +/// reconstructed particle +ROOT::VecOps::RVec +getRP2TRK_phi_cov(ROOT::VecOps::RVec in, + ROOT::VecOps::RVec tracks); + +/// Return the variance (not the sigma) of the omega of a track to a +/// reconstructed particle +ROOT::VecOps::RVec +getRP2TRK_omega_cov(ROOT::VecOps::RVec in, + ROOT::VecOps::RVec tracks); + +/// Return the variance (not the sigma) of the tanLambda of a track to a +/// reconstructed particle +ROOT::VecOps::RVec getRP2TRK_tanLambda_cov( + ROOT::VecOps::RVec in, + ROOT::VecOps::RVec tracks); + +/// Return the off-diag term (d0, phi0) of the covariance matrix +ROOT::VecOps::RVec +getRP2TRK_d0_phi0_cov(ROOT::VecOps::RVec in, + ROOT::VecOps::RVec tracks); + +/// Return the off-diag term (d0, omega) of the covariance matrix +ROOT::VecOps::RVec getRP2TRK_d0_omega_cov( + ROOT::VecOps::RVec in, + ROOT::VecOps::RVec tracks); + +/// Return the off-diag term (d0,z0) of the covariance matrix +ROOT::VecOps::RVec +getRP2TRK_d0_z0_cov(ROOT::VecOps::RVec in, + ROOT::VecOps::RVec tracks); + +/// Return the off-diag term (d0,tanlambda) of the covariance matrix +ROOT::VecOps::RVec getRP2TRK_d0_tanlambda_cov( + ROOT::VecOps::RVec in, + ROOT::VecOps::RVec tracks); + +/// Return the off-diag term (phi0,omega) of the covariance matrix +ROOT::VecOps::RVec getRP2TRK_phi0_omega_cov( + ROOT::VecOps::RVec in, + ROOT::VecOps::RVec tracks); + +/// Return the off-diag term (phi0,z0) of the covariance matrix +ROOT::VecOps::RVec +getRP2TRK_phi0_z0_cov(ROOT::VecOps::RVec in, + ROOT::VecOps::RVec tracks); + +/// Return the off-diag term (phi0,tanlambda) of the covariance matrix +ROOT::VecOps::RVec getRP2TRK_phi0_tanlambda_cov( + ROOT::VecOps::RVec in, + ROOT::VecOps::RVec tracks); + +/// Return the off-diag term (omega,z0) of the covariance matrix +ROOT::VecOps::RVec getRP2TRK_omega_z0_cov( + ROOT::VecOps::RVec in, + ROOT::VecOps::RVec tracks); + +/// Return the off-diag term (omega,tanlambda) of the covariance matrix +ROOT::VecOps::RVec getRP2TRK_omega_tanlambda_cov( + ROOT::VecOps::RVec in, + ROOT::VecOps::RVec tracks); + +/// Return the off-diag term (z0,tanlambda) of the covariance matrix +ROOT::VecOps::RVec getRP2TRK_z0_tanlambda_cov( + ROOT::VecOps::RVec in, + ROOT::VecOps::RVec tracks); + +/// Return the tracks associated to reco'ed particles +ROOT::VecOps::RVec +getRP2TRK(ROOT::VecOps::RVec in, + ROOT::VecOps::RVec tracks); + +/// Return the reco indices of particles that have tracks +ROOT::VecOps::RVec +get_recoindTRK(ROOT::VecOps::RVec in, + ROOT::VecOps::RVec tracks); + +/// Return the size of a collection of TrackStates +int getTK_n(ROOT::VecOps::RVec x); + +/// Return if a Reco particle have an associated track +ROOT::VecOps::RVec +hasTRK(ROOT::VecOps::RVec in); + +} // namespace ReconstructedParticle2Track }//end NS FCCAnalyses #endif diff --git a/analyzers/dataframe/FCCAnalyses/Utils.h b/analyzers/dataframe/FCCAnalyses/Utils.h index 78c5616d0d7..3739545e9bd 100644 --- a/analyzers/dataframe/FCCAnalyses/Utils.h +++ b/analyzers/dataframe/FCCAnalyses/Utils.h @@ -1,13 +1,188 @@ #ifndef UTILS_ANALYZERS_H #define UTILS_ANALYZERS_H +#include "ROOT/RVec.hxx" +#include #include namespace FCCAnalyses { namespace Utils { + /// @brief Helper struct to select entries matching a certain predicate. + /// When instantiated with a function that gives a "yes / no" decision + /// based on an input object, will support three functionalities: + /// - select accepted objects from a list of inputs, return as list + /// - select indices of selected objects from list of inputs and indices to + /// study, + /// return index ist + /// - select indicises of objects passing for a set of input lists (nested + /// selection) This allows to consistently apply the same selection in + /// multiple ways, and to easily write new cuts without having to write a lot + /// of boilerplate code. + /// @tparam thingToSelect: The type of object our decision will be based on + /// For example: edm4hep::MCParticleData (for selecting MC particles). + /// + /// One pattern is to create a selector class inheriting from + /// `selByPredicate`, where the selection function is specified within the + /// constructor. Constructor arguments can be passed via lambda capture or + /// `std::bind` to parametrise cuts where needed. This can look like + /// @code + /// struct myPtCut : selByPredicate{/// + /// myPtCut(double ptCutVal): selByPredicate( + /// [ptCutVal](const edm4hep::MCParticleData & mc){ + /// return mc.pt() > ptCutVal; + /// }){ + /// } + /// }; + /// @endcode + /// Or, assuming you have an existing selection function + /// `doesThisPass(const edm4hep::MCParticleData &, double, int)` + /// with some customizable cut parameters + /// @code + /// struct myDoesThisPass : selByPredicate{/// + /// myPtCut(double par1, int par2): + /// selByPredicate( + /// std::bind(doesThisPass,std::placeholders::_1, par1, par2)){ + /// } + /// }; + /// @endcode + /// Examples for practical usage can be found in `MCParticle.h`, see for + /// example #sel_pt. In the python steering, you would then typically + /// instantiate the derived cut classes, for example + /// @code + /// df = df.Define("particles_passing","myPtCut(10)(particles_to_check)") + /// @endcode + /// or (for indices) + /// @code + /// df = + /// df.Define("particle_indices_passing","myPtCut(10)(particle_indices_to_check, + /// allParticles)") + /// @endcode + template struct selByPredicate { + /// @brief constructor - instantiate by passing a selection function. + /// For implementing selections with configurable parameters + /// (e.g. cut values), look at lambda captures or `std::bind`. + selByPredicate(std::function thePredicate) + : m_predicate(thePredicate) {} + + /// @brief Select passing objects from a list of candidates, + /// return by (deep-)copy + /// @param in: List of candidates to run the selection on + /// @return A list of candidates for which `m_predicate` evaluates + /// as `true`. + ROOT::VecOps::RVec + operator()(const ROOT::VecOps::RVec &in) { + ROOT::VecOps::RVec result; + result.reserve(in.size()); + for (auto &p : in) { + if (m_predicate(p)) + result.emplace_back(p); + } + return result; + } + + /// @brief Select passing objects from a list of candidates, + /// return by index + /// @param indices: List of indices of candidates within `in` + // to consider for selection + /// @param in: List of candidates to which indices in `indices` refer. + /// @return A list of indices of candidates within the `in` + // vector for which `m_predicate` evaluates as `true`. + ROOT::VecOps::RVec + operator()(const ROOT::VecOps::RVec &indices, + const ROOT::VecOps::RVec &in) { + ROOT::VecOps::RVec result; + result.reserve(in.size()); + for (int index : indices) { + if (index < 0 || index >= in.size()) + continue; + if (m_predicate(in[index])) + result.emplace_back(index); + } + return result; + } + + /// @brief Select (several) lists of passing objects from (several) lists of + /// candidates, return by lists of indices. Used to select on several + /// (similar) lists at once. For example: Track selection within (many) + /// vertex candidates. + /// @param setsOfIndices: Multiple lists of indices referring to `in` + /// @param in: List of candidates to which indices in `setsOfIndices` refer. + /// @return A list of index lists of similar size as `setsOfIndices`. + /// For each element, the list contains the output of `operator()` + /// run on the indices of the given element of `setsOfIndices`. + ROOT::VecOps::RVec> + operator()(const ROOT::VecOps::RVec> &setsOfIndices, + const ROOT::VecOps::RVec &in) { + ROOT::VecOps::RVec> result(setsOfIndices.size()); + for (int elem = 0; elem < setsOfIndices.size(); ++elem) { + result[elem] = this->operator()(setsOfIndices[elem], in); + } + return result; + } + + private: + /// Stored selection function. + std::function m_predicate; + }; + template inline auto getsize( T& vec){ return vec.size();}; - template inline ROOT::VecOps::RVec > as_vector(const ROOT::VecOps::RVec& in){return ROOT::VecOps::RVec >(1, in);}; + template + inline ROOT::VecOps::RVec> + as_vector(const ROOT::VecOps::RVec &in) { + return ROOT::VecOps::RVec>(1, in); + }; + template + inline ROOT::VecOps::RVec + index_range(const ROOT::VecOps::RVec &in) { + ROOT::VecOps::RVec indices(in.size()); + std::iota(indices.begin(), indices.end(), 0); + return indices; + }; + /// @brief count the number of valid (>=0, < size of collection) indices in + /// an index list + /// @param in: index list + /// @param ref: particle vector to which the indices refer. + /// @return integer count of valid indices + template + inline int count_valid_indices(const ROOT::VecOps::RVec &in, + const ROOT::VecOps::RVec &ref) { + int maxSize = ref.size(); + return std::count_if(in.begin(), in.end(), [maxSize](const int &i) { + return (i >= 0 && i < maxSize); + }); + }; + + // @brief for a given list of indices, returns a list of particle (copies) + /// @param idx : Indices of desired particles within the full set ("in") + /// @param in : Full set of particles + /// @return A vector of particles with the desired indices. + template + inline ROOT::VecOps::RVec sel_byIndex(const ROOT::VecOps::RVec &idx, + const ROOT::VecOps::RVec &in) { + ROOT::VecOps::RVec found; + for (int index : idx) { + if (index < 0 || index >= in.size()) + continue; + found.push_back(in.at(index)); + } + return found; + } + // @brief merge (concatenate) two collections of arbitrary content + /// @param x : first collection - entries will be copied in-order + /// @param y : second collection - entries will be copied in-order after the + /// last element of the first + /// @return A combined collection of size (x.size()+y.size()), containing + /// the content of x followed by that of y + template + inline ROOT::VecOps::RVec merge(const ROOT::VecOps::RVec &x, + const ROOT::VecOps::RVec &y) { + ROOT::VecOps::RVec merged; + merged.reserve(x.size() + y.size()); + merged.insert(merged.end(), x.begin(), x.end()); + merged.insert(merged.end(), y.begin(), y.end()); + return merged; + } } } diff --git a/analyzers/dataframe/FCCAnalyses/VertexFitterSimple.h b/analyzers/dataframe/FCCAnalyses/VertexFitterSimple.h index a3a401aa86b..948a630eaa1 100644 --- a/analyzers/dataframe/FCCAnalyses/VertexFitterSimple.h +++ b/analyzers/dataframe/FCCAnalyses/VertexFitterSimple.h @@ -71,7 +71,6 @@ namespace VertexFitterSimple{ ROOT::VecOps::RVec IsPrimary_forTracks( ROOT::VecOps::RVec allTracks, ROOT::VecOps::RVec primaryTracks ) ; - /* Double_t FastRv(TVectorD p1, TVectorD p2) ; TMatrixDSym RegInv3(TMatrixDSym &Smat0) ; diff --git a/analyzers/dataframe/FCCAnalyses/VertexingUtils.h b/analyzers/dataframe/FCCAnalyses/VertexingUtils.h index 96721fde586..fdbbb456340 100644 --- a/analyzers/dataframe/FCCAnalyses/VertexingUtils.h +++ b/analyzers/dataframe/FCCAnalyses/VertexingUtils.h @@ -107,7 +107,30 @@ namespace VertexingUtils{ /// Retrieve the indices of the tracks fitted to that vertex, but now in the collection of RecoParticles ROOT::VecOps::RVec get_VertexRecoParticlesInd( FCCAnalysesVertex TheVertex, const ROOT::VecOps::RVec& reco ); - + + /// Retrieve the indices of the tracks fitted to a vector of vertices, but now + /// in the collection of RecoParticles + ROOT::VecOps::RVec get_VerticesRecoParticlesInd( + ROOT::VecOps::RVec vertices, + const ROOT::VecOps::RVec &reco); + + /// @brief Find (by index) a vertex that is made of reco particles from a + /// passed list. + /// @param vertices: List of possible vertices to match + /// @param recoParticleIndices: Indices of the reco particles to find in the + /// vertex + /// @param reco: list of all reco particles (ReconstructedParticles) + /// @param require_all: If set, require one vertex to contain all + /// recoParticleIndices (no 'missed' tracks). Else only require that all + /// tracks in the vertex come from the vector (but allow for unused tracks) + /// @return: Index within the list of the (first) vertex fulfilling the + /// criteria. If none are found, return -1 + int getVertex_matching_recoParticles( + const ROOT::VecOps::RVec &vertices, + const ROOT::VecOps::RVec &recoParticleIndices, + const ROOT::VecOps::RVec &reco, + bool require_all = false); + /// Return the number of tracks in a given track collection int get_nTracks(ROOT::VecOps::RVec tracks); diff --git a/analyzers/dataframe/FCCAnalyses/myUtils.h b/analyzers/dataframe/FCCAnalyses/myUtils.h index 2b8f3a6ded2..fbff7eaffb1 100644 --- a/analyzers/dataframe/FCCAnalyses/myUtils.h +++ b/analyzers/dataframe/FCCAnalyses/myUtils.h @@ -1,8 +1,8 @@ #ifndef MYUTILS_ANALYZERS_H #define MYUTILS_ANALYZERS_H #include "ROOT/RVec.hxx" -#include "edm4hep/ReconstructedParticleData.h" #include "edm4hep/MCParticleData.h" +#include "edm4hep/ReconstructedParticleData.h" #include "edm4hep/TrackState.h" #include "edm4hep/VertexData.h" @@ -77,7 +77,6 @@ namespace myUtils{ ROOT::VecOps::RVec operator() (ROOT::VecOps::RVec recop); }; - ROOT::VecOps::RVec get_pseudotrack(ROOT::VecOps::RVec vertex, ROOT::VecOps::RVec recop); diff --git a/analyzers/dataframe/src/MCParticle.cc b/analyzers/dataframe/src/MCParticle.cc index 009d0d4c255..0c75dff8653 100644 --- a/analyzers/dataframe/src/MCParticle.cc +++ b/analyzers/dataframe/src/MCParticle.cc @@ -8,36 +8,43 @@ namespace FCCAnalyses{ namespace MCParticle{ -sel_genStatus::sel_genStatus(int arg_status) : m_status(arg_status) {}; -ROOT::VecOps::RVec sel_genStatus::operator() (ROOT::VecOps::RVec in) { - ROOT::VecOps::RVec result; - result.reserve(in.size()); - for (size_t i = 0; i < in.size(); ++i) { - auto & p = in[i]; - if (p.generatorStatus == m_status) { - result.emplace_back(p); - } - } - return result; -} - -sel_pdgID::sel_pdgID(int arg_pdg, bool arg_chargeconjugate) : m_pdg(arg_pdg), m_chargeconjugate( arg_chargeconjugate ) {}; -ROOT::VecOps::RVec sel_pdgID::operator() (ROOT::VecOps::RVec in) { - ROOT::VecOps::RVec result; - result.reserve(in.size()); - for (size_t i = 0; i < in.size(); ++i) { - auto & p = in[i]; - if ( m_chargeconjugate ) { - if ( std::abs( p.PDG ) == std::abs( m_pdg) ) result.emplace_back(p); - } - else { - if ( p.PDG == m_pdg ) result.emplace_back(p); - } - } - return result; -} - - +sel_pt::sel_pt(float arg_min_pt) + : MCParticleSelection( + [arg_min_pt](const edm4hep::MCParticleData &p) -> bool { + return (p.momentum.x * p.momentum.x + p.momentum.y * p.momentum.y > + arg_min_pt * arg_min_pt); + }) {} + +sel_eta::sel_eta(float arg_max_eta) + : MCParticleSelection( + [arg_max_eta](const edm4hep::MCParticleData &p) -> bool { + ROOT::Math::PxPyPzM4D vec(p.momentum.x, p.momentum.y, p.momentum.z, + p.mass); + return std::abs(vec.Eta()) < std::abs(arg_max_eta); + }) {} + +sel_genStatus::sel_genStatus(int arg_status) + : MCParticleSelection( + [arg_status](const edm4hep::MCParticleData &p) -> bool { + return p.generatorStatus == arg_status; + }) {} + +sel_pdgID::sel_pdgID(int arg_pdg, bool arg_chargeconjugate) + : MCParticleSelection([arg_pdg, arg_chargeconjugate]( + const edm4hep::MCParticleData &p) -> bool { + if (arg_chargeconjugate) + return std::abs(p.PDG) == std::abs(arg_pdg); + else + return p.PDG == arg_pdg; + }) {} + +sel_charged::sel_charged() + : MCParticleSelection([](const edm4hep::MCParticleData &p) -> bool { + // try to avoid floating comp to zero + // and I am sad that some neutral particles carry a -999 + // charge! + return std::abs(p.charge) > 1e-6 && p.charge != -999; + }) {} get_decay::get_decay(int arg_mother, int arg_daughters, bool arg_inf){m_mother=arg_mother; m_daughters=arg_daughters; m_inf=arg_inf;}; bool get_decay::operator() (ROOT::VecOps::RVec in, ROOT::VecOps::RVec ind){ @@ -59,20 +66,6 @@ bool get_decay::operator() (ROOT::VecOps::RVec in, ROO return result; } -sel_pt::sel_pt(float arg_min_pt) : m_min_pt(arg_min_pt) {}; -ROOT::VecOps::RVec sel_pt::operator() (ROOT::VecOps::RVec in) { - ROOT::VecOps::RVec result; - result.reserve(in.size()); - for (size_t i = 0; i < in.size(); ++i) { - auto & p = in[i]; - if (std::sqrt(std::pow(p.momentum.x,2) + std::pow(p.momentum.y,2)) > m_min_pt) { - result.emplace_back(p); - } - } - return result; -} - - filter_pdgID::filter_pdgID(int arg_pdgid, bool arg_abs){m_pdgid = arg_pdgid; m_abs = arg_abs;}; bool filter_pdgID::operator() (ROOT::VecOps::RVec in) { for (size_t i = 0; i < in.size(); ++i) { @@ -502,6 +495,9 @@ std::vector get_list_of_stable_particles_from_decay( int i, ROOT::VecOps::R //for (int idaughter = d1; idaughter <= d2; idaughter++) { for (int id = db; id < de; id++) { int idaughter = ind[ id ]; + // prevent endless loop in case of looping MC record + if (idaughter == i) + continue; std::vector rr = get_list_of_stable_particles_from_decay( idaughter, in, ind) ; res.insert( res.end(), rr.begin(), rr.end() ); } @@ -540,6 +536,18 @@ std::vector get_list_of_particles_from_decay(int i, ROOT::VecOps::RVec> +get_lists_of_stable_particles_from_decays( + ROOT::VecOps::RVec i, ROOT::VecOps::RVec in, + ROOT::VecOps::RVec ind) { + + ROOT::VecOps::RVec> ret; + ret.reserve(i.size()); + for (int ix : i) { + ret.push_back(get_list_of_stable_particles_from_decay(ix, in, ind)); + } + return ret; +} // ---------------------------------------------------------------------------------------------------------------------------------- @@ -555,11 +563,6 @@ std::vector list_of_particles_from_decay(int i, ROOT::VecOps::RVec get_indices_MotherByIndex ( int imother, diff --git a/analyzers/dataframe/src/ReconstructedParticle2MC.cc b/analyzers/dataframe/src/ReconstructedParticle2MC.cc index 45abccdd4e8..6a6d510082a 100644 --- a/analyzers/dataframe/src/ReconstructedParticle2MC.cc +++ b/analyzers/dataframe/src/ReconstructedParticle2MC.cc @@ -337,26 +337,23 @@ selRP_ChargedHadrons (ROOT::VecOps::RVec recind, // ------------------------------------------------------------------------------------------------- -// -- select RecoParticles associated with a list of MC particles (passed by their index in the Particle block) +// -- select indices of RecoParticles associated with a list of MC particles +// (passed by their index in the Particle block) -ROOT::VecOps::RVec -selRP_matched_to_list( ROOT::VecOps::RVec mcParticles_indices, - ROOT::VecOps::RVec recind, - ROOT::VecOps::RVec mcind, - ROOT::VecOps::RVec reco, - ROOT::VecOps::RVec mc) { - - ROOT::VecOps::RVec results; - - edm4hep::ReconstructedParticleData dummy; - dummy.energy = -9999; - dummy.tracks_begin = -9999 ; +ROOT::VecOps::RVec selRP_indices_matched_to_list( + const ROOT::VecOps::RVec &mcParticles_indices, + const ROOT::VecOps::RVec &recind, const ROOT::VecOps::RVec &mcind, + const ROOT::VecOps::RVec &mc, + bool require_stable) { + ROOT::VecOps::RVec results; for ( auto & idx: mcParticles_indices ) { // exclude unstable particles - e.g. the list may contain the index of // the mother - if ( mc.at(idx).generatorStatus != 1 ) continue ; + // MG: Should we push_back -1 here for consistency? + if (require_stable && mc.at(idx).generatorStatus != 1) + continue; // is this MC particle associated with a Reco particle : bool found = false; @@ -365,23 +362,46 @@ selRP_matched_to_list( ROOT::VecOps::RVec mcParticles_indices, int mc_idx = mcind.at(i); if ( mc_idx == idx ) { found = true; - results.push_back( reco.at( reco_idx ) ); + results.push_back(reco_idx); break; } } // no Reco particle has been found for idx: add a dummy particle such that // one preserves the mapping with the input list - if ( ! found) results.push_back( dummy ); - - + // Note: This is inconsistent, the mapping will be + // lost whenever there is an unstable particle in the input list. + // Will keep current behaviour for backward compatibility... + if (!found) + results.push_back(-1); } // loop over the indices in the list return results; - } +// -- select RecoParticles associated with a list of MC particles (passed by +// their index in the Particle block) - +ROOT::VecOps::RVec selRP_matched_to_list( + const ROOT::VecOps::RVec &mcParticles_indices, + const ROOT::VecOps::RVec &recind, const ROOT::VecOps::RVec &mcind, + const ROOT::VecOps::RVec &reco, + const ROOT::VecOps::RVec &mc, + bool require_stable) { + edm4hep::ReconstructedParticleData dummy; + dummy.energy = -9999; + dummy.tracks_begin = -9999; + + ROOT::VecOps::RVec results; + ROOT::VecOps::RVec indices = selRP_indices_matched_to_list( + mcParticles_indices, recind, mcind, mc, require_stable); + for (int reco_idx : indices) { + if (reco_idx < 0) + results.push_back(dummy); + else + results.push_back(reco.at(reco_idx)); + } + return results; +} // ------------------------------------------------------------------------------------------------- diff --git a/analyzers/dataframe/src/VertexingUtils.cc b/analyzers/dataframe/src/VertexingUtils.cc index f60ccfba7e9..feac0b5af89 100644 --- a/analyzers/dataframe/src/VertexingUtils.cc +++ b/analyzers/dataframe/src/VertexingUtils.cc @@ -1,7 +1,7 @@ #include "FCCAnalyses/VertexingUtils.h" #include "FCCAnalyses/VertexFitterSimple.h" - #include "TrkUtil.h" // from delphes +#include namespace FCCAnalyses { @@ -369,6 +369,51 @@ ROOT::VecOps::RVec get_VertexRecoParticlesInd( return result; } +ROOT::VecOps::RVec get_VerticesRecoParticlesInd( + ROOT::VecOps::RVec vertices, + const ROOT::VecOps::RVec &reco) { + + ROOT::VecOps::RVec result; + for (int j = 0; j < vertices.size(); ++j) { + ROOT::VecOps::RVec indices_tracks = vertices[j].reco_ind; + for (int i = 0; i < indices_tracks.size(); i++) { + int tk_index = indices_tracks[i]; + for (int j = 0; j < reco.size(); j++) { + auto &p = reco[j]; + if (p.tracks_begin == p.tracks_end) + continue; + if (p.tracks_begin == tk_index) { + result.push_back(j); + break; + } + } + } + } + + return result; +} + +int getVertex_matching_recoParticles( + const ROOT::VecOps::RVec &vertices, + const ROOT::VecOps::RVec &recoParticleIndices, + const ROOT::VecOps::RVec &reco, + bool require_all) { + std::set indicesWeWant; + indicesWeWant.insert(recoParticleIndices.begin(), recoParticleIndices.end()); + // correct for "-1" representing missed tracks in the recoParticleIndices + int correctMissing = indicesWeWant.count(-1); + for (int iVX = 0; iVX < vertices.size(); ++iVX) { + auto vxParticleIndices = get_VertexRecoParticlesInd(vertices[iVX], reco); + int nFound = std::count_if( + vxParticleIndices.begin(), vxParticleIndices.end(), + [&](int recoIndex) { return indicesWeWant.count(recoIndex); }); + if (require_all && nFound == indicesWeWant.size() - correctMissing || + nFound == vxParticleIndices.size()) + return iVX; + } + return -1; +} + TVectorD ParToACTS(TVectorD Par) { TVectorD pACTS(6); // Return vector