Fork me on GitHub

source: git/modules/FastJetFinder.cc@ 45e58be

ImprovedOutputFile Timing llp
Last change on this file since 45e58be was 70bb4cb, checked in by Pavel Demin <pavel-demin@…>, 6 years ago

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