]>
Commit | Line | Data |
---|---|---|
370be031 | 1 | //STARTHEADER |
2 | // $Id: ClusterSequenceActiveArea.hh 1134 2008-03-15 17:05:16Z 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_CLUSTERSEQUENCEACTIVEAREA_HH__ | |
32 | #define __FASTJET_CLUSTERSEQUENCEACTIVEAREA_HH__ | |
33 | ||
34 | ||
35 | #include "fastjet/PseudoJet.hh" | |
36 | #include "fastjet/ClusterSequenceAreaBase.hh" | |
37 | #include "fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh" | |
38 | #include<iostream> | |
39 | #include<vector> | |
40 | ||
41 | //------------ backwards compatibility with version 2.1 ------------- | |
42 | // for backwards compatibility make ActiveAreaSpec name available | |
43 | #include "fastjet/ActiveAreaSpec.hh" | |
44 | #include "fastjet/ClusterSequenceWithArea.hh" | |
45 | //-------------------------------------------------------------------- | |
46 | ||
47 | ||
48 | namespace fastjet { // defined in fastjet/internal/base.hh | |
49 | ||
50 | //using namespace std; | |
51 | ||
52 | /// Class that behaves essentially like ClusterSequence except | |
53 | /// that it also provides access to the area of a jet (which | |
54 | /// will be a random quantity... Figure out what to do about seeds | |
55 | /// later...) | |
56 | class ClusterSequenceActiveArea : public ClusterSequenceAreaBase { | |
57 | public: | |
58 | ||
59 | /// default constructor | |
60 | ClusterSequenceActiveArea() {} | |
61 | ||
62 | /// constructor based on JetDefinition and GhostedAreaSpec | |
63 | template<class L> ClusterSequenceActiveArea | |
64 | (const std::vector<L> & pseudojets, | |
65 | const JetDefinition & jet_def, | |
66 | const GhostedAreaSpec & area_spec, | |
67 | const bool & writeout_combinations = false) ; | |
68 | ||
69 | virtual double area (const PseudoJet & jet) const { | |
70 | return _average_area[jet.cluster_hist_index()];}; | |
71 | virtual double area_error (const PseudoJet & jet) const { | |
72 | return _average_area2[jet.cluster_hist_index()];}; | |
73 | ||
74 | virtual PseudoJet area_4vector (const PseudoJet & jet) const { | |
75 | return _average_area_4vector[jet.cluster_hist_index()];}; | |
76 | ||
77 | /// enum providing a variety of tentative strategies for estimating | |
78 | /// the background (e.g. non-jet) activity in a highly populated event; the | |
79 | /// one that has been most extensively tested is median. | |
80 | enum mean_pt_strategies{median=0, non_ghost_median, pttot_over_areatot, | |
81 | pttot_over_areatot_cut, mean_ratio_cut, play, | |
82 | median_4vector}; | |
83 | ||
84 | /// return the transverse momentum per unit area according to one | |
85 | /// of the above strategies; for some strategies (those with "cut" | |
86 | /// in their name) the parameter "range" allows one to exclude a | |
87 | /// subset of the jets for the background estimation, those that | |
88 | /// have pt/area > median(pt/area)*range. | |
89 | /// | |
90 | /// NB: This call is OBSOLETE; use media_pt_per_unit_area from the | |
91 | // ClusterSequenceAreaBase class instead | |
92 | double pt_per_unit_area(mean_pt_strategies strat=median, | |
93 | double range=2.0 ) const; | |
94 | ||
95 | // following code removed -- now dealt with by AreaBase class (and | |
96 | // this definition here conflicts with it). | |
97 | // /// fits a form pt_per_unit_area(y) = a + b*y^2 in the range | |
98 | // /// abs(y)<raprange (for negative raprange, it defaults to | |
99 | // /// _safe_rap_for_area). | |
100 | // void parabolic_pt_per_unit_area(double & a,double & b, double raprange=-1.0, | |
101 | // double exclude_above=-1.0, | |
102 | // bool use_area_4vector=false ) const; | |
103 | // | |
104 | /// rewrite the empty area from the parent class, so as to use | |
105 | /// all info at our disposal | |
106 | /// return the total area, in the given y-phi range, that consists of ghost | |
107 | /// jets or unclustered ghosts | |
108 | virtual double empty_area(const RangeDefinition & range) const; | |
109 | ||
110 | /// return the true number of empty jets (replaces | |
111 | /// ClusterSequenceAreaBase::n_empty_jets(...)) | |
112 | virtual double n_empty_jets(const RangeDefinition & range) const; | |
113 | ||
114 | protected: | |
115 | void _resize_and_zero_AA (); | |
116 | void _initialise_AA(const JetDefinition & jet_def, | |
117 | const GhostedAreaSpec & area_spec, | |
118 | const bool & writeout_combinations, | |
119 | bool & continue_running); | |
120 | ||
121 | void _run_AA(const GhostedAreaSpec & area_spec); | |
122 | ||
123 | void _postprocess_AA(const GhostedAreaSpec & area_spec); | |
124 | ||
125 | /// does the initialisation and running specific to the active | |
126 | /// areas class | |
127 | void _initialise_and_run_AA (const JetDefinition & jet_def, | |
128 | const GhostedAreaSpec & area_spec, | |
129 | const bool & writeout_combinations = false); | |
130 | ||
131 | /// transfer the history (and jet-momenta) from clust_seq to our | |
132 | /// own internal structure while removing ghosts | |
133 | void _transfer_ghost_free_history( | |
134 | const ClusterSequenceActiveAreaExplicitGhosts & clust_seq); | |
135 | ||
136 | ||
137 | /// transfer areas from the ClusterSequenceActiveAreaExplicitGhosts | |
138 | /// object into our internal area bookkeeping... | |
139 | void _transfer_areas(const std::vector<int> & unique_hist_order, | |
140 | const ClusterSequenceActiveAreaExplicitGhosts & ); | |
141 | ||
142 | /// child classes benefit from having these at their disposal | |
143 | std::valarray<double> _average_area, _average_area2; | |
144 | std::valarray<PseudoJet> _average_area_4vector; | |
145 | ||
146 | /// returns true if there are any particles whose transverse momentum | |
147 | /// if so low that there's a risk of the ghosts having modified the | |
148 | /// clustering sequence | |
149 | bool has_dangerous_particles() const {return _has_dangerous_particles;} | |
150 | ||
151 | private: | |
152 | ||
153 | ||
154 | double _non_jet_area, _non_jet_area2, _non_jet_number; | |
155 | ||
156 | double _maxrap_for_area; // max rap where we put ghosts | |
157 | double _safe_rap_for_area; // max rap where we trust jet areas | |
158 | ||
159 | bool _has_dangerous_particles; | |
160 | ||
161 | ||
162 | /// routine for extracting the tree in an order that will be independent | |
163 | /// of any degeneracies in the recombination sequence that don't | |
164 | /// affect the composition of the final jets | |
165 | void _extract_tree(std::vector<int> &) const; | |
166 | /// do the part of the extraction associated with pos, working | |
167 | /// through its children and their parents | |
168 | void _extract_tree_children(int pos, std::valarray<bool> &, const std::valarray<int> &, std::vector<int> &) const; | |
169 | /// do the part of the extraction associated with the parents of pos. | |
170 | void _extract_tree_parents (int pos, std::valarray<bool> &, const std::valarray<int> &, std::vector<int> &) const; | |
171 | ||
172 | /// check if two jets have the same momentum to within the | |
173 | /// tolerance (and if pt's are not the same we're forgiving and | |
174 | /// look to see if the energy is the same) | |
175 | void _throw_unless_jets_have_same_perp_or_E(const PseudoJet & jet, | |
176 | const PseudoJet & refjet, | |
177 | double tolerance, | |
178 | const ClusterSequenceActiveAreaExplicitGhosts & jets_ghosted_seq | |
179 | ) const; | |
180 | ||
181 | /// since we are playing nasty games with seeds, we should warn | |
182 | /// the user a few times | |
183 | //static int _n_seed_warnings; | |
184 | //const static int _max_seed_warnings = 10; | |
185 | ||
186 | // record the number of repeats | |
187 | int _area_spec_repeat; | |
188 | ||
189 | /// a class for our internal storage of ghost jets | |
190 | class GhostJet : public PseudoJet { | |
191 | public: | |
192 | GhostJet(const PseudoJet & j, double a) : PseudoJet(j), area(a){} | |
193 | double area; | |
194 | }; | |
195 | ||
196 | std::vector<GhostJet> _ghost_jets; | |
197 | std::vector<GhostJet> _unclustered_ghosts; | |
198 | }; | |
199 | ||
200 | ||
201 | ||
202 | ||
203 | template<class L> ClusterSequenceActiveArea::ClusterSequenceActiveArea | |
204 | (const std::vector<L> & pseudojets, | |
205 | const JetDefinition & jet_def, | |
206 | const GhostedAreaSpec & area_spec, | |
207 | const bool & writeout_combinations) { | |
208 | ||
209 | // transfer the initial jets (type L) into our own array | |
210 | _transfer_input_jets(pseudojets); | |
211 | ||
212 | // run the clustering for active areas | |
213 | _initialise_and_run_AA(jet_def, area_spec, writeout_combinations); | |
214 | ||
215 | } | |
216 | ||
217 | ||
218 | ||
219 | } // fastjet namespace | |
220 | ||
221 | #endif // __FASTJET_CLUSTERSEQUENCEACTIVEAREA_HH__ |