include/fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh

0001 //FJSTARTHEADER
0002 // $Id$
0003 //
0004 // Copyright (c) 2005-2021, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
0005 //
0006 //----------------------------------------------------------------------
0007 // This file is part of FastJet.
0008 //
0009 //  FastJet is free software; you can redistribute it and/or modify
0010 //  it under the terms of the GNU General Public License as published by
0011 //  the Free Software Foundation; either version 2 of the License, or
0012 //  (at your option) any later version.
0013 //
0014 //  The algorithms that underlie FastJet have required considerable
0015 //  development. They are described in the original FastJet paper,
0016 //  hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
0017 //  FastJet as part of work towards a scientific publication, please
0018 //  quote the version you use and include a citation to the manual and
0019 //  optionally also to hep-ph/0512210.
0020 //
0021 //  FastJet is distributed in the hope that it will be useful,
0022 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
0023 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0024 //  GNU General Public License for more details.
0025 //
0026 //  You should have received a copy of the GNU General Public License
0027 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
0028 //----------------------------------------------------------------------
0029 //FJENDHEADER
0030
0031 #ifndef __FASTJET_CLUSTERSEQUENCEACTIVEAREAEXPLICITGHOSTS_HH_
0032 #define __FASTJET_CLUSTERSEQUENCEACTIVEAREAEXPLICITGHOSTS_HH_
0033
0034 #include "fastjet/PseudoJet.hh"
0035 #include "fastjet/ClusterSequenceAreaBase.hh"
0036 #include "fastjet/GhostedAreaSpec.hh"
0037 #include "fastjet/LimitedWarning.hh"
0038 #include<iostream>
0039 #include<vector>
0040 #include <cstdio>
0041
0042 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
0043
0044 //======================================================================
0045 /// @ingroup sec_area_classes
0046 /// \class ClusterSequenceActiveAreaExplicitGhosts
0047 /// Like ClusterSequence with computation of the active jet area with the
0048 /// addition of explicit ghosts
0049 ///
0050 /// Class that behaves essentially like ClusterSequence except
0051 /// that it also provides access to the area of a jet (which
0052 /// will be a random quantity... Figure out what to do about seeds
0053 /// later...)
0054 ///
0055 /// This class should not be used directly. Rather use
0056 /// ClusterSequenceArea with the appropriate AreaDefinition
0057 class ClusterSequenceActiveAreaExplicitGhosts :
0058   public ClusterSequenceAreaBase {
0059 public:
0060   /// constructor using a GhostedAreaSpec to specify how the area is
0061   /// to be measured
0062   template<class L> ClusterSequenceActiveAreaExplicitGhosts
0063          (const std::vector<L> & pseudojets,
0064           const JetDefinition & jet_def_in,
0065       const GhostedAreaSpec & ghost_spec,
0066       const bool & writeout_combinations = false)
0067        : ClusterSequenceAreaBase() {
0068            std::vector<L> * ghosts = NULL;
0069        _initialise(pseudojets,jet_def_in,&ghost_spec,ghosts,0.0,
0070                        writeout_combinations); }
0071
0072   template<class L> ClusterSequenceActiveAreaExplicitGhosts
0073          (const std::vector<L> & pseudojets,
0074           const JetDefinition & jet_def_in,
0075           const std::vector<L> & ghosts,
0076           double ghost_area,
0077       const bool & writeout_combinations = false)
0078        : ClusterSequenceAreaBase() {
0079            const GhostedAreaSpec * ghost_spec = NULL;
0080        _initialise(pseudojets,jet_def_in,ghost_spec,&ghosts,ghost_area,
0081                        writeout_combinations); }
0082
0083
0084   /// does the actual work of initialisation
0085   template<class L> void _initialise
0086          (const std::vector<L> & pseudojets,
0087           const JetDefinition & jet_def_in,
0088       const GhostedAreaSpec * ghost_spec,
0089       const std::vector<L> * ghosts,
0090       double                 ghost_area,
0091       const bool & writeout_combinations);
0092
0093   //vector<PseudoJet> constituents (const PseudoJet & jet) const;
0094
0095   /// returns the number of hard particles (i.e. those supplied by the user).
0096   unsigned int n_hard_particles() const;
0097
0098   /// returns the area of a jet
0099   virtual double area (const PseudoJet & jet) const FASTJET_OVERRIDE;
0100
0101   /// returns a four vector corresponding to the sum (E-scheme) of the
0102   /// ghost four-vectors composing the jet area, normalised such that
0103   /// for a small contiguous area the p_t of the extended_area jet is
0104   /// equal to area of the jet.
0105   virtual PseudoJet area_4vector (const PseudoJet & jet) const FASTJET_OVERRIDE;
0106
0107   /// true if a jet is made exclusively of ghosts
0108   virtual bool is_pure_ghost(const PseudoJet & jet) const FASTJET_OVERRIDE;
0109
0110   /// true if the entry in the history index corresponds to a
0111   /// ghost; if hist_ix does not correspond to an actual particle
0112   /// (i.e. hist_ix < 0), then the result is false.
0113   bool is_pure_ghost(int history_index) const;
0114
0115   /// this class does have explicit ghosts
0116   virtual bool has_explicit_ghosts() const FASTJET_OVERRIDE {return true;}
0117
0118   /// return the total area, corresponding to a given Selector, that
0119   /// consists of unclustered ghosts
0120   ///
0121   /// The selector needs to apply jet by jet
0122   virtual double empty_area(const Selector & selector) const FASTJET_OVERRIDE;
0123
0124   /// returns the total area under study
0125   double total_area () const;
0126
0127   /// returns the largest squared transverse momentum among
0128   /// all ghosts
0129   double max_ghost_perp2() const {return _max_ghost_perp2;}
0130
0131   /// returns true if there are any particles whose transverse momentum
0132   /// if so low that there's a risk of the ghosts having modified the
0133   /// clustering sequence
0134   bool has_dangerous_particles() const {return _has_dangerous_particles;}
0135
0136   /// get the area of the ghosts
0137   //double ghost_area() const{return _ghost_area;}
0138
0139 private:
0140
0141   int    _n_ghosts;
0142   double _ghost_area;
0143   std::vector<bool> _is_pure_ghost;
0144   std::vector<double> _areas;
0145   std::vector<PseudoJet> _area_4vectors;
0146
0147   // things related to checks for dangerous particles
0148   double _max_ghost_perp2;
0149   bool   _has_dangerous_particles;
0150   static LimitedWarning _warnings;
0151
0152   //static int _n_warn_dangerous_particles;
0153   //static const int _max_warn_dangerous_particles = 5;
0154
0155
0156   unsigned int _initial_hard_n;
0157
0158   /// adds the "ghost" momenta, which will be used to estimate
0159   /// the jet area
0160   void _add_ghosts(const GhostedAreaSpec & ghost_spec);
0161
0162   /// another way of adding ghosts
0163   template<class L> void _add_ghosts (
0164       const std::vector<L> & ghosts,
0165       double                 ghost_area);
0166
0167   /// routine to be called after the processing is done so as to
0168   /// establish summary information on all the jets (areas, whether
0169   /// pure ghost, etc.)
0170   void _post_process();
0171
0172 };
0173
0174
0175 //----------------------------------------------------------------------
0176 // initialise from some generic type... Has to be made available
0177 // here in order for the template aspect of it to work...
0178 template<class L> void ClusterSequenceActiveAreaExplicitGhosts::_initialise
0179          (const std::vector<L> & pseudojets,
0180           const JetDefinition & jet_def_in,
0181       const GhostedAreaSpec * ghost_spec,
0182       const std::vector<L> * ghosts,
0183       double                 ghost_area,
0184       const bool & writeout_combinations) {
0185   // don't reserve space yet -- will be done below
0186
0187   // insert initial jets this way so that any type L that can be
0188   // converted to a pseudojet will work fine (basically PseudoJet
0189   // and any type that has [] subscript access to the momentum
0190   // components, such as CLHEP HepLorentzVector).
0191   for (unsigned int i = 0; i < pseudojets.size(); i++) {
0192     PseudoJet mom(pseudojets[i]);
0193     //mom.set_user_index(0); // for user's particles (user index now lost...)
0194     _jets.push_back(mom);
0195     _is_pure_ghost.push_back(false);
0196   }
0197
0198   _initial_hard_n = _jets.size();
0199
0200   if (ghost_spec != NULL) {
0201     //std::cout << "about to reserve " << (_jets.size()+ghost_spec->n_ghosts())*2 << std::endl;
0202     _jets.reserve((_jets.size()+ghost_spec->n_ghosts()));
0203     _add_ghosts(*ghost_spec);
0204   } else {
0205     _jets.reserve(_jets.size()+ghosts->size());
0206     _add_ghosts(*ghosts, ghost_area);
0207   }
0208
0209   if (writeout_combinations) {
0210     std::cout << "# Printing particles including ghosts\n";
0211     for (unsigned j = 0; j < _jets.size(); j++) {
0212       printf("%5u %20.13f %20.13f %20.13e\n",
0213            j,_jets[j].rap(),_jets[j].phi_02pi(),_jets[j].kt2());
0214     }
0215     std::cout << "# Finished printing particles including ghosts\n";
0216   }
0217
0218   // this will ensure that we can still point to jets without
0219   // difficulties arising!
0220   //std::cout << _jets.size() << " " << _jets.size()*2 << " " << _jets.max_size() << std::endl;
0221   _jets.reserve(_jets.size()*2); //GPS tmp removed
0222
0223   // run the clustering
0224   _initialise_and_run(jet_def_in,writeout_combinations);
0225
0226   // set up all other information
0227   _post_process();
0228 }
0229
0230
0231 inline unsigned int ClusterSequenceActiveAreaExplicitGhosts::n_hard_particles() const {return _initial_hard_n;}
0232
0233
0234 //----------------------------------------------------------------------
0235 /// add an explicitly specified bunch of ghosts
0236 template<class L> void ClusterSequenceActiveAreaExplicitGhosts::_add_ghosts (
0237       const std::vector<L> & ghosts,
0238       double                 ghost_area) {
0239
0240
0241   for (unsigned i = 0; i < ghosts.size(); i++) {
0242     _is_pure_ghost.push_back(true);
0243     _jets.push_back(ghosts[i]);
0244   }
0245   // and record some info about ghosts
0246   _ghost_area = ghost_area;
0247   _n_ghosts   = ghosts.size();
0248 }
0249
0250
0251 FASTJET_END_NAMESPACE
0252
0253 #endif // __FASTJET_CLUSTERSEQUENCEACTIVEAREAEXPLICITGHOSTS_HH_