Fork me on GitHub

source: git/external/fastjet/contribs/Nsubjettiness/MeasureFunction.cc@ f14bd6a

ImprovedOutputFile Timing dual_readout llp
Last change on this file since f14bd6a was 35cdc46, checked in by Pavel Demin <demin@…>, 10 years ago

upgrade FastJet to version 3.1.0-beta.1, upgrade Nsubjettiness to version 2.1.0, add SoftKiller version 1.0.0

  • Property mode set to 100644
File size: 6.4 KB
Line 
1// Nsubjettiness Package
2// Questions/Comments? jthaler@jthaler.net
3//
4// Copyright (c) 2011-14
5// Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
6//
7// $Id: MeasureFunction.cc 670 2014-06-06 01:24:42Z jthaler $
8//----------------------------------------------------------------------
9// This file is part of FastJet contrib.
10//
11// It is free software; you can redistribute it and/or modify it under
12// the terms of the GNU General Public License as published by the
13// Free Software Foundation; either version 2 of the License, or (at
14// your option) any later version.
15//
16// It is distributed in the hope that it will be useful, but WITHOUT
17// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18// or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
19// License for more details.
20//
21// You should have received a copy of the GNU General Public License
22// along with this code. If not, see <http://www.gnu.org/licenses/>.
23//----------------------------------------------------------------------
24
25
26#include "MeasureFunction.hh"
27
28FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
29
30namespace contrib{
31
32///////
33//
34// Measure Function
35//
36///////
37
38// Return all of the necessary TauComponents for specific input particles and axes
39TauComponents MeasureFunction::result(const std::vector<fastjet::PseudoJet>& particles, const std::vector<fastjet::PseudoJet>& axes) const {
40
41 // first find partition
42 // this sets jetPartitionStorage and beamPartitionStorage
43 PseudoJet beamPartitionStorage;
44 std::vector<fastjet::PseudoJet> jetPartitionStorage = get_partition(particles,axes,&beamPartitionStorage);
45
46 // then return result calculated from partition
47 return result_from_partition(jetPartitionStorage,axes,&beamPartitionStorage);
48}
49
50std::vector<fastjet::PseudoJet> MeasureFunction::get_partition(const std::vector<fastjet::PseudoJet>& particles,
51 const std::vector<fastjet::PseudoJet>& axes,
52 PseudoJet * beamPartitionStorage) const {
53
54 std::vector<std::vector<PseudoJet> > jetPartition(axes.size());
55 std::vector<PseudoJet> beamPartition;
56
57 // Figures out the partiting of the input particles into the various jet pieces
58 // Based on which axis the parition is closest to
59 for (unsigned i = 0; i < particles.size(); i++) {
60
61 // find minimum distance; start with beam (-1) for reference
62 int j_min = -1;
63 double minRsq;
64 if (_has_beam) minRsq = beam_distance_squared(particles[i]);
65 else minRsq = std::numeric_limits<double>::max(); // make it large value
66
67 // check to see which axis the particle is closest to
68 for (unsigned j = 0; j < axes.size(); j++) {
69 double tempRsq = jet_distance_squared(particles[i],axes[j]); // delta R distance
70 if (tempRsq < minRsq) {
71 minRsq = tempRsq;
72 j_min = j;
73 }
74 }
75
76 if (j_min == -1) {
77 if (_has_beam) beamPartition.push_back(particles[i]);
78 else assert(_has_beam); // this should never happen.
79 } else {
80 jetPartition[j_min].push_back(particles[i]);
81 }
82 }
83
84 // Store beam partition
85 if (beamPartitionStorage) {
86 *beamPartitionStorage = join(beamPartition);
87 }
88
89 // Store jet partitions
90 std::vector<PseudoJet> jetPartitionStorage(axes.size(),PseudoJet(0,0,0,0));
91 for (unsigned j = 0; j < axes.size(); j++) {
92 jetPartitionStorage[j] = join(jetPartition[j]);
93 }
94
95 return jetPartitionStorage;
96}
97
98// does partition, but only stores index of PseudoJets
99std::vector<std::list<int> > MeasureFunction::get_partition_list(const std::vector<fastjet::PseudoJet>& particles,
100 const std::vector<fastjet::PseudoJet>& axes) const {
101
102 std::vector<std::list<int> > jetPartition(axes.size());
103
104 // Figures out the partiting of the input particles into the various jet pieces
105 // Based on which axis the parition is closest to
106 for (unsigned i = 0; i < particles.size(); i++) {
107
108 // find minimum distance; start with beam (-1) for reference
109 int j_min = -1;
110 double minRsq;
111 if (_has_beam) minRsq = beam_distance_squared(particles[i]);
112 else minRsq = std::numeric_limits<double>::max(); // make it large value
113
114 // check to see which axis the particle is closest to
115 for (unsigned j = 0; j < axes.size(); j++) {
116 double tempRsq = jet_distance_squared(particles[i],axes[j]); // delta R distance
117 if (tempRsq < minRsq) {
118 minRsq = tempRsq;
119 j_min = j;
120 }
121 }
122
123 if (j_min == -1) {
124 assert(_has_beam); // consistency check
125 } else {
126 jetPartition[j_min].push_back(i);
127 }
128 }
129
130 return jetPartition;
131}
132
133
134// Uses existing partition and calculates result
135// TODO: Can we cache this for speed up when doing area subtraction?
136TauComponents MeasureFunction::result_from_partition(const std::vector<fastjet::PseudoJet>& jet_partition,
137 const std::vector<fastjet::PseudoJet>& axes,
138 PseudoJet * beamPartitionStorage) const {
139
140 std::vector<double> jetPieces(axes.size(), 0.0);
141 double beamPiece = 0.0;
142
143 double tauDen = 0.0;
144 if (!_has_denominator) tauDen = 1.0; // if no denominator, then 1.0 for no normalization factor
145
146 // first find jet pieces
147 for (unsigned j = 0; j < axes.size(); j++) {
148 std::vector<PseudoJet> thisPartition = jet_partition[j].constituents();
149 for (unsigned i = 0; i < thisPartition.size(); i++) {
150 jetPieces[j] += jet_numerator(thisPartition[i],axes[j]); //numerator jet piece
151 if (_has_denominator) tauDen += denominator(thisPartition[i]); // denominator
152 }
153 }
154
155 // then find beam piece
156 if (_has_beam) {
157 assert(beamPartitionStorage); // make sure I have beam information
158 std::vector<PseudoJet> beamPartition = beamPartitionStorage->constituents();
159
160 for (unsigned i = 0; i < beamPartition.size(); i++) {
161 beamPiece += beam_numerator(beamPartition[i]); //numerator beam piece
162 if (_has_denominator) tauDen += denominator(beamPartition[i]); // denominator
163 }
164 }
165 return TauComponents(jetPieces, beamPiece, tauDen, _has_denominator, _has_beam);
166}
167
168
169
170
171
172} //namespace contrib
173
174FASTJET_END_NAMESPACE
Note: See TracBrowser for help on using the repository browser.