Index: contrib/contribs/SignalFreeBackgroundEstimator/trunk/SignalFreeBackgroundEstimator.cc
===================================================================
--- contrib/contribs/SignalFreeBackgroundEstimator/trunk/SignalFreeBackgroundEstimator.cc (revision 0)
+++ contrib/contribs/SignalFreeBackgroundEstimator/trunk/SignalFreeBackgroundEstimator.cc (revision 1377)
@@ -0,0 +1,280 @@
+//
+// 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 (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));
+ }
+ }
+ 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) {
+ 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) {
+ 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);
+ }
+
+} // contrib namespace
+
+
+FASTJET_END_NAMESPACE
Index: contrib/contribs/SignalFreeBackgroundEstimator/trunk/SignalFreeBackgroundEstimator.hh
===================================================================
--- contrib/contribs/SignalFreeBackgroundEstimator/trunk/SignalFreeBackgroundEstimator.hh (revision 0)
+++ contrib/contribs/SignalFreeBackgroundEstimator/trunk/SignalFreeBackgroundEstimator.hh (revision 1377)
@@ -0,0 +1,167 @@
+//
+// 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.
+ fastjet::BackgroundEstimate estimate() const override;
+
+ /// 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.
+ double rho() const override;
+
+ /// 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
+ 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:
+ /// 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;
+
+ };
+
+}
+
+FASTJET_END_NAMESPACE
+
+#endif // __SIGNAL_FREE_BACKGROUND_ESTIMATOR_H__
+