]>
Commit | Line | Data |
---|---|---|
370be031 | 1 | //STARTHEADER |
2 | // $Id: ClusterSequenceAreaBase.hh 1503 2009-04-06 11:32:52Z salam $ | |
3 | // | |
4 | // Copyright (c) 2005-2006, Matteo Cacciari and Gavin Salam | |
5 | // | |
6 | //---------------------------------------------------------------------- | |
7 | // This file is part of FastJet. | |
8 | // | |
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. | |
13 | // | |
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. | |
18 | // | |
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. | |
23 | // | |
24 | // You should have received a copy of the GNU General Public License | |
25 | // along with FastJet; if not, write to the Free Software | |
26 | // Foundation, Inc.: | |
27 | // 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA | |
28 | //---------------------------------------------------------------------- | |
29 | //ENDHEADER | |
30 | ||
31 | #ifndef __FASTJET_CLUSTERSEQUENCEAREABASE_HH__ | |
32 | #define __FASTJET_CLUSTERSEQUENCEAREABASE_HH__ | |
33 | ||
34 | #include "fastjet/ClusterSequence.hh" | |
35 | #include "fastjet/internal/LimitedWarning.hh" | |
36 | #include "fastjet/RangeDefinition.hh" | |
37 | ||
38 | namespace fastjet { | |
39 | ||
40 | /// base class that sets interface for extensions of ClusterSequence | |
41 | /// that provide information about the area of each jet; | |
42 | /// | |
43 | /// the virtual functions here all return 0, since no area determination | |
44 | /// is implemented. | |
45 | class ClusterSequenceAreaBase : public ClusterSequence { | |
46 | public: | |
47 | ||
48 | /// a constructor which just carries out the construction of the | |
49 | /// parent class | |
50 | template<class L> ClusterSequenceAreaBase | |
51 | (const std::vector<L> & pseudojets, | |
52 | const JetDefinition & jet_def, | |
53 | const bool & writeout_combinations = false) : | |
54 | ClusterSequence(pseudojets, jet_def, writeout_combinations) {} | |
55 | ||
56 | ||
57 | /// default constructor | |
58 | ClusterSequenceAreaBase() {} | |
59 | ||
60 | ||
61 | /// destructor | |
62 | virtual ~ClusterSequenceAreaBase() {} | |
63 | ||
64 | ||
65 | /// return the area associated with the given jet; this base class | |
66 | /// returns 0. | |
67 | virtual double area (const PseudoJet & jet) const {return 0.0;} | |
68 | ||
69 | /// return the error (uncertainty) associated with the determination | |
70 | /// of the area of this jet; this base class returns 0. | |
71 | virtual double area_error (const PseudoJet & jet) const {return 0.0;} | |
72 | ||
73 | /// return a PseudoJet whose 4-vector is defined by the following integral | |
74 | /// | |
75 | /// \int drap d\phi PseudoJet("rap,phi,pt=one") * | |
76 | /// * Theta("rap,phi inside jet boundary") | |
77 | /// | |
78 | /// where PseudoJet("rap,phi,pt=one") is a 4-vector with the given | |
79 | /// rapidity (rap), azimuth (phi) and pt=1, while Theta("rap,phi | |
80 | /// inside jet boundary") is a function that is 1 when rap,phi | |
81 | /// define a direction inside the jet boundary and 0 otherwise. | |
82 | /// | |
83 | /// This base class returns a null 4-vector. | |
84 | virtual PseudoJet area_4vector(const PseudoJet & jet) const { | |
85 | return PseudoJet(0.0,0.0,0.0,0.0);} | |
86 | ||
87 | /// true if a jet is made exclusively of ghosts | |
88 | /// | |
89 | /// NB: most area classes do not give any explicit ghost jets, but | |
90 | /// some do, and they should replace this function with their own | |
91 | /// version. | |
92 | virtual bool is_pure_ghost(const PseudoJet & jet) const { | |
93 | return false; | |
94 | } | |
95 | ||
96 | /// returns true if ghosts are explicitly included within | |
97 | /// jets for this ClusterSequence; | |
98 | /// | |
99 | /// Derived classes that do include explicit ghosts should provide | |
100 | /// an alternative version of this routine and set it properly. | |
101 | virtual bool has_explicit_ghosts() const { | |
102 | return false; | |
103 | } | |
104 | ||
105 | /// return the total area, within range, that is free of jets, in | |
106 | /// general based on the inclusive jets | |
107 | virtual double empty_area(const RangeDefinition & range) const; | |
108 | ||
109 | /// return the total area, within range, that is free of jets, based | |
110 | /// on the supplied all_jets | |
111 | double empty_area_from_jets(const std::vector<PseudoJet> & all_jets, | |
112 | const RangeDefinition & range) const; | |
113 | ||
114 | /// return something similar to the number of pure ghost jets | |
115 | /// in the given range in an active area case. | |
116 | /// For the local implementation we return empty_area/(0.55 pi R^2), | |
117 | /// based on measured properties of ghost jets with kt and cam. Note | |
118 | /// that the number returned is a double. | |
119 | virtual double n_empty_jets(const RangeDefinition & range) const { | |
120 | double R = jet_def().R(); | |
121 | return empty_area(range)/(0.55*pi*R*R); | |
122 | } | |
123 | ||
124 | /// the median of (pt/area) for jets contained within range, | |
125 | /// making use also of the info on n_empty_jets | |
126 | double median_pt_per_unit_area(const RangeDefinition & range) const; | |
127 | ||
128 | /// the median of (pt/area_4vector) for jets contained within | |
129 | /// making use also of the info on n_empty_jets | |
130 | double median_pt_per_unit_area_4vector(const RangeDefinition & range) const; | |
131 | ||
132 | /// the function that does the work for median_pt_per_unit_area and | |
133 | /// median_pt_per_unit_area_4vector: | |
134 | /// - something_is_area_4vect = false -> use plain area | |
135 | /// - something_is_area_4vect = true -> use 4-vector area | |
136 | double median_pt_per_unit_something( | |
137 | const RangeDefinition & range, bool use_area_4vector) const; | |
138 | ||
139 | /// using jets withing range (and with 4-vector areas if | |
140 | /// use_area_4vector), calculate the median pt/area, as well as an | |
141 | /// "error" (uncertainty), which is defined as the 1-sigma | |
142 | /// half-width of the distribution of pt/A, obtained by looking for | |
143 | /// the point below which we have (1-0.6827)/2 of the jets | |
144 | /// (including empty jets). | |
145 | /// | |
146 | /// The subtraction for a jet with uncorrected pt pt^U and area A is | |
147 | /// | |
148 | /// pt^S = pt^U - median*A +- sigma*sqrt(A) | |
149 | /// | |
150 | /// where the error is only that associated with the fluctuations | |
151 | /// in the noise and not that associated with the noise having | |
152 | /// caused changes in the hard-particle content of the jet. | |
153 | /// | |
154 | /// NB: subtraction may also be done with 4-vector area of course, | |
155 | /// and this is recommended for jets with larger values of R, as | |
156 | /// long as rho has also been determined with a 4-vector area; | |
157 | /// using a scalar area causes one to neglect terms of relative | |
158 | /// order $R^2/8$ in the jet $p_t$. | |
159 | virtual void get_median_rho_and_sigma(const RangeDefinition & range, | |
160 | bool use_area_4vector, | |
161 | double & median, double & sigma, | |
162 | double & mean_area) const; | |
163 | ||
164 | /// a more advanced version of get_median_rho_and_sigma, which allows | |
165 | /// one to use any "view" of the event containing all jets (so that, | |
166 | /// e.g. one might use Cam on a different resolution scale without | |
167 | /// have to rerun the algorithm). | |
168 | /// | |
169 | /// By default it will assume that "all" are not inclusive jets, | |
170 | /// so that in dealing with empty area it has to calculate | |
171 | /// the number of empty jets based on the empty area and the | |
172 | /// the observed <area> of jets rather than a surmised area | |
173 | /// | |
174 | /// Note that for small effective radii, this can cause problems | |
175 | /// because the harder jets get an area >> <ghost-jet-area> | |
176 | /// and so the estimate comes out all wrong. In these situations | |
177 | /// it is highly advisable to use an area with explicit ghosts, since | |
178 | /// then the "empty" jets are actually visible. | |
179 | virtual void get_median_rho_and_sigma(const std::vector<PseudoJet> & all_jets, | |
180 | const RangeDefinition & range, | |
181 | bool use_area_4vector, | |
182 | double & median, double & sigma, | |
183 | double & mean_area, | |
184 | bool all_are_inclusive = false) const; | |
185 | ||
186 | /// same as the full version of get_median_rho_and_error, but without | |
187 | /// access to the mean_area | |
188 | virtual void get_median_rho_and_sigma(const RangeDefinition & range, | |
189 | bool use_area_4vector, | |
190 | double & median, double & sigma) const { | |
191 | double mean_area; | |
192 | get_median_rho_and_sigma(range, use_area_4vector, | |
193 | median, sigma, mean_area); | |
194 | } | |
195 | ||
196 | ||
197 | /// fits a form pt_per_unit_area(y) = a + b*y^2 in the range "range". | |
198 | /// exclude_above allows one to exclude large values of pt/area from fit. | |
199 | /// (if negative, the cut is discarded) | |
200 | /// use_area_4vector = true uses the 4vector areas. | |
201 | virtual void parabolic_pt_per_unit_area(double & a, double & b, | |
202 | const RangeDefinition & range, | |
203 | double exclude_above=-1.0, | |
204 | bool use_area_4vector=false) const; | |
205 | ||
206 | /// return a vector of all subtracted jets, using area_4vector, given rho. | |
207 | /// Only inclusive_jets above ptmin are subtracted and returned. | |
208 | /// the ordering is the same as that of sorted_by_pt(cs.inclusive_jets()), | |
209 | /// i.e. not necessarily ordered in pt once subtracted | |
210 | std::vector<PseudoJet> subtracted_jets(const double rho, | |
211 | const double ptmin=0.0) const; | |
212 | ||
213 | /// return a vector of subtracted jets, using area_4vector. | |
214 | /// Only inclusive_jets above ptmin are subtracted and returned. | |
215 | /// the ordering is the same as that of sorted_by_pt(cs.inclusive_jets()), | |
216 | /// i.e. not necessarily ordered in pt once subtracted | |
217 | std::vector<PseudoJet> subtracted_jets(const RangeDefinition & range, | |
218 | const double ptmin=0.0) const; | |
219 | ||
220 | /// return a subtracted jet, using area_4vector, given rho | |
221 | PseudoJet subtracted_jet(const PseudoJet & jet, | |
222 | const double rho) const; | |
223 | ||
224 | /// return a subtracted jet, using area_4vector; note | |
225 | /// that this is potentially inefficient if repeatedly used for many | |
226 | /// different jets, because rho will be recalculated each time | |
227 | /// around. | |
228 | PseudoJet subtracted_jet(const PseudoJet & jet, | |
229 | const RangeDefinition & range) const; | |
230 | ||
231 | /// return the subtracted pt, given rho | |
232 | double subtracted_pt(const PseudoJet & jet, | |
233 | const double rho, | |
234 | bool use_area_4vector=false) const; | |
235 | ||
236 | /// return the subtracted pt; note that this is | |
237 | /// potentially inefficient if repeatedly used for many different | |
238 | /// jets, because rho will be recalculated each time around. | |
239 | double subtracted_pt(const PseudoJet & jet, | |
240 | const RangeDefinition & range, | |
241 | bool use_area_4vector=false) const; | |
242 | ||
243 | ||
244 | private: | |
245 | /// handle warning messages | |
246 | static LimitedWarning _warnings; | |
247 | ||
248 | /// check the jet algorithm is suitable (and if not issue a warning) | |
249 | void _check_jet_alg_good_for_median() const; | |
250 | ||
251 | }; | |
252 | ||
253 | ||
254 | ||
255 | } // fastjet namespace | |
256 | ||
257 | #endif // __FASTJET_CLUSTERSEQUENCEAREABASE_HH__ |