fastjet/contrib/CMPPlugin.hh

0001 /// Implementation of the flavoured anti-kt algorithm
0002 /// by Michal Czakon, Alexander Mitov and Rene Poncelet,
0003 /// as described in https://arxiv.org/pdf/2205.11879 (v2)
0004 ///
0005 /// Authors: Michal Czakon, Alexander Mitov, and Rene Poncelet
0006 ///
0007 /// based on initial version by
0008 ///
0009 /// Authors: Fabrizio Caola, Radoslaw Grabarczyk, Maxwell Hutt,
0010 ///          Gavin P. Salam, Ludovic Scyboz, and Jesse Thaler
0011 ///
0012
0013 #ifndef __CMPPLUGIN_HH__
0014 #define __CMPPLUGIN_HH__
0015
0016 #include "fastjet/contrib/FlavInfo.hh"
0017 // to facilitate use with fjcore
0018 #ifndef __FJC_FLAVINFO_USEFJCORE__
0019 #include "fastjet/NNH.hh"
0020 #endif
0021
0022 FASTJET_BEGIN_NAMESPACE  // defined in fastjet/internal/base.hh
0023 using namespace std;
0024 using namespace contrib;
0025
0026 /// Plugin for algorithm  proposed by Czakon, Mitov and Poncelet
0027 /// as described in arXiv:2205.11879 with the modification
0028 class CMPPlugin : public JetDefinition::Plugin {
0029 public:
0030
0031  /// correction to original CMP
0032  enum CorrectionType {
0033   /// this uses the CMP distance measure factor as written in 2205.11879v1
0034   /// Eq. 2.9:
0035   ///
0036   ///     Sij = 1 - Theta(1-κ)cos(πκ/2)
0037   NoCorrection,
0038   /// SqrtCoshyCosPhiArgument causes kappa in Eq.(2.9) to be
0039   /// mutiplied by
0040   ///
0041   ///    sqrt(M), with M = 2*(cosh(drap) - cos(dphi))/deltaR2
0042   ///
0043   /// which tends to 1 in the small deltaR^2 limit, and fixes
0044   /// a joint infrared/collinear safety issue when deltaR is large
0045   ///
0046   SqrtCoshyCosPhiArgument,
0047   /// Default: same as above, but with damping a = 2
0048   ///
0049   ///  sqrt(M), with M = 2*(1/a^2*[cosh(a*drap)-1] - [cos(dphi)-1])/deltaR2
0050   ///
0051   SqrtCoshyCosPhiArgument_a2,
0052   // NOT RECOMMENDED: with M from above, uses
0053   //
0054   //     Sij = 1 - Theta(1-κ)cos(πκ/2) * M
0055   CoshyCosPhi,
0056   // with M from above, uses
0057   //
0058   //     Sij = (1 - Theta(1-κ)cos(πκ/2)) * M
0059   OverAllCoshyCosPhi,
0060   // same as above, but with damping a = 2
0061   OverAllCoshyCosPhi_a2
0062  };
0063
0064  /// definition of ktmax
0065  ///   0 : dynamic ktmax (pT of the hardest pseudojet lying around, default)
0066  ///   1 : fixed ktmax   (pT of the hardest pseudojet clustered with anti-kt)
0067  enum ClusteringType {DynamicKtMax, FixedKtMax};
0068
0069   /// Main constructor for the class
0070   ///
0071   /// @param R: the jet radius parameter
0072   ///
0073   /// @param a: CMP parameter a
0074   ///
0075   /// @param correction_type: the type of correction we use in the CMP distance
0076   ///
0077   /// @param clustering_type: the definition of ktmax in the CMP distance
0078   ///
0079   /// @param spherical: true if using the e+e- version of the algorithm
0080  CMPPlugin(double R, double a,
0081            CorrectionType correction_type = SqrtCoshyCosPhiArgument_a2,
0082            ClusteringType clustering_type = DynamicKtMax, bool spherical=false)
0083      : _R(R), _a(a), _correction_type(correction_type),
0084        _clustering_type(clustering_type), _spherical(spherical) {}
0085
0086  /// copy constructor
0087  CMPPlugin(const CMPPlugin &plugin) { *this = plugin; }
0088
0089  /// CMP parameter
0090  double a() const { return _a; }
0091  /// Jet radius of antikt distance
0092  double R() const { return _R; }
0093
0094  /// whether to use the e+e- version of the algo
0095  bool is_spherical() const { return _spherical; }
0096
0097  // Required by base class:
0098  virtual std::string description() const;
0099
0100  virtual double precise_squared_distance(const PseudoJet & j1, const PseudoJet & j2) const;
0101
0102  virtual double distance_opposite_flavour(const PseudoJet & j1, const PseudoJet & j2,
0103                                           const double ktmax) const;
0104
0105  virtual void run_clustering(ClusterSequence &) const;
0106
0107  void cross_product(const PseudoJet & p1, const PseudoJet & p2,
0108                     PseudoJet & retp, bool lightlike=false) const {
0109    double px = p1.py() * p2.pz() - p2.py() * p1.pz();
0110    double py = p1.pz() * p2.px() - p2.pz() * p1.px();
0111    double pz = p1.px() * p2.py() - p2.px() * p1.py();
0112
0113    double E;
0114    if (lightlike) {
0115      E = sqrt(px*px + py*py + pz*pz);
0116    } else {
0117      E = 0.0;
0118    }
0119
0120    // is this exactly what we want? The masses seem to change, need to check!
0121    retp.reset(px, py, pz, E);
0122
0123    return;
0124  }
0125  double dot_product_3d(const PseudoJet & a, const PseudoJet & b) const {
0126    return a.px()*b.px() + a.py()*b.py() + a.pz()*b.pz();
0127  }
0128  /// Returns (1-cos theta) where theta is the angle between p1 and p2
0129  double one_minus_costheta(const PseudoJet & p1, const PseudoJet & p2) const {
0130
0131    if (p1.m2() == 0 && p2.m2() == 0) {
0132      // use the 4-vector dot product.
0133      // For massless particles it gives us E1*E2*(1-cos theta)
0134      double res = dot_product(p1,p2) / (p1.E() * p2.E());
0135      return res;
0136    } else {
0137      double p1mod = sqrt(p1.modp2());
0138      double p2mod = sqrt(p2.modp2());
0139      double p1p2mod = p1mod*p2mod;
0140      double dot = dot_product_3d(p1,p2);
0141
0142      if (dot > (1-std::numeric_limits<double>::epsilon()) * p1p2mod) {
0143        PseudoJet cross_result;
0144        cross_product(p1, p2, cross_result, false);
0145        // the mass^2 of cross_result is equal to
0146        // -(px^2 + py^2 + pz^2) = (p1mod*p2mod*sintheta_ab)^2
0147        // so we can get
0148        double res = -cross_result.m2()/(p1p2mod * (p1p2mod+dot));
0149
0150        return res;
0151      }
0152      return 1.0 - dot/p1p2mod;
0153
0154    }
0155  }
0156
0157 private:
0158   static const double _deltaR2_handover;
0159   double _R, _a;
0160   CorrectionType _correction_type;
0161   ClusteringType _clustering_type;
0162   bool _spherical;
0163   template <typename NN>
0164   void _NN_clustering(ClusterSequence &cs, NN &nn) const;
0165   //template <typename NN, typename CMPNNInfo>
0166   //void k_tmax_clustering(ClusterSequence &cs, NN &nn, CMPNNInfo info) const;
0167 };
0168
0169 FASTJET_END_NAMESPACE
0170
0171 #endif  // __CMPPLUGIN_HH__