PxConePlugin is a plugin for fastjet (v2.1 upwards) that provides an interface to the fortran pxcone iterative cone algorithm with midpoint seeds. More...

#include <PxConePlugin.hh>

Inheritance diagram for PxConePlugin:

[legend]

Collaboration diagram for PxConePlugin:

[legend]

List of all members.

Public Member Functions
	PxConePlugin (double cone_radius, double min_jet_energy=5.0, double overlap_threshold=0.5, bool E_scheme_jets=false)
	constructor for the PxConePlugin, whose arguments have the following meaning:
double	cone_radius () const
	the cone radius
double	min_jet_energy () const
	minimum jet energy (protojets below this are thrown own before merging/splitting) -- called epslon in pxcone
double	overlap_threshold () const
	Maximum fraction of overlap energy in a jet -- called ovlim in pxcone.
bool	E_scheme_jets () const
	if true then the final jets are returned as the E-scheme recombination of the particle momenta (by default, pxcone returns massless jets with a mean phi,eta type of recombination); regardless of what is returned, the internal pxcone jet-finding procedure is unaffected.
virtual std::string	description () const
	return a textual description of the jet-definition implemented in this plugin
virtual void	run_clustering (ClusterSequence &) const
	given a ClusterSequence that has been filled up with initial particles, the following function should fill up the rest of the ClusterSequence, using the following member functions of ClusterSequence: plugin_do_ij_recombination(...) plugin_do_iB_recombination(...)
virtual double	R () const
	the plugin mechanism's standard way of accessing the jet radius
Private Attributes
double	_cone_radius
double	_min_jet_energy
double	_overlap_threshold
bool	_E_scheme_jets

Detailed Description

PxConePlugin is a plugin for fastjet (v2.1 upwards) that provides an interface to the fortran pxcone iterative cone algorithm with midpoint seeds.

Pxcone was written by Luis del Pozo and Michael H. Seymour. It is not a "supported" program, so if you encounter problems, you are on your own...

Note that pxcone sometimes encounters non-stable iterations; in such cases it returns an error -- the plugin propagates this by throwing a fastjet::Error exception; if the user wishes to have robust code, they should catch this exception.

Pxcone has a hard-coded limit (by default 4000) on the maximum number of particles and protojets; if the number of particles or protojets exceeds this, again a fastjet::Error exception will be thrown.

The functionality of pxcone is described at http://www.hep.man.ac.uk/u/wplano/ConeJet.ps

Definition at line 65 of file PxConePlugin.hh.

Constructor & Destructor Documentation

PxConePlugin::PxConePlugin	(	double	cone_radius,
		double	min_jet_energy = `5.0`,
		double	overlap_threshold = `0.5`,
		bool	E_scheme_jets = `false`
	)		`[inline]`

constructor for the PxConePlugin, whose arguments have the following meaning:

the cone_radius is as usual in cone algorithms

stables cones (protojets) below min_jet_energy are discarded before calling the splitting procedure to resolve overlaps (called epslon in pxcone).

when two protojets overlap, if (overlapping_Et)/(Et_of_softer_protojet) < overlap_threshold the overlapping energy is split between the two protojets; otherwise the less energetic protojet is discarded. Called ovlim in pxcone.

pxcone carries out p-scheme recombination, and the resulting jets are massless; setting E_scheme_jets = true (default false) doesn't change the jet composition, but the final momentum sum for the jets is carried out by direct four-vector addition instead of p-scheme recombination.

Definition at line 89 of file PxConePlugin.hh.

                                               : 
    _cone_radius        (cone_radius        ),
    _min_jet_energy     (min_jet_energy     ),
    _overlap_threshold  (overlap_threshold  ),
    _E_scheme_jets      (E_scheme_jets      ) {}

Member Function Documentation

double PxConePlugin::cone_radius ( ) const [inline]

the cone radius

Definition at line 102 of file PxConePlugin.hh.

References _cone_radius.

Referenced by R().

{return _cone_radius        ;}

string PxConePlugin::description ( ) const [virtual]

return a textual description of the jet-definition implemented in this plugin

Implements JetDefinition::Plugin.

Definition at line 44 of file PxConePlugin.cc.

                                        {
  ostringstream desc;
  
  desc << "PxCone jet algorithm with " 
       << "cone_radius = "        << cone_radius        () << ", "
       << "min_jet_energy = "     << min_jet_energy     () << ", "
       << "overlap_threshold  = " << overlap_threshold  () << ", "
       << "E_scheme_jets  = "     << E_scheme_jets      () 
       << " (NB: non-standard version of PxCone, containing small bug fixes by Gavin Salam)";

  return desc.str();
}

bool PxConePlugin::E_scheme_jets ( ) const [inline]

if true then the final jets are returned as the E-scheme recombination of the particle momenta (by default, pxcone returns massless jets with a mean phi,eta type of recombination); regardless of what is returned, the internal pxcone jet-finding procedure is unaffected.

Definition at line 116 of file PxConePlugin.hh.

References _E_scheme_jets.

{return _E_scheme_jets      ;}

double PxConePlugin::min_jet_energy ( ) const [inline]

minimum jet energy (protojets below this are thrown own before merging/splitting) -- called epslon in pxcone

Definition at line 106 of file PxConePlugin.hh.

References _min_jet_energy.

{return _min_jet_energy     ;}

double PxConePlugin::overlap_threshold ( ) const [inline]

Maximum fraction of overlap energy in a jet -- called ovlim in pxcone.

Definition at line 109 of file PxConePlugin.hh.

References _overlap_threshold.

{return _overlap_threshold  ;}

virtual double PxConePlugin::R ( ) const [inline, virtual]

the plugin mechanism's standard way of accessing the jet radius

Implements JetDefinition::Plugin.

Definition at line 123 of file PxConePlugin.hh.

References cone_radius().

{return cone_radius();}

void PxConePlugin::run_clustering ( ClusterSequence & ) const [virtual]

given a ClusterSequence that has been filled up with initial particles, the following function should fill up the rest of the ClusterSequence, using the following member functions of ClusterSequence:

plugin_do_ij_recombination(...)
plugin_do_iB_recombination(...)

Implements JetDefinition::Plugin.

Definition at line 58 of file PxConePlugin.cc.

References ClusterSequence::jets(), ClusterSequence::plugin_record_iB_recombination(), ClusterSequence::plugin_record_ij_recombination(), and pxcone.

                                                                   {
 
  // only have hh mode
  int mode = 2;

  int    ntrak = clust_seq.jets().size(), itkdm = 4;
  double *ptrak = new double[ntrak*4+1];
  for (int i = 0; i < ntrak; i++) {
    ptrak[4*i+0] = clust_seq.jets()[i].px();
    ptrak[4*i+1] = clust_seq.jets()[i].py();
    ptrak[4*i+2] = clust_seq.jets()[i].pz();
    ptrak[4*i+3] = clust_seq.jets()[i].E();
  }  

  // max number of allowed jets
  int mxjet = ntrak;
  int njet;
  double *pjet  = new double[mxjet*5+1];
  int    *ipass = new int[ntrak+1];
  int    *ijmul = new int[mxjet+1];
  int ierr;

  // run pxcone
  pxcone(
    mode   ,    // 1=>e+e-, 2=>hadron-hadron
    ntrak  ,    // Number of particles
    itkdm  ,    // First dimension of PTRAK array: 
    ptrak  ,    // Array of particle 4-momenta (Px,Py,Pz,E)
    cone_radius()  ,    // Cone size (half angle) in radians
    min_jet_energy() ,    // Minimum Jet energy (GeV)
    overlap_threshold()  ,    // Maximum fraction of overlap energy in a jet
    mxjet  ,    // Maximum possible number of jets
    njet   ,    // Number of jets found
    pjet ,       // 5-vectors of jets
    ipass,      // Particle k belongs to jet number IPASS(k)-1
                // IPASS = -1 if not assosciated to a jet
    ijmul,      // Jet i contains IJMUL[i] particles
    ierr        // = 0 if all is OK ;   = -1 otherwise
    );

  if (ierr != 0) throw Error("An error occurred while running PXCONE");

  // now transfer information back 
  valarray<int> last_index_created(njet);

  vector<vector<int> > jet_particle_content(njet);

  // get a list of particles in each jet
  for (int itrak = 0; itrak < ntrak; itrak++) {
    int jet_i = ipass[itrak] - 1;
    if (jet_i >= 0) jet_particle_content[jet_i].push_back(itrak);
  }

  // now transfer the jets back into our own structure -- we will
  // mimic the cone code with a sequential recombination sequence in
  // which the jets are built up by adding one particle at a time
  for(int ipxjet = njet-1; ipxjet >= 0; ipxjet--) {
    const vector<int> & jet_trak_list = jet_particle_content[ipxjet];
    int jet_k = jet_trak_list[0];
  
    for (unsigned ilist = 1; ilist < jet_trak_list.size(); ilist++) {
      int jet_i = jet_k;
      // retrieve our misappropriated index for the jet
      int jet_j = jet_trak_list[ilist];
      // do a fake recombination step with dij=0
      double dij = 0.0;
      //clust_seq.plugin_record_ij_recombination(jet_i, jet_j, dij, jet_k);
      if (ilist != jet_trak_list.size()-1 || E_scheme_jets()) {
        // our E-scheme recombination in cases where it doesn't matter
        clust_seq.plugin_record_ij_recombination(jet_i, jet_j, dij, jet_k);
      } else {
        // put in pxcone's momentum for the last recombination so that the
        // final inclusive jet corresponds exactly to PXCONE's
        clust_seq.plugin_record_ij_recombination(jet_i, jet_j, dij, 
                                PseudoJet(pjet[5*ipxjet+0],pjet[5*ipxjet+1],
                                          pjet[5*ipxjet+2],pjet[5*ipxjet+3]),
                                                 jet_k);
      }
    }
  
    // NB: put a sensible looking d_iB just to be nice...
    double d_iB = clust_seq.jets()[jet_k].perp2();
    clust_seq.plugin_record_iB_recombination(jet_k, d_iB);
  }


  //cout << endl;
  //for (int ijet = 0; ijet < njet; ijet++) {
  //  PseudoJet jet(pjet[ijet][0],pjet[ijet][1],pjet[ijet][2],pjet[ijet][3]);
  //  cout << jet.perp() << " " << jet.rap() << endl;
  //}
  //cout << "-----------------------------------------------------\n";
  //vector<PseudoJet> ourjets(clust_seq.inclusive_jets());
  //for (vector<PseudoJet>::const_iterator ourjet = ourjets.begin();
  //     ourjet != ourjets.end(); ourjet++) {
  //  cout << ourjet->perp() << " " << ourjet->rap() << endl;
  //}

  delete[] ptrak;
  delete[] ipass;
  delete[] ijmul;
  delete[] pjet;
}