Index: contrib/contribs/Nsubjettiness/trunk/VERSION =================================================================== --- contrib/contribs/Nsubjettiness/trunk/VERSION (revision 1410) +++ contrib/contribs/Nsubjettiness/trunk/VERSION (revision 1411) @@ -1 +1 @@ -2.3.1 \ No newline at end of file +2.3.2 \ No newline at end of file Index: contrib/contribs/Nsubjettiness/trunk/NEWS =================================================================== --- contrib/contribs/Nsubjettiness/trunk/NEWS (revision 1410) +++ contrib/contribs/Nsubjettiness/trunk/NEWS (revision 1411) @@ -1,108 +1,110 @@ ------------------------- Version 2 ------------------------- This is a streamlining of the N-subjettiness code, developed mainly by TJ Wilkason. The core functionality is unchanged, but classes have been dramatically reorganized to allow for later expansion. Because the API for Njettiness has changed, we have called this v2 (http://semver.org). Note that we have maintain backwards compatibility for the typical ways that Nsubjettiness was used. In particular, all of the Nsubjettiness class code in the example file from v1.0.3 still compiles, as does the NjettinessPlugin class code that uses the default measure. The key new features are: * NsubjettinessRatio: Direct access to tau_N / tau_M (the most requested feature) * MeasureDefinition to allow access to normalized and unnormalized measures * AxesDefinition to allow for access to more general axes modes * Winner-Take-All recombination axes: a faster way to find axes than beta=1 minimization, but with comparable performance. * TauComponents to get access to the pieces of the N-(sub)jettiness calculation. * TauExtras to get complete access to get partitioning and axes information. * For clarity, split the example file into an example_basic_usage and example_advanced_usage. * In Nsubjettiness, access to seedAxes() and currentAxes() to figure out the axes used before and after minimization. * In Nsubjettiness, access to currentSubjets() to get the subjet fourvectors. * (v2.2) XConePlugin, which improves on the previous NjettinessPlugin to use N-jettiness as a jet finder using the new ConicalGeometric measure. +-- 2.3.2: (February 27, 2024) Fixing unsigned int warning in + example_advanced_usage.cc (thanks Gregory Soyez) -- 2.3.1: (February 27, 2024) Fixing unsigned int warning in AxesDefinition.hh (thanks Gregory Soyez) -- 2.3.0: (February 23, 2024) Changed recommended axes away from OnePass; added HalfKT axes options recommended for beta = 2; updated example files; removed ee testing since ee axes choices are not available -- 2.2.6: (June 13, 2022) Removed "static" for thread safety (thanks Tomek Procter and Andy Buckley) -- 2.2.5: (June 6, 2018) Fixed bug involved undefined pointer for in AxesDefinition (thanks Attila Krasznahorkay) -- 2.2.4: (Jun 14, 2016) Fixed bug where multi-pass minimization could yield pathological axes (thanks Gregory Soyez) -- 2.2.3: (Apr 4, 2016) Fixed bug where a jet with fewer than N constituents could give random value for tau_N (thanks Nathan Hartland) -- 2.2.2: (Mar 29, 2016) Updating SharedPtr interface for FJ 3.2 -- 2.2.1: (Sept 28, 2015) Fix of small Makefile bug -- 2.2.0: (Sept 7, 2015) Inclusion of the XCone jet algorithm, as well as a few new measures, including the "conical geometric" measure and options for e+e- colliders. Improvement of the Measure/AxesDefinition interface to allow for direct use in calculations. * Fixed bug where MultiPass_Axes did not actually minimize * Fixed floating point error with infinity^2 in various measures -- 2.1.0: (July 9, 2014) Inclusion of Measure/AxesDefinition interface. This was the first publicly available version of Nsubjettiness v2. -- 2.0.0: Initial release of v2.0. This was never officially made public. ------------------------- Version 1 ------------------------- This was a new release using FastJet contrib framework, primary developed by Jesse Thaler. -- 1.0.3: Added inlines to fix compile issue (thanks Matthew Low) -- 1.0.2: Fixed potential dependency issue (thanks FJ authors) -- 1.0.1: Fixed memory leak in Njettiness.hh (thanks FJ authors) -- 1.0.0: New release using FastJet contrib framework. This includes a new makefile and a simplified example program. ------------------------- Previous Versions ------------------------- The previous versions of this code were developed initially by Ken Van Tilburg, tweaked by Jesse Thaler, and made into a robust FastJet add on by Chris Vermilion. Previous versions available from: http://jthaler.net/jets/Njettiness-0.5.1.tar.gz (Experimental Version) http://jthaler.net/jets/Njettiness-0.4.1.tar.gz (Stable Version) Previous version history: -- 0.5.1: For Njettiness Plugin, added access to currentTau values and axes via ClusterSequence::Extras class. (thanks to Dinko Ferencek and John Paul Chou) -- 0.5.0: Corrected fatal error in ConstituentTauValue (TauValue unaffected). Started process of allowing for more general measures and alternative minimization schemes. Extremely preliminary inclusion of alternative "geometric" measure. -- 0.4.1: Corrected bug where a too-small value of Rcut would cause the minimization procedure to fail (thanks Marat Freytsis, Brian Shuve) -- 0.4.0: Adding Nsubjettiness FunctionOfPseudoJet. Re-organizing file structure and doing some re-naming to clarify Njettiness vs. Nsubjettiness. Some speedup in UpdateAxes code. (CKV) -- 0.3.2: Returns zero instead of a segmentation fault when the number of particles in a jet is smaller than the N value in tau_N (thanks Grigory Ovanesyan) -- 0.3.2: Fixed -Wshadow errors (thanks Grigory Ovanesyan) -- 0.3.1: Fixed stray comma/semicolon compiler error (thanks Grigory Ovanesyan) -- 0.3.1: Corrected tarbomb issue (thanks Jonathan Walsh) -- 0.3.1: Added anti-kT seeds as option -- 0.3.1: Fixed bug in minimization code with R_cutoff (thanks Chris Vermilion) -- 0.3.1: Added getPartition() and getJets() functions as helper functions for Chris Vermilion. (JT) Index: contrib/contribs/Nsubjettiness/trunk/ChangeLog =================================================================== --- contrib/contribs/Nsubjettiness/trunk/ChangeLog (revision 1410) +++ contrib/contribs/Nsubjettiness/trunk/ChangeLog (revision 1411) @@ -1,348 +1,351 @@ +2024-02-28 + Fixing unsigned int warning in example_advanced_usage.cc + Updated NEWS and VERSION getting ready for v2.3.2 release 2024-02-27 Fixing unsigned int warning in AxesDefinition.hh - Updated README and VERSION getting ready for v2.3.1 release + Updated NEWS and VERSION getting ready for v2.3.1 release 2024-02-22 Updated README and VERSION getting ready for v2.3.0 release Removed extraneous "-lm" from Makefile Added HalfKT options to AxesDefinition.hh Updated example files with new recommended usage and higher accuracy for OnePass minimization Removed "ee" example file, since it uses ee measure with pp axes (ee axes not available in current version) 2022-06-13 Removed -std=c++11 flag from makefile Updated MeasureDefinition.cc to remove "static thread_local" (since it doesn't really seem to help with timing) 2022-06-10 Updated makefile with -std=c++11 flag Updated MeasureDefinition.cc with thread_local for thread safety 2018-07-06 Updated comments in AxesDefinition.hh about role of JetDefinitionWrapper Updated AUTHORS with JHEP publication information for XCone Prepared VERSION and NEWS for 2.2.5 release 2018-07-05 Fixed bug in AxesDefinition.hh where _recomb was used before it was declared. 2016-06-08 Fixed bug in MeasureDefinition.cc where axes were not completely defined, leading to problems with multi-pass axes 2016-04-04 Fixed Njettiness.cc to give value of _current_tau_components even if less than N constituents Delete extraneous code in example_advanced_usage.cc 2016-03-29 Update for FJ 3.2.0 to deal with SharedPtr () deprecation 2015-09-28 Updated NEWS for 2.2.1 release. 2015-09-18 Fixed duplicate XConePlugin entry in Makefile. 2015-08-20 Trying to fix "abs" bug in ExtraRecombiners.cc 2015-08-19 Adding arXiv numbers to XCone papers Used this_jet in example_basic_usage. Fixed typo in example_advanced_usage header. Added copy/paste code in README file. 2015-08-13 Ready for 2.2.0 release 2015-07-23 Fixed typo in GenET_GenKT_Axes error message Added _too_few_axes_warning to ExclusiveJetAxes and ExclusiveCombinatorialJetAxes Switched to ../data/single_event_ee.dat for make check 2015-07-20 Renamed WinnerTakeAllRecombiner.hh/cc to ExtraRecombiners.hh/cc Added _too_few_axes_warning to HardestJetAxes Added GenKT_Axes and OnePass_GenKT_Axes and Comb_GenKT_Axes (using E-scheme recombination). Added warning about p < 0 or delta <=0 in GenKT axes finders. Added warning about beta <= 0 in all measures. 2015-07-10 Putting in small tweaks in documentation to get ready for 2.2 release candidate 1. 2015-06-15 Starting doxygen file for eventual improved documentation. Starting long process of improving documentation throughout. Made the basic usage file a bit easier to read. Adding in LimitedWarnings for old style constructors 2015-06-12 Synchronized definition of new measures with XCone paper. In MeasureDefinition, added default values of jet_distance_squared and beam_distance_squared for cases where we don't want to optimize specifically. Fixed bug in OriginalGeometricMeasure and ModifiedGeometric Measure Commented out DeprecatedGeometricMeasure and DeprecatedGeometricCutoffMeasure since they were only causing confusion 2015-05-26 Removed axis_scale_factor(), added bool to calculate this value if needed to save computation time Defined small offset in denominator of axis scaling according to accuracy of refinement Updated advanced examples to include tau values and number of jet constituents 2015-05-25 Clean up of AxesDefinition Splitting get_axes into get_starting_axes and get_refined axes Putting in proper noise ranges (hopefully) for MultiPass Clean up of MeasureDefinition, rename jet_gamma to beam_gamma Put in zero checking for jet_distance in ConicalGeometricMeasure Added in ConicalMeasure for consistency Changing OnePass Minimization to allow for temporary uphill 2015-05-24 Added Combinatorial GenET_GenKT_Axes and MultiPass_Manual_Axes Moved Axes refining information into MeasureDefinition, associated each measure with corresponding axes refiner Moved get_one_pass_axes into MeasureDefinition, removed any mention of Npass Moved all information on number of passes to AxesDefinition Made AxesRefiner.hh/.cc into defunct files 2015-05-22 Cleaning out commented text. Renaming classes to be consistent with recommended usage. 2015-05-22 Added XConePlugin as a specific implementation of NjettinessPlugin Added usage of XCone beta = 1.0 and beta = 2.0 to both basic and advanced example files Added OriginalGeometric, ModifiedGeometric, ConicalGeometric, and XCone measures to list of test measures Added OnePass_GenRecomb_GenKT_Axes to list of test axes Added description to XCone measure in MeasureDefinition 2015-05-21 Updated minimization scheme to avoid divide-by-zero errors Fixed various factors of 2 in the definition of the measures 2015-04-19 Fixed bug in minimization scheme for GeneralAxesRefiner Moved measure_type to DefaultMeasure, removed geometric measure from e+e- example file 2015-03-22 Added OriginalGeometricMeasure and ModifiedGeometricMeasure definitions Changed all instances of GeometricMeasure to DeprecatedGeometricMeasure, and added error statements Made GeneralAxesRefiner the default axes refiner for Measure Definition, overwritten by DefaultMeasure and GeometricMeasure Created DefaultMeasure class for all the conical measure subclasses Separated out e+e- and pp measures into separate example files 2015-03-09 Added ConicalGeometric measures with jet_beta and jet_gamma definitions Added XCone measures derived from ConicalGeometric with jet_gamma = 1.0 Added GeneralAxesRefiner for use with any measure (currently defined with XCone measure) Added axes_numerator in MeasureDefinition to define the momentum scaling for minimization (currently only defined for Conical Geometric measure) 2014-11-28 Minor change to default parameters in axes definition 2014-10-08 Updated example file with new e+e- measure definitions Added measure type to measure definition descriptions Changed order of parameters in new axes definitions Added standard C++ epsilon definition to GeneralERecombiner 2014-10-07 Updated example_advanced_usage with new axes choices Reversed inheritance of NormalizedMeasure and NormalizedCutoffMeasure (and Geometric) back to original Storing _RcutoffSq as separate variable, and recalculating it in NormalizedMeasure Cleaning up ExclusiveCombinatorialJetAxes and added comments to explain the process Fixed memory leaks using delete_recombiner_when_unused() Fixed manual axes bug in Njettiness Cleaned up enum definitions 2014-10-01 Added new parameterized recombination scheme to Winner-Take-All recombiner Created Winner-Take-All GenKT and general Recomb GenKT axes finders and onepass versions Created new N choose M minimization axis finder, created N choose M WTA GenKT axis finder as example Removed NPass as constructor argument in AxesDefinition, made it set through protected method Removed TauMode as constructor argument in MeasureDefinition, made it set through protected method Flipped inheritance of NormalizedMeasure and NormalizedCutoffMeasure (same for Geometric) to remove error of squaring the integer maximum Created new MeasureType enum to allow user to choose between pp and ee variables (ee variables need testing) Updated MeasureDefinition constructors to take in extra MeasureType parameter (but defaulted to pp variables) Added new Default TauMode argument Fixed unsigned integers in various places Added setAxes method to NjettinessPlugin 2014-08-26 Enhanced TauComponents to include more infomation NjettinessExtras now inherits from TauComponents Removed getPartition from Njettiness, to avoid code duplication Fixed double calculating issue in NjettinessPlugin::run_clustering() Now AxesDefinition can use measure information without running AxesRefiner Added TauStructure so the jets returned by TauComponents can know their tau value. 2014-08-25 Merged MeasureDefinition and MeasureFunction into new MeasureDefinition. Merged StartingAxesFinder and AxesDefinition into new AxesDefinition. Renamed AxesFinder.cc/hh to AxesRefiner.cc/hh Renamed NjettinessDefinition.cc/hh to AxesDefinition.cc/hh Renamed MeasureFunction.cc/hh to MeasureDefinition.cc/hh Renaming result() function in MeasureDefinition to be consistent with Nsubjettiness interface. Split off TauComponents into separate header Added TauPartition class for readability of partitioning Moved NjettinessExtras into TauComponents, as this will eventually be the logical location Added cross check of new MeasureDefinition and AxesDefinition in example_advanced_usage. Lots of comments updated. Changed version number to 2.2.0-alpha-dev, since this is going to be a bigger update than I had originally thought 2014-08-20 Incorporated code in NjettinessPlugin to handle FJ3.1 treatment of auto_ptr (thanks Gregory) Changed version number to 2.1.1-alpha-dev Split AxesFinder into StartingAxesFinder and RefiningAxesFinder for clarity. Manual axes mode now corresponds to a NULL StartingAxesFinder in Njettiness (so removed AxesFinderFromUserInput) Added AxesRefiningMode to make selection of minimization routine more transparent in Njettiness Moved sq() to more appropriate place in AxesFinder.hh Rearranged Nsubjettiness.hh to make the old code less visible. Renamed AxesFinderFromOnePassMinimization -> AxesFinderFromConicalMinimization Renamed DefaultUnnormalizedMeasureFunction -> ConicalUnnormalizedMeasureFunction Removed supportsMultiPassMinimization() from MeasureDefinition since any One Pass algorithm can be multipass. 2014-07-09 Changed version for 2.1.0 release. Updated NEWS to reflect 2.1.0 release 2014-07-07 Added forward declaration of options in NjettinessDefinition for readability. Updated README with some clarifications Added usage information in the example file Reran svn propset svn:keywords Id *.cc *.hh 2014-06-25 Declaring release candidate of 2.1 2014-06-11 Fixed virtual destructor issue in AxesFinder Changing copy() to create() in NjettinessDefinition for "new" clarity Converted some SharedPtr to regular pointers in NjettinessDefinition to be consistent on meaning of "create" commands. 2014-06-10 Slight modification of example_advanced_usage Fixed bug in GeometricCutoffMeasure (incorrect denominator setting) 2014-06-05 Moved public before private in the .hh files for readability Starting process of switching to SharedPtr internally Clean up of AxesFinderFromGeometricMinimization Simplified AxesFinder interface such that it doesn't know about starting axes finders (this is now handled in Njettiness). Added const qualifiers in Njettiness 2014-06-04 Implemented AxesDefinition class Added descriptions to AxesDefinition and MeasureDefinition Simplified example_advanced_usage with new Definitions Made copy constructor private for Njettiness, to avoid copying 2014-06-03 Implemented remaining suggestions from FJ authors (Thanks!) Fixed bug in example_advanced_usage where wrong beta value was used for NjettinessPlugin tests. Removed NANs as signals for number of parameters in Nsubjettiness and NjettinessPlugin Reduced the number of allowed parameters from 4 to 3. Wrapped NEWS to 80 characters Added MeasureDefinition as way to safely store information about the measures used Converted a few NANs to std::numeric_limits::quiet_NaN() when a parameter shouldn't be used. Added AxesStruct and MeasureStruct to simplify the form of example_advanced_usage Added example_v1p0p3 to check for backwards compatibility with v1.0.3 Changed the names of the MeasureFunctions in order to avoid conflicts with the new MeasureDefinitions Fixed bug in correlation between subjets and tau values in NjettinessPlugin Added currentTauComponents to Nsubjettiness Added subTau information to example_basic_usage Added file NjettinessDefinition to hold MeasureDefinition Changed Njettiness constructors to treat MeasureSpecification as primary object Fixed segmentation fault with ClusterSequenceAreas 2014-06-02 Implemented many suggestions from FJ authors (Thanks!) Removed FastJet 2 specific code Made sq() function into internal namespace (as "inline static" to avoid conflicts with other packages) Made setAxes() take const reference argument Rewrapped README to 80 characters and updated/improved some of the descriptions Clarified NEWS about what parts of the Nsubjettiness code is backwards compatible with v1.0.3 Clarified the para choices in Nsubjettiness constructor 2014-04-30 Added (void)(n_jets) in AxesFinder.hh to fix unused-parameter warning 2014-04-29 Added manual definition of NAN for compilers that don't have it. Removed a few more unused parameters for compilation 2014-04-22 Turned on -Wunused-parameter compiler flag to fix ATLAS compile issues. 2014-04-18 Tweaks to NEWS and README. Preparing for 2.0.0-rc1 release. 2014-04-16 Decided that enough has changed that this should be v2.0 Added Id tags 2014-04-14 Added get_partition_list to MeasureFunction Removed do_cluster from MeasureFunction (no longer needed) Fixed bug with NjettinessPlugin where jets were listed in backwards order from axes. Removed various commented out pieces of code. 2014-03-16 Added partitioning information to Nsubjettiness Partitioning is now calculated in MeasureFunction and stored by Njettiness. Rewrote MeasureFunction result() to call result_from_partition() Added subjet (and constituent counting) information to example_basic_usage Commented out redundant "getJets" function 2014-02-25 Added access to seedAxes used for one-pass minimization routines. Added axes print out to example_basic_usage, and fixed too many PrintJets declarations 2014-02-24 Fixed embarrassing bug with min_axes (error introduced after v1.0 to make it the same as onepass_kt) Simplified GeometricMeasure and added possibility of beta dependence Commented out WTA2 options, since those have not been fully tested (nor do they seem particularly useful at the moment). They can be reinstated if the physics case can be made to use them. Split example into example_basic_usage and example_advanced_usage 2014-01-28 Added new options in WinnerTakeAllRecombiner to use either pT or pT^2/E to recombine particles 2014-01-24 Added access to currentAxes from Nsubjettiness. 2014-01-18 Added beam regions to MeasureFunction, correspondingly renamed functions to have jet and beam regions Renamed functions in TauComponents for consistency with MeasureFunction Adding debugging code to AxesFinderFromOnePassMinimization::getAxes Worked extensively on example.cc to make sure that it tests all available options. Rewrote PrintJets command in example.cc for later improvements Converted some magic numbers to std::numeric_limits::max() 2014-01-17 Rewrote KMeansMinimization to call OnePassMinimization, adding noise explicitly. Removed any nothing of noise from OnePassMinimization Removed Double32_t for root usage is Nsubjettiness Clean up of many comments throughout the code, updating of README file Removed unnecessary establishAxes in Njettiness Removed bare constructor for Njettiness to avoid incompatibility with enum choices, may reinstate later. Also removed setMeasureFunction, setAxesFinder for same reason NjettinessExtras now calls TauComponents 2014-01-16 Moved minimization functions to OnePassMinimization, changed KMeansMinimization class to simply call OnePassMinimization a specified number of times Added extra tau function in TauComponents for users to get tau directly Changed radius parameter in AxesFinderFromExclusiveJet subclasses to use max_allowable_R Updated example.ref to account for changes due to change in radius parameter 2014-01-15 Changed NjettinessComponents to TauComponents Updated MeasureFunction with "result" function that returns TauComponents object TauComponents changed to calculate all tau components given subtaus_numerator and tau_denominator Njettiness updated to return TauComponents object rather than individual components Nsubjettiness and NjettinessPlugin updated to have option for 4th parameter 2014-01-14 Added NjettinessComponents class so Njettiness does not recalculate tau values Removed old Njettiness constructors, updated Nsubjettiness and NjettinessPlugin constructors to use new constructor Added geometric minimization to OnePassAxesFinders Created new Njettiness function to set OnePassAxesFinders to reduce code Updated LightLikeAxis with ConvertToPseudoJet function Updated README with new functionality of code 2014-01-12 Removed NsubGeometricParameters in all functions/constructors, replaced with Rcutoff double Added three new measure mode options where Rcutoff is declared explicitly in parameters Added checks so minimization axes finders are not used for geometric measures AxesFinderFromOnePassMinimization class created as child of AxesFinderFromKmeansMinimization Added new NsubjettinessRatio constructor to include MeasureMode option Moved AxesFinder and MeasureFunction declarations from AxesMode and MeasureMode into separate Njettiness function Removed R0 from AxesFinderFromKmeansMinimization Changed example.cc to get rid of use of NsubGeometricParameters 2014-01-9 Removed NsubParameters in all functions/constructors, replaced with three separate parameters Added checks for correct number of parameters in Njettiness constructor 2014-01-8 Removed normalization information from Nsubjettiness Added flag to MeasureFunction to give option of using the denominator Split DefaultMeasure into separate normalized and unnormalized classes 2014-01-7 Added capability of choosing a specific Measure in Njettiness Added new Nsubjettiness constructor to allow choice of both AxesMode and MeasureMode 2014-01-6 Updated copyright information Fixed bug in WinnerTakeAllRecombiner Moved KMeansParameters to AxesFinder Updated README with descriptions of new header files 2013-12-30 Changed name of MeasureFunctor to MeasureFunction Created separate .hh/.cc files for MeasureFunction, AxesFinder, and WinnerTakeAllRecombiner Updated Makefile to account for new files Removed getMinimumAxes in AxesFinderFromKMeansMinimization, consolidated with getAxes Updated comments on classes and major functions 2013-12-22 Created .cc files and moved all function definitions into .cc files Updated Makefile to account for new .cc files 2013-11-12 Added to fjcontrib svn 2013-11-12 Debugging svn 2013-11-11 Changed MeasureFunctor to separately treat tau numerator and denominator Changed some of the function names in MeasureFunctor. Should not affect users Added more informative function names to Njettiness. Njettiness now allows finding unnormalized tau values Added WTARecombiner to define winner-take-all axes Added various WTA options to AxesMode Added setAxes to Nsubjettiness Added NsubjettinessRatio function 2013-08-26 Added inlines to fix compile issue Put some of the minimization code inside of the AxesFinderFromKmeansMinimization class 2013-02-23 Fixed dependency issue (now using make depend) 2013-02-22 Fixed memory management and failed make check issues. 2013-02-21 First version submitted to fjcontrib 2013-02-20 Initial creation based on previous plugin hosted at http://www.jthaler.net/jets/ Index: contrib/contribs/Nsubjettiness/trunk/example_advanced_usage.cc =================================================================== --- contrib/contribs/Nsubjettiness/trunk/example_advanced_usage.cc (revision 1410) +++ contrib/contribs/Nsubjettiness/trunk/example_advanced_usage.cc (revision 1411) @@ -1,996 +1,996 @@ // Nsubjettiness Package // Questions/Comments? jthaler@jthaler.net // // Copyright (c) 2011-13 // Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason // // Run this example with // ./example_advanced_usage < ../data/single-event.dat // // $Id$ //---------------------------------------------------------------------- // 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/ClusterSequenceArea.hh" #include #include "Nsubjettiness.hh" // In external code, this should be fastjet/contrib/Nsubjettiness.hh #include "Njettiness.hh" #include "NjettinessPlugin.hh" #include "XConePlugin.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 Nsubjettiness contrib works analyze(event); return 0; } // Simple class to store Axes along with a name for display class AxesStruct { private: // Shared Ptr so it handles memory management SharedPtr _axes_def; public: AxesStruct(const AxesDefinition & axes_def) : _axes_def(axes_def.create()) {} // Need special copy constructor to make it possible to put in a std::vector AxesStruct(const AxesStruct& myStruct) : _axes_def(myStruct._axes_def->create()) {} const AxesDefinition & def() const {return *_axes_def;} AxesDefinition & def_for_change() {return *_axes_def;} string description() const {return _axes_def->description();} string short_description() const {return _axes_def->short_description();} }; // Simple class to store Measures to make it easier to put in std::vector class MeasureStruct { private: // Shared Ptr so it handles memory management SharedPtr _measure_def; public: MeasureStruct(const MeasureDefinition& measure_def) : _measure_def(measure_def.create()) {} // Need special copy constructor to make it possible to put in a std::vector MeasureStruct(const MeasureStruct& myStruct) : _measure_def(myStruct._measure_def->create()) {} const MeasureDefinition & def() const {return *_measure_def;} string description() const {return _measure_def->description();} }; // 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); } } // Helper Function for Printing out Jet Information void PrintJets(const vector & jets, bool commentOut = false); void PrintAxes(const vector & jets, bool commentOut = false); void PrintJetsWithComponents(const vector & jets, bool commentOut = false); //////// // // Main Routine for Analysis // /////// void analyze(const vector & input_particles) { //////// // // This code will check multiple axes/measure modes // First thing we do is establish the various modes we will check // /////// //Define characteristic test parameters to use here double p = 0.5; double delta = 10.0; // close to winner-take-all. TODO: Think about right value here. double R0 = 0.2; double Rcutoff = 0.5; double infinity = std::numeric_limits::max(); int nExtra = 2; int NPass = 10; // A list of all of the available axes modes vector _testAxes; _testAxes.push_back(HalfKT_Axes()); _testAxes.push_back(KT_Axes()); _testAxes.push_back(CA_Axes()); _testAxes.push_back(AntiKT_Axes(R0)); _testAxes.push_back(WTA_HalfKT_Axes()); _testAxes.push_back(WTA_KT_Axes()); _testAxes.push_back(WTA_CA_Axes()); _testAxes.push_back(GenKT_Axes(p, R0)); _testAxes.push_back(WTA_GenKT_Axes(p, R0)); _testAxes.push_back(GenET_GenKT_Axes(delta, p, R0)); _testAxes.push_back(OnePass_HalfKT_Axes()); _testAxes.push_back(OnePass_KT_Axes()); _testAxes.push_back(OnePass_AntiKT_Axes(R0)); _testAxes.push_back(OnePass_WTA_HalfKT_Axes()); _testAxes.push_back(OnePass_WTA_KT_Axes()); _testAxes.push_back(OnePass_GenKT_Axes(p, R0)); _testAxes.push_back(OnePass_WTA_GenKT_Axes(p, R0)); _testAxes.push_back(OnePass_GenET_GenKT_Axes(delta, p, R0)); _testAxes.push_back(Comb_GenKT_Axes(nExtra, p, R0)); _testAxes.push_back(Comb_WTA_GenKT_Axes(nExtra, p, R0)); _testAxes.push_back(Comb_GenET_GenKT_Axes(nExtra, delta, p, R0)); // manual axes (should be identical to kt axes) _testAxes.push_back(Manual_Axes()); _testAxes.push_back(OnePass_Manual_Axes()); // these axes are not checked during make check since they do not give reliable results _testAxes.push_back(OnePass_CA_Axes()); // not recommended _testAxes.push_back(OnePass_WTA_CA_Axes()); // not recommended _testAxes.push_back(MultiPass_Axes(NPass)); _testAxes.push_back(MultiPass_Manual_Axes(NPass)); int num_unchecked = 4; // number of unchecked axes // For testing purposes on different platforms, increase number of iterations // and accuracy for OnePass minimization to be less sensitive to floating point rounding - for (int i = 0; i < _testAxes.size(); i++) { + for (unsigned int i = 0; i < _testAxes.size(); i++) { if (_testAxes[i].def().nPass() == 1) { _testAxes[i].def_for_change().setNPass(1, 10000, 0.00001); // defaut is 1, 1000, 0.0001 } } // // Note: Njettiness::min_axes is not guarenteed to give a global // minimum, only a local minimum, and different choices of the random // number seed can give different results. For that reason, // the one-pass minimization are recommended over min_axes. // // Getting a smaller list of recommended axes modes // These are the ones that are more likely to give sensible results (and are all IRC safe) vector _testRecommendedAxes; _testRecommendedAxes.push_back(HalfKT_Axes()); _testRecommendedAxes.push_back(WTA_KT_Axes()); // new axes options added in most recent version of Nsubjettiness // these are separate from above since they should only be defined with a cutoff value for sensible results vector _testAlgorithmRecommendedAxes; _testAlgorithmRecommendedAxes.push_back(GenET_GenKT_Axes(1.0, 1.0, Rcutoff)); _testAlgorithmRecommendedAxes.push_back(GenET_GenKT_Axes(infinity, 1.0, Rcutoff)); _testAlgorithmRecommendedAxes.push_back(GenET_GenKT_Axes(1.0, 0.5, Rcutoff)); _testAlgorithmRecommendedAxes.push_back(OnePass_GenET_GenKT_Axes(1.0, 1.0, Rcutoff)); _testAlgorithmRecommendedAxes.push_back(OnePass_GenET_GenKT_Axes(infinity, 1.0, Rcutoff)); _testAlgorithmRecommendedAxes.push_back(OnePass_GenET_GenKT_Axes(1.0, 0.5, Rcutoff)); // Getting some of the measure modes to test // (When applied to a single jet we won't test the cutoff measures, // since cutoffs aren't typically helpful when applied to single jets) // Note that we are calling measures by their MeasureDefinition vector _testMeasures; _testMeasures.push_back( NormalizedMeasure(1.0, 1.0, pt_R)); _testMeasures.push_back(UnnormalizedMeasure(1.0 , pt_R)); _testMeasures.push_back( NormalizedMeasure(2.0, 1.0, pt_R)); _testMeasures.push_back(UnnormalizedMeasure(2.0 , pt_R)); // When doing Njettiness as a jet algorithm, want to test the cutoff measures. // (Since they are not senisible without a cutoff) vector _testCutoffMeasures; _testCutoffMeasures.push_back(UnnormalizedCutoffMeasure(1.0, Rcutoff, pt_R)); _testCutoffMeasures.push_back(UnnormalizedCutoffMeasure(2.0, Rcutoff, pt_R)); // new measures added in the most recent version of NSubjettiness _testCutoffMeasures.push_back(ConicalMeasure(1.0, Rcutoff)); _testCutoffMeasures.push_back(ConicalMeasure(2.0, Rcutoff)); _testCutoffMeasures.push_back(OriginalGeometricMeasure(Rcutoff)); _testCutoffMeasures.push_back(ModifiedGeometricMeasure(Rcutoff)); _testCutoffMeasures.push_back(ConicalGeometricMeasure(1.0, 1.0, Rcutoff)); _testCutoffMeasures.push_back(ConicalGeometricMeasure(2.0, 1.0, Rcutoff)); _testCutoffMeasures.push_back(XConeMeasure(1.0, Rcutoff)); // Should be identical to ConicalGeometric _testCutoffMeasures.push_back(XConeMeasure(2.0, Rcutoff)); /////// N-subjettiness ///////////////////////////// //////// // // Start of analysis. First find anti-kT jets, then find N-subjettiness values of those jets // /////// // 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()); // clust_seq.delete_self_when_unused(); // small number to show equivalence of doubles double epsilon = 0.0001; for (int j = 0; j < 2; j++) { // Two hardest jets per event if (antikt_jets[j].perp() < 200) continue; cout << "-----------------------------------------------------------------------------------------------" << endl; cout << "Analyzing Jet " << j + 1 << ":" << endl; cout << "-----------------------------------------------------------------------------------------------" << endl; //////// // // Basic checks of tau values first // // If you don't want to know the directions of the subjets, // then you can use the simple function Nsubjettiness. // // Recommended usage for Nsubjettiness: // AxesMode: kt_axes, wta_kt_axes, onepass_kt_axes, or onepass_wta_kt_axes // MeasureMode: unnormalized_measure // beta with kt_axes: 2.0 // beta with wta_kt_axes: anything greater than 0.0 (particularly good for 1.0) // beta with onepass_kt_axes or onepass_wta_kt_axes: between 1.0 and 3.0 // /////// cout << "-----------------------------------------------------------------------------------------------" << endl; cout << "Outputting N-subjettiness Values" << endl; cout << "-----------------------------------------------------------------------------------------------" << endl; // Now loop through all options cout << setprecision(6) << right << fixed; for (unsigned iM = 0; iM < _testMeasures.size(); iM++) { cout << "-----------------------------------------------------------------------------------------------" << endl; cout << _testMeasures[iM].description() << ":" << endl; cout << setw(25) << "AxisMode" << setw(14) << "tau1" << setw(14) << "tau2" << setw(14) << "tau3" << setw(14) << "tau2/tau1" << setw(14) << "tau3/tau2" << endl; for (unsigned iA = 0; iA < _testAxes.size(); iA++) { // Current axes/measure modes and particles const PseudoJet & my_jet = antikt_jets[j]; const vector particles = my_jet.constituents(); const AxesDefinition & axes_def = _testAxes[iA].def(); const MeasureDefinition & measure_def = _testMeasures[iM].def(); // This case doesn't work, so skip it. // if (axes_def.givesRandomizedResults()) continue; // define Nsubjettiness functions Nsubjettiness nSub1(1, axes_def, measure_def); Nsubjettiness nSub2(2, axes_def, measure_def); Nsubjettiness nSub3(3, axes_def, measure_def); // define manual axes when they are necessary (should be identical to KT_Axes) if (axes_def.needsManualAxes()) { JetDefinition manual_jetDef(fastjet::kt_algorithm, fastjet::JetDefinition::max_allowable_R, // new WinnerTakeAllRecombiner(), fastjet::E_scheme, fastjet::Best); fastjet::ClusterSequence manual_clustSeq(particles, manual_jetDef); nSub1.setAxes(manual_clustSeq.exclusive_jets(1)); nSub2.setAxes(manual_clustSeq.exclusive_jets(2)); nSub3.setAxes(manual_clustSeq.exclusive_jets(3)); } // calculate Nsubjettiness values double tau1 = nSub1(my_jet); double tau2 = nSub2(my_jet); double tau3 = nSub3(my_jet); //These should only happen if the axes are not manual and are not multipass double tau21, tau32; if (!axes_def.needsManualAxes() && !axes_def.givesRandomizedResults()) { // An entirely equivalent, but painful way to calculate is: double tau1alt = measure_def(particles,axes_def(1,particles,&measure_def)); double tau2alt = measure_def(particles,axes_def(2,particles,&measure_def)); double tau3alt = measure_def(particles,axes_def(3,particles,&measure_def)); assert(tau1alt == tau1); assert(tau2alt == tau2); assert(tau3alt == tau3); NsubjettinessRatio nSub21(2,1, axes_def, measure_def); NsubjettinessRatio nSub32(3,2, axes_def, measure_def); tau21 = nSub21(my_jet); tau32 = nSub32(my_jet); // Make sure calculations are consistent if (!_testAxes[iA].def().givesRandomizedResults()) { assert(abs(tau21 - tau2/tau1) < epsilon); assert(abs(tau32 - tau3/tau2) < epsilon); } } else { tau21 = tau2/tau1; tau32 = tau3/tau2; } string axesName = _testAxes[iA].short_description(); string left_hashtag; // comment out with # because MultiPass uses random number seed, or because axes do not give reliable results (those at the end of axes vector) if (_testAxes[iA].def().givesRandomizedResults() || iA >= (_testAxes.size() - num_unchecked)) left_hashtag = "#"; else left_hashtag = " "; // Output results: cout << std::right << left_hashtag << setw(23) << axesName << ":" << setw(14) << tau1 << setw(14) << tau2 << setw(14) << tau3 << setw(14) << tau21 << setw(14) << tau32 << endl; } } cout << "-----------------------------------------------------------------------------------------------" << endl; cout << "Done Outputting N-subjettiness Values" << endl; cout << "-----------------------------------------------------------------------------------------------" << endl; //////// // // Finding axes/jets found by N-subjettiness partitioning // // This uses the component_results function to get the subjet information // /////// cout << "-----------------------------------------------------------------------------------------------" << endl; cout << "Outputting N-subjettiness Subjets" << endl; cout << "-----------------------------------------------------------------------------------------------" << endl; // Loop through all options, this time setting up jet finding cout << setprecision(6) << left << fixed; for (unsigned iM = 0; iM < _testMeasures.size(); iM++) { for (unsigned iA = 0; iA < _testRecommendedAxes.size(); iA++) { const PseudoJet & my_jet = antikt_jets[j]; const AxesDefinition & axes_def = _testRecommendedAxes[iA].def(); const MeasureDefinition & measure_def = _testMeasures[iM].def(); // This case doesn't work, so skip it. if (axes_def.givesRandomizedResults()) continue; // define Nsubjettiness functions Nsubjettiness nSub1(1, axes_def, measure_def); Nsubjettiness nSub2(2, axes_def, measure_def); Nsubjettiness nSub3(3, axes_def, measure_def); // get component results TauComponents tau1comp = nSub1.component_result(my_jet); TauComponents tau2comp = nSub2.component_result(my_jet); TauComponents tau3comp = nSub3.component_result(my_jet); vector jets1 = tau1comp.jets(); vector jets2 = tau2comp.jets(); vector jets3 = tau3comp.jets(); vector axes1 = tau1comp.axes(); vector axes2 = tau2comp.axes(); vector axes3 = tau3comp.axes(); cout << "-----------------------------------------------------------------------------------------------" << endl; cout << measure_def.description() << ":" << endl; cout << axes_def.description() << ":" << endl; bool commentOut = false; if (axes_def.givesRandomizedResults()) commentOut = true; // have to comment out min_axes, because it has random values // This helper function tries to find out if the jets have tau information for printing PrintJetsWithComponents(jets1,commentOut); cout << "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -" << endl; PrintJetsWithComponents(jets2,commentOut); cout << "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -" << endl; PrintJetsWithComponents(jets3,commentOut); cout << "^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^" << endl; cout << "Axes Used for Above Subjets" << endl; PrintAxes(axes1,commentOut); cout << "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -" << endl; PrintAxes(axes2,commentOut); cout << "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -" << endl; PrintAxes(axes3,commentOut); } } cout << "-----------------------------------------------------------------------------------------------" << endl; cout << "Done Outputting N-subjettiness Subjets" << endl; cout << "-----------------------------------------------------------------------------------------------" << endl; } ////////// the XCone Jet Algorithm /////////////////////////// //////// // // We define a specific implementation of N-jettiness as a jet algorithm, which we call "XCone". // This is the recommended version for all users. // // Recommended usage of XConePlugin is with beta = 2.0 // Beta = 1.0 is also useful as a recoil-free variant in the face of pile-up. // /////// cout << "-----------------------------------------------------------------------------------------------" << endl; cout << "Using the XCone Jet Algorithm" << endl; cout << "-----------------------------------------------------------------------------------------------" << endl; //create list of various values of beta vector betalist; betalist.push_back(1.0); betalist.push_back(2.0); unsigned int n_betas = betalist.size(); for (unsigned iB = 0; iB < n_betas; iB++) { double beta = betalist[iB]; // define the plugins XConePlugin xcone_plugin2(2, Rcutoff, beta); XConePlugin xcone_plugin3(3, Rcutoff, beta); XConePlugin xcone_plugin4(4, Rcutoff, beta); // and the jet definitions JetDefinition xcone_jetDef2(&xcone_plugin2); JetDefinition xcone_jetDef3(&xcone_plugin3); JetDefinition xcone_jetDef4(&xcone_plugin4); // and the cluster sequences ClusterSequence xcone_seq2(input_particles, xcone_jetDef2); ClusterSequence xcone_seq3(input_particles, xcone_jetDef3); ClusterSequence xcone_seq4(input_particles, xcone_jetDef4); // and associated extras for more information const NjettinessExtras * extras2 = njettiness_extras(xcone_seq2); const NjettinessExtras * extras3 = njettiness_extras(xcone_seq3); const NjettinessExtras * extras4 = njettiness_extras(xcone_seq4); // and find the jets vector xcone_jets2 = xcone_seq2.inclusive_jets(); vector xcone_jets3 = xcone_seq3.inclusive_jets(); vector xcone_jets4 = xcone_seq4.inclusive_jets(); // (alternative way to find the jets) //vector xcone_jets2 = extras2->jets(); //vector xcone_jets3 = extras3->jets(); //vector xcone_jets4 = extras4->jets(); cout << "-----------------------------------------------------------------------------------------------" << endl; cout << "Using beta = " << setprecision(2) << beta << ", Rcut = " << setprecision(2) << Rcutoff << endl; cout << "-----------------------------------------------------------------------------------------------" << endl; PrintJets(xcone_jets2); cout << "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -" << endl; PrintJets(xcone_jets3); cout << "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -" << endl; PrintJets(xcone_jets4); // The axes might point in a different direction than the jets // Using the NjettinessExtras pointer (ClusterSequence::Extras) to access that information vector xcone_axes2 = extras2->axes(); vector xcone_axes3 = extras3->axes(); vector xcone_axes4 = extras4->axes(); cout << "^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^" << endl; cout << "Axes Used for Above Jets" << endl; PrintAxes(xcone_axes2); cout << "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -" << endl; PrintAxes(xcone_axes3); cout << "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -" << endl; PrintAxes(xcone_axes4); bool calculateArea = false; if (calculateArea) { cout << "^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^" << endl; cout << "Adding Area Information (quite slow)" << endl; double ghost_maxrap = 5.0; // e.g. if particles go up to y=5 AreaDefinition area_def(active_area_explicit_ghosts, GhostedAreaSpec(ghost_maxrap)); // Defining cluster sequences with area ClusterSequenceArea xcone_seq_area2(input_particles, xcone_jetDef2, area_def); ClusterSequenceArea xcone_seq_area3(input_particles, xcone_jetDef3, area_def); ClusterSequenceArea xcone_seq_area4(input_particles, xcone_jetDef4, area_def); vector xcone_jets_area2 = xcone_seq_area2.inclusive_jets(); vector xcone_jets_area3 = xcone_seq_area3.inclusive_jets(); vector xcone_jets_area4 = xcone_seq_area4.inclusive_jets(); PrintJets(xcone_jets_area2); cout << "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -" << endl; PrintJets(xcone_jets_area3); cout << "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -" << endl; PrintJets(xcone_jets_area4); } } cout << "-----------------------------------------------------------------------------------------------" << endl; cout << "Done Using the XCone Jet Algorithm" << endl; cout << "-----------------------------------------------------------------------------------------------" << endl; ////////// N-jettiness as a jet algorithm /////////////////////////// //////// // // The user can also defined N-jettiness as a jet algorithm more generally, using different choice // for measures and for axis finding. // // Recommended usage of NjettinessPlugin (event-wide) // AxesMode: wta_kt_axes or onepass_wta_kt_axes // MeasureMode: unnormalized_measure // beta with wta_kt_axes: anything greater than 0.0 (particularly good for 1.0) // beta with onepass_wta_kt_axes: between 1.0 and 3.0 // // Note that the user should find that the usage of Conical Geometric Measure beta = 1.0 with // GenET_GenKT_Axes(std::numeric_limits::max(), 1.0, Rcutoff) should be identical to XCone beta = 1.0, // and Conical Geometric Measure beta = 2.0 with GenET_GenKT_Axes(1.0, 0.5, Rcutoff) should be identical to // XCone beta = 2.0. // /////// cout << "-----------------------------------------------------------------------------------------------" << endl; cout << "Using N-jettiness as a Jet Algorithm" << endl; cout << "-----------------------------------------------------------------------------------------------" << endl; for (unsigned iM = 0; iM < _testCutoffMeasures.size(); iM++) { for (unsigned iA = 0; iA < _testAlgorithmRecommendedAxes.size(); iA++) { const AxesDefinition & axes_def = _testAlgorithmRecommendedAxes[iA].def(); const MeasureDefinition & measure_def = _testCutoffMeasures[iM].def(); // define the plugins NjettinessPlugin njet_plugin2(2, axes_def,measure_def); NjettinessPlugin njet_plugin3(3, axes_def,measure_def); NjettinessPlugin njet_plugin4(4, axes_def,measure_def); // and the jet definitions JetDefinition njet_jetDef2(&njet_plugin2); JetDefinition njet_jetDef3(&njet_plugin3); JetDefinition njet_jetDef4(&njet_plugin4); // and the cluster sequences ClusterSequence njet_seq2(input_particles, njet_jetDef2); ClusterSequence njet_seq3(input_particles, njet_jetDef3); ClusterSequence njet_seq4(input_particles, njet_jetDef4); // and associated extras for more information const NjettinessExtras * extras2 = njettiness_extras(njet_seq2); const NjettinessExtras * extras3 = njettiness_extras(njet_seq3); const NjettinessExtras * extras4 = njettiness_extras(njet_seq4); // and find the jets vector njet_jets2 = njet_seq2.inclusive_jets(); vector njet_jets3 = njet_seq3.inclusive_jets(); vector njet_jets4 = njet_seq4.inclusive_jets(); // (alternative way to find the jets) //vector njet_jets2 = extras2->jets(); //vector njet_jets3 = extras3->jets(); //vector njet_jets4 = extras4->jets(); cout << "-----------------------------------------------------------------------------------------------" << endl; cout << measure_def.description() << ":" << endl; cout << axes_def.description() << ":" << endl; PrintJets(njet_jets2); cout << "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -" << endl; PrintJets(njet_jets3); cout << "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -" << endl; PrintJets(njet_jets4); // The axes might point in a different direction than the jets // Using the NjettinessExtras pointer (ClusterSequence::Extras) to access that information vector njet_axes2 = extras2->axes(); vector njet_axes3 = extras3->axes(); vector njet_axes4 = extras4->axes(); cout << "^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^" << endl; cout << "Axes Used for Above Jets" << endl; PrintAxes(njet_axes2); cout << "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -" << endl; PrintAxes(njet_axes3); cout << "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -" << endl; PrintAxes(njet_axes4); bool calculateArea = false; if (calculateArea) { cout << "^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^" << endl; cout << "Adding Area Information (quite slow)" << endl; double ghost_maxrap = 5.0; // e.g. if particles go up to y=5 AreaDefinition area_def(active_area_explicit_ghosts, GhostedAreaSpec(ghost_maxrap)); // Defining cluster sequences with area ClusterSequenceArea njet_seq_area2(input_particles, njet_jetDef2, area_def); ClusterSequenceArea njet_seq_area3(input_particles, njet_jetDef3, area_def); ClusterSequenceArea njet_seq_area4(input_particles, njet_jetDef4, area_def); vector njet_jets_area2 = njet_seq_area2.inclusive_jets(); vector njet_jets_area3 = njet_seq_area3.inclusive_jets(); vector njet_jets_area4 = njet_seq_area4.inclusive_jets(); PrintJets(njet_jets_area2); cout << "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -" << endl; PrintJets(njet_jets_area3); cout << "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -" << endl; PrintJets(njet_jets_area4); } } } cout << "-----------------------------------------------------------------------------------------------" << endl; cout << "Done Using N-jettiness as a Jet Algorithm" << endl; cout << "-----------------------------------------------------------------------------------------------" << endl; // Below are timing tests for the developers double do_timing_test = false; if (do_timing_test) { clock_t clock_begin, clock_end; double num_iter; cout << setprecision(6); num_iter = 1000; double R0 = 0.5; double beta = 2.0; double N = 6; // AKT JetDefinition aktDef = JetDefinition(antikt_algorithm,R0,E_scheme,Best); // XC XConePlugin xconePlugin(N, R0, beta); JetDefinition xconeDef = JetDefinition(&xconePlugin); // pXC PseudoXConePlugin pseudoxconePlugin(N, R0, beta); JetDefinition pseudoxconeDef = JetDefinition(&pseudoxconePlugin); //AKT cout << "Timing for " << aktDef.description() << endl; clock_begin = clock(); for (int t = 0; t < num_iter; t++) { ClusterSequence clust_seq(input_particles,aktDef); clust_seq.inclusive_jets(); } clock_end = clock(); cout << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per AKT"<< endl; // XC cout << "Timing for " << xconeDef.description() << endl; clock_begin = clock(); for (int t = 0; t < num_iter; t++) { ClusterSequence clust_seq(input_particles,xconeDef); clust_seq.inclusive_jets(); } clock_end = clock(); cout << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per XCone"<< endl; // pXC cout << "Timing for " << pseudoxconePlugin.description() << endl; clock_begin = clock(); for (int t = 0; t < num_iter; t++) { ClusterSequence clust_seq(input_particles,pseudoxconeDef); clust_seq.inclusive_jets(); } clock_end = clock(); cout << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per PseudoXCone"<< endl; } } void PrintJets(const vector & jets, bool commentOut) { string commentStr = ""; if (commentOut) commentStr = "#"; // gets extras information if (jets.size() == 0) return; const NjettinessExtras * extras = njettiness_extras(jets[0]); // bool useExtras = true; bool useExtras = (extras != NULL); bool useArea = jets[0].has_area(); bool useConstit = jets[0].has_constituents(); // define nice tauN header int N = jets.size(); stringstream ss(""); ss << "tau" << N; string tauName = ss.str(); cout << fixed << right; cout << commentStr << setw(5) << "jet #" << " " << setw(10) << "rap" << setw(10) << "phi" << setw(11) << "pt" << setw(11) << "m" << setw(11) << "e"; if (useConstit) cout << setw(11) << "constit"; if (useExtras) cout << setw(14) << tauName; if (useArea) cout << setw(10) << "area"; cout << endl; fastjet::PseudoJet total(0,0,0,0); int total_constit = 0; // print out individual jet information for (unsigned i = 0; i < jets.size(); i++) { cout << commentStr << setw(5) << i+1 << " " << setprecision(4) << setw(10) << jets[i].rap() << setprecision(4) << setw(10) << jets[i].phi() << setprecision(4) << setw(11) << jets[i].perp() << setprecision(4) << setw(11) << max(jets[i].m(),0.0) // needed to fix -0.0 issue on some compilers. << setprecision(4) << setw(11) << jets[i].e(); if (useConstit) cout << setprecision(4) << setw(11) << jets[i].constituents().size(); if (useExtras) cout << setprecision(6) << setw(14) << max(extras->subTau(jets[i]),0.0); if (useArea) cout << setprecision(4) << setw(10) << (jets[i].has_area() ? jets[i].area() : 0.0 ); cout << endl; total += jets[i]; if (useConstit) total_constit += jets[i].constituents().size(); } // print out total jet if (useExtras) { double beamTau = extras->beamTau(); if (beamTau > 0.0) { cout << commentStr << setw(5) << " beam" << " " << setw(10) << "" << setw(10) << "" << setw(11) << "" << setw(11) << "" << setw(11) << "" << setw(11) << "" << setw(14) << setprecision(6) << beamTau << endl; } cout << commentStr << setw(5) << "total" << " " << setprecision(4) << setw(10) << total.rap() << setprecision(4) << setw(10) << total.phi() << setprecision(4) << setw(11) << total.perp() << setprecision(4) << setw(11) << max(total.m(),0.0) // needed to fix -0.0 issue on some compilers. << setprecision(4) << setw(11) << total.e(); if (useConstit) cout << setprecision(4) << setw(11) << total_constit; if (useExtras) cout << setprecision(6) << setw(14) << extras->totalTau(); if (useArea) cout << setprecision(4) << setw(10) << (total.has_area() ? total.area() : 0.0); cout << endl; } } void PrintAxes(const vector & jets, bool commentOut) { string commentStr = ""; if (commentOut) commentStr = "#"; // gets extras information if (jets.size() == 0) return; const NjettinessExtras * extras = njettiness_extras(jets[0]); // bool useExtras = true; bool useExtras = (extras != NULL); bool useArea = jets[0].has_area(); // define nice tauN header int N = jets.size(); stringstream ss(""); ss << "tau" << N; string tauName = ss.str(); cout << fixed << right; cout << commentStr << setw(5) << "jet #" << " " << setw(10) << "rap" << setw(10) << "phi" << setw(11) << "pt" << setw(11) << "m" << setw(11) << "e"; if (useExtras) cout << setw(14) << tauName; if (useArea) cout << setw(10) << "area"; cout << endl; fastjet::PseudoJet total(0,0,0,0); // print out individual jet information for (unsigned i = 0; i < jets.size(); i++) { cout << commentStr << setw(5) << i+1 << " " << setprecision(4) << setw(10) << jets[i].rap() << setprecision(4) << setw(10) << jets[i].phi() << setprecision(4) << setw(11) << jets[i].perp() << setprecision(4) << setw(11) << max(jets[i].m(),0.0) // needed to fix -0.0 issue on some compilers. << setprecision(4) << setw(11) << jets[i].e(); if (useExtras) cout << setprecision(6) << setw(14) << max(extras->subTau(jets[i]),0.0); if (useArea) cout << setprecision(4) << setw(10) << (jets[i].has_area() ? jets[i].area() : 0.0 ); cout << endl; total += jets[i]; } // print out total jet if (useExtras) { double beamTau = extras->beamTau(); if (beamTau > 0.0) { cout << commentStr << setw(5) << " beam" << " " << setw(10) << "" << setw(10) << "" << setw(11) << "" << setw(11) << "" << setw(11) << "" << setw(14) << setprecision(6) << beamTau << endl; } cout << commentStr << setw(5) << "total" << " " << setprecision(4) << setw(10) << total.rap() << setprecision(4) << setw(10) << total.phi() << setprecision(4) << setw(11) << total.perp() << setprecision(4) << setw(11) << max(total.m(),0.0) // needed to fix -0.0 issue on some compilers. << setprecision(4) << setw(11) << total.e() << setprecision(6) << setw(14) << extras->totalTau(); if (useArea) cout << setprecision(4) << setw(10) << (total.has_area() ? total.area() : 0.0); cout << endl; } } void PrintJetsWithComponents(const vector & jets, bool commentOut) { string commentStr = ""; if (commentOut) commentStr = "#"; bool useArea = jets[0].has_area(); // define nice tauN header int N = jets.size(); stringstream ss(""); ss << "tau" << N; string tauName = ss.str(); cout << fixed << right; cout << commentStr << setw(5) << "jet #" << " " << setw(10) << "rap" << setw(10) << "phi" << setw(11) << "pt" << setw(11) << "m" << setw(11) << "e"; if (jets[0].has_constituents()) cout << setw(11) << "constit"; cout << setw(14) << tauName; if (useArea) cout << setw(10) << "area"; cout << endl; fastjet::PseudoJet total(0,0,0,0); double total_tau = 0; int total_constit = 0; // print out individual jet information for (unsigned i = 0; i < jets.size(); i++) { double thisTau = jets[i].structure_of().tau(); cout << commentStr << setw(5) << i+1 << " " << setprecision(4) << setw(10) << jets[i].rap() << setprecision(4) << setw(10) << jets[i].phi() << setprecision(4) << setw(11) << jets[i].perp() << setprecision(4) << setw(11) << max(jets[i].m(),0.0) // needed to fix -0.0 issue on some compilers. << setprecision(4) << setw(11) << jets[i].e(); if (jets[i].has_constituents()) cout << setprecision(4) << setw(11) << jets[i].constituents().size(); cout << setprecision(6) << setw(14) << max(thisTau,0.0); if (useArea) cout << setprecision(4) << setw(10) << (jets[i].has_area() ? jets[i].area() : 0.0 ); cout << endl; total += jets[i]; total_tau += thisTau; if (jets[i].has_constituents()) total_constit += jets[i].constituents().size(); } cout << commentStr << setw(5) << "total" << " " << setprecision(4) << setw(10) << total.rap() << setprecision(4) << setw(10) << total.phi() << setprecision(4) << setw(11) << total.perp() << setprecision(4) << setw(11) << max(total.m(),0.0) // needed to fix -0.0 issue on some compilers. << setprecision(4) << setw(11) << total.e(); if (jets[0].has_constituents()) cout << setprecision(4) << setw(11) << total_constit; cout << setprecision(6) << setw(14) << total_tau; if (useArea) cout << setprecision(4) << setw(10) << (total.has_area() ? total.area() : 0.0); cout << endl; }