Fork me on GitHub

source: git/modules/FastJetFinder.cc@ 9da65a5

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 9da65a5 was ce4feac, checked in by Michele Selvaggi <michele.selvaggi@…>, 9 years ago

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