ClusterSequenceActiveArea.cc

Go to the documentation of this file.

//STARTHEADER
// $Id: ClusterSequenceActiveArea.cc 514 2007-03-10 17:50:06Z salam $
//
// 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

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

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


using namespace std;


//int ClusterSequenceActiveArea::_n_seed_warnings = 0;
//const int _max_seed_warnings = 10;

void ClusterSequenceActiveArea::_initialise_and_run_AA (
                const JetDefinition & jet_def,
                const ActiveAreaSpec & area_spec,
                const bool & writeout_combinations) 
{

  // initialize our local area information
  _average_area.resize(2*_jets.size());  _average_area  = 0.0;
  _average_area2.resize(2*_jets.size()); _average_area2 = 0.0;
  _average_area_4vector.resize(2*_jets.size()); 
  _average_area_4vector = PseudoJet(0.0,0.0,0.0,0.0);
  _non_jet_area = 0.0; _non_jet_area2 = 0.0; _non_jet_number=0.0;
     
  // for future reference...
  _maxrap_for_area = area_spec.ghost_maxrap();
  _safe_rap_for_area = _maxrap_for_area - jet_def.R();

  // Make sure we'll have at least one repetition -- then we can
  // deduce the unghosted clustering sequence from one of the ghosted
  // sequences. If we do not have any repetitions, then get the
  // unghosted sequence from the plain unghosted clustering.
  //
  // NB: all decanting and filling of initial history will then
  // be carried out by base-class routine
  if (area_spec.repeat() <= 0) {
    _initialise_and_run(jet_def, writeout_combinations);
    return;
  }

  // transfer all relevant info into internal variables
  _decant_options(jet_def, writeout_combinations);

  // set up the history entries for the initial particles (those
  // currently in _jets)
  _fill_initial_history();

  // record the input jets as they are currently
  vector<PseudoJet> input_jets(_jets);

  // code for testing the unique tree
  vector<int> unique_tree;

  
  

  // run the clustering multiple times so as to get areas of all the jets
  for (int irepeat = 0; irepeat < area_spec.repeat(); irepeat++) {

    ClusterSequenceActiveAreaExplicitGhosts clust_seq(input_jets, 
                                                      jet_def, area_spec);

    if (irepeat == 0) {
      // take the non-ghost part of the history and put into our own
      // history.
      _transfer_ghost_free_history(clust_seq);
      // get the "unique" order that will be used for transferring all areas. 
      unique_tree = unique_history_order();
    }

    // transfer areas from clust_seq into our object
    _transfer_areas(unique_tree, clust_seq);
  }
  
  _average_area  /= area_spec.repeat();
  _average_area2 /= area_spec.repeat();
  if (area_spec.repeat() > 1) {
    _average_area2 = sqrt(abs(_average_area2 - _average_area*_average_area)/
                          (area_spec.repeat()-1));
  } else {
    _average_area2 = 0.0;
  }

  _non_jet_area  /= area_spec.repeat();
  _non_jet_area2 /= area_spec.repeat();
  _non_jet_area2  = sqrt(abs(_non_jet_area2 - _non_jet_area*_non_jet_area)/
                         area_spec.repeat());
  _non_jet_number /= area_spec.repeat();

  // following bizarre way of writing things is related to 
  // poverty of operations on PseudoJet objects (as well as some confusion
  // in one or two places)
  for (unsigned i = 0; i < _average_area_4vector.size(); i++) {
    _average_area_4vector[i] = (1.0/area_spec.repeat()) * _average_area_4vector[i];
  }
  //cerr << "Non-jet area = " << _non_jet_area << " +- " << _non_jet_area2<<endl;

  
}


//----------------------------------------------------------------------
double ClusterSequenceActiveArea::pt_per_unit_area(
                       mean_pt_strategies strat, double range) const {
  
  vector<PseudoJet> incl_jets = inclusive_jets();
  vector<double> pt_over_areas;

  for (unsigned i = 0; i < incl_jets.size(); i++) {
    if (abs(incl_jets[i].rap()) < _safe_rap_for_area) {
      double this_area;
      if ( strat == median_4vector ) {
          this_area = area_4vector(incl_jets[i]).perp();
      } else {
          this_area = area(incl_jets[i]);
      }
      pt_over_areas.push_back(incl_jets[i].perp()/this_area);
    }
  }
  
  // there is nothing inside our region, so answer will always be zero
  if (pt_over_areas.size() == 0) {return 0.0;}
  
  // get median (pt/area) [this is the "old" median definition. It considers
  // only the "real" jets in calculating the median, i.e. excluding the
  // only-ghost ones]
  sort(pt_over_areas.begin(), pt_over_areas.end());
  double non_ghost_median_ratio = pt_over_areas[pt_over_areas.size()/2];

  // new median definition that takes into account non-jet area (i.e.
  // jets composed only of ghosts), and for fractional median position 
  // interpolates between the corresponding entries in the pt_over_areas array
  double nj_median_pos = (pt_over_areas.size()-1 - _non_jet_number)/2.0;
  double nj_median_ratio;
  if (nj_median_pos >= 0 && pt_over_areas.size() > 1) {
    int int_nj_median = int(nj_median_pos);
    nj_median_ratio = 
      pt_over_areas[int_nj_median] * (int_nj_median+1-nj_median_pos)
      + pt_over_areas[int_nj_median+1] * (nj_median_pos - int_nj_median);
  } else {
    nj_median_ratio = 0.0;
  }


  // get various forms of mean (pt/area)
  double pt_sum = 0.0, pt_sum_with_cut = 0.0;
  double area_sum = _non_jet_area, area_sum_with_cut = _non_jet_area;
  double ratio_sum = 0.0; 
  double ratio_n = _non_jet_number;
  for (unsigned i = 0; i < incl_jets.size(); i++) {
    if (abs(incl_jets[i].rap()) < _safe_rap_for_area) {
      double this_area;
      if ( strat == median_4vector ) {
          this_area = area_4vector(incl_jets[i]).perp();
      } else {
          this_area = area(incl_jets[i]);
      }
      pt_sum   += incl_jets[i].perp();
      area_sum += this_area;
      double ratio = incl_jets[i].perp()/this_area;
      if (ratio < range*nj_median_ratio) {
        pt_sum_with_cut   += incl_jets[i].perp();
        area_sum_with_cut += this_area;
        ratio_sum += ratio; ratio_n++;
      }
    }
  }
  
  if (strat == play) {
    double trunc_sum = 0, trunc_sumsqr = 0;
    vector<double> means(pt_over_areas.size()), sd(pt_over_areas.size());
    for (unsigned i = 0; i < pt_over_areas.size() ; i++ ) {
      double ratio = pt_over_areas[i];
      trunc_sum += ratio;
      trunc_sumsqr += ratio*ratio;
      means[i] = trunc_sum / (i+1);
      sd[i]    = sqrt(abs(means[i]*means[i]  - trunc_sumsqr/(i+1)));
      cerr << "i, means, sd: " <<i<<", "<< means[i] <<", "<<sd[i]<<", "<<
        sd[i]/sqrt(i+1.0)<<endl;
    }
    cout << "-----------------------------------"<<endl;
    for (unsigned i = 0; i <= pt_over_areas.size()/2 ; i++ ) {
      cout << "Median "<< i <<" = " << pt_over_areas[i]<<endl;
    }
    cout << "Number of non-jets: "<<_non_jet_number<<endl;
    cout << "Area of non-jets: "<<_non_jet_area<<endl;
    cout << "Default median position: " << (pt_over_areas.size()-1)/2.0<<endl;
    cout << "NJ median position: " << nj_median_pos <<endl;
    cout << "NJ median value: " << nj_median_ratio <<endl;
    return 0.0;
  }

  switch(strat) {
  case median:
  case median_4vector:
    return nj_median_ratio;
  case non_ghost_median:
    return non_ghost_median_ratio; 
  case pttot_over_areatot:
    return pt_sum / area_sum;
  case pttot_over_areatot_cut:
    return pt_sum_with_cut / area_sum_with_cut;
  case mean_ratio_cut:
    return ratio_sum/ratio_n;
  default:
    return nj_median_ratio;
  }

}


//----------------------------------------------------------------------
// fit a parabola to pt/area as a function of rapidity, using the
// formulae of CCN28-36 (which actually fits f = a+b*x^2)
void ClusterSequenceActiveArea::parabolic_pt_per_unit_area(
       double & a, double & b, double raprange, double exclude_above,
       bool use_area_4vector) const {
  
  double this_raprange;
  if (raprange <= 0) {this_raprange = _safe_rap_for_area;}
  else {this_raprange = raprange;}

  int n=0;
  int n_excluded = 0;
  double mean_f=0, mean_x2=0, mean_x4=0, mean_fx2=0; 

  vector<PseudoJet> incl_jets = inclusive_jets();

  for (unsigned i = 0; i < incl_jets.size(); i++) {
    if (abs(incl_jets[i].rap()) < this_raprange) {
      double this_area;
      if ( use_area_4vector ) {
          this_area = area_4vector(incl_jets[i]).perp();     
      } else {
          this_area = area(incl_jets[i]);
      }
      double f = incl_jets[i].perp()/this_area;
      if (exclude_above <= 0.0 || f < exclude_above) {
        double x = incl_jets[i].rap(); double x2 = x*x;
        mean_f   += f;
        mean_x2  += x2;
        mean_x4  += x2*x2;
        mean_fx2 += f*x2;
        n++;
      } else {
        n_excluded++;
      }
    }
  }

  if (n <= 1) {
    // meaningful results require at least two jets inside the
    // area -- mind you if there are empty jets we should be in 
    // any case doing something special...
    a = 0.0;
    b = 0.0;
  } else {
    mean_f   /= n;
    mean_x2  /= n;
    mean_x4  /= n;
    mean_fx2 /= n;
    
    b = (mean_f*mean_x2 - mean_fx2)/(mean_x2*mean_x2 - mean_x4);
    a = mean_f - b*mean_x2;
  }
  //cerr << "n_excluded = "<< n_excluded << endl;
}


//----------------------------------------------------------------------
void ClusterSequenceActiveArea::_transfer_ghost_free_history(
             const ClusterSequenceActiveAreaExplicitGhosts & ghosted_seq) {
  
  const vector<history_element> & gs_history  = ghosted_seq.history();
  vector<int> gs2self_hist_map(gs_history.size());

  // work our way through to first non-trivial combination
  unsigned igs = 0;
  unsigned iself = 0;
  while (gs_history[igs].parent1 == InexistentParent) {
    // record correspondence 
    if (!ghosted_seq.is_pure_ghost(igs)) {
      gs2self_hist_map[igs] = iself++; 
    } else {
      gs2self_hist_map[igs] = Invalid; 
    }
    igs++;
  };

  // make sure the count of non-ghost initial jets is equal to
  // what we already have in terms of initial jets
  assert(iself == _history.size());

  // now actually transfer things
  do  {
    // if we are a pure ghost, then go on to next round
    if (ghosted_seq.is_pure_ghost(igs)) {
      gs2self_hist_map[igs] = Invalid;
      continue;
    }

    const history_element & gs_hist_el = gs_history[igs];

    bool parent1_is_ghost = ghosted_seq.is_pure_ghost(gs_hist_el.parent1);
    bool parent2_is_ghost = ghosted_seq.is_pure_ghost(gs_hist_el.parent2);

    // if exactly one parent is a ghost then maintain info about the
    // non-ghost correspondence for this jet, and then go on to next
    // recombination in the ghosted sequence
    if (parent1_is_ghost && !parent2_is_ghost && gs_hist_el.parent2 >= 0) {
      gs2self_hist_map[igs] = gs2self_hist_map[gs_hist_el.parent2];
      continue;
    }
    if (!parent1_is_ghost && parent2_is_ghost) {
      gs2self_hist_map[igs] = gs2self_hist_map[gs_hist_el.parent1];
      continue;
    }

    // no parents are ghosts...
    if (gs_hist_el.parent2 >= 0) {
      // recombination of two non-ghosts
      gs2self_hist_map[igs] = _history.size();
      // record the recombination in our own sequence
      int newjet_k; // dummy var -- not used
      //cerr << igs << " " << gs_hist_el.parent1 << " " << gs_hist_el.parent2 << endl;
      //cerr << gs2self_hist_map[gs_hist_el.parent1] << " " << gs2self_hist_map[gs_hist_el.parent2] << endl;
      int jet_i = _history[gs2self_hist_map[gs_hist_el.parent1]].jetp_index;
      int jet_j = _history[gs2self_hist_map[gs_hist_el.parent2]].jetp_index;
      //cerr << "recombining "<< jet_i << " and "<< jet_j << endl;
      _do_ij_recombination_step(jet_i, jet_j, gs_hist_el.dij, newjet_k);
    } else {
      // we have a non-ghost that has become a beam-jet
      assert(gs_history[igs].parent2 == BeamJet);
      // record position
      gs2self_hist_map[igs] = _history.size();
      // record the recombination in our own sequence
      _do_iB_recombination_step(
             _history[gs2self_hist_map[gs_hist_el.parent1]].jetp_index,
             gs_hist_el.dij);
    }
  } while (++igs < gs_history.size());

  // finally transfer info about strategy used (which isn't necessarily
  // always the one that got asked for...)
  _strategy = ghosted_seq.strategy_used();
}

//----------------------------------------------------------------------
void ClusterSequenceActiveArea::_transfer_areas(
            const vector<int> & unique_hist_order,
            const ClusterSequenceActiveAreaExplicitGhosts & ghosted_seq  ) {

  const vector<history_element> & gs_history  = ghosted_seq.history();
  const vector<PseudoJet>       & gs_jets     = ghosted_seq.jets();
  vector<int>    gs_unique_hist_order = ghosted_seq.unique_history_order();

  const double tolerance = 1e-13; // to decide when two jets are the same

  int j = -1;
  int hist_index = -1;
  
  valarray<double> our_areas(_history.size());
  our_areas = 0.0;

  valarray<PseudoJet> our_area_4vectors(_history.size());
  our_area_4vectors = PseudoJet(0.0,0.0,0.0,0.0);

  for (unsigned i = 0; i < gs_history.size(); i++) {
    // only consider composite particles
    unsigned gs_hist_index = gs_unique_hist_order[i];
    if (gs_hist_index < ghosted_seq.n_particles()) continue;
    const history_element & gs_hist = gs_history[gs_unique_hist_order[i]];
    int parent1 = gs_hist.parent1;
    int parent2 = gs_hist.parent2;

    if (parent2 == BeamJet) {
      // need to look at parent to get the actual jet
      const PseudoJet & jet = 
          gs_jets[gs_history[parent1].jetp_index];
      double area = ghosted_seq.area(jet);
      PseudoJet ext_area = ghosted_seq.area_4vector(jet);

      if (ghosted_seq.is_pure_ghost(parent1)) {
        if (abs(jet.rap()) < _safe_rap_for_area) {
          _non_jet_area  += area;
          _non_jet_area2 += area*area;
          _non_jet_number += 1;
        }
      } else {

        // get next "combined-particle" index in our own history
        // making sure we don't go beyond it's bounds (if we do
        // then we're in big trouble anyway...)
        while (++j < static_cast<int>(_history.size())) {
          hist_index = unique_hist_order[j];
          if (hist_index >= _initial_n) break;}

        // sanity check 
        const PseudoJet & refjet = 
          _jets[_history[_history[hist_index].parent1].jetp_index];
        //if (jet.perp2() != refjet.perp2()) {
        //if (abs(jet.perp2()-refjet.perp2()) > 
        //            tolerance*max(jet.perp2(),refjet.perp2())) {

        // If pt disagrees check E; if they both disagree there's a
        // problem here... NB: a massive particle with zero pt may
        // have its pt changed when a ghost is added -- this is why we
        // also require the energy to be wrong before complaining
        if (abs(jet.perp2()-refjet.perp2()) > 
                    tolerance*max(jet.perp2(),refjet.perp2())
            && abs(jet.E()-refjet.E()) > tolerance*max(jet.E(),refjet.E())) {
          cerr << jet.perp() << " " << refjet.perp() << " "
               << jet.perp() - refjet.perp() << endl;
          cerr << refjet.px() << " " 
               << refjet.py() << " " 
               << refjet.pz() << " " 
               << refjet.E() << endl;
          throw Error("Could not match clustering sequence for an inclusive jet when reconstructing areas"); }

        // set the area at this clustering stage
        our_areas[hist_index]  = area; 
        our_area_4vectors[hist_index]  = ext_area; 

        // update the parent as well -- that way its area is the area
        // immediately before clustering (i.e. resolve an ambiguity in
        // the Cambridge case and ensure in the kt case that the original
        // particles get a correct area)
        our_areas[_history[hist_index].parent1] = area;
        our_area_4vectors[_history[hist_index].parent1] = ext_area;
        
      }
    }
    else if (!ghosted_seq.is_pure_ghost(parent1) && 
             !ghosted_seq.is_pure_ghost(parent2)) {

      // get next "combined-particle" index in our own history
      while (++j < static_cast<int>(_history.size())) {
        hist_index = unique_hist_order[j];
        if (hist_index >= _initial_n) break;}
      
      const PseudoJet & jet = gs_jets[gs_hist.jetp_index];
      const PseudoJet & refjet = _jets[_history[hist_index].jetp_index];

      // run sanity check 
      if (abs(jet.perp2()-refjet.perp2()) > 
          tolerance*max(jet.perp2(),refjet.perp2())) {
          cerr << jet.perp() << " " << refjet.perp() << " "<< jet.perp() - refjet.perp() << endl;
          throw Error("Could not match clustering sequence for an exclusive jet when reconstructing areas"); }

      // update area and our local index (maybe redundant since later
      // the descendants will reupdate it?)
      double area  = ghosted_seq.area(jet);
      our_areas[hist_index]  += area; 

      PseudoJet ext_area = ghosted_seq.area_4vector(jet);
      //our_area_4vectors[hist_index] = our_area_4vectors[hist_index] + ext_area; 
      _jet_def.recombiner()->plus_equal(our_area_4vectors[hist_index], ext_area);

      // now update areas of parents (so that they becomes areas
      // immediately before clustering occurred). This is of use
      // because it allows us to set the areas of the original hard
      // particles in the kt algorithm; for the Cambridge case it
      // means a jet's area will be the area just before it clusters
      // with another hard jet.
      const PseudoJet & jet1 = gs_jets[gs_history[parent1].jetp_index];
      int our_parent1 = _history[hist_index].parent1;
      our_areas[our_parent1] = ghosted_seq.area(jet1);
      our_area_4vectors[our_parent1] = ghosted_seq.area_4vector(jet1);

      const PseudoJet & jet2 = gs_jets[gs_history[parent2].jetp_index];
      int our_parent2 = _history[hist_index].parent2;
      our_areas[our_parent2] = ghosted_seq.area(jet2);
      our_area_4vectors[our_parent2] = ghosted_seq.area_4vector(jet2);
    }

  }

  _average_area  += our_areas; 
  _average_area2 += our_areas*our_areas; 

  //_average_area_4vector += our_area_4vectors;
  // Use the proper recombination scheme when averaging the area_4vectors
  // over multiple ghost runs (i.e. the repeat stage);
  for (unsigned i = 0; i < _average_area_4vector.size(); i++) {
    _jet_def.recombiner()->plus_equal(_average_area_4vector[i],
                                       our_area_4vectors[i]);
  }
}



FASTJET_END_NAMESPACE

---