Fork me on GitHub

source: git/modules/FastJetFinder.cc@ a3c2d52

Last change on this file since a3c2d52 was 7278220, checked in by Pavel Demin <pavel.demin@…>, 10 years ago

replace map with vector< TEstimatorStruct > and replace SelectorAbsRapRange with SelectorEtaRange

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