Fork me on GitHub

source: git/modules/FastJetFinder.cc@ c6f9311

Last change on this file since c6f9311 was f9aeea6, checked in by GitHub <noreply@…>, 4 years ago

Update FastJetFinder.cc

Removed comments, fixed issue with Trimmed Jet Constituent Count

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