Fork me on GitHub

source: git/modules/FastJetFinder.cc@ ededa33

ImprovedOutputFile Timing
Last change on this file since ededa33 was 341014c, checked in by Pavel Demin <pavel-demin@…>, 5 years ago

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