fastjet/contrib/FlavNeutraliser.hh

0001 #ifndef __FJCONTRIB_FLAVNEUTRALISER_HH__
0002 #define __FJCONTRIB_FLAVNEUTRALISER_HH__
0003
0004 #include "fastjet/contrib/FlavInfo.hh"
0005
0006 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
0007 namespace contrib{
0008
0009 // ----------------------------------------------------------------------
0010 /// Compare the net flavour of two vector<fastjet::PseudoJet> objects.
0011 bool jet_net_flavour_compare(std::vector<fastjet::PseudoJet> & j, std::vector<fastjet::PseudoJet> & k);
0012
0013 /// Compare the flavour of each jet of two vector<fastjet::PseudoJet> objects.
0014 /// If either vector has >3 jets, only the 3 hardest (i.e. largest energy) jets
0015 /// are compared.
0016 bool jet_flavour_compare(const std::vector<fastjet::PseudoJet> & j, const std::vector<fastjet::PseudoJet> & k,
0017                         const int max_jets=-1, bool sort_by_px = false);
0018
0019
0020 /// ----------------------------------------------------------------
0021 /// Class for the new flavour neutraliser
0022 class FlavNeutraliser {
0023 public:
0024   enum measure { sinh_delta_R, delta_R, jade_delta_R, maxscale_delta_R, phi2_coshy, cosphi_coshy,
0025                  aktlike_pair_refratio, aktlike_pair_dynrefratio,
0026                  jade, jadea2, maxscale, general };
0027
0028   /// Main constructor for the FlavNeutraliser class.
0029   FlavNeutraliser(bool modulo_2,
0030                   measure measure_in,
0031                   bool use_mass_flav = false,
0032                   bool spherical_algo = false,
0033                   double pp = 1.0,
0034                   bool recursive_neutralisation = true) :
0035                   _modulo_2(modulo_2), _measure(measure_in), _use_mass_flav(use_mass_flav),
0036                   _spherical_algo(spherical_algo), _pp(pp), _recursive(recursive_neutralisation)
0037                   {}
0038
0039   /// A constructor for the FlavNeutraliser class for
0040   /// the general measure of the form:
0041   ///
0042   ///   u_ij = max( pti^2 , ptj^2 )^p min( pti^2 , ptj^2 )^q
0043   ///          x 2[ 1/a^2 (cosh(a*Δy_ij) - 1) + (cos(Δφ_ij) - 1) ]
0044   FlavNeutraliser(double p, double q, double a,
0045                   bool modulo_2,
0046                   measure measure_in = general,
0047                   bool use_mass_flav = false,
0048                   bool spherical_algo = false,
0049                   double pp = 1.0,
0050                   bool recursive_neutralisation = true) :
0051                   _p(p), _q(q), _a(a), _modulo_2(modulo_2), _measure(measure_in),
0052                   _use_mass_flav(use_mass_flav), _spherical_algo(spherical_algo),
0053                   _pp(pp), _recursive(recursive_neutralisation)
0054                   {}
0055
0056   void set_recursive(bool val) {_recursive = val;}
0057   bool recursive() const {return _recursive;}
0058
0059   /// @brief Returns the neutralisation distance (u12) between the
0060   ///        two input jets
0061   /// @param p1 input pseudojet 1
0062   /// @param p2 input pseudojet 1
0063   /// @param ref_scale reference scale for distance measures that need such a scale
0064   /// @return u12
0065   double neutralisation_distance(const PseudoJet &p1,
0066                                  const PseudoJet &p2, double ref_scale) const;
0067
0068   /// Uses the ClusterSequence to build up neutralisation and returns a
0069   /// flavour-neutralised vector<PseudoJet> that is intended to replace
0070   /// cs.jets() (which contains the whole clustering history).
0071   std::vector<PseudoJet> neutralise(ClusterSequence &cs) const;
0072
0073   /// returns the inclusive jets as obtained after flavour
0074   /// neutralisation (i.e. with the neutralise function).
0075   std::vector<PseudoJet> inclusive_jets(ClusterSequence &cs,
0076                                    const double ptmin = 0.0) const {
0077     const std::vector<ClusterSequence::history_element> & hist = cs.history();
0078     std::vector<PseudoJet> jets = neutralise(cs);
0079     std::vector<PseudoJet> jets_local;
0080     int u = int(hist.size()) - 1;
0081     while (u >= 0) {
0082       if (hist[u].parent2 == ClusterSequence::BeamJet) {
0083         int parent1 = hist[u].parent1;
0084         const PseudoJet & jet = jets[hist[parent1].jetp_index];
0085         if (jet.pt() >= ptmin) {jets_local.push_back(jet);}
0086       }
0087       u--;
0088     }
0089     return jets_local;
0090   }
0091
0092   void use_neutralisation_candidates(
0093     PseudoJet & jet_i,
0094     double uij,
0095     int ih_step,
0096     std::vector<std::pair<PseudoJet*, double>> & flavour_candidates,
0097     double ref_scale) const;
0098
0099   void use_neutralisation_candidates_recursive(
0100     PseudoJet & jet_i,
0101     double uij,
0102     int ih_step,
0103     std::vector<std::pair<PseudoJet*, double>> & flavour_candidates,
0104     double ref_scale,
0105     const PseudoJet * exclude = nullptr
0106     ) const;
0107
0108   /// @brief  return a value for the jet's rapidity
0109   ///         that is FJ's default except at small rapidity
0110   ///         where a more accurate formula is used, allowing
0111   ///         representations of rapidities down to the
0112   ///         numeric_limits<double>::min()
0113   ///
0114   /// @param p the input jet
0115   /// @return the input jet's rapidity
0116   static double accurate_rap(const PseudoJet & p) {
0117     double rap = p.rap();
0118
0119     // switch to accurate rap if the momentum has moderately small
0120     // rapidity
0121     constexpr double rap_transition = 0.1;
0122     if (std::fabs(rap) < rap_transition) rap = 0.5 * log1p(2*p.pz()/(p.E() - p.pz()));
0123     return rap;
0124   }
0125
0126   /// @brief  return a value for the absdphi between p1 and p2
0127   /// @param p1 input jet 1
0128   /// @param p2 input jet 2
0129   /// @return abs(dphi_{12})
0130   ///
0131   /// the result is FJ's default, except at small dphi values, where a
0132   /// more accurate calculation is used based on a transverse cross
0133   /// product, which is more accurate in particular when both momenta
0134   /// are close to either the x or y axis
0135   static double accurate_absdphi(const PseudoJet & p1, const PseudoJet & p2) {
0136     double dphi = std::fabs(p1.delta_phi_to(p2));
0137     constexpr double phi_transition = 0.1;
0138     if (dphi < phi_transition) {
0139       // use a more accurate calculation based on a transverse cross product;
0140       // organise the cross product so as to avoid excessively large or small
0141       // numbers appearing in intermediate stages
0142       double invp1pt = 1.0/p1.pt();
0143       double invp2pt = 1.0/p2.pt();
0144       double cross = (p1.px()*invp1pt) * (p2.py()*invp2pt) - (p2.px()*invp2pt) * (p1.py()*invp1pt);
0145       //double cross = (p1.px() * p2.py() - p2.px() * p1.py())/sqrt(p1.pt2() * p2.pt2());
0146       assert(cross <= 1.0 && cross >= -1.0);
0147       dphi = std::fabs(asin(cross));
0148     }
0149     return dphi;
0150   }
0151
0152 private:
0153   /// the value of deltaR2 below which we replace 2*(cosh-cos) with
0154   /// (addition as of 2021-09-26, but not yet being used)
0155   double _p, _q, _a;
0156   static const double _deltaR2_handover;
0157   bool    _modulo_2;
0158   measure _measure;
0159   bool    _use_mass_flav;
0160   bool    _spherical_algo;
0161   double  _pp;
0162   bool    _recursive;
0163
0164   bool    _writeout_uijs = false;
0165
0166
0167 };
0168
0169 } // namespace contrib
0170
0171 FASTJET_END_NAMESPACE        // defined in fastjet/internal/base.hh
0172 #endif // __FJCONTRIB_FLAVNEUTRALISER_HH__