JETAN/fastjet/fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh

   1 //STARTHEADER
   2 // $Id: ClusterSequenceActiveAreaExplicitGhosts.hh 1303 2008-10-29 16:42:09Z salam $
   3 //
   4 // Copyright (c) 2005-2006, Matteo Cacciari and Gavin Salam
   5 //
   6 //----------------------------------------------------------------------
   7 // This file is part of FastJet.
   8 //
   9 //  FastJet is free software; you can redistribute it and/or modify
  10 //  it under the terms of the GNU General Public License as published by
  11 //  the Free Software Foundation; either version 2 of the License, or
  12 //  (at your option) any later version.
  13 //
  14 //  The algorithms that underlie FastJet have required considerable
  15 //  development and are described in hep-ph/0512210. If you use
  16 //  FastJet as part of work towards a scientific publication, please
  17 //  include a citation to the FastJet paper.
  18 //
  19 //  FastJet is distributed in the hope that it will be useful,
  20 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  21 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  22 //  GNU General Public License for more details.
  23 //
  24 //  You should have received a copy of the GNU General Public License
  25 //  along with FastJet; if not, write to the Free Software
  26 //  Foundation, Inc.:
  27 //      59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  28 //----------------------------------------------------------------------
  29 //ENDHEADER
  30
  31 #ifndef __FASTJET_CLUSTERSEQUENCEACTIVEAREAEXPLICITGHOSTS_HH_
  32 #define __FASTJET_CLUSTERSEQUENCEACTIVEAREAEXPLICITGHOSTS_HH_
  33
  34 #include "fastjet/PseudoJet.hh"
  35 #include "fastjet/ClusterSequenceAreaBase.hh"
  36 #include "fastjet/GhostedAreaSpec.hh"
  37 #include "fastjet/internal/LimitedWarning.hh"
  38 #include<iostream>
  39 #include<vector>
  40
  41 namespace fastjet {      // defined in fastjet/internal/base.hh
  42
  43 //======================================================================
  44 /// Class that behaves essentially like ClusterSequence except
  45 /// that it also provides access to the area of a jet (which
  46 /// will be a random quantity... Figure out what to do about seeds
  47 /// later...)
  48 class ClusterSequenceActiveAreaExplicitGhosts :
  49   public ClusterSequenceAreaBase {
  50 public:
  51   /// constructor using a GhostedAreaSpec to specify how the area is
  52   /// to be measured
  53   template<class L> ClusterSequenceActiveAreaExplicitGhosts
  54          (const std::vector<L> & pseudojets,
  55           const JetDefinition & jet_def,
  56           const GhostedAreaSpec & area_spec,
  57           const bool & writeout_combinations = false)
  58            : ClusterSequenceAreaBase() {
  59            std::vector<L> * ghosts = NULL;
  60            _initialise(pseudojets,jet_def,&area_spec,ghosts,0.0,
  61                        writeout_combinations); }
  62
  63   template<class L> ClusterSequenceActiveAreaExplicitGhosts
  64          (const std::vector<L> & pseudojets,
  65           const JetDefinition & jet_def,
  66           const std::vector<L> & ghosts,
  67           double ghost_area,
  68           const bool & writeout_combinations = false)
  69            : ClusterSequenceAreaBase() {
  70            const GhostedAreaSpec * area_spec = NULL;
  71            _initialise(pseudojets,jet_def,area_spec,&ghosts,ghost_area,
  72                        writeout_combinations); }
  73
  74
  75   /// does the actual work of initialisation
  76   template<class L> void _initialise
  77          (const std::vector<L> & pseudojets,
  78           const JetDefinition & jet_def,
  79           const GhostedAreaSpec * area_spec,
  80           const std::vector<L> * ghosts,
  81           double                 ghost_area,
  82           const bool & writeout_combinations);
  83
  84   //vector<PseudoJet> constituents (const PseudoJet & jet) const;
  85
  86   /// returns the number of hard particles (i.e. those supplied by the user).
  87   unsigned int n_hard_particles() const;
  88
  89   /// returns the area of a jet
  90   virtual double area (const PseudoJet & jet) const;
  91
  92   /// returns a four vector corresponding to the sum (E-scheme) of the
  93   /// ghost four-vectors composing the jet area, normalised such that
  94   /// for a small contiguous area the p_t of the extended_area jet is
  95   /// equal to area of the jet.
  96   virtual PseudoJet area_4vector (const PseudoJet & jet) const;
  97
  98   /// true if a jet is made exclusively of ghosts
  99   virtual bool is_pure_ghost(const PseudoJet & jet) const;
 100
 101   /// true if the entry in the history index corresponds to a
 102   /// ghost; if hist_ix does not correspond to an actual particle
 103   /// (i.e. hist_ix < 0), then the result is false.
 104   bool is_pure_ghost(int history_index) const;
 105
 106   /// this class does have explicit ghosts
 107   virtual bool has_explicit_ghosts() const {return true;}
 108
 109   /// return the total area, up to |y|<maxrap, that consists of
 110   /// unclustered ghosts
 111   virtual double empty_area(const RangeDefinition & range) const;
 112
 113   /// returns the total area under study
 114   double total_area () const;
 115
 116   /// returns the largest squared transverse momentum among
 117   /// all ghosts
 118   double max_ghost_perp2() const {return _max_ghost_perp2;}
 119
 120   /// returns true if there are any particles whose transverse momentum
 121   /// if so low that there's a risk of the ghosts having modified the
 122   /// clustering sequence
 123   bool has_dangerous_particles() const {return _has_dangerous_particles;}
 124
 125 private:
 126
 127   int    _n_ghosts;
 128   double _ghost_area;
 129   std::vector<bool> _is_pure_ghost;
 130   std::vector<double> _areas;
 131   std::vector<PseudoJet> _area_4vectors;
 132
 133   // things related to checks for dangerous particles
 134   double _max_ghost_perp2;
 135   bool   _has_dangerous_particles;
 136   static LimitedWarning _warnings;
 137
 138   //static int _n_warn_dangerous_particles;
 139   //static const int _max_warn_dangerous_particles = 5;
 140
 141
 142   unsigned int _initial_hard_n;
 143
 144   /// adds the "ghost" momenta, which will be used to estimate
 145   /// the jet area
 146   void _add_ghosts(const GhostedAreaSpec & area_spec);
 147
 148   /// another way of adding ghosts
 149   template<class L> void _add_ghosts (
 150           const std::vector<L> & ghosts,
 151           double                 ghost_area);
 152
 153   /// routine to be called after the processing is done so as to
 154   /// establish summary information on all the jets (areas, whether
 155   /// pure ghost, etc.)
 156   void _post_process();
 157
 158 };
 159
 160
 161 //----------------------------------------------------------------------
 162 // initialise from some generic type... Has to be made available
 163 // here in order for the template aspect of it to work...
 164 template<class L> void ClusterSequenceActiveAreaExplicitGhosts::_initialise
 165          (const std::vector<L> & pseudojets,
 166           const JetDefinition & jet_def,
 167           const GhostedAreaSpec * area_spec,
 168           const std::vector<L> * ghosts,
 169           double                 ghost_area,
 170           const bool & writeout_combinations) {
 171   // don't reserve space yet -- will be done below
 172
 173   // insert initial jets this way so that any type L that can be
 174   // converted to a pseudojet will work fine (basically PseudoJet
 175   // and any type that has [] subscript access to the momentum
 176   // components, such as CLHEP HepLorentzVector).
 177   for (unsigned int i = 0; i < pseudojets.size(); i++) {
 178     PseudoJet mom(pseudojets[i]);
 179     //mom.set_user_index(0); // for user's particles (user index now lost...)
 180     _jets.push_back(mom);
 181     _is_pure_ghost.push_back(false);
 182   }
 183
 184   _initial_hard_n = _jets.size();
 185
 186   if (area_spec != NULL) {
 187     _add_ghosts(*area_spec);
 188   } else {
 189     _add_ghosts(*ghosts, ghost_area);
 190   }
 191
 192   if (writeout_combinations) {
 193     std::cout << "# Printing particles including ghosts\n";
 194     for (unsigned j = 0; j < _jets.size(); j++) {
 195       printf("%5u %20.13f %20.13f %20.13e\n",
 196                j,_jets[j].rap(),_jets[j].phi_02pi(),_jets[j].kt2());
 197     }
 198     std::cout << "# Finished printing particles including ghosts\n";
 199   }
 200
 201   // this will ensure that we can still point to jets without
 202   // difficulties arising!
 203   _jets.reserve(_jets.size()*2);
 204
 205   // run the clustering
 206   _initialise_and_run(jet_def,writeout_combinations);
 207
 208   // set up all other information
 209   _post_process();
 210 }
 211
 212
 213 inline unsigned int ClusterSequenceActiveAreaExplicitGhosts::n_hard_particles() const {return _initial_hard_n;}
 214
 215
 216 //----------------------------------------------------------------------
 217 /// add an explicitly specified bunch of ghosts
 218 template<class L> void ClusterSequenceActiveAreaExplicitGhosts::_add_ghosts (
 219           const std::vector<L> & ghosts,
 220           double                 ghost_area) {
 221
 222
 223   for (unsigned i = 0; i < ghosts.size(); i++) {
 224     _is_pure_ghost.push_back(true);
 225     _jets.push_back(ghosts[i]);
 226   }
 227   // and record some info about ghosts
 228   _ghost_area = ghost_area;
 229   _n_ghosts   = ghosts.size();
 230 }
 231
 232
 233 } // fastjet namespace
 234
 235 #endif // __FASTJET_CLUSTERSEQUENCEACTIVEAREAEXPLICITGHOSTS_HH_