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 #include <stdio.h>
  41
  42 namespace fastjet {      // defined in fastjet/internal/base.hh
  43
  44 //======================================================================
  45 /// Class that behaves essentially like ClusterSequence except
  46 /// that it also provides access to the area of a jet (which
  47 /// will be a random quantity... Figure out what to do about seeds
  48 /// later...)
  49 class ClusterSequenceActiveAreaExplicitGhosts :
  50   public ClusterSequenceAreaBase {
  51 public:
  52   /// constructor using a GhostedAreaSpec to specify how the area is
  53   /// to be measured
  54   template<class L> ClusterSequenceActiveAreaExplicitGhosts
  55          (const std::vector<L> & pseudojets,
  56           const JetDefinition & jet_def,
  57           const GhostedAreaSpec & area_spec,
  58           const bool & writeout_combinations = false)
  59            : ClusterSequenceAreaBase() {
  60            std::vector<L> * ghosts = NULL;
  61            _initialise(pseudojets,jet_def,&area_spec,ghosts,0.0,
  62                        writeout_combinations); }
  63
  64   template<class L> ClusterSequenceActiveAreaExplicitGhosts
  65          (const std::vector<L> & pseudojets,
  66           const JetDefinition & jet_def,
  67           const std::vector<L> & ghosts,
  68           double ghost_area,
  69           const bool & writeout_combinations = false)
  70            : ClusterSequenceAreaBase() {
  71            const GhostedAreaSpec * area_spec = NULL;
  72            _initialise(pseudojets,jet_def,area_spec,&ghosts,ghost_area,
  73                        writeout_combinations); }
  74
  75
  76   /// does the actual work of initialisation
  77   template<class L> void _initialise
  78          (const std::vector<L> & pseudojets,
  79           const JetDefinition & jet_def,
  80           const GhostedAreaSpec * area_spec,
  81           const std::vector<L> * ghosts,
  82           double                 ghost_area,
  83           const bool & writeout_combinations);
  84
  85   //vector<PseudoJet> constituents (const PseudoJet & jet) const;
  86
  87   /// returns the number of hard particles (i.e. those supplied by the user).
  88   unsigned int n_hard_particles() const;
  89
  90   /// returns the area of a jet
  91   virtual double area (const PseudoJet & jet) const;
  92
  93   /// returns a four vector corresponding to the sum (E-scheme) of the
  94   /// ghost four-vectors composing the jet area, normalised such that
  95   /// for a small contiguous area the p_t of the extended_area jet is
  96   /// equal to area of the jet.
  97   virtual PseudoJet area_4vector (const PseudoJet & jet) const;
  98
  99   /// true if a jet is made exclusively of ghosts
 100   virtual bool is_pure_ghost(const PseudoJet & jet) const;
 101
 102   /// true if the entry in the history index corresponds to a
 103   /// ghost; if hist_ix does not correspond to an actual particle
 104   /// (i.e. hist_ix < 0), then the result is false.
 105   bool is_pure_ghost(int history_index) const;
 106
 107   /// this class does have explicit ghosts
 108   virtual bool has_explicit_ghosts() const {return true;}
 109
 110   /// return the total area, up to |y|<maxrap, that consists of
 111   /// unclustered ghosts
 112   virtual double empty_area(const RangeDefinition & range) const;
 113
 114   /// returns the total area under study
 115   double total_area () const;
 116
 117   /// returns the largest squared transverse momentum among
 118   /// all ghosts
 119   double max_ghost_perp2() const {return _max_ghost_perp2;}
 120
 121   /// returns true if there are any particles whose transverse momentum
 122   /// if so low that there's a risk of the ghosts having modified the
 123   /// clustering sequence
 124   bool has_dangerous_particles() const {return _has_dangerous_particles;}
 125
 126 private:
 127
 128   int    _n_ghosts;
 129   double _ghost_area;
 130   std::vector<bool> _is_pure_ghost;
 131   std::vector<double> _areas;
 132   std::vector<PseudoJet> _area_4vectors;
 133
 134   // things related to checks for dangerous particles
 135   double _max_ghost_perp2;
 136   bool   _has_dangerous_particles;
 137   static LimitedWarning _warnings;
 138
 139   //static int _n_warn_dangerous_particles;
 140   //static const int _max_warn_dangerous_particles = 5;
 141
 142
 143   unsigned int _initial_hard_n;
 144
 145   /// adds the "ghost" momenta, which will be used to estimate
 146   /// the jet area
 147   void _add_ghosts(const GhostedAreaSpec & area_spec);
 148
 149   /// another way of adding ghosts
 150   template<class L> void _add_ghosts (
 151           const std::vector<L> & ghosts,
 152           double                 ghost_area);
 153
 154   /// routine to be called after the processing is done so as to
 155   /// establish summary information on all the jets (areas, whether
 156   /// pure ghost, etc.)
 157   void _post_process();
 158
 159 };
 160
 161
 162 //----------------------------------------------------------------------
 163 // initialise from some generic type... Has to be made available
 164 // here in order for the template aspect of it to work...
 165 template<class L> void ClusterSequenceActiveAreaExplicitGhosts::_initialise
 166          (const std::vector<L> & pseudojets,
 167           const JetDefinition & jet_def,
 168           const GhostedAreaSpec * area_spec,
 169           const std::vector<L> * ghosts,
 170           double                 ghost_area,
 171           const bool & writeout_combinations) {
 172   // don't reserve space yet -- will be done below
 173
 174   // insert initial jets this way so that any type L that can be
 175   // converted to a pseudojet will work fine (basically PseudoJet
 176   // and any type that has [] subscript access to the momentum
 177   // components, such as CLHEP HepLorentzVector).
 178   for (unsigned int i = 0; i < pseudojets.size(); i++) {
 179     PseudoJet mom(pseudojets[i]);
 180     //mom.set_user_index(0); // for user's particles (user index now lost...)
 181     _jets.push_back(mom);
 182     _is_pure_ghost.push_back(false);
 183   }
 184
 185   _initial_hard_n = _jets.size();
 186
 187   if (area_spec != NULL) {
 188     _add_ghosts(*area_spec);
 189   } else {
 190     _add_ghosts(*ghosts, ghost_area);
 191   }
 192
 193   if (writeout_combinations) {
 194     std::cout << "# Printing particles including ghosts\n";
 195     for (unsigned j = 0; j < _jets.size(); j++) {
 196       printf("%5u %20.13f %20.13f %20.13e\n",
 197                j,_jets[j].rap(),_jets[j].phi_02pi(),_jets[j].kt2());
 198     }
 199     std::cout << "# Finished printing particles including ghosts\n";
 200   }
 201
 202   // this will ensure that we can still point to jets without
 203   // difficulties arising!
 204   _jets.reserve(_jets.size()*2);
 205
 206   // run the clustering
 207   _initialise_and_run(jet_def,writeout_combinations);
 208
 209   // set up all other information
 210   _post_process();
 211 }
 212
 213
 214 inline unsigned int ClusterSequenceActiveAreaExplicitGhosts::n_hard_particles() const {return _initial_hard_n;}
 215
 216
 217 //----------------------------------------------------------------------
 218 /// add an explicitly specified bunch of ghosts
 219 template<class L> void ClusterSequenceActiveAreaExplicitGhosts::_add_ghosts (
 220           const std::vector<L> & ghosts,
 221           double                 ghost_area) {
 222
 223
 224   for (unsigned i = 0; i < ghosts.size(); i++) {
 225     _is_pure_ghost.push_back(true);
 226     _jets.push_back(ghosts[i]);
 227   }
 228   // and record some info about ghosts
 229   _ghost_area = ghost_area;
 230   _n_ghosts   = ghosts.size();
 231 }
 232
 233
 234 } // fastjet namespace
 235
 236 #endif // __FASTJET_CLUSTERSEQUENCEACTIVEAREAEXPLICITGHOSTS_HH_