Fork me on GitHub

source: git/modules/FastJetFinder.cc@ ccea85b

ImprovedOutputFile
Last change on this file since ccea85b was c614dd7, checked in by Michele Selvaggi <michele.selvaggi@…>, 5 years ago

added neutral and charged energy fraction to jets
  • Property mode set to 100644
File size: 17.9 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 Double_t neutralEnergyFraction, chargedEnergyFraction;
319
320 Int_t number, ncharged, nneutrals;
321 Int_t charge;
322 Double_t rho = 0.0;
323 PseudoJet jet, area;
324 ClusterSequence *sequence;
325 vector<PseudoJet> inputList, outputList, subjets;
326 vector<PseudoJet>::iterator itInputList, itOutputList;
327 vector<TEstimatorStruct>::iterator itEstimators;
328 Double_t excl_ymerge23 = 0.0;
329 Double_t excl_ymerge34 = 0.0;
330 Double_t excl_ymerge45 = 0.0;
331 Double_t excl_ymerge56 = 0.0;
332
333 DelphesFactory *factory = GetFactory();
334
335 inputList.clear();
336
337 // loop over input objects
338 fItInputArray->Reset();
339 number = 0;
340 while((candidate = static_cast<Candidate *>(fItInputArray->Next())))
341 {
342 momentum = candidate->Momentum;
343 jet = PseudoJet(momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E());
344 jet.set_user_index(number);
345 inputList.push_back(jet);
346 ++number;
347 }
348
349 // construct jets
350 if(fAreaDefinition)
351 {
352 sequence = new ClusterSequenceArea(inputList, *fDefinition, *fAreaDefinition);
353 }
354 else
355 {
356 sequence = new ClusterSequence(inputList, *fDefinition);
357 }
358
359 // compute rho and store it
360 if(fComputeRho && fAreaDefinition)
361 {
362 for(itEstimators = fEstimators.begin(); itEstimators != fEstimators.end(); ++itEstimators)
363 {
364 itEstimators->estimator->set_particles(inputList);
365 rho = itEstimators->estimator->rho();
366
367 candidate = factory->NewCandidate();
368 candidate->Momentum.SetPtEtaPhiE(rho, 0.0, 0.0, rho);
369 candidate->Edges[0] = itEstimators->etaMin;
370 candidate->Edges[1] = itEstimators->etaMax;
371 fRhoOutputArray->Add(candidate);
372 }
373 }
374
375 outputList.clear();
376
377 if(fExclusiveClustering)
378 {
379 try
380 {
381 outputList = sorted_by_pt(sequence->exclusive_jets(fNJets));
382 }
383 catch(fastjet::Error)
384 {
385 outputList.clear();
386 }
387
388 excl_ymerge23 = sequence->exclusive_ymerge(2);
389 excl_ymerge34 = sequence->exclusive_ymerge(3);
390 excl_ymerge45 = sequence->exclusive_ymerge(4);
391 excl_ymerge56 = sequence->exclusive_ymerge(5);
392 }
393 else
394 {
395 outputList = sorted_by_pt(sequence->inclusive_jets(fJetPTMin));
396 }
397
398 // loop over all jets and export them
399 detaMax = 0.0;
400 dphiMax = 0.0;
401
402 for(itOutputList = outputList.begin(); itOutputList != outputList.end(); ++itOutputList)
403 {
404 jet = *itOutputList;
405 if(fJetAlgorithm == 7) jet = join(jet.constituents());
406
407 momentum.SetPxPyPzE(jet.px(), jet.py(), jet.pz(), jet.E());
408
409 area.reset(0.0, 0.0, 0.0, 0.0);
410 if(fAreaDefinition) area = itOutputList->area_4vector();
411
412 candidate = factory->NewCandidate();
413
414 time = 0.0;
415 timeWeight = 0.0;
416
417 charge = 0;
418
419 ncharged = 0;
420 nneutrals = 0;
421
422 neutralEnergyFraction =0.;
423 chargedEnergyFraction =0.;
424
425 inputList.clear();
426 inputList = sequence->constituents(*itOutputList);
427
428 for(itInputList = inputList.begin(); itInputList != inputList.end(); ++itInputList)
429 {
430 if(itInputList->user_index() < 0) continue;
431 constituent = static_cast<Candidate *>(fInputArray->At(itInputList->user_index()));
432
433 deta = TMath::Abs(momentum.Eta() - constituent->Momentum.Eta());
434 dphi = TMath::Abs(momentum.DeltaPhi(constituent->Momentum));
435 if(deta > detaMax) detaMax = deta;
436 if(dphi > dphiMax) dphiMax = dphi;
437
438 if(constituent->Charge == 0)
439 {
440 nneutrals++;
441 neutralEnergyFraction += constituent->Momentum.E();
442 }
443 else
444 {
445 ncharged++;
446 chargedEnergyFraction += constituent->Momentum.E();
447 }
448
449 time += TMath::Sqrt(constituent->Momentum.E()) * (constituent->Position.T());
450 timeWeight += TMath::Sqrt(constituent->Momentum.E());
451
452 charge += constituent->Charge;
453
454 fConstituentsOutputArray->Add(constituent);
455 candidate->AddCandidate(constituent);
456 }
457
458 candidate->Momentum = momentum;
459 candidate->Position.SetT(time / timeWeight);
460 candidate->Area.SetPxPyPzE(area.px(), area.py(), area.pz(), area.E());
461
462 candidate->DeltaEta = detaMax;
463 candidate->DeltaPhi = dphiMax;
464 candidate->Charge = charge;
465 candidate->NNeutrals = nneutrals;
466 candidate->NCharged = ncharged;
467
468 candidate->NeutralEnergyFraction = (momentum.E() > 0 ) ? neutralEnergyFraction/momentum.E() : 0.0;
469 candidate->ChargedEnergyFraction = (momentum.E() > 0 ) ? chargedEnergyFraction/momentum.E() : 0.0;
470
471 //for exclusive clustering, access y_n,n+1 as exclusive_ymerge (fNJets);
472 candidate->ExclYmerge23 = excl_ymerge23;
473 candidate->ExclYmerge34 = excl_ymerge34;
474 candidate->ExclYmerge45 = excl_ymerge45;
475 candidate->ExclYmerge56 = excl_ymerge56;
476
477 //------------------------------------
478 // Trimming
479 //------------------------------------
480
481 if(fComputeTrimming)
482 {
483
484 fastjet::Filter trimmer(fastjet::JetDefinition(fastjet::kt_algorithm, fRTrim), fastjet::SelectorPtFractionMin(fPtFracTrim));
485 fastjet::PseudoJet trimmed_jet = trimmer(*itOutputList);
486
487 trimmed_jet = join(trimmed_jet.constituents());
488
489 candidate->TrimmedP4[0].SetPtEtaPhiM(trimmed_jet.pt(), trimmed_jet.eta(), trimmed_jet.phi(), trimmed_jet.m());
490
491 // four hardest subjets
492 subjets.clear();
493 subjets = trimmed_jet.pieces();
494 subjets = sorted_by_pt(subjets);
495
496 candidate->NSubJetsTrimmed = subjets.size();
497
498 for(size_t i = 0; i < subjets.size() and i < 4; i++)
499 {
500 if(subjets.at(i).pt() < 0) continue;
501 candidate->TrimmedP4[i + 1].SetPtEtaPhiM(subjets.at(i).pt(), subjets.at(i).eta(), subjets.at(i).phi(), subjets.at(i).m());
502 }
503 }
504
505 //------------------------------------
506 // Pruning
507 //------------------------------------
508
509 if(fComputePruning)
510 {
511
512 fastjet::Pruner pruner(fastjet::JetDefinition(fastjet::cambridge_algorithm, fRPrun), fZcutPrun, fRcutPrun);
513 fastjet::PseudoJet pruned_jet = pruner(*itOutputList);
514
515 candidate->PrunedP4[0].SetPtEtaPhiM(pruned_jet.pt(), pruned_jet.eta(), pruned_jet.phi(), pruned_jet.m());
516
517 // four hardest subjet
518 subjets.clear();
519 subjets = pruned_jet.pieces();
520 subjets = sorted_by_pt(subjets);
521
522 candidate->NSubJetsPruned = subjets.size();
523
524 for(size_t i = 0; i < subjets.size() and i < 4; i++)
525 {
526 if(subjets.at(i).pt() < 0) continue;
527 candidate->PrunedP4[i + 1].SetPtEtaPhiM(subjets.at(i).pt(), subjets.at(i).eta(), subjets.at(i).phi(), subjets.at(i).m());
528 }
529 }
530
531 //------------------------------------
532 // SoftDrop
533 //------------------------------------
534
535 if(fComputeSoftDrop)
536 {
537
538 contrib::SoftDrop softDrop(fBetaSoftDrop, fSymmetryCutSoftDrop, fR0SoftDrop);
539 fastjet::PseudoJet softdrop_jet = softDrop(*itOutputList);
540
541 candidate->SoftDroppedP4[0].SetPtEtaPhiM(softdrop_jet.pt(), softdrop_jet.eta(), softdrop_jet.phi(), softdrop_jet.m());
542
543 // four hardest subjet
544
545 subjets.clear();
546 subjets = softdrop_jet.pieces();
547 subjets = sorted_by_pt(subjets);
548 candidate->NSubJetsSoftDropped = softdrop_jet.pieces().size();
549
550 candidate->SoftDroppedJet = candidate->SoftDroppedP4[0];
551
552 for(size_t i = 0; i < subjets.size() and i < 4; i++)
553 {
554 if(subjets.at(i).pt() < 0) continue;
555 candidate->SoftDroppedP4[i + 1].SetPtEtaPhiM(subjets.at(i).pt(), subjets.at(i).eta(), subjets.at(i).phi(), subjets.at(i).m());
556 if(i == 0) candidate->SoftDroppedSubJet1 = candidate->SoftDroppedP4[i + 1];
557 if(i == 1) candidate->SoftDroppedSubJet2 = candidate->SoftDroppedP4[i + 1];
558 }
559 }
560
561 // --- compute N-subjettiness with N = 1,2,3,4,5 ----
562
563 if(fComputeNsubjettiness)
564 {
565
566 Nsubjettiness nSub1(1, *fAxesDef, *fMeasureDef);
567 Nsubjettiness nSub2(2, *fAxesDef, *fMeasureDef);
568 Nsubjettiness nSub3(3, *fAxesDef, *fMeasureDef);
569 Nsubjettiness nSub4(4, *fAxesDef, *fMeasureDef);
570 Nsubjettiness nSub5(5, *fAxesDef, *fMeasureDef);
571
572 candidate->Tau[0] = nSub1(*itOutputList);
573 candidate->Tau[1] = nSub2(*itOutputList);
574 candidate->Tau[2] = nSub3(*itOutputList);
575 candidate->Tau[3] = nSub4(*itOutputList);
576 candidate->Tau[4] = nSub5(*itOutputList);
577 }
578
579 fOutputArray->Add(candidate);
580 }
581 delete sequence;
582}
Note: See TracBrowser for help on using the repository browser.