ClusterSequenceAreaBase.cc

Go to the documentation of this file.
00001 
00002 //STARTHEADER
00003 // $Id: ClusterSequenceAreaBase.cc 960 2007-11-29 17:03:47Z cacciari $
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 
00037 FASTJET_BEGIN_NAMESPACE
00038 
00039 using namespace std;
00040 
00041 
00043 LimitedWarning ClusterSequenceAreaBase::_warnings;
00044 
00045 //----------------------------------------------------------------------
00050 double ClusterSequenceAreaBase::empty_area(const RangeDefinition & range) const {
00051   double empty = range.area();
00052   vector<PseudoJet> incl_jets(inclusive_jets(0.0));
00053   for (unsigned i = 0; i < incl_jets.size(); i++) {
00054     if (range.is_in_range(incl_jets[i])) empty -= area(incl_jets[i]);
00055   }
00056   return empty;
00057 }
00058 
00059 double ClusterSequenceAreaBase::median_pt_per_unit_area(const RangeDefinition & range) const {
00060   return median_pt_per_unit_something(range,false);
00061 }
00062 
00063 double ClusterSequenceAreaBase::median_pt_per_unit_area_4vector(const RangeDefinition & range) const {
00064   return median_pt_per_unit_something(range,true);
00065 }
00066 
00067 
00068 //----------------------------------------------------------------------
00072 double ClusterSequenceAreaBase::median_pt_per_unit_something(
00073                 const RangeDefinition & range, bool use_area_4vector) const {
00074 
00075   double median, sigma, mean_area;
00076   get_median_rho_and_sigma(range, use_area_4vector, median, sigma, mean_area);
00077   return median;
00078 
00079 }
00080 
00081 
00082 //----------------------------------------------------------------------
00086 void ClusterSequenceAreaBase::parabolic_pt_per_unit_area(
00087        double & a, double & b, const RangeDefinition & range, 
00088        double exclude_above, bool use_area_4vector) const {
00089   
00090   int n=0;
00091   int n_excluded = 0;
00092   double mean_f=0, mean_x2=0, mean_x4=0, mean_fx2=0; 
00093 
00094   vector<PseudoJet> incl_jets = inclusive_jets();
00095 
00096   for (unsigned i = 0; i < incl_jets.size(); i++) {
00097     if (range.is_in_range(incl_jets[i])) {
00098       double this_area;
00099       if ( use_area_4vector ) {
00100           this_area = area_4vector(incl_jets[i]).perp();     
00101       } else {
00102           this_area = area(incl_jets[i]);
00103       }
00104       double f = incl_jets[i].perp()/this_area;
00105       if (exclude_above <= 0.0 || f < exclude_above) {
00106         double x = incl_jets[i].rap(); double x2 = x*x;
00107         mean_f   += f;
00108         mean_x2  += x2;
00109         mean_x4  += x2*x2;
00110         mean_fx2 += f*x2;
00111         n++;
00112       } else {
00113         n_excluded++;
00114       }
00115     }
00116   }
00117 
00118   if (n <= 1) {
00119     // meaningful results require at least two jets inside the
00120     // area -- mind you if there are empty jets we should be in 
00121     // any case doing something special...
00122     a = 0.0;
00123     b = 0.0;
00124   } else {
00125     mean_f   /= n;
00126     mean_x2  /= n;
00127     mean_x4  /= n;
00128     mean_fx2 /= n;
00129     
00130     b = (mean_f*mean_x2 - mean_fx2)/(mean_x2*mean_x2 - mean_x4);
00131     a = mean_f - b*mean_x2;
00132   }
00133   //cerr << "n_excluded = "<< n_excluded << endl;
00134 }
00135 
00136 
00137 
00138 
00139 void ClusterSequenceAreaBase::get_median_rho_and_sigma(
00140             const RangeDefinition & range, bool use_area_4vector,
00141             double & median, double & sigma, double & mean_area) const {
00142 
00143   _check_jet_alg_good_for_median();
00144 
00145   vector<double> pt_over_areas;
00146   vector<PseudoJet> incl_jets = inclusive_jets();
00147   double total_area  = 0.0;
00148   double total_njets = 0;
00149 
00150   for (unsigned i = 0; i < incl_jets.size(); i++) {
00151     if (range.is_in_range(incl_jets[i])) {
00152       double this_area;
00153       if (use_area_4vector) {
00154           this_area = area_4vector(incl_jets[i]).perp();
00155       } else {
00156           this_area = area(incl_jets[i]);
00157       }
00158       pt_over_areas.push_back(incl_jets[i].perp()/this_area);
00159       total_area  += this_area;
00160       total_njets += 1.0;
00161     }
00162   }
00163 
00164   // there is nothing inside our region, so answer will always be zero
00165   if (pt_over_areas.size() == 0) {
00166     median = 0.0;
00167     sigma  = 0.0;
00168     mean_area = 0.0;
00169     return;
00170   }
00171   
00172   // get median (pt/area) [this is the "old" median definition. It considers
00173   // only the "real" jets in calculating the median, i.e. excluding the
00174   // only-ghost ones]
00175   sort(pt_over_areas.begin(), pt_over_areas.end());
00176 
00177   // now get the median & error, accounting for empty jets
00178   // define the fractions of distribution at median, median-1sigma
00179   double posn[2] = {0.5, (1.0-0.6827)/2.0};
00180   double res[2];
00181   
00182   double n_empty = n_empty_jets(range);
00183   total_njets += n_empty;
00184   total_area  += empty_area(range);
00185 
00186   for (int i = 0; i < 2; i++) {
00187     double nj_median_pos = 
00188       (pt_over_areas.size()-1 + n_empty)*posn[i] - n_empty;
00189     double nj_median_ratio;
00190     if (nj_median_pos >= 0 && pt_over_areas.size() > 1) {
00191       int int_nj_median = int(nj_median_pos);
00192       nj_median_ratio = 
00193         pt_over_areas[int_nj_median] * (int_nj_median+1-nj_median_pos)
00194         + pt_over_areas[int_nj_median+1] * (nj_median_pos - int_nj_median);
00195     } else {
00196       nj_median_ratio = 0.0;
00197     }
00198     res[i] = nj_median_ratio;
00199   }
00200   median = res[0];
00201   double error  = res[0] - res[1];
00202   mean_area = total_area / total_njets;
00203   sigma  = error * sqrt(mean_area);
00204 }
00205 
00206 
00211 vector<PseudoJet> ClusterSequenceAreaBase::subtracted_jets(const double rho,
00212                                                            const double ptmin) 
00213                                                            const {
00214   vector<PseudoJet> sub_jets;
00215   vector<PseudoJet> jets = sorted_by_pt(inclusive_jets(ptmin));
00216   for (unsigned i=0; i<jets.size(); i++) {
00217      PseudoJet sub_jet = subtracted_jet(jets[i],rho);
00218      sub_jets.push_back(sub_jet);
00219   }
00220   return sub_jets;
00221 }
00222 
00227 vector<PseudoJet> ClusterSequenceAreaBase::subtracted_jets(
00228                                                  const RangeDefinition & range, 
00229                                                  const double ptmin)
00230                                                  const {
00231   double rho = median_pt_per_unit_area_4vector(range);
00232   return subtracted_jets(rho,ptmin);
00233 }
00234 
00235 
00237 PseudoJet ClusterSequenceAreaBase::subtracted_jet(const PseudoJet & jet,
00238                                                   const double rho) const {
00239   PseudoJet area4vect = area_4vector(jet);
00240   PseudoJet sub_jet;
00241   // sanity check
00242   if (rho*area4vect.perp() < jet.perp() ) { 
00243     sub_jet = jet - rho*area4vect;
00244   } else { sub_jet = PseudoJet(0.0,0.0,0.0,0.0); }
00245   return sub_jet;
00246 }
00247 
00248 
00252 PseudoJet ClusterSequenceAreaBase::subtracted_jet(const PseudoJet & jet,
00253                                        const RangeDefinition & range) const {
00254   double rho = median_pt_per_unit_area_4vector(range);
00255   PseudoJet sub_jet = subtracted_jet(jet, rho);
00256   return sub_jet;
00257 }
00258 
00259 
00261 double ClusterSequenceAreaBase::subtracted_pt(const PseudoJet & jet,
00262                                               const double rho,
00263                                               bool use_area_4vector) const {
00264   if ( use_area_4vector ) { 
00265      PseudoJet sub_jet = subtracted_jet(jet,rho);
00266      return sub_jet.perp();
00267   } else {
00268      return jet.perp() - rho*area(jet);
00269   }
00270 }  
00271 
00272 
00276 double ClusterSequenceAreaBase::subtracted_pt(const PseudoJet & jet,
00277                                               const RangeDefinition & range,
00278                                               bool use_area_4vector) const {
00279   if ( use_area_4vector ) { 
00280      PseudoJet sub_jet = subtracted_jet(jet,range);
00281      return sub_jet.perp();
00282   } else {
00283      double rho = median_pt_per_unit_area(range);
00284      return subtracted_pt(jet,rho,false);
00285   }
00286 }  
00287 
00288 
00290 void ClusterSequenceAreaBase::_check_jet_alg_good_for_median() const {
00291   if (jet_def().jet_algorithm() != kt_algorithm
00292       && jet_def().jet_algorithm() != cambridge_algorithm
00293       && jet_def().jet_algorithm() !=  cambridge_for_passive_algorithm) {
00294     _warnings.warn("ClusterSequenceAreaBase: jet_def being used may not be suitable for estimating diffuse backgrounds (good options are kt, cam)");
00295   }
00296 }
00297 
00298 
00299 
00300 FASTJET_END_NAMESPACE

Generated on Tue Dec 18 17:05:03 2007 for fastjet by  doxygen 1.5.2