packages/fjcontrib/centauro-euclidean-distance.patch

0001 Index: Centauro.cc
0002 ===================================================================
0003 --- Centauro.cc (revision 1555)
0004 +++ Centauro.cc (working copy)
0005 @@ -15,9 +15,10 @@
0006  //----------------------------------------------------------------------x
0007
0008  #include "Centauro.hh"
0009 -#include "fastjet/NNH.hh"
0010 +#include "fastjet/NNFJN2Plain.hh"
0011
0012  // strings and streams
0013 +#include <cmath>
0014  #include <sstream>
0015  #include <limits>
0016
0017 @@ -75,14 +76,20 @@
0018        }
0019        // beam distance
0020        diB = 1.0;
0021 +
0022 +      // Precompute 2D Cartesian coordinates in (etabar·cos φ, etabar·sin φ) space.
0023 +      // The Centauro distance is a pure 2D Euclidean squared distance:
0024 +      //   dij = ((cx_i-cx_j)^2 + (cy_i-cy_j)^2) / R^2
0025 +      // This eliminates cos() from the O(N^2) inner loop.
0026 +      _cx = etabar * std::cos(phi);
0027 +      _cy = etabar * std::sin(phi);
0028 +      _R2 = R * R;
0029      }
0030
0031      double distance(const CentauroBriefJet * jet) const {
0032 -
0033 -      double dij = pow(etabar - jet->etabar, 2.0) + 2*etabar*jet->etabar*(1-cos(phi- jet->phi));
0034 -      dij = dij/pow(R,2.0);
0035 -
0036 -      return dij;
0037 +      double dx = _cx - jet->_cx;
0038 +      double dy = _cy - jet->_cy;
0039 +      return (dx*dx + dy*dy) / _R2;
0040      }
0041
0042      double beam_distance() const {
0043 @@ -89,10 +96,21 @@
0044        return diB;
0045      }
0046
0047 +    // NNFJN2Plain interface: dij = min(mom_i, mom_j) * geom(i,j)
0048 +    // Since diB = 1.0 for all particles, momentum_factor = 1 trivially
0049 +    // satisfies the factorisation requirement.
0050 +    double geometrical_distance(const CentauroBriefJet * jet) const {
0051 +      return distance(jet);
0052 +    }
0053 +    double geometrical_beam_distance() const { return diB; }
0054 +    double momentum_factor()           const { return 1.0; }
0055 +
0056      double pT, phi, nx, ny, nz;
0057      double etabar;
0058      double diB;
0059      double R, gammaE, gammaPz, Q;
0060 +  private:
0061 +    double _cx, _cy, _R2;
0062    };
0063
0064
0065 @@ -109,7 +127,7 @@
0066      int njets = cs.jets().size();
0067      CentauroInfo vinfo(R(), gammaE(), gammaPz());
0068
0069 -    NNH<CentauroBriefJet,CentauroInfo> nnh(cs.jets(),&vinfo);
0070 +    NNFJN2Plain<CentauroBriefJet,CentauroInfo> nnh(cs.jets(),&vinfo);
0071
0072      while (njets > 0) {
0073        int i, j, k;