JadeDistanceTagger.cc

//STARTHEADER
// $Id: MassDropTagger.cc 2138 2011-05-16 18:56:46Z soyez $
//
// Copyright (c) 2005-2011, Matteo Cacciari, Gavin Salam and Gregory Soyez
//
//----------------------------------------------------------------------
// 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/tools/JadeDistanceTagger.hh>
#include <sstream>

FASTJET_BEGIN_NAMESPACE

using namespace std;

//----------------------------------------------------------------------
// JadeDistanceTagger implementation
//----------------------------------------------------------------------

//------------------------------------------------------------------------
// description of the tagger
string JadeDistanceTagger::description() const{ 
  ostringstream oss;
  oss << "JadeDistanceTagger with mdrop=" << _mdrop << " and mcut=" << _mcut;
  return oss.str();
}


//------------------------------------------------------------------------
// the tagging itself
//  - jet   the PseudoJet to tag
PseudoJet JadeDistanceTagger::apply(const PseudoJet & jet) const{
  // search the subjets with the mass drop and mass cut parameters
  vector<PseudoJet> subjets;
  _recursively_find_subjets(jet,subjets);

  // stupidly make sure we have at least 2!
  if (subjets.size()<2) return join(PseudoJet(0.0,0.0,0.0,0.0));

  // among all possible pairs, maximise the (modified) Jade distance
  double dmax=0.0;
  unsigned int i1=0, i2=0;
  for (unsigned int it1=0; it1<subjets.size(); it1++){
    const PseudoJet & j1 = subjets[it1];
    double pt1 = j1.perp();
    for (unsigned int it2=it1+1; it2<subjets.size(); it2++){
      const PseudoJet & j2 = subjets[it2];
      double d = j1.squared_distance(j2);
      d = pt1*j2.perp()*d*d; // could improve slightly by store the pt... or using the ^2

      if (d>dmax){ i1 = it1; i1=it2; dmax=d;}
    }
  }
  const PseudoJet & j1 = subjets[i1];
  const PseudoJet & j2 = subjets[i2];

  // apply the filter
  PseudoJet result = join<StructureType>(j1,j2);
  result.extra_properties<JadeDistanceTagger>()._mdrop 
    = max(j1.m(), j2.m())/jet.m();
  result.extra_properties<JadeDistanceTagger>()._jade_distance = dmax;
  return result;
}


//------------------------------------------------------------------------
void JadeDistanceTagger::_recursively_find_subjets(const PseudoJet &jet, 
                                                   vector<PseudoJet> &subjets) const{
  // check if we should simply keep it
  PseudoJet j1, j2;

  if ((jet.m()<_mcut) || (!jet.has_parents(j1, j2))){
    subjets.push_back(jet);
    return;
  }

  // make j1 the more massive jet
  if (j1.m() < j2.m()) std::swap(j1,j2);

  if (j1.m() < _mdrop*jet.m())
    _recursively_find_subjets(j2, subjets);

  _recursively_find_subjets(j1, subjets);
}

FASTJET_END_NAMESPACE