Index: contrib/contribs/SignalFreeBackgroundEstimator/trunk/NEWS
===================================================================
--- contrib/contribs/SignalFreeBackgroundEstimator/trunk/NEWS (revision 0)
+++ contrib/contribs/SignalFreeBackgroundEstimator/trunk/NEWS (revision 1376)
@@ -0,0 +1 @@
+2024/02/01: release of version 1.0.0
Index: contrib/contribs/SignalFreeBackgroundEstimator/trunk/Makefile
===================================================================
--- contrib/contribs/SignalFreeBackgroundEstimator/trunk/Makefile (revision 0)
+++ contrib/contribs/SignalFreeBackgroundEstimator/trunk/Makefile (revision 1376)
@@ -0,0 +1,78 @@
+# If you are using this Makefile standalone and fastjet-config is not
+# in your path, edit this line to specify the full path
+FASTJETCONFIG=fastjet-config
+PREFIX=`$(FASTJETCONFIG) --prefix`
+CXX=g++
+CXXFLAGS= -O3 -Wall -g -fPIC
+install_script = $(SHELL) ../utils/install-sh
+check_script = ../utils/check.sh
+
+# global contrib-wide Makefile include may override some of the above
+# variables (leading "-" means don't give an error if you can't find
+# the file)
+-include ../.Makefile.inc
+
+#------------------------------------------------------------------------
+# things that are specific to this contrib
+NAME=SignalFreeBackgroundEstimator
+SRCS=SignalFreeBackgroundEstimator.cc
+EXAMPLES=example
+INSTALLED_HEADERS=SignalFreeBackgroundEstimator.hh
+#------------------------------------------------------------------------
+
+CXXFLAGS+= $(shell $(FASTJETCONFIG) --cxxflags)
+LDFLAGS += -lm $(shell $(FASTJETCONFIG) --libs)
+
+OBJS = $(SRCS:.cc=.o)
+EXAMPLES_SRCS = $(EXAMPLES:=.cc)
+
+install_HEADER = $(install_script) -c -m 644
+install_LIB = $(install_script) -c -m 644
+install_DIR = $(install_script) -d
+install_DATA = $(install_script) -c -m 644
+install_PROGRAM = $(install_script) -c -s
+install_SCRIPT = $(install_script) -c
+
+.PHONY: clean distclean examples check install
+
+# compilation of the code (default target)
+all: lib$(NAME).a
+
+lib$(NAME).a: $(OBJS)
+ ar cru lib$(NAME).a $(OBJS)
+ ranlib lib$(NAME).a
+
+# building the examples
+examples: $(EXAMPLES)
+
+# the following construct makes it possible to automatically build
+# each of the examples listed in $EXAMPLES
+$(EXAMPLES): % : %.o all
+ $(CXX) $(CXXFLAGS) -o $@ $< -L. -l$(NAME) $(LDFLAGS)
+
+# check that everything went fine
+check: examples
+ @for prog in $(EXAMPLES); do\
+ $(check_script) $${prog} ../data/Pythia-Zp2jets-lhc-pileup-1ev.dat || echo "failed"; \
+ done
+# @echo "All tests successful"
+# $(check_script) $${prog} ../data/Pythia-Zp2jets-lhc-pileup-1ev.dat || exit 1; \
+
+# cleaning the directory
+clean:
+ rm -f *~ *.o
+
+distclean: clean
+ rm -f lib$(NAME).a $(EXAMPLES)
+
+# install things in PREFIX/...
+install: all
+ $(install_DIR) $(PREFIX)/include/fastjet/contrib
+ for header in $(INSTALLED_HEADERS); do\
+ $(install_HEADER) $$header $(PREFIX)/include/fastjet/contrib/;\
+ done
+ $(install_DIR) $(PREFIX)/lib
+ $(install_LIB) lib$(NAME).a $(PREFIX)/lib
+
+depend:
+ makedepend -Y -- -- $(SRCS) $(EXAMPLES_SRCS)
Index: contrib/contribs/SignalFreeBackgroundEstimator/trunk/README
===================================================================
--- contrib/contribs/SignalFreeBackgroundEstimator/trunk/README (revision 0)
+++ contrib/contribs/SignalFreeBackgroundEstimator/trunk/README (revision 1376)
@@ -0,0 +1,19 @@
+SignalFreeBackgroundEstimator is part of the fastjet-contrib project.
+
+The SignalFreeBackgroundEstimator class estimates pileup pT density (rho) with lower dependence on the signal particles, as described in
+
+ Pileup density estimate independent on jet multiplicity
+ Peter Berta, Juraj Smiesko, Martin Spousta
+ arXiv:2304.08383
+
+The code example.cc illustrates how to use the class on individual jets.
+It can be built using
+
+ make example
+
+example.cc makes use of an input datafile which is provided in the
+/data directory, and it should be run with
+
+ ./example < ../data/Pythia-Zp2jets-lhc-pileup-1ev.dat
+
+The expected output can be found in example.ref
\ No newline at end of file
Index: contrib/contribs/SignalFreeBackgroundEstimator/trunk/example.cc
===================================================================
--- contrib/contribs/SignalFreeBackgroundEstimator/trunk/example.cc (revision 0)
+++ contrib/contribs/SignalFreeBackgroundEstimator/trunk/example.cc (revision 1376)
@@ -0,0 +1,204 @@
+//
+//----------------------------------------------------------------------
+// Example on how to estimate background (rho) using the SignalFreeBackgroundEstimator. This background estimate is then used within the Iterative Constituent Subtraction.
+//
+// run it with
+// ./example_signalFreeBackgroundEstimator < ../data/Pythia-Zp2jets-lhc-pileup-1ev.dat
+//----------------------------------------------------------------------
+//
+//----------------------------------------------------------------------
+// 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 "fastjet/ClusterSequence.hh"
+#include "fastjet/Selector.hh"
+#include "SignalFreeBackgroundEstimator.hh"
+#include "fastjet/contrib/IterativeConstituentSubtractor.hh"
+
+#include "functions.hh"
+
+
+using namespace std;
+using namespace fastjet;
+
+//----------------------------------------------------------------------
+int main(){
+ // This should be set up before event loop:
+ double max_eta=4; // specify the maximal pseudorapidity for the input particles. It is used for the subtraction. Particles with eta>|max_eta| are removed and not used during the subtraction (they are not returned). The same parameter should be used for the SignalFreeBackgroundEstimator as it is demonstrated in this example.
+ double max_eta_jet=3; // the maximal pseudorapidity for selected jets. Not important for the subtraction.
+
+ // Signal-free background estimator
+ fastjet::contrib::SignalFreeBackgroundEstimator bge_rho(max_eta,0.55); // maximal pseudo-rapidity cut is used inside ConstituentSubtraction, but in SignalFreeBackgroundEstimator, the range is specified by maximal rapidity cut. Therefore, it is important to apply the same pseudo-rapidity cut also for particles used for background estimation (using function "set_particles") and also derive the rho dependence on rapidity using this max pseudo-rapidity cut to get the correct rescaling function!
+ bge_rho.set_signal_seed_parameters(0.3,0.4); // Set the parameters for signal seeds (the first is the distance parameter for anti-kt clustering of signal particles into jets, and the second is the deltaR distance around signal jets to exclude areas with signal).
+ bge_rho.set_jet_rho_min(5,8); // Set the parameters for the minimal jet rho when selecting jets from the signal jets (intercept and slope of the linear function as a function of some measure of pileup specified in the set_particles function)
+ bge_rho.set_jet_rho_min_charged(20); // If tracking info available, set the minimal jet rho when selecting jets obtained from charged signal particles
+ bge_rho.set_window_parameters(0.5,0.4,0.1); // Set parameters for the window to be used to compute the weighted median
+
+
+ contrib::IterativeConstituentSubtractor subtractor;
+ subtractor.set_max_eta(max_eta); // parameter for maximal |eta| cut. It is specified above.
+ subtractor.set_distance_type(contrib::ConstituentSubtractor::deltaR); // free parameter for the type of distance between particle i and ghost k. There are two options: "deltaR" or "angle" which are defined as deltaR=sqrt((y_i-y_k)^2+(phi_i-phi_k)^2) or Euclidean angle between the momenta
+ vector max_distances;
+ max_distances.push_back(0.15);
+ max_distances.push_back(0.2);
+ vector alphas;
+ alphas.push_back(1);
+ alphas.push_back(1);
+ subtractor.set_parameters(max_distances,alphas); // in this example, 2 CS corrections will be performed: 1. correction with max_distance of 0.1 and alpha of 0, 2. correction with max_distance of 0.15 and alpha of 0. After the first correction, the scalar sum of pt from remaining ghosts is evaluated and redistributed uniformly accross the event.
+ subtractor.set_ghost_removal(true); // set to true if the ghosts (proxies) which were not used in the previous CS procedure should be removed for the next CS procedure
+ subtractor.set_ghost_area(0.004); // parameter for the density of ghosts. The smaller, the better - but also the computation is slower.
+ subtractor.set_background_estimator(&bge_rho); // specify the background estimator to estimate rho. You can use also specify background estimator for the mass term as an additional parameter
+ subtractor.initialize(); // this is new compared to previous usages of ConstituentSubtractor! It should be used after specifying all the parameters and before event loop.
+ // print info (optional)
+ cout << subtractor.description() << endl;
+
+
+ // Grid-median background estimator and alternative ICS method to obtain input to the signal-free background estimator
+ GridMedianBackgroundEstimator bge_rho_grid(max_eta,0.5);
+ contrib::IterativeConstituentSubtractor subtractor_to_get_signal_seeds;
+ subtractor_to_get_signal_seeds.set_parameters(max_distances,alphas); // can be arbitrary
+ subtractor_to_get_signal_seeds.set_ghost_removal(true);
+ subtractor_to_get_signal_seeds.set_max_eta(max_eta);
+ subtractor_to_get_signal_seeds.set_background_estimator(&bge_rho_grid);
+
+
+ // in event loop
+ // read in input particles
+ vector hard_event, full_event;
+ vector *hard_event_charged=new vector;
+ vector *background_event_charged=new vector;
+ read_event(hard_event, full_event, hard_event_charged, background_event_charged);
+
+ hard_event = SelectorAbsEtaMax(max_eta)(hard_event);
+ full_event = SelectorAbsEtaMax(max_eta)(full_event);
+ *hard_event_charged = SelectorAbsEtaMax(max_eta)(*hard_event_charged);
+ *background_event_charged = SelectorAbsEtaMax(max_eta)(*background_event_charged);
+
+ cout << "# read an event with " << hard_event.size() << " signal particles and " << full_event.size() - hard_event.size() << " background particles with pseudo-rapidity |eta|<4" << endl;
+
+ // jet definition and selector for jets
+ JetDefinition jet_def(antikt_algorithm, 0.7);
+ Selector sel_jets = SelectorNHardest(3) * SelectorAbsEtaMax(max_eta_jet);
+
+ // clustering of the hard-scatter event and uncorrected event
+ ClusterSequence clust_seq_hard(hard_event, jet_def);
+ ClusterSequence clust_seq_full(full_event, jet_def);
+ vector hard_jets = sel_jets(clust_seq_hard.inclusive_jets());
+ vector full_jets = sel_jets(clust_seq_full.inclusive_jets());
+
+
+ // There are several options how to define signal seeds for the signal-free background estimator:
+ // 1. Perform an event-wide pileup correction with your preferred method to get estimated signal particles. No tracking information is needed. For example, you can use ICS with grid-median background estimator:
+ //bge_rho_grid.set_particles(full_event);
+ //vector estimated_signal=subtractor_to_get_signal_seeds.subtract_event(full_event);
+ //fastjet::PseudoJet jet_zeroEta(5,0,0,10); // px,py,pz,E
+ //double estimated_rho=bge_rho_grid.rho(jet_zeroEta);
+
+ //bge_rho.set_particles(full_event,estimated_signal,estimated_rho);
+ //vector corrected_event=subtractor.subtract_event(full_event);
+
+ // 2. An additional improvement can be done by doing an additional iteration: The corrected event from option 1 can be used as input to the signal-free background estimator:
+ // bge_rho_grid.set_particles(full_event);
+ // vector estimated_signal=subtractor_to_get_signal_seeds.subtract_event(full_event);
+ // fastjet::PseudoJet jet_zeroEta(5,0,0,10); // px,py,pz,E
+ // double estimated_rho=bge_rho_grid.rho(jet_zeroEta);
+ // bge_rho.set_particles(full_event,estimated_signal,estimated_rho); // first iteration
+ // vector estimated_signal2=subtractor.subtract_event(full_event);
+ // double estimated_rho2=bge_rho.rho(jet_zeroEta);
+
+ // bge_rho.set_particles(full_event,estimated_signal2,estimated_rho2); // second iteration
+ // vector corrected_event=subtractor.subtract_event(full_event);
+
+
+ // 3. If tracking information is available, then the charged signal particles can be used to find the signal seeds, in addition to the method described in 1 or 2.
+ bge_rho_grid.set_particles(full_event); // Alternatively, one can use only the neutral component, i.e. do Charged Hadron Subtraction and use the output here and on the next line.
+ vector estimated_signal=subtractor_to_get_signal_seeds.subtract_event(full_event);
+ fastjet::PseudoJet jet_zeroEta(5,0,0,10); // px,py,pz,E
+ double estimated_rho=bge_rho_grid.rho(jet_zeroEta);
+
+ bge_rho.set_particles(full_event,estimated_signal,estimated_rho,*hard_event_charged);
+ vector corrected_event=subtractor.subtract_event(full_event);
+
+
+ // 4. The user can provide his/her own set of signal seeds:
+
+ // 5. It is also possible to not provide any input for signal seeds. Then no areas will be excluded, and the algortithm is similar to GridMedianBackgroundEstimator, except the fact that the size and the position of the window to compute rho can be specified using the function SignalFreeBackgroundEstimator::set_window_parameters. In this case use:
+ // bge_rho.set_particles(full_event);
+ // vector corrected_event=subtractor.subtract_event(full_event);
+
+
+
+ // clustering of the corrected event
+ ClusterSequence clust_seq_corr(corrected_event, jet_def);
+ vector corrected_jets = sel_jets(clust_seq_corr.inclusive_jets());
+
+ ios::fmtflags f( cout.flags() );
+ cout << setprecision(4) << fixed;
+ cout << endl << "Corrected particles in the whole event:" << endl;
+ for (unsigned int i=0; i{
+public:
+
+ // action of the function
+ double result(const fastjet::PseudoJet &jet) const{
+ if (!jet.has_constituents()){
+ return -0.1;
+ }
+ double width = 1e-6;
+ double ptSum = 0;
+
+ if (jet.constituents().size() < 2) return width;
+
+ std::vector constituents = jet.constituents();
+
+ for (unsigned int iconst=0; iconst &hard_event, vector &full_event, vector *hard_event_charged=0, vector *pileup_event_charged=0){
+ string line;
+ int nsub = 0; // counter to keep track of which sub-event we're reading
+ bool format_ptRapPhi=false;
+ 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,9) == "#PTRAPPHI") format_ptRapPhi=true;
+
+ if (line.substr(0,4) == "#END") {break;}
+ if (line.substr(0,9) == "#SUBSTART") {
+ // if more sub events follow, make copy of first one (the hard one) here
+ if (nsub == 1) hard_event = full_event;
+ nsub += 1;
+ }
+ if (line.substr(0,1) == "#") {continue;}
+ PseudoJet particle(0,0,1,1);
+ double charge,pid;
+
+ if (format_ptRapPhi){
+ double pt,y,phi;
+ linestream >> pt >> y >> phi >> pid >> charge;
+ particle.reset_PtYPhiM(pt,y,phi);
+ }
+ else{
+ double px,py,pz,E;
+ linestream >> px >> py >> pz >> E >> pid >> charge;
+ particle.reset(px,py,pz,E);
+ }
+
+ // push event onto back of full_event vector
+
+ if ( nsub <= 1 ) {
+ if (hard_event_charged && fabs(charge)>0.99) hard_event_charged->push_back(particle);
+ }
+ else {
+ if (pileup_event_charged && fabs(charge)>0.99) pileup_event_charged->push_back(particle);
+ }
+
+ full_event.push_back(particle);
+ }
+
+ // if we have read in only one event, copy it across here...
+ if (nsub == 1) hard_event = full_event;
+
+ // if there was nothing in the event
+ if (nsub == 0) {
+ cerr << "Error: read empty event\n";
+ exit(-1);
+ }
+
+ cout << "# " << nsub-1 << " pileup events on top of the hard event" << endl;
+}