fastjet/contrib/KTClusCXX.hh

0001 // KTClusCXX Package
0002 //
0003 // Copyright (c) 2021 Andrii Verbytskyi <andrii.verbytskyi@mpp.mpg.de>
0004 //
0005 // $Id$
0006 //----------------------------------------------------------------------
0007 //    The package aims to implement the functionality of KTCLUS from the
0008 //    'ktclus.f' package by M.Seymour, 1992-2010, which was widely
0009 //    used in the past. Therefore, the KTClusCXXJetDefinition function
0010 //    attempts to mimic the interface of the KTCLUS and accepts a single
0011 //    4-digit integer to create a jet algorithm definition. The integer
0012 //    is decoded as follows (the definition is taken from the 'ktclus.f'):
0013 //
0014 //    The collision type and analysis type are indicated by the first
0015 //    argument of KTCLUS. IMODE=<TYPE><ANGLE><MONO><RECOM> where
0016 //    TYPE:  1=>ee, 2=>ep with p in -z direction, 3=>pe, 4=>pp
0017 //    ANGLE: 1=>angular kt def., 2=>DeltaR, 3=>f(DeltaEta,DeltaPhi)
0018 //           where f()=2(cosh(eta)-cos(phi)) is the QCD emission metric
0019 //    MONO:  1=>derive relative pseudoparticle angles from jets
0020 //           2=>monotonic definitions of relative angles
0021 //    RECOM: 1=>E recombination scheme, 2=>pt scheme, 3=>pt**2 scheme
0022 //
0023 //----------------------------------------------------------------------
0024 // This file is part of FastJet contrib.
0025 //
0026 // It is free software; you can redistribute it and/or modify it under
0027 // the terms of the GNU General Public License as published by the
0028 // Free Software Foundation; either version 2 of the License, or (at
0029 // your option) any later version.
0030 //
0031 // It is distributed in the hope that it will be useful, but WITHOUT
0032 // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
0033 // or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
0034 // License for more details.
0035 //
0036 // You should have received a copy of the GNU General Public License
0037 // along with this code. If not, see <http://www.gnu.org/licenses/>.
0038 //----------------------------------------------------------------------
0039
0040 #ifndef __FASTJET_CONTRIB_KTCLUSCXX_HH__
0041 #define __FASTJET_CONTRIB_KTCLUSCXX_HH__
0042
0043 #include <iostream>
0044 #include <limits>
0045 #include <vector>
0046 #include <sstream>
0047 #include <string>
0048 #include "fastjet/PseudoJet.hh"
0049 #include "fastjet/JetDefinition.hh"
0050 #include "fastjet/NNH.hh"
0051
0052 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
0053 // forward declaration to reduce includes
0054 class ClusterSequence;
0055
0056 namespace contrib {
0057
0058
0059 namespace KTClusCXX {
0060    enum Type {
0061        ee = 1000,
0062        ep = 2000,
0063        pe = 3000,
0064        pp = 4000
0065      };
0066    enum Angle {
0067        angular = 100,
0068        deltar = 200,
0069        f = 300
0070      };
0071    enum Mono {
0072        jets = 10,
0073        monolitic = 20
0074      };
0075    enum Recom {
0076        e = 1,
0077        pt =2,
0078        pt2 =3
0079      };
0080 };
0081
0082 class KTClusInfo {
0083 public:
0084     KTClusInfo(const double R = 1.0): _R(R) {}
0085     double R() const {
0086         return _R;
0087     }
0088 private:
0089     double _R;
0090 };
0091
0092 template<int TYPE, int ANGLE, int MONO, int RECOM>
0093 class KTClusBriefJet {
0094 public:
0095     inline double to_pi(const double& a1) const {
0096       double  phi = a1;
0097       while (phi>M_PI) phi -= (2*M_PI);
0098       while (phi<-M_PI) phi += (2*M_PI);
0099       return phi;
0100     }
0101     void init(const PseudoJet & jet, KTClusInfo* IN) {
0102         double norm = 1.0/std::sqrt(jet.modp2());
0103         nx = jet.px() * norm;
0104         ny = jet.py() * norm;
0105         nz = jet.pz() * norm;
0106         E = jet.E();
0107         pt = jet.pt();
0108         phi = jet.phi();
0109         eta = jet.eta();
0110         rap = jet.rap();
0111         R = IN->R();
0112     }
0113
0114
0115     double distance(const KTClusBriefJet * jet) const {
0116         /// For a better readability we have a single line for each set of template arguments.
0117         if (MONO != 10) throw Error("The algorithms with MONO!=1 are not implemented in the current version.");
0118
0119         if (TYPE == 1000 && ANGLE == 100) return 2*std::min(E*E, jet->E*jet->E)*(1.0 - (nx*jet->nx + ny*jet->ny + nz*jet->nz));
0120         if (TYPE == 2000 && ANGLE == 100) return 2*std::min(E*E, jet->E*jet->E)*(1.0 - (nx*jet->nx + ny*jet->ny + nz*jet->nz));
0121         if (TYPE == 3000 && ANGLE == 100) return 2*std::min(E*E, jet->E*jet->E)*(1.0 - (nx*jet->nx + ny*jet->ny + nz*jet->nz));
0122         if (TYPE == 4000 && ANGLE == 100) return 2*std::min(E*E, jet->E*jet->E)*(1.0 - (nx*jet->nx + ny*jet->ny + nz*jet->nz));
0123
0124         if (TYPE == 1000 && ANGLE == 200) return std::min(pt*pt, jet->pt*jet->pt)*(std::pow(to_pi(phi - jet->phi), 2)  + std::pow(rap - jet->rap, 2));
0125         if (TYPE == 2000 && ANGLE == 200) return std::min(pt*pt, jet->pt*jet->pt)*(std::pow(to_pi(phi - jet->phi), 2)  + std::pow(rap - jet->rap, 2));
0126         if (TYPE == 3000 && ANGLE == 200) return std::min(pt*pt, jet->pt*jet->pt)*(std::pow(to_pi(phi - jet->phi), 2)  + std::pow(rap - jet->rap, 2));
0127         if (TYPE == 4000 && ANGLE == 200) return std::min(pt*pt, jet->pt*jet->pt)*(std::pow(to_pi(phi - jet->phi), 2)  + std::pow(rap - jet->rap, 2));
0128
0129         if (TYPE == 1000 && ANGLE == 300) return std::min(pt*pt, jet->pt*jet->pt)*2*(std::cosh(rap - jet->rap) - std::cos(to_pi(phi - jet->phi)));
0130         if (TYPE == 2000 && ANGLE == 300) return std::min(pt*pt, jet->pt*jet->pt)*2*(std::cosh(rap - jet->rap) - std::cos(to_pi(phi - jet->phi)));
0131         if (TYPE == 3000 && ANGLE == 300) return std::min(pt*pt, jet->pt*jet->pt)*2*(std::cosh(rap - jet->rap) - std::cos(to_pi(phi - jet->phi)));
0132         if (TYPE == 4000 && ANGLE == 300) return std::min(pt*pt, jet->pt*jet->pt)*2*(std::cosh(rap - jet->rap) - std::cos(to_pi(phi - jet->phi)));
0133
0134         throw Error("Unknown distance function");
0135
0136     }
0137
0138     double beam_distance() const {
0139         /// For a better readability we have a single line for each set of template arguments.
0140         if (TYPE == 1000 && ANGLE == 100) return std::numeric_limits<double>::max()/256;
0141         if (TYPE == 2000 && ANGLE == 100) return 2*R*R*E*E*(1 + nz);
0142         if (TYPE == 3000 && ANGLE == 100) return 2*R*R*E*E*(1 - nz);
0143         if (TYPE == 4000 && ANGLE == 100) return 2*R*R*E*E*std::min(1 - nz, 1 + nz);
0144
0145         if (TYPE == 1000 && ANGLE == 200) return std::numeric_limits<double>::max()/256;
0146         if (TYPE == 2000 && ANGLE == 200) return R*R*pt*pt;
0147         if (TYPE == 3000 && ANGLE == 200) return R*R*pt*pt;
0148         if (TYPE == 4000 && ANGLE == 200) return R*R*pt*pt;
0149
0150         if (TYPE == 1000 && ANGLE == 300) return std::numeric_limits<double>::max()/256;
0151         if (TYPE == 2000 && ANGLE == 300) return R*R*pt*pt;
0152         if (TYPE == 3000 && ANGLE == 300) return R*R*pt*pt;
0153         if (TYPE == 4000 && ANGLE == 300) return R*R*pt*pt;
0154
0155         throw Error("Unknown beam distance function");
0156     }
0157
0158 private:
0159     double E, pt, eta, rap, phi, nx, ny, nz, R;
0160 };
0161
0162
0163 template<int TYPE, int ANGLE, int MONO, int RECOM>
0164 class KTClusCXXPlugin : public JetDefinition::Plugin {
0165 public:
0166     enum Strategy { strategy_NNH = 0};
0167
0168     KTClusCXXPlugin (const double radius = 1.0, Strategy strategy = strategy_NNH): _R(radius),  _strategy(strategy) {}
0169
0170     KTClusCXXPlugin (const KTClusCXXPlugin & plugin) {
0171         *this = plugin;
0172     }
0173
0174     virtual std::string description () const {
0175         std::ostringstream desc;
0176         desc << " KTCLUS algorithm plugin in mode " << TYPE << ANGLE << MONO << RECOM;
0177         switch(_strategy) {
0178         case strategy_NNH:
0179             desc << ", using NNH strategy";
0180             break;
0181         default:
0182             throw Error("Unrecognized strategy in KTClusCXX");
0183         }
0184
0185         return desc.str();
0186     }
0187
0188     virtual void run_clustering(ClusterSequence & cs) const {
0189         switch(_strategy) {
0190         case strategy_NNH:
0191             _actual_run_clustering<NNH<KTClusBriefJet < TYPE, ANGLE, MONO, 0 >,KTClusInfo> >(cs);
0192             break;
0193         default:
0194             throw Error("Unrecognized strategy in KTClusCXX");
0195         }
0196     }
0197
0198     virtual double R() const {
0199         return _R;
0200     }
0201
0202     virtual bool exclusive_sequence_meaningful() const {
0203         return true;
0204     }
0205
0206 private:
0207     double _R;
0208     Strategy _strategy;
0209     template<class N> void _actual_run_clustering(ClusterSequence & cs) const {
0210         int njets = cs.jets().size();
0211         KTClusInfo  IN(_R);
0212         N nn(cs.jets(),&IN);
0213         while (njets > 0) {
0214             int i, j, k;
0215             double dij = -1;
0216             dij= nn.dij_min(i, j);
0217             if (j >= 0) {
0218                 cs.plugin_record_ij_recombination(i, j, dij, k);
0219                 nn.merge_jets(i, j, cs.jets()[k], k);
0220             } else {
0221                 KTClusBriefJet < TYPE, ANGLE, MONO, 0 > jt;
0222                 jt.init(cs.jets()[i],&IN);
0223                 double diB = jt.beam_distance();
0224                 cs.plugin_record_iB_recombination(i, diB);
0225                 nn.remove_jet(i);
0226             }
0227             njets--;
0228         }
0229     }
0230 };
0231
0232 template<int MODE>
0233 JetDefinition KTClusCXXJetDefinition(const double Rad = 1.0) {
0234     RecombinationScheme rs;
0235     if (MODE%10 == 1) rs = E_scheme;
0236     if (MODE%10 == 2) rs = pt_scheme;
0237     if (MODE%10 == 3) rs = pt2_scheme;
0238     JetDefinition JD(new KTClusCXXPlugin<1000*(MODE/1000), 100*((MODE/100)%10), 10*((MODE/10)%10), MODE%10>(Rad));
0239     JD.set_recombination_scheme(rs);
0240     JD.delete_plugin_when_unused();
0241     return JD;
0242 }
0243
0244 inline JetDefinition KTClusCXXJetDefinition(const int m, const double x = 1.0) {
0245     switch (m) {
0246     case 1111: return KTClusCXXJetDefinition<1111>(x); break;
0247     case 2111: return KTClusCXXJetDefinition<2111>(x); break;
0248     case 3111: return KTClusCXXJetDefinition<3111>(x); break;
0249     case 4111: return KTClusCXXJetDefinition<4111>(x); break;
0250
0251     case 1211: return KTClusCXXJetDefinition<1211>(x); break;
0252     case 2211: return KTClusCXXJetDefinition<2211>(x); break;
0253     case 3211: return KTClusCXXJetDefinition<3211>(x); break;
0254     case 4211: return KTClusCXXJetDefinition<4211>(x); break;
0255
0256     case 1311: return KTClusCXXJetDefinition<1311>(x); break;
0257     case 2311: return KTClusCXXJetDefinition<2311>(x); break;
0258     case 3311: return KTClusCXXJetDefinition<3311>(x); break;
0259     case 4311: return KTClusCXXJetDefinition<4311>(x); break;
0260
0261     case 1121: return KTClusCXXJetDefinition<1121>(x); break;
0262     case 2121: return KTClusCXXJetDefinition<2121>(x); break;
0263     case 3121: return KTClusCXXJetDefinition<3121>(x); break;
0264     case 4121: return KTClusCXXJetDefinition<4121>(x); break;
0265
0266     case 1221: return KTClusCXXJetDefinition<1221>(x); break;
0267     case 2221: return KTClusCXXJetDefinition<2221>(x); break;
0268     case 3221: return KTClusCXXJetDefinition<3221>(x); break;
0269     case 4221: return KTClusCXXJetDefinition<4221>(x); break;
0270
0271     case 1321: return KTClusCXXJetDefinition<1321>(x); break;
0272     case 2321: return KTClusCXXJetDefinition<2321>(x); break;
0273     case 3321: return KTClusCXXJetDefinition<3321>(x); break;
0274     case 4321: return KTClusCXXJetDefinition<4321>(x); break;
0275
0276     case 1112: return KTClusCXXJetDefinition<1112>(x); break;
0277     case 2112: return KTClusCXXJetDefinition<2112>(x); break;
0278     case 3112: return KTClusCXXJetDefinition<3112>(x); break;
0279     case 4112: return KTClusCXXJetDefinition<4112>(x); break;
0280
0281     case 1212: return KTClusCXXJetDefinition<1212>(x); break;
0282     case 2212: return KTClusCXXJetDefinition<2212>(x); break;
0283     case 3212: return KTClusCXXJetDefinition<3212>(x); break;
0284     case 4212: return KTClusCXXJetDefinition<4212>(x); break;
0285
0286     case 1312: return KTClusCXXJetDefinition<1312>(x); break;
0287     case 2312: return KTClusCXXJetDefinition<2312>(x); break;
0288     case 3312: return KTClusCXXJetDefinition<3312>(x); break;
0289     case 4312: return KTClusCXXJetDefinition<4312>(x); break;
0290
0291     case 1122: return KTClusCXXJetDefinition<1122>(x); break;
0292     case 2122: return KTClusCXXJetDefinition<2122>(x); break;
0293     case 3122: return KTClusCXXJetDefinition<3122>(x); break;
0294     case 4122: return KTClusCXXJetDefinition<4122>(x); break;
0295
0296     case 1222: return KTClusCXXJetDefinition<1222>(x); break;
0297     case 2222: return KTClusCXXJetDefinition<2222>(x); break;
0298     case 3222: return KTClusCXXJetDefinition<3222>(x); break;
0299     case 4222: return KTClusCXXJetDefinition<4222>(x); break;
0300
0301     case 1322: return KTClusCXXJetDefinition<1322>(x); break;
0302     case 2322: return KTClusCXXJetDefinition<2322>(x); break;
0303     case 3322: return KTClusCXXJetDefinition<3322>(x); break;
0304     case 4322: return KTClusCXXJetDefinition<4322>(x); break;
0305
0306     case 1113: return KTClusCXXJetDefinition<1113>(x); break;
0307     case 2113: return KTClusCXXJetDefinition<2113>(x); break;
0308     case 3113: return KTClusCXXJetDefinition<3113>(x); break;
0309     case 4113: return KTClusCXXJetDefinition<4113>(x); break;
0310
0311     case 1213: return KTClusCXXJetDefinition<1213>(x); break;
0312     case 2213: return KTClusCXXJetDefinition<2213>(x); break;
0313     case 3213: return KTClusCXXJetDefinition<3213>(x); break;
0314     case 4213: return KTClusCXXJetDefinition<4213>(x); break;
0315
0316     case 1313: return KTClusCXXJetDefinition<1313>(x); break;
0317     case 2313: return KTClusCXXJetDefinition<2313>(x); break;
0318     case 3313: return KTClusCXXJetDefinition<3313>(x); break;
0319     case 4313: return KTClusCXXJetDefinition<4313>(x); break;
0320
0321     case 1123: return KTClusCXXJetDefinition<1123>(x); break;
0322     case 2123: return KTClusCXXJetDefinition<2123>(x); break;
0323     case 3123: return KTClusCXXJetDefinition<3123>(x); break;
0324     case 4123: return KTClusCXXJetDefinition<4123>(x); break;
0325
0326     case 1223: return KTClusCXXJetDefinition<1223>(x); break;
0327     case 2223: return KTClusCXXJetDefinition<2223>(x); break;
0328     case 3223: return KTClusCXXJetDefinition<3223>(x); break;
0329     case 4223: return KTClusCXXJetDefinition<4223>(x); break;
0330
0331     case 1323: return KTClusCXXJetDefinition<1323>(x); break;
0332     case 2323: return KTClusCXXJetDefinition<2323>(x); break;
0333     case 3323: return KTClusCXXJetDefinition<3323>(x); break;
0334     case 4323: return KTClusCXXJetDefinition<4323>(x); break;
0335     default:
0336     throw Error(std::string("In KTClusCXXJetDefinition, unrecognized mode =")+ std::to_string(m));
0337     }
0338     return JetDefinition();
0339 }
0340 inline JetDefinition KTClusCXXJetDefinition(const int m1,const int m2,const int m3,const int m4, const double x = 1.0)
0341 { return KTClusCXXJetDefinition(m1+m2+m3+m4, x);}
0342
0343 }  // namespace contrib
0344 FASTJET_END_NAMESPACE        // defined in fastjet/internal/base.hh
0345
0346 #endif // __KTCLUSCXX_HH__
0347