fastjet 2.4.5
ClusterSequenceAreaBase.cc
Go to the documentation of this file.
00001 
00002 //STARTHEADER
00003 // $Id: ClusterSequenceAreaBase.cc 2281 2011-06-26 08:32:47Z salam $
00004 //
00005 // Copyright (c) 2005-2006, Matteo Cacciari and Gavin Salam
00006 //
00007 //----------------------------------------------------------------------
00008 // This file is part of FastJet.
00009 //
00010 //  FastJet is free software; you can redistribute it and/or modify
00011 //  it under the terms of the GNU General Public License as published by
00012 //  the Free Software Foundation; either version 2 of the License, or
00013 //  (at your option) any later version.
00014 //
00015 //  The algorithms that underlie FastJet have required considerable
00016 //  development and are described in hep-ph/0512210. If you use
00017 //  FastJet as part of work towards a scientific publication, please
00018 //  include a citation to the FastJet paper.
00019 //
00020 //  FastJet is distributed in the hope that it will be useful,
00021 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00022 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00023 //  GNU General Public License for more details.
00024 //
00025 //  You should have received a copy of the GNU General Public License
00026 //  along with FastJet; if not, write to the Free Software
00027 //  Foundation, Inc.:
00028 //      59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00029 //----------------------------------------------------------------------
00030 //ENDHEADER
00031 
00032 
00033 
00034 
00035 #include "fastjet/ClusterSequenceAreaBase.hh"
00036 #include <algorithm>
00037 
00038 FASTJET_BEGIN_NAMESPACE
00039 
00040 using namespace std;
00041 
00042 
00044 LimitedWarning ClusterSequenceAreaBase::_warnings;
00045 LimitedWarning ClusterSequenceAreaBase::_warnings_zero_area;
00046 LimitedWarning ClusterSequenceAreaBase::_warnings_empty_area;
00047 
00048 //----------------------------------------------------------------------
00056 double ClusterSequenceAreaBase::empty_area(const RangeDefinition & range) const {
00057 
00058   if (has_explicit_ghosts()) {return 0.0;}
00059   else { return empty_area_from_jets(inclusive_jets(0.0), range);}
00060 
00061 }
00062 
00063 //----------------------------------------------------------------------
00068 double ClusterSequenceAreaBase::empty_area_from_jets(
00069                       const std::vector<PseudoJet> & all_jets,
00070                       const RangeDefinition & range) const {
00071 
00072   double empty = range.area();
00073   for (unsigned i = 0; i < all_jets.size(); i++) {
00074     if (range.is_in_range(all_jets[i])) empty -= area(all_jets[i]);
00075   }
00076   return empty;
00077 }
00078 
00079 double ClusterSequenceAreaBase::median_pt_per_unit_area(const RangeDefinition & range) const {
00080   return median_pt_per_unit_something(range,false);
00081 }
00082 
00083 double ClusterSequenceAreaBase::median_pt_per_unit_area_4vector(const RangeDefinition & range) const {
00084   return median_pt_per_unit_something(range,true);
00085 }
00086 
00087 
00088 //----------------------------------------------------------------------
00092 double ClusterSequenceAreaBase::median_pt_per_unit_something(
00093                 const RangeDefinition & range, bool use_area_4vector) const {
00094 
00095   double median, sigma, mean_area;
00096   get_median_rho_and_sigma(range, use_area_4vector, median, sigma, mean_area);
00097   return median;
00098 
00099 }
00100 
00101 
00102 //----------------------------------------------------------------------
00106 void ClusterSequenceAreaBase::parabolic_pt_per_unit_area(
00107        double & a, double & b, const RangeDefinition & range, 
00108        double exclude_above, bool use_area_4vector) const {
00109   
00110   int n=0;
00111   int n_excluded = 0;
00112   double mean_f=0, mean_x2=0, mean_x4=0, mean_fx2=0; 
00113 
00114   vector<PseudoJet> incl_jets = inclusive_jets();
00115 
00116   for (unsigned i = 0; i < incl_jets.size(); i++) {
00117     if (range.is_in_range(incl_jets[i])) {
00118       double this_area;
00119       if ( use_area_4vector ) {
00120           this_area = area_4vector(incl_jets[i]).perp();     
00121       } else {
00122           this_area = area(incl_jets[i]);
00123       }
00124       double f = incl_jets[i].perp()/this_area;
00125       if (exclude_above <= 0.0 || f < exclude_above) {
00126         double x = incl_jets[i].rap(); double x2 = x*x;
00127         mean_f   += f;
00128         mean_x2  += x2;
00129         mean_x4  += x2*x2;
00130         mean_fx2 += f*x2;
00131         n++;
00132       } else {
00133         n_excluded++;
00134       }
00135     }
00136   }
00137 
00138   if (n <= 1) {
00139     // meaningful results require at least two jets inside the
00140     // area -- mind you if there are empty jets we should be in 
00141     // any case doing something special...
00142     a = 0.0;
00143     b = 0.0;
00144   } else {
00145     mean_f   /= n;
00146     mean_x2  /= n;
00147     mean_x4  /= n;
00148     mean_fx2 /= n;
00149     
00150     b = (mean_f*mean_x2 - mean_fx2)/(mean_x2*mean_x2 - mean_x4);
00151     a = mean_f - b*mean_x2;
00152   }
00153   //cerr << "n_excluded = "<< n_excluded << endl;
00154 }
00155 
00156 
00157 
00158 void ClusterSequenceAreaBase::get_median_rho_and_sigma(
00159             const RangeDefinition & range, bool use_area_4vector,
00160             double & median, double & sigma, double & mean_area) const {
00161 
00162   vector<PseudoJet> incl_jets = inclusive_jets();
00163   get_median_rho_and_sigma(incl_jets, range, use_area_4vector,
00164                            median, sigma, mean_area, true);
00165 }
00166 
00167 
00168 void ClusterSequenceAreaBase::get_median_rho_and_sigma(
00169             const vector<PseudoJet> & all_jets,
00170             const RangeDefinition & range, bool use_area_4vector,
00171             double & median, double & sigma, double & mean_area,
00172             bool all_are_incl) const {
00173 
00174   _check_jet_alg_good_for_median();
00175 
00176   vector<double> pt_over_areas;
00177   double total_area  = 0.0;
00178   double total_njets = 0;
00179 
00180   for (unsigned i = 0; i < all_jets.size(); i++) {
00181     if (range.is_in_range(all_jets[i])) {
00182       double this_area;
00183       if (use_area_4vector) {
00184           this_area = area_4vector(all_jets[i]).perp();
00185       } else {
00186           this_area = area(all_jets[i]);
00187       }
00188 
00189       if (this_area>0) {
00190         pt_over_areas.push_back(all_jets[i].perp()/this_area);
00191       } else {
00192         _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).");
00193       }
00194 
00195       total_area  += this_area;
00196       total_njets += 1.0;
00197     }
00198   }
00199 
00200   // there is nothing inside our region, so answer will always be zero
00201   if (pt_over_areas.size() == 0) {
00202     median = 0.0;
00203     sigma  = 0.0;
00204     mean_area = 0.0;
00205     return;
00206   }
00207   
00208   // get median (pt/area) [this is the "old" median definition. It considers
00209   // only the "real" jets in calculating the median, i.e. excluding the
00210   // only-ghost ones; it will be supplemented with more info below]
00211   sort(pt_over_areas.begin(), pt_over_areas.end());
00212 
00213   // now get the median & error, accounting for empty jets
00214   // define the fractions of distribution at median, median-1sigma
00215   double posn[2] = {0.5, (1.0-0.6827)/2.0};
00216   double res[2];
00217   
00218   double n_empty, empty_a;
00219   if (has_explicit_ghosts()) {
00220     // NB: the following lines of code are potentially incorrect in cases
00221     //     where there are unclustered particles (empty_area would do a better job,
00222     //     at least for active areas). This is not an issue with kt or C/A, or other
00223     //     algorithms that cluster all particles (and the median estimation should in 
00224     //     any case only be done with kt or C/A!)
00225     empty_a = 0.0;
00226     n_empty = 0;
00227   } else if (all_are_incl) {
00228     // the default case
00229     empty_a = empty_area(range);
00230     n_empty = n_empty_jets(range);
00231   } else {
00232     // this one is intended to be used when e.g. one runs C/A, then looks at its
00233     // exclusive jets in order to get an effective smaller R value, and passes those
00234     // to this routine.
00235     empty_a = empty_area_from_jets(all_jets, range);
00236     mean_area = total_area / total_njets; // temporary value
00237     n_empty   = empty_a / mean_area;
00238   }
00239   //cout << "*** tot_area = " << total_area << ", empty_a = " << empty_a << endl;
00240   //cout << "*** n_empty = " << n_empty << ", ntotal =  " << total_njets << endl;
00241   total_njets += n_empty;
00242   total_area  += empty_a;
00243 
00244   // we need an int (rather than an unsigned int) with the size of the
00245   // pt_over_areas array, because we'll often be doing subtraction of
00246   // -1, negating it, etc. All of these operations go crazy with unsigned ints.
00247   int pt_over_areas_size = pt_over_areas.size();
00248   if (n_empty < -pt_over_areas_size/4)
00249     _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.");
00250 
00251   for (int i = 0; i < 2; i++) {
00252     double nj_median_pos = 
00253       (pt_over_areas_size-1.0 + n_empty)*posn[i] - n_empty;
00254     double nj_median_ratio;
00255     if (nj_median_pos >= 0 && pt_over_areas_size > 1) {
00256       int int_nj_median = int(nj_median_pos);
00257 
00258       // avoid potential overflow issues
00259       if (int_nj_median+1 > pt_over_areas_size-1){
00260         int_nj_median = pt_over_areas_size-2;
00261         nj_median_pos = pt_over_areas_size-1;
00262       }
00263 
00264       nj_median_ratio = 
00265         pt_over_areas[int_nj_median] * (int_nj_median+1-nj_median_pos)
00266         + pt_over_areas[int_nj_median+1] * (nj_median_pos - int_nj_median);
00267     } else {
00268       nj_median_ratio = 0.0;
00269     }
00270     res[i] = nj_median_ratio;
00271   }
00272   median = res[0];
00273   double error  = res[0] - res[1];
00274   mean_area = total_area / total_njets;
00275   sigma  = error * sqrt(mean_area);
00276 }
00277 
00278 
00283 vector<PseudoJet> ClusterSequenceAreaBase::subtracted_jets(const double rho,
00284                                                            const double ptmin) 
00285                                                            const {
00286   vector<PseudoJet> sub_jets;
00287   vector<PseudoJet> jets = sorted_by_pt(inclusive_jets(ptmin));
00288   for (unsigned i=0; i<jets.size(); i++) {
00289      PseudoJet sub_jet = subtracted_jet(jets[i],rho);
00290      sub_jets.push_back(sub_jet);
00291   }
00292   return sub_jets;
00293 }
00294 
00299 vector<PseudoJet> ClusterSequenceAreaBase::subtracted_jets(
00300                                                  const RangeDefinition & range, 
00301                                                  const double ptmin)
00302                                                  const {
00303   double rho = median_pt_per_unit_area_4vector(range);
00304   return subtracted_jets(rho,ptmin);
00305 }
00306 
00307 
00309 PseudoJet ClusterSequenceAreaBase::subtracted_jet(const PseudoJet & jet,
00310                                                   const double rho) const {
00311   PseudoJet area4vect = area_4vector(jet);
00312   PseudoJet sub_jet;
00313   // sanity check
00314   if (rho*area4vect.perp() < jet.perp() ) { 
00315     sub_jet = jet - rho*area4vect;
00316   } else { sub_jet = PseudoJet(0.0,0.0,0.0,0.0); }
00317   
00318   // make sure the subtracted jet has the same index (cluster and user)
00319   // (i.e. "looks like") the original jet
00320   sub_jet.set_cluster_hist_index(jet.cluster_hist_index());
00321   sub_jet.set_user_index(jet.user_index());
00322   
00323   return sub_jet;
00324 }
00325 
00326 
00330 PseudoJet ClusterSequenceAreaBase::subtracted_jet(const PseudoJet & jet,
00331                                        const RangeDefinition & range) const {
00332   double rho = median_pt_per_unit_area_4vector(range);
00333   PseudoJet sub_jet = subtracted_jet(jet, rho);
00334   return sub_jet;
00335 }
00336 
00337 
00339 double ClusterSequenceAreaBase::subtracted_pt(const PseudoJet & jet,
00340                                               const double rho,
00341                                               bool use_area_4vector) const {
00342   if ( use_area_4vector ) { 
00343      PseudoJet sub_jet = subtracted_jet(jet,rho);
00344      return sub_jet.perp();
00345   } else {
00346      return jet.perp() - rho*area(jet);
00347   }
00348 }  
00349 
00350 
00354 double ClusterSequenceAreaBase::subtracted_pt(const PseudoJet & jet,
00355                                               const RangeDefinition & range,
00356                                               bool use_area_4vector) const {
00357   if ( use_area_4vector ) { 
00358      PseudoJet sub_jet = subtracted_jet(jet,range);
00359      return sub_jet.perp();
00360   } else {
00361      double rho = median_pt_per_unit_area(range);
00362      return subtracted_pt(jet,rho,false);
00363   }
00364 }  
00365 
00366 
00368 void ClusterSequenceAreaBase::_check_jet_alg_good_for_median() const {
00369   if (jet_def().jet_algorithm() != kt_algorithm
00370       && jet_def().jet_algorithm() != cambridge_algorithm
00371       && jet_def().jet_algorithm() !=  cambridge_for_passive_algorithm) {
00372     _warnings.warn("ClusterSequenceAreaBase: jet_def being used may not be suitable for estimating diffuse backgrounds (good options are kt, cam)");
00373   }
00374 }
00375 
00376 
00377 
00378 FASTJET_END_NAMESPACE
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines