Index: contrib/contribs/LundPlane/trunk/example.cc =================================================================== --- contrib/contribs/LundPlane/trunk/example.cc (revision 1288) +++ contrib/contribs/LundPlane/trunk/example.cc (revision 1289) @@ -1,118 +1,119 @@ //---------------------------------------------------------------------- /// \file example.cc /// /// This example program is meant to illustrate how the /// fastjet::contrib::LundGenerator class is used. /// /// Run this example with /// /// \verbatim /// ./example < ../data/single-event.dat /// \endverbatim //---------------------------------------------------------------------- // $Id$ // -// Copyright (c) -, Frederic A. Dreyer, Gavin P. Salam, Gregory Soyez +// Copyright (c) 2018-, Frederic A. Dreyer, Keith Hamilton, Alexander Karlberg, +// Gavin P. Salam, Ludovic Scyboz, Gregory Soyez, Rob Verheyen // //---------------------------------------------------------------------- // 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 "fastjet/PseudoJet.hh" #include #include "LundGenerator.hh" // In external code, this should be fastjet/contrib/LundGenerator.hh #include "LundJSON.hh" // In external code, this should be fastjet/contrib/LundJSON.hh using namespace std; using namespace fastjet; // forward declaration to make things clearer void read_event(vector &event); //---------------------------------------------------------------------- int main(){ //---------------------------------------------------------- // read in input particles vector event; read_event(event); string filename = "jets.json"; cout << "# read an event with " << event.size() << " particles" << endl; cout << "# writing declusterings of primary and secondary plane to file " << filename << endl; ofstream outfile; outfile.open(filename.c_str()); // first get some anti-kt jets double R = 1.0, ptmin = 100.0; JetDefinition jet_def(antikt_algorithm, R); ClusterSequence cs(event, jet_def); vector jets = sorted_by_pt(cs.inclusive_jets(ptmin)); //---------------------------------------------------------- // create an instance of LundGenerator, with default options contrib::LundGenerator lund; cout << lund.description() << endl; for (unsigned ijet = 0; ijet < jets.size(); ijet++) { cout << endl << "Lund coordinates ( ln 1/Delta, ln kt ) of declusterings of jet " << ijet << " are:" << endl; vector declusts = lund(jets[ijet]); for (int idecl = 0; idecl < declusts.size(); idecl++) { pair coords = declusts[idecl].lund_coordinates(); cout << "(" << coords.first << ", " << coords.second << ")"; if (idecl < declusts.size() - 1) cout << "; "; } cout << endl; // outputs the primary Lund plane lund_to_json(outfile, declusts); outfile << endl; // outputs the full Lund tree //to_json(cout, lund_gen, jets[ijet]); cout << endl; } cout << endl << "File " << filename << " written." << endl; 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 is 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); } } Index: contrib/contribs/LundPlane/trunk/example.py =================================================================== --- contrib/contribs/LundPlane/trunk/example.py (revision 1288) +++ contrib/contribs/LundPlane/trunk/example.py (revision 1289) @@ -1,83 +1,84 @@ #!/usr/bin/env python3 # #---------------------------------------------------------------------- # $Id$ # -# Copyright (c) -, Frederic A. Dreyer, Gavin P. Salam, Gregory Soyez +# Copyright (c) 2018-, Frederic A. Dreyer, Keith Hamilton, Alexander Karlberg, +# Gavin P. Salam, Ludovic Scyboz, Gregory Soyez, Rob Verheyen # #---------------------------------------------------------------------- # 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 . #---------------------------------------------------------------------- # # Load a sample file and plot it. # # Usage: # python3 plot_lund.py [--file filename] [--bkg file_bkg] # [--njet njet] [--npxl npixels] # import read_lund_json as lund #from import LundImage from matplotlib.colors import LogNorm import numpy as np import matplotlib.pyplot as plt import argparse parser = argparse.ArgumentParser(description='Plot lund images') parser.add_argument('--file', action = 'store', default = 'jets.json') parser.add_argument('--njet', type = int, default = 2, help='Maximum number of jets to analyse') parser.add_argument('--npxl', type = int, default = 25, help="Number of pixels in each dimension of the image") args = parser.parse_args() # set up the reader and get array from file xval = [0.0, 3.0] yval = [-3.0, 5.0] # start by creating a reader for the json file produced by example.cc # (one json entry per line, correspond to one jet per json entry) reader = lund.Reader(args.file, args.njet) # Then examine the jets it contains print ("Contents of the file", args.file) for jet in reader: # jet is an array of declusterings. # The jet's pt can be obtained by looking at the first declustering (jet[0]) # and extracting the subjet-pair pt ("p_pt") print(" Jet with pt = {:6.1f} GeV with {:3d} primary Lund-plane declusterings".format(jet[0]["p_pt"], len(jet))) print() # Reset the reader to the start and use it with a helper # class to extract an image for each jet reader.reset() image_generator = lund.LundImage(reader, args.njet, args.npxl, xval, yval) images = image_generator.values() # Get the average of the images print("Now creating average lund image from the {} jets".format(len(images))) avg_img = np.average(images,axis=0) # Plot the result fig=plt.figure(figsize=(6, 4.5)) plt.title('Averaged Lund image') plt.xlabel('$\ln(R / \Delta)$') plt.ylabel('$\ln(k_t / \mathrm{GeV})$') plt.imshow(avg_img.transpose(), origin='lower', aspect='auto', extent=xval+yval, cmap=plt.get_cmap('BuPu')) plt.colorbar() print("Close the viewer window to exit") plt.show() Index: contrib/contribs/LundPlane/trunk/example_secondary.cc =================================================================== --- contrib/contribs/LundPlane/trunk/example_secondary.cc (revision 1288) +++ contrib/contribs/LundPlane/trunk/example_secondary.cc (revision 1289) @@ -1,116 +1,117 @@ //---------------------------------------------------------------------- /// \file example_secondary.cc /// /// This example program is meant to illustrate how the /// fastjet::contrib::LundWithSecondary class is used. /// /// Run this example with /// /// \verbatim /// ./example_secondary < ../data/single-event.dat /// \endverbatim //---------------------------------------------------------------------- // $Id$ // -// Copyright (c) -, Frederic A. Dreyer, Gavin P. Salam, Gregory Soyez +// Copyright (c) 2018-, Frederic A. Dreyer, Keith Hamilton, Alexander Karlberg, +// Gavin P. Salam, Ludovic Scyboz, Gregory Soyez, Rob Verheyen // //---------------------------------------------------------------------- // 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 "fastjet/PseudoJet.hh" #include #include "LundWithSecondary.hh" // In external code, this should be fastjet/contrib/LundWithSecondary.hh using namespace std; using namespace fastjet; // forward declaration to make things clearer void read_event(vector &event); //---------------------------------------------------------------------- int main(){ //---------------------------------------------------------- // read in input particles vector event; read_event(event); cout << "# read an event with " << event.size() << " particles" << endl; // first get some anti-kt jets double R = 1.0, ptmin = 100.0; JetDefinition jet_def(antikt_algorithm, R); ClusterSequence cs(event, jet_def); vector jets = sorted_by_pt(cs.inclusive_jets(ptmin)); //---------------------------------------------------------- // create an instance of LundWithSecondary, with default options contrib::SecondaryLund_mMDT secondary; contrib::LundWithSecondary lund(&secondary); cout << lund.description() << endl; for (unsigned ijet = 0; ijet < jets.size(); ijet++) { cout << endl << "Lund coordinates ( ln 1/Delta, ln kt ) of declusterings of jet " << ijet << " are:" << endl; vector declusts = lund.primary(jets[ijet]); for (int idecl = 0; idecl < declusts.size(); idecl++) { pair coords = declusts[idecl].lund_coordinates(); cout << "["<< idecl<< "](" << coords.first << ", " << coords.second << ")"; if (idecl < declusts.size() - 1) cout << "; "; } cout << endl; vector sec_declusts = lund.secondary(declusts); cout << endl << "with Lund coordinates for the secondary plane (from primary declustering [" << lund.secondary_index(declusts) <<"]):" << endl; for (int idecl = 0; idecl < sec_declusts.size(); ++idecl) { pair coords = sec_declusts[idecl].lund_coordinates(); cout << "["<< idecl<< "](" << coords.first << ", " << coords.second << ")"; if (idecl < sec_declusts.size() - 1) cout << "; "; } cout << endl; } 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 is 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); } } Index: contrib/contribs/LundPlane/trunk/LundWithSecondary.cc =================================================================== --- contrib/contribs/LundPlane/trunk/LundWithSecondary.cc (revision 1288) +++ contrib/contribs/LundPlane/trunk/LundWithSecondary.cc (revision 1289) @@ -1,77 +1,78 @@ // $Id$ // -// Copyright (c) -, Frederic A. Dreyer, Gavin P. Salam, Gregory Soyez +// Copyright (c) 2018-, Frederic A. Dreyer, Keith Hamilton, Alexander Karlberg, +// Gavin P. Salam, Ludovic Scyboz, Gregory Soyez, Rob Verheyen // //---------------------------------------------------------------------- // 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 "LundWithSecondary.hh" #include FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh namespace contrib{ //---------------------------------------------------------------------- /// return LundDeclustering sequence of primary plane std::vector LundWithSecondary::primary(const PseudoJet& jet) const { return lund_gen_(jet); } //---------------------------------------------------------------------- /// return LundDeclustering sequence of secondary plane (slow version) std::vector LundWithSecondary::secondary(const PseudoJet& jet) const { // this is not optimal as one is computing the primary plane twice. std::vector declusts = lund_gen_(jet); return secondary(declusts); } //---------------------------------------------------------------------- /// return LundDeclustering sequence of secondary plane with primary sequence as input std::vector LundWithSecondary::secondary( const std::vector & declusts) const { int sec_index = secondary_index(declusts); // if we found the index of secondary emission, return its declustering sequence if (sec_index >= 0) { return lund_gen_(declusts[sec_index].softer()); } else { return std::vector(); } } //---------------------------------------------------------------------- /// return the index associated of the primary declustering that is to be /// used for the secondary plane. int LundWithSecondary::secondary_index(const std::vector & declusts) const { if (secondary_def_ == 0) { throw Error("secondary class is a null pointer, cannot identify element to use for secondary plane"); } return (*secondary_def_)(declusts); } //---------------------------------------------------------------------- /// description std::string LundWithSecondary::description() const { std::ostringstream oss; oss << "LundWithSecondary using " << secondary_def_->description() << " and " << lund_gen_.description(); return oss.str(); } } // namespace contrib FASTJET_END_NAMESPACE Index: contrib/contribs/LundPlane/trunk/SecondaryLund.cc =================================================================== --- contrib/contribs/LundPlane/trunk/SecondaryLund.cc (revision 1288) +++ contrib/contribs/LundPlane/trunk/SecondaryLund.cc (revision 1289) @@ -1,123 +1,124 @@ // $Id$ // -// Copyright (c) 2018-, Frederic A. Dreyer, Gavin P. Salam, Gregory Soyez +// Copyright (c) 2018-, Frederic A. Dreyer, Keith Hamilton, Alexander Karlberg, +// Gavin P. Salam, Ludovic Scyboz, Gregory Soyez, Rob Verheyen // //---------------------------------------------------------------------- // 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 "SecondaryLund.hh" #include #include #include FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh namespace contrib{ //---------------------------------------------------------------------- /// retrieve the vector of declusterings of the secondary plane of a jet int SecondaryLund_mMDT::result(const std::vector& declusts) const { // iterate through primary branchings for (int i=0; i < declusts.size(); ++i) { // mMDTZ: find the first emission passing z>zcut if (declusts[i].z() > zcut_) return i; } return -1; } int SecondaryLund_dotmMDT::result(const std::vector& declusts) const { // set up leading emission pointer and bookkeeping variables int secondary_index = -1; double dot_prod_max = 0.0; // iterate through primary branchings for (int i=0; i < declusts.size(); ++i) { // dotmMDTZ: find emission passing z>zcut with largest dot product if (declusts[i].z() > zcut_) { double dot_prod = declusts[i].harder().pt()*declusts[i].softer().pt() *declusts[i].Delta()*declusts[i].Delta(); if (dot_prod > dot_prod_max) { dot_prod_max = dot_prod; secondary_index = i; } } } // return index of secondary emission return secondary_index; } int SecondaryLund_Mass::result(const std::vector& declusts) const { // set up leading emission pointer and bookkeeping variables int secondary_index = -1; double mass_diff = std::numeric_limits::max(); // iterate through primary branchings for (int i=0; i < declusts.size(); ++i) { // Mass: find emission that minimises the distance to reference mass double dist = std::abs(log(declusts[i].harder().pt()*declusts[i].softer().pt() * declusts[i].Delta()*declusts[i].Delta() / mref2_) * log(1.0/declusts[i].z())); if (dist < mass_diff) { mass_diff = dist; secondary_index = i; } } // return index of secondary emission return secondary_index; } //---------------------------------------------------------------------- /// description std::string SecondaryLund::description() const { std::ostringstream oss; oss << "SecondaryLund"; return oss.str(); } //---------------------------------------------------------------------- /// description std::string SecondaryLund_mMDT::description() const { std::ostringstream oss; oss << "SecondaryLund (mMDT selection of leading emission, zcut=" << zcut_<<")"; return oss.str(); } //---------------------------------------------------------------------- /// description std::string SecondaryLund_dotmMDT::description() const { std::ostringstream oss; oss << "SecondaryLund (dotmMDT selection of leading emission, zcut=" << zcut_<<")"; return oss.str(); } //---------------------------------------------------------------------- /// description std::string SecondaryLund_Mass::description() const { std::ostringstream oss; oss << " (Mass selection of leading emission, m="<. //---------------------------------------------------------------------- #ifndef __FASTJET_CONTRIB_LUNDWITHSECONDARY_HH__ #define __FASTJET_CONTRIB_LUNDWITHSECONDARY_HH__ #include "LundGenerator.hh" #include "SecondaryLund.hh" FASTJET_BEGIN_NAMESPACE namespace contrib{ //---------------------------------------------------------------------- /// \class LundWithSecondary /// Define a primary and secondary Lund plane /// /// \param secondary_def definition used for the leading emission. // class LundWithSecondary { public: /// LundWithSecondary constructor LundWithSecondary(SecondaryLund * secondary_def = 0) : secondary_def_(secondary_def) {} /// LundWithSecondary constructor with jet alg LundWithSecondary(JetAlgorithm jet_alg, SecondaryLund * secondary_def = 0) : lund_gen_(jet_alg), secondary_def_(secondary_def) {} /// LundWithSecondary constructor with jet def LundWithSecondary(const JetDefinition & jet_def, SecondaryLund * secondary_def = 0) : lund_gen_(jet_def), secondary_def_(secondary_def) {} /// destructor virtual ~LundWithSecondary() {} /// primary Lund declustering std::vector primary(const PseudoJet& jet) const; /// secondary Lund declustering (slow) std::vector secondary(const PseudoJet& jet) const; /// secondary Lund declustering with primary sequence as input std::vector secondary( const std::vector & declusts) const; /// return the index associated of the primary declustering that is to be /// used for the secondary plane. int secondary_index(const std::vector & declusts) const; /// description of the class std::string description() const; private: /// lund generator LundGenerator lund_gen_; /// secondary definition SecondaryLund * secondary_def_; }; // //---------------------------------------------------------------------- // /// \class LundWithSecondaryAndTertiary // /// Define a primary, secondary and tertiary Lund plane // class LundWithSecondaryAndTertiary : public LundWithSecondary { // public: // /// LundWithSecondaryAndTertiary constructor // LundWithSecondaryAndTertiary(SecondaryLund * secondary_def = 0, // SecondaryLund * tertiary_def = 0) // : LundWithSecondary(secondary_def), tertiary_def_(tertiary_def) {} // /// LundWithSecondaryAndTertiary constructor with jet alg // LundWithSecondaryAndTertiary(JetAlgorithm jet_alg, // SecondaryLund * secondary_def = 0, // SecondaryLund * tertiary_def = 0) // : LundWithSecondary(jet_alg, secondary_def), tertiary_def_(tertiary_def) {} // /// LundWithSecondaryAndTertiary constructor with jet def // LundWithSecondaryAndTertiary(const JetDefinition & jet_def, // SecondaryLund * secondary_def = 0, // SecondaryLund * tertiary_def = 0) // : LundWithSecondary(jet_def, secondary_def), tertiary_def_(tertiary_def) {} // /// destructor // virtual ~LundWithSecondaryAndTertiary() {} // /// tertiary Lund declustering // virtual std::vector tertiary(const PseudoJet& jet) const; // /// description of the class // virtual std::string description() const; // private: // /// tertiary definition // SecondaryLund * tertiary_def_; // }; } // namespace contrib FASTJET_END_NAMESPACE #endif // __FASTJET_CONTRIB_LUNDWITHSECONDARY_HH__ Index: contrib/contribs/LundPlane/trunk/LundGenerator.hh =================================================================== --- contrib/contribs/LundPlane/trunk/LundGenerator.hh (revision 1288) +++ contrib/contribs/LundPlane/trunk/LundGenerator.hh (revision 1289) @@ -1,137 +1,138 @@ // $Id$ // -// Copyright (c) 2018-, Frederic A. Dreyer, Gavin P. Salam, Gregory Soyez +// Copyright (c) 2018-, Frederic A. Dreyer, Keith Hamilton, Alexander Karlberg, +// Gavin P. Salam, Ludovic Scyboz, Gregory Soyez, Rob Verheyen // //---------------------------------------------------------------------- // This file is part of FastJet contrib. // // It is free software; you can redistribute it and/or modify it under // the terms of the GNU General Public License as published by the // Free Software Foundation; either version 2 of the License, or (at // your option) any later version. // // It is distributed in the hope that it will be useful, but WITHOUT // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY // or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public // License for more details. // // You should have received a copy of the GNU General Public License // along with this code. If not, see . //---------------------------------------------------------------------- #ifndef __FASTJET_CONTRIB_LUNDGENERATOR_HH__ #define __FASTJET_CONTRIB_LUNDGENERATOR_HH__ #include #include "fastjet/tools/Recluster.hh" #include "fastjet/JetDefinition.hh" #include "fastjet/PseudoJet.hh" #include #include #include // TODO: // - add interface to write declusterings to json files // [gps, possibly as a separate header, in order to factorise the json.hh dependence] // - something for pileup subtraction? // - do we want to update json.hh to latest? And handle // the precision issue more elegantly than the current // hack of editing json.hh // - what do we do about the fact that json.hh is c++11? FASTJET_BEGIN_NAMESPACE namespace contrib{ class LundGenerator; //---------------------------------------------------------------------- /// \class LundDeclustering /// Contains the declustering variables associated with a single qnode /// on the Lund plane class LundDeclustering { public: /// return the pair PseudoJet, i.e. sum of the two subjets const PseudoJet & pair() const {return pair_;} /// returns the subjet with larger transverse momentum const PseudoJet & harder() const {return harder_;} /// returns the subjet with smaller transverse momentum const PseudoJet & softer() const {return softer_;} /// returns pair().m() [cached] double m() const {return m_;} /// returns the rapidity-azimuth separation of the pair of subjets [cached] double Delta() const {return Delta_;} /// returns softer().pt() / (softer().pt() + harder().pt()) [cached] double z() const {return z_;} /// returns softer().pt() * Delta() [cached] double kt() const {return kt_;} /// returns z() * Delta() [cached] double kappa() const {return kappa_;} /// returns an azimuthal type angle of softer() around harder() double psi() const {return psi_;} /// returns the x,y coordinates that are used in the Lund-plane plots /// of arXiv:1807.04758: ln(1/Delta()), and ln(kt()) respectively std::pair const lund_coordinates() const { return std::pair(std::log(1.0/Delta()),std::log(kt())); } virtual ~LundDeclustering() {} private: double m_, Delta_, z_, kt_, kappa_, psi_; PseudoJet pair_, harder_, softer_; protected: /// the constructor is private, because users will not generally be /// constructing a LundDeclustering element themselves. LundDeclustering(const PseudoJet& pair, const PseudoJet& j1, const PseudoJet& j2); friend class LundGenerator; }; //---------------------------------------------------------------------- /// \class LundGenerator /// Generates vector of LundDeclustering for a given jet /// corresponding to its Lund plane. class LundGenerator : public FunctionOfPseudoJet< std::vector > { public: /// LundGenerator constructor LundGenerator(JetAlgorithm jet_alg = cambridge_algorithm) : recluster_(JetDefinition(jet_alg, JetDefinition::max_allowable_R)){} /// LundGenerator constructor with jet definition LundGenerator(const JetDefinition & jet_def) : recluster_(jet_def){} /// destructor virtual ~LundGenerator() {} /// obtain the declusterings of the primary plane of the jet virtual std::vector result(const PseudoJet& jet) const; /// description of the class virtual std::string description() const; private: /// recluster definition Recluster recluster_; }; } // namespace contrib FASTJET_END_NAMESPACE #endif // __FASTJET_CONTRIB_LUNDGENERATOR_HH__ Index: contrib/contribs/LundPlane/trunk/SecondaryLund.hh =================================================================== --- contrib/contribs/LundPlane/trunk/SecondaryLund.hh (revision 1288) +++ contrib/contribs/LundPlane/trunk/SecondaryLund.hh (revision 1289) @@ -1,121 +1,122 @@ // $Id$ // -// Copyright (c) 2018-, Frederic A. Dreyer, Gavin P. Salam, Gregory Soyez +// Copyright (c) 2018-, Frederic A. Dreyer, Keith Hamilton, Alexander Karlberg, +// Gavin P. Salam, Ludovic Scyboz, Gregory Soyez, Rob Verheyen // //---------------------------------------------------------------------- // This file is part of FastJet contrib. // // It is free software; you can redistribute it and/or modify it under // the terms of the GNU General Public License as published by the // Free Software Foundation; either version 2 of the License, or (at // your option) any later version. // // It is distributed in the hope that it will be useful, but WITHOUT // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY // or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public // License for more details. // // You should have received a copy of the GNU General Public License // along with this code. If not, see . //---------------------------------------------------------------------- #ifndef __FASTJET_CONTRIB_SECONDARYLUND_HH__ #define __FASTJET_CONTRIB_SECONDARYLUND_HH__ #include "LundGenerator.hh" FASTJET_BEGIN_NAMESPACE namespace contrib{ //---------------------------------------------------------------------- /// \class SecondaryLund /// Base class for definitions of the leading emission class SecondaryLund { public: /// SecondaryLund constructor SecondaryLund() {} /// destructor virtual ~SecondaryLund() {} /// returns the index of branch corresponding to the root of the secondary plane virtual int result(const std::vector & declusts) const = 0; int operator()(const std::vector & declusts) const { return result(declusts); } /// description of the class virtual std::string description() const; }; //---------------------------------------------------------------------- /// \class SecondaryLund_mMDT /// Contains a definition for the leading emission using mMDTZ class SecondaryLund_mMDT : public SecondaryLund { public: /// SecondaryLund_mMDT constructor SecondaryLund_mMDT(double zcut = 0.025) : zcut_(zcut) {} /// destructor virtual ~SecondaryLund_mMDT() {} /// returns the index of branch corresponding to the root of the secondary plane virtual int result(const std::vector & declusts) const; /// description of the class virtual std::string description() const; private: /// zcut parameter double zcut_; }; //---------------------------------------------------------------------- /// \class SecondaryLund_dotmMDT /// Contains a definition for the leading emission using dotmMDT class SecondaryLund_dotmMDT : public SecondaryLund { public: /// SecondaryLund_dotmMDT constructor SecondaryLund_dotmMDT(double zcut = 0.025) : zcut_(zcut) {} /// destructor virtual ~SecondaryLund_dotmMDT() {} /// returns the index of branch corresponding to the root of the secondary plane virtual int result(const std::vector & declusts) const; /// description of the class virtual std::string description() const; private: /// zcut parameter double zcut_; }; //---------------------------------------------------------------------- /// \class SecondaryLund_Mass /// Contains a definition for the leading emission using mass class SecondaryLund_Mass : public SecondaryLund { public: /// SecondaryLund_Mass constructor (default mass reference is W mass) SecondaryLund_Mass(double ref_mass = 80.4) : mref2_(ref_mass*ref_mass) {} /// destructor virtual ~SecondaryLund_Mass() {} /// returns the index of branch corresponding to the root of the secondary plane virtual int result(const std::vector & declusts) const; /// description of the class virtual std::string description() const; private: double mref2_; }; } // namespace contrib FASTJET_END_NAMESPACE #endif // __FASTJET_CONTRIB_SECONDARYLUND_HH__ Index: contrib/contribs/LundPlane/trunk/read_lund_json.py =================================================================== --- contrib/contribs/LundPlane/trunk/read_lund_json.py (revision 1288) +++ contrib/contribs/LundPlane/trunk/read_lund_json.py (revision 1289) @@ -1,224 +1,225 @@ #---------------------------------------------------------------------- # $Id$ # -# Copyright (c) -, Frederic A. Dreyer, Gavin P. Salam, Gregory Soyez +# Copyright (c) 2018-, Frederic A. Dreyer, Keith Hamilton, Alexander Karlberg, +# Gavin P. Salam, Ludovic Scyboz, Gregory Soyez, Rob Verheyen # #---------------------------------------------------------------------- # 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 . #---------------------------------------------------------------------- from abc import ABC, abstractmethod import numpy as np from math import log, ceil, cos, sin import json, gzip, sys #====================================================================== class Reader(object): """ Reader for files consisting of a sequence of json objects, one per line Any pure string object is considered to be part of a header (even if it appears at the end!). Once you have created a Reader, it is a standard Python iterable object. Each iteration gives the json entry for one jet's declustering. """ #---------------------------------------------------------------------- def __init__(self, infile, nmax = -1): self.infile = infile self.nmax = nmax self.reset() #---------------------------------------------------------------------- def set_nmax(self, nmax): self.nmax = nmax #---------------------------------------------------------------------- def reset(self): """ Reset the reader to the start of the file, clear the header and event count. """ if self.infile.endswith('.gz'): self.stream = gzip.open(self.infile,'r') else: self.stream = open(self.infile,'r') self.n = 0 self.header = [] #---------------------------------------------------------------------- def __iter__(self): # needed for iteration to work return self #---------------------------------------------------------------------- def __next__(self): ev = self.next_event() if (ev is None): raise StopIteration else : return ev #---------------------------------------------------------------------- def next(self): return self.__next__() #---------------------------------------------------------------------- def next_event(self): # we have hit the maximum number of events if (self.n == self.nmax): #print ("# Exiting after having read nmax jet declusterings") return None try: line = self.stream.readline() if self.infile.endswith('.gz'): j = json.loads(line.decode('utf-8')) else: j = json.loads(line) except IOError: print("# got to end with IOError (maybe gzip structure broken?) around event", self.n, file=sys.stderr) return None except EOFError: print("# got to end with EOFError (maybe gzip structure broken?) around event", self.n, file=sys.stderr) return None except ValueError: print("# got to end with ValueError (empty json entry?) around event", self.n, file=sys.stderr) return None # skip this if (type(j) is str): self.header.append(j) return self.next_event() self.n += 1 return j #====================================================================== class Image(ABC): """ Abstract base class for batch processing Lund declusterings """ def __init__(self, infile, nmax=-1): if (type(infile) is Reader): self.reader = infile self.reader.set_nmax(nmax) else: self.reader = Reader(infile, nmax) #---------------------------------------------------------------------- @abstractmethod def process(self, event): pass #---------------------------------------------------------------------- def values(self): """Return the values of each event in the reader.""" res = [] while True: event = self.reader.next_event() if event!=None: res.append(self.process(event)) else: break self.reader.reset() return res #====================================================================== class LundImage(Image): """ Class to take input file (or a reader) of jet declusterings in json format, one json entry per line, and transform it into lund images. - infile: a filename or a reader - nmax: the maximum number of jets to process - npxl: the number of bins (pixels) in each dimension - xval: the range of x (ln 1/Delta) values - yval: the range of y (ln kt) values Once you've created the class, call values() (inherited the abstract base class) and you will get a python list of images (formatted as 2d numpy arrays). """ #---------------------------------------------------------------------- def __init__(self, infile, nmax, npxl, xval = [0.0, 7.0], yval = [-3.0, 7.0]): Image.__init__(self, infile, nmax) self.npxl = npxl self.xmin = xval[0] self.ymin = yval[0] self.x_pxl_wdth = (xval[1] - xval[0])/npxl self.y_pxl_wdth = (yval[1] - yval[0])/npxl #---------------------------------------------------------------------- def process(self, event): """Process an event and return an image of the primary Lund plane.""" res = np.zeros((self.npxl,self.npxl)) for declust in event: x = log(1.0/declust['Delta']) y = log(declust['kt']) psi = declust['psi'] xind = ceil((x - self.xmin)/self.x_pxl_wdth - 1.0) yind = ceil((y - self.ymin)/self.y_pxl_wdth - 1.0) # print((x - self.xmin)/self.x_pxl_wdth,xind, # (y - self.ymin)/self.y_pxl_wdth,yind,':', # declust['delta_R'],declust['pt2']) if (max(xind,yind) < self.npxl and min(xind, yind) >= 0): res[xind,yind] += 1 return res #====================================================================== class LundDense(Image): """ Class to take input file (or a reader) of jet declusterings in json format, one json entry per line, and reduces them to the minimal information needed as an input to LSTM or dense network learning. - infile: a filename or a reader - nmax: the maximum number of jets to process - nlen: the size of the output array (for each jet), zero padded if the declustering sequence is shorter; if the declustering sequence is longer, entries beyond nlen are discarded. Calling values() returns a python list, each entry of which is a numpy array of dimension (nlen,2). values()[i,:] = (log(1/Delta),log(kt)) for declustering i. """ #---------------------------------------------------------------------- def __init__(self,infile, nmax, nlen): Image.__init__(self, infile, nmax) self.nlen = nlen #---------------------------------------------------------------------- def process(self, event): """Process an event and return an array of declusterings.""" res = np.zeros((self.nlen, 2)) # go over the declusterings and fill the res array # with the Lund coordinates for i in range(self.nlen): if (i >= len(event)): break res[i,:] = self.fill_declust(event[i]) return res #---------------------------------------------------------------------- def fill_declust(self,declust): """ Create an array of size two and fill it with the Lund coordinates (log(1/Delta),log(kt)). """ res = np.zeros(2) res[0] = log(1.0/declust['Delta']) res[1] = log(declust['kt']) return res