ClusterSequenceActiveArea.hh

//STARTHEADER
// $Id: ClusterSequenceActiveArea.hh 1761 2010-09-16 10:43:18Z soyez $
//
// Copyright (c) 2005-2006, Matteo Cacciari and Gavin Salam
//
//----------------------------------------------------------------------
// This file is part of FastJet.
//
//  FastJet is free software; you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation; either version 2 of the License, or
//  (at your option) any later version.
//
//  The algorithms that underlie FastJet have required considerable
//  development and are described in hep-ph/0512210. If you use
//  FastJet as part of work towards a scientific publication, please
//  include a citation to the FastJet paper.
//
//  FastJet is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.
//
//  You should have received a copy of the GNU General Public License
//  along with FastJet; if not, write to the Free Software
//  Foundation, Inc.:
//      59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//----------------------------------------------------------------------
//ENDHEADER

#ifndef __FASTJET_CLUSTERSEQUENCEACTIVEAREA_HH__
#define __FASTJET_CLUSTERSEQUENCEACTIVEAREA_HH__


#include "fastjet/PseudoJet.hh"
#include "fastjet/ClusterSequenceAreaBase.hh"
#include "fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh"
#include<iostream>
#include<vector>

//------------ backwards compatibility with version 2.1 -------------
// for backwards compatibility make ActiveAreaSpec name available
#include "fastjet/ActiveAreaSpec.hh"
#include "fastjet/ClusterSequenceWithArea.hh"
//--------------------------------------------------------------------


FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh

//using namespace std;

class ClusterSequenceActiveArea : public ClusterSequenceAreaBase {
public:

  ClusterSequenceActiveArea() {}

  template<class L> ClusterSequenceActiveArea
         (const std::vector<L> & pseudojets, 
          const JetDefinition & jet_def,
          const GhostedAreaSpec & ghost_spec,
          const bool & writeout_combinations = false) ;

  virtual double area (const PseudoJet & jet) const {
                             return _average_area[jet.cluster_hist_index()];};
  virtual double area_error (const PseudoJet & jet) const {
                             return _average_area2[jet.cluster_hist_index()];};

  virtual PseudoJet area_4vector (const PseudoJet & jet) const {
                    return _average_area_4vector[jet.cluster_hist_index()];};

  enum mean_pt_strategies{median=0, non_ghost_median, pttot_over_areatot, 
                          pttot_over_areatot_cut, mean_ratio_cut, play,
                          median_4vector};

  //      ClusterSequenceAreaBase class instead
  double pt_per_unit_area(mean_pt_strategies strat=median, 
                          double range=2.0 ) const;

  // following code removed -- now dealt with by AreaBase class (and
  // this definition here conflicts with it).
//   /// fits a form pt_per_unit_area(y) = a + b*y^2 in the range
//   /// abs(y)<raprange (for negative raprange, it defaults to
//   /// _safe_rap_for_area).
//   void parabolic_pt_per_unit_area(double & a,double & b, double raprange=-1.0,
//                                double exclude_above=-1.0, 
//                                bool use_area_4vector=false ) const;
// 
  virtual double empty_area(const RangeDefinition & range) const;

  virtual double n_empty_jets(const RangeDefinition & range) const;

protected:
  void _resize_and_zero_AA ();
  void _initialise_AA(const JetDefinition & jet_def,
                      const GhostedAreaSpec & ghost_spec,
                      const bool & writeout_combinations,
                      bool & continue_running);

  void _run_AA(const GhostedAreaSpec & ghost_spec);

  void _postprocess_AA(const GhostedAreaSpec & ghost_spec);

  void _initialise_and_run_AA (const JetDefinition & jet_def,
                               const GhostedAreaSpec & ghost_spec,
                               const bool & writeout_combinations = false);

  void _transfer_ghost_free_history(
           const ClusterSequenceActiveAreaExplicitGhosts & clust_seq);


  void _transfer_areas(const std::vector<int> & unique_hist_order, 
                       const ClusterSequenceActiveAreaExplicitGhosts & );

  std::valarray<double> _average_area, _average_area2;
  std::valarray<PseudoJet> _average_area_4vector;

  bool has_dangerous_particles() const {return _has_dangerous_particles;}

private:


  double           _non_jet_area, _non_jet_area2, _non_jet_number;

  double _maxrap_for_area; // max rap where we put ghosts
  double _safe_rap_for_area; // max rap where we trust jet areas

  bool   _has_dangerous_particles; 


  void _extract_tree(std::vector<int> &) const;
  void _extract_tree_children(int pos, std::valarray<bool> &, const std::valarray<int> &, std::vector<int> &) const;
  void _extract_tree_parents (int pos, std::valarray<bool> &, const std::valarray<int> &,  std::vector<int> &) const;

  void _throw_unless_jets_have_same_perp_or_E(const PseudoJet & jet, 
                                              const PseudoJet & refjet, 
                                              double tolerance,
             const ClusterSequenceActiveAreaExplicitGhosts & jets_ghosted_seq
                                              ) const;

  //static int _n_seed_warnings;
  //const static int _max_seed_warnings = 10;

  // record the number of repeats
  int _ghost_spec_repeat;

  class GhostJet : public PseudoJet {
  public:
    GhostJet(const PseudoJet & j, double a) : PseudoJet(j), area(a){}
    double area;
  };

  std::vector<GhostJet> _ghost_jets;
  std::vector<GhostJet> _unclustered_ghosts;
};




template<class L> ClusterSequenceActiveArea::ClusterSequenceActiveArea 
(const std::vector<L> & pseudojets, 
 const JetDefinition & jet_def,
 const GhostedAreaSpec & ghost_spec,
 const bool & writeout_combinations) {

  // transfer the initial jets (type L) into our own array
  _transfer_input_jets(pseudojets);

  // run the clustering for active areas
  _initialise_and_run_AA(jet_def, ghost_spec, writeout_combinations);

}


  
FASTJET_END_NAMESPACE

#endif // __FASTJET_CLUSTERSEQUENCEACTIVEAREA_HH__