diff --git a/configure.ac b/configure.ac --- a/configure.ac +++ b/configure.ac @@ -1,331 +1,333 @@ ## Process this file with autoconf to produce a configure script. AC_PREREQ(2.59) AC_INIT([Rivet],[2.6.1],[rivet@projects.hepforge.org],[Rivet]) ## Check and block installation into the src/build dir if test "$prefix" = "$PWD"; then AC_MSG_ERROR([Installation into the build directory is not supported: use a different --prefix argument]) fi ## Force default prefix to have a path value rather than NONE if test "$prefix" = "NONE"; then prefix=/usr/local fi AC_CONFIG_SRCDIR([src/Core/Analysis.cc]) AC_CONFIG_HEADERS([include/Rivet/Config/DummyConfig.hh include/Rivet/Config/RivetConfig.hh include/Rivet/Config/BuildOptions.hh]) AM_INIT_AUTOMAKE([dist-bzip2 -Wall 1.10]) m4_ifdef([AM_SILENT_RULES], [AM_SILENT_RULES([yes])]) m4_ifdef([AM_PROG_AR], [AM_PROG_AR]) AC_CONFIG_MACRO_DIR([m4]) AC_SUBST(LT_OBJDIR) ## Package-specific #defines AC_DEFINE_UNQUOTED(RIVET_VERSION, "$PACKAGE_VERSION", "Rivet version string") AC_DEFINE_UNQUOTED(RIVET_NAME, "$PACKAGE_NAME", "Rivet name string") AC_DEFINE_UNQUOTED(RIVET_STRING, "$PACKAGE_STRING", "Rivet name and version string") AC_DEFINE_UNQUOTED(RIVET_TARNAME, "$PACKAGE_TARNAME", "Rivet short name string") AC_DEFINE_UNQUOTED(RIVET_BUGREPORT, "$PACKAGE_BUGREPORT", "Rivet contact email address") ## OS X AC_CEDAR_OSX ## Work out the LCG platform tag AC_LCG_TAG ## Set default compiler flags if test "x$CXXFLAGS" == "x"; then CXXFLAGS="-O2"; fi ## Compiler setup AC_LANG(C++) AC_PROG_CXX AX_CXX_COMPILE_STDCXX([11], [noext], [mandatory]) ## Store and propagate the compiler identity and flags RIVETCXX="$CXX" AC_SUBST(RIVETCXX) RIVETCXXFLAGS="$CXXFLAGS" AC_SUBST(RIVETCXXFLAGS) ## Checks for programs. AC_PROG_INSTALL AC_PROG_LN_S AC_DISABLE_STATIC AC_LIBTOOL_DLOPEN AC_PROG_LIBTOOL AC_FUNC_STRERROR_R ## YODA histogramming library AC_CEDAR_LIBRARYANDHEADERS([YODA], , , [AC_MSG_ERROR([YODA is required])]) YODABINPATH=$YODALIBPATH/../bin AC_SUBST(YODABINPATH) AC_PATH_PROG(YODACONFIG, yoda-config, [], [$YODALIBPATH/../bin:$PATH]) YODA_PYTHONPATH="" if test -f "$YODACONFIG"; then AC_MSG_CHECKING([YODA version using yoda-config]) YODA_VERSION=`$YODACONFIG --version` AC_MSG_RESULT([$YODA_VERSION]) YODA_VERSION1=[`echo $YODA_VERSION | cut -d. -f1 | sed -e 's/\([0-9]*\).*/\1/g'`] YODA_VERSION2=[`echo $YODA_VERSION | cut -d. -f2 | sed -e 's/\([0-9]*\).*/\1/g'`] YODA_VERSION3=[`echo $YODA_VERSION | cut -d. -f3 | sed -e 's/\([0-9]*\).*/\1/g'`] let YODA_VERSION_INT=YODA_VERSION1*10000+YODA_VERSION2*100+YODA_VERSION3 if test $YODA_VERSION_INT -lt 10500; then AC_MSG_ERROR([YODA version isn't sufficient: at least version 1.5.0 required]) fi AC_MSG_CHECKING([YODA Python path using yoda-config]) YODA_PYTHONPATH=`$YODACONFIG --pythonpath` AC_MSG_RESULT([$YODA_PYTHONPATH]) fi AC_SUBST(YODA_PYTHONPATH) ## HepMC event record library AC_CEDAR_LIBRARYANDHEADERS([HepMC], , , [AC_MSG_ERROR([HepMC is required])]) oldCPPFLAGS=$CPPFLAGS CPPFLAGS="$CPPFLAGS -I$HEPMCINCPATH" if test -e "$HEPMCINCPATH/HepMC/HepMCDefs.h"; then AC_LANG_CONFTEST([AC_LANG_SOURCE([#include #include "HepMC/HepMCDefs.h" int main() { std::cout << HEPMC_VERSION << std::endl; return 0; }])]) else AC_LANG_CONFTEST([AC_LANG_SOURCE([#include #include "HepMC/defs.h" int main() { std::cout << VERSION << std::endl; return 0; }])]) fi if test -f conftest.cc; then $CXX $CPPFLAGS conftest.cc -o conftest 2>&1 1>&5 elif test -f conftest.C; then $CXX $CPPFLAGS conftest.C -o conftest 2>&1 1>&5 else $CXX $CPPFLAGS conftest.cpp -o conftest 2>&1 1>&5 fi hepmc_version=`./conftest` if test x$hepmc_version != x; then let hepmc_major=`echo "$hepmc_version" | cut -d. -f1` let hepmc_minor=`echo "$hepmc_version" | cut -d. -f2` fi rm -f conftest conftest.cpp conftest.cc conftest.C HEPMC_VERSION=$hepmc_major$hepmc_minor AC_MSG_NOTICE([HepMC version is $hepmc_version -> $HEPMC_VERSION]) AC_SUBST(HEPMC_VERSION) CPPFLAGS=$oldCPPFLAGS ## FastJet clustering library AC_CEDAR_LIBRARYANDHEADERS([fastjet], , , [AC_MSG_ERROR([FastJet is required])]) AC_PATH_PROG(FJCONFIG, fastjet-config, [], $FASTJETPATH/bin:$PATH) if test -f "$FJCONFIG"; then AC_MSG_CHECKING([FastJet version using fastjet-config]) fjversion=`$FJCONFIG --version` AC_MSG_RESULT([$fjversion]) fjmajor=$(echo $fjversion | cut -f1 -d.) fjminor=$(echo $fjversion | cut -f2 -d.) fjmicro=$(echo $fjversion | cut -f3 -d.) if test "$fjmajor" -lt 3; then AC_MSG_ERROR([FastJet version 3.0.0 or later is required]) fi FASTJETCONFIGLIBADD="$($FJCONFIG --plugins --shared --libs)" else FASTJETCONFIGLIBADD="-L$FASTJETLIBPATH -l$FASTJETLIBNAME" FASTJETCONFIGLIBADD="$FASTJETCONFIGLIBADD -lSISConePlugin -lsiscone -lsiscone_spherical" FASTJETCONFIGLIBADD="$FASTJETCONFIGLIBADD -lCDFConesPlugin -lD0RunIIConePlugin -lNestedDefsPlugin" FASTJETCONFIGLIBADD="$FASTJETCONFIGLIBADD -lTrackJetPlugin -lATLASConePlugin -lCMSIterativeConePlugin" FASTJETCONFIGLIBADD="$FASTJETCONFIGLIBADD -lEECambridgePlugin -lJadePlugin" fi; AC_MSG_NOTICE([FastJet LIBADD = $FASTJETCONFIGLIBADD]) AC_SUBST(FASTJETCONFIGLIBADD) # Check for FastJet headers that require the --enable-all(cxx)plugins option FASTJET_ERRMSG="Required FastJet plugin headers were not found: did you build FastJet with the --enable-allcxxplugins option?" oldCPPFLAGS=$CPPFLAGS CPPFLAGS="$CPPFLAGS -I$FASTJETINCPATH" AC_CHECK_HEADER([fastjet/D0RunIIConePlugin.hh], [], [AC_MSG_ERROR([$FASTJET_ERRMSG])]) AC_CHECK_HEADER([fastjet/TrackJetPlugin.hh], [], [AC_MSG_ERROR([$FASTJET_ERRMSG])]) CPPFLAGS=$oldCPPFLAGS # ## GNU Scientific Library # AC_SEARCH_GSL # AC_CEDAR_HEADERS([gsl], , , [AC_MSG_ERROR([GSL (GNU Scientific Library) is required])]) # oldCPPFLAGS=$CPPFLAGS # CPPFLAGS="$CPPFLAGS -I$GSLINCPATH" # AC_CHECK_HEADER([gsl/gsl_vector.h], [], [AC_MSG_ERROR([GSL vectors not found.])]) # CPPFLAGS=$oldCPPFLAGS ## Disable build/install of standard analyses AC_ARG_ENABLE([analyses], [AC_HELP_STRING(--disable-analyses, [don't try to build or install standard analyses])], [], [enable_analyses=yes]) if test x$enable_analyses != xyes; then AC_MSG_WARN([Not building standard Rivet analyses, by request]) fi AM_CONDITIONAL(ENABLE_ANALYSES, [test x$enable_analyses = xyes]) ## Build LaTeX docs if possible... AC_PATH_PROG(PDFLATEX, pdflatex) AM_CONDITIONAL(WITH_PDFLATEX, [test x$PDFLATEX != x]) ## ... unless told otherwise! AC_ARG_ENABLE([pdfmanual], [AC_HELP_STRING(--enable-pdfmanual, [build and install the manual])], [], [enable_pdfmanual=no]) if test x$enable_pdfmanual = xyes; then AC_MSG_WARN([Building Rivet manual, by request]) fi AM_CONDITIONAL(ENABLE_PDFMANUAL, [test x$enable_pdfmanual = xyes]) ## Build Doxygen documentation if possible AC_ARG_ENABLE([doxygen], [AC_HELP_STRING(--disable-doxygen, [don't try to make Doxygen documentation])], [], [enable_doxygen=yes]) if test x$enable_doxygen = xyes; then AC_PATH_PROG(DOXYGEN, doxygen) fi AM_CONDITIONAL(WITH_DOXYGEN, [test x$DOXYGEN != x]) ## Build asciidoc docs if possible AC_PATH_PROG(ASCIIDOC, asciidoc) AM_CONDITIONAL(WITH_ASCIIDOC, [test x$ASCIIDOC != x]) ## Python extension AC_ARG_ENABLE(pyext, [AC_HELP_STRING(--disable-pyext, [don't build Python module (default=build)])], [], [enable_pyext=yes]) ## Basic Python checks if test x$enable_pyext == xyes; then AX_PYTHON_DEVEL([>= '2.7.3']) AC_SUBST(PYTHON_VERSION) RIVET_PYTHONPATH=`$PYTHON -c "from __future__ import print_function; import distutils.sysconfig; print(distutils.sysconfig.get_python_lib(prefix='$prefix', plat_specific=True));"` AC_SUBST(RIVET_PYTHONPATH) if test -z "$PYTHON"; then AC_MSG_ERROR([Can't build Python extension since python can't be found]) enable_pyext=no fi if test -z "$PYTHON_CPPFLAGS"; then AC_MSG_ERROR([Can't build Python extension since Python.h header file cannot be found]) enable_pyext=no fi fi AM_CONDITIONAL(ENABLE_PYEXT, [test x$enable_pyext == xyes]) dnl dnl setup.py puts its build artifacts into a labelled path dnl this helps the test scripts to find them locally instead of dnl having to install first dnl RIVET_SETUP_PY_PATH=$(${PYTHON} -c 'from __future__ import print_function; import distutils.util as u, sys; vi=sys.version_info; print("lib.%s-%s.%s" % (u.get_platform(),vi.major, vi.minor))') AC_SUBST(RIVET_SETUP_PY_PATH) ## Cython checks if test x$enable_pyext == xyes; then AM_CHECK_CYTHON([0.24.0], [:], [:]) if test x$CYTHON_FOUND = xyes; then AC_MSG_NOTICE([Cython >= 0.24 found: Python extension source can be rebuilt (for developers)]) fi AC_CHECK_FILE([pyext/rivet/core.cpp], [], [if test "x$CYTHON_FOUND" != "xyes"; then AC_MSG_ERROR([Cython is required for --enable-pyext, no pre-built core.cpp was found.]) fi]) cython_compiler=$CXX ## Set extra Python extension build flags (to cope with Cython output code oddities) PYEXT_CXXFLAGS="$CXXFLAGS" AC_CEDAR_CHECKCXXFLAG([-Wno-unused-but-set-variable], [PYEXT_CXXFLAGS="$PYEXT_CXXFLAGS -Wno-unused-but-set-variable"]) AC_CEDAR_CHECKCXXFLAG([-Wno-sign-compare], [PYEXT_CXXFLAGS="$PYEXT_CXXFLAGS -Wno-sign-compare"]) AC_SUBST(PYEXT_CXXFLAGS) AC_MSG_NOTICE([All Python build checks successful: 'rivet' Python extension will be built]) fi AM_CONDITIONAL(WITH_CYTHON, [test x$CYTHON_FOUND = xyes]) ## Set default build flags AM_CPPFLAGS="-I\$(top_srcdir)/include -I\$(top_builddir)/include" #AM_CPPFLAGS="$AM_CPPFLAGS -I\$(top_srcdir)/include/eigen3" #AM_CPPFLAGS="$AM_CPPFLAGS \$(GSL_CPPFLAGS)" dnl AM_CPPFLAGS="$AM_CPPFLAGS \$(BOOST_CPPFLAGS)" AM_CPPFLAGS="$AM_CPPFLAGS -I\$(YODAINCPATH)" AM_CPPFLAGS="$AM_CPPFLAGS -I\$(HEPMCINCPATH)" AM_CPPFLAGS="$AM_CPPFLAGS -I\$(FASTJETINCPATH)" AC_CEDAR_CHECKCXXFLAG([-pedantic], [AM_CXXFLAGS="$AM_CXXFLAGS -pedantic"]) AC_CEDAR_CHECKCXXFLAG([-Wall], [AM_CXXFLAGS="$AM_CXXFLAGS -Wall"]) AC_CEDAR_CHECKCXXFLAG([-Wno-long-long], [AM_CXXFLAGS="$AM_CXXFLAGS -Wno-long-long"]) AC_CEDAR_CHECKCXXFLAG([-Wno-format], [AM_CXXFLAGS="$AM_CXXFLAGS -Wno-format"]) dnl AC_CEDAR_CHECKCXXFLAG([-Wno-unused-variable], [AM_CXXFLAGS="$AM_CXXFLAGS -Wno-unused-variable"]) AC_CEDAR_CHECKCXXFLAG([-Werror=uninitialized], [AM_CXXFLAGS="$AM_CXXFLAGS -Werror=uninitialized"]) AC_CEDAR_CHECKCXXFLAG([-Werror=delete-non-virtual-dtor], [AM_CXXFLAGS="$AM_CXXFLAGS -Werror=delete-non-virtual-dtor"]) ## Add OpenMP-enabling flags if possible AX_OPENMP([AM_CXXFLAGS="$AM_CXXFLAGS $OPENMP_CXXFLAGS"]) ## Optional zlib support for gzip-compressed data streams/files AX_CHECK_ZLIB ## Debug flag (default=-DNDEBUG, enabled=-g) AC_ARG_ENABLE([debug], [AC_HELP_STRING(--enable-debug, [build with debugging symbols @<:@default=no@:>@])], [], [enable_debug=no]) if test x$enable_debug == xyes; then AM_CXXFLAGS="$AM_CXXFLAGS -g" fi ## Extra warnings flag (default=none) AC_ARG_ENABLE([extra-warnings], [AC_HELP_STRING(--enable-extra-warnings, [build with extra compiler warnings (recommended for developers) @<:@default=no@:>@])], [], [enable_extra_warnings=no]) if test x$enable_extra_warnings == xyes; then AC_CEDAR_CHECKCXXFLAG([-Wextra], [AM_CXXFLAGS="$AM_CXXFLAGS -Wextra "]) fi AC_SUBST(AM_CPPFLAGS) AC_SUBST(AM_CXXFLAGS) AC_EMPTY_SUBST AC_CONFIG_FILES(Makefile Doxyfile) AC_CONFIG_FILES(include/Makefile include/Rivet/Makefile) AC_CONFIG_FILES(src/Makefile) AC_CONFIG_FILES(src/Core/Makefile src/Core/yamlcpp/Makefile) AC_CONFIG_FILES(src/Tools/Makefile) -AC_CONFIG_FILES(src/Tools/Nsubjettiness/Makefile) +AC_CONFIG_FILES(src/Tools/fjcontrib/Makefile) +AC_CONFIG_FILES(src/Tools/fjcontrib/EnergyCorrelator/Makefile) +AC_CONFIG_FILES(src/Tools/fjcontrib/Nsubjettiness/Makefile) AC_CONFIG_FILES(src/Projections/Makefile) AC_CONFIG_FILES(src/AnalysisTools/Makefile) AC_CONFIG_FILES(analyses/Makefile) AC_CONFIG_FILES(test/Makefile) AC_CONFIG_FILES(pyext/Makefile pyext/rivet/Makefile pyext/setup.py) AC_CONFIG_FILES(data/Makefile data/texmf/Makefile) AC_CONFIG_FILES(doc/Makefile) AC_CONFIG_FILES(doc/rivetversion.sty) AC_CONFIG_FILES(bin/Makefile bin/rivet-config bin/rivet-buildplugin) AC_CONFIG_FILES(rivetenv.sh rivetenv.csh rivet.pc) AC_OUTPUT if test x$enable_pyrivet == xyes; then cat < + Lina Necib + Gavin Salam + Jesse Thaler + +For physics details, see: + + Energy Correlation Functions for Jet Substructure. + Andrew J. Larkoski, Gavin P. Salam, and Jesse Thaler. + JHEP 1306, 108 (2013) + arXiv:1305.0007. + + Power Counting to Better Jet Observables. + Andrew J. Larkoski, Ian Moult, and Duff Neill. + JHEP 1412, 009 (2014) + arXiv:1409.6298. + + New Angles on Energy Correlation Functions. + Ian Moult, Lina Necib, and Jesse Thaler. + arXiv:1609.07483. +---------------------------------------------------------------------- diff --git a/src/Tools/fjcontrib/EnergyCorrelator/COPYING b/src/Tools/fjcontrib/EnergyCorrelator/COPYING new file mode 100644 --- /dev/null +++ b/src/Tools/fjcontrib/EnergyCorrelator/COPYING @@ -0,0 +1,360 @@ +The EnergyCorrelator contrib to FastJet is released +under the terms of the GNU General Public License v2 (GPLv2). + +A copy of the GPLv2 is to be found at the end of this file. + +While the GPL license grants you considerable freedom, please bear in +mind that this code's use falls under guidelines similar to those that +are standard for Monte Carlo event generators +(http://www.montecarlonet.org/GUIDELINES). In particular, if you use +this code as part of work towards a scientific publication, whether +directly or contained within another program you should include a citation +to + + arXiv:1305.0007 + arXiv:1409.6298 + arXiv:1609.07483 + +====================================================================== +====================================================================== +====================================================================== + GNU GENERAL PUBLIC LICENSE + Version 2, June 1991 + + Copyright (C) 1989, 1991 Free Software Foundation, Inc. + 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + Preamble + + The licenses for most software are designed to take away your +freedom to share and change it. By contrast, the GNU General Public +License is intended to guarantee your freedom to share and change free +software--to make sure the software is free for all its users. This +General Public License applies to most of the Free Software +Foundation's software and to any other program whose authors commit to +using it. (Some other Free Software Foundation software is covered by +the GNU Library General Public License instead.) You can apply it to +your programs, too. + + When we speak of free software, we are referring to freedom, not +price. Our General Public Licenses are designed to make sure that you +have the freedom to distribute copies of free software (and charge for +this service if you wish), that you receive source code or can get it +if you want it, that you can change the software or use pieces of it +in new free programs; and that you know you can do these things. + + To protect your rights, we need to make restrictions that forbid +anyone to deny you these rights or to ask you to surrender the rights. +These restrictions translate to certain responsibilities for you if you +distribute copies of the software, or if you modify it. + + For example, if you distribute copies of such a program, whether +gratis or for a fee, you must give the recipients all the rights that +you have. You must make sure that they, too, receive or can get the +source code. And you must show them these terms so they know their +rights. + + We protect your rights with two steps: (1) copyright the software, and +(2) offer you this license which gives you legal permission to copy, +distribute and/or modify the software. + + Also, for each author's protection and ours, we want to make certain +that everyone understands that there is no warranty for this free +software. If the software is modified by someone else and passed on, we +want its recipients to know that what they have is not the original, so +that any problems introduced by others will not reflect on the original +authors' reputations. + + Finally, any free program is threatened constantly by software +patents. We wish to avoid the danger that redistributors of a free +program will individually obtain patent licenses, in effect making the +program proprietary. To prevent this, we have made it clear that any +patent must be licensed for everyone's free use or not licensed at all. + + The precise terms and conditions for copying, distribution and +modification follow. + + GNU GENERAL PUBLIC LICENSE + TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION + + 0. This License applies to any program or other work which contains +a notice placed by the copyright holder saying it may be distributed +under the terms of this General Public License. The "Program", below, +refers to any such program or work, and a "work based on the Program" +means either the Program or any derivative work under copyright law: +that is to say, a work containing the Program or a portion of it, +either verbatim or with modifications and/or translated into another +language. (Hereinafter, translation is included without limitation in +the term "modification".) Each licensee is addressed as "you". + +Activities other than copying, distribution and modification are not +covered by this License; they are outside its scope. The act of +running the Program is not restricted, and the output from the Program +is covered only if its contents constitute a work based on the +Program (independent of having been made by running the Program). +Whether that is true depends on what the Program does. + + 1. You may copy and distribute verbatim copies of the Program's +source code as you receive it, in any medium, provided that you +conspicuously and appropriately publish on each copy an appropriate +copyright notice and disclaimer of warranty; keep intact all the +notices that refer to this License and to the absence of any warranty; +and give any other recipients of the Program a copy of this License +along with the Program. + +You may charge a fee for the physical act of transferring a copy, and +you may at your option offer warranty protection in exchange for a fee. + + 2. You may modify your copy or copies of the Program or any portion +of it, thus forming a work based on the Program, and copy and +distribute such modifications or work under the terms of Section 1 +above, provided that you also meet all of these conditions: + + a) You must cause the modified files to carry prominent notices + stating that you changed the files and the date of any change. + + b) You must cause any work that you distribute or publish, that in + whole or in part contains or is derived from the Program or any + part thereof, to be licensed as a whole at no charge to all third + parties under the terms of this License. + + c) If the modified program normally reads commands interactively + when run, you must cause it, when started running for such + interactive use in the most ordinary way, to print or display an + announcement including an appropriate copyright notice and a + notice that there is no warranty (or else, saying that you provide + a warranty) and that users may redistribute the program under + these conditions, and telling the user how to view a copy of this + License. (Exception: if the Program itself is interactive but + does not normally print such an announcement, your work based on + the Program is not required to print an announcement.) + +These requirements apply to the modified work as a whole. If +identifiable sections of that work are not derived from the Program, +and can be reasonably considered independent and separate works in +themselves, then this License, and its terms, do not apply to those +sections when you distribute them as separate works. But when you +distribute the same sections as part of a whole which is a work based +on the Program, the distribution of the whole must be on the terms of +this License, whose permissions for other licensees extend to the +entire whole, and thus to each and every part regardless of who wrote it. + +Thus, it is not the intent of this section to claim rights or contest +your rights to work written entirely by you; rather, the intent is to +exercise the right to control the distribution of derivative or +collective works based on the Program. + +In addition, mere aggregation of another work not based on the Program +with the Program (or with a work based on the Program) on a volume of +a storage or distribution medium does not bring the other work under +the scope of this License. + + 3. You may copy and distribute the Program (or a work based on it, +under Section 2) in object code or executable form under the terms of +Sections 1 and 2 above provided that you also do one of the following: + + a) Accompany it with the complete corresponding machine-readable + source code, which must be distributed under the terms of Sections + 1 and 2 above on a medium customarily used for software interchange; or, + + b) Accompany it with a written offer, valid for at least three + years, to give any third party, for a charge no more than your + cost of physically performing source distribution, a complete + machine-readable copy of the corresponding source code, to be + distributed under the terms of Sections 1 and 2 above on a medium + customarily used for software interchange; or, + + c) Accompany it with the information you received as to the offer + to distribute corresponding source code. (This alternative is + allowed only for noncommercial distribution and only if you + received the program in object code or executable form with such + an offer, in accord with Subsection b above.) + +The source code for a work means the preferred form of the work for +making modifications to it. For an executable work, complete source +code means all the source code for all modules it contains, plus any +associated interface definition files, plus the scripts used to +control compilation and installation of the executable. However, as a +special exception, the source code distributed need not include +anything that is normally distributed (in either source or binary +form) with the major components (compiler, kernel, and so on) of the +operating system on which the executable runs, unless that component +itself accompanies the executable. + +If distribution of executable or object code is made by offering +access to copy from a designated place, then offering equivalent +access to copy the source code from the same place counts as +distribution of the source code, even though third parties are not +compelled to copy the source along with the object code. + + 4. You may not copy, modify, sublicense, or distribute the Program +except as expressly provided under this License. Any attempt +otherwise to copy, modify, sublicense or distribute the Program is +void, and will automatically terminate your rights under this License. +However, parties who have received copies, or rights, from you under +this License will not have their licenses terminated so long as such +parties remain in full compliance. + + 5. You are not required to accept this License, since you have not +signed it. However, nothing else grants you permission to modify or +distribute the Program or its derivative works. These actions are +prohibited by law if you do not accept this License. Therefore, by +modifying or distributing the Program (or any work based on the +Program), you indicate your acceptance of this License to do so, and +all its terms and conditions for copying, distributing or modifying +the Program or works based on it. + + 6. Each time you redistribute the Program (or any work based on the +Program), the recipient automatically receives a license from the +original licensor to copy, distribute or modify the Program subject to +these terms and conditions. You may not impose any further +restrictions on the recipients' exercise of the rights granted herein. +You are not responsible for enforcing compliance by third parties to +this License. + + 7. If, as a consequence of a court judgment or allegation of patent +infringement or for any other reason (not limited to patent issues), +conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot +distribute so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you +may not distribute the Program at all. For example, if a patent +license would not permit royalty-free redistribution of the Program by +all those who receive copies directly or indirectly through you, then +the only way you could satisfy both it and this License would be to +refrain entirely from distribution of the Program. + +If any portion of this section is held invalid or unenforceable under +any particular circumstance, the balance of the section is intended to +apply and the section as a whole is intended to apply in other +circumstances. + +It is not the purpose of this section to induce you to infringe any +patents or other property right claims or to contest validity of any +such claims; this section has the sole purpose of protecting the +integrity of the free software distribution system, which is +implemented by public license practices. Many people have made +generous contributions to the wide range of software distributed +through that system in reliance on consistent application of that +system; it is up to the author/donor to decide if he or she is willing +to distribute software through any other system and a licensee cannot +impose that choice. + +This section is intended to make thoroughly clear what is believed to +be a consequence of the rest of this License. + + 8. If the distribution and/or use of the Program is restricted in +certain countries either by patents or by copyrighted interfaces, the +original copyright holder who places the Program under this License +may add an explicit geographical distribution limitation excluding +those countries, so that distribution is permitted only in or among +countries not thus excluded. In such case, this License incorporates +the limitation as if written in the body of this License. + + 9. The Free Software Foundation may publish revised and/or new versions +of the General Public License from time to time. Such new versions will +be similar in spirit to the present version, but may differ in detail to +address new problems or concerns. + +Each version is given a distinguishing version number. If the Program +specifies a version number of this License which applies to it and "any +later version", you have the option of following the terms and conditions +either of that version or of any later version published by the Free +Software Foundation. If the Program does not specify a version number of +this License, you may choose any version ever published by the Free Software +Foundation. + + 10. If you wish to incorporate parts of the Program into other free +programs whose distribution conditions are different, write to the author +to ask for permission. For software which is copyrighted by the Free +Software Foundation, write to the Free Software Foundation; we sometimes +make exceptions for this. Our decision will be guided by the two goals +of preserving the free status of all derivatives of our free software and +of promoting the sharing and reuse of software generally. + + NO WARRANTY + + 11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY +FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN +OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES +PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED +OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF +MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS +TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE +PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING, +REPAIR OR CORRECTION. + + 12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING +WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR +REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, +INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING +OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED +TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY +YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER +PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE +POSSIBILITY OF SUCH DAMAGES. + + END OF TERMS AND CONDITIONS + + How to Apply These Terms to Your New Programs + + If you develop a new program, and you want it to be of the greatest +possible use to the public, the best way to achieve this is to make it +free software which everyone can redistribute and change under these terms. + + To do so, attach the following notices to the program. It is safest +to attach them to the start of each source file to most effectively +convey the exclusion of warranty; and each file should have at least +the "copyright" line and a pointer to where the full notice is found. + + + Copyright (C) + + This program 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. + + This program 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 program; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + +Also add information on how to contact you by electronic and paper mail. + +If the program is interactive, make it output a short notice like this +when it starts in an interactive mode: + + Gnomovision version 69, Copyright (C) year name of author + Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'. + This is free software, and you are welcome to redistribute it + under certain conditions; type `show c' for details. + +The hypothetical commands `show w' and `show c' should show the appropriate +parts of the General Public License. Of course, the commands you use may +be called something other than `show w' and `show c'; they could even be +mouse-clicks or menu items--whatever suits your program. + +You should also get your employer (if you work as a programmer) or your +school, if any, to sign a "copyright disclaimer" for the program, if +necessary. Here is a sample; alter the names: + + Yoyodyne, Inc., hereby disclaims all copyright interest in the program + `Gnomovision' (which makes passes at compilers) written by James Hacker. + + , 1 April 1989 + Ty Coon, President of Vice + +This General Public License does not permit incorporating your program into +proprietary programs. If your program is a subroutine library, you may +consider it more useful to permit linking proprietary applications with the +library. If this is what you want to do, use the GNU Library General +Public License instead of this License. diff --git a/src/Tools/fjcontrib/EnergyCorrelator/ChangeLog b/src/Tools/fjcontrib/EnergyCorrelator/ChangeLog new file mode 100644 --- /dev/null +++ b/src/Tools/fjcontrib/EnergyCorrelator/ChangeLog @@ -0,0 +1,307 @@ +2018-02-08 Jesse Thaler + + * NEWS + * VERSION + Making 1.3.1 + + +2018-02-08 Lina Necib + + * EnergyCorrelator.cc + Debugging memory allocation, deleting energyStore/angleStore properly this time + + +2018-02-04 Lina Necib + + * EnergyCorrelator.cc + Debugging memory allocation, deleting energyStore/angleStore + + * VERSION + Making 1.3.1-devel + +2018-01-09 Jesse Thaler + + * NEWS + * VERSION + Finalizing 1.3.0 + +2018-01-09 Lina Necib + + * NEWS + Preparing for version candidate + + * VERSION + Making 1.3.0-rc1 + +2018-01-07 Lina Necib + + * EnergyCorrelator.hh/cc + Sped up evaluations of ECFs. + Functions are defined in ECFs and called to evaluate ECFGs. + + * VERSION + Making 1.2.3-devel + + +2018-01-04 Lina Necib + + * EnergyCorrelator.hh/cc + Sped up evaluations of ECFG by a factor of 4. + Refactored some of the code in functions + + * VERSION + Making 1.2.2-devel + + * example.cc + Added some timings tests for ECFGs, N2, and N3. + + +2017-01-25 Jesse Thaler + + * EnergyCorrelator.hh/cc + Converting all _N to unsigned integers, removing _N < 0 warning + Added warning to constructors for N < 1 in some cases. + + * VERSION + Making 1.2.1-devel + + * Makefile + Added -Wextra to compiler flags + +2016-10-07 Jesse Thaler + + * AUTHORS/COPYING: + Updated the publication arXiv number to 1609.07483. + + * EnergyCorrelator.hh + Fixed typo in EnergyCorrelatorGeneralized description + + * NEWS + Added "USeries" to news. + + * VERSION + Changed version to 1.2.0 + +2016-09-27 Lina Necib + + * EnergyCorrelator.hh/README + Updated the publication arXiv number to 1609.07483. + +2016-09-23 Lina Necib + + * EnergyCorrelator.cc + Made the evaluation of ECFG faster by replacing sort -> partial_sort and simplified multiplication syntax + * example.cc/ref + Fixed minor typo + * VERSION + Changed version to 1.2.0-rc3 + +2016-09-23 Lina Necib + + * EnergyCorrelator.hh + Fixed minor doxygen typo + * example.cc/ref + Changed EnergyCorrelatorNormalized -> EnergyCorrelatorGeneralized in function explanations + * VERSION + Changed version to 1.2.0-rc2 + +2016-09-22 Jesse Thaler + + * EnergyCorrelator.cc/hh + Renamed EnergyCorrelatorNormalized -> EnergyCorrelatorGeneralized + Changed order of arguments for EnergyCorrelatorGeneralized + Updated doxygen documentation + * example_basic_usage.cc + Changed EnergyCorrelatorNormalized -> EnergyCorrelatorGeneralized + * README + Updated for EnergyCorrelatorGeneralized + +2016-09-15 Lina Necib + + *VERSION: + 1.2.0-rc1 + +2016-08-24 Jesse Thaler + + Minor comment fixes throughout. + + * EnergyCorrelator.cc/hh + Put in _helper_correlator when defining overall normalization + Removed angle error for result_all_angles() + *NEWS: + Made this more readable. + *README: + Clarified some of the documentation. + +2016-08-23 Lina Necib + + *VERSION: + *NEWS: + *AUTHORS: + *COPYING: + *README: + *EnergyCorrelator.cc + Added if statement that the ECF and ECFNs return exactly zero + if the number of particles is less than N of the ECF. + *EnergyCorrelator.hh + *example.cc + *example_basic_usage.cc + Updated this example. + +2016-08-21 Lina Necib + + *VERSION: + *NEWS: + *AUTHORS: + *COPYING: + *README: + *EnergyCorrelator.cc + Added Cseries. + *EnergyCorrelator.hh + Added Cseries. + *example.cc + Added example of Cseries + *example_basic_usage.cc + Simplified examples. + +2016-08-17 Lina Necib + + *VERSION: + *NEWS: + *AUTHORS: + Added placeholder for new reference + *COPYING: + Added placeholder for new reference + *README: + added information on different measures E_inv + *EnergyCorrelator.cc + Minor text edits + added comments + *EnergyCorrelator.hh + Minor text edits + added comments + *example_basic_usage.cc + added a simplified example file that shows the use of the + different observables. + + +2016-06-23 Lina Necib + + *VERSION: + 1.2.0-alpha0 + *NEWS: + *AUTHORS: + Lina Necib + *COPYING: + *README: + added descriptions of functions that will appear shortly + arXiv:160X.XXXXX + + *EnergyCorrelator.cc + added code to calculate normalized ECFS. + *EnergyCorrelator.hh + added code to calculate normalized ECFS, Nseries, generalized + D2, N2, N3, and M2. + + *example.cc + added calculation of normalized ECFS, Nseries, generalized + D2, N2, N3, and M2 to example file. + + +2014-11-20 Jesse Thaler + + * README: + Typos fixed + +2014-11-13 Andrew Larkoski + + *VERSION: + *NEWS: + release of version 1.1.0 + + *AUTHORS: + *COPYING: + *README: + added reference to arXiv:1409.6298. + + *EnergyCorrelator.cc + *EnergyCorrelator.hh + added code to calculate C1, C2, and D2 observables. + + *example.cc + added calculation of C1, C2, and D2 to example file. + +2013-05-01 Gavin Salam + + * VERSION: + * NEWS: + release of version 1.0.1 + + * README: + updated to reflect v1.0 interface. + +2013-05-01 Jesse Thaler + + * EnergyCorrelator.cc + Switched from NAN to std::numeric_limits::quiet_NaN() + +2013-04-30 Jesse Thaler + + * EnergyCorrelator.cc + Implemented Gregory Soyez's suggestions on errors/asserts + +2013-04-30 Gavin Salam + + * VERSION: + * NEWS: + released v. 1.0.0 + + * EnergyCorrelator.hh: + * example.cc + small changes to documentation to reflect the change below + an + gave an explicit command-line in the example. + +2013-04-29 Jesse Thaler + + * EnergyCorrelator.cc + - Added support for N = 5 + + * example.cc|.ref + - Added N = 5 to test suite. + +2013-04-29 Gavin Salam + + * EnergyCorrelator.hh|cc: + - reduced memory usage from roughly 8n^2 to 4n^2 bytes (note that + sums are now performed in a different order, so results may + change to within rounding error). + + - Implemented description() for all three classes. + + - Worked on doxygen-style comments and moved some bits of code into + the .cc file. + + * Doxyfile: *** ADDED *** + + * example.cc: + developers' timing component now uses clock(), to get + finer-grained timing, and also outputs a description for some of + the correlators being used. + +2013-04-26 + 04-27 Jesse Thaler + + * EnergyCorrelator.hh|cc: + added temporary storage array for interparticle angles, to speed + up the energy correlator calculation for N>2 + + * example.cc + has internal option for printing out timing info. + +2013-04-26 Gavin Salam + Jesse & Andrew + + * Makefile: + added explicit dependencies + + * example.cc (analyze): + added a little bit of commented code for timing tests. + +2013-03-10 + Initial creation diff --git a/src/Tools/fjcontrib/EnergyCorrelator/EnergyCorrelator.cc b/src/Tools/fjcontrib/EnergyCorrelator/EnergyCorrelator.cc new file mode 100644 --- /dev/null +++ b/src/Tools/fjcontrib/EnergyCorrelator/EnergyCorrelator.cc @@ -0,0 +1,1091 @@ +// EnergyCorrelator Package +// Questions/Comments? Email the authors: +// larkoski@mit.edu, lnecib@mit.edu, +// gavin.salam@cern.ch jthaler@jthaler.net +// +// Copyright (c) 2013-2016 +// Andrew Larkoski, Lina Necib, Gavin Salam, and Jesse Thaler +// +// $Id: EnergyCorrelator.cc 1106 2018-02-09 01:47:28Z linoush $ +//---------------------------------------------------------------------- +// 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 . +//---------------------------------------------------------------------- + +#include "EnergyCorrelator.hh" +#include +#include +using namespace std; + +FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh + +namespace contrib { + + double EnergyCorrelator::result(const PseudoJet& jet) const { + + // if jet does not have constituents, throw error + if (!jet.has_constituents()) throw Error("EnergyCorrelator called on jet with no constituents."); + + // get N = 0 case out of the way + if (_N == 0) return 1.0; + + // find constituents + std::vector particles = jet.constituents(); + + // return zero if the number of constituents is less than _N + if (particles.size() < _N) return 0.0 ; + + double answer = 0.0; + + // take care of N = 1 case. + if (_N == 1) { + for (unsigned int i = 0; i < particles.size(); i++) { + answer += energy(particles[i]); + } + return answer; + } + + double half_beta = _beta/2.0; + + // take care of N = 2 case. + if (_N == 2) { + for (unsigned int i = 0; i < particles.size(); i++) { + for (unsigned int j = i + 1; j < particles.size(); j++) { //note offset by one so that angle is never called on identical pairs + answer += energy(particles[i]) + * energy(particles[j]) + * pow(angleSquared(particles[i],particles[j]), half_beta); + } + } + return answer; + } + + + // if N > 5, then throw error + if (_N > 5) { + throw Error("EnergyCorrelator is only hard coded for N = 0,1,2,3,4,5"); + } + + + // Now deal with N = 3,4,5. Different options if storage array is used or not. + if (_strategy == storage_array) { + + // For N > 2, fill static storage array to save computation time. + unsigned int nC = particles.size(); + // Make energy storage + double *energyStore = new double[nC]; + + // Make angular storage + double **angleStore = new double*[nC]; + + precompute_energies_and_angles(particles, energyStore, angleStore); + + // Define n_angles so it is the same function for ECFs and ECFGs + unsigned int n_angles = _N * (_N - 1) / 2; + // now do recursion + if (_N == 3) { + answer = evaluate_n3(nC, n_angles, energyStore, angleStore); + } else if (_N == 4) { + answer = evaluate_n4(nC, n_angles, energyStore, angleStore); + } else if (_N == 5) { + answer = evaluate_n5(nC, n_angles, energyStore, angleStore); + } else { + assert(_N <= 5); + } + // Deleting arrays + delete[] energyStore; + + for (unsigned int i = 0; i < particles.size(); i++) { + delete[] angleStore[i]; + } + delete[] angleStore; + + } else if (_strategy == slow) { + if (_N == 3) { + for (unsigned int i = 0; i < particles.size(); i++) { + for (unsigned int j = i + 1; j < particles.size(); j++) { + double ans_ij = energy(particles[i]) + * energy(particles[j]) + * pow(angleSquared(particles[i],particles[j]), half_beta); + for (unsigned int k = j + 1; k < particles.size(); k++) { + answer += ans_ij + * energy(particles[k]) + * pow(angleSquared(particles[i],particles[k]), half_beta) + * pow(angleSquared(particles[j],particles[k]), half_beta); + } + } + } + } else if (_N == 4) { + for (unsigned int i = 0; i < particles.size(); i++) { + for (unsigned int j = i + 1; j < particles.size(); j++) { + double ans_ij = energy(particles[i]) + * energy(particles[j]) + * pow(angleSquared(particles[i],particles[j]), half_beta); + for (unsigned int k = j + 1; k < particles.size(); k++) { + double ans_ijk = ans_ij + * energy(particles[k]) + * pow(angleSquared(particles[i],particles[k]), half_beta) + * pow(angleSquared(particles[j],particles[k]), half_beta); + for (unsigned int l = k + 1; l < particles.size(); l++) { + answer += ans_ijk + * energy(particles[l]) + * pow(angleSquared(particles[i],particles[l]), half_beta) + * pow(angleSquared(particles[j],particles[l]), half_beta) + * pow(angleSquared(particles[k],particles[l]), half_beta); + } + } + } + } + } else if (_N == 5) { + for (unsigned int i = 0; i < particles.size(); i++) { + for (unsigned int j = i + 1; j < particles.size(); j++) { + double ans_ij = energy(particles[i]) + * energy(particles[j]) + * pow(angleSquared(particles[i],particles[j]), half_beta); + for (unsigned int k = j + 1; k < particles.size(); k++) { + double ans_ijk = ans_ij + * energy(particles[k]) + * pow(angleSquared(particles[i],particles[k]), half_beta) + * pow(angleSquared(particles[j],particles[k]), half_beta); + for (unsigned int l = k + 1; l < particles.size(); l++) { + double ans_ijkl = ans_ijk + * energy(particles[l]) + * pow(angleSquared(particles[i],particles[l]), half_beta) + * pow(angleSquared(particles[j],particles[l]), half_beta) + * pow(angleSquared(particles[k],particles[l]), half_beta); + for (unsigned int m = l + 1; m < particles.size(); m++) { + answer += ans_ijkl + * energy(particles[m]) + * pow(angleSquared(particles[i],particles[m]), half_beta) + * pow(angleSquared(particles[j],particles[m]), half_beta) + * pow(angleSquared(particles[k],particles[m]), half_beta) + * pow(angleSquared(particles[l],particles[m]), half_beta); + } + } + } + } + } + } else { + assert(_N <= 5); + } + } else { + assert(_strategy == slow || _strategy == storage_array); + } + + return answer; + } + + double EnergyCorrelator::energy(const PseudoJet& jet) const { + if (_measure == pt_R) { + return jet.perp(); + } else if (_measure == E_theta || _measure == E_inv) { + return jet.e(); + } else { + assert(_measure==pt_R || _measure==E_theta || _measure==E_inv); + return std::numeric_limits::quiet_NaN(); + } + } + + double EnergyCorrelator::angleSquared(const PseudoJet& jet1, const PseudoJet& jet2) const { + if (_measure == pt_R) { + return jet1.squared_distance(jet2); + } else if (_measure == E_theta) { + // doesn't seem to be a fastjet built in for this + double dot = jet1.px()*jet2.px() + jet1.py()*jet2.py() + jet1.pz()*jet2.pz(); + double norm1 = jet1.px()*jet1.px() + jet1.py()*jet1.py() + jet1.pz()*jet1.pz(); + double norm2 = jet2.px()*jet2.px() + jet2.py()*jet2.py() + jet2.pz()*jet2.pz(); + + double costheta = dot/(sqrt(norm1 * norm2)); + if (costheta > 1.0) costheta = 1.0; // Need to handle case of numerical overflow + double theta = acos(costheta); + return theta*theta; + + } else if (_measure == E_inv) { + if (jet1.E() < 0.0000001 || jet2.E() < 0.0000001) return 0.0; + else { + double dot4 = max(jet1.E()*jet2.E() - jet1.px()*jet2.px() - jet1.py()*jet2.py() - jet1.pz()*jet2.pz(),0.0); + return 2.0 * dot4 / jet1.E() / jet2.E(); + } + } else { + assert(_measure==pt_R || _measure==E_theta || _measure==E_inv); + return std::numeric_limits::quiet_NaN(); + } + } + + + double EnergyCorrelator::multiply_angles(double angle_list[], int n_angles, unsigned int N_total) const { + // Compute the product of the n_angles smallest angles. + // std::partial_sort could also work, but since angle_list contains + // less than 10 elements, this way is usually faster. + double product = 1; + + for (int a = 0; a < n_angles; a++) { + double cur_min = angle_list[0]; + int cur_min_pos = 0; + for (unsigned int b = 1; b < N_total; b++) { + if (angle_list[b] < cur_min) { + cur_min = angle_list[b]; + cur_min_pos = b; + } + } + + // multiply it by the next smallest + product *= cur_min; + angle_list[cur_min_pos] = INT_MAX; + } + return product; + } + + void EnergyCorrelator::precompute_energies_and_angles(std::vector const &particles, double* energyStore, double** angleStore) const { + // Fill storage with energy/angle information + unsigned int nC = particles.size(); + for (unsigned int i = 0; i < nC; i++) { + angleStore[i] = new double[i]; + } + + double half_beta = _beta/2.0; + for (unsigned int i = 0; i < particles.size(); i++) { + energyStore[i] = energy(particles[i]); + for (unsigned int j = 0; j < i; j++) { + if (half_beta == 1.0){ + angleStore[i][j] = angleSquared(particles[i], particles[j]); + } else { + angleStore[i][j] = pow(angleSquared(particles[i], particles[j]), half_beta); + } + } + } + } + + double EnergyCorrelator::evaluate_n3(unsigned int nC, unsigned int n_angles, double* energyStore, double** angleStore) const { + unsigned int N_total = 3; + double angle1, angle2, angle3; + double angle; + double answer = 0; + + for (unsigned int i = 2; i < nC; i++) { + for (unsigned int j = 1; j < i; j++) { + double mult_energy_i_j = energyStore[i] * energyStore[j]; + + for (unsigned int k = 0; k < j; k++) { + angle1 = angleStore[i][j]; + angle2 = angleStore[i][k]; + angle3 = angleStore[j][k]; + + double angle_list[] = {angle1, angle2, angle3}; + + if (n_angles == N_total) { + angle = angle1 * angle2 * angle3; + } else { + angle = multiply_angles(angle_list, n_angles, N_total); + } + + answer += mult_energy_i_j + * energyStore[k] + * angle; + } + } + } + return answer; + } + + double EnergyCorrelator::evaluate_n4(unsigned int nC, unsigned int n_angles, double* energyStore, double** angleStore) const { + double answer = 0; + double angle1, angle2, angle3, angle4, angle5, angle6; + unsigned int N_total = 6; + double angle; + + for (unsigned int i = 3; i < nC; i++) { + for (unsigned int j = 2; j < i; j++) { + for (unsigned int k = 1; k < j; k++) { + for (unsigned int l = 0; l < k; l++) { + + angle1 = angleStore[i][j]; + angle2 = angleStore[i][k]; + angle3 = angleStore[i][l]; + angle4 = angleStore[j][k]; + angle5 = angleStore[j][l]; + angle6 = angleStore[k][l]; + + double angle_list[] = {angle1, angle2, angle3, angle4, angle5, angle6}; + + if (n_angles == N_total) { + angle = angle1 * angle2 * angle3 * angle4 * angle5 * angle6; + } else { + angle = multiply_angles(angle_list, n_angles, N_total); + } + + answer += energyStore[i] + * energyStore[j] + * energyStore[k] + * energyStore[l] + * angle; + } + } + } + } + return answer; + } + + double EnergyCorrelator::evaluate_n5(unsigned int nC, unsigned int n_angles, double* energyStore, double** angleStore) const { + + double answer = 0; + double angle1, angle2, angle3, angle4, angle5, angle6, angle7, angle8, angle9, angle10; + unsigned int N_total = 10; + double angle; + + for (unsigned int i = 4; i < nC; i++) { + for (unsigned int j = 3; j < i; j++) { + for (unsigned int k = 2; k < j; k++) { + for (unsigned int l = 1; l < k; l++) { + for (unsigned int m = 0; m < l; m++) { + + angle1 = angleStore[i][j]; + angle2 = angleStore[i][k]; + angle3 = angleStore[i][l]; + angle4 = angleStore[i][m]; + angle5 = angleStore[j][k]; + angle6 = angleStore[j][l]; + angle7 = angleStore[j][m]; + angle8 = angleStore[k][l]; + angle9 = angleStore[k][m]; + angle10 = angleStore[l][m]; + + double angle_list[] = {angle1, angle2, angle3, angle4, angle5, angle6, angle7, angle8, + angle9, angle10}; + + angle = multiply_angles(angle_list, n_angles, N_total); + + answer += energyStore[i] + * energyStore[j] + * energyStore[k] + * energyStore[l] + * energyStore[m] + * angle; + } + } + } + } + } + return answer; + } + + + + double EnergyCorrelatorGeneralized::multiply_angles(double angle_list[], int n_angles, unsigned int N_total) const { + + return _helper_correlator.multiply_angles(angle_list, n_angles, N_total); + } + + void EnergyCorrelatorGeneralized::precompute_energies_and_angles(std::vector const &particles, double* energyStore, double** angleStore) const { + + return _helper_correlator.precompute_energies_and_angles(particles, energyStore, angleStore); + } + + double EnergyCorrelatorGeneralized::evaluate_n3(unsigned int nC, unsigned int n_angles, double* energyStore, double** angleStore) const { + + return _helper_correlator.evaluate_n3(nC, n_angles, energyStore, angleStore); + } + + double EnergyCorrelatorGeneralized::evaluate_n4(unsigned int nC, unsigned int n_angles, double* energyStore, double** angleStore) const { + + return _helper_correlator.evaluate_n4(nC, n_angles, energyStore, angleStore); + } + + double EnergyCorrelatorGeneralized::evaluate_n5(unsigned int nC, unsigned int n_angles, double* energyStore, double** angleStore) const { + + return _helper_correlator.evaluate_n5(nC, n_angles, energyStore, angleStore); + } + + double EnergyCorrelatorGeneralized::result(const PseudoJet& jet) const { + + // if jet does not have constituents, throw error + if (!jet.has_constituents()) throw Error("EnergyCorrelator called on jet with no constituents."); + + // Throw an error if N < 0 + // Not needed if N is unsigned integer + //if (_N < 0 ) throw Error("N cannot be negative"); + // get N = 0 case out of the way + if (_N == 0) return 1.0; + + // take care of N = 1 case. + if (_N == 1) return 1.0; + + // find constituents + std::vector particles = jet.constituents(); + double answer = 0.0; + + // return zero if the number of constituents is less than _N for the ECFG + if (particles.size() < _N) return 0.0 ; + + // The normalization is the energy or pt of the jet, which is also ECF(1, beta) + double EJ = _helper_correlator.result(jet); + + // The overall normalization + double norm = pow(EJ, _N); + + // Find the max number of angles and throw an error if unsuitable + int N_total = int(_N*(_N-1)/2); + if (_angles > N_total) throw Error("Requested number of angles for EnergyCorrelatorGeneralized is larger than number of angles available"); + if (_angles < -1) throw Error("Negative number of angles called for EnergyCorrelatorGeneralized"); + + double half_beta = _beta/2.0; + + // take care of N = 2 case. + if (_N == 2) { + for (unsigned int i = 0; i < particles.size(); i++) { + for (unsigned int j = i + 1; j < particles.size(); j++) { //note offset by one so that angle is never called on identical pairs + answer += energy(particles[i]) + * energy(particles[j]) + * pow(angleSquared(particles[i],particles[j]), half_beta)/norm; + } + } + return answer; + } + + + // if N > 4, then throw error + if (_N > 5) { + throw Error("EnergyCorrelatorGeneralized is only hard coded for N = 0,1,2,3,4,5"); + } + + // Now deal with N = 3,4,5. Different options if storage array is used or not. + if (_strategy == EnergyCorrelator::storage_array) { + + // For N > 2, fill static storage array to save computation time. + + unsigned int nC = particles.size(); + // Make energy storage +// double energyStore[nC]; + double *energyStore = new double[nC]; + + // Make angular storage +// double angleStore[nC][nC]; + double **angleStore = new double*[nC]; + + precompute_energies_and_angles(particles, energyStore, angleStore); + + unsigned int n_angles = _angles; + if (_angles < 0) { + n_angles = N_total; + } + + // now do recursion + if (_N == 3) { + answer = evaluate_n3(nC, n_angles, energyStore, angleStore) / norm; + } else if (_N == 4) { + answer = evaluate_n4(nC, n_angles, energyStore, angleStore) / norm; + } else if (_N == 5) { + answer = evaluate_n5(nC, n_angles, energyStore, angleStore) / norm; + } else { + assert(_N <= 5); + } + // Deleting arrays + delete[] energyStore; + + for (unsigned int i = 0; i < particles.size(); i++) { + delete[] angleStore[i]; + } + delete[] angleStore; + } else if (_strategy == EnergyCorrelator::slow) { + if (_N == 3) { + unsigned int N_total = 3; + double angle1, angle2, angle3; + double angle; + + for (unsigned int i = 0; i < particles.size(); i++) { + for (unsigned int j = i + 1; j < particles.size(); j++) { + for (unsigned int k = j + 1; k < particles.size(); k++) { + + angle1 = angleSquared(particles[i], particles[j]); + angle2 = angleSquared(particles[i], particles[k]); + angle3 = angleSquared(particles[j], particles[k]); + + if (_angles == -1){ + angle = angle1*angle2*angle3; + } else { + double angle_list[] = {angle1, angle2, angle3}; + std::vector angle_vector(angle_list, angle_list + N_total); + std::sort(angle_vector.begin(), angle_vector.begin() + N_total); + + angle = angle_vector[0]; + for ( int l = 1; l < _angles; l++) { angle = angle * angle_vector[l]; } + } + answer += energy(particles[i]) + * energy(particles[j]) + * energy(particles[k]) + * pow(angle, half_beta) /norm; + } + } + } + } else if (_N == 4) { + double angle1, angle2, angle3, angle4, angle5, angle6; + unsigned int N_total = 6; + double angle; + + for (unsigned int i = 0; i < particles.size(); i++) { + for (unsigned int j = i + 1; j < particles.size(); j++) { + for (unsigned int k = j + 1; k < particles.size(); k++) { + for (unsigned int l = k + 1; l < particles.size(); l++) { + + angle1 = angleSquared(particles[i], particles[j]); + angle2 = angleSquared(particles[i], particles[k]); + angle3 = angleSquared(particles[i], particles[l]); + angle4 = angleSquared(particles[j], particles[k]); + angle5 = angleSquared(particles[j], particles[l]); + angle6 = angleSquared(particles[k], particles[l]); + + if(_angles == -1) { + angle = angle1*angle2*angle3*angle4*angle5*angle6; + } else { + + double angle_list[] = {angle1, angle2, angle3, angle4, angle5, angle6}; + std::vector angle_vector(angle_list, angle_list + N_total); + std::sort(angle_vector.begin(), angle_vector.begin() + N_total); + + angle = angle_vector[0]; + for ( int s = 1; s < _angles; s++) { angle = angle * angle_vector[s]; } + + } + answer += energy(particles[i]) + * energy(particles[j]) + * energy(particles[k]) + * energy(particles[l]) + * pow(angle, half_beta)/norm; + } + } + } + } + } else if (_N == 5) { + double angle1, angle2, angle3, angle4, angle5, angle6, angle7, angle8, angle9, angle10; + unsigned int N_total = 10; + double angle; + + for (unsigned int i = 0; i < particles.size(); i++) { + for (unsigned int j = i + 1; j < particles.size(); j++) { + for (unsigned int k = j + 1; k < particles.size(); k++) { + for (unsigned int l = k + 1; l < particles.size(); l++) { + for (unsigned int m = l + 1; m < particles.size(); m++) { + + angle1 = angleSquared(particles[i], particles[j]); + angle2 = angleSquared(particles[i], particles[k]); + angle3 = angleSquared(particles[i], particles[l]); + angle4 = angleSquared(particles[j], particles[k]); + angle5 = angleSquared(particles[j], particles[l]); + angle6 = angleSquared(particles[k], particles[l]); + angle7 = angleSquared(particles[m], particles[i]); + angle8 = angleSquared(particles[m], particles[j]); + angle9 = angleSquared(particles[m], particles[k]); + angle10 = angleSquared(particles[m], particles[l]); + + if (_angles == -1){ + angle = angle1*angle2*angle3*angle4*angle5*angle6*angle7*angle8*angle9*angle10; + } else { + double angle_list[] = {angle1, angle2, angle3, angle4, angle5, angle6, + angle7, angle8, angle9, angle10}; + std::vector angle_vector(angle_list, angle_list + N_total); + std::sort(angle_vector.begin(), angle_vector.begin() + N_total); + + angle = angle_vector[0]; + for ( int s = 1; s < _angles; s++) { angle = angle * angle_vector[s]; } + } + answer += energy(particles[i]) + * energy(particles[j]) + * energy(particles[k]) + * energy(particles[l]) + * energy(particles[m]) + * pow(angle, half_beta) /norm; + } + } + } + } + } + } else { + assert(_N <= 5); + } + } else { + assert(_strategy == EnergyCorrelator::slow || _strategy == EnergyCorrelator::storage_array); + } + return answer; + } + + + std::vector EnergyCorrelatorGeneralized::result_all_angles(const PseudoJet& jet) const { + + // if jet does not have constituents, throw error + if (!jet.has_constituents()) throw Error("EnergyCorrelator called on jet with no constituents."); + + // Throw an error if N < 1 + if (_N < 1 ) throw Error("N cannot be negative or zero"); + + // get the N = 1 case out of the way + if (_N == 1) { + std::vector ans (1, 1.0); + return ans; + } + + // find constituents + std::vector particles = jet.constituents(); + + // return zero if the number of constituents is less than _N for the ECFG + if (particles.size() < _N) { + std::vector ans (_N, 0.0); + return ans; + } + + // The normalization is the energy or pt of the jet, which is also ECF(1, beta) + double EJ = _helper_correlator.result(jet); + + // The overall normalization + double norm = pow(EJ, _N); + + // Find the max number of angles and throw an error if it unsuitable + int N_total = _N * (_N - 1)/2; + + double half_beta = _beta/2.0; + + // take care of N = 2 case. + if (_N == 2) { + double answer = 0.0; + for (unsigned int i = 0; i < particles.size(); i++) { + for (unsigned int j = i + 1; j < particles.size(); j++) { //note offset by one so that angle is never called on identical pairs + answer += energy(particles[i]) + * energy(particles[j]) + * pow(angleSquared(particles[i],particles[j]), half_beta)/norm; + } + } + std::vector ans(N_total, answer); + return ans; + } + + // Prepare the answer vector + std::vector ans (N_total, 0.0); + // if N > 4, then throw error + if (_N > 5) { + throw Error("EnergyCorrelatorGeneralized is only hard coded for N = 0,1,2,3,4,5"); + } + + // Now deal with N = 3,4,5. Different options if storage array is used or not. + if (_strategy == EnergyCorrelator::storage_array) { + + // For N > 2, fill static storage array to save computation time. + + // Make energy storage + std::vector energyStore; + energyStore.resize(particles.size()); + + // Make angular storage + std::vector < std::vector > angleStore; + angleStore.resize(particles.size()); + for (unsigned int i = 0; i < angleStore.size(); i++) { + angleStore[i].resize(i); + } + + // Fill storage with energy/angle information + for (unsigned int i = 0; i < particles.size(); i++) { + energyStore[i] = energy(particles[i]); + for (unsigned int j = 0; j < i; j++) { + if (half_beta == 1){ + angleStore[i][j] = angleSquared(particles[i], particles[j]); + } else { + angleStore[i][j] = pow(angleSquared(particles[i], particles[j]), half_beta); + } + } + } + + // now do recursion + if (_N == 3) { + double angle1, angle2, angle3; + + for (unsigned int i = 2; i < particles.size(); i++) { + for (unsigned int j = 1; j < i; j++) { + for (unsigned int k = 0; k < j; k++) { + + angle1 = angleStore[i][j]; + angle2 = angleStore[i][k]; + angle3 = angleStore[j][k]; + + double angle_list[] = {angle1, angle2, angle3}; + std::vector angle_vector(angle_list, angle_list + N_total); + std::sort(angle_vector.begin(), angle_vector.begin() + N_total); + + std::vector final_angles (N_total, angle_vector[0]); + + double z_product = energyStore[i] * energyStore[j] * energyStore[k]/norm; + ans[0] += z_product * final_angles[0]; + for ( int s=1 ; s angle_vector(angle_list, angle_list + N_total); + std::sort(angle_vector.begin(), angle_vector.begin() + N_total); + + std::vector final_angles (N_total, angle_vector[0]); + + double z_product = energyStore[i] * energyStore[j] * energyStore[k] * energyStore [l]/norm; + ans[0] += z_product * final_angles[0]; + for ( int s=1 ; s angle_vector(angle_list, angle_list + N_total); + std::sort(angle_vector.begin(), angle_vector.begin() + N_total); + + std::vector final_angles (N_total, angle_vector[0]); + + double z_product = energyStore[i] * energyStore[j] * energyStore[k] + * energyStore[l] * energyStore[m]/norm; + ans[0] += z_product * final_angles[0]; + for ( int s=1 ; s angle_vector(angle_list, angle_list + N_total); + std::sort(angle_vector.begin(), angle_vector.begin() + N_total); + + std::vector final_angles (N_total, angle_vector[0]); + + double z_product = energy(particles[i]) + * energy(particles[j]) + * energy(particles[k])/norm; + + ans[0] += z_product * final_angles[0]; + for ( int s=1 ; s angle_vector(angle_list, angle_list + N_total); + std::sort(angle_vector.begin(), angle_vector.begin() + N_total); + + std::vector final_angles (N_total, angle_vector[0]); + + double z_product = energy(particles[i]) + * energy(particles[j]) + * energy(particles[k]) + * energy(particles[l])/norm; + ans[0] += z_product * final_angles[0]; + for ( int s=1 ; s angle_vector(angle_list, angle_list + N_total); + std::sort(angle_vector.begin(), angle_vector.begin() + N_total); + + std::vector final_angles (N_total, angle_vector[0]); + + double z_product = energy(particles[i]) + * energy(particles[j]) + * energy(particles[k]) + * energy(particles[l]) + * energy(particles[m])/norm; + + ans[0] += z_product * final_angles[0]; + for ( int s=1 ; s. +//---------------------------------------------------------------------- + +#include +#include "fastjet/FunctionOfPseudoJet.hh" + +#include +#include + +FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh + +namespace contrib{ + +/// \mainpage EnergyCorrelator contrib +/// +/// The EnergyCorrelator contrib provides an implementation of energy +/// correlators and their ratios as described in arXiv:1305.0007 by +/// Larkoski, Salam and Thaler. Additionally, the ratio observable +/// D2 described in arXiv:1409.6298 by Larkoski, Moult and Neill +/// is also included in this contrib. Finally, a generalized version of +/// the energy correlation functions is added, defined in +/// arXiv:1609.07483 by Moult, Necib and Thaler, which allow the +/// definition of the M series, N series, and U series observables. +/// There is also a generalized version of D2. +/// +/// +///

There are 4 main classes: +/// +/// - EnergyCorrelator +/// - EnergyCorrelatorRatio +/// - EnergyCorrelatorDoubleRatio +/// - EnergyCorrelatorGeneralized +/// +///

There are five classes that define useful combinations of the ECFs. +/// +/// - EnergyCorrelatorNseries +/// - EnergyCorrelatorMseries +/// - EnergyCorrelatorUseries +/// - EnergyCorrelatorD2 +/// - EnergyCorrelatorGeneralizedD2 +/// +///

There are also aliases for easier access: +/// - EnergyCorrelatorCseries (same as EnergyCorrelatorDoubleRatio) +/// - EnergyCorrelatorC1 (EnergyCorrelatorCseries with i=1) +/// - EnergyCorrelatorC2 (EnergyCorrelatorCseries with i=2) +/// - EnergyCorrelatorN2 (EnergyCorrelatorNseries with i=2) +/// - EnergyCorrelatorN3 (EnergyCorrelatorNseries with i=3) +/// - EnergyCorrelatorM2 (EnergyCorrelatorMseries with i=2) +/// - EnergyCorrelatorU1 (EnergyCorrelatorUseries with i=1) +/// - EnergyCorrelatorU2 (EnergyCorrelatorUseries with i=2) +/// - EnergyCorrelatorU3 (EnergyCorrelatorUseries with i=3) +/// +/// Each of these classes is a FastJet FunctionOfPseudoJet. +/// EnergyCorrelatorDoubleRatio (which is equivalent to EnergyCorrelatorCseries) +/// is in particular is useful for quark/gluon discrimination and boosted +/// object tagging. +/// +/// Using the original 2- and 3-point correlators, EnergyCorrelationD2 has +/// been shown to be the optimal combination for boosted 2-prong tagging. +/// +/// The EnergyCorrelatorNseries and EnergyCorrelatorMseries use +/// generalized correlation functions with different angular scaling, +/// and are intended for use on 2-prong and 3-prong jets. +/// The EnergyCorrelatorUseries is useful for quark/gluon discrimimation. +/// +/// See the file example.cc for an illustration of usage and +/// example_basic_usage.cc for the most commonly used functions. + +//------------------------------------------------------------------------ +/// \class EnergyCorrelator +/// ECF(N,beta) is the N-point energy correlation function, with an angular exponent beta. +/// +/// It is defined as follows +/// +/// - \f$ \mathrm{ECF}(1,\beta) = \sum_i E_i \f$ +/// - \f$ \mathrm{ECF}(2,\beta) = \sum_{i { + friend class EnergyCorrelatorGeneralized; ///< This allow ECFG to access the energy and angle definitions + ///< of this class, which are otherwise private. + public: + + enum Measure { + pt_R, ///< use transverse momenta and boost-invariant angles, + ///< eg \f$\mathrm{ECF}(2,\beta) = \sum_{i=3 this leads to many expensive recomputations, + ///< but has only O(n) memory usage for n particles + + storage_array /// the interparticle angles are cached. This gives a significant speed + /// improvement for N>=3, but has a memory requirement of (4n^2) bytes. + }; + + public: + + /// constructs an N-point correlator with angular exponent beta, + /// using the specified choice of energy and angular measure as well + /// one of two possible underlying computational Strategy + EnergyCorrelator(unsigned int N, + double beta, + Measure measure = pt_R, + Strategy strategy = storage_array) : + _N(N), _beta(beta), _measure(measure), _strategy(strategy) {}; + + /// destructor + virtual ~EnergyCorrelator(){} + + /// returns the value of the energy correlator for a jet's + /// constituents. (Normally accessed by the parent class's + /// operator()). + double result(const PseudoJet& jet) const; + + std::string description() const; + + /// returns the the part of the description related to the parameters + std::string description_parameters() const; + std::string description_no_N() const; + + private: + + unsigned int _N; + double _beta; + Measure _measure; + Strategy _strategy; + + double energy(const PseudoJet& jet) const; + double angleSquared(const PseudoJet& jet1, const PseudoJet& jet2) const; + double multiply_angles(double angles[], int n_angles, unsigned int N_total) const; + void precompute_energies_and_angles(std::vector const &particles, double* energyStore, double** angleStore) const; + double evaluate_n3(unsigned int nC, unsigned int n_angles, double* energyStore, double** angleStore) const; + double evaluate_n4(unsigned int nC, unsigned int n_angles, double* energyStore, double** angleStore) const; + double evaluate_n5(unsigned int nC, unsigned int n_angles, double* energyStore, double** angleStore) const; + }; + +// core EnergyCorrelator::result code in .cc file. + + +//------------------------------------------------------------------------ +/// \class EnergyCorrelatorRatio +/// A class to calculate the ratio of (N+1)-point to N-point energy correlators, +/// ECF(N+1,beta)/ECF(N,beta), +/// called \f$ r_N^{(\beta)} \f$ in the publication. + class EnergyCorrelatorRatio : public FunctionOfPseudoJet { + + public: + + /// constructs an (N+1)-point to N-point correlator ratio with + /// angular exponent beta, using the specified choice of energy and + /// angular measure as well one of two possible underlying + /// computational strategies + EnergyCorrelatorRatio(unsigned int N, + double beta, + EnergyCorrelator::Measure measure = EnergyCorrelator::pt_R, + EnergyCorrelator::Strategy strategy = EnergyCorrelator::storage_array) + : _N(N), _beta(beta), _measure(measure), _strategy(strategy) {}; + + virtual ~EnergyCorrelatorRatio() {} + + /// returns the value of the energy correlator ratio for a jet's + /// constituents. (Normally accessed by the parent class's + /// operator()). + double result(const PseudoJet& jet) const; + + std::string description() const; + + private: + + unsigned int _N; + double _beta; + + EnergyCorrelator::Measure _measure; + EnergyCorrelator::Strategy _strategy; + + + }; + + inline double EnergyCorrelatorRatio::result(const PseudoJet& jet) const { + + double numerator = EnergyCorrelator(_N + 1, _beta, _measure, _strategy).result(jet); + double denominator = EnergyCorrelator(_N, _beta, _measure, _strategy).result(jet); + + return numerator/denominator; + + } + + +//------------------------------------------------------------------------ +/// \class EnergyCorrelatorDoubleRatio +/// Calculates the double ratio of energy correlators, ECF(N-1,beta)*ECF(N+1)/ECF(N,beta)^2. +/// +/// A class to calculate a double ratio of energy correlators, +/// ECF(N-1,beta)*ECF(N+1,beta)/ECF(N,beta)^2, +/// called \f$C_N^{(\beta)}\f$ in the publication, and equal to +/// \f$ r_N^{(\beta)}/r_{N-1}^{(\beta)} \f$. +/// + + class EnergyCorrelatorDoubleRatio : public FunctionOfPseudoJet { + + public: + + EnergyCorrelatorDoubleRatio(unsigned int N, + double beta, + EnergyCorrelator::Measure measure = EnergyCorrelator::pt_R, + EnergyCorrelator::Strategy strategy = EnergyCorrelator::storage_array) + : _N(N), _beta(beta), _measure(measure), _strategy(strategy) { + + if (_N < 1) throw Error("EnergyCorrelatorDoubleRatio: N must be 1 or greater."); + + }; + + virtual ~EnergyCorrelatorDoubleRatio() {} + + + /// returns the value of the energy correlator double-ratio for a + /// jet's constituents. (Normally accessed by the parent class's + /// operator()). + double result(const PseudoJet& jet) const; + + std::string description() const; + + private: + + unsigned int _N; + double _beta; + EnergyCorrelator::Measure _measure; + EnergyCorrelator::Strategy _strategy; + + + }; + + + inline double EnergyCorrelatorDoubleRatio::result(const PseudoJet& jet) const { + + double numerator = EnergyCorrelator(_N - 1, _beta, _measure, _strategy).result(jet) * EnergyCorrelator(_N + 1, _beta, _measure, _strategy).result(jet); + double denominator = pow(EnergyCorrelator(_N, _beta, _measure, _strategy).result(jet), 2.0); + + return numerator/denominator; + + } + +//------------------------------------------------------------------------ +/// \class EnergyCorrelatorC1 +/// A class to calculate the normalized 2-point energy correlators, +/// ECF(2,beta)/ECF(1,beta)^2, +/// called \f$ C_1^{(\beta)} \f$ in the publication. + class EnergyCorrelatorC1 : public FunctionOfPseudoJet { + + public: + + /// constructs a 2-point correlator ratio with + /// angular exponent beta, using the specified choice of energy and + /// angular measure as well one of two possible underlying + /// computational strategies + EnergyCorrelatorC1(double beta, + EnergyCorrelator::Measure measure = EnergyCorrelator::pt_R, + EnergyCorrelator::Strategy strategy = EnergyCorrelator::storage_array) + : _beta(beta), _measure(measure), _strategy(strategy) {}; + + virtual ~EnergyCorrelatorC1() {} + + /// returns the value of the energy correlator ratio for a jet's + /// constituents. (Normally accessed by the parent class's + /// operator()). + double result(const PseudoJet& jet) const; + + std::string description() const; + + private: + + double _beta; + + EnergyCorrelator::Measure _measure; + EnergyCorrelator::Strategy _strategy; + + + }; + + + inline double EnergyCorrelatorC1::result(const PseudoJet& jet) const { + + double numerator = EnergyCorrelator(2, _beta, _measure, _strategy).result(jet); + double denominator = EnergyCorrelator(1, _beta, _measure, _strategy).result(jet); + + return numerator/denominator/denominator; + + } + + +//------------------------------------------------------------------------ +/// \class EnergyCorrelatorC2 +/// A class to calculate the double ratio of 3-point to 2-point +/// energy correlators, +/// ECF(3,beta)*ECF(1,beta)/ECF(2,beta)^2, +/// called \f$ C_2^{(\beta)} \f$ in the publication. + class EnergyCorrelatorC2 : public FunctionOfPseudoJet { + + public: + + /// constructs a 3-point to 2-point correlator double ratio with + /// angular exponent beta, using the specified choice of energy and + /// angular measure as well one of two possible underlying + /// computational strategies + EnergyCorrelatorC2(double beta, + EnergyCorrelator::Measure measure = EnergyCorrelator::pt_R, + EnergyCorrelator::Strategy strategy = EnergyCorrelator::storage_array) + : _beta(beta), _measure(measure), _strategy(strategy) {}; + + virtual ~EnergyCorrelatorC2() {} + + /// returns the value of the energy correlator ratio for a jet's + /// constituents. (Normally accessed by the parent class's + /// operator()). + double result(const PseudoJet& jet) const; + + std::string description() const; + + private: + + double _beta; + + EnergyCorrelator::Measure _measure; + EnergyCorrelator::Strategy _strategy; + + + }; + + + inline double EnergyCorrelatorC2::result(const PseudoJet& jet) const { + + double numerator3 = EnergyCorrelator(3, _beta, _measure, _strategy).result(jet); + double numerator1 = EnergyCorrelator(1, _beta, _measure, _strategy).result(jet); + double denominator = EnergyCorrelator(2, _beta, _measure, _strategy).result(jet); + + return numerator3*numerator1/denominator/denominator; + + } + + +//------------------------------------------------------------------------ +/// \class EnergyCorrelatorD2 +/// A class to calculate the observable formed from the ratio of the +/// 3-point and 2-point energy correlators, +/// ECF(3,beta)*ECF(1,beta)^3/ECF(2,beta)^3, +/// called \f$ D_2^{(\beta)} \f$ in the publication. + class EnergyCorrelatorD2 : public FunctionOfPseudoJet { + + public: + + /// constructs an 3-point to 2-point correlator ratio with + /// angular exponent beta, using the specified choice of energy and + /// angular measure as well one of two possible underlying + /// computational strategies + EnergyCorrelatorD2(double beta, + EnergyCorrelator::Measure measure = EnergyCorrelator::pt_R, + EnergyCorrelator::Strategy strategy = EnergyCorrelator::storage_array) + : _beta(beta), _measure(measure), _strategy(strategy) {}; + + virtual ~EnergyCorrelatorD2() {} + + /// returns the value of the energy correlator ratio for a jet's + /// constituents. (Normally accessed by the parent class's + /// operator()). + double result(const PseudoJet& jet) const; + + std::string description() const; + + private: + + double _beta; + + EnergyCorrelator::Measure _measure; + EnergyCorrelator::Strategy _strategy; + + + }; + + + inline double EnergyCorrelatorD2::result(const PseudoJet& jet) const { + + double numerator3 = EnergyCorrelator(3, _beta, _measure, _strategy).result(jet); + double numerator1 = EnergyCorrelator(1, _beta, _measure, _strategy).result(jet); + double denominator2 = EnergyCorrelator(2, _beta, _measure, _strategy).result(jet); + + return numerator3*numerator1*numerator1*numerator1/denominator2/denominator2/denominator2; + + } + + +//------------------------------------------------------------------------ +/// \class EnergyCorrelatorGeneralized +/// A generalized and normalized version of the N-point energy correlators, with +/// angular exponent beta and v number of pairwise angles. When \f$v = {N \choose 2}\f$ +/// (or, for convenience, \f$v = -1\f$), EnergyCorrelatorGeneralized just gives normalized +/// versions of EnergyCorrelator: +/// - \f$ \mathrm{ECFG}(-1,1,\beta) = \mathrm{ECFN}(N,\beta) = \mathrm{ECF}(N,\beta)/\mathrm{ECF}(1,\beta)\f$ +/// +/// Note that there is no separate class that implements ECFN, though it is a +/// notation that we will use in this documentation. Examples of the low-point normalized +/// correlators are: +/// - \f$\mathrm{ECFN}(1,\beta) = 1\f$ +/// - \f$\mathrm{ECFN}(2,\beta) = \sum_{i { + public: + + /// constructs an N-point correlator with v_angles pairwise angles + /// and angular exponent beta, + /// using the specified choice of energy and angular measure as well + /// one of two possible underlying computational Strategy + EnergyCorrelatorGeneralized(int v_angles, + unsigned int N, + double beta, + EnergyCorrelator::Measure measure = EnergyCorrelator::pt_R, + EnergyCorrelator::Strategy strategy = EnergyCorrelator::storage_array) + : _angles(v_angles), _N(N), _beta(beta), _measure(measure), _strategy(strategy), _helper_correlator(1,_beta, _measure, _strategy) {}; + + /// destructor + virtual ~EnergyCorrelatorGeneralized(){} + + /// returns the value of the normalized energy correlator for a jet's + /// constituents. (Normally accessed by the parent class's + /// operator()). + + double result(const PseudoJet& jet) const; + std::vector result_all_angles(const PseudoJet& jet) const; + + private: + + int _angles; + unsigned int _N; + double _beta; + EnergyCorrelator::Measure _measure; + EnergyCorrelator::Strategy _strategy; + EnergyCorrelator _helper_correlator; + + double energy(const PseudoJet& jet) const; + double angleSquared(const PseudoJet& jet1, const PseudoJet& jet2) const; + double multiply_angles(double angles[], int n_angles, unsigned int N_total) const; + void precompute_energies_and_angles(std::vector const &particles, double* energyStore, double** angleStore) const; + double evaluate_n3(unsigned int nC, unsigned int n_angles, double* energyStore, double** angleStore) const; + double evaluate_n4(unsigned int nC, unsigned int n_angles, double* energyStore, double** angleStore) const; + double evaluate_n5(unsigned int nC, unsigned int n_angles, double* energyStore, double** angleStore) const; + }; + + + +//------------------------------------------------------------------------ +/// \class EnergyCorrelatorGeneralizedD2 +/// A class to calculate the observable formed from the ratio of the +/// 3-point and 2-point energy correlators, +/// ECFN(3,alpha)/ECFN(2,beta)^3 alpha/beta, +/// called \f$ D_2^{(\alpha, \beta)} \f$ in the publication. + class EnergyCorrelatorGeneralizedD2 : public FunctionOfPseudoJet { + + public: + + /// constructs an 3-point to 2-point correlator ratio with + /// angular exponent beta, using the specified choice of energy and + /// angular measure as well one of two possible underlying + /// computational strategies + EnergyCorrelatorGeneralizedD2( + double alpha, + double beta, + EnergyCorrelator::Measure measure = EnergyCorrelator::pt_R, + EnergyCorrelator::Strategy strategy = EnergyCorrelator::storage_array) + : _alpha(alpha), _beta(beta), _measure(measure), _strategy(strategy) {}; + + virtual ~EnergyCorrelatorGeneralizedD2() {} + + /// returns the value of the energy correlator ratio for a jet's + /// constituents. (Normally accessed by the parent class's + /// operator()). + double result(const PseudoJet& jet) const; + + std::string description() const; + + private: + + double _alpha; + double _beta; + + EnergyCorrelator::Measure _measure; + EnergyCorrelator::Strategy _strategy; + + + }; + + + inline double EnergyCorrelatorGeneralizedD2::result(const PseudoJet& jet) const { + + double numerator = EnergyCorrelatorGeneralized(-1, 3, _alpha, _measure, _strategy).result(jet); + double denominator = EnergyCorrelatorGeneralized(-1, 2, _beta, _measure, _strategy).result(jet); + + return numerator/pow(denominator, 3.0*_alpha/_beta); + + } + + +//------------------------------------------------------------------------ +/// \class EnergyCorrelatorNseries +/// A class to calculate the observable formed from the ratio of the +/// 3-point and 2-point energy correlators, +/// N_n = ECFG(2,n+1,beta)/ECFG(1,n,beta)^2, +/// called \f$ N_i^{(\alpha, \beta)} \f$ in the publication. +/// By definition, N_1^{beta} = ECFG(1, 2, 2*beta), where the angular exponent +/// is twice as big since the N series should involve two pairwise angles. + class EnergyCorrelatorNseries : public FunctionOfPseudoJet { + + public: + + /// constructs a n 3-point to 2-point correlator ratio with + /// angular exponent beta, using the specified choice of energy and + /// angular measure as well one of two possible underlying + /// computational strategies + EnergyCorrelatorNseries( + unsigned int n, + double beta, + EnergyCorrelator::Measure measure = EnergyCorrelator::pt_R, + EnergyCorrelator::Strategy strategy = EnergyCorrelator::storage_array) + : _n(n), _beta(beta), _measure(measure), _strategy(strategy) { + + if (_n < 1) throw Error("EnergyCorrelatorNseries: n must be 1 or greater."); + + }; + + virtual ~EnergyCorrelatorNseries() {} + + /// returns the value of the energy correlator ratio for a jet's + /// constituents. (Normally accessed by the parent class's + /// operator()). + double result(const PseudoJet& jet) const; + + std::string description() const; + + private: + + unsigned int _n; + double _beta; + EnergyCorrelator::Measure _measure; + EnergyCorrelator::Strategy _strategy; + + }; + + + inline double EnergyCorrelatorNseries::result(const PseudoJet& jet) const { + + if (_n == 1) return EnergyCorrelatorGeneralized(1, 2, 2*_beta, _measure, _strategy).result(jet); + // By definition, N1 = ECFN(2, 2 beta) + double numerator = EnergyCorrelatorGeneralized(2, _n + 1, _beta, _measure, _strategy).result(jet); + double denominator = EnergyCorrelatorGeneralized(1, _n, _beta, _measure, _strategy).result(jet); + + return numerator/denominator/denominator; + + } + + + +//------------------------------------------------------------------------ +/// \class EnergyCorrelatorN2 +/// A class to calculate the observable formed from the ratio of the +/// 3-point and 2-point energy correlators, +/// ECFG(2,3,beta)/ECFG(1,2,beta)^2, +/// called \f$ N_2^{(\beta)} \f$ in the publication. + class EnergyCorrelatorN2 : public FunctionOfPseudoJet { + + public: + + /// constructs an 3-point to 2-point correlator ratio with + /// angular exponent beta, using the specified choice of energy and + /// angular measure as well one of two possible underlying + /// computational strategies + EnergyCorrelatorN2(double beta, + EnergyCorrelator::Measure measure = EnergyCorrelator::pt_R, + EnergyCorrelator::Strategy strategy = EnergyCorrelator::storage_array) + : _beta(beta), _measure(measure), _strategy(strategy) {}; + + virtual ~EnergyCorrelatorN2() {} + + /// returns the value of the energy correlator ratio for a jet's + /// constituents. (Normally accessed by the parent class's + /// operator()). + double result(const PseudoJet& jet) const; + + std::string description() const; + + private: + + double _beta; + + EnergyCorrelator::Measure _measure; + EnergyCorrelator::Strategy _strategy; + + + }; + + + inline double EnergyCorrelatorN2::result(const PseudoJet& jet) const { + + double numerator = EnergyCorrelatorGeneralized(2, 3, _beta, _measure, _strategy).result(jet); + double denominator = EnergyCorrelatorGeneralized(1, 2, _beta, _measure, _strategy).result(jet); + + return numerator/denominator/denominator; + + } + + +//------------------------------------------------------------------------ +/// \class EnergyCorrelatorN3 +/// A class to calculate the observable formed from the ratio of the +/// 3-point and 2-point energy correlators, +/// ECFG(2,4,beta)/ECFG(1,3,beta)^2, +/// called \f$ N_3^{(\beta)} \f$ in the publication. + class EnergyCorrelatorN3 : public FunctionOfPseudoJet { + + public: + + /// constructs an 3-point to 2-point correlator ratio with + /// angular exponent beta, using the specified choice of energy and + /// angular measure as well one of two possible underlying + /// computational strategies + EnergyCorrelatorN3(double beta, + EnergyCorrelator::Measure measure = EnergyCorrelator::pt_R, + EnergyCorrelator::Strategy strategy = EnergyCorrelator::storage_array) + : _beta(beta), _measure(measure), _strategy(strategy) {}; + + virtual ~EnergyCorrelatorN3() {} + + /// returns the value of the energy correlator ratio for a jet's + /// constituents. (Normally accessed by the parent class's + /// operator()). + double result(const PseudoJet& jet) const; + + std::string description() const; + + private: + + double _beta; + + EnergyCorrelator::Measure _measure; + EnergyCorrelator::Strategy _strategy; + + + }; + + + inline double EnergyCorrelatorN3::result(const PseudoJet& jet) const { + + double numerator = EnergyCorrelatorGeneralized(2, 4, _beta, _measure, _strategy).result(jet); + double denominator = EnergyCorrelatorGeneralized(1, 3, _beta, _measure, _strategy).result(jet); + + return numerator/denominator/denominator; + + } + + +//------------------------------------------------------------------------ +/// \class EnergyCorrelatorMseries +/// A class to calculate the observable formed from the ratio of the +/// 3-point and 2-point energy correlators, +/// M_n = ECFG(1,n+1,beta)/ECFG(1,n,beta), +/// called \f$ M_i^{(\alpha, \beta)} \f$ in the publication. +/// By definition, M_1^{beta} = ECFG(1,2,beta) + class EnergyCorrelatorMseries : public FunctionOfPseudoJet { + + public: + + /// constructs a n 3-point to 2-point correlator ratio with + /// angular exponent beta, using the specified choice of energy and + /// angular measure as well one of two possible underlying + /// computational strategies + EnergyCorrelatorMseries( + unsigned int n, + double beta, + EnergyCorrelator::Measure measure = EnergyCorrelator::pt_R, + EnergyCorrelator::Strategy strategy = EnergyCorrelator::storage_array) + : _n(n), _beta(beta), _measure(measure), _strategy(strategy) { + + if (_n < 1) throw Error("EnergyCorrelatorMseries: n must be 1 or greater."); + + }; + + virtual ~EnergyCorrelatorMseries() {} + + /// returns the value of the energy correlator ratio for a jet's + /// constituents. (Normally accessed by the parent class's + /// operator()). + double result(const PseudoJet& jet) const; + + std::string description() const; + + private: + + unsigned int _n; + double _beta; + EnergyCorrelator::Measure _measure; + EnergyCorrelator::Strategy _strategy; + + }; + + + inline double EnergyCorrelatorMseries::result(const PseudoJet& jet) const { + + if (_n == 1) return EnergyCorrelatorGeneralized(1, 2, _beta, _measure, _strategy).result(jet); + + double numerator = EnergyCorrelatorGeneralized(1, _n + 1, _beta, _measure, _strategy).result(jet); + double denominator = EnergyCorrelatorGeneralized(1, _n, _beta, _measure, _strategy).result(jet); + + return numerator/denominator; + + } + +//------------------------------------------------------------------------ +/// \class EnergyCorrelatorM2 +/// A class to calculate the observable formed from the ratio of the +/// 3-point and 2-point energy correlators, +/// ECFG(1,3,beta)/ECFG(1,2,beta), +/// called \f$ M_2^{(\beta)} \f$ in the publication. + class EnergyCorrelatorM2 : public FunctionOfPseudoJet { + + public: + + /// constructs an 3-point to 2-point correlator ratio with + /// angular exponent beta, using the specified choice of energy and + /// angular measure as well one of two possible underlying + /// computational strategies + EnergyCorrelatorM2(double beta, + EnergyCorrelator::Measure measure = EnergyCorrelator::pt_R, + EnergyCorrelator::Strategy strategy = EnergyCorrelator::storage_array) + : _beta(beta), _measure(measure), _strategy(strategy) {}; + + virtual ~EnergyCorrelatorM2() {} + + /// returns the value of the energy correlator ratio for a jet's + /// constituents. (Normally accessed by the parent class's + /// operator()). + double result(const PseudoJet& jet) const; + + std::string description() const; + + private: + + double _beta; + + EnergyCorrelator::Measure _measure; + EnergyCorrelator::Strategy _strategy; + + + }; + + + inline double EnergyCorrelatorM2::result(const PseudoJet& jet) const { + + double numerator = EnergyCorrelatorGeneralized(1, 3, _beta, _measure, _strategy).result(jet); + double denominator = EnergyCorrelatorGeneralized(1, 2, _beta, _measure, _strategy).result(jet); + + return numerator/denominator; + + } + + +//------------------------------------------------------------------------ +/// \class EnergyCorrelatorCseries +/// Calculates the C series energy correlators, ECFN(N-1,beta)*ECFN(N+1,beta)/ECFN(N,beta)^2. +/// This is equivalent to EnergyCorrelatorDoubleRatio +/// +/// A class to calculate a double ratio of energy correlators, +/// ECFN(N-1,beta)*ECFN(N+1,beta)/ECFN(N,beta)^2, +/// called \f$C_N^{(\beta)}\f$ in the publication, and equal to +/// \f$ r_N^{(\beta)}/r_{N-1}^{(\beta)} \f$. +/// + + class EnergyCorrelatorCseries : public FunctionOfPseudoJet { + + public: + + EnergyCorrelatorCseries(unsigned int N, + double beta, + EnergyCorrelator::Measure measure = EnergyCorrelator::pt_R, + EnergyCorrelator::Strategy strategy = EnergyCorrelator::storage_array) + : _N(N), _beta(beta), _measure(measure), _strategy(strategy) { + + if (_N < 1) throw Error("EnergyCorrelatorCseries: N must be 1 or greater."); + + }; + + virtual ~EnergyCorrelatorCseries() {} + + + /// returns the value of the energy correlator double-ratio for a + /// jet's constituents. (Normally accessed by the parent class's + /// operator()). + double result(const PseudoJet& jet) const; + + std::string description() const; + + private: + + unsigned int _N; + double _beta; + EnergyCorrelator::Measure _measure; + EnergyCorrelator::Strategy _strategy; + + + }; + + + inline double EnergyCorrelatorCseries::result(const PseudoJet& jet) const { + + double numerator = EnergyCorrelatorGeneralized(-1, _N - 1, _beta, _measure, _strategy).result(jet) * EnergyCorrelatorGeneralized(-1, _N + 1, _beta, _measure, _strategy).result(jet); + double denominator = pow(EnergyCorrelatorGeneralized(-1, _N, _beta, _measure, _strategy).result(jet), 2.0); + + return numerator/denominator; + + } + +//------------------------------------------------------------------------ +/// \class EnergyCorrelatorUseries +/// A class to calculate the observable used for quark versus gluon discrimination +/// U_n = ECFG(1,n+1,beta), +/// called \f$ U_i^{(\beta)} \f$ in the publication. + + class EnergyCorrelatorUseries : public FunctionOfPseudoJet { + + public: + + /// constructs a n 3-point to 2-point correlator ratio with + /// angular exponent beta, using the specified choice of energy and + /// angular measure as well one of two possible underlying + /// computational strategies + EnergyCorrelatorUseries( + unsigned int n, + double beta, + EnergyCorrelator::Measure measure = EnergyCorrelator::pt_R, + EnergyCorrelator::Strategy strategy = EnergyCorrelator::storage_array) + : _n(n), _beta(beta), _measure(measure), _strategy(strategy) { + + if (_n < 1) throw Error("EnergyCorrelatorUseries: n must be 1 or greater."); + + }; + + virtual ~EnergyCorrelatorUseries() {} + + /// returns the value of the energy correlator ratio for a jet's + /// constituents. (Normally accessed by the parent class's + /// operator()). + double result(const PseudoJet& jet) const; + + std::string description() const; + + private: + + unsigned int _n; + double _beta; + EnergyCorrelator::Measure _measure; + EnergyCorrelator::Strategy _strategy; + + }; + + + inline double EnergyCorrelatorUseries::result(const PseudoJet& jet) const { + + double answer = EnergyCorrelatorGeneralized(1, _n + 1, _beta, _measure, _strategy).result(jet); + return answer; + + } + + +//------------------------------------------------------------------------ +/// \class EnergyCorrelatorU1 +/// A class to calculate the observable formed from +/// ECFG(1,2,beta), +/// called \f$ U_1^{(\beta)} \f$ in the publication. + class EnergyCorrelatorU1 : public FunctionOfPseudoJet { + + public: + + /// constructs a 2-point correlator with + /// angular exponent beta, using the specified choice of energy and + /// angular measure as well one of two possible underlying + /// computational strategies + EnergyCorrelatorU1(double beta, + EnergyCorrelator::Measure measure = EnergyCorrelator::pt_R, + EnergyCorrelator::Strategy strategy = EnergyCorrelator::storage_array) + : _beta(beta), _measure(measure), _strategy(strategy) {}; + + virtual ~EnergyCorrelatorU1() {} + + /// returns the value of the energy correlator ratio for a jet's + /// constituents. (Normally accessed by the parent class's + /// operator()). + double result(const PseudoJet& jet) const; + + std::string description() const; + + private: + + double _beta; + + EnergyCorrelator::Measure _measure; + EnergyCorrelator::Strategy _strategy; + + + }; + + + inline double EnergyCorrelatorU1::result(const PseudoJet& jet) const { + + double answer = EnergyCorrelatorGeneralized(1, 2, _beta, _measure, _strategy).result(jet); + + return answer; + + } + + + //------------------------------------------------------------------------ + /// \class EnergyCorrelatorU2 + /// A class to calculate the observable formed from + /// ECFG(1,3,beta), + /// called \f$ U_2^{(\beta)} \f$ in the publication. + class EnergyCorrelatorU2 : public FunctionOfPseudoJet { + + public: + + /// constructs a 3-point correlator with + /// angular exponent beta, using the specified choice of energy and + /// angular measure as well one of two possible underlying + /// computational strategies + EnergyCorrelatorU2(double beta, + EnergyCorrelator::Measure measure = EnergyCorrelator::pt_R, + EnergyCorrelator::Strategy strategy = EnergyCorrelator::storage_array) + : _beta(beta), _measure(measure), _strategy(strategy) {}; + + virtual ~EnergyCorrelatorU2() {} + + /// returns the value of the energy correlator ratio for a jet's + /// constituents. (Normally accessed by the parent class's + /// operator()). + double result(const PseudoJet& jet) const; + + std::string description() const; + + private: + + double _beta; + + EnergyCorrelator::Measure _measure; + EnergyCorrelator::Strategy _strategy; + + + }; + + + inline double EnergyCorrelatorU2::result(const PseudoJet& jet) const { + + double answer = EnergyCorrelatorGeneralized(1, 3, _beta, _measure, _strategy).result(jet); + + return answer; + + } + + + //------------------------------------------------------------------------ + /// \class EnergyCorrelatorU3 + /// A class to calculate the observable formed from + /// ECFG(1,4,beta), + /// called \f$ U_3^{(\beta)} \f$ in the publication. + class EnergyCorrelatorU3 : public FunctionOfPseudoJet { + + public: + + /// constructs a 4-point correlator with + /// angular exponent beta, using the specified choice of energy and + /// angular measure as well one of two possible underlying + /// computational strategies + EnergyCorrelatorU3(double beta, + EnergyCorrelator::Measure measure = EnergyCorrelator::pt_R, + EnergyCorrelator::Strategy strategy = EnergyCorrelator::storage_array) + : _beta(beta), _measure(measure), _strategy(strategy) {}; + + virtual ~EnergyCorrelatorU3() {} + + /// returns the value of the energy correlator ratio for a jet's + /// constituents. (Normally accessed by the parent class's + /// operator()). + double result(const PseudoJet& jet) const; + + std::string description() const; + + private: + + double _beta; + + EnergyCorrelator::Measure _measure; + EnergyCorrelator::Strategy _strategy; + + + }; + + + inline double EnergyCorrelatorU3::result(const PseudoJet& jet) const { + + double answer = EnergyCorrelatorGeneralized(1, 4, _beta, _measure, _strategy).result(jet); + + return answer; + + } + + + +} // namespace contrib + +FASTJET_END_NAMESPACE + +#endif // __FASTJET_CONTRIB_ENERGYCORRELATOR_HH__ diff --git a/src/Tools/fjcontrib/EnergyCorrelator/NEWS b/src/Tools/fjcontrib/EnergyCorrelator/NEWS new file mode 100644 --- /dev/null +++ b/src/Tools/fjcontrib/EnergyCorrelator/NEWS @@ -0,0 +1,14 @@ +2018-02-08: Release of version 1.3.1: + Fixing memory leak, functionality unchanged. +2018-01-09: Release of version 1.3.0: + Speed up of core ECF code, functionality unchanged. +2016-10-07: Release of version 1.2.0: + Incorporates the generalized energy correlation functions + and adds the NSeries, MSeries, USeries, and GeneralizedD2 + observables. +2014-11-13: Release of version 1.1.0: + Added the C1, C2, and D2 observables +2013-05-01: Release of version 1.0.1: + Improved error reporting +2013-04-30: Release of version 1.0.0: + Initial release diff --git a/src/Tools/fjcontrib/EnergyCorrelator/README b/src/Tools/fjcontrib/EnergyCorrelator/README new file mode 100644 --- /dev/null +++ b/src/Tools/fjcontrib/EnergyCorrelator/README @@ -0,0 +1,190 @@ +The EnergyCorrelator package is based on the physics described in: + + Energy Correlation Functions for Jet Substructure. + Andrew J. Larkoski, Gavin Salam, and Jesse Thaler. + JHEP 1306, 108 (2013) + arXiv:1305.0007. + +Additional information and a new observable formed from the +energy correlation functions was described in + + Power Counting to Better Jet Observables. + Andrew J. Larkoski, Ian Moult, and Duff Neill. + JHEP 1412, 009 (2014) + arXiv:1409.6298. + +Additional observables based on generalizations of the energy +correlation functions are described in + + New Angles on Energy Correlation Functions. + Ian Moult, Lina Necib, and Jesse Thaler. + arXiv:1609.07483. + +This FastJet-contrib package contains a number of classes derived from +FunctionOfPseudoJet. + +---------------------------------------------------------------------------- + +The core classes from 1305.0007, and defined since version 1.0, are: + +EnergyCorrelator(int N, double beta, Measure measure) + + Called ECF(N,beta) in arXiv:1305.0007. Corresponds to the N-point + correlation function, with beta the angular exponent, while measure + = pt_R (default) or E_theta sets how energies and angles are + determined. + +EnergyCorrelatorRatio(int N, double beta, Measure measure) + + Called r_N^(beta) in arXiv:1305.0007. + Equals ECF(N+1,beta)/ECF(N,beta). + +EnergyCorrelatorDoubleRatio(int N, double beta, Measure measure) + + Called C_N^(beta) in arXiv:1305.0007. Equals r_N/r_{N-1}. This + observable provides good boosted N-prong object discrimination. + (N=1 for quark/gluon, N=2 for boosted W/Z/H, N=3 for boosted top) + Also given in EnergyCorrelatorCseries as of version 1.2. + +---------------------------------------------------------------------------- + +The D2 observable from 1409.6298, as well as C1 and C2 alias classes, were +added in version 1.1: + +EnergyCorrelatorC1(double beta, Measure measure) + + This calculates the double ratio observable C_1^(beta) which is + useful for quark versus gluon discrimination. + +EnergyCorrelatorC2(double beta, Measure measure) + + This calculates the double ratio observable C_2^(beta) which is + useful for boosted W/Z/H identification. + +EnergyCorrelatorD2(double beta, Measure measure) + + Called D_2^(beta) in arXiv:1409.6298. + Equals ECF(3,beta)*ECF(1,beta)^3/ECF(2,beta)^3. + This is the recommended function for boosted 2-prong object + discrimination (boosted W/Z/H). + +---------------------------------------------------------------------------- + +Generalized energy correlators were introduced in 1609.07483 and appear in +version 1.2. They are defined in the class: + +EnergyCorrelatorGeneralized(int angles, int N, double beta, Measure measure) + + Called {}_v e_n^{(beta)} in 1609.07483, but will be denoted here as + ECFG(angles,N,beta), where v=angles and n=N. As for EnergyCorrelator, + beta is the angular exponent, while measure = pt_R (default) or E_theta + sets how energies and angles are determined. The integer angles + determines the number of angles in the observable. The choice angles=-1 + sets angles = N choose 2, which corresponds to the N-point + normalized (dimensionless) correlation function, with + ECFN(N,beta) = ECFG(N choose 2,N,beta) = ECF(N,beta)/ECF(1,beta)^N + +From the generalized correlators, a variety of useful ratios are defined. +They are mainly organized by series, with special values highlighted for +recommended usage. + +---------------------------------------------------------------------------- + +EnergyCorrelatorGeneralizedD2(double alpha, double beta, Measure measure) + + Called D_2^(alpha, beta) in arXiv:1609.07483 + Equals ECFN(3,alpha)/ECFN(2,beta)^(3 alpha/beta). + Useful for groomed 2-prong object tagging. We recommend the use of alpha=1 + and beta=2. + +---------------------------------------------------------------------------- + +EnergyCorrelatorNseries(int i, double beta, Measure measure) + + Called N_i^(beta) in arXiv:1609.07483 + Equals ECFG(2,n+1,beta)/ECFN(1,n,beta)^2. + +EnergyCorrelatorN2(double beta, Measure measure) + + Called N_2^(beta) in arXiv:1609.07483 + Equals ECFG(2,3,beta)/ECFG(1,2,beta)^2. + Useful for groomed and ungroomed 2-prong object tagging. + +EnergyCorrelatorN3(double beta, Measure measure) + + Called N_3^(beta) in arXiv:1609.07483 + Equals ECFG(2,4,beta)/ECFG(1,3,beta)^2. + Useful for groomed 3-prong object tagging. + +---------------------------------------------------------------------------- + +EnergyCorrelatorMseries(int i, double beta, Measure measure) + + Called M_i^(beta) in arXiv:1609.07483 + Equals ECFG(1,n+1,beta)/ECFG(1,n,beta). + +EnergyCorrelatorM2(double beta, Measure measure) + + Called M_2^(beta) in arXiv:1609.07483 + Equals ECFG(1,3,beta)/ECFG(1,2,beta). + Useful for groomed 2-prong object tagging. + + +---------------------------------------------------------------------------- + +EnergyCorrelatorUseries(int i, double beta, Measure measure) + + Called U_i^(beta) in arXiv:1609.07483 + Equals ECFG(1,n+1,beta). + +EnergyCorrelatorU1(double beta, Measure measure) + + Called U_1^(beta) in arXiv:1609.07483 + Equals ECFG(1,2,beta). + Useful for quark vs. gluon discrimination. + +EnergyCorrelatorU2(double beta, Measure measure) + + Called U_2^(beta) in arXiv:1609.07483 + Equals ECFG(1,3,beta). + Useful for quark vs. gluon discrimination. + +EnergyCorrelatorU3(double beta, Measure measure) + + Called U_3^(beta) in arXiv:1609.07483 + Equals ECFG(1,4,beta). + Useful for quark vs. gluon discrimination. +---------------------------------------------------------------------------- + +The argument Measure in each of the above functions sets how energies +and angles are defined in the observable. The measure + + pt_R + +uses hadron collider coordinates (transverse momenta and boost-invariant +angles). The "energy" in this case is defined as the pT of the jet, +and the "angle" is the distance between the jets in phi, eta space. + +The measure + + E_theta + +uses particle energies and angles and is appropriate for e+e- +collider applications. The "energy" is the jet energy and the angle +between 2 jets is computed from the dot product of the 3 vectors p1 and p2. + +The measure + + E_inv + +uses particle energies and angles and is also appropriate for e+e- +collider applications. In this case “theta” is replaced by Mandelstam +invariants with the same behavior in the collinear limits, leading to a more +calculation friendly observables. The "energy" is defined as the jet energy +and the "angle squared" is defined as (2p_i \cdot p_j/E_i E_j), +where p_i,p_j are the momenta of the jets i adn j, and E_i, E_j are their +respective energies. + +General usage is shown in the example.cc program, and recommended usage +is shown in example_basic_usage.cc. + diff --git a/src/Tools/fjcontrib/EnergyCorrelator/VERSION b/src/Tools/fjcontrib/EnergyCorrelator/VERSION new file mode 100644 --- /dev/null +++ b/src/Tools/fjcontrib/EnergyCorrelator/VERSION @@ -0,0 +1,1 @@ +1.3.1 \ No newline at end of file diff --git a/src/Tools/fjcontrib/EnergyCorrelator/example.cc b/src/Tools/fjcontrib/EnergyCorrelator/example.cc new file mode 100644 --- /dev/null +++ b/src/Tools/fjcontrib/EnergyCorrelator/example.cc @@ -0,0 +1,686 @@ +// Example showing usage of energy correlator classes. +// +// Compile it with "make example" and run it with +// +// ./example < ../data/single-event.dat +// +// Copyright (c) 2013-2016 +// Andrew Larkoski, Lina Necib, Gavin Salam, and Jesse Thaler +// +// $Id: example.cc 1097 2018-01-05 00:04:20Z linoush $ +//---------------------------------------------------------------------- +// 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 . +//---------------------------------------------------------------------- + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "fastjet/PseudoJet.hh" +#include "fastjet/ClusterSequence.hh" +#include "fastjet/JetDefinition.hh" + +#include +#include "EnergyCorrelator.hh" // In external code, this should be fastjet/contrib/EnergyCorrelator.hh + +using namespace std; +using namespace fastjet; +using namespace fastjet::contrib; + +// forward declaration to make things clearer +void read_event(vector &event); +void analyze(const vector & input_particles); + +//---------------------------------------------------------------------- +int main(){ + + //---------------------------------------------------------- + // read in input particles + vector event; + read_event(event); + cout << "# read an event with " << event.size() << " particles" << endl; + + //---------------------------------------------------------- + // illustrate how this EnergyCorrelator contrib works + + analyze(event); + + return 0; +} + +// read in input particles +void read_event(vector &event){ + string line; + while (getline(cin, line)) { + istringstream linestream(line); + // take substrings to avoid problems when there are extra "pollution" + // characters (e.g. line-feed). + if (line.substr(0,4) == "#END") {return;} + if (line.substr(0,1) == "#") {continue;} + double px,py,pz,E; + linestream >> px >> py >> pz >> E; + PseudoJet particle(px,py,pz,E); + + // push event onto back of full_event vector + event.push_back(particle); + } +} + +//////// +// +// Main Routine for Analysis +// +/////// + +void analyze(const vector & input_particles) { + + /////// EnergyCorrelator ///////////////////////////// + + // Initial clustering with anti-kt algorithm + JetAlgorithm algorithm = antikt_algorithm; + double jet_rad = 1.00; // jet radius for anti-kt algorithm + JetDefinition jetDef = JetDefinition(algorithm,jet_rad,E_scheme,Best); + ClusterSequence clust_seq(input_particles,jetDef); + vector antikt_jets = sorted_by_pt(clust_seq.inclusive_jets()); + + for (int j = 0; j < 2; j++) { // Two hardest jets per event + if (antikt_jets[j].perp() > 200) { + + PseudoJet myJet = antikt_jets[j]; + + // various values of beta + vector betalist; + betalist.push_back(0.1); + betalist.push_back(0.2); + betalist.push_back(0.5); + betalist.push_back(1.0); + betalist.push_back(1.5); + betalist.push_back(2.0); + + // various values of alpha + vector alphalist; + alphalist.push_back(0.1); + alphalist.push_back(0.2); + alphalist.push_back(0.5); + alphalist.push_back(1.0); + + + // checking the two energy/angle modes + vector measurelist; + measurelist.push_back(EnergyCorrelator::pt_R); + measurelist.push_back(EnergyCorrelator::E_theta); + //measurelist.push_back(EnergyCorrelator::E_inv); + + vector modename; + modename.push_back("pt_R"); + modename.push_back("E_theta"); + //modename.push_back("E_inv"); + + for (unsigned int M = 0; M < measurelist.size(); M++) { + + cout << "-------------------------------------------------------------------------------------" << endl; + cout << "EnergyCorrelator: ECF(N,beta) with " << modename[M] << endl; + cout << "-------------------------------------------------------------------------------------" << endl; + printf("%7s %14s %14s %14s %14s %15s\n","beta", "N=1 (GeV)", "N=2 (GeV^2)", "N=3 (GeV^3)", "N=4 (GeV^4)", "N=5 (GeV^5)"); + + for (unsigned int B = 0; B < betalist.size(); B++) { + double beta = betalist[B]; + + EnergyCorrelator ECF0(0,beta,measurelist[M]); + EnergyCorrelator ECF1(1,beta,measurelist[M]); + EnergyCorrelator ECF2(2,beta,measurelist[M]); + EnergyCorrelator ECF3(3,beta,measurelist[M]); + EnergyCorrelator ECF4(4,beta,measurelist[M]); + EnergyCorrelator ECF5(5,beta,measurelist[M]); + + printf("%7.3f %14.2f %14.2f %14.2f %14.2f %15.2f \n",beta,ECF1(myJet),ECF2(myJet),ECF3(myJet),ECF4(myJet),ECF5(myJet)); + } + cout << "-------------------------------------------------------------------------------------" << endl << endl; + + cout << "-------------------------------------------------------------------------------------" << endl; + cout << "EnergyCorrelatorRatio: r_N^(beta) = ECF(N+1,beta)/ECF(N,beta) with " << modename[M] << endl; + cout << "-------------------------------------------------------------------------------------" << endl; + printf("%7s %14s %14s %14s %14s %15s \n","beta", "N=0 (GeV)", "N=1 (GeV)", "N=2 (GeV)", "N=3 (GeV)","N=4 (GeV)"); + + for (unsigned int B = 0; B < betalist.size(); B++) { + double beta = betalist[B]; + + EnergyCorrelatorRatio r0(0,beta,measurelist[M]); + EnergyCorrelatorRatio r1(1,beta,measurelist[M]); + EnergyCorrelatorRatio r2(2,beta,measurelist[M]); + EnergyCorrelatorRatio r3(3,beta,measurelist[M]); + EnergyCorrelatorRatio r4(4,beta,measurelist[M]); + + printf("%7.3f %14.4f %14.4f %14.4f %14.4f %15.4f \n",beta,r0(myJet),r1(myJet),r2(myJet),r3(myJet),r4(myJet)); + } + cout << "-------------------------------------------------------------------------------------" << endl << endl; + + cout << "-------------------------------------------------------------------------------------" << endl; + cout << "EnergyCorrelatorDoubleRatio: C_N^(beta) = r_N^(beta)/r_{N-1}^(beta) with " << modename[M] << endl; + cout << "-------------------------------------------------------------------------------------" << endl; + printf("%7s %14s %14s %14s %14s \n","beta", "N=1", "N=2", "N=3", "N=4"); + + for (unsigned int B = 0; B < betalist.size(); B++) { + double beta = betalist[B]; + + EnergyCorrelatorDoubleRatio C1(1,beta,measurelist[M]); + EnergyCorrelatorDoubleRatio C2(2,beta,measurelist[M]); + EnergyCorrelatorDoubleRatio C3(3,beta,measurelist[M]); + EnergyCorrelatorDoubleRatio C4(4,beta,measurelist[M]); + + printf("%7.3f %14.6f %14.6f %14.6f %14.6f \n",beta,C1(myJet),C2(myJet),C3(myJet),C4(myJet)); + } + cout << "-------------------------------------------------------------------------------------" << endl << endl; + + cout << "-------------------------------------------------------------------------------------" << endl; + cout << "EnergyCorrelatorC1: C_1^(beta) = ECF(2,beta)/ECF(1,beta)^2 with " << modename[M] << endl; + cout << "-------------------------------------------------------------------------------------" << endl; + printf("%7s %14s \n","beta","C1 obs"); + + for (unsigned int B = 0; B < betalist.size(); B++) { + double beta = betalist[B]; + + EnergyCorrelatorC1 c1(beta,measurelist[M]); + + printf("%7.3f %14.6f \n",beta,c1(myJet)); + } + cout << "-------------------------------------------------------------------------------------" << endl << endl; + + cout << "-------------------------------------------------------------------------------------" << endl; + cout << "EnergyCorrelatorC2: C_2^(beta) = ECF(3,beta)*ECF(1,beta)/ECF(2,beta)^2 with " << modename[M] << endl; + cout << "-------------------------------------------------------------------------------------" << endl; + printf("%7s %14s \n","beta","C2 obs"); + + for (unsigned int B = 0; B < betalist.size(); B++) { + double beta = betalist[B]; + + EnergyCorrelatorC2 c2(beta,measurelist[M]); + + printf("%7.3f %14.6f \n",beta,c2(myJet)); + } + cout << "-------------------------------------------------------------------------------------" << endl << endl; + + + cout << "-------------------------------------------------------------------------------------" << endl; + cout << "EnergyCorrelatorD2: D_2^(beta) = ECF(3,beta)*ECF(1,beta)^3/ECF(2,beta)^3 with " << modename[M] << endl; + cout << "-------------------------------------------------------------------------------------" << endl; + printf("%7s %14s \n","beta","D2 obs"); + + for (unsigned int B = 0; B < betalist.size(); B++) { + double beta = betalist[B]; + + EnergyCorrelatorD2 d2(beta,measurelist[M]); + + printf("%7.3f %14.6f \n",beta,d2(myJet)); + } + cout << "-------------------------------------------------------------------------------------" << endl << endl; + + cout << "-------------------------------------------------------------------------------------" << endl; + cout << "EnergyCorrelatorGeneralizedD2: D_2^(alpha, beta) = ECFN(3,alpha)/ECFN(2,beta)^(3*alpha/beta) with " << modename[M] << endl; + cout << "-------------------------------------------------------------------------------------" << endl; + printf("%7s %18s %18s %18s %18s\n","beta","alpha = 0.100","alpha = 0.200","alpha = 0.500","alpha = 1.000"); + + for (unsigned int B = 1; B < betalist.size(); B++) { + double beta = betalist[B]; + + printf("%7.3f ", beta); + for (unsigned int A = 0; A < alphalist.size(); A++) { + double alpha = alphalist[A]; + + EnergyCorrelatorGeneralizedD2 d2(alpha, beta, measurelist[M]); + + printf("%18.4f ", d2(myJet)); + } + printf("\n"); + } + cout << "-------------------------------------------------------------------------------------" << endl << endl; + + + + cout << "-------------------------------------------------------------------------------------" << endl; + cout << "EnergyCorrelatorGeneralized (angles = N Choose 2): ECFN(N, beta) with " << modename[M] << endl; + cout << "-------------------------------------------------------------------------------------" << endl; + printf("%7s %7s %14s %14s %14s\n","beta", "N=1", "N=2", "N=3", "N=4"); + + + for (unsigned int B = 0; B < betalist.size(); B++) { + double beta = betalist[B]; + + EnergyCorrelatorGeneralized ECF1(-1,1, beta, measurelist[M]); + EnergyCorrelatorGeneralized ECF2(-1,2, beta, measurelist[M]); + EnergyCorrelatorGeneralized ECF3(-1,3, beta, measurelist[M]); + EnergyCorrelatorGeneralized ECF4(-1,4, beta, measurelist[M]); + //EnergyCorrelatorGeneralized ECF5(-1, 5, beta, measurelist[M]); + + printf("%7.3f %7.2f %14.10f %14.10f %14.10f \n", beta, ECF1(myJet), ECF2(myJet), ECF3(myJet), + ECF4(myJet)); + } + cout << "-------------------------------------------------------------------------------------" << + endl << endl; + + cout << "-------------------------------------------------------------------------------------" << endl; + cout << "EnergyCorrelatorGeneralized: ECFG(angles, N, beta=1) with " << modename[M] << endl; + cout << "-------------------------------------------------------------------------------------" << endl; + printf("%7s %7s %14s %14s %14s\n","angles", "N=1", "N=2", "N=3", "N=4"); + + double beta = 1.0; + for (unsigned int A = 1; A < 2; A++) { + double angle = A; + + EnergyCorrelatorGeneralized ECF1(angle, 1, beta, measurelist[M]); + EnergyCorrelatorGeneralized ECF2(angle, 2, beta, measurelist[M]); + EnergyCorrelatorGeneralized ECF3(angle, 3, beta, measurelist[M]); + EnergyCorrelatorGeneralized ECF4(angle, 4, beta, measurelist[M], EnergyCorrelator::slow); + + printf("%7.0f %7.2f %14.10f %14.10f %14.10f \n", angle, ECF1(myJet), ECF2(myJet), ECF3(myJet), + ECF4(myJet)); + + } + + for (unsigned int A = 2; A < 4; A++) { + double angle = A; + + EnergyCorrelatorGeneralized ECF3(angle, 3, beta, measurelist[M]); + EnergyCorrelatorGeneralized ECF4(angle, 4, beta, measurelist[M]); + + printf("%7.0f %7s %14s %14.10f %14.10f \n", angle, " " , " " ,ECF3(myJet), ECF4(myJet)); + } + + for (unsigned int A = 4; A < 7; A++) { + double angle = A; + + EnergyCorrelatorGeneralized ECF4(angle, 4, beta, measurelist[M]); + printf("%7.0f %7s %14s %14s %14.10f \n", angle, " ", " ", " ", ECF4(myJet) ); + } + cout << "-------------------------------------------------------------------------------------" << + endl << endl; + + cout << "-------------------------------------------------------------------------------------" << endl; + cout << "EnergyCorrelatorNseries: N_i(beta) = ECFG(i+1, 2, beta)/ECFG(i, 1, beta)^2 with " << modename[M] << endl; + cout << "-------------------------------------------------------------------------------------" << endl; + printf("%7s %14s %14s %14s \n","beta", "N=1", "N=2", "N=3"); + + for (unsigned int B = 0; B < betalist.size(); B++) { + double beta = betalist[B]; + + EnergyCorrelatorNseries N1(1,beta,measurelist[M]); + EnergyCorrelatorNseries N2(2,beta,measurelist[M]); + EnergyCorrelatorNseries N3(3,beta,measurelist[M]); + + printf("%7.3f %14.6f %14.6f %14.6f \n",beta,N1(myJet),N2(myJet),N3(myJet)); + + } + cout << "-------------------------------------------------------------------------------------" << endl << endl; + + + cout << "-------------------------------------------------------------------------------------" << endl; + cout << "EnergyCorrelatorN2: N2(beta) = ECFG(3, 2, beta)/ECFG(2, 1, beta)^2 with " << modename[M] << endl; + cout << "-------------------------------------------------------------------------------------" << endl; + printf("%7s %14s \n","beta", "N2 obs"); + + for (unsigned int B = 0; B < betalist.size(); B++) { + double beta = betalist[B]; + + EnergyCorrelatorN2 N2(beta,measurelist[M]); + + printf("%7.3f %14.6f \n",beta,N2(myJet)); + } + cout << "-------------------------------------------------------------------------------------" << endl << endl; + + + cout << "-------------------------------------------------------------------------------------" << endl; + cout << "EnergyCorrelatorN3: N3(beta) = ECFG(4, 2, beta)/ECFG(3, 1, beta)^2 with " << modename[M] << endl; + cout << "-------------------------------------------------------------------------------------" << endl; + printf("%7s %14s \n","beta", "N3 obs"); + + for (unsigned int B = 0; B < betalist.size(); B++) { + double beta = betalist[B]; + + EnergyCorrelatorN3 N3(beta,measurelist[M]); + + printf("%7.3f %14.6f \n",beta,N3(myJet)); + } + cout << "-------------------------------------------------------------------------------------" << endl << endl; + + cout << "-------------------------------------------------------------------------------------" << endl; + cout << "EnergyCorrelatorMseries: M_i(beta) = ECFG(i+1, 1, beta)/ECFN(i, 1, beta) with " << modename[M] << endl; + cout << "-------------------------------------------------------------------------------------" << endl; + printf("%7s %14s %14s %14s \n","beta", "N=1", "N=2", "N=3"); + + for (unsigned int B = 0; B < betalist.size(); B++) { + double beta = betalist[B]; + + EnergyCorrelatorMseries M1(1,beta,measurelist[M]); + EnergyCorrelatorMseries M2(2,beta,measurelist[M]); + EnergyCorrelatorMseries M3(3,beta,measurelist[M]); + + + printf("%7.3f %14.6f %14.6f %14.6f \n",beta,M1(myJet),M2(myJet),M3(myJet)); + } + cout << "-------------------------------------------------------------------------------------" << endl << endl; + + + cout << "-------------------------------------------------------------------------------------" << endl; + cout << "EnergyCorrelatorM2: M2(beta) = ECFG(3, 1, beta)/ECFG(3, 1, beta) with " << modename[M] << endl; + cout << "-------------------------------------------------------------------------------------" << endl; + printf("%7s %14s \n","beta", "M2 obs"); + + for (unsigned int B = 0; B < betalist.size(); B++) { + double beta = betalist[B]; + + EnergyCorrelatorM2 M2(beta,measurelist[M]); + + printf("%7.3f %14.6f \n",beta,M2(myJet)); + } + cout << "-------------------------------------------------------------------------------------" << endl << endl; + + cout << "-------------------------------------------------------------------------------------" << endl; + cout << "EnergyCorrelatorCseries: C_i(beta) = ECFN(i-1, beta)*ECFN(i+1, beta)/ECFN(i, beta)^2 with " << modename[M] << endl; + cout << "-------------------------------------------------------------------------------------" << endl; + printf("%7s %20s %20s %20s \n","beta", "N=1", "N=2", "N=3"); + + for (unsigned int B = 0; B < betalist.size(); B++) { + double beta = betalist[B]; + + EnergyCorrelatorCseries C1(1,beta,measurelist[M]); + EnergyCorrelatorCseries C2(2,beta,measurelist[M]); + EnergyCorrelatorCseries C3(3,beta,measurelist[M]); + + + printf("%7.3f %20.10f %20.10f %20.10f \n",beta,C1(myJet),C2(myJet),C3(myJet)); + } + cout << "-------------------------------------------------------------------------------------" << endl << endl; + + + cout << "-------------------------------------------------------------------------------------" << endl; + cout << "EnergyCorrelatorUseries: U_i(beta) = ECFG(i+1, 1, beta) with " << modename[M] << endl; + cout << "-------------------------------------------------------------------------------------" << endl; + printf("%7s %20s %20s %20s \n","beta", "N=1", "N=2", "N=3"); + + for (unsigned int B = 0; B < betalist.size(); B++) { + double beta = betalist[B]; + + EnergyCorrelatorUseries U1(1,beta,measurelist[M]); + EnergyCorrelatorUseries U2(2,beta,measurelist[M]); + EnergyCorrelatorUseries U3(3,beta,measurelist[M]); + + + printf("%7.3f %20.10f %20.10f %20.10f \n",beta,U1(myJet),U2(myJet),U3(myJet)); + } + cout << "-------------------------------------------------------------------------------------" << endl << endl; + + cout << "-------------------------------------------------------------------------------------" << endl; + cout << "EnergyCorrelatorU1: U1(beta) = ECFG(2, 1, beta) with " << modename[M] << endl; + cout << "-------------------------------------------------------------------------------------" << endl; + printf("%7s %14s \n","beta", "U1 obs"); + + for (unsigned int B = 0; B < betalist.size(); B++) { + double beta = betalist[B]; + + EnergyCorrelatorU1 U1(beta,measurelist[M]); + + printf("%7.3f %14.10f \n",beta,U1(myJet)); + } + cout << "-------------------------------------------------------------------------------------" << endl << endl; + + cout << "-------------------------------------------------------------------------------------" << endl; + cout << "EnergyCorrelatorU2: U2(beta) = ECFG(3, 1, beta) with " << modename[M] << endl; + cout << "-------------------------------------------------------------------------------------" << endl; + printf("%7s %14s \n","beta", "U2 obs"); + + for (unsigned int B = 0; B < betalist.size(); B++) { + double beta = betalist[B]; + + EnergyCorrelatorU2 U2(beta,measurelist[M]); + + printf("%7.3f %14.10f \n",beta,U2(myJet)); + } + cout << "-------------------------------------------------------------------------------------" << endl << endl; + + cout << "-------------------------------------------------------------------------------------" << endl; + cout << "EnergyCorrelatorU3: U3(beta) = ECFG(4, 1, beta) with " << modename[M] << endl; + cout << "-------------------------------------------------------------------------------------" << endl; + printf("%7s %14s \n","beta", "U3 obs"); + + for (unsigned int B = 0; B < betalist.size(); B++) { + double beta = betalist[B]; + + EnergyCorrelatorU3 U3(beta,measurelist[M]); + + printf("%7.3f %14.10f \n",beta,U3(myJet)); + } + cout << "-------------------------------------------------------------------------------------" << endl << endl; + + + // timing tests for the developers + double do_timing_test = false; + if (do_timing_test) { + + cout << "jet with pt = " << myJet.pt() << " and " << myJet.constituents().size() << " constituents" << endl; + + clock_t clock_begin, clock_end; + double num_iter; + double beta = 0.5; + + cout << setprecision(6); + + // test C1 + num_iter = 20000; + clock_begin = clock(); + EnergyCorrelatorDoubleRatio C1s(1,beta,measurelist[M],EnergyCorrelator::slow); + EnergyCorrelatorDoubleRatio C1f(1,beta,measurelist[M],EnergyCorrelator::storage_array); + cout << "timing " << C1s.description() << endl; + cout << "timing " << C1f.description() << endl; + for (int t = 0; t < num_iter; t++) { + C1s(myJet); + } + clock_end = clock(); + cout << "Slow method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per C1"<< endl; + + num_iter = 20000; + clock_begin = clock(); + for (int t = 0; t < num_iter; t++) { + C1f(myJet); + } + clock_end = clock(); + cout << "Storage array method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per C1"<< endl; + + + // test C2 + num_iter = 1000; + clock_begin = clock(); + for (int t = 0; t < num_iter; t++) { + EnergyCorrelatorDoubleRatio C2(2,beta,measurelist[M],EnergyCorrelator::slow); + C2(myJet); + } + clock_end = clock(); + cout << "Slow method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per C2"<< endl; + + num_iter = 10000; + clock_begin = clock(); + for (int t = 0; t < num_iter; t++) { + EnergyCorrelatorDoubleRatio C2(2,beta,measurelist[M],EnergyCorrelator::storage_array); + C2(myJet); + } + clock_end = clock(); + cout << "Storage array method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per C2"<< endl; + + // test C3 + num_iter = 100; + clock_begin = clock(); + + for (int t = 0; t < num_iter; t++) { + EnergyCorrelatorDoubleRatio C3(3,beta,measurelist[M],EnergyCorrelator::slow); + C3(myJet); + } + clock_end = clock(); + cout << "Slow method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per C3"<< endl; + + num_iter = 3000; + clock_begin = clock(); + for (int t = 0; t < num_iter; t++) { + EnergyCorrelatorDoubleRatio C3(3,beta,measurelist[M],EnergyCorrelator::storage_array); + C3(myJet); + } + clock_end = clock(); + cout << "Storage array method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per C3"<< endl; + + // test C4 + num_iter = 10; + clock_begin = clock(); + + for (int t = 0; t < num_iter; t++) { + EnergyCorrelatorDoubleRatio C4(4,beta,measurelist[M],EnergyCorrelator::slow); + C4(myJet); + } + clock_end = clock(); + cout << "Slow method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per C4"<< endl; + + num_iter = 300; + clock_begin = clock(); + for (int t = 0; t < num_iter; t++) { + EnergyCorrelatorDoubleRatio C4(4,beta,measurelist[M],EnergyCorrelator::storage_array); + C4(myJet); + } + clock_end = clock(); + cout << "Storage array method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per C4"<< endl; + + // test N2 + num_iter = 10; + clock_begin = clock(); + + num_iter = 300; + clock_begin = clock(); + for (int t = 0; t < num_iter; t++) { + EnergyCorrelatorN2 N2(beta,measurelist[M],EnergyCorrelator::storage_array); + N2(myJet); + } + clock_end = clock(); + EnergyCorrelatorN2 N2test(beta,measurelist[M],EnergyCorrelator::storage_array); + cout << "Beta is: "<< beta << endl; + cout << "Result of N2: "<< N2test(myJet) << endl; + cout << "Storage array method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per N2"<< endl; + + + num_iter = 300; + clock_begin = clock(); + for (int t = 0; t < num_iter; t++) { + EnergyCorrelatorN3 N3(beta,measurelist[M],EnergyCorrelator::storage_array); + N3(myJet); + } + clock_end = clock(); + EnergyCorrelatorN3 N3test(beta,measurelist[M],EnergyCorrelator::storage_array); + cout << "Beta is: "<< beta << endl; + cout << "Result of N3: "<< N3test(myJet) << endl; + cout << "Storage array method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per N3"<< endl; + + + + + num_iter = 300; + clock_begin = clock(); + for (int t = 0; t < num_iter; t++) { + EnergyCorrelatorGeneralized ECF1(2,3, beta, measurelist[M]); + ECF1(myJet); + } + clock_end = clock(); + EnergyCorrelatorGeneralized ECF1test(2,3, beta, measurelist[M]); + cout << "Beta is: "<< beta << endl; + cout << "Result of 2e3: "<< ECF1test(myJet) << endl; + cout << "Storage array method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per 2e3"<< endl; + + + num_iter = 300; + clock_begin = clock(); + for (int t = 0; t < num_iter; t++) { + EnergyCorrelatorGeneralized ECF3(2,4, beta, measurelist[M]); + ECF3(myJet); + } + clock_end = clock(); + EnergyCorrelatorGeneralized ECF2test(2,4, beta, measurelist[M]); + cout << "Beta is: "<< beta << endl; + cout << "Result of 2e4: "<< ECF2test(myJet) << endl; + cout << "Storage array method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per 2e4"<< endl; + + +// num_iter = 300; +// clock_begin = clock(); +// for (int t = 0; t < num_iter; t++) { +// EnergyCorrelatorGeneralized ECF5(2,5, beta, measurelist[M]); +// ECF5(myJet); +// } +// clock_end = clock(); +// EnergyCorrelatorGeneralized ECF5test(2,5, beta, measurelist[M]); +// cout << "Beta is: "<< beta << endl; +// cout << "Result of 2e5: "<< ECF5test(myJet) << endl; +// cout << "Storage array method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per 2e5"<< endl; +// + + // test M2 + num_iter = 10; + clock_begin = clock(); + + for (int t = 0; t < num_iter; t++) { + EnergyCorrelatorM2 M2(beta,measurelist[M],EnergyCorrelator::slow); + M2(myJet); + } + clock_end = clock(); + cout << "Slow method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per M2"<< endl; + + num_iter = 300; + clock_begin = clock(); + for (int t = 0; t < num_iter; t++) { + EnergyCorrelatorM2 M2(beta,measurelist[M],EnergyCorrelator::storage_array); + M2(myJet); + } + clock_end = clock(); + cout << "Storage array method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per M2"<< endl; + + // test M3 + num_iter = 10; + clock_begin = clock(); + + for (int t = 0; t < num_iter; t++) { + EnergyCorrelatorMseries M3(3,beta,measurelist[M],EnergyCorrelator::slow); + M3(myJet); + } + clock_end = clock(); + cout << "Slow method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per M3"<< endl; + + num_iter = 300; + clock_begin = clock(); + for (int t = 0; t < num_iter; t++) { + EnergyCorrelatorMseries M3(3,beta,measurelist[M],EnergyCorrelator::storage_array); + M3(myJet); + } + clock_end = clock(); + cout << "Storage array method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per M3"<< endl; + + + } + } + } + } +} + + + diff --git a/src/Tools/fjcontrib/EnergyCorrelator/example.ref b/src/Tools/fjcontrib/EnergyCorrelator/example.ref new file mode 100644 --- /dev/null +++ b/src/Tools/fjcontrib/EnergyCorrelator/example.ref @@ -0,0 +1,922 @@ +# read an event with 354 particles +#-------------------------------------------------------------------------- +# FastJet release 3.2.0 +# M. Cacciari, G.P. Salam and G. Soyez +# A software package for jet finding and analysis at colliders +# http://fastjet.fr +# +# Please cite EPJC72(2012)1896 [arXiv:1111.6097] if you use this package +# for scientific work and optionally PLB641(2006)57 [hep-ph/0512210]. +# +# FastJet is provided without warranty under the terms of the GNU GPLv2. +# It uses T. Chan's closest pair algorithm, S. Fortune's Voronoi code +# and 3rd party plugin jet algorithms. See COPYING file for details. +#-------------------------------------------------------------------------- +------------------------------------------------------------------------------------- +EnergyCorrelator: ECF(N,beta) with pt_R +------------------------------------------------------------------------------------- + beta N=1 (GeV) N=2 (GeV^2) N=3 (GeV^3) N=4 (GeV^4) N=5 (GeV^5) + 0.100 983.64 265690.56 27540663.76 1287213926.38 30113720909.66 + 0.200 983.64 172787.11 8222752.45 134175168.40 888560847.26 + 0.500 983.64 52039.98 364021.49 608244.06 356259.05 + 1.000 983.64 10006.49 9934.06 2670.10 271.98 + 1.500 983.64 3001.20 1066.15 118.17 4.24 + 2.000 983.64 1272.64 260.18 14.88 0.22 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorRatio: r_N^(beta) = ECF(N+1,beta)/ECF(N,beta) with pt_R +------------------------------------------------------------------------------------- + beta N=0 (GeV) N=1 (GeV) N=2 (GeV) N=3 (GeV) N=4 (GeV) + 0.100 983.6369 270.1104 103.6569 46.7387 23.3945 + 0.200 983.6369 175.6615 47.5889 16.3175 6.6224 + 0.500 983.6369 52.9057 6.9950 1.6709 0.5857 + 1.000 983.6369 10.1730 0.9928 0.2688 0.1019 + 1.500 983.6369 3.0511 0.3552 0.1108 0.0359 + 2.000 983.6369 1.2938 0.2044 0.0572 0.0151 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorDoubleRatio: C_N^(beta) = r_N^(beta)/r_{N-1}^(beta) with pt_R +------------------------------------------------------------------------------------- + beta N=1 N=2 N=3 N=4 + 0.100 0.274604 0.383758 0.450898 0.500538 + 0.200 0.178584 0.270913 0.342885 0.405845 + 0.500 0.053786 0.132217 0.238870 0.350540 + 1.000 0.010342 0.097588 0.270742 0.378972 + 1.500 0.003102 0.116430 0.312016 0.323906 + 2.000 0.001315 0.158011 0.279675 0.263310 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorC1: C_1^(beta) = ECF(2,beta)/ECF(1,beta)^2 with pt_R +------------------------------------------------------------------------------------- + beta C1 obs + 0.100 0.274604 + 0.200 0.178584 + 0.500 0.053786 + 1.000 0.010342 + 1.500 0.003102 + 2.000 0.001315 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorC2: C_2^(beta) = ECF(3,beta)*ECF(1,beta)/ECF(2,beta)^2 with pt_R +------------------------------------------------------------------------------------- + beta C2 obs + 0.100 0.383758 + 0.200 0.270913 + 0.500 0.132217 + 1.000 0.097588 + 1.500 0.116430 + 2.000 0.158011 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorD2: D_2^(beta) = ECF(3,beta)*ECF(1,beta)^3/ECF(2,beta)^3 with pt_R +------------------------------------------------------------------------------------- + beta D2 obs + 0.100 1.397496 + 0.200 1.517007 + 0.500 2.458216 + 1.000 9.435950 + 1.500 37.535182 + 2.000 120.129760 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorGeneralizedD2: D_2^(alpha, beta) = ECFN(3,alpha)/ECFN(2,beta)^(3*alpha/beta) with pt_R +------------------------------------------------------------------------------------- + beta alpha = 0.100 alpha = 0.200 alpha = 0.500 alpha = 1.000 + 0.200 0.3834 1.5170 156.2463 1741793.7615 + 0.500 0.1671 0.2882 2.4582 431.1389 + 1.000 0.1140 0.1342 0.3637 9.4359 + 1.500 0.0919 0.0871 0.1233 1.0849 + 2.000 0.0783 0.0632 0.0554 0.2188 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorGeneralized (angles = N Choose 2): ECFN(N, beta) with pt_R +------------------------------------------------------------------------------------- + beta N=1 N=2 N=3 N=4 + 0.100 1.00 0.2746037664 0.0289380967 0.0013750278 + 0.200 1.00 0.1785836568 0.0086399808 0.0001433286 + 0.500 1.00 0.0537857869 0.0003824922 0.0000006497 + 1.000 1.00 0.0103421840 0.0000104381 0.0000000029 + 1.500 1.00 0.0031018827 0.0000011202 0.0000000001 + 2.000 1.00 0.0013153375 0.0000002734 0.0000000000 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorGeneralized: ECFG(angles, N, beta=1) with pt_R +------------------------------------------------------------------------------------- + angles N=1 N=2 N=3 N=4 + 1 1.00 0.0103421840 0.0007353402 0.0000733472 + 2 0.0000460939 0.0000010704 + 3 0.0000104381 0.0000000552 + 4 0.0000000114 + 5 0.0000000045 + 6 0.0000000029 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorNseries: N_i(beta) = ECFG(i+1, 2, beta)/ECFG(i, 1, beta)^2 with pt_R +------------------------------------------------------------------------------------- + beta N=1 N=2 N=3 + 0.100 0.178584 0.551898 1.506618 + 0.200 0.078954 0.534774 1.521976 + 0.500 0.010342 0.487547 1.601014 + 1.000 0.001315 0.430942 1.979487 + 1.500 0.000437 0.391935 3.250084 + 2.000 0.000220 0.365859 6.700247 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorN2: N2(beta) = ECFG(3, 2, beta)/ECFG(2, 1, beta)^2 with pt_R +------------------------------------------------------------------------------------- + beta N2 obs + 0.100 0.551898 + 0.200 0.534774 + 0.500 0.487547 + 1.000 0.430942 + 1.500 0.391935 + 2.000 0.365859 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorN3: N3(beta) = ECFG(4, 2, beta)/ECFG(3, 1, beta)^2 with pt_R +------------------------------------------------------------------------------------- + beta N3 obs + 0.100 1.506618 + 0.200 1.521976 + 0.500 1.601014 + 1.000 1.979487 + 1.500 3.250084 + 2.000 6.700247 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorMseries: M_i(beta) = ECFG(i+1, 1, beta)/ECFN(i, 1, beta) with pt_R +------------------------------------------------------------------------------------- + beta N=1 N=2 N=3 + 0.100 0.274604 0.225158 0.151396 + 0.200 0.178584 0.206357 0.146715 + 0.500 0.053786 0.149827 0.130985 + 1.000 0.010342 0.071101 0.099746 + 1.500 0.003102 0.027603 0.065806 + 2.000 0.001315 0.010448 0.036511 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorM2: M2(beta) = ECFG(3, 1, beta)/ECFG(3, 1, beta) with pt_R +------------------------------------------------------------------------------------- + beta M2 obs + 0.100 0.225158 + 0.200 0.206357 + 0.500 0.149827 + 1.000 0.071101 + 1.500 0.027603 + 2.000 0.010448 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorCseries: C_i(beta) = ECFN(i-1, beta)*ECFN(i+1, beta)/ECFN(i, beta)^2 with pt_R +------------------------------------------------------------------------------------- + beta N=1 N=2 N=3 + 0.100 0.2746037664 0.3837575950 0.4508977207 + 0.200 0.1785836568 0.2709126926 0.3428854490 + 0.500 0.0537857869 0.1322170711 0.2388696533 + 1.000 0.0103421840 0.0975883303 0.2707417692 + 1.500 0.0031018827 0.1164297328 0.3120164994 + 2.000 0.0013153375 0.1580111821 0.2796746959 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorUseries: U_i(beta) = ECFG(i+1, 1, beta) with pt_R +------------------------------------------------------------------------------------- + beta N=1 N=2 N=3 + 0.100 0.2746037664 0.0618292258 0.0093606928 + 0.200 0.1785836568 0.0368519186 0.0054067344 + 0.500 0.0537857869 0.0080585540 0.0010555469 + 1.000 0.0103421840 0.0007353402 0.0000733472 + 1.500 0.0031018827 0.0000856204 0.0000056343 + 2.000 0.0013153375 0.0000137425 0.0000005017 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorU1: U1(beta) = ECFG(2, 1, beta) with pt_R +------------------------------------------------------------------------------------- + beta U1 obs + 0.100 0.2746037664 + 0.200 0.1785836568 + 0.500 0.0537857869 + 1.000 0.0103421840 + 1.500 0.0031018827 + 2.000 0.0013153375 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorU2: U2(beta) = ECFG(3, 1, beta) with pt_R +------------------------------------------------------------------------------------- + beta U2 obs + 0.100 0.0618292258 + 0.200 0.0368519186 + 0.500 0.0080585540 + 1.000 0.0007353402 + 1.500 0.0000856204 + 2.000 0.0000137425 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorU3: U3(beta) = ECFG(4, 1, beta) with pt_R +------------------------------------------------------------------------------------- + beta U3 obs + 0.100 0.0093606928 + 0.200 0.0054067344 + 0.500 0.0010555469 + 1.000 0.0000733472 + 1.500 0.0000056343 + 2.000 0.0000005017 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelator: ECF(N,beta) with E_theta +------------------------------------------------------------------------------------- + beta N=1 (GeV) N=2 (GeV^2) N=3 (GeV^3) N=4 (GeV^4) N=5 (GeV^5) + 0.100 1378.16 504056.94 68362746.10 4040636656.75 115642116882.67 + 0.200 1378.16 316828.68 18441161.70 345185274.86 2475861856.72 + 0.500 1378.16 86243.22 614652.00 941844.70 459093.24 + 1.000 1378.16 14186.79 11682.63 2060.77 92.98 + 1.500 1378.16 3751.90 907.11 39.34 0.27 + 2.000 1378.16 1439.57 152.71 1.84 0.00 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorRatio: r_N^(beta) = ECF(N+1,beta)/ECF(N,beta) with E_theta +------------------------------------------------------------------------------------- + beta N=0 (GeV) N=1 (GeV) N=2 (GeV) N=3 (GeV) N=4 (GeV) + 0.100 1378.1622 365.7457 135.6250 59.1058 28.6198 + 0.200 1378.1622 229.8922 58.2055 18.7182 7.1726 + 0.500 1378.1622 62.5784 7.1270 1.5323 0.4874 + 1.000 1378.1622 10.2940 0.8235 0.1764 0.0451 + 1.500 1378.1622 2.7224 0.2418 0.0434 0.0069 + 2.000 1378.1622 1.0446 0.1061 0.0121 0.0013 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorDoubleRatio: C_N^(beta) = r_N^(beta)/r_{N-1}^(beta) with E_theta +------------------------------------------------------------------------------------- + beta N=1 N=2 N=3 N=4 + 0.100 0.265387 0.370818 0.435803 0.484212 + 0.200 0.166811 0.253186 0.321588 0.383186 + 0.500 0.045407 0.113888 0.215004 0.318106 + 1.000 0.007469 0.079997 0.214207 0.255786 + 1.500 0.001975 0.088810 0.179360 0.160180 + 2.000 0.000758 0.101555 0.113867 0.106401 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorC1: C_1^(beta) = ECF(2,beta)/ECF(1,beta)^2 with E_theta +------------------------------------------------------------------------------------- + beta C1 obs + 0.100 0.265387 + 0.200 0.166811 + 0.500 0.045407 + 1.000 0.007469 + 1.500 0.001975 + 2.000 0.000758 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorC2: C_2^(beta) = ECF(3,beta)*ECF(1,beta)/ECF(2,beta)^2 with E_theta +------------------------------------------------------------------------------------- + beta C2 obs + 0.100 0.370818 + 0.200 0.253186 + 0.500 0.113888 + 1.000 0.079997 + 1.500 0.088810 + 2.000 0.101555 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorD2: D_2^(beta) = ECF(3,beta)*ECF(1,beta)^3/ECF(2,beta)^3 with E_theta +------------------------------------------------------------------------------------- + beta D2 obs + 0.100 1.397274 + 0.200 1.517804 + 0.500 2.508161 + 1.000 10.709981 + 1.500 44.958245 + 2.000 133.988418 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorGeneralizedD2: D_2^(alpha, beta) = ECFN(3,alpha)/ECFN(2,beta)^(3*alpha/beta) with E_theta +------------------------------------------------------------------------------------- + beta alpha = 0.100 alpha = 0.200 alpha = 0.500 alpha = 1.000 + 0.200 0.3833 1.5178 159.9740 2071485.6064 + 0.500 0.1670 0.2880 2.5082 509.2064 + 1.000 0.1135 0.1330 0.3637 10.7100 + 1.500 0.0907 0.0850 0.1189 1.1438 + 2.000 0.0767 0.0608 0.0514 0.2139 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorGeneralized (angles = N Choose 2): ECFN(N, beta) with E_theta +------------------------------------------------------------------------------------- + beta N=1 N=2 N=3 N=4 + 0.100 1.00 0.2653865755 0.0261167137 0.0011200786 + 0.200 1.00 0.1668106769 0.0070451023 0.0000956866 + 0.500 1.00 0.0454071552 0.0002348163 0.0000002611 + 1.000 1.00 0.0074693640 0.0000044631 0.0000000006 + 1.500 1.00 0.0019753776 0.0000003465 0.0000000000 + 2.000 1.00 0.0007579352 0.0000000583 0.0000000000 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorGeneralized: ECFG(angles, N, beta=1) with E_theta +------------------------------------------------------------------------------------- + angles N=1 N=2 N=3 N=4 + 1 1.00 0.0074693640 0.0005219446 0.0000521319 + 2 0.0000242551 0.0000005388 + 3 0.0000044631 0.0000000205 + 4 0.0000000035 + 5 0.0000000011 + 6 0.0000000006 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorNseries: N_i(beta) = ECFG(i+1, 2, beta)/ECFG(i, 1, beta)^2 with E_theta +------------------------------------------------------------------------------------- + beta N=1 N=2 N=3 + 0.100 0.166811 0.551800 1.506192 + 0.200 0.068920 0.534652 1.521253 + 0.500 0.007469 0.487710 1.598907 + 1.000 0.000758 0.434746 1.977674 + 1.500 0.000207 0.401351 3.333270 + 2.000 0.000082 0.372570 7.518334 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorN2: N2(beta) = ECFG(3, 2, beta)/ECFG(2, 1, beta)^2 with E_theta +------------------------------------------------------------------------------------- + beta N2 obs + 0.100 0.551800 + 0.200 0.534652 + 0.500 0.487710 + 1.000 0.434746 + 1.500 0.401351 + 2.000 0.372570 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorN3: N3(beta) = ECFG(4, 2, beta)/ECFG(3, 1, beta)^2 with E_theta +------------------------------------------------------------------------------------- + beta N3 obs + 0.100 1.506192 + 0.200 1.521253 + 0.500 1.598907 + 1.000 1.977674 + 1.500 3.333270 + 2.000 7.518334 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorMseries: M_i(beta) = ECFG(i+1, 1, beta)/ECFN(i, 1, beta) with E_theta +------------------------------------------------------------------------------------- + beta N=1 N=2 N=3 + 0.100 0.265387 0.225134 0.151350 + 0.200 0.166811 0.206335 0.146685 + 0.500 0.045407 0.149607 0.131015 + 1.000 0.007469 0.069878 0.099880 + 1.500 0.001975 0.025930 0.065823 + 2.000 0.000758 0.009248 0.036027 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorM2: M2(beta) = ECFG(3, 1, beta)/ECFG(3, 1, beta) with E_theta +------------------------------------------------------------------------------------- + beta M2 obs + 0.100 0.225134 + 0.200 0.206335 + 0.500 0.149607 + 1.000 0.069878 + 1.500 0.025930 + 2.000 0.009248 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorCseries: C_i(beta) = ECFN(i-1, beta)*ECFN(i+1, beta)/ECFN(i, beta)^2 with E_theta +------------------------------------------------------------------------------------- + beta N=1 N=2 N=3 + 0.100 0.2653865755 0.3708178471 0.4358031826 + 0.200 0.1668106769 0.2531859580 0.3215882713 + 0.500 0.0454071552 0.1138884647 0.2150035566 + 1.000 0.0074693640 0.0799967491 0.2142068062 + 1.500 0.0019753776 0.0888095113 0.1793603915 + 2.000 0.0007579352 0.1015545382 0.1138668084 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorUseries: U_i(beta) = ECFG(i+1, 1, beta) with E_theta +------------------------------------------------------------------------------------- + beta N=1 N=2 N=3 + 0.100 0.2653865755 0.0597476192 0.0090428096 + 0.200 0.1668106769 0.0344188111 0.0050487063 + 0.500 0.0454071552 0.0067932369 0.0008900155 + 1.000 0.0074693640 0.0005219446 0.0000521319 + 1.500 0.0019753776 0.0000512213 0.0000033715 + 2.000 0.0007579352 0.0000070092 0.0000002525 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorU1: U1(beta) = ECFG(2, 1, beta) with E_theta +------------------------------------------------------------------------------------- + beta U1 obs + 0.100 0.2653865755 + 0.200 0.1668106769 + 0.500 0.0454071552 + 1.000 0.0074693640 + 1.500 0.0019753776 + 2.000 0.0007579352 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorU2: U2(beta) = ECFG(3, 1, beta) with E_theta +------------------------------------------------------------------------------------- + beta U2 obs + 0.100 0.0597476192 + 0.200 0.0344188111 + 0.500 0.0067932369 + 1.000 0.0005219446 + 1.500 0.0000512213 + 2.000 0.0000070092 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorU3: U3(beta) = ECFG(4, 1, beta) with E_theta +------------------------------------------------------------------------------------- + beta U3 obs + 0.100 0.0090428096 + 0.200 0.0050487063 + 0.500 0.0008900155 + 1.000 0.0000521319 + 1.500 0.0000033715 + 2.000 0.0000002525 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelator: ECF(N,beta) with pt_R +------------------------------------------------------------------------------------- + beta N=1 (GeV) N=2 (GeV^2) N=3 (GeV^3) N=4 (GeV^4) N=5 (GeV^5) + 0.100 910.03 163661.63 10858223.95 340994615.02 5907941830.38 + 0.200 910.03 111307.49 3950123.18 59638399.31 478368362.42 + 0.500 910.03 40470.38 446453.20 2147680.90 5564873.72 + 1.000 910.03 13863.90 65435.19 142222.94 166409.02 + 1.500 910.03 8417.59 26548.97 39560.12 31484.21 + 2.000 910.03 6328.24 16871.11 21081.97 13743.58 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorRatio: r_N^(beta) = ECF(N+1,beta)/ECF(N,beta) with pt_R +------------------------------------------------------------------------------------- + beta N=0 (GeV) N=1 (GeV) N=2 (GeV) N=3 (GeV) N=4 (GeV) + 0.100 910.0320 179.8416 66.3456 31.4043 17.3256 + 0.200 910.0320 122.3116 35.4884 15.0979 8.0211 + 0.500 910.0320 44.4714 11.0316 4.8105 2.5911 + 1.000 910.0320 15.2345 4.7198 2.1735 1.1701 + 1.500 910.0320 9.2498 3.1540 1.4901 0.7959 + 2.000 910.0320 6.9539 2.6660 1.2496 0.6519 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorDoubleRatio: C_N^(beta) = r_N^(beta)/r_{N-1}^(beta) with pt_R +------------------------------------------------------------------------------------- + beta N=1 N=2 N=3 N=4 + 0.100 0.197621 0.368911 0.473344 0.551696 + 0.200 0.134404 0.290147 0.425431 0.531277 + 0.500 0.048868 0.248061 0.436069 0.538632 + 1.000 0.016741 0.309811 0.460503 0.538330 + 1.500 0.010164 0.340980 0.472444 0.534103 + 2.000 0.007641 0.383385 0.468712 0.521701 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorC1: C_1^(beta) = ECF(2,beta)/ECF(1,beta)^2 with pt_R +------------------------------------------------------------------------------------- + beta C1 obs + 0.100 0.197621 + 0.200 0.134404 + 0.500 0.048868 + 1.000 0.016741 + 1.500 0.010164 + 2.000 0.007641 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorC2: C_2^(beta) = ECF(3,beta)*ECF(1,beta)/ECF(2,beta)^2 with pt_R +------------------------------------------------------------------------------------- + beta C2 obs + 0.100 0.368911 + 0.200 0.290147 + 0.500 0.248061 + 1.000 0.309811 + 1.500 0.340980 + 2.000 0.383385 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorD2: D_2^(beta) = ECF(3,beta)*ECF(1,beta)^3/ECF(2,beta)^3 with pt_R +------------------------------------------------------------------------------------- + beta D2 obs + 0.100 1.866758 + 0.200 2.158775 + 0.500 5.076145 + 1.000 18.506536 + 1.500 33.547049 + 2.000 50.172493 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorGeneralizedD2: D_2^(alpha, beta) = ECFN(3,alpha)/ECFN(2,beta)^(3*alpha/beta) with pt_R +------------------------------------------------------------------------------------- + beta alpha = 0.100 alpha = 0.200 alpha = 0.500 alpha = 1.000 + 0.200 0.2924 2.1588 2039.4962 1029141861.7962 + 0.500 0.0881 0.1962 5.0761 6375.2544 + 1.000 0.0491 0.0610 0.2735 18.5065 + 1.500 0.0361 0.0329 0.0583 0.8404 + 2.000 0.0299 0.0226 0.0229 0.1300 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorGeneralized (angles = N Choose 2): ECFN(N, beta) with pt_R +------------------------------------------------------------------------------------- + beta N=1 N=2 N=3 N=4 + 0.100 1.00 0.1976212201 0.0144075072 0.0004971883 + 0.200 1.00 0.1344036550 0.0052413202 0.0000869560 + 0.500 1.00 0.0488679347 0.0005923876 0.0000031314 + 1.000 1.00 0.0167406417 0.0000868243 0.0000002074 + 1.500 1.00 0.0101642297 0.0000352272 0.0000000577 + 2.000 1.00 0.0076413374 0.0000223859 0.0000000307 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorGeneralized: ECFG(angles, N, beta=1) with pt_R +------------------------------------------------------------------------------------- + angles N=1 N=2 N=3 N=4 + 1 1.00 0.0167406417 0.0005149099 0.0000275956 + 2 0.0001201519 0.0000010521 + 3 0.0000868243 0.0000003239 + 4 0.0000002207 + 5 0.0000001864 + 6 0.0000002074 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorNseries: N_i(beta) = ECFG(i+1, 2, beta)/ECFG(i, 1, beta)^2 with pt_R +------------------------------------------------------------------------------------- + beta N=1 N=2 N=3 + 0.100 0.134404 0.501515 2.093230 + 0.200 0.066590 0.490962 2.107355 + 0.500 0.016741 0.472128 2.270271 + 1.000 0.007641 0.428733 3.968204 + 1.500 0.005197 0.367850 11.456776 + 2.000 0.003938 0.329248 20.513687 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorN2: N2(beta) = ECFG(3, 2, beta)/ECFG(2, 1, beta)^2 with pt_R +------------------------------------------------------------------------------------- + beta N2 obs + 0.100 0.501515 + 0.200 0.490962 + 0.500 0.472128 + 1.000 0.428733 + 1.500 0.367850 + 2.000 0.329248 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorN3: N3(beta) = ECFG(4, 2, beta)/ECFG(3, 1, beta)^2 with pt_R +------------------------------------------------------------------------------------- + beta N3 obs + 0.100 2.093230 + 0.200 2.107355 + 0.500 2.270271 + 1.000 3.968204 + 1.500 11.456776 + 2.000 20.513687 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorMseries: M_i(beta) = ECFG(i+1, 1, beta)/ECFN(i, 1, beta) with pt_R +------------------------------------------------------------------------------------- + beta N=1 N=2 N=3 + 0.100 0.197621 0.140440 0.090587 + 0.200 0.134404 0.127907 0.086996 + 0.500 0.048868 0.087042 0.076041 + 1.000 0.016741 0.030758 0.053593 + 1.500 0.010164 0.010637 0.025871 + 2.000 0.007641 0.005839 0.009492 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorM2: M2(beta) = ECFG(3, 1, beta)/ECFG(3, 1, beta) with pt_R +------------------------------------------------------------------------------------- + beta M2 obs + 0.100 0.140440 + 0.200 0.127907 + 0.500 0.087042 + 1.000 0.030758 + 1.500 0.010637 + 2.000 0.005839 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorCseries: C_i(beta) = ECFN(i-1, beta)*ECFN(i+1, beta)/ECFN(i, beta)^2 with pt_R +------------------------------------------------------------------------------------- + beta N=1 N=2 N=3 + 0.100 0.1976212201 0.3689110739 0.4733439247 + 0.200 0.1344036550 0.2901472976 0.4254309496 + 0.500 0.0488679347 0.2480607340 0.4360689366 + 1.000 0.0167406417 0.3098112813 0.4605028221 + 1.500 0.0101642297 0.3409799143 0.4724436009 + 2.000 0.0076413374 0.3833849510 0.4687122507 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorUseries: U_i(beta) = ECFG(i+1, 1, beta) with pt_R +------------------------------------------------------------------------------------- + beta N=1 N=2 N=3 + 0.100 0.1976212201 0.0277540195 0.0025141656 + 0.200 0.1344036550 0.0171911179 0.0014955586 + 0.500 0.0488679347 0.0042535497 0.0003234433 + 1.000 0.0167406417 0.0005149099 0.0000275956 + 1.500 0.0101642297 0.0001081211 0.0000027972 + 2.000 0.0076413374 0.0000446167 0.0000004235 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorU1: U1(beta) = ECFG(2, 1, beta) with pt_R +------------------------------------------------------------------------------------- + beta U1 obs + 0.100 0.1976212201 + 0.200 0.1344036550 + 0.500 0.0488679347 + 1.000 0.0167406417 + 1.500 0.0101642297 + 2.000 0.0076413374 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorU2: U2(beta) = ECFG(3, 1, beta) with pt_R +------------------------------------------------------------------------------------- + beta U2 obs + 0.100 0.0277540195 + 0.200 0.0171911179 + 0.500 0.0042535497 + 1.000 0.0005149099 + 1.500 0.0001081211 + 2.000 0.0000446167 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorU3: U3(beta) = ECFG(4, 1, beta) with pt_R +------------------------------------------------------------------------------------- + beta U3 obs + 0.100 0.0025141656 + 0.200 0.0014955586 + 0.500 0.0003234433 + 1.000 0.0000275956 + 1.500 0.0000027972 + 2.000 0.0000004235 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelator: ECF(N,beta) with E_theta +------------------------------------------------------------------------------------- + beta N=1 (GeV) N=2 (GeV^2) N=3 (GeV^3) N=4 (GeV^4) N=5 (GeV^5) + 0.100 934.39 173035.93 11869881.64 387687952.45 7030672399.84 + 0.200 934.39 117791.85 4351824.23 69148951.24 587495287.84 + 0.500 934.39 43220.34 507378.92 2589996.65 7002091.93 + 1.000 934.39 15182.02 75656.66 165115.95 178354.07 + 1.500 934.39 9296.03 30143.02 40417.01 24563.71 + 2.000 934.39 6943.30 18259.83 17679.34 7116.38 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorRatio: r_N^(beta) = ECF(N+1,beta)/ECF(N,beta) with E_theta +------------------------------------------------------------------------------------- + beta N=0 (GeV) N=1 (GeV) N=2 (GeV) N=3 (GeV) N=4 (GeV) + 0.100 934.3868 185.1866 68.5978 32.6615 18.1349 + 0.200 934.3868 126.0633 36.9450 15.8896 8.4961 + 0.500 934.3868 46.2553 11.7394 5.1047 2.7035 + 1.000 934.3868 16.2481 4.9833 2.1824 1.0802 + 1.500 934.3868 9.9488 3.2426 1.3408 0.6078 + 2.000 934.3868 7.4309 2.6299 0.9682 0.4025 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorDoubleRatio: C_N^(beta) = r_N^(beta)/r_{N-1}^(beta) with E_theta +------------------------------------------------------------------------------------- + beta N=1 N=2 N=3 N=4 + 0.100 0.198191 0.370425 0.476130 0.555237 + 0.200 0.134916 0.293067 0.430089 0.534693 + 0.500 0.049503 0.253795 0.434833 0.529617 + 1.000 0.017389 0.306701 0.437950 0.494940 + 1.500 0.010647 0.325926 0.413512 0.453265 + 2.000 0.007953 0.353909 0.368161 0.415742 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorC1: C_1^(beta) = ECF(2,beta)/ECF(1,beta)^2 with E_theta +------------------------------------------------------------------------------------- + beta C1 obs + 0.100 0.198191 + 0.200 0.134916 + 0.500 0.049503 + 1.000 0.017389 + 1.500 0.010647 + 2.000 0.007953 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorC2: C_2^(beta) = ECF(3,beta)*ECF(1,beta)/ECF(2,beta)^2 with E_theta +------------------------------------------------------------------------------------- + beta C2 obs + 0.100 0.370425 + 0.200 0.293067 + 0.500 0.253795 + 1.000 0.306701 + 1.500 0.325926 + 2.000 0.353909 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorD2: D_2^(beta) = ECF(3,beta)*ECF(1,beta)^3/ECF(2,beta)^3 with E_theta +------------------------------------------------------------------------------------- + beta D2 obs + 0.100 1.869035 + 0.200 2.172229 + 0.500 5.126818 + 1.000 17.637552 + 1.500 30.610782 + 2.000 44.502015 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorGeneralizedD2: D_2^(alpha, beta) = ECFN(3,alpha)/ECFN(2,beta)^(3*alpha/beta) with E_theta +------------------------------------------------------------------------------------- + beta alpha = 0.100 alpha = 0.200 alpha = 0.500 alpha = 1.000 + 0.200 0.2936 2.1722 2081.0777 1038338617.9694 + 0.500 0.0883 0.1966 5.1268 6301.7080 + 1.000 0.0491 0.0607 0.2712 17.6376 + 1.500 0.0361 0.0328 0.0584 0.8180 + 2.000 0.0300 0.0227 0.0234 0.1308 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorGeneralized (angles = N Choose 2): ECFN(N, beta) with E_theta +------------------------------------------------------------------------------------- + beta N=1 N=2 N=3 N=4 + 0.100 1.00 0.1981905437 0.0145501117 0.0005085991 + 0.200 1.00 0.1349155085 0.0053344701 0.0000907149 + 0.500 1.00 0.0495033749 0.0006219455 0.0000033978 + 1.000 1.00 0.0173890665 0.0000927400 0.0000002166 + 1.500 1.00 0.0106474132 0.0000369493 0.0000000530 + 2.000 1.00 0.0079526591 0.0000223829 0.0000000232 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorGeneralized: ECFG(angles, N, beta=1) with E_theta +------------------------------------------------------------------------------------- + angles N=1 N=2 N=3 N=4 + 1 1.00 0.0173890665 0.0005243822 0.0000276026 + 2 0.0001290539 0.0000011436 + 3 0.0000927400 0.0000003720 + 4 0.0000002479 + 5 0.0000002021 + 6 0.0000002166 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorNseries: N_i(beta) = ECFG(i+1, 2, beta)/ECFG(i, 1, beta)^2 with E_theta +------------------------------------------------------------------------------------- + beta N=1 N=2 N=3 + 0.100 0.134916 0.502292 2.091587 + 0.200 0.067180 0.491827 2.107891 + 0.500 0.017389 0.472830 2.289141 + 1.000 0.007953 0.426794 4.159019 + 1.500 0.005251 0.366194 11.630367 + 2.000 0.003828 0.328039 19.240170 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorN2: N2(beta) = ECFG(3, 2, beta)/ECFG(2, 1, beta)^2 with E_theta +------------------------------------------------------------------------------------- + beta N2 obs + 0.100 0.502292 + 0.200 0.491827 + 0.500 0.472830 + 1.000 0.426794 + 1.500 0.366194 + 2.000 0.328039 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorN3: N3(beta) = ECFG(4, 2, beta)/ECFG(3, 1, beta)^2 with E_theta +------------------------------------------------------------------------------------- + beta N3 obs + 0.100 2.091587 + 0.200 2.107891 + 0.500 2.289141 + 1.000 4.159019 + 1.500 11.630367 + 2.000 19.240170 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorMseries: M_i(beta) = ECFG(i+1, 1, beta)/ECFN(i, 1, beta) with E_theta +------------------------------------------------------------------------------------- + beta N=1 N=2 N=3 + 0.100 0.198191 0.140807 0.091081 + 0.200 0.134916 0.127947 0.087436 + 0.500 0.049503 0.086121 0.076197 + 1.000 0.017389 0.030156 0.052638 + 1.500 0.010647 0.010921 0.024416 + 2.000 0.007953 0.006325 0.009110 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorM2: M2(beta) = ECFG(3, 1, beta)/ECFG(3, 1, beta) with E_theta +------------------------------------------------------------------------------------- + beta M2 obs + 0.100 0.140807 + 0.200 0.127947 + 0.500 0.086121 + 1.000 0.030156 + 1.500 0.010921 + 2.000 0.006325 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorCseries: C_i(beta) = ECFN(i-1, beta)*ECFN(i+1, beta)/ECFN(i, beta)^2 with E_theta +------------------------------------------------------------------------------------- + beta N=1 N=2 N=3 + 0.100 0.1981905437 0.3704251544 0.4761303115 + 0.200 0.1349155085 0.2930674174 0.4300888197 + 0.500 0.0495033749 0.2537948134 0.4348330516 + 1.000 0.0173890665 0.3067005727 0.4379497669 + 1.500 0.0106474132 0.3259256469 0.4135119649 + 2.000 0.0079526591 0.3539093501 0.3681613681 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorUseries: U_i(beta) = ECFG(i+1, 1, beta) with E_theta +------------------------------------------------------------------------------------- + beta N=1 N=2 N=3 + 0.100 0.1981905437 0.0279066279 0.0025417550 + 0.200 0.1349155085 0.0172619856 0.0015093221 + 0.500 0.0495033749 0.0042632685 0.0003248489 + 1.000 0.0173890665 0.0005243822 0.0000276026 + 1.500 0.0106474132 0.0001162849 0.0000028392 + 2.000 0.0079526591 0.0000503042 0.0000004582 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorU1: U1(beta) = ECFG(2, 1, beta) with E_theta +------------------------------------------------------------------------------------- + beta U1 obs + 0.100 0.1981905437 + 0.200 0.1349155085 + 0.500 0.0495033749 + 1.000 0.0173890665 + 1.500 0.0106474132 + 2.000 0.0079526591 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorU2: U2(beta) = ECFG(3, 1, beta) with E_theta +------------------------------------------------------------------------------------- + beta U2 obs + 0.100 0.0279066279 + 0.200 0.0172619856 + 0.500 0.0042632685 + 1.000 0.0005243822 + 1.500 0.0001162849 + 2.000 0.0000503042 +------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------------- +EnergyCorrelatorU3: U3(beta) = ECFG(4, 1, beta) with E_theta +------------------------------------------------------------------------------------- + beta U3 obs + 0.100 0.0025417550 + 0.200 0.0015093221 + 0.500 0.0003248489 + 1.000 0.0000276026 + 1.500 0.0000028392 + 2.000 0.0000004582 +------------------------------------------------------------------------------------- + diff --git a/src/Tools/fjcontrib/EnergyCorrelator/example_basic_usage.cc b/src/Tools/fjcontrib/EnergyCorrelator/example_basic_usage.cc new file mode 100644 --- /dev/null +++ b/src/Tools/fjcontrib/EnergyCorrelator/example_basic_usage.cc @@ -0,0 +1,234 @@ +// Example showing basic usage of energy correlator classes. +// +// Compile it with "make example" and run it with +// +// ./example_basic_usage < ../data/single-event.dat +// +// Copyright (c) 2013-2016 +// Andrew Larkoski, Lina Necib, Gavin Salam, and Jesse Thaler +// +// $Id: example.cc 958 2016-08-17 00:25:14Z linoush $ +//---------------------------------------------------------------------- +// 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 . +//---------------------------------------------------------------------- + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "fastjet/PseudoJet.hh" +#include "fastjet/ClusterSequence.hh" +#include "fastjet/JetDefinition.hh" + +#include +#include "EnergyCorrelator.hh" // In external code, this should be fastjet/contrib/EnergyCorrelator.hh + +using namespace std; +using namespace fastjet; +using namespace fastjet::contrib; + +// forward declaration to make things clearer +void read_event(vector &event); +void analyze(const vector & input_particles); + +//---------------------------------------------------------------------- +int main(){ + + //---------------------------------------------------------- + // read in input particles + vector event; + read_event(event); + cout << "# read an event with " << event.size() << " particles" << endl; + + //---------------------------------------------------------- + // illustrate how this EnergyCorrelator contrib works + + analyze(event); + + return 0; +} + +// read in input particles +void read_event(vector &event){ + string line; + while (getline(cin, line)) { + istringstream linestream(line); + // take substrings to avoid problems when there are extra "pollution" + // characters (e.g. line-feed). + if (line.substr(0,4) == "#END") {return;} + if (line.substr(0,1) == "#") {continue;} + double px,py,pz,E; + linestream >> px >> py >> pz >> E; + PseudoJet particle(px,py,pz,E); + + // push event onto back of full_event vector + event.push_back(particle); + } +} + +//////// +// +// Main Routine for Analysis +// +/////// + +void analyze(const vector & input_particles) { + + /////// EnergyCorrelator ///////////////////////////// + + // Initial clustering with anti-kt algorithm + JetAlgorithm algorithm = antikt_algorithm; + double jet_rad = 1.00; // jet radius for anti-kt algorithm + JetDefinition jetDef = JetDefinition(algorithm,jet_rad,E_scheme,Best); + ClusterSequence clust_seq(input_particles,jetDef); + vector antikt_jets = sorted_by_pt(clust_seq.inclusive_jets()); + + for (int j = 0; j < 1; j++) { // Hardest jet per event + if (antikt_jets[j].perp() > 200) { + + PseudoJet myJet = antikt_jets[j]; + + + // The angularity is set by the value of beta + double beta; + + cout << "-------------------------------------------------------------------------------------" << endl; + cout << "EnergyCorrelator: C series " << endl; + cout << "-------------------------------------------------------------------------------------" << endl; + printf("%7s %14s %14s %14s\n","beta", "C1", "C2", "C3"); + + + beta = 1.0; + //Defining the Cseries for beta= 1.0 + EnergyCorrelatorC1 C1s(beta); + EnergyCorrelatorC2 C2s(beta); + EnergyCorrelatorCseries C3s(3, beta); + + + printf("%7.3f %14.6f %14.6f %14.6f\n",beta,C1s(myJet),C2s(myJet),C3s(myJet)); + + beta = 2.0; + //Defining the Cseries for beta= 2.0 + EnergyCorrelatorC1 C1s_2(beta); + EnergyCorrelatorC2 C2s_2(beta); + EnergyCorrelatorCseries C3s_2(3, beta); + + + printf("%7.3f %14.6f %14.6f %14.6f\n",beta,C1s_2(myJet),C2s_2(myJet),C3s_2(myJet)); + + + cout << "-------------------------------------------------------------------------------------" << endl; + cout << "EnergyCorrelator: D2, orignal (alpha=beta) and generalized " << endl; + cout << "-------------------------------------------------------------------------------------" << endl; + printf("%7s %14s %14s %14s\n","beta","D2", "D2(alpha=1)", "D2(alpha=2)"); + + + beta = 1.0; + EnergyCorrelatorD2 d2(beta); + double alpha = 1.0; + EnergyCorrelatorGeneralizedD2 d2_generalized(alpha, beta); + alpha = 2.0; + EnergyCorrelatorGeneralizedD2 d2_generalized_2(alpha, beta); + + printf("%7.3f %14.6f %14.6f %14.6f\n",beta,d2(myJet), d2_generalized(myJet), d2_generalized_2(myJet)); + beta = 2.0; + EnergyCorrelatorD2 d2_2(beta); + alpha = 1.0; + EnergyCorrelatorGeneralizedD2 d2_generalized_3(alpha, beta); + alpha = 2.0; + EnergyCorrelatorGeneralizedD2 d2_generalized_4(alpha, beta); + printf("%7.3f %14.6f %14.6f %14.6f\n",beta,d2_2(myJet), d2_generalized_3(myJet), d2_generalized_4(myJet)); + + + cout << "-------------------------------------------------------------------------------------" << endl; + cout << "EnergyCorrelator: N series " << endl; + cout << "-------------------------------------------------------------------------------------" << endl; + printf("%7s %14s %14s\n","beta", "N2", "N3"); + + + beta = 1.0; + // Directly defining the EnergyCorrelator N2 and N3 + EnergyCorrelatorN2 N2(beta); + EnergyCorrelatorN3 N3(beta); + printf("%7.3f %14.6f %14.6f\n",beta, N2(myJet), N3(myJet)); + + beta = 2.0; + EnergyCorrelatorN2 N2_2(beta); + EnergyCorrelatorN3 N3_2(beta); + printf("%7.3f %14.6f %14.6f\n",beta, N2_2(myJet), N3_2(myJet)); + + + cout << "-------------------------------------------------------------------------------------" << endl; + cout << "EnergyCorrelator: M series " << endl; + cout << "-------------------------------------------------------------------------------------" << endl; + printf("%7s %14s %14s\n","beta", "M2", "M3"); + + beta = 1.0; + //Directly defining M2 + EnergyCorrelatorM2 M2(beta); + EnergyCorrelatorMseries M3(3, beta); + printf("%7.3f %14.6f %14.6f\n", beta, M2(myJet), M3(myJet)); + + beta = 2.0; + EnergyCorrelatorM2 M2_2(beta); + EnergyCorrelatorMseries M3_2(3, beta); + printf("%7.3f %14.6f %14.6f\n", beta, M2_2(myJet), M3_2(myJet)); + + cout << "-------------------------------------------------------------------------------------" << endl; + cout << "EnergyCorrelator: U series " << endl; + cout << "-------------------------------------------------------------------------------------" << endl; + printf("%7s %14s %14s %14s\n","beta", "U1", "U2", "U3"); + + + beta = 0.5; + //Defining the Useries for beta= 0.5 + EnergyCorrelatorU1 U1s(beta); + EnergyCorrelatorU2 U2s(beta); + EnergyCorrelatorU3 U3s(beta); + + + printf("%7.3f %14.8f %14.8f %14.8f\n",beta,U1s(myJet),U2s(myJet),U3s(myJet)); + + beta = 1.0; + //Defining the Useries for beta= 1.0 + EnergyCorrelatorU1 U1s_2(beta); + EnergyCorrelatorU2 U2s_2(beta); + EnergyCorrelatorU3 U3s_2(beta); + + + printf("%7.3f %14.8f %14.8f %14.8f\n",beta,U1s_2(myJet),U2s_2(myJet),U3s_2(myJet)); + + beta = 2.0; + //Defining the Useries for beta= 2.0 + EnergyCorrelatorU1 U1s_3(beta); + EnergyCorrelatorU2 U2s_3(beta); + EnergyCorrelatorU3 U3s_3(beta); + + + printf("%7.3f %14.8f %14.8f %14.8f\n",beta,U1s_3(myJet),U2s_3(myJet),U3s_3(myJet)); + } + } +} + + + + diff --git a/src/Tools/fjcontrib/EnergyCorrelator/example_basic_usage.ref b/src/Tools/fjcontrib/EnergyCorrelator/example_basic_usage.ref new file mode 100644 --- /dev/null +++ b/src/Tools/fjcontrib/EnergyCorrelator/example_basic_usage.ref @@ -0,0 +1,45 @@ +# read an event with 354 particles +#-------------------------------------------------------------------------- +# FastJet release 3.2.0 +# M. Cacciari, G.P. Salam and G. Soyez +# A software package for jet finding and analysis at colliders +# http://fastjet.fr +# +# Please cite EPJC72(2012)1896 [arXiv:1111.6097] if you use this package +# for scientific work and optionally PLB641(2006)57 [hep-ph/0512210]. +# +# FastJet is provided without warranty under the terms of the GNU GPLv2. +# It uses T. Chan's closest pair algorithm, S. Fortune's Voronoi code +# and 3rd party plugin jet algorithms. See COPYING file for details. +#-------------------------------------------------------------------------- +------------------------------------------------------------------------------------- +EnergyCorrelator: C series +------------------------------------------------------------------------------------- + beta C1 C2 C3 + 1.000 0.010342 0.097588 0.270742 + 2.000 0.001315 0.158011 0.279675 +------------------------------------------------------------------------------------- +EnergyCorrelator: D2, orignal (alpha=beta) and generalized +------------------------------------------------------------------------------------- + beta D2 D2(alpha=1) D2(alpha=2) + 1.000 9.435950 9.435950 223402.841831 + 2.000 120.129760 0.218810 120.129760 +------------------------------------------------------------------------------------- +EnergyCorrelator: N series +------------------------------------------------------------------------------------- + beta N2 N3 + 1.000 0.430942 1.979487 + 2.000 0.365859 6.700247 +------------------------------------------------------------------------------------- +EnergyCorrelator: M series +------------------------------------------------------------------------------------- + beta M2 M3 + 1.000 0.071101 0.099746 + 2.000 0.010448 0.036511 +------------------------------------------------------------------------------------- +EnergyCorrelator: U series +------------------------------------------------------------------------------------- + beta U1 U2 U3 + 0.500 0.05378579 0.00805855 0.00105555 + 1.000 0.01034218 0.00073534 0.00007335 + 2.000 0.00131534 0.00001374 0.00000050 diff --git a/src/Tools/fjcontrib/Makefile.am b/src/Tools/fjcontrib/Makefile.am new file mode 100644 --- /dev/null +++ b/src/Tools/fjcontrib/Makefile.am @@ -0,0 +1,1 @@ +SUBDIRS = EnergyCorrelator Nsubjettiness