JETAN/fastjet/fastjet/ClusterSequenceActiveArea.hh

   1 //STARTHEADER
   2 // $Id: ClusterSequenceActiveArea.hh 1134 2008-03-15 17:05:16Z 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_CLUSTERSEQUENCEACTIVEAREA_HH__
  32 #define __FASTJET_CLUSTERSEQUENCEACTIVEAREA_HH__
  33
  34
  35 #include "fastjet/PseudoJet.hh"
  36 #include "fastjet/ClusterSequenceAreaBase.hh"
  37 #include "fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh"
  38 #include<iostream>
  39 #include<vector>
  40
  41 //------------ backwards compatibility with version 2.1 -------------
  42 // for backwards compatibility make ActiveAreaSpec name available
  43 #include "fastjet/ActiveAreaSpec.hh"
  44 #include "fastjet/ClusterSequenceWithArea.hh"
  45 //--------------------------------------------------------------------
  46
  47
  48 namespace fastjet {      // defined in fastjet/internal/base.hh
  49
  50 //using namespace std;
  51
  52 /// Class that behaves essentially like ClusterSequence except
  53 /// that it also provides access to the area of a jet (which
  54 /// will be a random quantity... Figure out what to do about seeds
  55 /// later...)
  56 class ClusterSequenceActiveArea : public ClusterSequenceAreaBase {
  57 public:
  58
  59   /// default constructor
  60   ClusterSequenceActiveArea() {}
  61
  62   /// constructor based on JetDefinition and GhostedAreaSpec
  63   template<class L> ClusterSequenceActiveArea
  64          (const std::vector<L> & pseudojets,
  65           const JetDefinition & jet_def,
  66           const GhostedAreaSpec & area_spec,
  67           const bool & writeout_combinations = false) ;
  68
  69   virtual double area (const PseudoJet & jet) const {
  70                              return _average_area[jet.cluster_hist_index()];};
  71   virtual double area_error (const PseudoJet & jet) const {
  72                              return _average_area2[jet.cluster_hist_index()];};
  73
  74   virtual PseudoJet area_4vector (const PseudoJet & jet) const {
  75                     return _average_area_4vector[jet.cluster_hist_index()];};
  76
  77   /// enum providing a variety of tentative strategies for estimating
  78   /// the background (e.g. non-jet) activity in a highly populated event; the
  79   /// one that has been most extensively tested is median.
  80   enum mean_pt_strategies{median=0, non_ghost_median, pttot_over_areatot,
  81                           pttot_over_areatot_cut, mean_ratio_cut, play,
  82                           median_4vector};
  83
  84   /// return the transverse momentum per unit area according to one
  85   /// of the above strategies; for some strategies (those with "cut"
  86   /// in their name) the parameter "range" allows one to exclude a
  87   /// subset of the jets for the background estimation, those that
  88   /// have pt/area > median(pt/area)*range.
  89   ///
  90   /// NB: This call is OBSOLETE; use media_pt_per_unit_area from the
  91   //      ClusterSequenceAreaBase class instead
  92   double pt_per_unit_area(mean_pt_strategies strat=median,
  93                           double range=2.0 ) const;
  94
  95   // following code removed -- now dealt with by AreaBase class (and
  96   // this definition here conflicts with it).
  97 //   /// fits a form pt_per_unit_area(y) = a + b*y^2 in the range
  98 //   /// abs(y)<raprange (for negative raprange, it defaults to
  99 //   /// _safe_rap_for_area).
 100 //   void parabolic_pt_per_unit_area(double & a,double & b, double raprange=-1.0,
 101 //                                double exclude_above=-1.0,
 102 //                                bool use_area_4vector=false ) const;
 103 //
 104   /// rewrite the empty area from the parent class, so as to use
 105   /// all info at our disposal
 106   /// return the total area, in the given y-phi range, that consists of ghost
 107   /// jets or unclustered ghosts
 108   virtual double empty_area(const RangeDefinition & range) const;
 109
 110   /// return the true number of empty jets (replaces
 111   /// ClusterSequenceAreaBase::n_empty_jets(...))
 112   virtual double n_empty_jets(const RangeDefinition & range) const;
 113
 114 protected:
 115   void _resize_and_zero_AA ();
 116   void _initialise_AA(const JetDefinition & jet_def,
 117                       const GhostedAreaSpec & area_spec,
 118                       const bool & writeout_combinations,
 119                       bool & continue_running);
 120
 121   void _run_AA(const GhostedAreaSpec & area_spec);
 122
 123   void _postprocess_AA(const GhostedAreaSpec & area_spec);
 124
 125   /// does the initialisation and running specific to the active
 126   /// areas class
 127   void _initialise_and_run_AA (const JetDefinition & jet_def,
 128                                const GhostedAreaSpec & area_spec,
 129                                const bool & writeout_combinations = false);
 130
 131   /// transfer the history (and jet-momenta) from clust_seq to our
 132   /// own internal structure while removing ghosts
 133   void _transfer_ghost_free_history(
 134            const ClusterSequenceActiveAreaExplicitGhosts & clust_seq);
 135
 136
 137   /// transfer areas from the ClusterSequenceActiveAreaExplicitGhosts
 138   /// object into our internal area bookkeeping...
 139   void _transfer_areas(const std::vector<int> & unique_hist_order,
 140                        const ClusterSequenceActiveAreaExplicitGhosts & );
 141
 142   /// child classes benefit from having these at their disposal
 143   std::valarray<double> _average_area, _average_area2;
 144   std::valarray<PseudoJet> _average_area_4vector;
 145
 146   /// returns true if there are any particles whose transverse momentum
 147   /// if so low that there's a risk of the ghosts having modified the
 148   /// clustering sequence
 149   bool has_dangerous_particles() const {return _has_dangerous_particles;}
 150
 151 private:
 152
 153
 154   double           _non_jet_area, _non_jet_area2, _non_jet_number;
 155
 156   double _maxrap_for_area; // max rap where we put ghosts
 157   double _safe_rap_for_area; // max rap where we trust jet areas
 158
 159   bool   _has_dangerous_particles;
 160
 161
 162   /// routine for extracting the tree in an order that will be independent
 163   /// of any degeneracies in the recombination sequence that don't
 164   /// affect the composition of the final jets
 165   void _extract_tree(std::vector<int> &) const;
 166   /// do the part of the extraction associated with pos, working
 167   /// through its children and their parents
 168   void _extract_tree_children(int pos, std::valarray<bool> &, const std::valarray<int> &, std::vector<int> &) const;
 169   /// do the part of the extraction associated with the parents of pos.
 170   void _extract_tree_parents (int pos, std::valarray<bool> &, const std::valarray<int> &,  std::vector<int> &) const;
 171
 172   /// check if two jets have the same momentum to within the
 173   /// tolerance (and if pt's are not the same we're forgiving and
 174   /// look to see if the energy is the same)
 175   void _throw_unless_jets_have_same_perp_or_E(const PseudoJet & jet,
 176                                               const PseudoJet & refjet,
 177                                               double tolerance,
 178              const ClusterSequenceActiveAreaExplicitGhosts & jets_ghosted_seq
 179                                               ) const;
 180
 181   /// since we are playing nasty games with seeds, we should warn
 182   /// the user a few times
 183   //static int _n_seed_warnings;
 184   //const static int _max_seed_warnings = 10;
 185
 186   // record the number of repeats
 187   int _area_spec_repeat;
 188
 189   /// a class for our internal storage of ghost jets
 190   class GhostJet : public PseudoJet {
 191   public:
 192     GhostJet(const PseudoJet & j, double a) : PseudoJet(j), area(a){}
 193     double area;
 194   };
 195
 196   std::vector<GhostJet> _ghost_jets;
 197   std::vector<GhostJet> _unclustered_ghosts;
 198 };
 199
 200
 201
 202
 203 template<class L> ClusterSequenceActiveArea::ClusterSequenceActiveArea
 204 (const std::vector<L> & pseudojets,
 205  const JetDefinition & jet_def,
 206  const GhostedAreaSpec & area_spec,
 207  const bool & writeout_combinations) {
 208
 209   // transfer the initial jets (type L) into our own array
 210   _transfer_input_jets(pseudojets);
 211
 212   // run the clustering for active areas
 213   _initialise_and_run_AA(jet_def, area_spec, writeout_combinations);
 214
 215 }
 216
 217
 218
 219 } // fastjet namespace
 220
 221 #endif // __FASTJET_CLUSTERSEQUENCEACTIVEAREA_HH__