FastJet 3.0.5
ClusterSequenceActiveAreaExplicitGhosts.hh
00001 //STARTHEADER
00002 // $Id: ClusterSequenceActiveAreaExplicitGhosts.hh 2687 2011-11-14 11:17:51Z 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 #ifndef __FASTJET_CLUSTERSEQUENCEACTIVEAREAEXPLICITGHOSTS_HH_
00030 #define __FASTJET_CLUSTERSEQUENCEACTIVEAREAEXPLICITGHOSTS_HH_ 
00031 
00032 #include "fastjet/PseudoJet.hh"
00033 #include "fastjet/ClusterSequenceAreaBase.hh"
00034 #include "fastjet/GhostedAreaSpec.hh"
00035 #include "fastjet/LimitedWarning.hh"
00036 #include<iostream>
00037 #include<vector>
00038 #include <cstdio>
00039 
00040 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
00041 
00042 //======================================================================
00043 /// @ingroup sec_area_classes
00044 /// \class ClusterSequenceActiveAreaExplicitGhosts
00045 /// Like ClusterSequence with computation of the active jet area with the
00046 /// addition of explicit ghosts
00047 ///
00048 /// Class that behaves essentially like ClusterSequence except
00049 /// that it also provides access to the area of a jet (which
00050 /// will be a random quantity... Figure out what to do about seeds 
00051 /// later...)
00052 ///
00053 /// This class should not be used directly. Rather use
00054 /// ClusterSequenceArea with the appropriate AreaDefinition
00055 class ClusterSequenceActiveAreaExplicitGhosts : 
00056   public ClusterSequenceAreaBase {
00057 public:
00058   /// constructor using a GhostedAreaSpec to specify how the area is
00059   /// to be measured
00060   template<class L> ClusterSequenceActiveAreaExplicitGhosts
00061          (const std::vector<L> & pseudojets, 
00062           const JetDefinition & jet_def_in,
00063           const GhostedAreaSpec & ghost_spec,
00064           const bool & writeout_combinations = false) 
00065            : ClusterSequenceAreaBase() {
00066            std::vector<L> * ghosts = NULL;
00067            _initialise(pseudojets,jet_def_in,&ghost_spec,ghosts,0.0,
00068                        writeout_combinations); }
00069 
00070   template<class L> ClusterSequenceActiveAreaExplicitGhosts
00071          (const std::vector<L> & pseudojets, 
00072           const JetDefinition & jet_def_in,
00073           const std::vector<L> & ghosts,
00074           double ghost_area,
00075           const bool & writeout_combinations = false) 
00076            : ClusterSequenceAreaBase() {
00077            const GhostedAreaSpec * ghost_spec = NULL;
00078            _initialise(pseudojets,jet_def_in,ghost_spec,&ghosts,ghost_area,
00079                        writeout_combinations); }
00080 
00081 
00082   /// does the actual work of initialisation
00083   template<class L> void _initialise
00084          (const std::vector<L> & pseudojets, 
00085           const JetDefinition & jet_def_in,
00086           const GhostedAreaSpec * ghost_spec,
00087           const std::vector<L> * ghosts,
00088           double                 ghost_area,
00089           const bool & writeout_combinations); 
00090 
00091   //vector<PseudoJet> constituents (const PseudoJet & jet) const;
00092 
00093   /// returns the number of hard particles (i.e. those supplied by the user).
00094   unsigned int n_hard_particles() const;
00095 
00096   /// returns the area of a jet
00097   virtual double area (const PseudoJet & jet) const;
00098 
00099   /// returns a four vector corresponding to the sum (E-scheme) of the
00100   /// ghost four-vectors composing the jet area, normalised such that
00101   /// for a small contiguous area the p_t of the extended_area jet is
00102   /// equal to area of the jet.
00103   virtual PseudoJet area_4vector (const PseudoJet & jet) const;
00104 
00105   /// true if a jet is made exclusively of ghosts
00106   virtual bool is_pure_ghost(const PseudoJet & jet) const;
00107 
00108   /// true if the entry in the history index corresponds to a
00109   /// ghost; if hist_ix does not correspond to an actual particle
00110   /// (i.e. hist_ix < 0), then the result is false.
00111   bool is_pure_ghost(int history_index) const;
00112 
00113   /// this class does have explicit ghosts
00114   virtual bool has_explicit_ghosts() const {return true;}
00115 
00116   /// return the total area, corresponding to a given Selector, that
00117   /// consists of unclustered ghosts
00118   ///
00119   /// The selector needs to apply jet by jet
00120   virtual double empty_area(const Selector & selector) const;
00121 
00122   /// returns the total area under study
00123   double total_area () const;
00124   
00125   /// returns the largest squared transverse momentum among
00126   /// all ghosts
00127   double max_ghost_perp2() const {return _max_ghost_perp2;}
00128 
00129   /// returns true if there are any particles whose transverse momentum
00130   /// if so low that there's a risk of the ghosts having modified the
00131   /// clustering sequence
00132   bool has_dangerous_particles() const {return _has_dangerous_particles;}
00133 
00134   /// get the area of the ghosts
00135   //double ghost_area() const{return _ghost_area;}
00136 
00137 private:
00138 
00139   int    _n_ghosts;
00140   double _ghost_area;
00141   std::vector<bool> _is_pure_ghost;
00142   std::vector<double> _areas;
00143   std::vector<PseudoJet> _area_4vectors;
00144   
00145   // things related to checks for dangerous particles
00146   double _max_ghost_perp2;
00147   bool   _has_dangerous_particles; 
00148   static LimitedWarning _warnings;
00149 
00150   //static int _n_warn_dangerous_particles;
00151   //static const int _max_warn_dangerous_particles = 5;
00152 
00153   
00154   unsigned int _initial_hard_n;
00155 
00156   /// adds the "ghost" momenta, which will be used to estimate
00157   /// the jet area
00158   void _add_ghosts(const GhostedAreaSpec & ghost_spec); 
00159 
00160   /// another way of adding ghosts
00161   template<class L> void _add_ghosts (
00162           const std::vector<L> & ghosts,
00163           double                 ghost_area);
00164 
00165   /// routine to be called after the processing is done so as to
00166   /// establish summary information on all the jets (areas, whether
00167   /// pure ghost, etc.)
00168   void _post_process();
00169 
00170 };
00171 
00172 
00173 //----------------------------------------------------------------------
00174 // initialise from some generic type... Has to be made available
00175 // here in order for the template aspect of it to work...
00176 template<class L> void ClusterSequenceActiveAreaExplicitGhosts::_initialise
00177          (const std::vector<L> & pseudojets, 
00178           const JetDefinition & jet_def_in,
00179           const GhostedAreaSpec * ghost_spec,
00180           const std::vector<L> * ghosts,
00181           double                 ghost_area,
00182           const bool & writeout_combinations) {
00183   // don't reserve space yet -- will be done below
00184 
00185   // insert initial jets this way so that any type L that can be
00186   // converted to a pseudojet will work fine (basically PseudoJet
00187   // and any type that has [] subscript access to the momentum
00188   // components, such as CLHEP HepLorentzVector).
00189   for (unsigned int i = 0; i < pseudojets.size(); i++) {
00190     PseudoJet mom(pseudojets[i]);
00191     //mom.set_user_index(0); // for user's particles (user index now lost...)
00192     _jets.push_back(mom);
00193     _is_pure_ghost.push_back(false);
00194   }
00195 
00196   _initial_hard_n = _jets.size();
00197 
00198   if (ghost_spec != NULL) {
00199     //std::cout << "about to reserve " << (_jets.size()+ghost_spec->n_ghosts())*2 << std::endl;
00200     _jets.reserve((_jets.size()+ghost_spec->n_ghosts()));
00201     _add_ghosts(*ghost_spec);
00202   } else {
00203     _jets.reserve(_jets.size()+ghosts->size());
00204     _add_ghosts(*ghosts, ghost_area);
00205   }
00206 
00207   if (writeout_combinations) {
00208     std::cout << "# Printing particles including ghosts\n";
00209     for (unsigned j = 0; j < _jets.size(); j++) {
00210       printf("%5u %20.13f %20.13f %20.13e\n",
00211                j,_jets[j].rap(),_jets[j].phi_02pi(),_jets[j].kt2());
00212     }
00213     std::cout << "# Finished printing particles including ghosts\n";
00214   }
00215 
00216   // this will ensure that we can still point to jets without
00217   // difficulties arising!
00218   //std::cout << _jets.size() << " " << _jets.size()*2 << " " << _jets.max_size() << std::endl;
00219   _jets.reserve(_jets.size()*2); //GPS tmp removed
00220 
00221   // run the clustering
00222   _initialise_and_run(jet_def_in,writeout_combinations);
00223 
00224   // set up all other information
00225   _post_process();
00226 }
00227 
00228 
00229 inline unsigned int ClusterSequenceActiveAreaExplicitGhosts::n_hard_particles() const {return _initial_hard_n;}
00230 
00231 
00232 //----------------------------------------------------------------------
00233 /// add an explicitly specified bunch of ghosts
00234 template<class L> void ClusterSequenceActiveAreaExplicitGhosts::_add_ghosts (
00235           const std::vector<L> & ghosts,
00236           double                 ghost_area) {
00237 
00238   
00239   for (unsigned i = 0; i < ghosts.size(); i++) {
00240     _is_pure_ghost.push_back(true);
00241     _jets.push_back(ghosts[i]);
00242   }
00243   // and record some info about ghosts
00244   _ghost_area = ghost_area;
00245   _n_ghosts   = ghosts.size();
00246 }
00247 
00248 
00249 FASTJET_END_NAMESPACE
00250 
00251 #endif // __FASTJET_CLUSTERSEQUENCEACTIVEAREAEXPLICITGHOSTS_HH_
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends