FastJet 3.0.2
|
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