Index: contrib/contribs/SignalFreeBackgroundEstimator/trunk/SignalFreeBackgroundEstimator.hh
===================================================================
--- contrib/contribs/SignalFreeBackgroundEstimator/trunk/SignalFreeBackgroundEstimator.hh (revision 1377)
+++ contrib/contribs/SignalFreeBackgroundEstimator/trunk/SignalFreeBackgroundEstimator.hh (revision 1378)
@@ -1,167 +1,170 @@
//
// SignalFreeBackgroundEstimator package
// Questions/comments: peter.berta@cern.ch
//
// Copyright (c) 2024-, Peter Berta
//
//----------------------------------------------------------------------
// 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 __SIGNAL_FREE_BACKGROUND_ESTIMATOR_H__
#define __SIGNAL_FREE_BACKGROUND_ESTIMATOR_H__
#include "fastjet/tools/BackgroundEstimatorBase.hh"
#include "fastjet/RectangularGrid.hh"
#include "fastjet/tools/GridMedianBackgroundEstimator.hh"
#include "fastjet/ClusterSequenceArea.hh"
FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
namespace contrib{
/// \class SignalFreeBackgroundEstimator
/// Background estimator class based on the method described in arXiv:2304.08383. The design of this class is based on the GridMedianBackgroundEstimator class located in the official FastJet software: ./tools/fastjet/tools/GridMedianBackgroundEstimator.hh
class SignalFreeBackgroundEstimator : public fastjet::BackgroundEstimatorBase,
public fastjet::RectangularGrid {
public:
SignalFreeBackgroundEstimator(double rapidity_max,
double requested_grid_spacing): fastjet::RectangularGrid(rapidity_max, requested_grid_spacing), _rapidity_max(rapidity_max) {}
/// Return a pointer to a copy of this estimator
SignalFreeBackgroundEstimator * copy() const override {
return new SignalFreeBackgroundEstimator(*this);
};
/// Tell the background estimator that it has a new event, composed of the specified particles.
void set_particles(const std::vector& particles) override {
std::vector empty;
this->set_particles(particles,empty);
}
/// Tell the background estimator that it has a new event, composed of the specified particles and also include estimated signal particles. If the variable measureOfPileup is set to -1, then the minimal jet rho threshold is constant and set to the parameter _jet_rho_min_intercept. If it is >0, then it is not constant and varies linearly with measureOfPileup. If signalChargedParticles are also provided, then the signal jets are obtained also from clustering of these particles (separate jet clustering is used).
void set_particles(const std::vector& particles, const std::vector& signalParticles, const double& measureOfPileup=-1 /* can be nPU, nPV, mu, rho*/, const std::vector& signalChargedParticles=std::vector());
- /// Get the full set of background properties.
+ /// Add additional signal seeds from the user
+ void add_additional_seeds(const std::vector &additional_seeds);
+
+ /// Get the full set of background properties
fastjet::BackgroundEstimate estimate() const override;
- /// Get the full set of background properties for a given reference jet.
+ /// Get the full set of background properties for a given reference jet
fastjet::BackgroundEstimate estimate(const fastjet::PseudoJet& jet) const override;
- /// Returns rho, the average background density per unit area.
+ /// Returns rho, the average background density per unit area
double rho() const override;
- /// Returns rho, the average background density per unit area, locally at
- /// the position of a given jet.
+ /// Returns rho, the average background density per unit area, locally at the position of a given jet
double rho(const fastjet::PseudoJet & jet) override;
- /// Returns rho_m, the purely longitudinal, particle-mass-induced
- /// component of the background density per unit area
+ /// Returns rho_m, the purely longitudinal, particle-mass-induced component of the background density per unit area
double rho_m() const FASTJET_OVERRIDE;
/// Returns rho_m locally at the jet position. As for rho(jet), it is non-const.
double rho_m(const PseudoJet & jet) FASTJET_OVERRIDE;
/// determine whether the automatic calculation of rho_m (by default true)
void set_compute_rho_m(bool enable){ _enable_rho_m = enable; }
/// Returns true if this background estimator has support for determination of rho_m.
bool has_rho_m() const FASTJET_OVERRIDE {return _enable_rho_m;}
/// Returns a textual description of the background estimator.
std::string description() const override;
/// Set a pointer to a rescaling function
virtual void set_rescaling_class(const fastjet::FunctionOfPseudoJet* rescaling_class) override;
/// Set the parameters for signal seeds (distance parameter for clustering of signal particles into jets, and deltaR distance around signal jets to exclude areas with signal)
void set_signal_seed_parameters(const double &signal_jets_R, const double &exclusion_deltaR){ _signal_jets_R = signal_jets_R; _exclusion_deltaR = exclusion_deltaR; }
/// Set the parameters for the minimal jet rho when selecting jets from the signal jets
void set_jet_rho_min(const double &jet_rho_min_intercept, const double &jet_rho_min_slope){ _jet_rho_min_intercept = jet_rho_min_intercept; _jet_rho_min_slope=jet_rho_min_slope;}
/// Set the parameter for the minimal jet rho when selecting jets obtained from charged signal particles
void set_jet_rho_min_charged(const double &jet_rho_min_charged){ _jet_rho_min_charged = jet_rho_min_charged;}
/// Set the parameter for the spacing of ghosts
void set_ghost_size(const double &ghost_size){ _ghost_size = ghost_size; }
/// Set the parameter for the minimal tile area
void set_tile_area_min(const double &tile_area_min){ _tile_area_min = tile_area_min; }
/// Set parameters for the window to be used to compute the weighted median
void set_window_parameters(const double ¢er, const double &half_window, const double ®ulator_for_floating_center=-1, const double &signal_fraction_min=0.005, const double &shift_max=0.15){ _center = center; _half_window = half_window; _regulator_for_floating_center=regulator_for_floating_center; _signal_fraction_min=signal_fraction_min; _shift_max=shift_max; }
private:
+ // function to compute the weighted median
+ double _compute_weighted_median(std::vector > &tileRho) const;
+
/// Store eta extend.
double _rapidity_max=-1;
/// Distance parameter for jet clustering with anti-kt algorithm to obtain the signal seeds
double _signal_jets_R=0.3;
/// Parameter for the exclusion delta R around the signal seeds
double _exclusion_deltaR=0.4;
/// Parameters to compute the minimal jet rho threshold: intercept and slope of the linear function as a function of some measure of pileup specified in the set_particles function:
double _jet_rho_min_intercept=5;
double _jet_rho_min_slope=8;
/// Parameter for the minimal jet rho threshold for the jets clustered from charged particles. The default value is 20 GeV/(unit area), if the input charged particles are provided in GeV.
double _jet_rho_min_charged=20;
/// Parameter for the minimal tile area after excluding signal-contaminated areas. Tiles with smaller areas are excluded from the weighted median computation. Can be modified with function set_tile_area_min
double _tile_area_min=0.00001;
/// Parameter for the spacing of ghosts to evaluate tile areas after excluding signal-contaminated areas. The smaller ghostSize, the more precise is the computation, but also slower. The default is 0.04. Can be modified with function set_ghost_size
double _ghost_size=0.04;
/// Window parameters:
/// Parameters for the window center and the size of half window used to compute rho
double _center=0.5;
double _half_window=0.1;
/// Parameters for the floating center of the window to compute rho
double _regulator_for_floating_center=-1;
double _signal_fraction_min=0.005;
double _shift_max=0.15;
/// Fraction of scalar pT from signal particles used to move the window to get the median. It is computed based on the provided signal particles.
double _signal_fraction;
/// Sum of all tile areas which are used as weights to compute the weighted median
double _weights_sum;
/// verify that particles have been set and throw an error if not
void verify_particles_set() const;
/// Rescaling warning
fastjet::LimitedWarning _warning_rescaling;
bool _enable_rho_m=true;
-
- double _compute_weighted_median(std::vector > &tileRho) const;
+ std::vector _signal_seeds;
};
}
FASTJET_END_NAMESPACE
#endif // __SIGNAL_FREE_BACKGROUND_ESTIMATOR_H__
Index: contrib/contribs/SignalFreeBackgroundEstimator/trunk/example.ref
===================================================================
--- contrib/contribs/SignalFreeBackgroundEstimator/trunk/example.ref (revision 1377)
+++ contrib/contribs/SignalFreeBackgroundEstimator/trunk/example.ref (revision 1378)
@@ -1,404 +1,16 @@
-
-Description of fastjet::IterativeConstituentSubtractor:
- Using rho estimation: SignalFreeBackgroundEstimator, with rectangular grid with rapidity extent -4 < rap < 4, tile size drap x dphi = 0.533333 x 0.571199
- The masses of all particles will be set to zero.
- The rapidity of the particles will be kept unchanged (not pseudo-rapidity).
- The information about nearby hard proxies will not be used.
- IterativeConstituentSubtractor parameters:
- Iteration 1: max_distance = 0.15 alpha = 1
- Iteration 2: max_distance = 0.2 alpha = 1
-
# 20 pileup events on top of the hard event
# read an event with 203 signal particles and 1714 background particles with pseudo-rapidity |eta|<4
#--------------------------------------------------------------------------
# FastJet release 3.4.2
# M. Cacciari, G.P. Salam and G. Soyez
# A software package for jet finding and analysis at colliders
# http://fastjet.fr
#
# Please cite EPJC72(2012)1896 [arXiv:1111.6097] if you use this package
# for scientific work and optionally PLB641(2006)57 [hep-ph/0512210].
#
# FastJet is provided without warranty under the GNU GPL v2 or higher.
# It uses T. Chan's closest pair algorithm, S. Fortune's Voronoi code
# and 3rd party plugin jet algorithms. See COPYING file for details.
#--------------------------------------------------------------------------
-
-Corrected particles in the whole event:
-pt = 0.1207, phi = 4.5820, rap = -3.9140, |mass| = 0.0000
-pt = 1.1780, phi = 5.2599, rap = -3.6581, |mass| = 0.0000
-pt = 0.5552, phi = 4.9980, rap = -3.6347, |mass| = 0.0000
-pt = 0.3851, phi = 4.2556, rap = -3.5526, |mass| = 0.0000
-pt = 0.6025, phi = 5.6861, rap = -3.5493, |mass| = 0.0000
-pt = 0.6117, phi = 2.6409, rap = -3.4932, |mass| = 0.0000
-pt = 0.7114, phi = 5.5855, rap = -3.4828, |mass| = 0.0000
-pt = 2.2578, phi = 6.0415, rap = -3.4827, |mass| = 0.0000
-pt = 0.7888, phi = 4.9418, rap = -3.4609, |mass| = 0.0000
-pt = 0.7376, phi = 5.4112, rap = -3.4232, |mass| = 0.0000
-pt = 0.3035, phi = 6.0661, rap = -3.4124, |mass| = 0.0000
-pt = 0.1848, phi = 3.2413, rap = -3.3861, |mass| = 0.0000
-pt = 0.6880, phi = 0.2742, rap = -3.3640, |mass| = 0.0000
-pt = 0.0913, phi = 2.2892, rap = -3.3465, |mass| = 0.0000
-pt = 1.6549, phi = 5.8966, rap = -3.3435, |mass| = 0.0000
-pt = 0.3217, phi = 1.6309, rap = -3.3372, |mass| = 0.0000
-pt = 1.2491, phi = 2.2916, rap = -3.3330, |mass| = 0.0000
-pt = 0.5767, phi = 0.2181, rap = -3.3282, |mass| = 0.0000
-pt = 0.5802, phi = 2.6996, rap = -3.3022, |mass| = 0.0000
-pt = 3.1930, phi = 3.0096, rap = -3.2831, |mass| = 0.0000
-pt = 0.0436, phi = 2.3489, rap = -3.2400, |mass| = 0.0000
-pt = 0.2150, phi = 5.4556, rap = -3.2297, |mass| = 0.0000
-pt = 1.1087, phi = 1.4148, rap = -3.2284, |mass| = 0.0000
-pt = 0.1610, phi = 4.5544, rap = -3.1881, |mass| = 0.0000
-pt = 0.2646, phi = 5.3974, rap = -3.1863, |mass| = 0.0000
-pt = 1.3092, phi = 5.9358, rap = -3.1711, |mass| = 0.0000
-pt = 0.4090, phi = 3.1890, rap = -3.1614, |mass| = 0.0000
-pt = 1.4384, phi = 5.7996, rap = -3.1614, |mass| = 0.0000
-pt = 2.1327, phi = 4.4399, rap = -3.1188, |mass| = 0.0000
-pt = 0.3384, phi = 3.0074, rap = -3.1125, |mass| = 0.0000
-pt = 3.3494, phi = 3.1301, rap = -3.0963, |mass| = 0.0000
-pt = 0.6148, phi = 4.5479, rap = -3.0839, |mass| = 0.0000
-pt = 0.4756, phi = 2.7204, rap = -3.0545, |mass| = 0.0000
-pt = 1.4137, phi = 5.6871, rap = -3.0433, |mass| = 0.0000
-pt = 0.3270, phi = 4.5763, rap = -3.0201, |mass| = 0.0000
-pt = 0.3096, phi = 4.1676, rap = -3.0161, |mass| = 0.0000
-pt = 0.8174, phi = 5.6925, rap = -3.0015, |mass| = 0.0000
-pt = 0.5059, phi = 0.7339, rap = -2.9925, |mass| = 0.0000
-pt = 0.6720, phi = 5.2605, rap = -2.9511, |mass| = 0.0000
-pt = 1.2715, phi = 5.4574, rap = -2.9496, |mass| = 0.0000
-pt = 0.4487, phi = 5.7954, rap = -2.8763, |mass| = 0.0000
-pt = 0.3482, phi = 5.1832, rap = -2.8757, |mass| = 0.0000
-pt = 0.3969, phi = 4.9577, rap = -2.8737, |mass| = 0.0000
-pt = 2.5949, phi = 5.0940, rap = -2.8440, |mass| = 0.0000
-pt = 0.7022, phi = 5.1006, rap = -2.8366, |mass| = 0.0000
-pt = 0.1212, phi = 0.2435, rap = -2.8043, |mass| = 0.0000
-pt = 0.5727, phi = 5.0895, rap = -2.7945, |mass| = 0.0000
-pt = 0.2731, phi = 5.7579, rap = -2.7939, |mass| = 0.0000
-pt = 0.2301, phi = 5.2141, rap = -2.7867, |mass| = 0.0000
-pt = 1.2636, phi = 3.2932, rap = -2.7708, |mass| = 0.0000
-pt = 0.0675, phi = 4.3068, rap = -2.7459, |mass| = 0.0000
-pt = 0.1727, phi = 5.0971, rap = -2.7245, |mass| = 0.0000
-pt = 1.3920, phi = 4.3266, rap = -2.7104, |mass| = 0.0000
-pt = 0.2323, phi = 4.9647, rap = -2.6730, |mass| = 0.0000
-pt = 0.3236, phi = 0.7724, rap = -2.6691, |mass| = 0.0000
-pt = 1.2201, phi = 3.7925, rap = -2.6181, |mass| = 0.0000
-pt = 0.1799, phi = 4.2586, rap = -2.6133, |mass| = 0.0000
-pt = 0.7033, phi = 5.2473, rap = -2.5999, |mass| = 0.0000
-pt = 0.5591, phi = 0.8228, rap = -2.5632, |mass| = 0.0000
-pt = 0.5210, phi = 5.4743, rap = -2.5477, |mass| = 0.0000
-pt = 1.4704, phi = 3.7936, rap = -2.5437, |mass| = 0.0000
-pt = 0.7513, phi = 0.7076, rap = -2.5349, |mass| = 0.0000
-pt = 0.0043, phi = 0.2840, rap = -2.5265, |mass| = 0.0000
-pt = 0.4688, phi = 3.8120, rap = -2.5092, |mass| = 0.0000
-pt = 0.7452, phi = 5.7415, rap = -2.5082, |mass| = 0.0000
-pt = 0.0776, phi = 5.7518, rap = -2.4951, |mass| = 0.0000
-pt = 0.2097, phi = 5.5739, rap = -2.4928, |mass| = 0.0000
-pt = 0.5182, phi = 5.7067, rap = -2.4847, |mass| = 0.0000
-pt = 0.6459, phi = 2.2857, rap = -2.4781, |mass| = 0.0000
-pt = 0.0002, phi = 4.4585, rap = -2.4677, |mass| = 0.0000
-pt = 0.7019, phi = 2.3820, rap = -2.4577, |mass| = 0.0000
-pt = 0.0115, phi = 4.7296, rap = -2.4105, |mass| = 0.0000
-pt = 0.0417, phi = 0.2745, rap = -2.4075, |mass| = 0.0000
-pt = 0.4900, phi = 2.2732, rap = -2.3606, |mass| = 0.0000
-pt = 1.0378, phi = 5.6589, rap = -2.3513, |mass| = 0.0000
-pt = 0.3025, phi = 2.6826, rap = -2.3242, |mass| = 0.0000
-pt = 0.3434, phi = 4.2981, rap = -2.3130, |mass| = 0.0000
-pt = 0.3171, phi = 2.3724, rap = -2.3003, |mass| = 0.0000
-pt = 0.3051, phi = 4.0933, rap = -2.2952, |mass| = 0.0000
-pt = 0.3025, phi = 4.0379, rap = -2.2709, |mass| = 0.0000
-pt = 0.6777, phi = 6.1109, rap = -2.2503, |mass| = 0.0000
-pt = 0.3981, phi = 2.0559, rap = -2.2459, |mass| = 0.0000
-pt = 0.0453, phi = 5.8451, rap = -2.1326, |mass| = 0.0000
-pt = 0.1202, phi = 5.2157, rap = -2.0325, |mass| = 0.0000
-pt = 0.0039, phi = 2.5354, rap = -1.9378, |mass| = 0.0000
-pt = 1.1229, phi = 1.7065, rap = -1.9364, |mass| = 0.0000
-pt = 0.3434, phi = 4.6298, rap = -1.8117, |mass| = 0.0000
-pt = 0.5442, phi = 4.1310, rap = -1.7928, |mass| = 0.0000
-pt = 0.7649, phi = 4.0208, rap = -1.6984, |mass| = 0.0000
-pt = 0.1120, phi = 1.6670, rap = -1.6500, |mass| = 0.0000
-pt = 0.3269, phi = 4.3785, rap = -1.5446, |mass| = 0.0000
-pt = 0.4704, phi = 4.4554, rap = -1.4610, |mass| = 0.0000
-pt = 0.6879, phi = 0.6859, rap = -1.3568, |mass| = 0.0000
-pt = 0.1825, phi = 1.3556, rap = -1.3363, |mass| = 0.0000
-pt = 0.9195, phi = 0.0872, rap = -1.1921, |mass| = 0.0000
-pt = 24.4459, phi = 0.0579, rap = -1.0933, |mass| = 0.0000
-pt = 10.6629, phi = 0.0731, rap = -1.0859, |mass| = 0.0000
-pt = 22.7084, phi = 0.0773, rap = -1.0854, |mass| = 0.0000
-pt = 20.5108, phi = 0.0834, rap = -1.0851, |mass| = 0.0000
-pt = 70.5496, phi = 0.0743, rap = -1.0811, |mass| = 0.0000
-pt = 12.6654, phi = 0.0757, rap = -1.0660, |mass| = 0.0000
-pt = 8.1149, phi = 0.0626, rap = -1.0644, |mass| = 0.0000
-pt = 4.3234, phi = 0.0345, rap = -1.0580, |mass| = 0.0000
-pt = 1.1803, phi = 0.0917, rap = -1.0549, |mass| = 0.0000
-pt = 1.2027, phi = 0.0712, rap = -1.0548, |mass| = 0.0000
-pt = 6.0000, phi = 0.0565, rap = -1.0434, |mass| = 0.0000
-pt = 4.5794, phi = 0.0586, rap = -1.0426, |mass| = 0.0000
-pt = 12.7933, phi = 0.1181, rap = -1.0377, |mass| = 0.0000
-pt = 4.0644, phi = 6.2190, rap = -1.0143, |mass| = 0.0000
-pt = 4.5256, phi = 0.0865, rap = -1.0003, |mass| = 0.0000
-pt = 2.0608, phi = 0.0711, rap = -0.9980, |mass| = 0.0000
-pt = 0.2533, phi = 0.9472, rap = -0.9624, |mass| = 0.0000
-pt = 3.5843, phi = 6.2357, rap = -0.9621, |mass| = 0.0000
-pt = 0.5732, phi = 4.0146, rap = -0.9112, |mass| = 0.0000
-pt = 0.4447, phi = 4.3827, rap = -0.8746, |mass| = 0.0000
-pt = 0.0692, phi = 0.8457, rap = -0.8653, |mass| = 0.0000
-pt = 0.5969, phi = 0.1357, rap = -0.8338, |mass| = 0.0000
-pt = 0.3044, phi = 2.7181, rap = -0.7632, |mass| = 0.0000
-pt = 0.3242, phi = 1.5393, rap = -0.7199, |mass| = 0.0000
-pt = 0.0034, phi = 1.6905, rap = -0.6892, |mass| = 0.0000
-pt = 0.2333, phi = 0.6640, rap = -0.6877, |mass| = 0.0000
-pt = 0.5739, phi = 2.7094, rap = -0.6874, |mass| = 0.0000
-pt = 0.1001, phi = 1.5224, rap = -0.6862, |mass| = 0.0000
-pt = 0.2955, phi = 0.0850, rap = -0.6858, |mass| = 0.0000
-pt = 0.1202, phi = 2.6770, rap = -0.6664, |mass| = 0.0000
-pt = 0.1183, phi = 3.5134, rap = -0.6446, |mass| = 0.0000
-pt = 0.1013, phi = 2.8407, rap = -0.6369, |mass| = 0.0000
-pt = 0.2240, phi = 5.5749, rap = -0.6292, |mass| = 0.0000
-pt = 0.7548, phi = 1.5600, rap = -0.5723, |mass| = 0.0000
-pt = 0.8892, phi = 0.7584, rap = -0.5473, |mass| = 0.0000
-pt = 0.2254, phi = 1.3913, rap = -0.5288, |mass| = 0.0000
-pt = 1.7859, phi = 2.3648, rap = -0.4499, |mass| = 0.0000
-pt = 0.0341, phi = 2.9395, rap = -0.4264, |mass| = 0.0000
-pt = 1.2646, phi = 0.4795, rap = -0.3768, |mass| = 0.0000
-pt = 1.2770, phi = 0.3563, rap = -0.3342, |mass| = 0.0000
-pt = 0.0114, phi = 0.5804, rap = -0.2940, |mass| = 0.0000
-pt = 0.9964, phi = 2.7340, rap = -0.2597, |mass| = 0.0000
-pt = 0.1239, phi = 1.7579, rap = -0.2434, |mass| = 0.0000
-pt = 0.8876, phi = 1.8261, rap = -0.2425, |mass| = 0.0000
-pt = 0.2459, phi = 1.2481, rap = -0.2241, |mass| = 0.0000
-pt = 0.6269, phi = 0.5956, rap = -0.1828, |mass| = 0.0000
-pt = 0.0134, phi = 2.2445, rap = -0.1654, |mass| = 0.0000
-pt = 0.0417, phi = 0.6866, rap = -0.1372, |mass| = 0.0000
-pt = 0.5334, phi = 3.7742, rap = -0.1045, |mass| = 0.0000
-pt = 0.1726, phi = 2.9150, rap = -0.0938, |mass| = 0.0000
-pt = 0.2286, phi = 2.8847, rap = -0.0887, |mass| = 0.0000
-pt = 0.2015, phi = 3.6537, rap = -0.0631, |mass| = 0.0000
-pt = 1.5339, phi = 2.9835, rap = -0.0366, |mass| = 0.0000
-pt = 0.9162, phi = 2.5866, rap = -0.0260, |mass| = 0.0000
-pt = 0.1328, phi = 1.2679, rap = -0.0125, |mass| = 0.0000
-pt = 0.0909, phi = 3.0993, rap = -0.0100, |mass| = 0.0000
-pt = 0.8049, phi = 2.3563, rap = 0.0151, |mass| = 0.0000
-pt = 0.7226, phi = 0.9729, rap = 0.0474, |mass| = 0.0000
-pt = 1.8253, phi = 2.9511, rap = 0.0737, |mass| = 0.0000
-pt = 1.2129, phi = 0.9594, rap = 0.0939, |mass| = 0.0000
-pt = 1.3432, phi = 2.9165, rap = 0.1176, |mass| = 0.0000
-pt = 1.5255, phi = 2.9748, rap = 0.1206, |mass| = 0.0000
-pt = 0.3297, phi = 2.1295, rap = 0.1566, |mass| = 0.0000
-pt = 0.1920, phi = 4.5570, rap = 0.1722, |mass| = 0.0000
-pt = 0.2387, phi = 5.7815, rap = 0.1729, |mass| = 0.0000
-pt = 0.2428, phi = 4.8472, rap = 0.1892, |mass| = 0.0000
-pt = 0.6449, phi = 4.4850, rap = 0.1935, |mass| = 0.0000
-pt = 0.2817, phi = 0.2256, rap = 0.2082, |mass| = 0.0000
-pt = 0.1058, phi = 1.7689, rap = 0.2374, |mass| = 0.0000
-pt = 0.2601, phi = 0.0974, rap = 0.2805, |mass| = 0.0000
-pt = 0.4149, phi = 0.2089, rap = 0.2861, |mass| = 0.0000
-pt = 0.0644, phi = 0.3520, rap = 0.2948, |mass| = 0.0000
-pt = 0.8727, phi = 1.8427, rap = 0.2980, |mass| = 0.0000
-pt = 0.4952, phi = 0.8566, rap = 0.3299, |mass| = 0.0000
-pt = 0.7580, phi = 3.4631, rap = 0.3500, |mass| = 0.0000
-pt = 0.0390, phi = 0.8614, rap = 0.3626, |mass| = 0.0000
-pt = 0.2393, phi = 0.6328, rap = 0.4070, |mass| = 0.0000
-pt = 0.4407, phi = 0.9692, rap = 0.4178, |mass| = 0.0000
-pt = 1.7884, phi = 0.7335, rap = 0.4255, |mass| = 0.0000
-pt = 0.1680, phi = 5.4320, rap = 0.4910, |mass| = 0.0000
-pt = 0.2032, phi = 4.5484, rap = 0.4993, |mass| = 0.0000
-pt = 0.5485, phi = 2.9835, rap = 0.5077, |mass| = 0.0000
-pt = 0.2371, phi = 0.6706, rap = 0.5158, |mass| = 0.0000
-pt = 2.2708, phi = 0.8053, rap = 0.5295, |mass| = 0.0000
-pt = 0.7310, phi = 2.2633, rap = 0.5336, |mass| = 0.0000
-pt = 4.1556, phi = 3.4258, rap = 0.5413, |mass| = 0.0000
-pt = 1.4952, phi = 3.3710, rap = 0.5773, |mass| = 0.0000
-pt = 4.5495, phi = 3.4551, rap = 0.5839, |mass| = 0.0000
-pt = 1.1263, phi = 1.2803, rap = 0.5848, |mass| = 0.0000
-pt = 3.4646, phi = 3.4360, rap = 0.5910, |mass| = 0.0000
-pt = 2.0260, phi = 3.5194, rap = 0.5915, |mass| = 0.0000
-pt = 1.0205, phi = 3.6197, rap = 0.5925, |mass| = 0.0000
-pt = 2.3729, phi = 3.5442, rap = 0.6111, |mass| = 0.0000
-pt = 5.1586, phi = 3.4756, rap = 0.6166, |mass| = 0.0000
-pt = 0.6532, phi = 5.4022, rap = 0.6202, |mass| = 0.0000
-pt = 4.3862, phi = 3.4572, rap = 0.6241, |mass| = 0.0000
-pt = 0.3941, phi = 0.6815, rap = 0.6313, |mass| = 0.0000
-pt = 0.8450, phi = 2.1577, rap = 0.6338, |mass| = 0.0000
-pt = 8.3256, phi = 3.4533, rap = 0.6359, |mass| = 0.0000
-pt = 1.8902, phi = 3.4727, rap = 0.6368, |mass| = 0.0000
-pt = 10.1706, phi = 3.4104, rap = 0.6385, |mass| = 0.0000
-pt = 1.9006, phi = 3.4744, rap = 0.6420, |mass| = 0.0000
-pt = 1.8977, phi = 3.4329, rap = 0.6449, |mass| = 0.0000
-pt = 1.9560, phi = 3.3996, rap = 0.6533, |mass| = 0.0000
-pt = 0.5467, phi = 3.5051, rap = 0.6559, |mass| = 0.0000
-pt = 1.6145, phi = 3.4793, rap = 0.6577, |mass| = 0.0000
-pt = 6.7198, phi = 3.4263, rap = 0.6685, |mass| = 0.0000
-pt = 10.7997, phi = 3.4717, rap = 0.6706, |mass| = 0.0000
-pt = 0.9452, phi = 3.4751, rap = 0.6708, |mass| = 0.0000
-pt = 2.2925, phi = 3.3178, rap = 0.6756, |mass| = 0.0000
-pt = 5.6895, phi = 3.4260, rap = 0.6768, |mass| = 0.0000
-pt = 0.4761, phi = 6.0878, rap = 0.6805, |mass| = 0.0000
-pt = 1.8819, phi = 3.4712, rap = 0.6823, |mass| = 0.0000
-pt = 5.5091, phi = 3.4245, rap = 0.6858, |mass| = 0.0000
-pt = 1.8233, phi = 3.4087, rap = 0.6874, |mass| = 0.0000
-pt = 2.8740, phi = 3.4624, rap = 0.6911, |mass| = 0.0000
-pt = 13.7509, phi = 3.4120, rap = 0.7074, |mass| = 0.0000
-pt = 8.8154, phi = 3.4972, rap = 0.7090, |mass| = 0.0000
-pt = 5.3372, phi = 3.4490, rap = 0.7170, |mass| = 0.0000
-pt = 6.5059, phi = 3.4530, rap = 0.7207, |mass| = 0.0000
-pt = 0.7166, phi = 2.4215, rap = 0.7229, |mass| = 0.0000
-pt = 0.5389, phi = 4.0178, rap = 0.7319, |mass| = 0.0000
-pt = 0.2500, phi = 2.5919, rap = 0.7321, |mass| = 0.0000
-pt = 9.3849, phi = 3.4142, rap = 0.7452, |mass| = 0.0000
-pt = 0.0304, phi = 2.0804, rap = 0.7527, |mass| = 0.0000
-pt = 0.6741, phi = 1.9752, rap = 0.7563, |mass| = 0.0000
-pt = 1.1696, phi = 1.2971, rap = 0.7623, |mass| = 0.0000
-pt = 0.8436, phi = 3.7466, rap = 0.7690, |mass| = 0.0000
-pt = 0.5235, phi = 1.4775, rap = 0.8066, |mass| = 0.0000
-pt = 0.3310, phi = 3.0398, rap = 0.8188, |mass| = 0.0000
-pt = 1.2395, phi = 3.6233, rap = 0.8190, |mass| = 0.0000
-pt = 1.8261, phi = 2.4207, rap = 0.8200, |mass| = 0.0000
-pt = 0.5166, phi = 5.9957, rap = 0.8491, |mass| = 0.0000
-pt = 0.9150, phi = 6.0737, rap = 0.8541, |mass| = 0.0000
-pt = 0.2201, phi = 0.4515, rap = 0.8594, |mass| = 0.0000
-pt = 1.5799, phi = 5.8830, rap = 0.8746, |mass| = 0.0000
-pt = 0.6406, phi = 1.4696, rap = 0.8781, |mass| = 0.0000
-pt = 4.0321, phi = 5.8902, rap = 0.8995, |mass| = 0.0000
-pt = 0.8265, phi = 5.9910, rap = 0.9133, |mass| = 0.0000
-pt = 0.3163, phi = 4.6516, rap = 0.9169, |mass| = 0.0000
-pt = 0.1635, phi = 0.3048, rap = 0.9385, |mass| = 0.0000
-pt = 1.8945, phi = 5.9042, rap = 0.9416, |mass| = 0.0000
-pt = 0.1435, phi = 3.0645, rap = 0.9518, |mass| = 0.0000
-pt = 0.2727, phi = 2.6248, rap = 0.9718, |mass| = 0.0000
-pt = 0.0373, phi = 0.7729, rap = 0.9754, |mass| = 0.0000
-pt = 0.4325, phi = 6.2096, rap = 0.9917, |mass| = 0.0000
-pt = 0.8909, phi = 1.8193, rap = 0.9987, |mass| = 0.0000
-pt = 1.8756, phi = 1.9405, rap = 1.0063, |mass| = 0.0000
-pt = 0.3476, phi = 3.2114, rap = 1.0078, |mass| = 0.0000
-pt = 0.1939, phi = 3.8936, rap = 1.0109, |mass| = 0.0000
-pt = 0.4405, phi = 6.0103, rap = 1.0350, |mass| = 0.0000
-pt = 0.1149, phi = 1.5879, rap = 1.0489, |mass| = 0.0000
-pt = 0.7168, phi = 2.4110, rap = 1.0565, |mass| = 0.0000
-pt = 1.7619, phi = 2.1993, rap = 1.0685, |mass| = 0.0000
-pt = 0.7826, phi = 2.4491, rap = 1.0845, |mass| = 0.0000
-pt = 0.1973, phi = 4.5479, rap = 1.0946, |mass| = 0.0000
-pt = 0.0339, phi = 0.0435, rap = 1.0969, |mass| = 0.0000
-pt = 1.2336, phi = 3.7716, rap = 1.1266, |mass| = 0.0000
-pt = 0.7577, phi = 3.8805, rap = 1.1281, |mass| = 0.0000
-pt = 1.2612, phi = 2.1470, rap = 1.1299, |mass| = 0.0000
-pt = 0.2082, phi = 6.2627, rap = 1.1410, |mass| = 0.0000
-pt = 0.3986, phi = 3.6682, rap = 1.1423, |mass| = 0.0000
-pt = 0.3208, phi = 6.2192, rap = 1.1465, |mass| = 0.0000
-pt = 1.5281, phi = 2.2491, rap = 1.1824, |mass| = 0.0000
-pt = 0.5662, phi = 1.5378, rap = 1.2240, |mass| = 0.0000
-pt = 2.1756, phi = 2.0002, rap = 1.2478, |mass| = 0.0000
-pt = 0.0473, phi = 2.1465, rap = 1.2830, |mass| = 0.0000
-pt = 1.3244, phi = 2.4614, rap = 1.2917, |mass| = 0.0000
-pt = 0.8752, phi = 1.3460, rap = 1.3049, |mass| = 0.0000
-pt = 1.0254, phi = 2.1008, rap = 1.3221, |mass| = 0.0000
-pt = 0.7200, phi = 3.4732, rap = 1.3242, |mass| = 0.0000
-pt = 0.5594, phi = 0.4431, rap = 1.3560, |mass| = 0.0000
-pt = 0.1076, phi = 4.2127, rap = 1.3914, |mass| = 0.0000
-pt = 0.1558, phi = 3.5671, rap = 1.3936, |mass| = 0.0000
-pt = 0.6409, phi = 4.9412, rap = 1.4212, |mass| = 0.0000
-pt = 1.6360, phi = 2.6652, rap = 1.4574, |mass| = 0.0000
-pt = 0.6448, phi = 3.8357, rap = 1.4605, |mass| = 0.0000
-pt = 1.5019, phi = 4.7690, rap = 1.4731, |mass| = 0.0000
-pt = 0.7338, phi = 2.2565, rap = 1.5072, |mass| = 0.0000
-pt = 1.0588, phi = 2.6610, rap = 1.5637, |mass| = 0.0000
-pt = 1.6028, phi = 5.5777, rap = 1.5689, |mass| = 0.0000
-pt = 0.1342, phi = 4.1969, rap = 1.5827, |mass| = 0.0000
-pt = 0.2035, phi = 1.8101, rap = 1.6003, |mass| = 0.0000
-pt = 0.6901, phi = 0.9558, rap = 1.6013, |mass| = 0.0000
-pt = 2.0083, phi = 2.7272, rap = 1.6053, |mass| = 0.0000
-pt = 0.1648, phi = 5.0173, rap = 1.6058, |mass| = 0.0000
-pt = 1.8892, phi = 5.2340, rap = 1.6192, |mass| = 0.0000
-pt = 0.5032, phi = 4.9272, rap = 1.6319, |mass| = 0.0000
-pt = 0.1631, phi = 3.1377, rap = 1.6448, |mass| = 0.0000
-pt = 0.3457, phi = 3.9398, rap = 1.6559, |mass| = 0.0000
-pt = 1.4796, phi = 5.7800, rap = 1.6586, |mass| = 0.0000
-pt = 0.1031, phi = 2.9014, rap = 1.6933, |mass| = 0.0000
-pt = 0.3576, phi = 2.6133, rap = 1.7359, |mass| = 0.0000
-pt = 0.8553, phi = 5.0883, rap = 1.7570, |mass| = 0.0000
-pt = 1.3839, phi = 2.4587, rap = 1.7761, |mass| = 0.0000
-pt = 5.8307, phi = 2.7960, rap = 1.8007, |mass| = 0.0000
-pt = 3.1964, phi = 0.9863, rap = 1.8074, |mass| = 0.0000
-pt = 0.5766, phi = 3.2073, rap = 1.8452, |mass| = 0.0000
-pt = 2.2420, phi = 2.7201, rap = 1.8501, |mass| = 0.0000
-pt = 5.8164, phi = 2.7803, rap = 1.8695, |mass| = 0.0000
-pt = 2.3878, phi = 2.6740, rap = 1.8754, |mass| = 0.0000
-pt = 1.3726, phi = 1.8314, rap = 1.9129, |mass| = 0.0000
-pt = 0.3131, phi = 2.7014, rap = 1.9332, |mass| = 0.0000
-pt = 0.5653, phi = 1.0835, rap = 1.9387, |mass| = 0.0000
-pt = 0.2684, phi = 2.0417, rap = 2.0056, |mass| = 0.0000
-pt = 0.0512, phi = 3.0374, rap = 2.0565, |mass| = 0.0000
-pt = 2.3945, phi = 3.0890, rap = 2.0668, |mass| = 0.0000
-pt = 0.1842, phi = 2.8988, rap = 2.0761, |mass| = 0.0000
-pt = 1.4991, phi = 2.8476, rap = 2.0854, |mass| = 0.0000
-pt = 0.7112, phi = 2.8541, rap = 2.0864, |mass| = 0.0000
-pt = 3.2188, phi = 2.9009, rap = 2.0914, |mass| = 0.0000
-pt = 8.1435, phi = 2.9158, rap = 2.1664, |mass| = 0.0000
-pt = 1.8454, phi = 2.6275, rap = 2.1787, |mass| = 0.0000
-pt = 14.2462, phi = 2.9648, rap = 2.1793, |mass| = 0.0000
-pt = 7.7212, phi = 2.8981, rap = 2.1901, |mass| = 0.0000
-pt = 8.3485, phi = 2.9077, rap = 2.1933, |mass| = 0.0000
-pt = 2.7375, phi = 2.8874, rap = 2.1989, |mass| = 0.0000
-pt = 5.9625, phi = 2.9472, rap = 2.2211, |mass| = 0.0000
-pt = 0.8056, phi = 2.3070, rap = 2.2353, |mass| = 0.0000
-pt = 2.0749, phi = 2.9942, rap = 2.2515, |mass| = 0.0000
-pt = 0.6627, phi = 2.6010, rap = 2.2538, |mass| = 0.0000
-pt = 7.0869, phi = 2.8741, rap = 2.2979, |mass| = 0.0000
-pt = 0.1211, phi = 4.2786, rap = 2.3323, |mass| = 0.0000
-pt = 1.4674, phi = 2.4512, rap = 2.3893, |mass| = 0.0000
-pt = 0.8952, phi = 4.9860, rap = 2.4013, |mass| = 0.0000
-pt = 0.6324, phi = 3.2281, rap = 2.4204, |mass| = 0.0000
-pt = 0.1931, phi = 2.3345, rap = 2.4268, |mass| = 0.0000
-pt = 0.5642, phi = 5.0863, rap = 2.4376, |mass| = 0.0000
-pt = 0.0375, phi = 4.6126, rap = 2.4678, |mass| = 0.0000
-pt = 0.7009, phi = 5.6329, rap = 2.4704, |mass| = 0.0000
-pt = 0.7562, phi = 5.3862, rap = 2.4942, |mass| = 0.0000
-pt = 0.3816, phi = 0.2073, rap = 2.5013, |mass| = 0.0000
-pt = 0.1342, phi = 5.5875, rap = 2.5304, |mass| = 0.0000
-pt = 0.2490, phi = 5.9251, rap = 2.5418, |mass| = 0.0000
-pt = 0.2596, phi = 1.3474, rap = 2.5588, |mass| = 0.0000
-pt = 0.5431, phi = 5.3477, rap = 2.6706, |mass| = 0.0000
-pt = 0.3763, phi = 5.2658, rap = 2.6781, |mass| = 0.0000
-pt = 0.0923, phi = 0.7093, rap = 2.7173, |mass| = 0.0000
-pt = 0.5464, phi = 5.2734, rap = 2.7487, |mass| = 0.0000
-pt = 0.4436, phi = 5.2945, rap = 2.7975, |mass| = 0.0000
-pt = 0.4076, phi = 4.8803, rap = 2.8057, |mass| = 0.0000
-pt = 1.2479, phi = 5.3496, rap = 2.9503, |mass| = 0.0000
-pt = 0.4169, phi = 4.5024, rap = 2.9632, |mass| = 0.0000
-pt = 0.1712, phi = 5.4471, rap = 2.9999, |mass| = 0.0000
-pt = 0.2471, phi = 2.5911, rap = 3.1385, |mass| = 0.0000
-pt = 0.0762, phi = 2.1888, rap = 3.1689, |mass| = 0.0000
-pt = 0.1484, phi = 0.8686, rap = 3.2117, |mass| = 0.0000
-pt = 0.2602, phi = 4.3488, rap = 3.2525, |mass| = 0.0000
-pt = 0.1832, phi = 3.8185, rap = 3.2569, |mass| = 0.0000
-pt = 0.5142, phi = 3.8074, rap = 3.3040, |mass| = 0.0000
-pt = 0.2522, phi = 4.8649, rap = 3.3606, |mass| = 0.0000
-pt = 0.0264, phi = 3.7510, rap = 3.5445, |mass| = 0.0000
-pt = 0.2938, phi = 0.6887, rap = 3.5684, |mass| = 0.0000
-pt = 0.2841, phi = 5.1156, rap = 3.6265, |mass| = 0.0000
-pt = 0.1426, phi = 0.3094, rap = 3.6296, |mass| = 0.0000
-pt = 0.2145, phi = 5.1615, rap = 3.6548, |mass| = 0.0000
-pt = 0.0539, phi = 3.4492, rap = 3.6786, |mass| = 0.0000
-pt = 0.3386, phi = 3.5986, rap = 3.6873, |mass| = 0.0000
-pt = 0.0885, phi = 5.6649, rap = 3.7001, |mass| = 0.0000
-pt = 0.1187, phi = 5.7991, rap = 3.7410, |mass| = 0.0000
-pt = 3.3433, phi = 5.8328, rap = 3.8261, |mass| = 0.0000
-pt = 1.1742, phi = 5.9007, rap = 3.8434, |mass| = 0.0000
-pt = 0.9101, phi = 3.7877, rap = 3.8660, |mass| = 0.0000
-pt = 0.4457, phi = 3.9673, rap = 3.9384, |mass| = 0.0000
-pt = 0.3687, phi = 3.4089, rap = 3.9813, |mass| = 0.0000
-pt = 0.0535, phi = 0.2902, rap = 3.9936, |mass| = 0.0000
-
-# original hard jets
-pt = 89.57, rap = 2.088, mass = 21.55, width = 0.1866
-pt = 144.5, rap = 0.6689, mass = 19.34, width = 0.06909
-pt = 220.1, rap = -1.068, mass = 18.09, width = 0.0337
-
-# unsubtracted full jets
-pt = 115.3, rap = 2.08, mass = 45.33, width = 0.2569
-pt = 167.3, rap = 0.6891, mass = 58.31, width = 0.1268
-pt = 230.5, rap = -1.065, mass = 46.2, width = 0.05647
-
-# subtracted full jets
-pt = 92.94, rap = 2.088, mass = 23.19, width = 0.2037
-pt = 146.8, rap = 0.6785, mass = 18.79, width = 0.07731
-pt = 216.3, rap = -1.071, mass = 12.47, width = 0.02855
-
+obtained rho with SignalFreeBackgroundEstimator: 13.2552
Index: contrib/contribs/SignalFreeBackgroundEstimator/trunk/SignalFreeBackgroundEstimator.cc
===================================================================
--- contrib/contribs/SignalFreeBackgroundEstimator/trunk/SignalFreeBackgroundEstimator.cc (revision 1377)
+++ contrib/contribs/SignalFreeBackgroundEstimator/trunk/SignalFreeBackgroundEstimator.cc (revision 1378)
@@ -1,280 +1,286 @@
//
// SignalFreeBackgroundEstimator package
// Questions/comments: peter.berta@cern.ch
//
// Copyright (c) 2024-, Peter Berta
//
//----------------------------------------------------------------------
// 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 "SignalFreeBackgroundEstimator.hh"
FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
namespace contrib{
void SignalFreeBackgroundEstimator::set_particles(const std::vector& particles, const std::vector& signalParticles, const double& measureOfPileup, const std::vector& chargedSignalParticles){
assert(all_tiles_equal_area());
_cache_available = false;
_cached_estimate.reset();
_cached_estimate.set_has_sigma(false);
_cached_estimate.set_mean_area(mean_tile_area());
// Cluster signal particles
fastjet::AreaDefinition areaDef(fastjet::active_area, fastjet::GhostedAreaSpec(_rapidity_max));
fastjet::JetDefinition jetDef(fastjet::antikt_algorithm, _signal_jets_R);
fastjet::ClusterSequenceArea clusSeq(signalParticles, jetDef, areaDef);
std::vector allSignalJets = sorted_by_pt(clusSeq.inclusive_jets());
- std::vector signalPoints;
double jet_rho_min=_jet_rho_min_intercept;
if (measureOfPileup>0) jet_rho_min=_jet_rho_min_intercept+sqrt(measureOfPileup)*_jet_rho_min_slope;
for (size_t i = 0; i < allSignalJets.size(); ++i) {
- if (allSignalJets.at(i).pt()/allSignalJets.at(i).area() > jet_rho_min) signalPoints.emplace_back(allSignalJets.at(i));
+ if (allSignalJets.at(i).pt()/allSignalJets.at(i).area() > jet_rho_min) _signal_seeds.emplace_back(allSignalJets.at(i));
}
if (chargedSignalParticles.size()>0){
fastjet::ClusterSequenceArea clusSeq_charged(chargedSignalParticles, jetDef, areaDef);
std::vector chargedJets = sorted_by_pt(clusSeq_charged.inclusive_jets());
for (size_t i = 0; i < chargedJets.size(); ++i) {
- if (chargedJets.at(i).pt()/chargedJets.at(i).area() > _jet_rho_min_charged) signalPoints.emplace_back(chargedJets.at(i));
+ if (chargedJets.at(i).pt()/chargedJets.at(i).area() > _jet_rho_min_charged) _signal_seeds.emplace_back(chargedJets.at(i));
}
}
double scalarPtSumSignal=0;
for (size_t i = 0; i < signalParticles.size(); ++i) scalarPtSumSignal+=signalParticles[i].pt();
_signal_fraction=scalarPtSumSignal/14000.;
// Determine tile area outside the exclusion area
std::vector tileStates(n_tiles(),false);
std::vector tileAreas(n_tiles());
int nGhosts_phi=dphi()/_ghost_size;
int nGhosts_rap=drap()/_ghost_size;
fastjet::PseudoJet ghost(0,0,0,1);
int nPhi = std::round(fastjet::twopi / dphi());
for (int i = 0; i < n_tiles(); ++i) {
int tileIndexPhi = i % nPhi;
int tileIndexRap = i / nPhi;
int nGhostsCloseToHardScatter=0;
for (int iphi = 0; iphi < nGhosts_phi; ++iphi) {
double phi=(tileIndexPhi + (iphi+0.5)/(double)nGhosts_phi)*dphi();
for (int irap = 0; irap < nGhosts_rap; ++irap) {
double rap=(tileIndexRap + (irap+0.5)/(double)nGhosts_rap)*drap() + rapmin();
ghost.reset_momentum_PtYPhiM(1,rap,phi,1e-200);
- for (size_t isignal = 0; isignal < signalPoints.size(); ++isignal) {
- if (signalPoints.at(isignal).delta_R(ghost) < _exclusion_deltaR) {
+ for (size_t isignal = 0; isignal < _signal_seeds.size(); ++isignal) {
+ if (_signal_seeds.at(isignal).delta_R(ghost) < _exclusion_deltaR) {
nGhostsCloseToHardScatter++;
break;
}
}
}
}
double areaFraction=1-nGhostsCloseToHardScatter/(double)(nGhosts_phi*nGhosts_rap);
tileAreas[i]=mean_tile_area()*areaFraction;
if (tileAreas[i]>_tile_area_min) tileStates[i]=true;
}
std::vector tilePt(n_tiles(), 0.);
std::vector tileMtMinusPt(n_tiles(), 0.);
for (size_t i = 0; i < particles.size(); ++i) {
// Find index of the tile the particle belongs to
int tileIndex = tile_index(particles[i]);
if (tileIndex < 0) {
std::cout << "WARNING SignalFreeBackgroundEstimator::set_particles: Tile index is negative. Particle (pt,rap,phi,m): " << particles[i].pt() << " " << particles[i].rap() << " " << particles[i].phi() << " " << particles[i].m() << std::endl;
continue;
}
bool closeToHardScatter = false;
- for (size_t j = 0; j < signalPoints.size(); ++j) {
- if (signalPoints.at(j).delta_R(particles[i]) < _exclusion_deltaR) {
+ for (size_t j = 0; j < _signal_seeds.size(); ++j) {
+ if (_signal_seeds.at(j).delta_R(particles[i]) < _exclusion_deltaR) {
closeToHardScatter = true;
break;
}
}
if (closeToHardScatter) continue;
if (_rescaling_class) {
double r = (*_rescaling_class)(particles[i]);
tilePt[tileIndex] += particles[i].pt() / r;
if (_enable_rho_m) tileMtMinusPt[tileIndex] += (particles[i].mt()-particles[i].pt()) / r;
}
else{
tilePt[tileIndex] += particles[i].pt();
if (_enable_rho_m) tileMtMinusPt[tileIndex] += particles[i].mt()-particles[i].pt();
}
}
/// Compute rho and rho_m for each tile, and assign weights based on tile area
std::vector > tileRho, tileRhom;
_weights_sum=0;
for (int i = 0; i < n_tiles(); ++i) {
if (!tileStates[i]) continue;
tileRho.emplace_back(std::make_pair(tilePt[i]/tileAreas[i],tileAreas[i]));
_weights_sum+=tileAreas[i];
if (_enable_rho_m) tileRhom.emplace_back(std::make_pair(tileMtMinusPt[i]/tileAreas[i],tileAreas[i]));
}
if (tileRho.size()==0){
std::cout << "WARNING SignalFreeBackgroundEstimator::set_particles: All tiles are excluded in this event, and it is not possible to estimate rho with SignalFreeBackgroundEstimator! Estimating it with GridMedianBackgroundEstimator instead. Crosscheck whether it is expected that all tiles are excluded in this event and try to optimize the parameters of SignalFreeBackgroundEstimator, if needed." << std::endl;
std::unique_ptr gridMedian=std::make_unique(_rapidity_max, 0.55);
if (_rescaling_class) gridMedian->set_rescaling_class(_rescaling_class);
gridMedian->set_particles(particles);
_cached_estimate=gridMedian->estimate();
_cache_available = true;
return;
}
/// Compute and set the final rho and rho_m
double rho=_compute_weighted_median(tileRho);
_cached_estimate.set_rho(rho);
if (_enable_rho_m){
double rhom=_compute_weighted_median(tileRhom);;
_cached_estimate.set_has_rho_m(true);
_cached_estimate.set_rho_m(rhom);
}
_cache_available = true;
}
/// Compute weighted median based on the weights corresponding to tile areas.
double SignalFreeBackgroundEstimator::_compute_weighted_median(std::vector > &tileRho) const{
sort(tileRho.begin(),tileRho.end());
// Get the actual center in case the user requested to use floating center:
double shift=0;
if (_regulator_for_floating_center>=0) shift=(_signal_fraction-_signal_fraction_min)*_regulator_for_floating_center;
if (shift<0) shift=0;
if (shift>_shift_max) shift=_shift_max;
if (shift>_center-_half_window) shift=_center-_half_window;
double center_shifted=_center-shift;
// Compute rho from the window specified by the user:
double rho_nominator=0;
double rho_denominator=0;
bool insideWindow=false;
double fractions_partialSum=0;
for (size_t i = 0; i < tileRho.size(); ++i) {
double fraction_i=tileRho[i].second/_weights_sum;
double remaining_to_start=center_shifted-_half_window-fractions_partialSum;
double remaining_to_stop=center_shifted+_half_window-fractions_partialSum;
if (remaining_to_start* rescaling_class_in) {
if (!_cache_available) {
_warning_rescaling.warn("SignalFreeBackgroundEstimator::set_rescaling_class: Found cached result. Set particles again to obtain correct calculation!");
}
BackgroundEstimatorBase::set_rescaling_class(rescaling_class_in);
}
+ /// User can add his/her own signal seeds (either additional seeds to those obtained from jet clustering or only these user-defined seeds can be used)
+ void SignalFreeBackgroundEstimator::add_additional_seeds(const std::vector &additional_seeds){
+ for (size_t iseed = 0; iseed < additional_seeds.size(); ++iseed) {
+ _signal_seeds.emplace_back(additional_seeds[iseed]);
+ }
+ }
+
} // contrib namespace
FASTJET_END_NAMESPACE
Index: contrib/contribs/SignalFreeBackgroundEstimator/trunk/example.cc
===================================================================
--- contrib/contribs/SignalFreeBackgroundEstimator/trunk/example.cc (revision 1377)
+++ contrib/contribs/SignalFreeBackgroundEstimator/trunk/example.cc (revision 1378)
@@ -1,204 +1,212 @@
//
//----------------------------------------------------------------------
// 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
+// ./example < ../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 "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
-
+ /* Commenting event-wide pileup corrector (ICS) in this example, because currently no other contribs can be used.
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;
+ /* 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);
+ 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);
+ // 3. If tracking information is available, then the charged signal particles can be used to find the signal seeds, without the need to do any event-wide pileup correction in advance.
+ vector estimated_signal; // empty vector should be provided in this case
+ bge_rho.set_particles(full_event,estimated_signal,-1,*hard_event_charged);
+ cout << "obtained rho with SignalFreeBackgroundEstimator: " << bge_rho.rho() << endl;
+ // vector corrected_event=subtractor.subtract_event(full_event); // use this to get the corrected event using ICS
+
+ // 4. The startegy to find signal seeds described in the above steps 1, 2, and 3 can be arbitrarily combined. For example, the combination of steps 1 and 3 can be done as follows:
+ // bge_rho_grid.set_particles(full_event); // Alternatively, one can use only the neutral component, i.e. perform Charged Hadron Subtraction in adavnace, and then 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. The user can provide his/her own set of signal seeds:
+ // vector additional_seeds=....; // the user obtaines the signal seeds based on his/her preferred method
+ // bge_rho.add_additional_seeds(additional_seeds);
+ // bge_rho.set_particles(full_event); // no pileup-corrected particles or charged signal particles are provided here. Alternatively, the user can also provide them.
- // 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:
+ // 6. 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 and jet variable printing (in case pileup correction is applied)
// 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