Index: contrib/contribs/KTClusCXX/trunk/KTClusCXX.hh =================================================================== --- contrib/contribs/KTClusCXX/trunk/KTClusCXX.hh (revision 1333) +++ contrib/contribs/KTClusCXX/trunk/KTClusCXX.hh (revision 1334) @@ -1,296 +1,320 @@ // KTClusCXX Package // // Copyright (c) 2021 Andrii Verbytskyi // // $Id$ +//---------------------------------------------------------------------- +// The package aims to implement the functionality of KTCLUS from the +// 'ktclus.f' package by M.Seymour, 1992-2010, which was widely +// used in the past. Therefore, the KTClusCXXJetDefinition function +// attempts to mimic the interface of the KTCLUS and accepts a single +// 4-digin integer to create a jet algorithm definition. The integer +// is decoded as follows (the definition is taken from the 'ktclus.f'): +// +// The collision type and analysis type are indicated by the first +// argument of KTCLUS. IMODE= where +// TYPE: 1=>ee, 2=>ep with p in -z direction, 3=>pe, 4=>pp +// ANGLE: 1=>angular kt def., 2=>DeltaR, 3=>f(DeltaEta,DeltaPhi) +// where f()=2(cosh(eta)-cos(phi)) is the QCD emission metric +// MONO: 1=>derive relative pseudoparticle angles from jets +// 2=>monotonic definitions of relative angles +// RECOM: 1=>E recombination scheme, 2=>pt scheme, 3=>pt**2 scheme // //---------------------------------------------------------------------- // This file is part of FastJet contrib. // // It is free software; you can redistribute it and/or modify it under // the terms of the GNU General Public License as published by the // Free Software Foundation; either version 2 of the License, or (at // your option) any later version. // // It is distributed in the hope that it will be useful, but WITHOUT // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY // or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public // License for more details. // // You should have received a copy of the GNU General Public License // along with this code. If not, see . //---------------------------------------------------------------------- #ifndef __FASTJET_CONTRIB_KTCLUSCXX_HH__ #define __FASTJET_CONTRIB_KTCLUSCXX_HH__ -#include + #include -#include +#include #include #include +#include #include "fastjet/PseudoJet.hh" #include "fastjet/JetDefinition.hh" #include "fastjet/NNH.hh" FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh // forward declaration to reduce includes class ClusterSequence; namespace contrib { -inline double to_pi(const double& a1) { - double phi = a1; - while (phi>M_PI) phi -= (2*M_PI); - while (phi<-M_PI) phi += (2*M_PI); - return phi; -} + class KTClusInfo { public: KTClusInfo(const double R = 1.0): _R(R) {} double R() const { return _R; } private: double _R; }; template class KTClusBriefJet { public: + inline double to_pi(const double& a1) const { + double phi = a1; + while (phi>M_PI) phi -= (2*M_PI); + while (phi<-M_PI) phi += (2*M_PI); + return phi; + } void init(const PseudoJet & jet, KTClusInfo* IN) { double norm = 1.0/std::sqrt(jet.modp2()); nx = jet.px() * norm; ny = jet.py() * norm; nz = jet.pz() * norm; E = jet.E(); pt = jet.pt(); phi = jet.phi(); eta = jet.eta(); rap = jet.rap(); R = IN->R(); } double distance(const KTClusBriefJet * jet) const { + /// For a better readability we have a single line for each set of template arguments. + if (MONO != 1) throw Error("The algorithms with MONO!=1 are not implemented in the current version."); + if (TYPE == 1 && ANGLE == 1) return 2*std::min(E*E, jet->E*jet->E)*(1.0 - (nx*jet->nx + ny*jet->ny + nz*jet->nz)); if (TYPE == 2 && ANGLE == 1) return 2*std::min(E*E, jet->E*jet->E)*(1.0 - (nx*jet->nx + ny*jet->ny + nz*jet->nz)); if (TYPE == 3 && ANGLE == 1) return 2*std::min(E*E, jet->E*jet->E)*(1.0 - (nx*jet->nx + ny*jet->ny + nz*jet->nz)); if (TYPE == 4 && ANGLE == 1) return 2*std::min(E*E, jet->E*jet->E)*(1.0 - (nx*jet->nx + ny*jet->ny + nz*jet->nz)); if (TYPE == 1 && ANGLE == 2) return std::min(pt*pt, jet->pt*jet->pt)*(std::pow(to_pi(phi - jet->phi), 2) + std::pow(rap - jet->rap, 2)); if (TYPE == 2 && ANGLE == 2) return std::min(pt*pt, jet->pt*jet->pt)*(std::pow(to_pi(phi - jet->phi), 2) + std::pow(rap - jet->rap, 2)); if (TYPE == 3 && ANGLE == 2) return std::min(pt*pt, jet->pt*jet->pt)*(std::pow(to_pi(phi - jet->phi), 2) + std::pow(rap - jet->rap, 2)); if (TYPE == 4 && ANGLE == 2) return std::min(pt*pt, jet->pt*jet->pt)*(std::pow(to_pi(phi - jet->phi), 2) + std::pow(rap - jet->rap, 2)); if (TYPE == 1 && ANGLE == 3) return std::min(pt*pt, jet->pt*jet->pt)*2*(std::cosh(rap - jet->rap) - std::cos(to_pi(phi - jet->phi))); if (TYPE == 2 && ANGLE == 3) return std::min(pt*pt, jet->pt*jet->pt)*2*(std::cosh(rap - jet->rap) - std::cos(to_pi(phi - jet->phi))); if (TYPE == 3 && ANGLE == 3) return std::min(pt*pt, jet->pt*jet->pt)*2*(std::cosh(rap - jet->rap) - std::cos(to_pi(phi - jet->phi))); if (TYPE == 4 && ANGLE == 3) return std::min(pt*pt, jet->pt*jet->pt)*2*(std::cosh(rap - jet->rap) - std::cos(to_pi(phi - jet->phi))); throw Error("Unknown distance function"); + } double beam_distance() const { + /// For a better readability we have a single line for each set of template arguments. if (TYPE == 1 && ANGLE == 1) return std::numeric_limits::max()/256; if (TYPE == 2 && ANGLE == 1) return 2*R*R*E*E*(1 + nz); if (TYPE == 3 && ANGLE == 1) return 2*R*R*E*E*(1 - nz); if (TYPE == 4 && ANGLE == 1) return 2*R*R*E*E*std::min(1 - nz, 1 + nz); if (TYPE == 1 && ANGLE == 2) return std::numeric_limits::max()/256; if (TYPE == 2 && ANGLE == 2) return R*R*pt*pt; if (TYPE == 3 && ANGLE == 2) return R*R*pt*pt; if (TYPE == 4 && ANGLE == 2) return R*R*pt*pt; if (TYPE == 1 && ANGLE == 3) return std::numeric_limits::max()/256; if (TYPE == 2 && ANGLE == 3) return R*R*pt*pt; if (TYPE == 3 && ANGLE == 3) return R*R*pt*pt; if (TYPE == 4 && ANGLE == 3) return R*R*pt*pt; throw Error("Unknown beam distance function"); } private: double E, pt, eta, rap, phi, nx, ny, nz, R; }; template -class KTClusCXX : public JetDefinition::Plugin { +class KTClusCXXPlugin : public JetDefinition::Plugin { public: enum Strategy { strategy_NNH = 0}; - KTClusCXX (const double radius = 1.0, Strategy strategy = strategy_NNH): _R(radius), _strategy(strategy) {} + KTClusCXXPlugin (const double radius = 1.0, Strategy strategy = strategy_NNH): _R(radius), _strategy(strategy) {} - KTClusCXX (const KTClusCXX & plugin) { + KTClusCXXPlugin (const KTClusCXXPlugin & plugin) { *this = plugin; } virtual std::string description () const { std::ostringstream desc; desc << " KTCLUS algorithm plugin in mode " << TYPE << ANGLE << MONO << RECOM; switch(_strategy) { case strategy_NNH: desc << ", using NNH strategy"; break; default: throw Error("Unrecognized strategy in KTClusCXX"); } return desc.str(); } virtual void run_clustering(ClusterSequence & cs) const { switch(_strategy) { case strategy_NNH: _actual_run_clustering,KTClusInfo> >(cs); break; default: throw Error("Unrecognized strategy in KTClusCXX"); } } virtual double R() const { return _R; } virtual bool exclusive_sequence_meaningful() const { return true; } private: double _R; Strategy _strategy; template void _actual_run_clustering(ClusterSequence & cs) const { int njets = cs.jets().size(); KTClusInfo IN(_R); N nn(cs.jets(),&IN); while (njets > 0) { int i, j, k; double dij = -1; dij= nn.dij_min(i, j); if (j >= 0) { cs.plugin_record_ij_recombination(i, j, dij, k); nn.merge_jets(i, j, cs.jets()[k], k); } else { KTClusBriefJet < TYPE, ANGLE, MONO, 0 > jt; jt.init(cs.jets()[i],&IN); double diB = jt.beam_distance(); cs.plugin_record_iB_recombination(i, diB); nn.remove_jet(i); } njets--; } } }; template JetDefinition KTClusCXXJetDefinition(const double Rad = 1.0) { RecombinationScheme rs; if (MODE%10 == 1) rs = E_scheme; if (MODE%10 == 2) rs = pt_scheme; if (MODE%10 == 3) rs = pt2_scheme; - JetDefinition JD(new KTClusCXX(Rad)); + JetDefinition JD(new KTClusCXXPlugin(Rad)); JD.set_recombination_scheme(rs); + JD.delete_plugin_when_unused(); return JD; } inline JetDefinition KTClusCXXJetDefinition(const int m, const double x = 1.0) { switch (m) { case 1111: return KTClusCXXJetDefinition<1111>(x); break; case 2111: return KTClusCXXJetDefinition<2111>(x); break; case 3111: return KTClusCXXJetDefinition<3111>(x); break; case 4111: return KTClusCXXJetDefinition<4111>(x); break; case 1211: return KTClusCXXJetDefinition<1211>(x); break; case 2211: return KTClusCXXJetDefinition<2211>(x); break; case 3211: return KTClusCXXJetDefinition<3211>(x); break; case 4211: return KTClusCXXJetDefinition<4211>(x); break; case 1311: return KTClusCXXJetDefinition<1311>(x); break; case 2311: return KTClusCXXJetDefinition<2311>(x); break; case 3311: return KTClusCXXJetDefinition<3311>(x); break; case 4311: return KTClusCXXJetDefinition<4311>(x); break; case 1121: return KTClusCXXJetDefinition<1121>(x); break; case 2121: return KTClusCXXJetDefinition<2121>(x); break; case 3121: return KTClusCXXJetDefinition<3121>(x); break; case 4121: return KTClusCXXJetDefinition<4121>(x); break; case 1221: return KTClusCXXJetDefinition<1221>(x); break; case 2221: return KTClusCXXJetDefinition<2221>(x); break; case 3221: return KTClusCXXJetDefinition<3221>(x); break; case 4221: return KTClusCXXJetDefinition<4221>(x); break; case 1321: return KTClusCXXJetDefinition<1321>(x); break; case 2321: return KTClusCXXJetDefinition<2321>(x); break; case 3321: return KTClusCXXJetDefinition<3321>(x); break; case 4321: return KTClusCXXJetDefinition<4321>(x); break; case 1112: return KTClusCXXJetDefinition<1112>(x); break; case 2112: return KTClusCXXJetDefinition<2112>(x); break; case 3112: return KTClusCXXJetDefinition<3112>(x); break; case 4112: return KTClusCXXJetDefinition<4112>(x); break; case 1212: return KTClusCXXJetDefinition<1212>(x); break; case 2212: return KTClusCXXJetDefinition<2212>(x); break; case 3212: return KTClusCXXJetDefinition<3212>(x); break; case 4212: return KTClusCXXJetDefinition<4212>(x); break; case 1312: return KTClusCXXJetDefinition<1312>(x); break; case 2312: return KTClusCXXJetDefinition<2312>(x); break; case 3312: return KTClusCXXJetDefinition<3312>(x); break; case 4312: return KTClusCXXJetDefinition<4312>(x); break; case 1122: return KTClusCXXJetDefinition<1122>(x); break; case 2122: return KTClusCXXJetDefinition<2122>(x); break; case 3122: return KTClusCXXJetDefinition<3122>(x); break; case 4122: return KTClusCXXJetDefinition<4122>(x); break; case 1222: return KTClusCXXJetDefinition<1222>(x); break; case 2222: return KTClusCXXJetDefinition<2222>(x); break; case 3222: return KTClusCXXJetDefinition<3222>(x); break; case 4222: return KTClusCXXJetDefinition<4222>(x); break; case 1322: return KTClusCXXJetDefinition<1322>(x); break; case 2322: return KTClusCXXJetDefinition<2322>(x); break; case 3322: return KTClusCXXJetDefinition<3322>(x); break; case 4322: return KTClusCXXJetDefinition<4322>(x); break; case 1113: return KTClusCXXJetDefinition<1113>(x); break; case 2113: return KTClusCXXJetDefinition<2113>(x); break; case 3113: return KTClusCXXJetDefinition<3113>(x); break; case 4113: return KTClusCXXJetDefinition<4113>(x); break; case 1213: return KTClusCXXJetDefinition<1213>(x); break; case 2213: return KTClusCXXJetDefinition<2213>(x); break; case 3213: return KTClusCXXJetDefinition<3213>(x); break; case 4213: return KTClusCXXJetDefinition<4213>(x); break; case 1313: return KTClusCXXJetDefinition<1313>(x); break; case 2313: return KTClusCXXJetDefinition<2313>(x); break; case 3313: return KTClusCXXJetDefinition<3313>(x); break; case 4313: return KTClusCXXJetDefinition<4313>(x); break; case 1123: return KTClusCXXJetDefinition<1123>(x); break; case 2123: return KTClusCXXJetDefinition<2123>(x); break; case 3123: return KTClusCXXJetDefinition<3123>(x); break; case 4123: return KTClusCXXJetDefinition<4123>(x); break; case 1223: return KTClusCXXJetDefinition<1223>(x); break; case 2223: return KTClusCXXJetDefinition<2223>(x); break; case 3223: return KTClusCXXJetDefinition<3223>(x); break; case 4223: return KTClusCXXJetDefinition<4223>(x); break; case 1323: return KTClusCXXJetDefinition<1323>(x); break; case 2323: return KTClusCXXJetDefinition<2323>(x); break; case 3323: return KTClusCXXJetDefinition<3323>(x); break; case 4323: return KTClusCXXJetDefinition<4323>(x); break; default: - throw Error("Unrecognized option in KTCLUS"); + throw Error(std::string("In KTClusCXXJetDefinition, unrecognized mode =")+ std::to_string(m)); } return JetDefinition(); } } // namespace contrib FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh #endif // __KTCLUSCXX_HH__ Index: contrib/contribs/KTClusCXX/trunk/exampleall.cc =================================================================== --- contrib/contribs/KTClusCXX/trunk/exampleall.cc (revision 0) +++ contrib/contribs/KTClusCXX/trunk/exampleall.cc (revision 1334) @@ -0,0 +1,96 @@ +#include +#include +#include +#include +#include +#include +#include "KTClusCXX.hh" + +using namespace fastjet; + +class Input { +public: + Input() {}; + ~Input() {}; + const double PIMASS = 0.139; + std::vector input_particles; + std::size_t N; + void fill(const std::string& inp="") + { + input_particles.clear(); + std::size_t i = 0; + double px, py, pz, E; + + if (inp.size()) { + std::ifstream myfile; + myfile.open (inp.c_str()); + while (myfile >> px >> py >> pz >> E) { + input_particles.push_back(PseudoJet(px, py, pz, E)); + i++; + } + myfile.close(); + } else { + + while (std::cin >> px >> py >> pz >> E) { + input_particles.push_back(PseudoJet(px, py, pz, E)); + i++; + } + } + N = i; + } +}; +inline std::ostream& operator<<(std::ostream& os, const PseudoJet &j) { + os << j.px() << " " << j.py() << " " << j.pz() << " " << j.e(); + return os; +} +int main( int argc, char** argv) +{ + if (argc<2) std::cout <<"#Usage: "< " < modes { + 1111,2111,3111,4111,1211,2211,3211,4211, + 1311,2311,3311,4311,1112,2112,3112,4112, + 1212,2212,3212,4212,1312,2312,3312,4312, + 1113,2113,3113,4113,1213,2213,3213,4213, + 1313,2313,3313,4313 + }; + std::vector inclusives{true,false}; + double Radius = argc>2?std::atof(argv[2]):1.0; + double YCUT = argc>3?std::atof(argv[3]):0.0001; + double ECUT =argc>4?std::atof(argv[4]):1.0; + Input IN; + IN.fill(argc>1?argv[1]:""); + std::cout << std::fixed << std::setw(11) << std::setprecision(6); + std::cout<<"Input particles:"< fastjets; + fastjets.reserve(IN.N); + + for (const auto& inclusive: inclusives) { + for (const auto& mode: modes) { + ClusterSequence cs(IN.input_particles, contrib::KTClusCXXJetDefinition(mode, Radius)); + if (inclusive) { + fastjets = cs.inclusive_jets(); + } else { + fastjets = cs.exclusive_jets(ECUT*ECUT*YCUT); + } + std::vector < double > fastjets_dmerge; + std::sort(fastjets.begin(), fastjets.end(), [](const PseudoJet & a, const PseudoJet & b) -> bool { return a.E() > b.E(); }); + for (std::size_t i = 0; i < IN.N; i++) { + fastjets_dmerge.push_back(cs.exclusive_dmerge_max(i)); + } + + std::cout<<"Parameters: "<< mode << (inclusive?std::string("i"):std::string(""))<< " " << Radius << " " << YCUT<< " " << ECUT < +#include #include #include -#include #include #include #include "KTClusCXX.hh" using namespace fastjet; class Input { public: Input() {}; ~Input() {}; const double PIMASS = 0.139; std::vector input_particles; - size_t N; - void fill(const size_t N_ = 23) { - srand(1); - N = std::max(size_t{2},std::min(size_t{100},N_)); - for (size_t i = 0; i < N; i++) { - double px = 10 * (drand48() - 0.5); - double py = 10 * (drand48() - 0.5); - double pz = 10 * (drand48() - 0.3); - double e = std::sqrt(px*px + py*py + pz*pz + PIMASS*PIMASS); - auto pj = PseudoJet(px, py, pz, e); - input_particles.push_back(pj); - } - } + std::size_t N; void fill(const std::string& inp="") { - input_particles.clear(); - size_t i = 0; + std::size_t i = 0; double px, py, pz, E; if (inp.size()) { std::ifstream myfile; myfile.open (inp.c_str()); while (myfile >> px >> py >> pz >> E) { input_particles.push_back(PseudoJet(px, py, pz, E)); i++; } myfile.close(); } else { while (std::cin >> px >> py >> pz >> E) { input_particles.push_back(PseudoJet(px, py, pz, E)); i++; } } N = i; } - - void dump(const std:: string& out) { - std::ofstream myfile; - myfile.open (out.c_str()); - for (const auto& p: input_particles) { - myfile << " " << p.px() << " " << p.py() << " " << p.pz() << " " << p.e() << std::endl; - } - myfile.close(); - } - }; inline std::ostream& operator<<(std::ostream& os, const PseudoJet &j) { os << j.px() << " " << j.py() << " " << j.pz() << " " << j.e(); return os; } int main( int argc, char** argv) { ///2111 1.0 0.0001 1.0 - + if (argc<2) std::cout <<"#Usage: "< " <2?std::string(argv[2]):"2111"; bool inclusive = (smode.length() > 4 && smode[4] == 'i'); int mode = std::atoi(smode.substr(0, 4).c_str()); double Radius = argc>3?std::atof(argv[3]):1.0; double YCUT = argc>4?std::atof(argv[4]):0.0001; double ECUT =argc>5?std::atof(argv[5]):1.0; Input IN; IN.fill(argc>1?argv[1]:""); std::vector < PseudoJet > fastjets; fastjets.reserve(IN.N); ClusterSequence cs(IN.input_particles, contrib::KTClusCXXJetDefinition(mode, Radius)); if (inclusive) { fastjets = cs.inclusive_jets(); } else { fastjets = cs.exclusive_jets(ECUT*ECUT*YCUT); } std::vector < double > fastjets_dmerge; std::sort(fastjets.begin(), fastjets.end(), [](const PseudoJet & a, const PseudoJet & b) -> bool { return a.E() > b.E(); }); - for (size_t i = 0; i < IN.N; i++) { + for (std::size_t i = 0; i < IN.N; i++) { fastjets_dmerge.push_back(cs.exclusive_dmerge_max(i)); } std::cout << std::fixed << std::setw(11) << std::setprecision(6); std::cout<<"Input particles:"<