Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-12-13 10:26:15

0001 #ifndef RIVET_PARTICLEBASEUTILS_HH
0002 #define RIVET_PARTICLEBASEUTILS_HH
0003 
0004 #include "Rivet/Tools/Utils.hh"
0005 #include "Rivet/ParticleBase.hh"
0006 
0007 namespace Rivet {
0008 
0009 
0010   /// @defgroup particlebaseutils Functions for Particles and Jets
0011   /// @{
0012 
0013   /// @defgroup particlebasetutils_pb2bool ParticleBase classifier -> bool functors
0014   /// @todo Move to FourMomentum functions
0015   ///
0016   /// To be passed to any() or all() e.g. any(jets, DeltaRLess(electron, 0.4))
0017   /// @{
0018 
0019   /// std::function instantiation for functors taking a ParticleBase and returning a bool
0020   using ParticleBaseSelector = function<bool(const ParticleBase&)>;
0021   /// std::function instantiation for functors taking two ParticleBase and returning a bool
0022   using ParticleBaseSorter = function<bool(const ParticleBase&, const ParticleBase&)>;
0023 
0024 
0025   /// Base type for Particle -> bool functors
0026   struct BoolParticleBaseFunctor {
0027     virtual bool operator()(const ParticleBase& p) const = 0;
0028     virtual ~BoolParticleBaseFunctor() {}
0029   };
0030 
0031 
0032   /// Transverse momentum greater-than functor
0033   struct PtGtr : public BoolParticleBaseFunctor {
0034     PtGtr(double pt) : ptcut(pt) { }
0035     PtGtr(const FourMomentum& p) : ptcut(p.pT()) { }
0036     bool operator()(const ParticleBase& p) const { return p.pT() > ptcut; }
0037     double ptcut;
0038   };
0039   using pTGtr = PtGtr;
0040   using ptGtr = PtGtr;
0041 
0042   /// Transverse momentum less-than functor
0043   struct PtLess : public BoolParticleBaseFunctor {
0044     PtLess(const FourMomentum& p) : ptcut(p.pT()) { }
0045     PtLess(double pt) : ptcut(pt) { }
0046     bool operator()(const ParticleBase& p) const { return p.pT() < ptcut; }
0047     double ptcut;
0048   };
0049   using pTLess = PtLess;
0050   using ptLess = PtLess;
0051 
0052   /// Transverse momentum in-range functor
0053   struct PtInRange : public BoolParticleBaseFunctor {
0054     PtInRange(pair<double, double> ptcuts) : ptcut(ptcuts) { }
0055     PtInRange(double ptlow, double pthigh) : PtInRange(make_pair(ptlow, pthigh)) { }
0056     PtInRange(const FourMomentum& p1, const FourMomentum& p2) : PtInRange(p1.pT(), p2.pT()) { }
0057     bool operator()(const ParticleBase& p) const { return p.pT() >= ptcut.first && p.pT() < ptcut.second; }
0058     pair<double,double> ptcut;
0059   };
0060   using pTInRange = PtInRange;
0061   using ptInRange = PtInRange;
0062 
0063 
0064   /// Pseudorapidity greater-than functor
0065   struct EtaGtr : public BoolParticleBaseFunctor {
0066     EtaGtr(double eta) : etacut(eta) { }
0067     EtaGtr(const FourMomentum& p) : etacut(p.eta()) { }
0068     bool operator()(const ParticleBase& p) const { return p.eta() > etacut; }
0069     double etacut;
0070   };
0071   using etaGtr = EtaGtr;
0072 
0073   /// Pseudorapidity less-than functor
0074   struct EtaLess : public BoolParticleBaseFunctor {
0075     EtaLess(double eta) : etacut(eta) { }
0076     EtaLess(const FourMomentum& p) : etacut(p.eta()) { }
0077     bool operator()(const ParticleBase& p) const { return p.eta() < etacut; }
0078     double etacut;
0079   };
0080   using etaLess = EtaLess;
0081 
0082   /// Pseudorapidity in-range functor
0083   struct EtaInRange : public BoolParticleBaseFunctor {
0084     EtaInRange(pair<double, double> etacuts) : etacut(etacuts) { }
0085     EtaInRange(double etalow, double etahigh) : EtaInRange(make_pair(etalow, etahigh)) { }
0086     EtaInRange(const FourMomentum& p1, const FourMomentum& p2) : EtaInRange(p1.eta(), p2.eta()) { }
0087     bool operator()(const ParticleBase& p) const { return p.eta() >= etacut.first && p.eta() < etacut.second; }
0088     pair<double,double> etacut;
0089   };
0090   using etaInRange = EtaInRange;
0091 
0092 
0093   /// Abs pseudorapidity greater-than functor
0094   struct AbsEtaGtr : public BoolParticleBaseFunctor {
0095     AbsEtaGtr(double abseta) : absetacut(abseta) { }
0096     AbsEtaGtr(const FourMomentum& p) : absetacut(p.abseta()) { }
0097     bool operator()(const ParticleBase& p) const { return p.abseta() > absetacut; }
0098     double absetacut;
0099   };
0100   using absEtaGtr = AbsEtaGtr;
0101   using absetaGtr = AbsEtaGtr;
0102 
0103   /// Abs pseudorapidity momentum less-than functor
0104   struct AbsEtaLess : public BoolParticleBaseFunctor {
0105     AbsEtaLess(double abseta) : absetacut(abseta) { }
0106     AbsEtaLess(const FourMomentum& p) : absetacut(p.abseta()) { }
0107     bool operator()(const ParticleBase& p) const { return p.abseta() < absetacut; }
0108     double absetacut;
0109   };
0110   using absEtaLess = AbsEtaLess;
0111   using absetaLess = AbsEtaLess;
0112 
0113   /// Abs pseudorapidity in-range functor
0114   struct AbsEtaInRange : public BoolParticleBaseFunctor {
0115     AbsEtaInRange(const pair<double, double>& absetacuts) : absetacut(absetacuts) { }
0116     AbsEtaInRange(double absetalow, double absetahigh) : AbsEtaInRange(make_pair(absetalow, absetahigh)) { }
0117     AbsEtaInRange(const FourMomentum& p1, const FourMomentum& p2) : AbsEtaInRange(p1.abseta(), p2.abseta()) { }
0118     bool operator()(const ParticleBase& p) const { return p.abseta() >= absetacut.first && p.abseta() < absetacut.second; }
0119     pair<double,double> absetacut;
0120   };
0121   using absEtaInRange = AbsEtaInRange;
0122   using absetaInRange = AbsEtaInRange;
0123 
0124 
0125   /// Rapidity greater-than functor
0126   struct RapGtr : public BoolParticleBaseFunctor {
0127     RapGtr(double rap) : rapcut(rap) { }
0128     RapGtr(const FourMomentum& p) : rapcut(p.rap()) { }
0129     bool operator()(const ParticleBase& p) const { return p.rap() > rapcut; }
0130     double rapcut;
0131   };
0132   using rapGtr = RapGtr;
0133 
0134   /// Rapidity momentum less-than functor
0135   struct RapLess : public BoolParticleBaseFunctor {
0136     RapLess(double rap) : rapcut(rap) { }
0137     RapLess(const FourMomentum& p) : rapcut(p.rap()) { }
0138     bool operator()(const ParticleBase& p) const { return p.rap() < rapcut; }
0139     double rapcut;
0140   };
0141   using rapLess = RapLess;
0142 
0143   /// Rapidity in-range functor
0144   struct RapInRange : public BoolParticleBaseFunctor {
0145     RapInRange(const pair<double, double>& rapcuts) : rapcut(rapcuts) { }
0146     RapInRange(double raplow, double raphigh) : RapInRange(make_pair(raplow, raphigh)) { }
0147     RapInRange(const FourMomentum& p1, const FourMomentum& p2) : RapInRange(p1.rap(), p2.rap()) { }
0148     bool operator()(const ParticleBase& p) const { return p.rap() >= rapcut.first && p.rap() < rapcut.second; }
0149     pair<double,double> rapcut;
0150   };
0151   using rapInRange = RapInRange;
0152 
0153 
0154   /// Abs rapidity greater-than functor
0155   struct AbsRapGtr : public BoolParticleBaseFunctor {
0156     AbsRapGtr(double absrap) : absrapcut(absrap) { }
0157     AbsRapGtr(const FourMomentum& p) : absrapcut(p.absrap()) { }
0158     bool operator()(const ParticleBase& p) const { return p.absrap() > absrapcut; }
0159     double absrapcut;
0160   };
0161   using absRapGtr = AbsRapGtr;
0162   using absrapGtr = AbsRapGtr;
0163 
0164   /// Abs rapidity momentum less-than functor
0165   struct AbsRapLess : public BoolParticleBaseFunctor {
0166     AbsRapLess(double absrap) : absrapcut(absrap) { }
0167     AbsRapLess(const FourMomentum& p) : absrapcut(p.absrap()) { }
0168     bool operator()(const ParticleBase& p) const { return p.absrap() < absrapcut; }
0169     double absrapcut;
0170   };
0171   using absRapLess = AbsRapLess;
0172   using absrapLess = AbsRapLess;
0173 
0174   /// Abs rapidity in-range functor
0175   struct AbsRapInRange : public BoolParticleBaseFunctor {
0176     AbsRapInRange(const pair<double, double>& absrapcuts) : absrapcut(absrapcuts) { }
0177     AbsRapInRange(double absraplow, double absraphigh) : AbsRapInRange(make_pair(absraplow, absraphigh)) { }
0178     AbsRapInRange(const FourMomentum& p1, const FourMomentum& p2) : AbsRapInRange(p1.absrap(), p2.absrap()) { }
0179     bool operator()(const ParticleBase& p) const { return p.absrap() >= absrapcut.first && p.absrap() < absrapcut.second; }
0180     pair<double,double> absrapcut;
0181   };
0182   using absRapInRange = AbsRapInRange;
0183   using absrapInRange = AbsRapInRange;
0184 
0185 
0186   /// @todo Define dR and dphi functors w.r.t. *multiple* ref vectors, with "all" or "any" semantics
0187 
0188 
0189   /// @f$ \Delta R @f$ (with respect to another 4-momentum, @a vec) greater-than functor
0190   struct DeltaRGtr : public BoolParticleBaseFunctor {
0191     DeltaRGtr(const ParticleBase& vec, double dr, RapScheme scheme=PSEUDORAPIDITY)
0192       : refvec(vec.mom()), drcut(dr), rapscheme(scheme) { }
0193     DeltaRGtr(const FourMomentum& vec, double dr, RapScheme scheme=PSEUDORAPIDITY)
0194       : refvec(vec), drcut(dr), rapscheme(scheme) { }
0195     DeltaRGtr(const Vector3& vec, double dr)
0196       : drcut(dr), rapscheme(PSEUDORAPIDITY) { refvec.setPx(vec.x()); refvec.setPy(vec.y()); refvec.setPz(vec.z()); }
0197     bool operator()(const ParticleBase& p) const { return deltaR(p, refvec, rapscheme) > drcut; }
0198     FourMomentum refvec;
0199     double drcut;
0200     RapScheme rapscheme;
0201   };
0202   using deltaRGtr = DeltaRGtr;
0203 
0204   /// @f$ \Delta R @f$ (with respect to another 4-momentum, @a vec) less-than functor
0205   struct DeltaRLess : public BoolParticleBaseFunctor {
0206     DeltaRLess(const ParticleBase& vec, double dr, RapScheme scheme=PSEUDORAPIDITY)
0207       : refvec(vec.mom()), drcut(dr), rapscheme(scheme) { }
0208     DeltaRLess(const FourMomentum& vec, double dr, RapScheme scheme=PSEUDORAPIDITY)
0209       : refvec(vec), drcut(dr), rapscheme(scheme) { }
0210     DeltaRLess(const Vector3& vec, double dr)
0211       : drcut(dr), rapscheme(PSEUDORAPIDITY) { refvec.setPx(vec.x()); refvec.setPy(vec.y()); refvec.setPz(vec.z()); }
0212     bool operator()(const ParticleBase& p) const { return deltaR(p, refvec, rapscheme) < drcut; }
0213     FourMomentum refvec;
0214     double drcut;
0215     RapScheme rapscheme;
0216   };
0217   using deltaRLess = DeltaRLess;
0218 
0219   /// @f$ \Delta R @f$ (with respect to another 4-momentum, @a vec) in-range functor
0220   struct DeltaRInRange : public BoolParticleBaseFunctor {
0221     DeltaRInRange(const ParticleBase& vec, const pair<double,double>& dr, RapScheme scheme=PSEUDORAPIDITY)
0222       : refvec(vec.mom()), drcut(dr), rapscheme(scheme) { }
0223     DeltaRInRange(const ParticleBase& vec, double drmin, double drmax, RapScheme scheme=PSEUDORAPIDITY)
0224       : DeltaRInRange(vec, make_pair(drmin, drmax), scheme) { }
0225     DeltaRInRange(const FourMomentum& vec, const pair<double,double>& dr, RapScheme scheme=PSEUDORAPIDITY)
0226       : refvec(vec), drcut(dr), rapscheme(scheme) { }
0227     DeltaRInRange(const FourMomentum& vec, double drmin, double drmax, RapScheme scheme=PSEUDORAPIDITY)
0228       : DeltaRInRange(vec, make_pair(drmin, drmax), scheme) { }
0229     DeltaRInRange(const Vector3& vec, const pair<double,double>& dr)
0230       : drcut(dr), rapscheme(PSEUDORAPIDITY) { refvec.setPx(vec.x()); refvec.setPy(vec.y()); refvec.setPz(vec.z()); }
0231     DeltaRInRange(const Vector3& vec, double drmin, double drmax)
0232       : DeltaRInRange(vec, make_pair(drmin, drmax)) { }
0233     bool operator()(const ParticleBase& p) const {
0234       const double dR = deltaR(p, refvec, rapscheme);
0235       return dR >= drcut.first && dR < drcut.second;
0236     }
0237     FourMomentum refvec;
0238     pair<double,double> drcut;
0239     RapScheme rapscheme;
0240   };
0241   using deltaRInRange = DeltaRInRange;
0242 
0243 
0244   /// @f$ |\Delta \phi| @f$ (with respect to another momentum, @a vec) greater-than functor
0245   struct DeltaPhiGtr : public BoolParticleBaseFunctor {
0246     DeltaPhiGtr(const ParticleBase& vec, double dphi)
0247       : refvec(vec.p3()), dphicut(dphi) { }
0248     DeltaPhiGtr(const FourMomentum& vec, double dphi)
0249       : refvec(vec.p3()), dphicut(dphi) { }
0250     DeltaPhiGtr(const Vector3& vec, double dphi)
0251       : refvec(vec), dphicut(dphi) { }
0252     bool operator()(const ParticleBase& p) const { return deltaPhi(p, refvec) > dphicut; }
0253     Vector3 refvec;
0254     double dphicut;
0255   };
0256   using deltaPhiGtr = DeltaPhiGtr;
0257 
0258   /// @f$ |\Delta \phi| @f$ (with respect to another momentum, @a vec) less-than functor
0259   struct DeltaPhiLess : public BoolParticleBaseFunctor {
0260     DeltaPhiLess(const ParticleBase& vec, double dphi)
0261       : refvec(vec.p3()), dphicut(dphi) { }
0262     DeltaPhiLess(const FourMomentum& vec, double dphi)
0263       : refvec(vec.p3()), dphicut(dphi) { }
0264     DeltaPhiLess(const Vector3& vec, double dphi)
0265       : refvec(vec), dphicut(dphi) { }
0266     bool operator()(const ParticleBase& p) const { return deltaPhi(p, refvec) < dphicut; }
0267     Vector3 refvec;
0268     double dphicut;
0269   };
0270   using deltaPhiLess = DeltaPhiLess;
0271 
0272   /// @f$ \Delta \phi @f$ (with respect to another 4-momentum, @a vec) in-range functor
0273   struct DeltaPhiInRange : public BoolParticleBaseFunctor {
0274     DeltaPhiInRange(const ParticleBase& vec, const pair<double,double>& dphi)
0275       : refvec(vec.mom()), dphicut(dphi) { }
0276     DeltaPhiInRange(const ParticleBase& vec, double dphimin, double dphimax)
0277       : DeltaPhiInRange(vec, make_pair(dphimin, dphimax)) { }
0278     DeltaPhiInRange(const FourMomentum& vec, const pair<double,double>& dphi)
0279       : refvec(vec), dphicut(dphi) { }
0280     DeltaPhiInRange(const FourMomentum& vec, double dphimin, double dphimax)
0281       : DeltaPhiInRange(vec, make_pair(dphimin, dphimax)) { }
0282     DeltaPhiInRange(const Vector3& vec, const pair<double,double>& dphi)
0283       : refvec(vec), dphicut(dphi) { }
0284     DeltaPhiInRange(const Vector3& vec, double dphimin, double dphimax)
0285       : DeltaPhiInRange(vec, make_pair(dphimin, dphimax)) { }
0286     bool operator()(const ParticleBase& p) const {
0287       const double dphi = deltaPhi(p, refvec);
0288       return dphi >= dphicut.first && dphi < dphicut.second;
0289     }
0290     Vector3 refvec;
0291     pair<double,double> dphicut;
0292   };
0293   using deltaPhiInRange = DeltaPhiInRange;
0294 
0295 
0296   /// @f$ |\Delta \eta| @f$ (with respect to another momentum, @a vec) greater-than functor
0297   struct DeltaEtaGtr : public BoolParticleBaseFunctor {
0298     DeltaEtaGtr(const ParticleBase& vec, double deta)
0299       : refvec(vec.p3()), detacut(deta) { }
0300     DeltaEtaGtr(const FourMomentum& vec, double deta)
0301       : refvec(vec.p3()), detacut(deta) { }
0302     DeltaEtaGtr(const Vector3& vec, double deta)
0303       : refvec(vec), detacut(deta) { }
0304     bool operator()(const ParticleBase& p) const { return std::abs(deltaEta(p, refvec)) > detacut; }
0305     Vector3 refvec;
0306     double detacut;
0307   };
0308   using deltaEtaGtr = DeltaEtaGtr;
0309 
0310   /// @f$ |\Delta \eta| @f$ (with respect to another momentum, @a vec) less-than functor
0311   struct DeltaEtaLess : public BoolParticleBaseFunctor {
0312     DeltaEtaLess(const ParticleBase& vec, double deta)
0313       : refvec(vec.p3()), detacut(deta) { }
0314     DeltaEtaLess(const FourMomentum& vec, double deta)
0315       : refvec(vec.p3()), detacut(deta) { }
0316     DeltaEtaLess(const Vector3& vec, double deta)
0317       : refvec(vec), detacut(deta) { }
0318     bool operator()(const ParticleBase& p) const { return std::abs(deltaEta(p, refvec)) < detacut; }
0319     Vector3 refvec;
0320     double detacut;
0321   };
0322   using deltaEtaLess = DeltaEtaLess;
0323 
0324   /// @f$ \Delta \eta @f$ (with respect to another 4-momentum, @a vec) in-range functor
0325   struct DeltaEtaInRange : public BoolParticleBaseFunctor {
0326     DeltaEtaInRange(const ParticleBase& vec, const pair<double,double>& deta)
0327       : refvec(vec.mom()), detacut(deta) { }
0328     DeltaEtaInRange(const ParticleBase& vec, double detamin, double detamax)
0329       : DeltaEtaInRange(vec, make_pair(detamin, detamax)) { }
0330     DeltaEtaInRange(const FourMomentum& vec, const pair<double,double>& deta)
0331       : refvec(vec), detacut(deta) { }
0332     DeltaEtaInRange(const FourMomentum& vec, double detamin, double detamax)
0333       : DeltaEtaInRange(vec, make_pair(detamin, detamax)) { }
0334     DeltaEtaInRange(const Vector3& vec, const pair<double,double>& deta)
0335       : refvec(vec), detacut(deta) { }
0336     DeltaEtaInRange(const Vector3& vec, double detamin, double detamax)
0337       : DeltaEtaInRange(vec, make_pair(detamin, detamax)) { }
0338     bool operator()(const ParticleBase& p) const {
0339       const double deta = deltaEta(p, refvec);
0340       return deta >= detacut.first && deta < detacut.second;
0341     }
0342     Vector3 refvec;
0343     pair<double,double> detacut;
0344   };
0345   using deltaEtaInRange = DeltaEtaInRange;
0346 
0347 
0348   /// @f$ |\Delta y| @f$ (with respect to another momentum, @a vec) greater-than functor
0349   struct DeltaRapGtr : public BoolParticleBaseFunctor {
0350     DeltaRapGtr(const ParticleBase& vec, double drap)
0351       : refvec(vec.mom()), drapcut(drap) { }
0352     DeltaRapGtr(const FourMomentum& vec, double drap)
0353       : refvec(vec), drapcut(drap) { }
0354     bool operator()(const ParticleBase& p) const { return std::abs(deltaRap(p, refvec)) > drapcut; }
0355     FourMomentum refvec;
0356     double drapcut;
0357   };
0358   using deltaRapGtr = DeltaRapGtr;
0359 
0360   /// @f$ |\Delta y| @f$ (with respect to another momentum, @a vec) less-than functor
0361   struct DeltaRapLess : public BoolParticleBaseFunctor {
0362     DeltaRapLess(const ParticleBase& vec, double drap)
0363       : refvec(vec.mom()), drapcut(drap) { }
0364     DeltaRapLess(const FourMomentum& vec, double drap)
0365       : refvec(vec), drapcut(drap) { }
0366     bool operator()(const ParticleBase& p) const { return std::abs(deltaRap(p, refvec)) < drapcut; }
0367     FourMomentum refvec;
0368     double drapcut;
0369   };
0370   using deltaRapLess = DeltaRapLess;
0371 
0372   /// @f$ \Delta y @f$ (with respect to another 4-momentum, @a vec) in-range functor
0373   struct DeltaRapInRange : public BoolParticleBaseFunctor {
0374     DeltaRapInRange(const ParticleBase& vec, const pair<double,double>& drap)
0375       : refvec(vec.mom()), drapcut(drap) { }
0376     DeltaRapInRange(const ParticleBase& vec, double drapmin, double drapmax)
0377       : DeltaRapInRange(vec, make_pair(drapmin, drapmax)) { }
0378     DeltaRapInRange(const FourMomentum& vec, const pair<double,double>& drap)
0379       : refvec(vec), drapcut(drap) { }
0380     DeltaRapInRange(const FourMomentum& vec, double drapmin, double drapmax)
0381       : DeltaRapInRange(vec, make_pair(drapmin, drapmax)) { }
0382     bool operator()(const ParticleBase& p) const {
0383       const double drap = deltaRap(p, refvec);
0384       return drap >= drapcut.first && drap < drapcut.second;
0385     }
0386     FourMomentum refvec;
0387     pair<double,double> drapcut;
0388   };
0389   using deltaRapInRange = DeltaRapInRange;
0390 
0391   /// @}
0392 
0393 
0394   /// @defgroup particlebaseutils_pb2dbl ParticleBase comparison -> double functors
0395   /// @todo Move to FourMomentum functions
0396   ///
0397   /// To be passed to transform()any(jets, DeltaRLess(electron, 0.4))
0398   /// @{
0399 
0400   /// Base type for Particle -> double functors
0401   struct DoubleParticleBaseFunctor {
0402     virtual double operator()(const ParticleBase& p) const = 0;
0403     virtual ~DoubleParticleBaseFunctor() {}
0404   };
0405 
0406   /// Calculator of @f$ \Delta R @f$ with respect to a given momentum
0407   struct DeltaRWRT : public DoubleParticleBaseFunctor {
0408     DeltaRWRT(const ParticleBase& pb, RapScheme scheme=PSEUDORAPIDITY) : p(pb.mom()), rapscheme(scheme) {}
0409     DeltaRWRT(const FourMomentum& p4, RapScheme scheme=PSEUDORAPIDITY) : p(p4), rapscheme(scheme) {}
0410     DeltaRWRT(const Vector3& p3) : p(p3.mod(), p3.x(), p3.y(), p3.z()), rapscheme(PSEUDORAPIDITY) {}
0411     double operator()(const ParticleBase& pb) const { return deltaR(p, pb, rapscheme); }
0412     double operator()(const FourMomentum& p4) const { return deltaR(p, p4, rapscheme); }
0413     double operator()(const Vector3& p3) const { return deltaR(p, p3); }
0414     const FourMomentum p;
0415     RapScheme rapscheme;
0416   };
0417   using deltaRWRT = DeltaRWRT;
0418 
0419   /// Calculator of @f$ \Delta \phi @f$ with respect to a given momentum
0420   struct DeltaPhiWRT : public DoubleParticleBaseFunctor {
0421     DeltaPhiWRT(const ParticleBase& pb) : p(pb.mom().vector3()) {}
0422     DeltaPhiWRT(const FourMomentum& p4) : p(p4.vector3()) {}
0423     DeltaPhiWRT(const Vector3& p3) : p(p3) {}
0424     double operator()(const ParticleBase& pb) const { return deltaPhi(p, pb); }
0425     double operator()(const FourMomentum& p4) const { return deltaPhi(p, p4); }
0426     double operator()(const Vector3& p3) const { return deltaPhi(p, p3); }
0427     const Vector3 p;
0428   };
0429   using deltaPhiWRT = DeltaPhiWRT;
0430 
0431   /// Calculator of @f$ \Delta \eta @f$ with respect to a given momentum
0432   struct DeltaEtaWRT : public DoubleParticleBaseFunctor {
0433     DeltaEtaWRT(const ParticleBase& pb) : p(pb.mom().vector3()) {}
0434     DeltaEtaWRT(const FourMomentum& p4) : p(p4.vector3()) {}
0435     DeltaEtaWRT(const Vector3& p3) : p(p3) {}
0436     double operator()(const ParticleBase& pb) const { return deltaEta(p, pb); }
0437     double operator()(const FourMomentum& p4) const { return deltaEta(p, p4); }
0438     double operator()(const Vector3& p3) const { return deltaEta(p, p3); }
0439     const Vector3 p;
0440   };
0441   using deltaEtaWRT = DeltaEtaWRT;
0442 
0443   /// Calculator of @f$ |\Delta \eta| @f$ with respect to a given momentum
0444   struct AbsDeltaEtaWRT : public DoubleParticleBaseFunctor {
0445     AbsDeltaEtaWRT(const ParticleBase& pb) : p(pb.mom().vector3()) {}
0446     AbsDeltaEtaWRT(const FourMomentum& p4) : p(p4.vector3()) {}
0447     AbsDeltaEtaWRT(const Vector3& p3) : p(p3) {}
0448     double operator()(const ParticleBase& pb) const { return fabs(deltaEta(p, pb)); }
0449     double operator()(const FourMomentum& p4) const { return fabs(deltaEta(p, p4)); }
0450     double operator()(const Vector3& p3) const { return fabs(deltaEta(p, p3)); }
0451     const Vector3 p;
0452   };
0453   using absDeltaEtaWRT = AbsDeltaEtaWRT;
0454 
0455   /// Calculator of @f$ \Delta y @f$ with respect to a given momentum
0456   struct DeltaRapWRT : public DoubleParticleBaseFunctor {
0457     DeltaRapWRT(const ParticleBase& pb) : p(pb.mom()) {}
0458     DeltaRapWRT(const FourMomentum& p4) : p(p4) {}
0459     double operator()(const ParticleBase& pb) const { return deltaRap(p, pb); }
0460     double operator()(const FourMomentum& p4) const { return deltaRap(p, p4); }
0461     const FourMomentum p;
0462   };
0463   using deltaRapWRT = DeltaRapWRT;
0464 
0465   /// Calculator of @f$ |\Delta y| @f$ with respect to a given momentum
0466   struct AbsDeltaRapWRT : public DoubleParticleBaseFunctor {
0467     AbsDeltaRapWRT(const ParticleBase& pb) : p(pb.mom()) {}
0468     AbsDeltaRapWRT(const FourMomentum& p4) : p(p4) {}
0469     double operator()(const ParticleBase& pb) const { return fabs(deltaRap(p, pb)); }
0470     double operator()(const FourMomentum& p4) const { return fabs(deltaRap(p, p4)); }
0471     const FourMomentum p;
0472   };
0473   using absDeltaRapWRT = AbsDeltaRapWRT;
0474 
0475   /// @}
0476 
0477 
0478   /// @defgroup particlebaseutils_uberfilt Next-level filtering
0479   /// @{
0480 
0481   template<typename PBCONTAINER1, typename PBCONTAINER2>
0482   inline void idiscardIfAny(PBCONTAINER1& tofilter, const PBCONTAINER2& tocompare,
0483                             typename std::function<bool(const typename PBCONTAINER1::value_type&,
0484                                                         const typename PBCONTAINER2::value_type&)> fn) {
0485     for (const auto& pbcmp : tocompare) {
0486       idiscard(tofilter, [&](const typename PBCONTAINER1::value_type& pbfilt){ return fn(pbfilt, pbcmp); });
0487     }
0488   }
0489 
0490   template<typename PBCONTAINER1, typename PBCONTAINER2>
0491   inline PBCONTAINER1 discardIfAny(const PBCONTAINER1& tofilter, const PBCONTAINER2& tocompare,
0492                                    typename std::function<bool(const typename PBCONTAINER1::value_type&,
0493                                                                const typename PBCONTAINER2::value_type&)> fn) {
0494     PBCONTAINER1 tmp{tofilter};
0495     idiscardIfAny(tmp, tocompare, fn);
0496     return tmp;
0497   }
0498 
0499 
0500   template<typename PBCONTAINER1, typename PBCONTAINER2>
0501   inline PBCONTAINER1 selectIfAny(const PBCONTAINER1& tofilter, const PBCONTAINER2& tocompare,
0502                                   typename std::function<bool(const typename PBCONTAINER1::value_type&,
0503                                                               const typename PBCONTAINER2::value_type&)> fn) {
0504     PBCONTAINER1 selected;
0505     for (const auto& pbfilt : tofilter) {
0506       if (any(tocompare, [&](const typename PBCONTAINER2::value_type& pbcmp){ return fn(pbfilt, pbcmp); })) {
0507         selected += pbfilt;
0508       }
0509     }
0510     return selected;
0511   }
0512 
0513   template<typename PBCONTAINER1, typename PBCONTAINER2>
0514   inline void iselectIfAny(PBCONTAINER1& tofilter, const PBCONTAINER2& tocompare,
0515                            typename std::function<bool(const typename PBCONTAINER1::value_type&,
0516                                                        const typename PBCONTAINER2::value_type&)> fn) {
0517     tofilter = selectIfAny(tofilter, tocompare, fn);
0518   }
0519 
0520 
0521 
0522   template<typename PBCONTAINER1, typename PBCONTAINER2>
0523   inline PBCONTAINER1 discardIfAll(const PBCONTAINER1& tofilter, const PBCONTAINER2& tocompare,
0524                                    typename std::function<bool(const typename PBCONTAINER1::value_type&,
0525                                                                const typename PBCONTAINER2::value_type&)> fn) {
0526     PBCONTAINER1 selected;
0527     for (const auto& pbfilt : tofilter) {
0528       if (!all(tocompare, [&](const typename PBCONTAINER2::value_type& pbcmp){ return fn(pbfilt, pbcmp); })) {
0529         selected += pbfilt;
0530       }
0531     }
0532     return selected;
0533   }
0534 
0535   template<typename PBCONTAINER1, typename PBCONTAINER2>
0536   inline void idiscardIfAll(PBCONTAINER1& tofilter, const PBCONTAINER2& tocompare,
0537                             typename std::function<bool(const typename PBCONTAINER1::value_type&,
0538                                                         const typename PBCONTAINER2::value_type&)> fn) {
0539     tofilter = discardIfAll(tofilter, tocompare, fn);
0540   }
0541 
0542 
0543   template<typename PBCONTAINER1, typename PBCONTAINER2>
0544   inline PBCONTAINER1 selectIfAll(const PBCONTAINER1& tofilter, const PBCONTAINER2& tocompare,
0545                                   typename std::function<bool(const typename PBCONTAINER1::value_type&,
0546                                                               const typename PBCONTAINER2::value_type&)> fn) {
0547     PBCONTAINER1 selected;
0548     for (const auto& pbfilt : tofilter) {
0549       if (all(tocompare, [&](const typename PBCONTAINER2::value_type& pbcmp){ return fn(pbfilt, pbcmp); })) {
0550         selected += pbfilt;
0551       }
0552     }
0553     return selected;
0554   }
0555 
0556   template<typename PBCONTAINER1, typename PBCONTAINER2>
0557   inline void iselectIfAll(PBCONTAINER1& tofilter, const PBCONTAINER2& tocompare,
0558                            typename std::function<bool(const typename PBCONTAINER1::value_type&,
0559                                                        const typename PBCONTAINER2::value_type&)> fn) {
0560     tofilter = selectIfAll(tofilter, tocompare, fn);
0561   }
0562 
0563   /// @}
0564 
0565 
0566   /// @defgroup particlebaseutils_iso Isolation helpers
0567   /// @{
0568 
0569   template<typename PBCONTAINER1, typename PBCONTAINER2>
0570   inline void idiscardIfAnyDeltaRLess(PBCONTAINER1& tofilter, const PBCONTAINER2& tocompare, double dR) {
0571     for (const typename PBCONTAINER2::value_type& pb : tocompare) {
0572       idiscard(tofilter, deltaRLess(pb, dR));
0573     }
0574   }
0575 
0576   template<typename PBCONTAINER1, typename PBCONTAINER2>
0577   inline PBCONTAINER1 discardIfAnyDeltaRLess(const PBCONTAINER1& tofilter, const PBCONTAINER2& tocompare, double dR) {
0578     PBCONTAINER1 tmp{tofilter};
0579     idiscardIfAnyDeltaRLess(tmp, tocompare, dR);
0580     return tmp;
0581   }
0582 
0583   template<typename PBCONTAINER1, typename PBCONTAINER2>
0584   inline void idiscardIfAnyDeltaPhiLess(PBCONTAINER1& tofilter, const PBCONTAINER2& tocompare, double dphi) {
0585     for (const typename PBCONTAINER2::value_type& pb : tocompare) {
0586       idiscard(tofilter, deltaPhiLess(pb, dphi));
0587     }
0588   }
0589 
0590   template<typename PBCONTAINER1, typename PBCONTAINER2>
0591   inline PBCONTAINER1 discardIfAnyDeltaPhiLess(const PBCONTAINER1& tofilter, const PBCONTAINER2& tocompare, double dphi) {
0592     PBCONTAINER1 tmp{tofilter};
0593     idiscardIfAnyDeltaPhiLess(tmp, tocompare, dphi);
0594     return tmp;
0595   }
0596 
0597 
0598 
0599   template<typename PBCONTAINER1, typename PBCONTAINER2>
0600   inline PBCONTAINER1 selectIfAnyDeltaRLess(const PBCONTAINER1& tofilter, const PBCONTAINER2& tocompare, double dR) {
0601     PBCONTAINER1 selected;
0602     for (const typename PBCONTAINER1::value_type& f : tofilter) {
0603       if (any(tocompare, deltaRLess(f, dR))) selected.push_back(f);
0604     }
0605     return selected;
0606   }
0607 
0608   template<typename PBCONTAINER1, typename PBCONTAINER2>
0609   inline void iselectIfAnyDeltaRLess(PBCONTAINER1& tofilter, const PBCONTAINER2& tocompare, double dR) {
0610     tofilter = selectIfAnyDeltaRLess(tofilter, tocompare, dR);
0611   }
0612 
0613 
0614   template<typename PBCONTAINER1, typename PBCONTAINER2>
0615   inline PBCONTAINER1 selectIfAnyDeltaPhiLess(const PBCONTAINER1& tofilter, const PBCONTAINER2& tocompare, double dphi) {
0616     PBCONTAINER1 selected;
0617     for (const typename PBCONTAINER1::value_type& f : tofilter) {
0618       if (any(tocompare, deltaPhiLess(f, dphi))) selected.push_back(f);
0619     }
0620     return selected;
0621   }
0622 
0623   template<typename PBCONTAINER1, typename PBCONTAINER2>
0624   inline void iselectIfAnyDeltaPhiLess(PBCONTAINER1& tofilter, const PBCONTAINER2& tocompare, double dphi) {
0625     tofilter = selectIfAnyDeltaPhiLess(tofilter, tocompare, dphi);
0626   }
0627 
0628 
0629   /// @todo Add 'all' variants
0630 
0631   /// @}
0632 
0633 
0634   /// @defgroup particlebaseutils_kin Unbound functions for kinematic properties
0635   ///
0636   /// @todo Mostly move to functions on FourMomentum
0637   /// @note In a sub-namespace (imported by default) for protection
0638   /// @{
0639   namespace Kin {
0640 
0641     /// Unbound function access to momentum
0642     inline FourMomentum mom(const ParticleBase& p) { return p.mom(); }
0643     /// Unbound function access to momentum
0644     inline FourMomentum p4(const ParticleBase& p) { return p.mom(); }
0645 
0646     /// Unbound function access to p3
0647     inline Vector3 p3(const ParticleBase& p) { return p.p3(); }
0648 
0649     /// Unbound function access to pTvec
0650     inline Vector3 pTvec(const ParticleBase& p) { return p.pTvec(); }
0651 
0652     /// Unbound function access to p
0653     inline double p(const ParticleBase& p) { return p.p(); }
0654 
0655     /// Unbound function access to pT
0656     inline double pT(const ParticleBase& p) { return p.pT(); }
0657 
0658     /// Unbound function access to E
0659     inline double E(const ParticleBase& p) { return p.E(); }
0660 
0661     /// Unbound function access to ET
0662     inline double Et(const ParticleBase& p) { return p.Et(); }
0663 
0664     /// Unbound function access to eta
0665     inline double eta(const ParticleBase& p) { return p.eta(); }
0666 
0667     /// Unbound function access to abseta
0668     inline double abseta(const ParticleBase& p) { return p.abseta(); }
0669 
0670     /// Unbound function access to rapidity
0671     inline double rap(const ParticleBase& p) { return p.rap(); }
0672 
0673     /// Unbound function access to abs rapidity
0674     inline double absrap(const ParticleBase& p) { return p.absrap(); }
0675 
0676     /// Unbound function access to mass
0677     inline double mass(const ParticleBase& p) { return p.mass(); }
0678 
0679 
0680     // /// Unbound function access to pair pT
0681     // inline double pairPt(const ParticleBase& p1, const ParticleBase& p2) { return (p1.mom() + p2.mom()).pT(); }
0682 
0683     // /// Unbound function access to pair mass
0684     // inline double pairMass(const ParticleBase& p1, const ParticleBase& p2) { return (p1.mom() + p2.mom()).mass(); }
0685 
0686 
0687     /// @brief Get the mass of a ParticleBase and a P4
0688     inline double mass(const ParticleBase& p, const FourMomentum& p4) {
0689       return mass(p.mom(), p4);
0690     }
0691 
0692     /// @brief Get the mass of a ParticleBase and a P4
0693     inline double mass(const FourMomentum& p4, const ParticleBase& p) {
0694       return mass(p4, p.mom());
0695     }
0696 
0697     /// @brief Get the mass of a pair of ParticleBase (as separate args)
0698     inline double mass(const ParticleBase& p1, const ParticleBase& p2) {
0699       return mass(p1.mom(), p2.mom());
0700     }
0701 
0702     /// @brief Get the mass^2 of a ParticleBase and a P4
0703     inline double mass2(const ParticleBase& p, const P4& p4) {
0704       return mass2(p.mom(), p4);
0705     }
0706 
0707     /// @brief Get the mass^2 of a ParticleBase and a P4
0708     inline double mass2(const P4& p4, const ParticleBase& p) {
0709       return mass2(p4, p.mom());
0710     }
0711 
0712     /// @brief Get the mass^2 of a pair of ParticleBase (as separate args)
0713     inline double mass2(const ParticleBase& p1, const ParticleBase& p2) {
0714       return mass2(p1.mom(), p2.mom());
0715     }
0716 
0717     /// @brief Get the transverse mass of a ParticleBase and a P4
0718     ///
0719     /// @note This ignores the particle mass and just computes mT from the 3-vectors
0720     /// @todo Fix!!
0721     inline double mT(const ParticleBase& p, const P4& p4) {
0722       return mT(p.mom(), p4);
0723     }
0724 
0725     /// @brief Get the transverse mass of a ParticleBase and a P4
0726     ///
0727     /// @note This ignores the particle mass and just computes mT from the 3-vectors
0728     /// @todo Fix!!
0729     inline double mT(const P4& p4, const ParticleBase& p) {
0730       return mT(p4, p.mom());
0731     }
0732 
0733     /// @brief Get the transverse mass of a pair of ParticleBase (as separate args)
0734     ///
0735     /// @note This ignores the particle mass and just computes mT from the 3-vectors
0736     /// @todo Fix!!
0737     inline double mT(const ParticleBase& p1, const ParticleBase& p2) {
0738       return mT(p1.mom(), p2.mom());
0739     }
0740 
0741     /// @brief Get the transverse momentum of a ParticleBase and a P4
0742     inline double pT(const ParticleBase& p, const P4& p4) {
0743       return pT(p.mom(), p4);
0744     }
0745 
0746     /// @brief Get the transverse momentum of a ParticleBase and a P4
0747     inline double pT(const P4& p4, const ParticleBase& p) {
0748       return pT(p4, p.mom());
0749     }
0750 
0751     /// @brief Get the transverse momentum of a pair of ParticleBase (as separate args)
0752     inline double pT(const ParticleBase& p1, const ParticleBase& p2) {
0753       return pT(p1.mom(), p2.mom());
0754     }
0755 
0756   }
0757 
0758   // Import Kin namespace into Rivet
0759   using namespace Kin;
0760 
0761   /// @}
0762 
0763 
0764   /// @defgroup pbcontcombutils ParticleBase container-combinatorics utils
0765   /// @{
0766 
0767   /// @brief Return the index from a vector which best matches mass(c[i]) to the target value
0768   ///
0769   /// A specialisation of closestMatchIndex from Utils.hh, with the
0770   /// function bound to Kin::mass as a common use-case.
0771   // auto closestMassIndex = std::bind(closestMatchIndex, std::placeholders::_1, Kin::mass, std::placeholders::_2, std::placeholders::_3);
0772   template <typename CONTAINER, typename = isCIterable<CONTAINER>>
0773   inline int closestMassIndex(CONTAINER&& c, double mtarget, double mmin=-DBL_MAX, double mmax=DBL_MAX) {
0774     return closestMatchIndex(std::forward<CONTAINER>(c), Kin::mass, mtarget, mmin, mmax);
0775   }
0776 
0777   /// @brief Return the indices from two vectors which best match mass(c1[i], c2[j]) to the target value
0778   ///
0779   /// A specialisation of closestMatchIndex from Utils.hh, with the
0780   /// function bound to Kin::mass as a common use-case.
0781   template <typename CONTAINER1, typename CONTAINER2, typename = isCIterable<CONTAINER1, CONTAINER2>>
0782   inline pair<int,int> closestMassIndices(CONTAINER1&& c1, CONTAINER2&& c2,
0783                            double mtarget, double mmin=-DBL_MAX, double mmax=DBL_MAX) {
0784     return closestMatchIndices(std::forward<CONTAINER1>(c1),
0785                                std::forward<CONTAINER2>(c2), Kin::mass, mtarget, mmin, mmax);
0786   }
0787 
0788   /// @brief Return the index from a vector which best matches mass(c[i], x) to the target value
0789   ///
0790   /// A specialisation of closestMatchIndex from Utils.hh, with the
0791   /// function bound to Kin::mass as a common use-case.
0792   template <typename CONTAINER, typename T, typename = isCIterable<CONTAINER>>
0793   inline int closestMassIndex(CONTAINER&& c, const T& x,
0794                               double mtarget, double mmin=-DBL_MAX, double mmax=DBL_MAX) {
0795     using FN = double(const ParticleBase&,const T&);
0796     return closestMatchIndex<CONTAINER, T, FN>(std::forward<CONTAINER>(c), x, Kin::mass, mtarget, mmin, mmax);
0797   }
0798 
0799 
0800   /// @brief Return the index from a vector which best matches mass(x, c[j]) to the target value
0801   ///
0802   /// A specialisation of closestMatchIndex from Utils.hh, with the
0803   /// function bound to Kin::mass as a common use-case.
0804   template <typename CONTAINER, typename T, typename = isCIterable<CONTAINER>>
0805   inline int closestMassIndex(T&& x, CONTAINER&& c,
0806                               double mtarget, double mmin=-DBL_MAX, double mmax=DBL_MAX) {
0807     return closestMatchIndex(std::forward<T>(x), std::forward<CONTAINER>(c), Kin::mass, mtarget, mmin, mmax);
0808   }
0809 
0810   /// @}
0811 
0812 
0813   /// @defgroup pbcontphysutils ParticleBase-container physics utils
0814   /// @{
0815 
0816   /// @brief Return the pT sum of a collection of particle-like objects
0817   template <typename CONTAINER, typename = isCIterable<CONTAINER>>
0818   inline double sumPt(CONTAINER&& c) {
0819     return sum(c, Kin::pT, 0.0);
0820   }
0821 
0822   /// @brief Return the four-momentum sum of a collection of particle-like objects
0823   template <typename CONTAINER, typename = isCIterable<CONTAINER>>
0824   inline FourMomentum sumP4(CONTAINER&& c) {
0825     return sum(c, Kin::p4, FourMomentum());
0826   }
0827 
0828   /// @brief Return the four-momentum sum of a collection of particle-like objects
0829   template <typename CONTAINER, typename = isCIterable<CONTAINER>>
0830   inline Vector3 sumP3(CONTAINER&& c) {
0831     return sum(c, Kin::p3, Vector3());
0832   }
0833 
0834   /// @brief Return the four-momentum sum of a collection of particle-like objects
0835   template <typename CONTAINER, typename = isCIterable<CONTAINER>>
0836   inline Vector3 sumPtVec(CONTAINER&& c) {
0837     return sum(c, Kin::pTvec, Vector3());
0838   }
0839 
0840 
0841   /// @brief Calculate the minimum phi separation in a collection of particle-like objects
0842   ///
0843   /// @note Returns in the range 0..pi
0844   ///
0845   /// @todo Generalise to min of any 2-ParticleBase function?
0846   template <typename CONTAINER, typename = isCIterable<CONTAINER>>
0847   inline double deltaPhiMin(CONTAINER&& c) {
0848     double rtn = DBL_MAX;
0849     for (size_t i = 0; i < c.size(); ++i) {
0850       for (size_t j = i; j < c.size(); ++j) {
0851         if (i == j) continue;
0852         const double dphi = fabs(deltaPhi(c[i], c[j]));
0853         if (dphi < rtn) rtn = dphi;
0854       }
0855     }
0856     return rtn;
0857   }
0858 
0859   /// @brief Calculate the minimum phi separation between two collections of particle-like objects
0860   ///
0861   /// @note Returns in the range 0..pi
0862   ///
0863   /// @warning Assumes the two vectors are mutually exclusive, so does not protect against self-comparisons.
0864   ///
0865   /// @todo Generalise to min of any 2-ParticleBase function?
0866   template <typename CONTAINER1, typename CONTAINER2, typename = isCIterable<CONTAINER1, CONTAINER2>>
0867   inline double deltaPhiMin(CONTAINER1&& c1, CONTAINER2&& c2) {
0868     double rtn = DBL_MAX;
0869     for (size_t i = 0; i < c1.size(); ++i) {
0870       for (size_t j = 0; j < c2.size(); ++j) {
0871         const double dphi = deltaPhi(c1[i], c2[j]);
0872         if (dphi < rtn) rtn = dphi;
0873       }
0874     }
0875     return rtn;
0876   }
0877 
0878 
0879 
0880   /// @brief Calculate the minimum pseudorapidity separation in a collection of particle-like objects
0881   ///
0882   /// @note Returns positive values
0883   ///
0884   /// @todo Generalise to min of any 2-ParticleBase function?
0885   template <typename CONTAINER, typename = isCIterable<CONTAINER>>
0886   inline double deltaEtaMin(CONTAINER&& c) {
0887     double rtn = DBL_MAX;
0888     for (size_t i = 0; i < c.size(); ++i) {
0889       for (size_t j = i; j < c.size(); ++j) {
0890         if (i == j) continue;
0891         const double deta = fabs(deltaEta(c[i], c[j]));
0892         if (deta < rtn) rtn = deta;
0893       }
0894     }
0895     return rtn;
0896   }
0897 
0898   /// @brief Calculate the minimum pseudorapidity separation between two collections of particle-like objects
0899   ///
0900   /// @note Returns positive value
0901   ///
0902   /// @warning Assumes the two vectors are mutually exclusive, so does not protect against self-comparisons.
0903   ///
0904   /// @todo Generalise to min of any 2-ParticleBase function?
0905   template <typename CONTAINER1, typename CONTAINER2, typename = isCIterable<CONTAINER1, CONTAINER2>>
0906   inline double deltaEtaMin(CONTAINER1&& c1, CONTAINER2&& c2) {
0907     double rtn = DBL_MAX;
0908     for (size_t i = 0; i < c1.size(); ++i) {
0909       for (size_t j = 0; j < c2.size(); ++j) {
0910         const double deta = deltaEta(c1[i], c2[j]);
0911         if (deta < rtn) rtn = deta;
0912       }
0913     }
0914     return rtn;
0915   }
0916 
0917 
0918   /// @brief Calculate the minimum rapidity separation in a collection of particle-like objects
0919   ///
0920   /// @note Returns positive values
0921   ///
0922   /// @todo Generalise to min of any 2-ParticleBase function?
0923   template <typename CONTAINER, typename = isCIterable<CONTAINER>>
0924   inline double deltaRapMin(CONTAINER&& c) {
0925     double rtn = DBL_MAX;
0926     for (size_t i = 0; i < c.size(); ++i) {
0927       for (size_t j = i; j < c.size(); ++j) {
0928         if (i == j) continue;
0929         const double dy = fabs(deltaRap(c[i], c[j]));
0930         if (dy < rtn) rtn = dy;
0931       }
0932     }
0933     return rtn;
0934   }
0935 
0936   /// @brief Calculate the minimum pseudorapidity separation between two collections of particle-like objects
0937   ///
0938   /// @note Returns positive values
0939   ///
0940   /// @warning Assumes the two vectors are mutually exclusive, so does not protect against self-comparisons.
0941   ///
0942   /// @todo Generalise to min of any 2-ParticleBase function?
0943   template <typename CONTAINER1, typename CONTAINER2, typename = isCIterable<CONTAINER1, CONTAINER2>>
0944   inline double deltaRapMin(CONTAINER1&& c1, CONTAINER2&& c2) {
0945     double rtn = DBL_MAX;
0946     for (size_t i = 0; i < c1.size(); ++i) {
0947       for (size_t j = 0; j < c2.size(); ++j) {
0948         const double dy = deltaRap(c1[i], c2[j]);
0949         if (dy < rtn) rtn = dy;
0950       }
0951     }
0952     return rtn;
0953   }
0954 
0955 
0956 
0957   /// @brief Calculate the minimum rapidity separation in a collection of particle-like objects
0958   ///
0959   /// @note Returns positive values
0960   ///
0961   /// @todo Generalise to min of any 2-ParticleBase function?
0962   template <typename CONTAINER, typename = isCIterable<CONTAINER>>
0963   inline double deltaRMin(CONTAINER&& c) {
0964     double rtn = DBL_MAX;
0965     for (size_t i = 0; i < c.size(); ++i) {
0966       for (size_t j = i; j < c.size(); ++j) {
0967         if (i == j) continue;
0968         const double dr = fabs(deltaR(c[i], c[j]));
0969         if (dr < rtn) rtn = dr;
0970       }
0971     }
0972     return rtn;
0973   }
0974 
0975   /// @brief Calculate the minimum pseudorapidity separation between two collections of particle-like objects
0976   ///
0977   /// @note Returns positive values
0978   ///
0979   /// @warning Assumes the two vectors are mutually exclusive, so does not protect against self-comparisons.
0980   ///
0981   /// @todo Generalise to min of any 2-ParticleBase function?
0982   template <typename CONTAINER1, typename CONTAINER2, typename = isCIterable<CONTAINER1, CONTAINER2>>
0983   inline double deltaRMin(CONTAINER1&& c1, CONTAINER2&& c2) {
0984     double rtn = DBL_MAX;
0985     for (size_t i = 0; i < c1.size(); ++i) {
0986       for (size_t j = 0; j < c2.size(); ++j) {
0987         const double dr = deltaR(c1[i], c2[j]);
0988         if (dr < rtn) rtn = dr;
0989       }
0990     }
0991     return rtn;
0992   }
0993 
0994   /// @}
0995 
0996 
0997   /// @}
0998 
0999 }
1000 
1001 #endif