Index: contrib/contribs/ConstituentSubtractor/tags/1.4.7/RescalingClasses.cc
===================================================================
--- contrib/contribs/ConstituentSubtractor/tags/1.4.7/RescalingClasses.cc (revision 0)
+++ contrib/contribs/ConstituentSubtractor/tags/1.4.7/RescalingClasses.cc (revision 1373)
@@ -0,0 +1,215 @@
+// ConstituentSubtractor package
+// Questions/comments: peter.berta@cern.ch, Martin.Spousta@cern.ch, David.W.Miller@uchicago.edu, Rupert.Leitner@mff.cuni.cz
+
+// Copyright (c) 2014-, Peter Berta, Martin Spousta, David W. Miller, Rupert Leitner
+//
+//----------------------------------------------------------------------
+// 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 "RescalingClasses.hh"
+
+FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
+
+
+namespace contrib{
+
+ ///-----------------------------------------------------
+ /// BackgroundRescalingYPhi
+ BackgroundRescalingYPhi::BackgroundRescalingYPhi(double v2, double v3, double v4, double psi, double a1, double sigma1, double a2, double sigma2){
+ _v2=v2;
+ _v3=v3;
+ _v4=v4;
+ _psi=psi;
+ _a1=a1;
+ _sigma1=sigma1;
+ _a2=a2;
+ _sigma2=sigma2;
+ _use_rap=true;
+ _use_phi=true;
+ }
+
+ void BackgroundRescalingYPhi::use_rap_term(bool use_rap){
+ _use_rap=use_rap;
+ }
+
+ void BackgroundRescalingYPhi::use_phi_term(bool use_phi){
+ _use_phi=use_phi;
+ }
+
+ double BackgroundRescalingYPhi::result(const PseudoJet & particle) const {
+ double phi_term=1;
+ if (_use_phi){
+ double phi=particle.phi();
+ phi_term=1 + 2*_v2*_v2*cos(2*(phi-_psi)) + 2*_v3*_v3*cos(3*(phi-_psi)) + 2*_v4*_v4*cos(4*(phi-_psi));
+ }
+ double rap_term=1;
+ if (_use_rap){
+ double y=particle.rap();
+ rap_term=_a1*exp(-y*y/(2*_sigma1*_sigma1)) + _a2*exp(-y*y/(2*_sigma2*_sigma2));
+ }
+
+ return phi_term*rap_term;
+ }
+
+
+
+
+
+
+ ///----------------------------------------------------
+ /// BackgroundRescalingYPhiUsingVectorForY
+ BackgroundRescalingYPhiUsingVectorForY::BackgroundRescalingYPhiUsingVectorForY(double v2, double v3, double v4, double psi, std::vector values, std::vector rap_binning, bool interpolate){
+ _v2=v2;
+ _v3=v3;
+ _v4=v4;
+ _psi=psi;
+ _values=values;
+ _rap_binning=rap_binning;
+ _use_phi=true;
+ _interpolate=interpolate;
+ if (_rap_binning.size()>=2){
+ _use_rap=true;
+ if (_values.size()!=_rap_binning.size()-1) throw Error("BackgroundRescalingYPhiUsingVectorForY (from ConstituentSubtractor) The input vectors have wrong dimension. The vector with binning shuld have the size by one higher than the vector with values.");
+ }
+ else _use_rap=false;
+ }
+
+ void BackgroundRescalingYPhiUsingVectorForY::use_rap_term(bool use_rap){
+ _use_rap=use_rap;
+ if (use_rap && _rap_binning.size()<2) throw Error("BackgroundRescalingYPhiUsingVectorForY (from ConstituentSubtractor) Requested rapidity rescaling, but the vector with binning has less than two elements!");
+ }
+
+ void BackgroundRescalingYPhiUsingVectorForY::use_phi_term(bool use_phi){
+ _use_phi=use_phi;
+ }
+
+ double BackgroundRescalingYPhiUsingVectorForY::result(const PseudoJet & particle) const {
+ double phi_term=1;
+ if (_use_phi){
+ double phi=particle.phi();
+ phi_term=1 + 2*_v2*_v2*cos(2*(phi-_psi)) + 2*_v3*_v3*cos(3*(phi-_psi)) + 2*_v4*_v4*cos(4*(phi-_psi));
+ }
+ unsigned int nBins=_rap_binning.size()-1;
+ double rap_term=1;
+ if (_use_rap){
+ int rap_index=0;
+ double rap=particle.rap();
+ if (rap<_rap_binning[0]) rap_index=0; // take the lowest rapidity bin in this case
+ else if (rap>=_rap_binning[nBins]) rap_index=nBins-1; // take the highest rapidity bin in this case
+ else{
+ for (unsigned int i=1;i=(_rap_binning[nBins-1]+_rap_binning[nBins])/2.) rap_term=_values[nBins-1];
+ else{
+ double value_low,value_high,rap_bin_low,rap_bin_high;
+ double rap_bin_center=(_rap_binning[rap_index]+_rap_binning[rap_index+1])/2.;
+ if (rap > values, std::vector rap_binning, std::vector phi_binning){
+ _values=values;
+ _rap_binning=rap_binning;
+ _phi_binning=phi_binning;
+ if (_rap_binning.size()>=2) _use_rap=true;
+ else _use_rap=false;
+ if (_phi_binning.size()>=2) _use_phi=true;
+ else _use_phi=false;
+ }
+
+ void BackgroundRescalingYPhiUsingVectors::use_rap_term(bool use_rap){
+ _use_rap=use_rap;
+ if (use_rap && _rap_binning.size()<2) throw Error("BackgroundRescalingYPhiUsingVectors (from ConstituentSubtractor) Requested rapidity rescaling, but the vector with binning has less than two elements!");
+ }
+
+ void BackgroundRescalingYPhiUsingVectors::use_phi_term(bool use_phi){
+ _use_phi=use_phi;
+ if (use_phi && _phi_binning.size()<2) throw Error("BackgroundRescalingYPhiUsingVectors (from ConstituentSubtractor) Requested azimuth rescaling, but the vector with binning has less than two elements!");
+ }
+
+ double BackgroundRescalingYPhiUsingVectors::result(const PseudoJet & particle) const {
+ unsigned int phi_index=0;
+ if (_use_phi){
+ double phi=particle.phi();
+ if (phi<_phi_binning[0] || phi>=_phi_binning[_phi_binning.size()-1]) throw Error("BackgroundRescalingYPhiUsingVectors (from ConstituentSubtractor) The phi binning does not correspond to the phi binning of the particles.");
+ for (unsigned int i=1;i<_phi_binning.size();++i){
+ if (phi<_phi_binning[i]){
+ phi_index=i-1;
+ break;
+ }
+ }
+ }
+ unsigned int rap_index=0;
+ if (_use_rap){
+ double rap=particle.rap();
+ if (rap<_rap_binning[0]) rap_index=0; // take the lowest rapidity bin in this case
+ else if (rap>=_rap_binning[_rap_binning.size()-1]) rap_index=_rap_binning.size()-2; // take the highest rapidity bin in this case
+ else{
+ for (unsigned int i=1;i<_rap_binning.size();++i){
+ if (rap<_rap_binning[i]){
+ rap_index=i-1;
+ break;
+ }
+ }
+ }
+ }
+ if (_values.size()<=rap_index) throw Error("BackgroundRescalingYPhiUsingVectors (from ConstituentSubtractor) The input vector > with values has wrong size.");
+ else if (_values[rap_index].size()<=phi_index) throw Error("BackgroundRescalingYPhiUsingVectors (from ConstituentSubtractor) The input vector > with values has wrong size.");
+
+ return _values[rap_index][phi_index];
+ }
+
+
+} // namespace contrib
+
+
+FASTJET_END_NAMESPACE
+
+
+
Index: contrib/contribs/ConstituentSubtractor/tags/1.4.7/example_whole_event_using_charged_info.cc
===================================================================
--- contrib/contribs/ConstituentSubtractor/tags/1.4.7/example_whole_event_using_charged_info.cc (revision 0)
+++ contrib/contribs/ConstituentSubtractor/tags/1.4.7/example_whole_event_using_charged_info.cc (revision 1373)
@@ -0,0 +1,146 @@
+// $Id$
+//
+//----------------------------------------------------------------------
+// Example on how to do pileup correction on the whole event
+//
+// run it with
+// ./example_whole_event_using_charged_info < ../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/Selector.hh"
+#include "fastjet/tools/GridMedianBackgroundEstimator.hh"
+#include "fastjet/ClusterSequence.hh"
+#include "ConstituentSubtractor.hh" // In external code, this should be fastjet/contrib/ConstituentSubtractor.hh
+
+#include "functions.hh"
+
+
+using namespace std;
+using namespace fastjet;
+
+//----------------------------------------------------------------------
+int main(){
+ // set up before event loop:
+ contrib::ConstituentSubtractor subtractor; // no need to provide background estimator in this case
+ 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
+ subtractor.set_max_distance(0.3); // free parameter for the maximal allowed distance between particle i and ghost k
+ subtractor.set_alpha(1); // free parameter for the distance measure (the exponent of particle pt). The larger the parameter alpha, the more are favoured the lower pt particles in the subtraction process
+ subtractor.set_ghost_area(0.01); // free parameter for the density of ghosts. The smaller, the better - but also the computation is slower.
+ subtractor.set_do_mass_subtraction(); // use this line if also the mass term sqrt(pT^2+m^2)-pT should be corrected or not. It is necessary to specify it like this because the function set_common_bge_for_rho_and_rhom cannot be used in this case.
+ double CBS=1.0; // choose the scale for scaling the background charged particles
+ double CSS=1.0; // choose the scale for scaling the signal charged particles
+ subtractor.set_remove_particles_with_zero_pt_and_mass(false); // set to false if you want to have also the zero pt and mtMinuspt particles in the output. Set to true, if not. The choice has no effect on the performance. By deafult, these particles are removed - this is the recommended way since then the output contains much less particles, and therefore the next step (e.g. clustering) is faster. In this example, it is set to false to make sure that this test is successful on all systems (mac, linux).
+ subtractor.set_grid_size_background_estimator(0.6); // set the grid size (not area) for the background estimation with GridMedianBackgroundEstimation which is used within CS correction using charged info
+ cout << subtractor.description() << endl;
+
+
+
+ // 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);
+
+ double max_eta=4; // specify the maximal pseudorapidity for the particles used in the subtraction
+ double max_eta_jet=3; // the maximal pseudorapidity for selected jets
+
+ // keep the particles up to 4 units in rapidity
+ 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, " << full_event.size() - hard_event.size() << " background particles, " << hard_event_charged->size() << " signal charged particles, and " << background_event_charged->size() << " background charged particles" << " with pseudo-rapidity |eta|<4" << endl;
+
+
+ // do the clustering with ghosts and get the jets
+ //----------------------------------------------------------
+ JetDefinition jet_def(antikt_algorithm, 0.7);
+ Selector sel_jets = SelectorNHardest(3) * SelectorAbsEtaMax(max_eta_jet);
+
+ 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());
+
+ // correction of the whole event using charged info
+ vector corrected_event=subtractor.subtract_event_using_charged_info(full_event,CBS,*background_event_charged,CSS,*hard_event_charged,max_eta);
+
+ ios::fmtflags f( cout.flags() );
+ cout << setprecision(5) << fixed;
+ cout << endl << "Corrected particles in the whole event:" << endl;
+ for (unsigned int i=0; i corrected_jets = sel_jets(clust_seq_corr.inclusive_jets());
+
+ // shape variables:
+ //----------------------------------------------------------
+ JetWidth width;
+
+ // subtract and print the result
+ //----------------------------------------------------------
+ cout.flags( f );
+ cout << setprecision(6);
+ cout << "Original hard jets" << endl;
+ for (unsigned int i=0; i0) _use_max_distance=true;
+ else _use_max_distance=false;
+ }
+
+
+ ///
+ /// Constructor that takes an externally supplied value for rho and, optionally, for rho_m.
+ ConstituentSubtractor::ConstituentSubtractor(double rho, double rhom, double alpha, double max_distance, Distance distance):
+ _bge_rho(0), _bge_rhom(0), _common_bge(false), _rhom_from_bge_rhom(false), _rho(rho), _rhom(rhom), _externally_supplied_rho_rhom(true), _alpha(alpha), _distance(distance), _max_distance(max_distance) {
+ if (_max_distance>0) _use_max_distance=true;
+ else _use_max_distance=false;
+ assert(_rho >= 0);
+ assert(_rhom >= 0);
+ }
+
+
+ void ConstituentSubtractor::_initialize_common(){
+ if (_max_eta<=0) throw Error("ConstituentSubtractor::initialize_common: The value for eta cut was not set or it is negative. It needs to be set before using the function initialize");
+ if (_masses_to_zero && _do_mass_subtraction) throw Error("ConstituentSubtractor::initialize_common: It is specified to do mass subtraction and also to keep the masses at zero. Something is wrong.");
+ if (_masses_to_zero && _scale_fourmomentum) throw Error("ConstituentSubtractor::initialize_common: It is specified to do scaling of fourmomenta and also to keep the masses at zero. Something is wrong.");
+ if (_do_mass_subtraction && _scale_fourmomentum) throw Error("ConstituentSubtractor::initialize_common: It is specified to do mass subtraction and also to do scaling of fourmomenta. Something is wrong.");
+ this->construct_ghosts_uniformly(_max_eta);
+ }
+
+
+ void ConstituentSubtractor::initialize(){
+ this->_initialize_common();
+ }
+
+
+
+ void ConstituentSubtractor::set_background_estimator(fastjet::BackgroundEstimatorBase *bge_rho, fastjet::BackgroundEstimatorBase *bge_rhom){
+ _bge_rho=bge_rho;
+ _bge_rhom=bge_rhom;
+ }
+
+
+ void ConstituentSubtractor::set_scalar_background_density(double rho, double rhom){
+ _rho=rho;
+ _rhom=rhom;
+ assert(_rho >= 0);
+ assert(_rhom >= 0);
+ _externally_supplied_rho_rhom=true;
+ _common_bge=false;
+ }
+
+ ///----------------------------------------------------------------------
+ /// the action on a given jet
+ fastjet::PseudoJet ConstituentSubtractor::result(const PseudoJet &jet) const{
+ // make sure we have a BGE or a rho value
+ if (!_bge_rho && !_externally_supplied_rho_rhom){
+ throw Error("ConstituentSubtractor::result() constituent subtraction needs a BackgroundEstimator or a value for rho.");
+ }
+ if (_ghosts_constructed) throw Error("ConstituentSubtractor::result() The ghosts are constructed, but they are not needed when using this function. When you want to perform jet-by-jet correction, initialize a new ConstituentSubtractor without construction of ghosts.");
+
+ ///----------------------------------------------------------------------
+ /// sift ghosts and particles in the input jet
+ std::vector particles, ghosts;
+ SelectorIsPureGhost().sift(jet.constituents(), ghosts, particles);
+
+ std::vector selected_particles,unselected_particles;
+ if (_particle_selector) _particle_selector->sift(particles, selected_particles, unselected_particles);
+ else selected_particles=particles;
+
+
+ std::vector ghosts_area;
+ unsigned long nGhosts=ghosts.size();
+ for (unsigned int j=0;j backgroundProxies=this->get_background_proxies_from_ghosts(ghosts,ghosts_area);
+ std::vector subtracted_particles=this->do_subtraction(selected_particles,backgroundProxies);
+ if (_particle_selector) subtracted_particles.insert(subtracted_particles.end(), unselected_particles.begin(), unselected_particles.end());
+ fastjet::PseudoJet subtracted_jet=join(subtracted_particles);
+ subtracted_jet.set_user_index(jet.user_index());
+ subtracted_jet.user_info_shared_ptr() = jet.user_info_shared_ptr();
+
+ return subtracted_jet;
+ }
+
+
+
+
+ std::vector ConstituentSubtractor::subtract_event(std::vector const &particles, double max_eta){
+ if (fabs(_max_eta/max_eta-1)>1e-5 && max_eta>0){
+ _ghosts_constructed=false;
+ _max_eta=max_eta;
+ }
+ if (!_ghosts_constructed) this->construct_ghosts_uniformly(_max_eta);
+ return subtract_event(particles);
+ }
+
+
+ std::vector ConstituentSubtractor::subtract_event(std::vector const &particles, std::vector const *hard_proxies){
+ std::vector backgroundProxies=this->get_background_proxies_from_ghosts(_ghosts,_ghosts_area);
+ std::vector selected_particles,unselected_particles;
+ for (unsigned int iparticle=0;iparticle_max_eta) continue;
+ if (particles[iparticle].pt()pass(particles[iparticle])) selected_particles.push_back(particles[iparticle]);
+ else unselected_particles.push_back(particles[iparticle]);
+ }
+ else selected_particles.push_back(particles[iparticle]);
+ }
+ if (_use_nearby_hard){
+ if (hard_proxies) _hard_proxies=hard_proxies;
+ else throw Error("ConstituentSubtractor::subtract_event: It was requested to use closeby hard proxies but they were not provided in this function!");
+ }
+ else if (hard_proxies) throw Error("ConstituentSubtractor::subtract_event: Hard proxies were provided but the set_use_hard_proxies function was not used before initialization. It needs to be called before initialization!");
+ std::vector subtracted_particles=this->do_subtraction(selected_particles,backgroundProxies);
+ if (_particle_selector) subtracted_particles.insert(subtracted_particles.end(), unselected_particles.begin(), unselected_particles.end());
+ return subtracted_particles;
+ }
+
+
+
+
+ std::vector ConstituentSubtractor::get_background_proxies_from_ghosts(std::vector const &ghosts,std::vector const &ghosts_area) const{
+ unsigned long nGhosts=ghosts.size();
+ std::vector proxies;
+ std::vector pt;
+ std::vector mtMinusPt;
+ proxies.reserve(nGhosts);
+ pt.reserve(nGhosts);
+ mtMinusPt.reserve(nGhosts);
+ if (_externally_supplied_rho_rhom){
+ for (unsigned int j=0;jrho(ghosts[j])*ghosts_area[j]);
+ if (_bge_rhom){
+ /* if (_rhom_from_bge_rhom){
+#if FASTJET_VERSION_NUMBER >= 30100
+ for (unsigned int j=0;jrho_m(ghosts[j])*ghosts_area[j]);
+#else
+ throw(Error("ConstituentSubtractor:: _rhom_from_bge_rhom not allowed for FJ<3.1"));
+#endif // end of code specific to FJ >= 3.1
+ } else {
+ for (unsigned int j=0;jrho(ghosts[j])*ghosts_area[j]);
+ }*/
+
+
+ // since FJ 3.1.0, some background estimators have an automatic internal calculation of rho_m
+#if FASTJET_VERSION_NUMBER >= 30100
+ // check if the BGE has internal support for rho_m
+ if (_bge_rhom->has_rho_m()){
+ for (unsigned int j=0;jrho_m(ghosts[j])*ghosts_area[j]);
+ } else {
+ throw Error("ConstituentSubtractor: The provided background estimator for rho_m has no support to compute rho_m, and other option to get it is not available in ConstituentSubtractor.");
+ }
+#endif // end of code specific to FJ >= 3.1
+#if FASTJET_VERSION_NUMBER < 30100
+ throw Error("ConstituentSubtractor: You provided a background estimator for rho_m, but this is not supported in ConstituentSubtractor when using fastjet version smaller than 30100.");
+#endif
+ }
+ else if (_common_bge){
+ // since FJ 3.1.0, some background estimators have an automatic internal calculation of rho_m
+#if FASTJET_VERSION_NUMBER >= 30100
+ // check if the BGE has internal support for rho_m
+ if (_bge_rho->has_rho_m()){
+ for (unsigned int j=0;jrho_m(ghosts[j])*ghosts_area[j]);
+ } else {
+#endif // end of code specific to FJ >= 3.1
+ BackgroundJetPtMDensity m_density;
+ JetMedianBackgroundEstimator *jmbge = dynamic_cast(_bge_rho);
+ const FunctionOfPseudoJet * orig_density = jmbge->jet_density_class();
+ jmbge->set_jet_density_class(&m_density);
+ for (unsigned int j=0;jrho(ghosts[j])*ghosts_area[j]);
+ jmbge->set_jet_density_class(orig_density);
+#if FASTJET_VERSION_NUMBER >= 30100
+ }
+#endif
+ }
+ else { // a single bge, only rho requested
+ for (unsigned int j=0;j= 30100
+ // In FJ3.1 and BGE with rho_m support, add a warning, similar to that in Subtractor
+ double const rho_m_warning_threshold = 1e-5;
+ if (_bge_rho->has_rho_m() && _bge_rho->rho_m()>rho_m_warning_threshold*_bge_rho->rho() && !_masses_to_zero && !_scale_fourmomentum){
+ _warning_unused_rhom.warn("ConstituentSubtractor:: Background estimator indicates non-zero rho_m, but the ConstituentSubtractor does not use rho_m information, nor the masses are set to zero, nor the 4-momentum is scaled. Consider calling set_common_bge_for_rho_and_rhom() to include the rho_m information; or call set_keep_original_masses(false) to set masses for all particles to zero; or call set_scale_fourmomentum to scale the fourmomentum.");
+ }
+#endif
+ }
+ }
+
+
+ fastjet::PseudoJet proxy(0,0,0,1);
+ for (unsigned int j=0;j0) mass=sqrt(mass_squared);
+ proxy.reset_momentum_PtYPhiM(pt[j],ghosts[j].rap(),ghosts[j].phi(),mass);
+ proxies.push_back(proxy);
+ }
+ return proxies;
+ }
+
+
+
+
+ double ConstituentSubtractor::_get_transformed_distance(double const &distance) const{
+ double max_distance_transformed=-1;
+ if (_distance==ConstituentSubtractor::deltaR) max_distance_transformed=pow(distance,2);
+ if (_distance==ConstituentSubtractor::angle) max_distance_transformed=-cos(distance);
+ return max_distance_transformed;
+ }
+
+
+
+ std::vector ConstituentSubtractor::do_subtraction(std::vector const &particles, std::vector const &backgroundProxies,std::vector *remaining_backgroundProxies) const{
+ unsigned int nBackgroundProxies=backgroundProxies.size();
+ unsigned int nParticles=particles.size();
+ double max_distance_transformed=this->_get_transformed_distance(_max_distance);
+
+ // sort particles according to rapidity to achieve faster performance for the whole event subtraction
+ std::vector particles_sorted=particles;
+ std::sort(particles_sorted.begin(),particles_sorted.end(),ConstituentSubtractor::_rapidity_sorting);
+
+ // get the kinematic variables for particles and background proxies in advance to achieve faster performance
+ std::vector particles_phi,particles_rap,particles_pt,particles_mt,particles_factor,particles_px_normalized,particles_py_normalized,particles_pz_normalized;
+ particles_phi.reserve(nParticles);
+ particles_rap.reserve(nParticles);
+ particles_pt.reserve(nParticles);
+ particles_mt.reserve(nParticles);
+ particles_factor.reserve(nParticles);
+ particles_px_normalized.reserve(nParticles);
+ particles_py_normalized.reserve(nParticles);
+ particles_pz_normalized.reserve(nParticles);
+
+ double pt_factor=1;
+ double polarAngle_factor=1;
+ double nearby_hard_factor=1;
+ double max_distance_from_hard_transformed=this->_get_transformed_distance(_nearby_hard_radius);
+ for (unsigned int i=0;i1e-5) pt_factor=pow(particles_pt[i],_alpha);
+ double momentum=sqrt(particles_sorted[i].pt2()+particles_sorted[i].pz()*particles_sorted[i].pz());
+ if (fabs(_polarAngleExp)>1e-5) polarAngle_factor=pow(particles_pt[i]/momentum,_polarAngleExp);
+ if (_distance==ConstituentSubtractor::angle){
+ particles_px_normalized.push_back(particles_sorted[i].px()/momentum);
+ particles_py_normalized.push_back(particles_sorted[i].py()/momentum);
+ particles_pz_normalized.push_back(particles_sorted[i].pz()/momentum);
+ }
+ if (_use_nearby_hard){
+ nearby_hard_factor=1;
+ double distance_from_hard_transformed=-1;
+ for (unsigned int ihard=0;ihard<_hard_proxies->size();++ihard){
+ if (_distance==ConstituentSubtractor::deltaR){
+ double deltaPhi=fabs(_hard_proxies->at(ihard).phi()-particles_phi[i]);
+ if (deltaPhi>pi) deltaPhi=twopi-deltaPhi;
+ double deltaRap=_hard_proxies->at(ihard).rap()-particles_rap[i];
+ distance_from_hard_transformed = deltaPhi*deltaPhi+deltaRap*deltaRap;
+ }
+ if (_distance==ConstituentSubtractor::angle){
+ distance_from_hard_transformed=-(particles_px_normalized[i]*_hard_proxies->at(ihard).px()+particles_py_normalized[i]*_hard_proxies->at(ihard).py()+particles_pz_normalized[i]*_hard_proxies->at(ihard).pz())/sqrt(_hard_proxies->at(ihard).pt2()*_hard_proxies->at(ihard).pt2()+_hard_proxies->at(ihard).pz()*_hard_proxies->at(ihard).pz());
+ }
+ if (distance_from_hard_transformed<=max_distance_from_hard_transformed){
+ nearby_hard_factor=_nearby_hard_factor;
+ break;
+ }
+ }
+ }
+ particles_factor.push_back(pt_factor*polarAngle_factor*nearby_hard_factor);
+ }
+
+ std::vector backgroundProxies_phi,backgroundProxies_rap,backgroundProxies_pt,backgroundProxies_mt,backgroundProxies_px_normalized,backgroundProxies_py_normalized,backgroundProxies_pz_normalized;
+ backgroundProxies_phi.reserve(nBackgroundProxies);
+ backgroundProxies_rap.reserve(nBackgroundProxies);
+ backgroundProxies_pt.reserve(nBackgroundProxies);
+ backgroundProxies_mt.reserve(nBackgroundProxies);
+ backgroundProxies_px_normalized.reserve(nBackgroundProxies);
+ backgroundProxies_py_normalized.reserve(nBackgroundProxies);
+ backgroundProxies_pz_normalized.reserve(nBackgroundProxies);
+ for (unsigned int j=0;j backgroundProxies_minParticleIndex,backgroundProxies_maxParticleIndex;
+ if (_use_max_distance && _distance==ConstituentSubtractor::deltaR && _ghosts_rapidity_sorted && !_ghost_selector){
+ for (unsigned int j=0;j<_ghosts_rapidities.size();++j){
+ unsigned int min=this->_find_index_after(_ghosts_rapidities[j]-_max_distance,particles_rap);
+ unsigned int max=this->_find_index_before(_ghosts_rapidities[j]+_max_distance,particles_rap);
+ for (int k=0;k<_n_ghosts_phi;++k){
+ backgroundProxies_minParticleIndex.push_back(min);
+ backgroundProxies_maxParticleIndex.push_back(max);
+ }
+ max_number_pairs+=max-min;
+ }
+ max_number_pairs=max_number_pairs*_n_ghosts_phi*_max_distance/3.1415;
+ }
+ else{
+ for (unsigned int j=0;j > > CS_distances; // storing three elements: the CS distance, and corresponding particle and proxy indexes
+ CS_distances.reserve(max_number_pairs);
+
+ bool skip_particles_outside_phi_range=false; // used for speed optimization
+ if (_distance==ConstituentSubtractor::deltaR && _use_max_distance && _max_distancetwopi){
+ particle_phi_min=particle_phi_max-twopi;
+ particle_phi_max=backgroundProxies_phi[j]-_max_distance;
+ switched=true;
+ }
+ if (particle_phi_min<0){
+ particle_phi_max=particle_phi_min+twopi;
+ particle_phi_min=backgroundProxies_phi[j]+_max_distance;
+ switched=true;
+ }
+ }
+ for (unsigned int i=backgroundProxies_minParticleIndex[j];iparticle_phi_min && particles_phi[i]particle_phi_max))) continue; // speed optimization only, this line has no effect on the subtraction performance
+ double deltaPhi=fabs(backgroundProxies_phi[j]-particles_phi[i]);
+ if (deltaPhi>pi) deltaPhi=twopi-deltaPhi;
+ double deltaRap=backgroundProxies_rap[j]-particles_rap[i];
+ distance_transformed = deltaPhi*deltaPhi+deltaRap*deltaRap;
+ }
+ if (_distance==ConstituentSubtractor::angle){
+ distance_transformed=-(particles_px_normalized[i]*backgroundProxies_px_normalized[j]+particles_py_normalized[i]*backgroundProxies_py_normalized[j]+particles_pz_normalized[i]*backgroundProxies_pz_normalized[j]);
+ }
+ if (!_use_max_distance || distance_transformed<=max_distance_transformed){
+ double CS_distance = distance_transformed*particles_factor[i];
+ CS_distances.push_back(std::make_pair(CS_distance,std::make_pair(i,j))); // have tried to use emplace_back and tuple here - did not lead to any speed improvement
+ }
+ }
+ }
+
+
+ // Sorting of the CS distances
+ std::sort(CS_distances.begin(),CS_distances.end(),ConstituentSubtractor::_function_used_for_sorting);
+
+
+ // The iterative process. Here, only finding the fractions of pt or delta_m=mt-pt to be corrected. The actual correction of particles is done later.
+ unsigned long nStoredPairs=CS_distances.size();
+
+ std::vector backgroundProxies_fraction_of_pt(nBackgroundProxies,1.);
+ std::vector particles_fraction_of_pt(nParticles,1.);
+ std::vector backgroundProxies_fraction_of_mtMinusPt(nBackgroundProxies,1.);
+ std::vector particles_fraction_of_mtMinusPt(nParticles,1.);
+ for (unsigned long iindices=0;iindices0 && particles_fraction_of_pt[particle_index]>0 && particles_pt[particle_index]>0 && backgroundProxies[proxy_index].pt()>0){
+ double ratio_pt=particles_pt[particle_index]*particles_fraction_of_pt[particle_index]/backgroundProxies_pt[proxy_index]/backgroundProxies_fraction_of_pt[proxy_index];
+ if (ratio_pt>1){
+ particles_fraction_of_pt[particle_index]*=1-1./ratio_pt;
+ backgroundProxies_fraction_of_pt[proxy_index]=-1;
+ }
+ else {
+ backgroundProxies_fraction_of_pt[proxy_index]*=1-ratio_pt;
+ particles_fraction_of_pt[particle_index]=-1;
+ }
+ }
+ if (_do_mass_subtraction && backgroundProxies_fraction_of_mtMinusPt[proxy_index]>0 && particles_fraction_of_mtMinusPt[particle_index]>0 && particles_mt[particle_index]>particles_pt[particle_index] && backgroundProxies_mt[proxy_index]>backgroundProxies_pt[proxy_index]){
+ double ratio_mtMinusPt=(particles_mt[particle_index]-particles_pt[particle_index])*particles_fraction_of_mtMinusPt[particle_index]/(backgroundProxies_mt[proxy_index]-backgroundProxies_pt[proxy_index])/backgroundProxies_fraction_of_mtMinusPt[proxy_index];
+ if (ratio_mtMinusPt>1){
+ particles_fraction_of_mtMinusPt[particle_index]*=1-1./ratio_mtMinusPt;
+ backgroundProxies_fraction_of_mtMinusPt[proxy_index]=-1;
+ }
+ else{
+ backgroundProxies_fraction_of_mtMinusPt[proxy_index]*=1-ratio_mtMinusPt;
+ particles_fraction_of_mtMinusPt[particle_index]=-1;
+ }
+ }
+ }
+
+
+ // do the actual correction for particles:
+ std::vector subtracted_particles;
+ PseudoJet subtracted_const;
+ for (unsigned int i=0;i0) corrected_pt=particles_pt[i]*particles_fraction_of_pt[i];
+
+ if (corrected_pt<=ConstituentSubtractorConstants::zero_pt){
+ if (_remove_all_zero_pt_particles) continue;
+ else{
+ corrected_pt=ConstituentSubtractorConstants::zero_pt;
+ particle_pt_larger_than_zero=false;
+ }
+ }
+
+
+ if (_scale_fourmomentum){
+ if (particle_pt_larger_than_zero) subtracted_const=particles_sorted[i]*particles_fraction_of_pt[i];
+ else{
+ double scale=1;
+ if (corrected_ptConstituentSubtractorConstants::zero_mass);
+ if (!particle_pt_larger_than_zero && !particle_mass_larger_than_zero && _remove_particles_with_zero_pt_and_mass) continue;
+
+ double new_mass=ConstituentSubtractorConstants::zero_mass;
+ if (_do_mass_subtraction){
+ if (particles_fraction_of_mtMinusPt[i]>0){
+ double subtracted_mtMinusPt=(particles_mt[i]-particles_pt[i])*particles_fraction_of_mtMinusPt[i];
+ double mass_squared=pow(corrected_pt+subtracted_mtMinusPt,2)-pow(corrected_pt,2);
+ if (mass_squared>0) new_mass=sqrt(mass_squared);
+ }
+ }
+ else if (!_masses_to_zero) new_mass=particles_sorted[i].m();
+
+ if (new_mass<=ConstituentSubtractorConstants::zero_mass){
+ if (!particle_pt_larger_than_zero && _remove_particles_with_zero_pt_and_mass) continue;
+ else new_mass=ConstituentSubtractorConstants::zero_mass;
+ }
+
+
+ // std::cout << "before correction: " << particles_sorted[i].pt() << " " << particles_sorted[i].eta() << " " << particles_sorted[i].rap() << " " << particles_sorted[i].m() << std::endl;
+ if (_fix_pseudorapidity){
+ double scale=corrected_pt/particles_pt[i];
+ subtracted_const.reset(particles_sorted[i].px()*scale,particles_sorted[i].py()*scale,particles_sorted[i].pz()*scale,sqrt(pow(corrected_pt,2)+pow(scale,2)*pow(particles_sorted[i].pz(),2)+pow(new_mass,2)));
+ }
+ else subtracted_const.reset_PtYPhiM(corrected_pt,particles_rap[i],particles_phi[i],new_mass);
+ // std::cout << "after correction: " << subtracted_const.pt() << " " << subtracted_const.eta() << " " << subtracted_const.rap() << " " << subtracted_const.m() << std::endl;
+ }
+ subtracted_const.set_user_index(particles_sorted[i].user_index());
+ subtracted_const.user_info_shared_ptr() = particles_sorted[i].user_info_shared_ptr();
+
+ subtracted_particles.push_back(subtracted_const);
+ }
+
+ // get the remaining background proxies if requested:
+ if (remaining_backgroundProxies){
+ PseudoJet subtracted_proxy;
+ for (unsigned int i=0;i0 && backgroundProxies_pt[i]>0);
+ bool proxy_mtMinusPt_larger_than_zero=(!_masses_to_zero && backgroundProxies_fraction_of_mtMinusPt[i]>0 && backgroundProxies_mt[i]>backgroundProxies_pt[i]);
+
+ double scale=1e-100;
+ if (proxy_pt_larger_than_zero) scale=backgroundProxies_fraction_of_pt[i];
+ if (_scale_fourmomentum) subtracted_proxy=backgroundProxies[i]*scale;
+ else{
+ double new_mass=1e-150;
+ if (proxy_mtMinusPt_larger_than_zero){
+ if (_do_mass_subtraction){
+ double subtracted_mtMinusPt=(backgroundProxies_mt[i]-backgroundProxies_pt[i])*backgroundProxies_fraction_of_mtMinusPt[i];
+ double mass_squared=pow(scale*backgroundProxies_pt[i]+subtracted_mtMinusPt,2)-pow(scale*backgroundProxies_pt[i],2);
+ if (mass_squared>0) new_mass=sqrt(mass_squared);
+ }
+ else new_mass=backgroundProxies[i].m();
+ }
+ if (_fix_pseudorapidity) subtracted_proxy.reset(backgroundProxies[i].px()*scale,backgroundProxies[i].py()*scale,backgroundProxies[i].pz()*scale,sqrt(pow(scale,2)*(backgroundProxies[i].pt2()+pow(backgroundProxies[i].pz(),2))+pow(new_mass,2)));
+ else subtracted_proxy.reset_PtYPhiM(scale*backgroundProxies_pt[i],backgroundProxies_rap[i],backgroundProxies_phi[i],new_mass);
+ }
+ remaining_backgroundProxies->push_back(subtracted_proxy);
+ }
+ }
+
+ return subtracted_particles;
+ }
+
+
+
+ void ConstituentSubtractor::set_keep_original_masses(){
+ _masses_to_zero=false;
+ }
+
+
+
+
+ std::vector ConstituentSubtractor::subtract_event_using_charged_info(std::vector const &particles, double charged_background_scale, std::vector const &charged_background, double charged_signal_scale, std::vector const &charged_signal, double max_eta){
+ if (fabs(_max_eta/max_eta-1)>1e-5) _ghosts_constructed=false;
+ if (!_ghosts_constructed) this->construct_ghosts_uniformly(max_eta);
+ _ghosts_rapidity_sorted=false; // no speed optimization implemented for this function yet
+
+ std::vector scaled_charged_all;
+ std::vector scaled_charged_signal;
+ std::vector scaled_charged_background;
+ for (unsigned int i=0;imax_eta) continue;
+ scaled_charged_background.push_back(charged_background_scale*charged_background[i]);
+ scaled_charged_all.push_back(scaled_charged_background[scaled_charged_background.size()-1]);
+ }
+ for (unsigned int i=0;imax_eta) continue;
+ scaled_charged_signal.push_back(charged_signal_scale*charged_signal[i]);
+ scaled_charged_all.push_back(scaled_charged_signal[scaled_charged_signal.size()-1]);
+ }
+ std::vector selected_particles;
+ for (unsigned int i=0;i *remaining_charged_background= new std::vector;
+ double maxDeltaR=this->get_max_distance();
+ if (maxDeltaR<=0) maxDeltaR=0.5;
+ this->set_max_distance(0.2);
+ std::vector subtracted_particles_using_scaled_charged_signal=this->do_subtraction(selected_particles,scaled_charged_signal);
+ std::vector subtracted_particles_using_scaled_charged_all=this->do_subtraction(subtracted_particles_using_scaled_charged_signal,scaled_charged_background,remaining_charged_background); // remaining neutral background particles
+ std::vector scaled_charged_background_used_for_subtraction=this->do_subtraction(scaled_charged_background,*remaining_charged_background);
+ _bge_rho= new fastjet::GridMedianBackgroundEstimator(max_eta, _grid_size_background_estimator);
+ if (_do_mass_subtraction) this->set_common_bge_for_rho_and_rhom();
+ _bge_rho->set_rescaling_class(_rescaling);
+ _bge_rho->set_particles(subtracted_particles_using_scaled_charged_all);
+
+ std::vector backgroundProxies=this->get_background_proxies_from_ghosts(_ghosts,_ghosts_area);
+ backgroundProxies.insert(backgroundProxies.end(), scaled_charged_background_used_for_subtraction.begin(), scaled_charged_background_used_for_subtraction.end());
+
+ this->set_max_distance(maxDeltaR);
+ std::vector subtracted_particles=this->do_subtraction(selected_particles,backgroundProxies);
+ delete remaining_charged_background;
+ delete _bge_rho;
+ return subtracted_particles;
+ }
+
+
+
+ void ConstituentSubtractor::set_grid_size_background_estimator(double const &grid_size_background_estimator){
+ _grid_size_background_estimator=grid_size_background_estimator;
+ }
+
+
+
+ void ConstituentSubtractor::set_common_bge_for_rho_and_rhom(bool value){
+ if (value) this->set_common_bge_for_rho_and_rhom();
+ else throw Error("ConstituentSubtractor::set_common_bge_for_rho_and_rhom: This function should be not used with false! Read the instructions for mass subtraction in the header file.");
+ }
+
+
+ void ConstituentSubtractor::set_common_bge_for_rho_and_rhom(){
+ if (!_bge_rho) throw Error("ConstituentSubtractor::set_common_bge_for_rho_and_rhom() is not allowed when _bge_rho is not set!");
+ if (_bge_rhom) throw Error("ConstituentSubtractor::set_common_bge_for_rho_and_rhom() is not allowed in the presence of an existing background estimator for rho_m.");
+ if (_externally_supplied_rho_rhom) throw Error("ConstituentSubtractor::set_common_bge_for_rho_and_rhom() is not allowed when supplying externally the values for rho and rho_m.");
+
+#if FASTJET_VERSION_NUMBER >= 30100
+ if (!_bge_rho->has_rho_m()){
+#endif
+ JetMedianBackgroundEstimator *jmbge = dynamic_cast(_bge_rho);
+ if (!jmbge) throw Error("ConstituentSubtractor::set_common_bge_for_rho_and_rhom() is currently only allowed for background estimators of JetMedianBackgroundEstimator type.");
+#if FASTJET_VERSION_NUMBER >= 30100
+ }
+#endif
+ _common_bge=true;
+ }
+
+
+ // setting this to true will result in rho_m being estimated using bge_rhom->rho_m() instead of bge_rhom->rho()
+ /* void ConstituentSubtractor::set_use_bge_rhom_rhom(bool value){
+ if (!value){
+ _rhom_from_bge_rhom=false;
+ return;
+ }
+
+#if FASTJET_VERSION_NUMBER < 30100
+ throw Error("ConnstituentSubtractor::use_rhom_from_bge_rhom() can only be used with FastJet >=3.1.");
+#else
+ if (!_bge_rhom) throw Error("ConstituentSubtractor::use_rhom_from_bge_rhom() requires a background estimator for rho_m.");
+
+ if (!(_bge_rhom->has_rho_m())) throw Error("ConstituentSubtractor::use_rhom_from_bge_rhom() requires rho_m support for the background estimator for rho_m.");
+#endif
+ _rhom_from_bge_rhom=true;
+ _do_mass_subtraction=true;
+ }*/
+
+
+ void ConstituentSubtractor::set_do_mass_subtraction(){
+ _do_mass_subtraction=true;
+ _masses_to_zero=false;
+ }
+
+
+
+ void ConstituentSubtractor::description_common(std::ostringstream &descr) const{
+ if ( _externally_supplied_rho_rhom){
+ descr << " Using externally supplied rho = " << _rho << " and rho_m = " << _rhom << std::endl;
+ } else {
+ if (_bge_rhom && _bge_rho) {
+ descr << " Using rho estimation: " << _bge_rho->description() << std::endl;
+ descr << " Using rho_m estimation: " << _bge_rhom->description() << std::endl;
+ } else {
+ if (_bge_rho) descr << " Using rho estimation: " << _bge_rho->description() << std::endl;
+ else descr << " No externally supplied rho, nor background estimator" << std::endl;
+ }
+ }
+
+ if (_do_mass_subtraction){
+ descr << " The mass part (delta_m) will be also corrected." << std::endl;
+ if (_common_bge) descr << " using the same background estimator for rho_m as for rho" << std::endl;
+ else descr << " using different background estimator for rho_m as for rho" << std::endl;
+ }
+ else if (_masses_to_zero) descr << " The masses of all particles will be set to zero." << std::endl;
+ else if (_scale_fourmomentum) descr << " The masses will be corrected by scaling the whole 4-momentum." << std::endl;
+ else descr << " The original mass of the particles will be kept." << std::endl;
+
+ if (!_scale_fourmomentum){
+ if (_fix_pseudorapidity) descr << " The pseudo-rapidity of the particles will be kept unchanged (not rapidity)." << std::endl;
+ else descr << " The rapidity of the particles will be kept unchanged (not pseudo-rapidity)." << std::endl;
+ }
+
+ if (_use_nearby_hard) descr << " Using information about nearby hard proxies with parameters _nearby_hard_radius=" << _nearby_hard_radius << " and _nearby_hard_factor=" << _nearby_hard_factor << std::endl;
+ else descr << " The information about nearby hard proxies will not be used." << std::endl;
+ }
+
+
+
+
+
+ std::string ConstituentSubtractor::description() const{
+ std::ostringstream descr;
+ descr << std::endl << "Description of fastjet::ConstituentSubtractor which can be used for event-wide or jet-by-jet correction:" << std::endl;
+ this->description_common(descr);
+ descr << " Using parameters: max_distance = " << _max_distance << " alpha = " << _alpha << std::endl;
+ return descr.str();
+ }
+
+
+
+ void ConstituentSubtractor::set_distance_type(Distance distance){
+ _distance=distance;
+ }
+
+
+ void ConstituentSubtractor::set_max_distance(double max_distance){
+ if (max_distance>0){
+ _use_max_distance=true;
+ _max_distance=max_distance;
+ }
+ else _use_max_distance=false;
+ }
+
+
+
+ void ConstituentSubtractor::set_max_standardDeltaR(double max_distance){
+ this->set_max_distance(max_distance);
+ }
+
+
+
+ double ConstituentSubtractor::get_max_distance(){
+ return _max_distance;
+ }
+
+
+ void ConstituentSubtractor::set_alpha(double alpha){
+ _alpha=alpha;
+ }
+
+ void ConstituentSubtractor::set_polarAngleExp(double polarAngleExp){
+ _polarAngleExp=polarAngleExp;
+ }
+
+
+ void ConstituentSubtractor::set_ghost_area(double ghost_area){
+ _ghost_area=ghost_area;
+ this->clear_ghosts();
+ }
+
+
+ void ConstituentSubtractor::set_max_eta(double max_eta){
+ _max_eta=max_eta;
+ }
+
+
+ void ConstituentSubtractor::set_fix_pseudorapidity(){
+ _fix_pseudorapidity=true;
+ }
+
+
+ void ConstituentSubtractor::set_scale_fourmomentum(){
+ _scale_fourmomentum=true;
+ _masses_to_zero=false;
+ }
+
+
+ void ConstituentSubtractor::set_remove_particles_with_zero_pt_and_mass(bool value){
+ _remove_particles_with_zero_pt_and_mass=value;
+ }
+
+
+ void ConstituentSubtractor::set_remove_all_zero_pt_particles(bool value){
+ _remove_all_zero_pt_particles=value;
+ }
+
+
+ bool ConstituentSubtractor::_function_used_for_sorting(std::pair > const &i,std::pair > const &j){
+ return (i.first < j.first);
+ }
+
+ bool ConstituentSubtractor::_rapidity_sorting(fastjet::PseudoJet const &i,fastjet::PseudoJet const &j){
+ return (i.rap() < j.rap());
+ }
+
+ unsigned int ConstituentSubtractor::_find_index_after(double const &value, std::vector const &vec) const{
+ int size=vec.size();
+ if (size==0) return -1;
+ int nIterations=log(size)/log(2)+2;
+ unsigned int lowerBound=0;
+ unsigned int upperBound=size-1;
+ if (value<=vec[0]) return 0;
+ if (value>vec[size-1]) return size; // if the value is larger than the maximal possible rapidity, than setting min to size - then also max is set to size, and nothing will be run in the loop
+ for (int i=0;ivec[test]){
+ if (value<=vec[test+1]) return test+1;
+ lowerBound=test;
+ }
+ else{
+ if (value>vec[test-1]) return test;
+ upperBound=test;
+ }
+ }
+ return lowerBound;
+ }
+
+
+ unsigned int ConstituentSubtractor::_find_index_before(double const &value, std::vector const &vec) const{
+ int size=vec.size();
+ if (size==0) return -1;
+ int nIterations=log(size)/log(2)+1;
+ unsigned int lowerBound=0;
+ unsigned int upperBound=size-1;
+ if (value=vec[size-1]) return size; // it is higher by one to account for the "<" comparison in the for loop
+ for (int i=0;i=vec[test]){
+ if (value=vec[test-1]) return test; // it is higher by one to account for the "<" comparison in the for loop
+ upperBound=test;
+ }
+ }
+ return upperBound+1;
+ }
+
+
+
+
+ void ConstituentSubtractor::clear_ghosts(){
+ _ghosts.clear();
+ _ghosts_rapidities.clear();
+ _ghosts_area.clear();
+ _ghosts_rapidity_sorted=false;
+ _ghosts_constructed=false;
+ }
+
+
+ void ConstituentSubtractor::construct_ghosts_uniformly(double max_eta){
+ this->clear_ghosts();
+ _max_eta=max_eta;
+ double a=sqrt(_ghost_area);
+ _n_ghosts_phi=(2*3.14159265/a+0.5); // rounding
+ int n_ghosts_rap=(2*max_eta/a+0.5); // rounding
+ _grid_size_phi=2*3.14159265/(double)_n_ghosts_phi;
+ _grid_size_rap=2*max_eta/(double)n_ghosts_rap;
+ double used_ghost_area=_grid_size_phi*_grid_size_rap;
+ fastjet::PseudoJet ghost(0,0,0,1);
+ for (int iRap=0;iRappass(ghost)) continue;
+ }
+ _ghosts.push_back(ghost);
+ _ghosts_area.push_back(used_ghost_area);
+ }
+ }
+ _ghosts_rapidity_sorted=true; // the ghosts are now sorted according to rapidity. This variable needs to be set to true to be able to use faster algorithm in "do_subtraction".
+ _ghosts_constructed=true;
+ }
+
+
+ std::vector ConstituentSubtractor::get_ghosts(){
+ return _ghosts;
+ }
+
+
+ std::vector ConstituentSubtractor::get_ghosts_area(){
+ return _ghosts_area;
+ }
+
+
+
+ void ConstituentSubtractor::set_ghost_selector(fastjet::Selector* selector){
+ _ghost_selector=selector;
+ this->clear_ghosts();
+ }
+
+ void ConstituentSubtractor::set_particle_selector(fastjet::Selector* selector){
+ _particle_selector=selector;
+ }
+
+ void ConstituentSubtractor::set_rescaling(fastjet::FunctionOfPseudoJet *rescaling){
+ _rescaling=rescaling;
+ }
+
+
+ void ConstituentSubtractor::set_use_nearby_hard(double const &nearby_hard_radius, double const &nearby_hard_factor){
+ _nearby_hard_radius=nearby_hard_radius;
+ _nearby_hard_factor=nearby_hard_factor;
+ if (_nearby_hard_radius>0) _use_nearby_hard=true;
+ else _use_nearby_hard=false;
+ }
+
+} // namespace contrib
+
+
+FASTJET_END_NAMESPACE
Property changes on: contrib/contribs/ConstituentSubtractor/tags/1.4.7/ConstituentSubtractor.cc
___________________________________________________________________
Added: svn:keywords
## -0,0 +1 ##
+Id
\ No newline at end of property
Index: contrib/contribs/ConstituentSubtractor/tags/1.4.7/NEWS
===================================================================
--- contrib/contribs/ConstituentSubtractor/tags/1.4.7/NEWS (revision 0)
+++ contrib/contribs/ConstituentSubtractor/tags/1.4.7/NEWS (revision 1373)
@@ -0,0 +1,123 @@
+2024/02/01: release of version 1.4.7:
+ - Lowering the number of printed digits in example example_whole_event_using_charged_info.cc to avoid differences when running on Mac with Apple clang 15.0.0 and an M2 (Arm) chip (reported by Gavin Salam on 2024/01/21)
+ - Adding the possibility to use linear interpolation for the rescaling functions defined in RescalingClasses. Using interpolation for rescaling class BackgroundRescalingYPhiUsingVectorForY by default. This should avoid to have differences for example example_background_rescaling.cc when running on Mac with Apple clang 15.0.0 and an M2 (Arm) chip (reported by Gavin Salam on 2024/01/21)
+ - Updating the description in example example_background_rescaling.cc to make clear the possibility to use interpolation.
+
+
+2023/05/11: release of version 1.4.6:
+ - Fixed a bug for jet-by-jet Constituent Subtractor when using particle selector
+ - Copying also user_info for the corrected jet and corrected constituents
+ - Initializing previously not initialized variables (caused a FPE warning in ATLAS)
+
+2020/02/23: release of version 1.4.5:
+ - Further removal of C++11 extensions which were missed previously
+
+
+2019/07/15: release of version 1.4.4:
+ - Removed two more C++11 extensions
+
+
+2019/07/09: release of version 1.4.3:
+ - Removed C++11 extensions
+ - Changed the way how the user should specify that he/she wants to do delta_m correction. Now, additionally the "set_do_mass_correction" function must be called
+ - Updated information about the options for the treatment of massive inputs.
+
+
+
+2019/07/04: release of version 1.4.2:
+ - Updated information about the treatment of massive inputs. No changes in the code.
+
+
+2019/06/21: release of version 1.4.1:
+ - Updated the treatment of additional background estimator for rho_m. Currently in ConstituentSubtractor, only such rho_m background estimators are supported which have member function rho_m(), i.e. the result from member function has_rho_m() is true (which is available since fastjet version 30100).
+
+
+2019/06/19: release of version 1.4.0:
+ - Updated naming of functions and examples to be consistent with arXiv:1905.03470
+ - Updated treatment of massive inputs - added new options and clarification in the examples
+ - Updated values of pt and mass for particles with zero pt and mass
+ - Added possibility to use selector for particles - only particles passing the selector are corrected, the other particles are unchanged
+ - Updated COPYING
+
+
+2019/03/28: release of version 1.3.2:
+ - Started to use C++11 features
+ - Added possibility to keep the masses unchanged. The possible options how to treat masses are explained at the beggining of ConstituentSubtractor.hh
+ - Added possibility to use information about hard particles (e.g. charged tracks from primary vertex) in the definition of the CS distance used for the correction procedure.
+ - Updated examples
+
+
+2018/06/28: release of version 1.3.1:
+ - Updated examples
+
+
+2018/06/25: release of version 1.3.0:
+ - Replaced the "sequential" CS by "iterative CS". Besides the change of name, also the usage changed, see example file example_whole_event_iterative.cc.
+ - The usage of the ConstituentSubtractor for whole event correction is changes, see updated example file example_whole_event.cc. Currently the new and old usage give the same behaviour, but it will change in the future.
+
+
+2018/06/12: release of version 1.2.8:
+ - Now removing all ghosts when setting new selector
+ - Updated sequential CS - possibility to do more iterations than 2 and possibility to remove ghosts which were not used.
+
+2018/05/04: release of version 1.2.7:
+ - Added possibility to correct only certain phase space defined by an input fastjet::Selector
+ - Added possibility to change the grid size for background estimation when using function subtract_event_using_charged_info
+ - Update example for whole event correction.
+
+
+2018/05/03: release of version 1.2.6:
+ - Updated examples to be more clear.
+
+
+2018/04/11: release of version 1.2.5:
+ - Fixed bug causing seg fault in ATLAS.
+
+
+2017/12/18: release of version 1.2.4:
+ - Added back the test example_whole_event_using_charged_info. The problem with this test was fixed by printing all the output particles (also with zero pt and zeo mtMinusPt values).
+
+
+2017/09/08: release of version 1.2.3:
+ - Fixed warning for too low pt particles in case of active area computation
+ - Temporarily removed the test example_whole_event_using_charged_info
+
+2017/08/21: release of version 1.2.2:
+ - Added function set_max_standardDeltaR for backward compatibility. It was replaced in version 1.1.3 by function set_max_distance. In this version, it does the same as the function set_max_distance.
+ - Further speed improvement for whole event subtraction (especially for small _max_distance parameters, such as 0.25)
+
+
+2017/08/07: release of version 1.2.1:
+ - Fixed bug for distance enum ConstituentSubtractor::angle
+ - Added check for the jet-by-jet correction to not use the speed optimization
+
+2017/07/20: release of version 1.2.0:
+ - Important change for the output corrected particles in case of correction of massive particles. The corrected zero pt particles with non-zero corrected delta_m (massive particles in rest) are now included in the output by default. This is the recommended usage. The user can change this behaviour by using function "set_remove_all_zero_pt_particles". By using argument true, all the zero pt corrected particles are removed.
+ - Significant change in the speed of the whole event CS procedure. For small values of the _max_distance parameter (~0.25), the speed is now ~10-times greater.
+ - Added new rescaling classes: rapidity vs azimuth dependence in TH2 Root histogram or in vector >, and rapidity vs azimuth dependence for heavy ions where the rapidity dependence is taken from a vector.
+
+2017/07/11: release of version 1.1.6:
+ - added new functions "get_ghosts" and "get_ghosts_area".
+ - now the user can modify the ghosts outside the CS after constructing them with function "construct_ghosts_uniformly" and getting them with "get_ghosts" and "get_ghosts_area"
+ - added new rescaling class useful for heavy ion events. Now the rapidity dependence can be taken from a Root TH1 object.
+
+2017/01/20: release of version 1.1.4:
+ - fixed few bugs concerning the correction of masses
+
+2016/11/14: release of version 1.1.3:
+ - updated rho rescaling functions
+ - new function for pileup subtraction: sequential_subtraction
+ - possibility to change the computation of distance between particles and background proxies
+
+2016/03/11: release of version 1.1.2:
+ - added two FunctionOfPseudojet