Fork me on GitHub

source: git/modules/FastJetFinder.cc@ e5ea42e

Last change on this file since e5ea42e was 58dcd2e, checked in by Michele Selvaggi <michele.selvaggi@…>, 9 years ago

ensure bwd compatibility (for real this time?)
  • Property mode set to 100644
File size: 16.0 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/ExtraRecombiners.hh"
67
68#include "fastjet/tools/Filter.hh"
69#include "fastjet/tools/Pruner.hh"
70#include "fastjet/contribs/RecursiveTools/SoftDrop.hh"
71
72using namespace std;
73using namespace fastjet;
74using namespace fastjet::contrib;
75
76
77//------------------------------------------------------------------------------
78
79FastJetFinder::FastJetFinder() :
80 fPlugin(0), fRecomb(0), fAxesDef(0), fMeasureDef(0), fNjettinessPlugin(0),
81 fDefinition(0), fAreaDefinition(0), fItInputArray(0)
82{
83
84}
85
86//------------------------------------------------------------------------------
87
88FastJetFinder::~FastJetFinder()
89{
90
91}
92
93//------------------------------------------------------------------------------
94
95void FastJetFinder::Init()
96{
97 JetDefinition::Plugin *plugin = 0;
98 JetDefinition::Recombiner *recomb = 0;
99 ExRootConfParam param;
100 Long_t i, size;
101 Double_t etaMin, etaMax;
102 TEstimatorStruct estimatorStruct;
103
104 // define algorithm
105
106 fJetAlgorithm = GetInt("JetAlgorithm", 6);
107 fParameterR = GetDouble("ParameterR", 0.5);
108
109 fConeRadius = GetDouble("ConeRadius", 0.5);
110 fSeedThreshold = GetDouble("SeedThreshold", 1.0);
111 fConeAreaFraction = GetDouble("ConeAreaFraction", 1.0);
112 fMaxIterations = GetInt("MaxIterations", 100);
113 fMaxPairSize = GetInt("MaxPairSize", 2);
114 fIratch = GetInt("Iratch", 1);
115 fAdjacencyCut = GetInt("AdjacencyCut", 2);
116 fOverlapThreshold = GetDouble("OverlapThreshold", 0.75);
117
118 fJetPTMin = GetDouble("JetPTMin", 10.0);
119
120 //-- N(sub)jettiness parameters --
121
122 fComputeNsubjettiness = GetBool("ComputeNsubjettiness", false);
123 fBeta = GetDouble("Beta", 1.0);
124 fAxisMode = GetInt("AxisMode", 1);
125 fRcutOff = GetDouble("RcutOff", 0.8); // used only if Njettiness is used as jet clustering algo (case 8)
126 fN = GetInt("N", 2); // used only if Njettiness is used as jet clustering algo (case 8)
127
128 fMeasureDef = new UnnormalizedMeasure(fBeta);
129
130 switch(fAxisMode)
131 {
132 default:
133 case 1:
134 fAxesDef = new WTA_KT_Axes();
135 break;
136 case 2:
137 fAxesDef = new OnePass_WTA_KT_Axes();
138 break;
139 case 3:
140 fAxesDef = new KT_Axes();
141 break;
142 case 4:
143 fAxesDef = new OnePass_KT_Axes();
144 }
145
146 //-- Trimming parameters --
147
148 fComputeTrimming = GetBool("ComputeTrimming", false);
149 fRTrim = GetDouble("RTrim", 0.2);
150 fPtFracTrim = GetDouble("PtFracTrim", 0.05);
151
152
153 //-- Pruning parameters --
154
155 fComputePruning = GetBool("ComputePruning", false);
156 fZcutPrun = GetDouble("ZcutPrun", 0.1);
157 fRcutPrun = GetDouble("RcutPrun", 0.5);
158 fRPrun = GetDouble("RPrun", 0.8);
159
160 //-- SoftDrop parameters --
161
162 fComputeSoftDrop = GetBool("ComputeSoftDrop", false);
163 fBetaSoftDrop = GetDouble("BetaSoftDrop", 0.0);
164 fSymmetryCutSoftDrop = GetDouble("SymmetryCutSoftDrop", 0.1);
165 fR0SoftDrop= GetDouble("R0SoftDrop=", 0.8);
166
167
168 // --- Jet Area Parameters ---
169 fAreaAlgorithm = GetInt("AreaAlgorithm", 0);
170 fComputeRho = GetBool("ComputeRho", false);
171
172 // - ghost based areas -
173 fGhostEtaMax = GetDouble("GhostEtaMax", 5.0);
174 fRepeat = GetInt("Repeat", 1);
175 fGhostArea = GetDouble("GhostArea", 0.01);
176 fGridScatter = GetDouble("GridScatter", 1.0);
177 fPtScatter = GetDouble("PtScatter", 0.1);
178 fMeanGhostPt = GetDouble("MeanGhostPt", 1.0E-100);
179
180 // - voronoi based areas -
181 fEffectiveRfact = GetDouble("EffectiveRfact", 1.0);
182
183 switch(fAreaAlgorithm)
184 {
185 case 1:
186 fAreaDefinition = new AreaDefinition(active_area_explicit_ghosts, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
187 break;
188 case 2:
189 fAreaDefinition = new AreaDefinition(one_ghost_passive_area, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
190 break;
191 case 3:
192 fAreaDefinition = new AreaDefinition(passive_area, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
193 break;
194 case 4:
195 fAreaDefinition = new AreaDefinition(VoronoiAreaSpec(fEffectiveRfact));
196 break;
197 case 5:
198 fAreaDefinition = new AreaDefinition(active_area, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
199 break;
200 default:
201 case 0:
202 fAreaDefinition = 0;
203 break;
204 }
205
206 switch(fJetAlgorithm)
207 {
208 case 1:
209 plugin = new CDFJetCluPlugin(fSeedThreshold, fConeRadius, fAdjacencyCut, fMaxIterations, fIratch, fOverlapThreshold);
210 fDefinition = new JetDefinition(plugin);
211 break;
212 case 2:
213 plugin = new CDFMidPointPlugin(fSeedThreshold, fConeRadius, fConeAreaFraction, fMaxPairSize, fMaxIterations, fOverlapThreshold);
214 fDefinition = new JetDefinition(plugin);
215 break;
216 case 3:
217 plugin = new SISConePlugin(fConeRadius, fOverlapThreshold, fMaxIterations, fJetPTMin);
218 fDefinition = new JetDefinition(plugin);
219 break;
220 case 4:
221 fDefinition = new JetDefinition(kt_algorithm, fParameterR);
222 break;
223 case 5:
224 fDefinition = new JetDefinition(cambridge_algorithm, fParameterR);
225 break;
226 default:
227 case 6:
228 fDefinition = new JetDefinition(antikt_algorithm, fParameterR);
229 break;
230 case 7:
231 recomb = new WinnerTakeAllRecombiner();
232 fDefinition = new JetDefinition(antikt_algorithm, fParameterR, recomb, Best);
233 break;
234 case 8:
235 fNjettinessPlugin = new NjettinessPlugin(fN, Njettiness::wta_kt_axes, Njettiness::unnormalized_cutoff_measure, fBeta, fRcutOff);
236 fDefinition = new JetDefinition(fNjettinessPlugin);
237 break;
238 }
239
240
241
242 fPlugin = plugin;
243 fRecomb = recomb;
244
245 ClusterSequence::print_banner();
246
247 if(fComputeRho && fAreaDefinition)
248 {
249 // read eta ranges
250
251 param = GetParam("RhoEtaRange");
252 size = param.GetSize();
253
254 fEstimators.clear();
255 for(i = 0; i < size/2; ++i)
256 {
257 etaMin = param[i*2].GetDouble();
258 etaMax = param[i*2 + 1].GetDouble();
259 estimatorStruct.estimator = new JetMedianBackgroundEstimator(SelectorRapRange(etaMin, etaMax), *fDefinition, *fAreaDefinition);
260 estimatorStruct.etaMin = etaMin;
261 estimatorStruct.etaMax = etaMax;
262 fEstimators.push_back(estimatorStruct);
263 }
264 }
265
266 // import input array
267
268 fInputArray = ImportArray(GetString("InputArray", "Calorimeter/towers"));
269 fItInputArray = fInputArray->MakeIterator();
270
271 // create output arrays
272
273 fOutputArray = ExportArray(GetString("OutputArray", "jets"));
274 fRhoOutputArray = ExportArray(GetString("RhoOutputArray", "rho"));
275}
276
277//------------------------------------------------------------------------------
278
279void FastJetFinder::Finish()
280{
281 vector< TEstimatorStruct >::iterator itEstimators;
282
283 for(itEstimators = fEstimators.begin(); itEstimators != fEstimators.end(); ++itEstimators)
284 {
285 if(itEstimators->estimator) delete itEstimators->estimator;
286 }
287
288 if(fItInputArray) delete fItInputArray;
289 if(fDefinition) delete fDefinition;
290 if(fAreaDefinition) delete fAreaDefinition;
291 if(fPlugin) delete static_cast<JetDefinition::Plugin*>(fPlugin);
292 if(fRecomb) delete static_cast<JetDefinition::Recombiner*>(fRecomb);
293 if(fNjettinessPlugin) delete static_cast<JetDefinition::Plugin*>(fNjettinessPlugin);
294 if(fAxesDef) delete fAxesDef;
295 if(fMeasureDef) delete fMeasureDef;
296}
297
298//------------------------------------------------------------------------------
299
300void FastJetFinder::Process()
301{
302 Candidate *candidate, *constituent;
303 TLorentzVector momentum;
304
305 Double_t deta, dphi, detaMax, dphiMax;
306 Double_t time, timeWeight;
307 Int_t number, ncharged, nneutrals;
308 Int_t charge;
309 Double_t rho = 0.0;
310 PseudoJet jet, area;
311 ClusterSequence *sequence;
312 vector< PseudoJet > inputList, outputList, subjets;
313 vector< PseudoJet >::iterator itInputList, itOutputList;
314 vector< TEstimatorStruct >::iterator itEstimators;
315
316 DelphesFactory *factory = GetFactory();
317
318 inputList.clear();
319
320 // loop over input objects
321 fItInputArray->Reset();
322 number = 0;
323 while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
324 {
325 momentum = candidate->Momentum;
326 jet = PseudoJet(momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E());
327 jet.set_user_index(number);
328 inputList.push_back(jet);
329 ++number;
330 }
331
332 // construct jets
333 if(fAreaDefinition)
334 {
335 sequence = new ClusterSequenceArea(inputList, *fDefinition, *fAreaDefinition);
336 }
337 else
338 {
339 sequence = new ClusterSequence(inputList, *fDefinition);
340 }
341
342 // compute rho and store it
343 if(fComputeRho && fAreaDefinition)
344 {
345 for(itEstimators = fEstimators.begin(); itEstimators != fEstimators.end(); ++itEstimators)
346 {
347 itEstimators->estimator->set_particles(inputList);
348 rho = itEstimators->estimator->rho();
349
350 candidate = factory->NewCandidate();
351 candidate->Momentum.SetPtEtaPhiE(rho, 0.0, 0.0, rho);
352 candidate->Edges[0] = itEstimators->etaMin;
353 candidate->Edges[1] = itEstimators->etaMax;
354 fRhoOutputArray->Add(candidate);
355 }
356 }
357
358 outputList.clear();
359 outputList = sorted_by_pt(sequence->inclusive_jets(fJetPTMin));
360
361
362 // loop over all jets and export them
363 detaMax = 0.0;
364 dphiMax = 0.0;
365
366 for(itOutputList = outputList.begin(); itOutputList != outputList.end(); ++itOutputList)
367 {
368 jet = *itOutputList;
369 if(fJetAlgorithm == 7) jet = join(jet.constituents());
370
371 momentum.SetPxPyPzE(jet.px(), jet.py(), jet.pz(), jet.E());
372
373 area.reset(0.0, 0.0, 0.0, 0.0);
374 if(fAreaDefinition) area = itOutputList->area_4vector();
375
376 candidate = factory->NewCandidate();
377
378 time = 0.0;
379 timeWeight = 0.0;
380
381 charge = 0;
382
383 ncharged = 0;
384 nneutrals = 0;
385
386 inputList.clear();
387 inputList = sequence->constituents(*itOutputList);
388
389 for(itInputList = inputList.begin(); itInputList != inputList.end(); ++itInputList)
390 {
391 if(itInputList->user_index() < 0) continue;
392 constituent = static_cast<Candidate*>(fInputArray->At(itInputList->user_index()));
393
394 deta = TMath::Abs(momentum.Eta() - constituent->Momentum.Eta());
395 dphi = TMath::Abs(momentum.DeltaPhi(constituent->Momentum));
396 if(deta > detaMax) detaMax = deta;
397 if(dphi > dphiMax) dphiMax = dphi;
398
399 if(constituent->Charge == 0) nneutrals++;
400 else ncharged++;
401
402 time += TMath::Sqrt(constituent->Momentum.E())*(constituent->Position.T());
403 timeWeight += TMath::Sqrt(constituent->Momentum.E());
404
405 charge += constituent->Charge;
406
407 candidate->AddCandidate(constituent);
408 }
409
410 candidate->Momentum = momentum;
411 candidate->Position.SetT(time/timeWeight);
412 candidate->Area.SetPxPyPzE(area.px(), area.py(), area.pz(), area.E());
413
414 candidate->DeltaEta = detaMax;
415 candidate->DeltaPhi = dphiMax;
416 candidate->Charge = charge;
417 candidate->NNeutrals = nneutrals;
418 candidate->NCharged = ncharged;
419
420 //------------------------------------
421 // Trimming
422 //------------------------------------
423
424 if(fComputeTrimming)
425 {
426
427 fastjet::Filter trimmer(fastjet::JetDefinition(fastjet::kt_algorithm,fRTrim),fastjet::SelectorPtFractionMin(fPtFracTrim));
428 fastjet::PseudoJet trimmed_jet = trimmer(*itOutputList);
429
430 trimmed_jet = join(trimmed_jet.constituents());
431
432 candidate->TrimmedP4[0].SetPtEtaPhiM(trimmed_jet.pt(), trimmed_jet.eta(), trimmed_jet.phi(), trimmed_jet.m());
433
434 // four hardest subjets
435 subjets.clear();
436 subjets = trimmed_jet.pieces();
437 subjets = sorted_by_pt(subjets);
438
439 candidate->NSubJetsTrimmed = subjets.size();
440
441 for (size_t i = 0; i < subjets.size() and i < 4; i++)
442 {
443 if(subjets.at(i).pt() < 0) continue ;
444 candidate->TrimmedP4[i+1].SetPtEtaPhiM(subjets.at(i).pt(), subjets.at(i).eta(), subjets.at(i).phi(), subjets.at(i).m());
445 }
446 }
447
448
449 //------------------------------------
450 // Pruning
451 //------------------------------------
452
453
454 if(fComputePruning)
455 {
456
457 fastjet::Pruner pruner(fastjet::JetDefinition(fastjet::cambridge_algorithm,fRPrun),fZcutPrun,fRcutPrun);
458 fastjet::PseudoJet pruned_jet = pruner(*itOutputList);
459
460 candidate->PrunedP4[0].SetPtEtaPhiM(pruned_jet.pt(), pruned_jet.eta(), pruned_jet.phi(), pruned_jet.m());
461
462 // four hardest subjet
463 subjets.clear();
464 subjets = pruned_jet.pieces();
465 subjets = sorted_by_pt(subjets);
466
467 candidate->NSubJetsPruned = subjets.size();
468
469 for (size_t i = 0; i < subjets.size() and i < 4; i++)
470 {
471 if(subjets.at(i).pt() < 0) continue ;
472 candidate->PrunedP4[i+1].SetPtEtaPhiM(subjets.at(i).pt(), subjets.at(i).eta(), subjets.at(i).phi(), subjets.at(i).m());
473 }
474
475 }
476
477 //------------------------------------
478 // SoftDrop
479 //------------------------------------
480
481 if(fComputeSoftDrop)
482 {
483
484 contrib::SoftDrop softDrop(fBetaSoftDrop,fSymmetryCutSoftDrop,fR0SoftDrop);
485 fastjet::PseudoJet softdrop_jet = softDrop(*itOutputList);
486
487 candidate->SoftDroppedP4[0].SetPtEtaPhiM(softdrop_jet.pt(), softdrop_jet.eta(), softdrop_jet.phi(), softdrop_jet.m());
488
489 // four hardest subjet
490
491 subjets.clear();
492 subjets = softdrop_jet.pieces();
493 subjets = sorted_by_pt(subjets);
494 candidate->NSubJetsSoftDropped = softdrop_jet.pieces().size();
495
496 for (size_t i = 0; i < subjets.size() and i < 4; i++)
497 {
498 if(subjets.at(i).pt() < 0) continue ;
499 candidate->SoftDroppedP4[i+1].SetPtEtaPhiM(subjets.at(i).pt(), subjets.at(i).eta(), subjets.at(i).phi(), subjets.at(i).m());
500 }
501 }
502
503 // --- compute N-subjettiness with N = 1,2,3,4,5 ----
504
505 if(fComputeNsubjettiness)
506 {
507
508 Nsubjettiness nSub1(1, *fAxesDef, *fMeasureDef);
509 Nsubjettiness nSub2(2, *fAxesDef, *fMeasureDef);
510 Nsubjettiness nSub3(3, *fAxesDef, *fMeasureDef);
511 Nsubjettiness nSub4(4, *fAxesDef, *fMeasureDef);
512 Nsubjettiness nSub5(5, *fAxesDef, *fMeasureDef);
513
514 candidate->Tau[0] = nSub1(*itOutputList);
515 candidate->Tau[1] = nSub2(*itOutputList);
516 candidate->Tau[2] = nSub3(*itOutputList);
517 candidate->Tau[3] = nSub4(*itOutputList);
518 candidate->Tau[4] = nSub5(*itOutputList);
519
520 }
521
522 fOutputArray->Add(candidate);
523 }
524 delete sequence;
525}
Note: See TracBrowser for help on using the repository browser.