1 |
|
---|
2 | /** \class FastJetFinder
|
---|
3 | *
|
---|
4 | * Finds jets using FastJet library.
|
---|
5 | *
|
---|
6 | * $Date: 2014-04-17 10:28:09 +0000 (Thu, 17 Apr 2014) $
|
---|
7 | * $Revision: 1383 $
|
---|
8 | *
|
---|
9 | *
|
---|
10 | * \author P. Demin - UCL, Louvain-la-Neuve
|
---|
11 | *
|
---|
12 | */
|
---|
13 |
|
---|
14 | #include "modules/FastJetFinder.h"
|
---|
15 |
|
---|
16 | #include "classes/DelphesClasses.h"
|
---|
17 | #include "classes/DelphesFactory.h"
|
---|
18 | #include "classes/DelphesFormula.h"
|
---|
19 |
|
---|
20 | #include "ExRootAnalysis/ExRootResult.h"
|
---|
21 | #include "ExRootAnalysis/ExRootFilter.h"
|
---|
22 | #include "ExRootAnalysis/ExRootClassifier.h"
|
---|
23 |
|
---|
24 | #include "TMath.h"
|
---|
25 | #include "TString.h"
|
---|
26 | #include "TFormula.h"
|
---|
27 | #include "TRandom3.h"
|
---|
28 | #include "TObjArray.h"
|
---|
29 | #include "TDatabasePDG.h"
|
---|
30 | #include "TLorentzVector.h"
|
---|
31 |
|
---|
32 | #include <algorithm>
|
---|
33 | #include <stdexcept>
|
---|
34 | #include <iostream>
|
---|
35 | #include <sstream>
|
---|
36 | #include <vector>
|
---|
37 |
|
---|
38 | #include "fastjet/PseudoJet.hh"
|
---|
39 | #include "fastjet/JetDefinition.hh"
|
---|
40 | #include "fastjet/ClusterSequence.hh"
|
---|
41 | #include "fastjet/Selector.hh"
|
---|
42 | #include "fastjet/ClusterSequenceArea.hh"
|
---|
43 | #include "fastjet/tools/JetMedianBackgroundEstimator.hh"
|
---|
44 |
|
---|
45 | #include "fastjet/plugins/SISCone/fastjet/SISConePlugin.hh"
|
---|
46 | #include "fastjet/plugins/CDFCones/fastjet/CDFMidPointPlugin.hh"
|
---|
47 | #include "fastjet/plugins/CDFCones/fastjet/CDFJetCluPlugin.hh"
|
---|
48 |
|
---|
49 | #include "fastjet/contribs/Nsubjettiness/Nsubjettiness.hh"
|
---|
50 | #include "fastjet/contribs/Nsubjettiness/Njettiness.hh"
|
---|
51 | #include "fastjet/contribs/Nsubjettiness/NjettinessPlugin.hh"
|
---|
52 | #include "fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.hh"
|
---|
53 |
|
---|
54 | using namespace std;
|
---|
55 | using namespace fastjet;
|
---|
56 | using namespace fastjet::contrib;
|
---|
57 |
|
---|
58 |
|
---|
59 | //------------------------------------------------------------------------------
|
---|
60 |
|
---|
61 | FastJetFinder::FastJetFinder() :
|
---|
62 | fPlugin(0), fRecomb(0), fNjettinessPlugin(0), fDefinition(0), fAreaDefinition(0), fItInputArray(0)
|
---|
63 | {
|
---|
64 |
|
---|
65 | }
|
---|
66 |
|
---|
67 | //------------------------------------------------------------------------------
|
---|
68 |
|
---|
69 | FastJetFinder::~FastJetFinder()
|
---|
70 | {
|
---|
71 |
|
---|
72 | }
|
---|
73 |
|
---|
74 | //------------------------------------------------------------------------------
|
---|
75 |
|
---|
76 | void FastJetFinder::Init()
|
---|
77 | {
|
---|
78 | JetDefinition::Plugin *plugin = 0;
|
---|
79 | JetDefinition::Recombiner *recomb = 0;
|
---|
80 | NjettinessPlugin *njetPlugin = 0;
|
---|
81 |
|
---|
82 | // read eta ranges
|
---|
83 |
|
---|
84 | ExRootConfParam param = GetParam("RhoEtaRange");
|
---|
85 | Long_t i, size;
|
---|
86 |
|
---|
87 | fEtaRangeMap.clear();
|
---|
88 | size = param.GetSize();
|
---|
89 | for(i = 0; i < size/2; ++i)
|
---|
90 | {
|
---|
91 | fEtaRangeMap[param[i*2].GetDouble()] = param[i*2 + 1].GetDouble();
|
---|
92 | }
|
---|
93 |
|
---|
94 | // define algorithm
|
---|
95 |
|
---|
96 | fJetAlgorithm = GetInt("JetAlgorithm", 6);
|
---|
97 | fParameterR = GetDouble("ParameterR", 0.5);
|
---|
98 |
|
---|
99 | fConeRadius = GetDouble("ConeRadius", 0.5);
|
---|
100 | fSeedThreshold = GetDouble("SeedThreshold", 1.0);
|
---|
101 | fConeAreaFraction = GetDouble("ConeAreaFraction", 1.0);
|
---|
102 | fMaxIterations = GetInt("MaxIterations", 100);
|
---|
103 | fMaxPairSize = GetInt("MaxPairSize", 2);
|
---|
104 | fIratch = GetInt("Iratch", 1);
|
---|
105 | fAdjacencyCut = GetInt("AdjacencyCut", 2);
|
---|
106 | fOverlapThreshold = GetDouble("OverlapThreshold", 0.75);
|
---|
107 |
|
---|
108 | fJetPTMin = GetDouble("JetPTMin", 10.0);
|
---|
109 |
|
---|
110 | //-- N(sub)jettiness parameters --
|
---|
111 |
|
---|
112 | fComputeNsubjettiness = GetBool("ComputeNsubjettiness", false);
|
---|
113 | fBeta = GetDouble("Beta", 1.0);
|
---|
114 | fAxisMode = GetInt("AxisMode", 1);
|
---|
115 | fRcutOff = GetDouble("RcutOff", 0.8); // used only if Njettiness is used as jet clustering algo (case 8)
|
---|
116 | fN = GetInt("N", 2); // used only if Njettiness is used as jet clustering algo (case 8)
|
---|
117 |
|
---|
118 | // --- Jet Area Parameters ---
|
---|
119 | fAreaAlgorithm = GetInt("AreaAlgorithm", 0);
|
---|
120 | fComputeRho = GetBool("ComputeRho", false);
|
---|
121 |
|
---|
122 | // - ghost based areas -
|
---|
123 | fGhostEtaMax = GetDouble("GhostEtaMax", 5.0);
|
---|
124 | fRepeat = GetInt("Repeat", 1);
|
---|
125 | fGhostArea = GetDouble("GhostArea", 0.01);
|
---|
126 | fGridScatter = GetDouble("GridScatter", 1.0);
|
---|
127 | fPtScatter = GetDouble("PtScatter", 0.1);
|
---|
128 | fMeanGhostPt = GetDouble("MeanGhostPt", 1.0E-100);
|
---|
129 |
|
---|
130 | // - voronoi based areas -
|
---|
131 | fEffectiveRfact = GetDouble("EffectiveRfact", 1.0);
|
---|
132 |
|
---|
133 | switch(fAreaAlgorithm)
|
---|
134 | {
|
---|
135 | case 1:
|
---|
136 | fAreaDefinition = new AreaDefinition(active_area_explicit_ghosts, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
|
---|
137 | break;
|
---|
138 | case 2:
|
---|
139 | fAreaDefinition = new AreaDefinition(one_ghost_passive_area, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
|
---|
140 | break;
|
---|
141 | case 3:
|
---|
142 | fAreaDefinition = new AreaDefinition(passive_area, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
|
---|
143 | break;
|
---|
144 | case 4:
|
---|
145 | fAreaDefinition = new AreaDefinition(VoronoiAreaSpec(fEffectiveRfact));
|
---|
146 | break;
|
---|
147 | case 5:
|
---|
148 | fAreaDefinition = new AreaDefinition(active_area, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
|
---|
149 | break;
|
---|
150 | default:
|
---|
151 | case 0:
|
---|
152 | fAreaDefinition = 0;
|
---|
153 | break;
|
---|
154 | }
|
---|
155 |
|
---|
156 | switch(fJetAlgorithm)
|
---|
157 | {
|
---|
158 | case 1:
|
---|
159 | plugin = new CDFJetCluPlugin(fSeedThreshold, fConeRadius, fAdjacencyCut, fMaxIterations, fIratch, fOverlapThreshold);
|
---|
160 | fDefinition = new JetDefinition(plugin);
|
---|
161 | break;
|
---|
162 | case 2:
|
---|
163 | plugin = new CDFMidPointPlugin(fSeedThreshold, fConeRadius, fConeAreaFraction, fMaxPairSize, fMaxIterations, fOverlapThreshold);
|
---|
164 | fDefinition = new JetDefinition(plugin);
|
---|
165 | break;
|
---|
166 | case 3:
|
---|
167 | plugin = new SISConePlugin(fConeRadius, fOverlapThreshold, fMaxIterations, fJetPTMin);
|
---|
168 | fDefinition = new JetDefinition(plugin);
|
---|
169 | break;
|
---|
170 | case 4:
|
---|
171 | fDefinition = new JetDefinition(kt_algorithm, fParameterR);
|
---|
172 | break;
|
---|
173 | case 5:
|
---|
174 | fDefinition = new JetDefinition(cambridge_algorithm, fParameterR);
|
---|
175 | break;
|
---|
176 | default:
|
---|
177 | case 6:
|
---|
178 | fDefinition = new JetDefinition(antikt_algorithm, fParameterR);
|
---|
179 | break;
|
---|
180 | case 7:
|
---|
181 | recomb = new WinnerTakeAllRecombiner();
|
---|
182 | fDefinition = new JetDefinition(antikt_algorithm, fParameterR, recomb, Best);
|
---|
183 | break;
|
---|
184 | case 8:
|
---|
185 | njetPlugin = new NjettinessPlugin(fN, Njettiness::wta_kt_axes, Njettiness::unnormalized_cutoff_measure, fBeta, fRcutOff);
|
---|
186 | fDefinition = new JetDefinition(njetPlugin);
|
---|
187 | break;
|
---|
188 | }
|
---|
189 |
|
---|
190 | fPlugin = plugin;
|
---|
191 | fRecomb = recomb;
|
---|
192 | fNjettinessPlugin = njetPlugin;
|
---|
193 |
|
---|
194 | ClusterSequence::print_banner();
|
---|
195 |
|
---|
196 | // import input array
|
---|
197 |
|
---|
198 | fInputArray = ImportArray(GetString("InputArray", "Calorimeter/towers"));
|
---|
199 | fItInputArray = fInputArray->MakeIterator();
|
---|
200 |
|
---|
201 | // create output arrays
|
---|
202 |
|
---|
203 | fOutputArray = ExportArray(GetString("OutputArray", "jets"));
|
---|
204 | fRhoOutputArray = ExportArray(GetString("RhoOutputArray", "rho"));
|
---|
205 | }
|
---|
206 |
|
---|
207 | //------------------------------------------------------------------------------
|
---|
208 |
|
---|
209 | void FastJetFinder::Finish()
|
---|
210 | {
|
---|
211 | if(fItInputArray) delete fItInputArray;
|
---|
212 | if(fDefinition) delete fDefinition;
|
---|
213 | if(fAreaDefinition) delete fAreaDefinition;
|
---|
214 | if(fPlugin) delete static_cast<JetDefinition::Plugin*>(fPlugin);
|
---|
215 | if(fRecomb) delete static_cast<JetDefinition::Recombiner*>(fRecomb);
|
---|
216 | if(fNjettinessPlugin) delete static_cast<JetDefinition::Plugin*>(fNjettinessPlugin);
|
---|
217 | }
|
---|
218 |
|
---|
219 | //------------------------------------------------------------------------------
|
---|
220 |
|
---|
221 | void FastJetFinder::Process()
|
---|
222 | {
|
---|
223 | Candidate *candidate, *constituent;
|
---|
224 | TLorentzVector momentum;
|
---|
225 | Double_t deta, dphi, detaMax, dphiMax;
|
---|
226 | Double_t time, weightTime, avTime;
|
---|
227 | Int_t number;
|
---|
228 | Double_t rho = 0;
|
---|
229 | PseudoJet jet, area;
|
---|
230 | vector<PseudoJet> inputList, outputList;
|
---|
231 | ClusterSequence *sequence;
|
---|
232 | map< Double_t, Double_t >::iterator itEtaRangeMap;
|
---|
233 |
|
---|
234 | DelphesFactory *factory = GetFactory();
|
---|
235 |
|
---|
236 | inputList.clear();
|
---|
237 |
|
---|
238 | // loop over input objects
|
---|
239 | fItInputArray->Reset();
|
---|
240 | number = 0;
|
---|
241 | while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
|
---|
242 | {
|
---|
243 | momentum = candidate->Momentum;
|
---|
244 | jet = PseudoJet(momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E());
|
---|
245 | jet.set_user_index(number);
|
---|
246 | inputList.push_back(jet);
|
---|
247 | ++number;
|
---|
248 | }
|
---|
249 |
|
---|
250 | // construct jets
|
---|
251 | if(fAreaDefinition)
|
---|
252 | {
|
---|
253 | sequence = new ClusterSequenceArea(inputList, *fDefinition, *fAreaDefinition);
|
---|
254 | }
|
---|
255 | else
|
---|
256 | {
|
---|
257 | sequence = new ClusterSequence(inputList, *fDefinition);
|
---|
258 | }
|
---|
259 |
|
---|
260 | // compute rho and store it
|
---|
261 | if(fComputeRho && fAreaDefinition)
|
---|
262 | {
|
---|
263 | for(itEtaRangeMap = fEtaRangeMap.begin(); itEtaRangeMap != fEtaRangeMap.end(); ++itEtaRangeMap)
|
---|
264 | {
|
---|
265 | Selector select_rapidity = SelectorAbsRapRange(itEtaRangeMap->first, itEtaRangeMap->second);
|
---|
266 | JetMedianBackgroundEstimator estimator(select_rapidity, *fDefinition, *fAreaDefinition);
|
---|
267 | estimator.set_particles(inputList);
|
---|
268 | rho = estimator.rho();
|
---|
269 |
|
---|
270 | candidate = factory->NewCandidate();
|
---|
271 | candidate->Momentum.SetPtEtaPhiE(rho, 0.0, 0.0, rho);
|
---|
272 | candidate->Edges[0] = itEtaRangeMap->first;
|
---|
273 | candidate->Edges[1] = itEtaRangeMap->second;
|
---|
274 | fRhoOutputArray->Add(candidate);
|
---|
275 | }
|
---|
276 | }
|
---|
277 |
|
---|
278 | outputList.clear();
|
---|
279 | outputList = sorted_by_pt(sequence->inclusive_jets(fJetPTMin));
|
---|
280 |
|
---|
281 |
|
---|
282 | // loop over all jets and export them
|
---|
283 | detaMax = 0.0;
|
---|
284 | dphiMax = 0.0;
|
---|
285 | vector<PseudoJet>::iterator itInputList, itOutputList;
|
---|
286 | for(itOutputList = outputList.begin(); itOutputList != outputList.end(); ++itOutputList)
|
---|
287 | {
|
---|
288 | momentum.SetPxPyPzE(itOutputList->px(), itOutputList->py(), itOutputList->pz(), itOutputList->E());
|
---|
289 | area.reset(0.0, 0.0, 0.0, 0.0);
|
---|
290 | if(fAreaDefinition) area = itOutputList->area_4vector();
|
---|
291 |
|
---|
292 | candidate = factory->NewCandidate();
|
---|
293 |
|
---|
294 | time=0;
|
---|
295 | weightTime=0;
|
---|
296 |
|
---|
297 | inputList.clear();
|
---|
298 | inputList = sequence->constituents(*itOutputList);
|
---|
299 |
|
---|
300 | for(itInputList = inputList.begin(); itInputList != inputList.end(); ++itInputList)
|
---|
301 | {
|
---|
302 | constituent = static_cast<Candidate*>(fInputArray->At(itInputList->user_index()));
|
---|
303 |
|
---|
304 | deta = TMath::Abs(momentum.Eta() - constituent->Momentum.Eta());
|
---|
305 | dphi = TMath::Abs(momentum.DeltaPhi(constituent->Momentum));
|
---|
306 | if(deta > detaMax) detaMax = deta;
|
---|
307 | if(dphi > dphiMax) dphiMax = dphi;
|
---|
308 |
|
---|
309 | time += TMath::Sqrt(constituent->Momentum.E())*(constituent->Position.T());
|
---|
310 | weightTime += TMath::Sqrt(constituent->Momentum.E());
|
---|
311 |
|
---|
312 | candidate->AddCandidate(constituent);
|
---|
313 | }
|
---|
314 |
|
---|
315 | avTime = time/weightTime;
|
---|
316 |
|
---|
317 | candidate->Momentum = momentum;
|
---|
318 | candidate->Position.SetT(avTime);
|
---|
319 | candidate->Area.SetPxPyPzE(area.px(), area.py(), area.pz(), area.E());
|
---|
320 |
|
---|
321 | candidate->DeltaEta = detaMax;
|
---|
322 | candidate->DeltaPhi = dphiMax;
|
---|
323 |
|
---|
324 | // --- compute N-subjettiness with N = 1,2,3,4,5 ----
|
---|
325 |
|
---|
326 | if(fComputeNsubjettiness)
|
---|
327 | {
|
---|
328 | Njettiness::AxesMode axisMode;
|
---|
329 |
|
---|
330 | switch(fAxisMode)
|
---|
331 | {
|
---|
332 | default:
|
---|
333 | case 1:
|
---|
334 | axisMode = Njettiness::wta_kt_axes;
|
---|
335 | break;
|
---|
336 | case 2:
|
---|
337 | axisMode = Njettiness::onepass_wta_kt_axes;
|
---|
338 | break;
|
---|
339 | case 3:
|
---|
340 | axisMode = Njettiness::kt_axes;
|
---|
341 | break;
|
---|
342 | case 4:
|
---|
343 | axisMode = Njettiness::onepass_kt_axes;
|
---|
344 | break;
|
---|
345 | }
|
---|
346 |
|
---|
347 | Njettiness::MeasureMode measureMode = Njettiness::unnormalized_measure;
|
---|
348 |
|
---|
349 | Nsubjettiness nSub1(1, axisMode, measureMode, fBeta);
|
---|
350 | Nsubjettiness nSub2(2, axisMode, measureMode, fBeta);
|
---|
351 | Nsubjettiness nSub3(3, axisMode, measureMode, fBeta);
|
---|
352 | Nsubjettiness nSub4(4, axisMode, measureMode, fBeta);
|
---|
353 | Nsubjettiness nSub5(5, axisMode, measureMode, fBeta);
|
---|
354 |
|
---|
355 | candidate->Tau[0] = nSub1(*itOutputList);
|
---|
356 | candidate->Tau[1] = nSub2(*itOutputList);
|
---|
357 | candidate->Tau[2] = nSub3(*itOutputList);
|
---|
358 | candidate->Tau[3] = nSub4(*itOutputList);
|
---|
359 | candidate->Tau[4] = nSub5(*itOutputList);
|
---|
360 | }
|
---|
361 |
|
---|
362 |
|
---|
363 | fOutputArray->Add(candidate);
|
---|
364 | }
|
---|
365 | delete sequence;
|
---|
366 | }
|
---|