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