ClusterSequenceAreaBase.cc

Go to the documentation of this file.

//STARTHEADER
// $Id: ClusterSequenceAreaBase.cc 2281 2011-06-26 08:32:47Z salam $
//
// Copyright (c) 2005-2006, Matteo Cacciari and Gavin Salam
//
//----------------------------------------------------------------------
// This file is part of FastJet.
//
//  FastJet 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.
//
//  The algorithms that underlie FastJet have required considerable
//  development and are described in hep-ph/0512210. If you use
//  FastJet as part of work towards a scientific publication, please
//  include a citation to the FastJet paper.
//
//  FastJet 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 FastJet; if not, write to the Free Software
//  Foundation, Inc.:
//      59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//----------------------------------------------------------------------
//ENDHEADER




#include "fastjet/ClusterSequenceAreaBase.hh"
#include <algorithm>

FASTJET_BEGIN_NAMESPACE

using namespace std;


LimitedWarning ClusterSequenceAreaBase::_warnings;
LimitedWarning ClusterSequenceAreaBase::_warnings_zero_area;
LimitedWarning ClusterSequenceAreaBase::_warnings_empty_area;

//----------------------------------------------------------------------
double ClusterSequenceAreaBase::empty_area(const RangeDefinition & range) const {

  if (has_explicit_ghosts()) {return 0.0;}
  else { return empty_area_from_jets(inclusive_jets(0.0), range);}

}

//----------------------------------------------------------------------
double ClusterSequenceAreaBase::empty_area_from_jets(
                      const std::vector<PseudoJet> & all_jets,
                      const RangeDefinition & range) const {

  double empty = range.area();
  for (unsigned i = 0; i < all_jets.size(); i++) {
    if (range.is_in_range(all_jets[i])) empty -= area(all_jets[i]);
  }
  return empty;
}

double ClusterSequenceAreaBase::median_pt_per_unit_area(const RangeDefinition & range) const {
  return median_pt_per_unit_something(range,false);
}

double ClusterSequenceAreaBase::median_pt_per_unit_area_4vector(const RangeDefinition & range) const {
  return median_pt_per_unit_something(range,true);
}


//----------------------------------------------------------------------
double ClusterSequenceAreaBase::median_pt_per_unit_something(
                const RangeDefinition & range, bool use_area_4vector) const {

  double median, sigma, mean_area;
  get_median_rho_and_sigma(range, use_area_4vector, median, sigma, mean_area);
  return median;

}


//----------------------------------------------------------------------
void ClusterSequenceAreaBase::parabolic_pt_per_unit_area(
       double & a, double & b, const RangeDefinition & range, 
       double exclude_above, bool use_area_4vector) const {
  
  int n=0;
  int n_excluded = 0;
  double mean_f=0, mean_x2=0, mean_x4=0, mean_fx2=0; 

  vector<PseudoJet> incl_jets = inclusive_jets();

  for (unsigned i = 0; i < incl_jets.size(); i++) {
    if (range.is_in_range(incl_jets[i])) {
      double this_area;
      if ( use_area_4vector ) {
          this_area = area_4vector(incl_jets[i]).perp();     
      } else {
          this_area = area(incl_jets[i]);
      }
      double f = incl_jets[i].perp()/this_area;
      if (exclude_above <= 0.0 || f < exclude_above) {
        double x = incl_jets[i].rap(); double x2 = x*x;
        mean_f   += f;
        mean_x2  += x2;
        mean_x4  += x2*x2;
        mean_fx2 += f*x2;
        n++;
      } else {
        n_excluded++;
      }
    }
  }

  if (n <= 1) {
    // meaningful results require at least two jets inside the
    // area -- mind you if there are empty jets we should be in 
    // any case doing something special...
    a = 0.0;
    b = 0.0;
  } else {
    mean_f   /= n;
    mean_x2  /= n;
    mean_x4  /= n;
    mean_fx2 /= n;
    
    b = (mean_f*mean_x2 - mean_fx2)/(mean_x2*mean_x2 - mean_x4);
    a = mean_f - b*mean_x2;
  }
  //cerr << "n_excluded = "<< n_excluded << endl;
}



void ClusterSequenceAreaBase::get_median_rho_and_sigma(
            const RangeDefinition & range, bool use_area_4vector,
            double & median, double & sigma, double & mean_area) const {

  vector<PseudoJet> incl_jets = inclusive_jets();
  get_median_rho_and_sigma(incl_jets, range, use_area_4vector,
                           median, sigma, mean_area, true);
}


void ClusterSequenceAreaBase::get_median_rho_and_sigma(
            const vector<PseudoJet> & all_jets,
            const RangeDefinition & range, bool use_area_4vector,
            double & median, double & sigma, double & mean_area,
            bool all_are_incl) const {

  _check_jet_alg_good_for_median();

  vector<double> pt_over_areas;
  double total_area  = 0.0;
  double total_njets = 0;

  for (unsigned i = 0; i < all_jets.size(); i++) {
    if (range.is_in_range(all_jets[i])) {
      double this_area;
      if (use_area_4vector) {
          this_area = area_4vector(all_jets[i]).perp();
      } else {
          this_area = area(all_jets[i]);
      }

      if (this_area>0) {
        pt_over_areas.push_back(all_jets[i].perp()/this_area);
      } else {
        _warnings_zero_area.warn("ClusterSequenceAreaBase::get_median_rho_and_sigma(...): discarded jet with zero area. Zero-area jets may be due to (i) too large a ghost area (ii) a jet being outside the ghost range (iii) the computation not being done using an appropriate algorithm (kt;C/A).");
      }

      total_area  += this_area;
      total_njets += 1.0;
    }
  }

  // there is nothing inside our region, so answer will always be zero
  if (pt_over_areas.size() == 0) {
    median = 0.0;
    sigma  = 0.0;
    mean_area = 0.0;
    return;
  }
  
  // get median (pt/area) [this is the "old" median definition. It considers
  // only the "real" jets in calculating the median, i.e. excluding the
  // only-ghost ones; it will be supplemented with more info below]
  sort(pt_over_areas.begin(), pt_over_areas.end());

  // now get the median & error, accounting for empty jets
  // define the fractions of distribution at median, median-1sigma
  double posn[2] = {0.5, (1.0-0.6827)/2.0};
  double res[2];
  
  double n_empty, empty_a;
  if (has_explicit_ghosts()) {
    // NB: the following lines of code are potentially incorrect in cases
    //     where there are unclustered particles (empty_area would do a better job,
    //     at least for active areas). This is not an issue with kt or C/A, or other
    //     algorithms that cluster all particles (and the median estimation should in 
    //     any case only be done with kt or C/A!)
    empty_a = 0.0;
    n_empty = 0;
  } else if (all_are_incl) {
    // the default case
    empty_a = empty_area(range);
    n_empty = n_empty_jets(range);
  } else {
    // this one is intended to be used when e.g. one runs C/A, then looks at its
    // exclusive jets in order to get an effective smaller R value, and passes those
    // to this routine.
    empty_a = empty_area_from_jets(all_jets, range);
    mean_area = total_area / total_njets; // temporary value
    n_empty   = empty_a / mean_area;
  }
  //cout << "*** tot_area = " << total_area << ", empty_a = " << empty_a << endl;
  //cout << "*** n_empty = " << n_empty << ", ntotal =  " << total_njets << endl;
  total_njets += n_empty;
  total_area  += empty_a;

  // we need an int (rather than an unsigned int) with the size of the
  // pt_over_areas array, because we'll often be doing subtraction of
  // -1, negating it, etc. All of these operations go crazy with unsigned ints.
  int pt_over_areas_size = pt_over_areas.size();
  if (n_empty < -pt_over_areas_size/4)
    _warnings_empty_area.warn("ClusterSequenceAreaBase::get_median_rho_and_sigma(...): the estimated empty area is suspiciously large and may lead to an over-estimation of rho. This may be due to (i) a rare statistical fluctuation or (ii) too small a range used to estimate the background properties.");

  for (int i = 0; i < 2; i++) {
    double nj_median_pos = 
      (pt_over_areas_size-1.0 + n_empty)*posn[i] - n_empty;
    double nj_median_ratio;
    if (nj_median_pos >= 0 && pt_over_areas_size > 1) {
      int int_nj_median = int(nj_median_pos);

      // avoid potential overflow issues
      if (int_nj_median+1 > pt_over_areas_size-1){
        int_nj_median = pt_over_areas_size-2;
        nj_median_pos = pt_over_areas_size-1;
      }

      nj_median_ratio = 
        pt_over_areas[int_nj_median] * (int_nj_median+1-nj_median_pos)
        + pt_over_areas[int_nj_median+1] * (nj_median_pos - int_nj_median);
    } else {
      nj_median_ratio = 0.0;
    }
    res[i] = nj_median_ratio;
  }
  median = res[0];
  double error  = res[0] - res[1];
  mean_area = total_area / total_njets;
  sigma  = error * sqrt(mean_area);
}


vector<PseudoJet> ClusterSequenceAreaBase::subtracted_jets(const double rho,
                                                           const double ptmin) 
                                                           const {
  vector<PseudoJet> sub_jets;
  vector<PseudoJet> jets = sorted_by_pt(inclusive_jets(ptmin));
  for (unsigned i=0; i<jets.size(); i++) {
     PseudoJet sub_jet = subtracted_jet(jets[i],rho);
     sub_jets.push_back(sub_jet);
  }
  return sub_jets;
}

vector<PseudoJet> ClusterSequenceAreaBase::subtracted_jets(
                                                 const RangeDefinition & range, 
                                                 const double ptmin)
                                                 const {
  double rho = median_pt_per_unit_area_4vector(range);
  return subtracted_jets(rho,ptmin);
}


PseudoJet ClusterSequenceAreaBase::subtracted_jet(const PseudoJet & jet,
                                                  const double rho) const {
  PseudoJet area4vect = area_4vector(jet);
  PseudoJet sub_jet;
  // sanity check
  if (rho*area4vect.perp() < jet.perp() ) { 
    sub_jet = jet - rho*area4vect;
  } else { sub_jet = PseudoJet(0.0,0.0,0.0,0.0); }
  
  // make sure the subtracted jet has the same index (cluster and user)
  // (i.e. "looks like") the original jet
  sub_jet.set_cluster_hist_index(jet.cluster_hist_index());
  sub_jet.set_user_index(jet.user_index());
  
  return sub_jet;
}


PseudoJet ClusterSequenceAreaBase::subtracted_jet(const PseudoJet & jet,
                                       const RangeDefinition & range) const {
  double rho = median_pt_per_unit_area_4vector(range);
  PseudoJet sub_jet = subtracted_jet(jet, rho);
  return sub_jet;
}


double ClusterSequenceAreaBase::subtracted_pt(const PseudoJet & jet,
                                              const double rho,
                                              bool use_area_4vector) const {
  if ( use_area_4vector ) { 
     PseudoJet sub_jet = subtracted_jet(jet,rho);
     return sub_jet.perp();
  } else {
     return jet.perp() - rho*area(jet);
  }
}  


double ClusterSequenceAreaBase::subtracted_pt(const PseudoJet & jet,
                                              const RangeDefinition & range,
                                              bool use_area_4vector) const {
  if ( use_area_4vector ) { 
     PseudoJet sub_jet = subtracted_jet(jet,range);
     return sub_jet.perp();
  } else {
     double rho = median_pt_per_unit_area(range);
     return subtracted_pt(jet,rho,false);
  }
}  


void ClusterSequenceAreaBase::_check_jet_alg_good_for_median() const {
  if (jet_def().jet_algorithm() != kt_algorithm
      && jet_def().jet_algorithm() != cambridge_algorithm
      && jet_def().jet_algorithm() !=  cambridge_for_passive_algorithm) {
    _warnings.warn("ClusterSequenceAreaBase: jet_def being used may not be suitable for estimating diffuse backgrounds (good options are kt, cam)");
  }
}



FASTJET_END_NAMESPACE