FastJet 3.0.2
JetMedianBackgroundEstimator.cc
00001 //STARTHEADER
00002 // $Id: JetMedianBackgroundEstimator.cc 2689 2011-11-14 14:51:06Z soyez $
00003 //
00004 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
00005 //
00006 //----------------------------------------------------------------------
00007 // This file is part of FastJet.
00008 //
00009 //  FastJet is free software; you can redistribute it and/or modify
00010 //  it under the terms of the GNU General Public License as published by
00011 //  the Free Software Foundation; either version 2 of the License, or
00012 //  (at your option) any later version.
00013 //
00014 //  The algorithms that underlie FastJet have required considerable
00015 //  development and are described in hep-ph/0512210. If you use
00016 //  FastJet as part of work towards a scientific publication, please
00017 //  include a citation to the FastJet paper.
00018 //
00019 //  FastJet is distributed in the hope that it will be useful,
00020 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00021 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00022 //  GNU General Public License for more details.
00023 //
00024 //  You should have received a copy of the GNU General Public License
00025 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
00026 //----------------------------------------------------------------------
00027 //ENDHEADER
00028 
00029 #include "fastjet/tools/JetMedianBackgroundEstimator.hh"
00030 #include <fastjet/ClusterSequenceArea.hh>
00031 #include <fastjet/ClusterSequenceStructure.hh>
00032 #include <iostream>
00033 
00034 FASTJET_BEGIN_NAMESPACE     // defined in fastjet/internal/base.hh
00035 
00036 using namespace std;
00037 
00038 double BackgroundJetScalarPtDensity::result(const PseudoJet & jet) const {
00039   std::vector<PseudoJet> constituents = jet.constituents();
00040   double scalar_pt = 0;
00041   for (unsigned i = 0; i < constituents.size(); i++) {
00042     scalar_pt += pow(constituents[i].perp(), _pt_power);
00043   }
00044   return scalar_pt / jet.area();
00045 }
00046 
00047 
00048 //----------------------------------------------------------------------
00049 double BackgroundRescalingYPolynomial::result(const PseudoJet & jet) const {
00050   double y = jet.rap();
00051   double y2 = y*y;
00052   double rescaling = _a0 + _a1*y + _a2*y2 + _a3*y2*y + _a4*y2*y2;
00053   return rescaling;
00054 }
00055 
00056 /// allow for warnings
00057 LimitedWarning JetMedianBackgroundEstimator::_warnings;
00058 LimitedWarning JetMedianBackgroundEstimator::_warnings_zero_area;
00059 LimitedWarning JetMedianBackgroundEstimator::_warnings_preliminary;
00060 
00061 
00062 //---------------------------------------------------------------------
00063 // class JetMedianBackgroundEstimator
00064 // Class to estimate the density of the background per unit area
00065 //---------------------------------------------------------------------
00066 
00067 //----------------------------------------------------------------------
00068 // ctors and dtors
00069 //----------------------------------------------------------------------
00070 // ctor that allows to set only the particles later on
00071 JetMedianBackgroundEstimator::JetMedianBackgroundEstimator(const Selector &rho_range,
00072                                          const JetDefinition &jet_def,
00073                                          const AreaDefinition &area_def)
00074   : _rho_range(rho_range), _jet_def(jet_def), _area_def(area_def) {
00075 
00076   // initialise things decently
00077   reset();
00078 
00079   // make a few checks
00080   _check_jet_alg_good_for_median();
00081 }
00082 
00083 
00084 //----------------------------------------------------------------------
00085 // ctor from a cluster sequence
00086 //  - csa        the ClusterSequenceArea to use
00087 //  - rho_range  the Selector specifying which jets will be considered
00088 JetMedianBackgroundEstimator::JetMedianBackgroundEstimator( const Selector &rho_range, const ClusterSequenceAreaBase &csa)
00089   : _rho_range(rho_range), _jet_def(JetDefinition()){
00090 
00091   // initialise things properly
00092   reset();
00093 
00094   // tell the BGE about the cluster sequence
00095   set_cluster_sequence(csa);
00096 }
00097 
00098 
00099 //----------------------------------------------------------------------
00100 // setting a new event
00101 //----------------------------------------------------------------------
00102 // tell the background estimator that it has a new event, composed
00103 // of the specified particles.
00104 void JetMedianBackgroundEstimator::set_particles(const vector<PseudoJet> & particles) {
00105   // make sure that we have been provided a genuine jet definition 
00106   if (_jet_def.jet_algorithm() == undefined_jet_algorithm)
00107     throw Error("JetMedianBackgroundEstimator::set_particles can only be called if you set the jet (and area) definition explicitly through the class constructor");
00108 
00109   // initialise things decently (including setting uptodate to false!)
00110   //reset();
00111   _uptodate=false;
00112 
00113   // cluster the particles
00114   // 
00115   // One may argue that it is better to cache the particles and only
00116   // do the clustering later but clustering the particles now has 2
00117   // practical advantages:
00118   //  - it allows to une only '_included_jets' in all that follows
00119   //  - it avoids adding another flag to ensure particles are 
00120   //    clustered only once
00121   ClusterSequenceArea *csa = new ClusterSequenceArea(particles, _jet_def, _area_def);
00122   _included_jets = csa->inclusive_jets();
00123 
00124   // store the CS for later on
00125   _csi = csa->structure_shared_ptr();
00126   csa->delete_self_when_unused();
00127 }
00128 
00129 //----------------------------------------------------------------------
00130 // (re)set the cluster sequence (with area support) to be used by
00131 // future calls to rho() etc. 
00132 //
00133 // \param csa  the cluster sequence area
00134 //
00135 // Pre-conditions: 
00136 //  - one should be able to estimate the "empty area" (i.e. the area
00137 //    not occupied by jets). This is feasible if at least one of the following
00138 //    conditions is satisfied:
00139 //     ( i) the ClusterSequence has explicit ghosts
00140 //     (ii) the range selected has a computable area.
00141 //  - the jet algorithm must be suited for median computation
00142 //    (otherwise a warning will be issues)
00143 //
00144 // Note that selectors with e.g. hardest-jets exclusion do not have
00145 // a well-defined area. For this reasons, it is STRONGLY advised to
00146 // use an area with explicit ghosts.
00147 void JetMedianBackgroundEstimator::set_cluster_sequence(const ClusterSequenceAreaBase & csa) {
00148   _csi = csa.structure_shared_ptr();
00149 
00150   // sanity checks
00151   //---------------
00152   //  (i) check the alg is appropriate
00153   _check_jet_alg_good_for_median();
00154 
00155   //  (ii) check that, if there are no explicit ghosts, the selector has a finite area
00156   if ((!csa.has_explicit_ghosts()) && (!_rho_range.has_finite_area())){
00157     throw Error("JetMedianBackgroundEstimator: either an area with explicit ghosts (recommended) or a Selector with finite area is needed (to allow for the computation of the empty area)");
00158   }
00159 
00160   // get the initial list of jets
00161   _included_jets = csa.inclusive_jets();
00162 
00163   _uptodate = false;
00164 }
00165 
00166 
00167 //----------------------------------------------------------------------
00168 // (re)set the jets (which must have area support) to be used by future
00169 // calls to rho() etc.; for the conditions that must be satisfied
00170 // by the jets, see the Constructor that takes jets.
00171 void JetMedianBackgroundEstimator::set_jets(const vector<PseudoJet> &jets) {
00172   
00173   if (! jets.size())
00174     throw Error("JetMedianBackgroundEstimator::JetMedianBackgroundEstimator: At least one jet is needed to compute the background properties");
00175 
00176   // sanity checks
00177   //---------------
00178   //  (o) check that there is an underlying CS shared by all the jets
00179   if (! (jets[0].has_associated_cluster_sequence()) && (jets[0].has_area()))
00180     throw Error("JetMedianBackgroundEstimator::JetMedianBackgroundEstimator: the jets used to estimate the background properties must be associated with a valid ClusterSequenceAreaBase");
00181 
00182   _csi = jets[0].structure_shared_ptr();
00183   ClusterSequenceStructure * csi = dynamic_cast<ClusterSequenceStructure*>(_csi());
00184   const ClusterSequenceAreaBase * csab = csi->validated_csab();
00185 
00186   for (unsigned int i=1;i<jets.size(); i++){
00187     if (! jets[i].has_associated_cluster_sequence()) // area automatic if the next test succeeds
00188       throw Error("JetMedianBackgroundEstimator::set_jets(...): the jets used to estimate the background properties must be associated with a valid ClusterSequenceAreaBase");
00189 
00190     if (jets[i].structure_shared_ptr().get() != _csi.get())
00191       throw Error("JetMedianBackgroundEstimator::set_jets(...): all the jets used to estimate the background properties must share the same ClusterSequence");
00192   }
00193 
00194   //  (i) check the alg is appropriate
00195   _check_jet_alg_good_for_median();
00196 
00197   //  (ii) check that, if there are no explicit ghosts, the selector has a finite area
00198   if ((!csab->has_explicit_ghosts()) && (!_rho_range.has_finite_area())){
00199     throw Error("JetMedianBackgroundEstimator: either an area with explicit ghosts (recommended) or a Selector with finite area is needed (to allow for the computation of the empty area)");
00200   }
00201 
00202 
00203   // get the initial list of jets
00204   _included_jets = jets;
00205 
00206   // ensure recalculation of quantities that need it
00207   _uptodate = false;
00208 }
00209 
00210 
00211 //----------------------------------------------------------------------
00212 // retrieving fundamental information
00213 //----------------------------------------------------------------
00214 
00215 // get rho, the median background density per unit area
00216 double JetMedianBackgroundEstimator::rho() const {
00217   if (_rho_range.takes_reference())
00218     throw Error("The background estimation is obtained from a selector that takes a reference jet. rho(PseudoJet) should be used in that case");
00219   _recompute_if_needed();
00220   return _rho;
00221 }
00222 
00223 // get sigma, the background fluctuations per unit area
00224 double JetMedianBackgroundEstimator::sigma() const {
00225   if (_rho_range.takes_reference())
00226     throw Error("The background estimation is obtained from a selector that takes a reference jet. rho(PseudoJet) should be used in that case");
00227   _recompute_if_needed();
00228   return _sigma;
00229 }
00230 
00231 // get rho, the median background density per unit area, locally at
00232 // the position of a given jet.
00233 //
00234 // If the Selector associated with the range takes a reference jet
00235 // (i.e. is relocatable), then for subsequent operations the
00236 // Selector has that jet set as its reference.
00237 double JetMedianBackgroundEstimator::rho(const PseudoJet & jet) {
00238   _recompute_if_needed(jet);
00239   double our_rho = _rho;
00240   if (_rescaling_class != 0) { 
00241     our_rho *= (*_rescaling_class)(jet);
00242   }
00243   return our_rho;
00244 }
00245 
00246 // get sigma, the background fluctuations per unit area,
00247 // locally at the position of a given jet.
00248 //
00249 // If the Selector associated with the range takes a reference jet
00250 // (i.e. is relocatable), then for subsequent operations the
00251 // Selector has that jet set as its reference.
00252 double JetMedianBackgroundEstimator::sigma(const PseudoJet &jet) {
00253   _recompute_if_needed(jet);
00254   double our_sigma = _sigma;
00255   if (_rescaling_class != 0) { 
00256     our_sigma *= (*_rescaling_class)(jet);
00257   }
00258   return our_sigma;
00259 }
00260 
00261 
00262 //----------------------------------------------------------------------
00263 // configuring behaviour
00264 //----------------------------------------------------------------------
00265 // reset to default values
00266 // 
00267 // set the variou options to their default values
00268 void JetMedianBackgroundEstimator::reset(){
00269   // set the remaining default parameters
00270   set_use_area_4vector();  // true by default
00271   set_provide_fj2_sigma(false);
00272 
00273   // reset the computed values
00274   _rho = _sigma = 0.0;
00275   _n_jets_used = _n_empty_jets = 0;
00276   _empty_area = _mean_area = 0.0;
00277 
00278   _included_jets.clear();
00279 
00280   _jet_density_class = 0; // null pointer
00281   _rescaling_class = 0;   // null pointer
00282 
00283   _uptodate = false;
00284 }
00285 
00286 
00287 // Set a pointer to a class that calculates the quantity whose
00288 // median will be calculated; if the pointer is null then pt/area
00289 // is used (as occurs also if this function is not called).
00290 void JetMedianBackgroundEstimator::set_jet_density_class(const FunctionOfPseudoJet<double> * jet_density_class_in) {
00291   _warnings_preliminary.warn("JetMedianBackgroundEstimator::set_jet_density_class: density classes are still preliminary in FastJet 3.0. Their interface may differ in future releases (without guaranteeing backward compatibility).");
00292   _jet_density_class = jet_density_class_in;
00293   _uptodate = false;
00294 }
00295 
00296 
00297 
00298 //----------------------------------------------------------------------
00299 // description
00300 //----------------------------------------------------------------------
00301 string JetMedianBackgroundEstimator::description() const { 
00302   ostringstream desc;
00303   desc << "JetMedianBackgroundEstimator, using " << _jet_def.description() 
00304        << " with " << _area_def.description() 
00305        << " and selecting jets with " << _rho_range.description();
00306   return desc.str();
00307 }       
00308 
00309 
00310 
00311 //----------------------------------------------------------------------
00312 // computation of the background properties
00313 //----------------------------------------------------------------------
00314 // for estimation using a relocatable selector (i.e. local range)
00315 // this allows to set its position. Note that this HAS to be called
00316 // before any attempt to compute the background properties
00317 void JetMedianBackgroundEstimator::_recompute_if_needed(const PseudoJet &jet){
00318   // if the range is relocatable, handles its relocation
00319   if (_rho_range.takes_reference()){
00320     // check that the reference is not the same as the previous one
00321     // (would avoid an unnecessary recomputation)
00322     if (jet == _current_reference) return;
00323 
00324     // relocate the range and make sure things get recomputed the next
00325     // time one tries to get some information
00326     _rho_range.set_reference(jet);
00327     _uptodate=false;
00328   }
00329 
00330   _recompute_if_needed();
00331 }
00332 
00333 // do the actual job
00334 void JetMedianBackgroundEstimator::_compute() const {
00335   // check if the clustersequence is still valid
00336   _check_csa_alive();
00337 
00338   // fill the vector of pt/area (or the quantity from the jet density class) 
00339   //  - in the range
00340   vector<double> vector_for_median;
00341   double total_area  = 0.0;
00342   _n_jets_used = 0;
00343 
00344   // apply the selector to the included jets
00345   vector<PseudoJet> selected_jets = _rho_range(_included_jets);
00346 
00347   // compute the pt/area for the selected jets
00348   for (unsigned i = 0; i < selected_jets.size(); i++) {
00349     const PseudoJet & current_jet = selected_jets[i];
00350 
00351     double this_area = (_use_area_4vector) ? current_jet.area_4vector().perp() : current_jet.area(); 
00352 
00353     if (this_area>0){
00354       double median_input;
00355       if (_jet_density_class == 0) {
00356         median_input = current_jet.perp()/this_area;
00357       } else {
00358         median_input = (*_jet_density_class)(current_jet);
00359       }
00360       if (_rescaling_class != 0) {
00361         median_input /= (*_rescaling_class)(current_jet);
00362       }
00363       vector_for_median.push_back(median_input);
00364       total_area  += this_area;
00365       _n_jets_used++;
00366     } else {
00367       _warnings_zero_area.warn("JetMedianBackgroundEstimator::_compute(...): 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).");
00368     }
00369       
00370   }
00371   
00372   // there is nothing inside our region, so answer will always be zero
00373   if (vector_for_median.size() == 0) {
00374     _rho        = 0.0;
00375     _sigma      = 0.0;
00376     _mean_area  = 0.0;
00377     return;
00378   }
00379 
00380   // determine the number of empty jets
00381   const ClusterSequenceAreaBase * csab = (dynamic_cast<ClusterSequenceStructure*>(_csi()))->validated_csab();
00382   if (csab->has_explicit_ghosts()) {
00383     _empty_area = 0.0;
00384     _n_empty_jets = 0;
00385   } else {
00386     _empty_area = csab->empty_area(_rho_range);
00387     _n_empty_jets = csab->n_empty_jets(_rho_range);
00388   }
00389 
00390   double total_njets = _n_jets_used + _n_empty_jets;
00391   total_area  += _empty_area;
00392 
00393   double stand_dev;
00394   _median_and_stddev(vector_for_median, _n_empty_jets, _rho, stand_dev, 
00395                      _provide_fj2_sigma);
00396 
00397   // process and store the results (_rho was already stored above)
00398   _mean_area  = total_area / total_njets;
00399   _sigma      = stand_dev * sqrt(_mean_area);
00400 
00401   // record that the computation has been performed  
00402   _uptodate = true;
00403 }
00404 
00405 
00406 
00407 // check that the underlying structure is still alive;
00408 // throw an error otherwise
00409 void JetMedianBackgroundEstimator::_check_csa_alive() const{
00410   ClusterSequenceStructure* csa = dynamic_cast<ClusterSequenceStructure*>(_csi());
00411   if (csa == 0) {
00412     throw Error("JetMedianBackgroundEstimator: there is no cluster sequence associated with the JetMedianBackgroundEstimator");
00413   }
00414   if (! dynamic_cast<ClusterSequenceStructure*>(_csi())->has_associated_cluster_sequence())
00415     throw Error("JetMedianBackgroundEstimator: modifications are no longer possible as the underlying ClusterSequence has gone out of scope");
00416 }
00417 
00418 
00419 // check that the algorithm used for the clustering is suitable for
00420 // background estimation (i.e. either kt or C/A).
00421 // Issue a warning otherwise
00422 void JetMedianBackgroundEstimator::_check_jet_alg_good_for_median() const{
00423   const JetDefinition * jet_def = &_jet_def;
00424 
00425   // if no explicit jet def has been provided, fall back on the
00426   // cluster sequence
00427   if (_jet_def.jet_algorithm() == undefined_jet_algorithm){
00428     const ClusterSequence * cs = dynamic_cast<ClusterSequenceStructure*>(_csi())->validated_cs();
00429     jet_def = &(cs->jet_def());
00430   }
00431 
00432   if (jet_def->jet_algorithm() != kt_algorithm
00433       && jet_def->jet_algorithm() != cambridge_algorithm
00434       && jet_def->jet_algorithm() != cambridge_for_passive_algorithm) {
00435     _warnings.warn("JetMedianBackgroundEstimator: jet_def being used may not be suitable for estimating diffuse backgrounds (good alternatives are kt, cam)");
00436   }
00437 }
00438 
00439 
00440 
00441 
00442 FASTJET_END_NAMESPACE
00443 
00444 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends