fastjetfortran_madfks.cc

Go to the documentation of this file.

//STARTHEADER
//
// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. 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, see <http://www.gnu.org/licenses/>.
//----------------------------------------------------------------------
//ENDHEADER
 
#include <iostream>
#include <memory>
 
#include "fastjet/ClusterSequence.hh"
#include "fastjet/ClusterSequenceArea.hh"
#include "fastjet/Selector.hh"
#include "fastjet/SISConePlugin.hh"
 
using namespace std;
using namespace fastjet;
 
FASTJET_BEGIN_NAMESPACE  // defined in fastjet/internal/base.hh
 
    namespace fwrapper {
  vector<PseudoJet> input_particles, jets;
  unique_ptr<JetDefinition::Plugin> plugin;
  JetDefinition jet_def;
  unique_ptr<ClusterSequence> cs;
 
  void transfer_input_particles(const double *p, const int &npart) {
    input_particles.resize(0);
    input_particles.reserve(npart);
    for (int i = 0; i < npart; i++) {
      valarray<double> mom(4);  // mom[0..3]
      for (int j = 0; j <= 3; j++) {
        mom[j] = *(p++);
        // RF-MZ: reorder the arguments because in madfks energy goes first; in fastjet last.
        // mom[(j+3) % 4] = *(p++);
      }
      PseudoJet psjet(mom);
      psjet.set_user_index(i);
      input_particles.push_back(psjet);
    }
  }
 
  void transfer_jets(double *f77jets, int &njets) {
    njets = jets.size();
    for (int i = 0; i < njets; i++) {
      for (int j = 0; j <= 3; j++) {
        *f77jets = jets[i][j];
        // RF-MZ: reorder the arguments because in madfks energy goes first; in fastjet last.
        //*f77jets = jets[i][(j+3) % 4];
        f77jets++;
      }
    }
  }
 
  void transfer_cluster_transfer(const double *p,
                                 const int &npart,
                                 const JetDefinition &jet_def,
                                 const double &ptmin,
                                 double *f77jets,
                                 int &njets,
                                 int *whichjet,
                                 const double &ghost_maxrap = 0.0,
                                 const int &nrepeat = 0,
                                 const double &ghost_area = 0.0) {
    // transfer p[4*ipart+0..3] -> input_particles[i]
    transfer_input_particles(p, npart);
 
    // perform the clustering
    if (ghost_maxrap == 0.0) {
      // cluster without areas
      cs = std::make_unique<ClusterSequence>(input_particles, jet_def);
    } else {
      // cluster with areas
      GhostedAreaSpec area_spec(ghost_maxrap, nrepeat, ghost_area);
      AreaDefinition area_def(active_area, area_spec);
      cs = std::make_unique<ClusterSequenceArea>(input_particles, jet_def, area_def);
    }
    // extract jets (pt-ordered)
    jets = sorted_by_pt(cs->inclusive_jets(ptmin));
 
    // transfer jets -> f77jets[4*ijet+0..3]
    transfer_jets(f77jets, njets);
 
    // Determine which parton/particle ended-up in which jet
    // set all jet entrie to zero first
    for (int ii = 0; ii < npart; ++ii)
      whichjet[ii] = 0;
 
    // Loop over jets and find constituents
    for (int kk = 0; kk < njets; ++kk) {
      vector<PseudoJet> constit = cs->constituents(jets[kk]);
      for (unsigned int ll = 0; ll < constit.size(); ++ll)
        whichjet[constit[ll].user_index()] = kk + 1;
    }
  }
}
FASTJET_END_NAMESPACE
 
using namespace fastjet::fwrapper;
 
extern "C" {
 
//
// Corresponds to the following Fortran subroutine
// interface structure:
//
//   SUBROUTINE FASTJETSISCONE(P,NPART,R,F,F77JETS,NJETS,WHICHJET)
//   DOUBLE PRECISION P(4,*), R, F, F77JETS(4,*)
//   INTEGER          NPART, NJETS, WHICHJET(*)
//
// where on input
//
//   P        the input particle 4-momenta
//   NPART    the number of input momenta
//   R        the radius parameter
//   F        the overlap threshold
//
// and on output
//
//   F77JETS  the output jet momenta (whose second dim should be >= NPART)
//            sorted in order of decreasing p_t.
//   NJETS    the number of output jets
//   WHICHJET(i) the jet of parton/particle 'i'
//
// NOTE: if you are interfacing fastjet to Pythia 6, Pythia stores its
// momenta as a matrix of the form P(4000,5), whereas this fortran
// interface to fastjet expects them as P(4,NPART), i.e. you must take
// the transpose of the Pythia array and drop the fifth component
// (particle mass).
//
void fastjetsiscone_(
    const double *p, const int &npart, const double &R, const double &f, double *f77jets, int &njets, int *whichjet) {
  // prepare jet def
  plugin = std::make_unique<SISConePlugin>(R, f);
  jet_def = plugin.get();
 
  // do everything
  //    transfer_cluster_transfer(p,npart,jet_def,f77jets,njets,whichjet);
}
 
//
// Corresponds to the following Fortran subroutine
// interface structure:
//
//   SUBROUTINE FASTJETSISCONEWITHAREA(P,NPART,R,F,GHMAXRAP,NREP,GHAREA,F77JETS,NJETS,WHICHJET)
//   DOUBLE PRECISION P(4,*), R, F, F77JETS(4,*), GHMAXRAP, GHAREA
//   INTEGER          NPART, NJETS, NREP, WHICHJET(*)
//
// where on input
//
//   P        the input particle 4-momenta
//   NPART    the number of input momenta
//   R        the radius parameter
//   F        the overlap threshold
//   GHMAXRAP the maximum (abs) rapidity covered by ghosts (FastJet default 6.0)
//   NREP     the number of repetitions used to evaluate the area (FastJet default 1)
//   GHAREA   the area of a single ghost (FastJet default 0.01)
//
// and on output
//
//   F77JETS  the output jet momenta (whose second dim should be >= NPART)
//            sorted in order of decreasing p_t.
//   NJETS    the number of output jets
//   WHICHJET(i) the jet of parton/particle 'i'
//
// NOTE: if you are interfacing fastjet to Pythia 6, Pythia stores its
// momenta as a matrix of the form P(4000,5), whereas this fortran
// interface to fastjet expects them as P(4,NPART), i.e. you must take
// the transpose of the Pythia array and drop the fifth component
// (particle mass).
//
void fastjetsisconewitharea_(const double *p,
                             const int &npart,
                             const double &R,
                             const double &f,
                             const double &ghost_rapmax,
                             const int &nrepeat,
                             const double &ghost_area,
                             double *f77jets,
                             int &njets,
                             int *whichjet) {
  // prepare jet def
  plugin = std::make_unique<SISConePlugin>(R, f);
  jet_def = plugin.get();
 
  // do everything
  //    transfer_cluster_transfer(p,npart,jet_def,f77jets,njets,whichjet,ghost_rapmax,nrepeat,ghost_area);
}
 
//
// Corresponds to the following Fortran subroutine
// interface structure:
//
//   SUBROUTINE FASTJETPPGENKT(P,NPART,R,PALG,F77JETS,NJETS,WHICHJET)
//   DOUBLE PRECISION P(4,*), R, PALG, F, F77JETS(4,*)
//   INTEGER          NPART, NJETS, WHICHJET(*)
//
// where on input
//
//   P        the input particle 4-momenta
//   NPART    the number of input momenta
//   R        the radius parameter
//   PALG     the power for the generalised kt alg
//            (1.0=kt, 0.0=C/A,  -1.0 = anti-kt)
//
// and on output
//
//   F77JETS  the output jet momenta (whose second dim should be >= NPART)
//            sorted in order of decreasing p_t.
//   NJETS    the number of output jets
//   WHICHJET(i) the jet of parton/particle 'i'
//
// For the values of PALG that correspond to "standard" cases (1.0=kt,
// 0.0=C/A, -1.0 = anti-kt) this routine actually calls the direct
// implementation of those algorithms, whereas for other values of
// PALG it calls the generalised kt implementation.
//
// NOTE: if you are interfacing fastjet to Pythia 6, Pythia stores its
// momenta as a matrix of the form P(4000,5), whereas this fortran
// interface to fastjet expects them as P(4,NPART), i.e. you must take
// the transpose of the Pythia array and drop the fifth component
// (particle mass).
//
void fastjetppgenkt_(const double *p,
                     const int &npart,
                     const double &R,
                     const double &ptjetmin,
                     const double &palg,
                     double *f77jets,
                     int &njets,
                     int *whichjet) {
  // prepare jet def
  if (palg == 1.0) {
    jet_def = JetDefinition(kt_algorithm, R);
  } else if (palg == 0.0) {
    jet_def = JetDefinition(cambridge_algorithm, R);
  } else if (palg == -1.0) {
    jet_def = JetDefinition(antikt_algorithm, R);
  } else {
    jet_def = JetDefinition(genkt_algorithm, R, palg);
  }
 
  // do everything
  transfer_cluster_transfer(p, npart, jet_def, ptjetmin, f77jets, njets, whichjet);
}
 
//
// Corresponds to the following Fortran subroutine
// interface structure:
//
//   SUBROUTINE FASTJETPPGENKTWITHAREA(P,NPART,R,PALG,GHMAXRAP,NREP,GHAREA,F77JETS,NJETS,WHICHJET)
//   DOUBLE PRECISION P(4,*), R, PALG, GHMAXRAP, GHAREA,  F77JETS(4,*)
//   INTEGER          NPART, NREP, NJETS, WHICHJET(*)
//
// where on input
//
//   P        the input particle 4-momenta
//   NPART    the number of input momenta
//   R        the radius parameter
//   PALG     the power for the generalised kt alg
//            (1.0=kt, 0.0=C/A,  -1.0 = anti-kt)
//   GHMAXRAP the maximum (abs) rapidity covered by ghosts (FastJet default 6.0)
//   NREP     the number of repetitions used to evaluate the area (FastJet default 1)
//   GHAREA   the area of a single ghost (FastJet default 0.01)
//
// and on output
//
//   F77JETS  the output jet momenta (whose second dim should be >= NPART)
//            sorted in order of decreasing p_t.
//   NJETS    the number of output jets
//   WHICHJET(i) the jet of parton/particle 'i'
//
// For the values of PALG that correspond to "standard" cases (1.0=kt,
// 0.0=C/A, -1.0 = anti-kt) this routine actually calls the direct
// implementation of those algorithms, whereas for other values of
// PALG it calls the generalised kt implementation.
//
// NOTE: if you are interfacing fastjet to Pythia 6, Pythia stores its
// momenta as a matrix of the form P(4000,5), whereas this fortran
// interface to fastjet expects them as P(4,NPART), i.e. you must take
// the transpose of the Pythia array and drop the fifth component
// (particle mass).
//
void fastjetppgenktwitharea_(const double *p,
                             const int &npart,
                             const double &R,
                             const double &palg,
                             const double &ghost_rapmax,
                             const int &nrepeat,
                             const double &ghost_area,
                             double *f77jets,
                             int &njets,
                             int *whichjet) {
  // prepare jet def
  if (palg == 1.0) {
    jet_def = JetDefinition(kt_algorithm, R);
  } else if (palg == 0.0) {
    jet_def = JetDefinition(cambridge_algorithm, R);
  } else if (palg == -1.0) {
    jet_def = JetDefinition(antikt_algorithm, R);
  } else {
    jet_def = JetDefinition(genkt_algorithm, R, palg);
  }
 
  // do everything
  //    transfer_cluster_transfer(p,npart,jet_def,f77jets,njets,whichjet,ghost_rapmax,nrepeat,ghost_area);
}
 
//
// Corresponds to the following Fortran subroutine
// interface structure:
//
//   SUBROUTINE FASTJETCONSTITUENTS(IJET,CONSTITUENT_INDICES,NCONSTITUENTS)
//   INTEGER    IJET
//   INTEGER    CONSTITUENT_INDICES(*)
//   INTEGER    nconstituents
//
void fastjetconstituents_(const int &ijet, int *constituent_indices, int &nconstituents) {
  assert(cs.get() != nullptr);
  assert(ijet > 0 && ijet <= int(jets.size()));
 
  vector<PseudoJet> constituents = cs->constituents(jets[ijet - 1]);
 
  nconstituents = constituents.size();
  for (int i = 0; i < nconstituents; i++) {
    constituent_indices[i] = constituents[i].cluster_hist_index() + 1;
  }
}
 
//
// Corresponds to the following Fortran subroutine
// interface structure:
//
//   FUNCTION FASTJETAREA(IJET)
//   DOUBLE PRECISION FASTJETAREA
//   INTEGER    IJET
//
double fastjetarea_(const int &ijet) {
  assert(ijet > 0 && ijet <= int(jets.size()));
  const ClusterSequenceAreaBase *csab = dynamic_cast<const ClusterSequenceAreaBase *>(cs.get());
  if (csab != nullptr) {
    // we have areas and can use csab to access all the area-related info
    return csab->area(jets[ijet - 1]);
  } else {
    return 0.;
    //  Error("No area information associated to this jet.");
  }
}
 
//
// Corresponds to the following Fortran interface
//
//   FUNCTION FASTJETDMERGE(N)
//   DOUBLE PRECISION FASTJETDMERGE
//   INTEGER N
//
double fastjetdmerge_(const int &n) {
  assert(cs.get() != nullptr);
  return cs->exclusive_dmerge(n);
}
 
//
// Corresponds to the following Fortran interface
//
//   FUNCTION FASTJETDMERGEMAX(N)
//   DOUBLE PRECISION FASTJETDMERGEMAX
//   INTEGER N
//
double fastjetdmergemax_(const int &n) {
  assert(cs.get() != nullptr);
  return cs->exclusive_dmerge_max(n);
}
 
//
// Corresponds to the following Fortran interface
//
//   SUBROUTINE FASTJETGLOBALRHOANDSIGMA(RAPMIN,RAPMAX,PHIMIN,PHIMAX,RHO,SIGMA,MEANAREA)
//   DOUBLE PRECISION RAPMIN,RAPMAX,PHIMIN,PHIMAX
//   DOUBLE PRECISION RHO,SIGMA,MEANAREA
//
void fastjetglobalrhoandsigma_(const double &rapmin,
                               const double &rapmax,
                               const double &phimin,
                               const double &phimax,
                               double &rho,
                               double &sigma,
                               double &meanarea) {
  const ClusterSequenceAreaBase *csab = dynamic_cast<const ClusterSequenceAreaBase *>(cs.get());
  if (csab != nullptr) {
    // we have areas and can use csab to access all the area-related info
    Selector range = SelectorRapRange(rapmin, rapmax) && SelectorPhiRange(phimin, phimax);
    bool use_area_4vector = false;
    csab->get_median_rho_and_sigma(range, use_area_4vector, rho, sigma, meanarea);
  } else {
    Error("Clustering with area is necessary in order to be able to evaluate rho.");
  }
}
}