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 |
|
---|
73 | using namespace std;
|
---|
74 | using namespace fastjet;
|
---|
75 | using namespace fastjet::contrib;
|
---|
76 |
|
---|
77 | //------------------------------------------------------------------------------
|
---|
78 |
|
---|
79 | FastJetFinder::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 |
|
---|
87 | FastJetFinder::~FastJetFinder()
|
---|
88 | {
|
---|
89 | }
|
---|
90 |
|
---|
91 | //------------------------------------------------------------------------------
|
---|
92 |
|
---|
93 | void 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 | fParameterP = GetDouble("ParameterP", -1.0);
|
---|
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);
|
---|
114 | fAdjacencyCut = GetInt("AdjacencyCut", 2);
|
---|
115 | fOverlapThreshold = GetDouble("OverlapThreshold", 0.75);
|
---|
116 |
|
---|
117 | fJetPTMin = GetDouble("JetPTMin", 10.0);
|
---|
118 |
|
---|
119 | //-- N(sub)jettiness parameters --
|
---|
120 |
|
---|
121 | fComputeNsubjettiness = GetBool("ComputeNsubjettiness", false);
|
---|
122 | fBeta = GetDouble("Beta", 1.0);
|
---|
123 | fAxisMode = GetInt("AxisMode", 1);
|
---|
124 | fRcutOff = GetDouble("RcutOff", 0.8); // used only if Njettiness is used as jet clustering algo (case 8)
|
---|
125 | fN = GetInt("N", 2); // used only if Njettiness is used as jet clustering algo (case 8)
|
---|
126 |
|
---|
127 | //-- Exclusive clustering for e+e- collisions --
|
---|
128 |
|
---|
129 | fNJets = GetInt("NJets", 2);
|
---|
130 | fExclusiveClustering = GetBool("ExclusiveClustering", false);
|
---|
131 |
|
---|
132 | //-- Valencia Linear Collider algorithm
|
---|
133 |
|
---|
134 | fGamma = GetDouble("Gamma", 1.0);
|
---|
135 | //fBeta parameter see above
|
---|
136 |
|
---|
137 | fMeasureDef = new NormalizedMeasure(fBeta, fParameterR);
|
---|
138 |
|
---|
139 | switch(fAxisMode)
|
---|
140 | {
|
---|
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();
|
---|
153 | }
|
---|
154 |
|
---|
155 | //-- Trimming parameters --
|
---|
156 |
|
---|
157 | fComputeTrimming = GetBool("ComputeTrimming", false);
|
---|
158 | fRTrim = GetDouble("RTrim", 0.2);
|
---|
159 | fPtFracTrim = GetDouble("PtFracTrim", 0.05);
|
---|
160 |
|
---|
161 | //-- Pruning parameters --
|
---|
162 |
|
---|
163 | fComputePruning = GetBool("ComputePruning", false);
|
---|
164 | fZcutPrun = GetDouble("ZcutPrun", 0.1);
|
---|
165 | fRcutPrun = GetDouble("RcutPrun", 0.5);
|
---|
166 | fRPrun = GetDouble("RPrun", 0.8);
|
---|
167 |
|
---|
168 | //-- SoftDrop parameters --
|
---|
169 |
|
---|
170 | fComputeSoftDrop = GetBool("ComputeSoftDrop", false);
|
---|
171 | fBetaSoftDrop = GetDouble("BetaSoftDrop", 0.0);
|
---|
172 | fSymmetryCutSoftDrop = GetDouble("SymmetryCutSoftDrop", 0.1);
|
---|
173 | fR0SoftDrop = GetDouble("R0SoftDrop=", 0.8);
|
---|
174 |
|
---|
175 | // --- Jet Area Parameters ---
|
---|
176 |
|
---|
177 | fAreaAlgorithm = GetInt("AreaAlgorithm", 0);
|
---|
178 | fComputeRho = GetBool("ComputeRho", false);
|
---|
179 |
|
---|
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);
|
---|
187 |
|
---|
188 | // - voronoi based areas -
|
---|
189 | fEffectiveRfact = GetDouble("EffectiveRfact", 1.0);
|
---|
190 |
|
---|
191 | switch(fAreaAlgorithm)
|
---|
192 | {
|
---|
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;
|
---|
212 | }
|
---|
213 |
|
---|
214 | switch(fJetAlgorithm)
|
---|
215 | {
|
---|
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;
|
---|
250 | case 10:
|
---|
251 | fDefinition = new JetDefinition(ee_genkt_algorithm,fParameterR,fParameterP);
|
---|
252 | break;
|
---|
253 |
|
---|
254 | }
|
---|
255 |
|
---|
256 | fPlugin = plugin;
|
---|
257 | fRecomb = recomb;
|
---|
258 |
|
---|
259 | ClusterSequence::print_banner();
|
---|
260 |
|
---|
261 | if(fComputeRho && fAreaDefinition)
|
---|
262 | {
|
---|
263 | // read eta ranges
|
---|
264 |
|
---|
265 | param = GetParam("RhoEtaRange");
|
---|
266 | size = param.GetSize();
|
---|
267 |
|
---|
268 | fEstimators.clear();
|
---|
269 | for(i = 0; i < size / 2; ++i)
|
---|
270 | {
|
---|
271 | etaMin = param[i * 2].GetDouble();
|
---|
272 | etaMax = param[i * 2 + 1].GetDouble();
|
---|
273 | estimatorStruct.estimator = new JetMedianBackgroundEstimator(SelectorRapRange(etaMin, etaMax), *fDefinition, *fAreaDefinition);
|
---|
274 | estimatorStruct.etaMin = etaMin;
|
---|
275 | estimatorStruct.etaMax = etaMax;
|
---|
276 | fEstimators.push_back(estimatorStruct);
|
---|
277 | }
|
---|
278 | }
|
---|
279 |
|
---|
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"));
|
---|
289 | fConstituentsOutputArray = ExportArray(GetString("ConstituentsOutputArray", "constituents"));
|
---|
290 | }
|
---|
291 |
|
---|
292 | //------------------------------------------------------------------------------
|
---|
293 |
|
---|
294 | void FastJetFinder::Finish()
|
---|
295 | {
|
---|
296 | vector<TEstimatorStruct>::iterator itEstimators;
|
---|
297 |
|
---|
298 | for(itEstimators = fEstimators.begin(); itEstimators != fEstimators.end(); ++itEstimators)
|
---|
299 | {
|
---|
300 | if(itEstimators->estimator) delete itEstimators->estimator;
|
---|
301 | }
|
---|
302 |
|
---|
303 | if(fItInputArray) delete fItInputArray;
|
---|
304 | if(fDefinition) delete fDefinition;
|
---|
305 | if(fAreaDefinition) delete fAreaDefinition;
|
---|
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);
|
---|
309 | if(fAxesDef) delete fAxesDef;
|
---|
310 | if(fMeasureDef) delete fMeasureDef;
|
---|
311 | if(fValenciaPlugin) delete static_cast<JetDefinition::Plugin *>(fValenciaPlugin);
|
---|
312 | }
|
---|
313 |
|
---|
314 | //------------------------------------------------------------------------------
|
---|
315 |
|
---|
316 | void FastJetFinder::Process()
|
---|
317 | {
|
---|
318 | Candidate *candidate, *constituent;
|
---|
319 | TLorentzVector momentum;
|
---|
320 |
|
---|
321 | Double_t deta, dphi, detaMax, dphiMax;
|
---|
322 | Double_t time, timeWeight;
|
---|
323 | Double_t neutralEnergyFraction, chargedEnergyFraction;
|
---|
324 |
|
---|
325 | Int_t number, ncharged, nneutrals;
|
---|
326 | Int_t charge;
|
---|
327 | Double_t rho = 0.0;
|
---|
328 | PseudoJet jet, area;
|
---|
329 | ClusterSequence *sequence;
|
---|
330 | vector<PseudoJet> inputList, outputList, subjets;
|
---|
331 | vector<PseudoJet>::iterator itInputList, itOutputList;
|
---|
332 | vector<TEstimatorStruct>::iterator itEstimators;
|
---|
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;
|
---|
337 |
|
---|
338 | DelphesFactory *factory = GetFactory();
|
---|
339 |
|
---|
340 | inputList.clear();
|
---|
341 |
|
---|
342 | // loop over input objects
|
---|
343 | fItInputArray->Reset();
|
---|
344 | number = 0;
|
---|
345 | while((candidate = static_cast<Candidate *>(fItInputArray->Next())))
|
---|
346 | {
|
---|
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 | }
|
---|
353 |
|
---|
354 | // construct jets
|
---|
355 | if(fAreaDefinition)
|
---|
356 | {
|
---|
357 | sequence = new ClusterSequenceArea(inputList, *fDefinition, *fAreaDefinition);
|
---|
358 | }
|
---|
359 | else
|
---|
360 | {
|
---|
361 | sequence = new ClusterSequence(inputList, *fDefinition);
|
---|
362 | }
|
---|
363 |
|
---|
364 | // compute rho and store it
|
---|
365 | if(fComputeRho && fAreaDefinition)
|
---|
366 | {
|
---|
367 | for(itEstimators = fEstimators.begin(); itEstimators != fEstimators.end(); ++itEstimators)
|
---|
368 | {
|
---|
369 | itEstimators->estimator->set_particles(inputList);
|
---|
370 | rho = itEstimators->estimator->rho();
|
---|
371 |
|
---|
372 | candidate = factory->NewCandidate();
|
---|
373 | candidate->Momentum.SetPtEtaPhiE(rho, 0.0, 0.0, rho);
|
---|
374 | candidate->Edges[0] = itEstimators->etaMin;
|
---|
375 | candidate->Edges[1] = itEstimators->etaMax;
|
---|
376 | fRhoOutputArray->Add(candidate);
|
---|
377 | }
|
---|
378 | }
|
---|
379 |
|
---|
380 | outputList.clear();
|
---|
381 |
|
---|
382 | if(fExclusiveClustering)
|
---|
383 | {
|
---|
384 | try
|
---|
385 | {
|
---|
386 | outputList = sorted_by_pt(sequence->exclusive_jets(fNJets));
|
---|
387 | }
|
---|
388 | catch(fastjet::Error)
|
---|
389 | {
|
---|
390 | outputList.clear();
|
---|
391 | }
|
---|
392 |
|
---|
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 |
|
---|
403 | // loop over all jets and export them
|
---|
404 | detaMax = 0.0;
|
---|
405 | dphiMax = 0.0;
|
---|
406 |
|
---|
407 | for(itOutputList = outputList.begin(); itOutputList != outputList.end(); ++itOutputList)
|
---|
408 | {
|
---|
409 | jet = *itOutputList;
|
---|
410 | if(fJetAlgorithm == 7) jet = join(jet.constituents());
|
---|
411 |
|
---|
412 | momentum.SetPxPyPzE(jet.px(), jet.py(), jet.pz(), jet.E());
|
---|
413 |
|
---|
414 | area.reset(0.0, 0.0, 0.0, 0.0);
|
---|
415 | if(fAreaDefinition) area = itOutputList->area_4vector();
|
---|
416 |
|
---|
417 | candidate = factory->NewCandidate();
|
---|
418 |
|
---|
419 | time = 0.0;
|
---|
420 | timeWeight = 0.0;
|
---|
421 |
|
---|
422 | charge = 0;
|
---|
423 |
|
---|
424 | ncharged = 0;
|
---|
425 | nneutrals = 0;
|
---|
426 |
|
---|
427 | neutralEnergyFraction =0.;
|
---|
428 | chargedEnergyFraction =0.;
|
---|
429 |
|
---|
430 | inputList.clear();
|
---|
431 | inputList = sequence->constituents(*itOutputList);
|
---|
432 |
|
---|
433 | for(itInputList = inputList.begin(); itInputList != inputList.end(); ++itInputList)
|
---|
434 | {
|
---|
435 | if(itInputList->user_index() < 0) continue;
|
---|
436 | constituent = static_cast<Candidate *>(fInputArray->At(itInputList->user_index()));
|
---|
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;
|
---|
442 |
|
---|
443 | if(constituent->Charge == 0)
|
---|
444 | {
|
---|
445 | nneutrals++;
|
---|
446 | neutralEnergyFraction += constituent->Momentum.E();
|
---|
447 | }
|
---|
448 | else
|
---|
449 | {
|
---|
450 | ncharged++;
|
---|
451 | chargedEnergyFraction += constituent->Momentum.E();
|
---|
452 | }
|
---|
453 |
|
---|
454 | time += TMath::Sqrt(constituent->Momentum.E()) * (constituent->Position.T());
|
---|
455 | timeWeight += TMath::Sqrt(constituent->Momentum.E());
|
---|
456 |
|
---|
457 | charge += constituent->Charge;
|
---|
458 |
|
---|
459 | fConstituentsOutputArray->Add(constituent);
|
---|
460 | candidate->AddCandidate(constituent);
|
---|
461 | }
|
---|
462 |
|
---|
463 | candidate->Momentum = momentum;
|
---|
464 | candidate->Position.SetT(time / timeWeight);
|
---|
465 | candidate->Area.SetPxPyPzE(area.px(), area.py(), area.pz(), area.E());
|
---|
466 |
|
---|
467 | candidate->DeltaEta = detaMax;
|
---|
468 | candidate->DeltaPhi = dphiMax;
|
---|
469 | candidate->Charge = charge;
|
---|
470 | candidate->NNeutrals = nneutrals;
|
---|
471 | candidate->NCharged = ncharged;
|
---|
472 |
|
---|
473 | candidate->NeutralEnergyFraction = (momentum.E() > 0 ) ? neutralEnergyFraction/momentum.E() : 0.0;
|
---|
474 | candidate->ChargedEnergyFraction = (momentum.E() > 0 ) ? chargedEnergyFraction/momentum.E() : 0.0;
|
---|
475 |
|
---|
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;
|
---|
481 |
|
---|
482 | //------------------------------------
|
---|
483 | // Trimming
|
---|
484 | //------------------------------------
|
---|
485 |
|
---|
486 | if(fComputeTrimming)
|
---|
487 | {
|
---|
488 |
|
---|
489 | fastjet::Filter trimmer(fastjet::JetDefinition(fastjet::kt_algorithm, fRTrim), fastjet::SelectorPtFractionMin(fPtFracTrim));
|
---|
490 | fastjet::PseudoJet trimmed_jet = trimmer(*itOutputList);
|
---|
491 |
|
---|
492 | candidate->TrimmedP4[0].SetPtEtaPhiM(trimmed_jet.pt(), trimmed_jet.eta(), trimmed_jet.phi(), trimmed_jet.m());
|
---|
493 |
|
---|
494 | // four hardest subjets
|
---|
495 | subjets.clear();
|
---|
496 | subjets = trimmed_jet.pieces();
|
---|
497 | subjets = sorted_by_pt(subjets);
|
---|
498 |
|
---|
499 | candidate->NSubJetsTrimmed = subjets.size();
|
---|
500 |
|
---|
501 | for(size_t i = 0; i < subjets.size() and i < 4; i++)
|
---|
502 | {
|
---|
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());
|
---|
505 | }
|
---|
506 | }
|
---|
507 |
|
---|
508 | //------------------------------------
|
---|
509 | // Pruning
|
---|
510 | //------------------------------------
|
---|
511 |
|
---|
512 | if(fComputePruning)
|
---|
513 | {
|
---|
514 |
|
---|
515 | fastjet::Pruner pruner(fastjet::JetDefinition(fastjet::cambridge_algorithm, fRPrun), fZcutPrun, fRcutPrun);
|
---|
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());
|
---|
519 |
|
---|
520 | // four hardest subjet
|
---|
521 | subjets.clear();
|
---|
522 | subjets = pruned_jet.pieces();
|
---|
523 | subjets = sorted_by_pt(subjets);
|
---|
524 |
|
---|
525 | candidate->NSubJetsPruned = subjets.size();
|
---|
526 |
|
---|
527 | for(size_t i = 0; i < subjets.size() and i < 4; i++)
|
---|
528 | {
|
---|
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());
|
---|
531 | }
|
---|
532 | }
|
---|
533 |
|
---|
534 | //------------------------------------
|
---|
535 | // SoftDrop
|
---|
536 | //------------------------------------
|
---|
537 |
|
---|
538 | if(fComputeSoftDrop)
|
---|
539 | {
|
---|
540 |
|
---|
541 | contrib::SoftDrop softDrop(fBetaSoftDrop, fSymmetryCutSoftDrop, fR0SoftDrop);
|
---|
542 | fastjet::PseudoJet softdrop_jet = softDrop(*itOutputList);
|
---|
543 |
|
---|
544 | candidate->SoftDroppedP4[0].SetPtEtaPhiM(softdrop_jet.pt(), softdrop_jet.eta(), softdrop_jet.phi(), softdrop_jet.m());
|
---|
545 |
|
---|
546 | // four hardest subjet
|
---|
547 |
|
---|
548 | subjets.clear();
|
---|
549 | subjets = softdrop_jet.pieces();
|
---|
550 | subjets = sorted_by_pt(subjets);
|
---|
551 | candidate->NSubJetsSoftDropped = softdrop_jet.pieces().size();
|
---|
552 |
|
---|
553 | candidate->SoftDroppedJet = candidate->SoftDroppedP4[0];
|
---|
554 |
|
---|
555 | for(size_t i = 0; i < subjets.size() and i < 4; i++)
|
---|
556 | {
|
---|
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];
|
---|
561 | }
|
---|
562 | }
|
---|
563 |
|
---|
564 | // --- compute N-subjettiness with N = 1,2,3,4,5 ----
|
---|
565 |
|
---|
566 | if(fComputeNsubjettiness)
|
---|
567 | {
|
---|
568 |
|
---|
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);
|
---|
574 |
|
---|
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);
|
---|
580 | }
|
---|
581 |
|
---|
582 | fOutputArray->Add(candidate);
|
---|
583 | }
|
---|
584 | delete sequence;
|
---|
585 | }
|
---|