Fork me on GitHub

source: git/external/fastjet/internal/ClusterSequence_N2.icc@ 1a4956e

ImprovedOutputFile Timing
Last change on this file since 1a4956e was b7b836a, checked in by Pavel Demin <pavel-demin@…>, 6 years ago

update FastJet library to 3.3.1 and FastJet Contrib library to 1.036

  • Property mode set to 100644
File size: 6.5 KB
Line 
1// -*- C++ -*-
2#ifndef __FASTJET_CLUSTERQUENCE_N2_ICC__
3#define __FASTJET_CLUSTERQUENCE_N2_ICC__
4#include "fastjet/ClusterSequence.hh"
5
6//FJSTARTHEADER
7// $Id: ClusterSequence_N2.cc 1351 2009-01-09 18:03:03Z salam $
8//
9// Copyright (c) 2005-2018, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
10//
11//----------------------------------------------------------------------
12// This file is part of FastJet.
13//
14// FastJet is free software; you can redistribute it and/or modify
15// it under the terms of the GNU General Public License as published by
16// the Free Software Foundation; either version 2 of the License, or
17// (at your option) any later version.
18//
19// The algorithms that underlie FastJet have required considerable
20// development. They are described in the original FastJet paper,
21// hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
22// FastJet as part of work towards a scientific publication, please
23// quote the version you use and include a citation to the manual and
24// optionally also to hep-ph/0512210.
25//
26// FastJet is distributed in the hope that it will be useful,
27// but WITHOUT ANY WARRANTY; without even the implied warranty of
28// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
29// GNU General Public License for more details.
30//
31// You should have received a copy of the GNU General Public License
32// along with FastJet. If not, see <http://www.gnu.org/licenses/>.
33//----------------------------------------------------------------------
34//FJENDHEADER
35
36//----------------------------------------------------------------------
37/// Order(N^2) clustering
38///
39/// Works for any class BJ that satisfies certain minimal
40/// requirements (which are ...?)
41///
42/// - need to have _bj_set_jetinfo
43/// - need to have _bj_dist
44/// - should contain members kt2 (=energy^2), NN, NN_dist, _jets_index
45
46FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
47
48// this should not normally appear in an include file, but we'll make an
49// exception seeing as this is
50//using namespace std;
51
52
53template<class BJ> void ClusterSequence::_simple_N2_cluster() {
54 int n = _jets.size();
55 BJ * briefjets = new BJ[n];
56 BJ * jetA = briefjets, * jetB;
57
58 // initialise the basic jet info
59 for (int i = 0; i< n; i++) {
60 _bj_set_jetinfo(jetA, i);
61 jetA++; // move on to next entry of briefjets
62 }
63 BJ * tail = jetA; // a semaphore for the end of briefjets
64 BJ * head = briefjets; // a nicer way of naming start
65
66 // now initialise the NN distances: jetA will run from 1..n-1; and
67 // jetB from 0..jetA-1
68 for (jetA = head + 1; jetA != tail; jetA++) {
69 // set NN info for jetA based on jets running from head..jetA-1,
70 // checking in the process whether jetA itself is an undiscovered
71 // NN of one of those jets.
72 _bj_set_NN_crosscheck(jetA, head, jetA);
73 }
74
75
76 // now create the diJ (where J is i's NN) table -- remember that
77 // we differ from standard normalisation here by a factor of R2
78 double * diJ = new double[n];
79 jetA = head;
80 for (int i = 0; i < n; i++) {
81 diJ[i] = _bj_diJ(jetA);
82 jetA++; // have jetA follow i
83 }
84
85 // now run the recombination loop
86 int history_location = n-1;
87 while (tail != head) {
88
89 // find the minimum of the diJ on this round
90 double diJ_min = diJ[0];
91 int diJ_min_jet = 0;
92 for (int i = 1; i < n; i++) {
93 if (diJ[i] < diJ_min) {diJ_min_jet = i; diJ_min = diJ[i];}
94 }
95
96 // do the recombination between A and B
97 history_location++;
98 jetA = & briefjets[diJ_min_jet];
99 // GPS mod 2009-02-11
100 //jetB = jetA->NN;
101 jetB = static_cast<BJ *>(jetA->NN);
102 // put the normalisation back in
103 diJ_min *= _invR2;
104 if (jetB != NULL) {
105 // jet-jet recombination
106 // If necessary relabel A & B to ensure jetB < jetA, that way if
107 // the larger of them == newtail then that ends up being jetA and
108 // the new jet that is added as jetB is inserted in a position that
109 // has a future!
110 if (jetA < jetB) {std::swap(jetA,jetB);}
111
112 int nn; // new jet index
113 _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
114
115 // what was jetB will now become the new jet
116 _bj_set_jetinfo(jetB, nn);
117
118 } else {
119 // jet-beam recombination
120 _do_iB_recombination_step(jetA->_jets_index, diJ_min);
121 }
122
123 // now update our nearest neighbour info and diJ table
124 // first reduce size of table
125 tail--; n--;
126 // Copy last jet contents and diJ info into position of jetA
127 *jetA = *tail;
128 diJ[jetA - head] = diJ[tail-head];
129
130 // Initialise jetB's NN distance as well as updating it for
131 // other particles.
132 // NB: by having different loops for jetB == or != NULL we could
133 // perhaps save a few percent (usually avoid one if inside loop),
134 // but will not do it for now because on laptop fluctuations are
135 // too large to reliably measure a few percent difference...
136 for (BJ * jetI = head; jetI != tail; jetI++) {
137 // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
138 if (jetI->NN == jetA || jetI->NN == jetB) {
139 _bj_set_NN_nocross(jetI, head, tail);
140 diJ[jetI-head] = _bj_diJ(jetI); // update diJ
141 }
142 // check whether new jetB is closer than jetI's current NN and
143 // if need be update things
144 if (jetB != NULL) {
145 double dist = _bj_dist(jetI,jetB);
146 if (dist < jetI->NN_dist) {
147 if (jetI != jetB) {
148 jetI->NN_dist = dist;
149 jetI->NN = jetB;
150 diJ[jetI-head] = _bj_diJ(jetI); // update diJ...
151 }
152 }
153 if (dist < jetB->NN_dist) {
154 if (jetI != jetB) {
155 jetB->NN_dist = dist;
156 jetB->NN = jetI;}
157 }
158 }
159 // if jetI's NN is the new tail then relabel it so that it becomes jetA
160 if (jetI->NN == tail) {jetI->NN = jetA;}
161 }
162
163
164 if (jetB != NULL) {diJ[jetB-head] = _bj_diJ(jetB);}
165
166
167 }
168
169
170 // final cleaning up;
171 delete[] diJ;
172 delete[] briefjets;
173}
174
175
176
177// //----------------------------------------------------------------------
178// // initialises a GenBriefJet
179// template<> inline void ClusterSequence::_bj_set_jetinfo(
180// GenBriefJet * const jetA, const int _jets_index) const {
181//
182// jetA->init(_jets[_jets_index]);
183// jetA->_jets_index = _jets_index;
184//
185// }
186//
187//
188// //----------------------------------------------------------------------
189// // returns the distance between two GenBriefJets
190// template<> double ClusterSequence::_bj_dist(
191// const GenBriefJet * const jeta,
192// const GenBriefJet * const jetb) const {
193// return jeta->geom_ij(jetb);
194// }
195
196FASTJET_END_NAMESPACE
197
198#endif // __FASTJET_CLUSTERQUENCE_N2_ICC__
199
200
Note: See TracBrowser for help on using the repository browser.