2
// $Id: fastjetfortran.cc 1570 2009-05-25 10:45:18Z salam $
4
// Copyright (c) 2005-2007, Matteo Cacciari, Gavin Salam and Gregory Soyez
6
//----------------------------------------------------------------------
7
// This file is part of FastJet.
9
// FastJet is free software; you can redistribute it and/or modify
10
// it under the terms of the GNU General Public License as published by
11
// the Free Software Foundation; either version 2 of the License, or
12
// (at your option) any later version.
14
// The algorithms that underlie FastJet have required considerable
15
// development and are described in hep-ph/0512210. If you use
16
// FastJet as part of work towards a scientific publication, please
17
// include a citation to the FastJet paper.
19
// FastJet is distributed in the hope that it will be useful,
20
// but WITHOUT ANY WARRANTY; without even the implied warranty of
21
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22
// GNU General Public License for more details.
24
// You should have received a copy of the GNU General Public License
25
// along with FastJet; if not, write to the Free Software
27
// 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
28
//----------------------------------------------------------------------
32
#include "fastjet/ClusterSequence.hh"
34
namespace fj = fastjet;
36
bool _first_time = true;
37
auto_ptr<fj::ClusterSequence> cs;
41
// f77 interface to the pp generalised-kt (sequential recombination)
42
// algorithms, as defined in arXiv.org:0802.1189, which includes
43
// kt, Cambridge/Aachen and anti-kt as special cases.
45
// Corresponds to the following Fortran subroutine
46
// interface structure:
48
// SUBROUTINE FASTJETPPSEQREC(P,NPART,R,PALG,F77JETS,NJETS)
49
// DOUBLE PRECISION P(4,*), R, PALG, F, F77JETS(4,*)
50
// INTEGER NPART, NJETS
54
// P the input particle 4-momenta
55
// NPART the number of input momenta
56
// R the radius parameter
57
// PALG the power for the generalised kt alg
58
// (1.0=kt, 0.0=C/A, -1.0 = anti-kt)
62
// F77JETS the output jet momenta (whose second dim should be >= NPART)
63
// sorted in order of decreasing p_t.
64
// NJETS the number of output jets
66
// For the values of PALG that correspond to "standard" cases (1.0=kt,
67
// 0.0=C/A, -1.0 = anti-kt) this routine actually calls the direct
68
// implementation of those algorithms, whereas for other values of
69
// PALG it calls the generalised kt implementation.
71
void fastjetppgenkt_(const double * p, const int & npart,
72
const double & R, const double & jetptmin, const double & palg,
73
double * f77jets, int & njets, int * whichjet) {
75
// transfer p[4*ipart+0..3] -> input_particles[i]
76
vector<fj::PseudoJet> input_particles;
77
for (int i=0; i<npart; i++) {
78
valarray<double> mom(4); // mom[0..3]
79
for (int j=0;j<=3; j++) {
81
// RF: reorder the arguments because in madfks energy goes first; in fastjet last.
82
mom[(j+3) % 4] = *(p++);
84
fj::PseudoJet psjet(mom);
85
psjet.set_user_index(i);
86
input_particles.push_back(psjet);
89
// prepare jet def and run fastjet
90
fj::JetDefinition jet_def;
92
jet_def = fj::JetDefinition(fj::kt_algorithm, R);
93
} else if (palg == 0.0) {
94
jet_def = fj::JetDefinition(fj::cambridge_algorithm, R);
95
} else if (palg == -1.0) {
96
jet_def = fj::JetDefinition(fj::antikt_algorithm, R);
98
jet_def = fj::JetDefinition(fj::genkt_algorithm, R, palg);
102
// perform clustering
103
cs.reset(new fj::ClusterSequence(input_particles,jet_def));
104
// extract jets (pt-ordered)
105
vector<fj::PseudoJet> jets = sorted_by_pt(cs->inclusive_jets(jetptmin));
108
// tell the user what was done
110
cout << "# " << jet_def.description() << endl;
111
cout <<("# and minimum jet pt of %E \n",jetptmin);
112
cout <<("#------------------------------------------------------------------------- \n");
116
// transfer jets -> f77jets[4*ijet+0..3]
117
for (int i=0; i<njets; i++) {
118
for (int j=0;j<=3; j++) {
119
// *f77jets = jets[i][j];
120
// RF: reorder the arguments because in madfks energy goes first; in fastjet last.
121
*f77jets = jets[i][(j+3) % 4];
128
// set all jet entrie to zero first
129
for(unsigned int ii=0; ii<npart; ++ii) whichjet[ii]=0;
131
for (unsigned int kk=0; kk<jets.size(); ++kk) {
132
vector<fj::PseudoJet> constit = cs->constituents(jets[kk]);
133
for(unsigned int ll=0; ll<constit.size(); ++ll)
134
whichjet[constit[ll].user_index()]=kk+1;
139
/// return the dmin corresponding to the recombination that went from
140
/// n+1 to n jets (sometimes known as d_{n n+1}).
142
// Corresponds to the following Fortran interface
144
// FUNCTION FASTJETDMERGE(N)
145
// DOUBLE PRECISION FASTJETDMERGE
148
double fastjetdmerge_(const int & n) {
149
assert(cs.get() != 0);
150
return cs->exclusive_dmerge(n);
153
/// return the maximum of the dmin encountered during all recombinations
154
/// up to the one that led to an n-jet final state; identical to
155
/// exclusive_dmerge, except in cases where the dmin do not increase
158
// Corresponds to the following Fortran interface
160
// FUNCTION FASTJETDMERGEMAX(N)
161
// DOUBLE PRECISION FASTJETDMERGEMAX
164
double fastjetdmergemax_(const int & n) {
165
assert(cs.get() != 0);
166
return cs->exclusive_dmerge_max(n);
added
removed
Template/FKS-born/SubProcesses/fastjetfortran_madfks.cc.fj-2.4.xx
