Fork me on GitHub

source: git/modules/FastJetFinder.cc@ 0f8390e

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 0f8390e was 1fa50c2, checked in by Pavel Demin <pavel.demin@…>, 10 years ago

fix GPLv3 header
  • Property mode set to 100644
File size: 12.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 * $Date$
25 * $Revision$
26 *
27 *
28 * \author P. Demin - UCL, Louvain-la-Neuve
29 *
30 */
31
32#include "modules/FastJetFinder.h"
33
34#include "classes/DelphesClasses.h"
35#include "classes/DelphesFactory.h"
36#include "classes/DelphesFormula.h"
37
38#include "ExRootAnalysis/ExRootResult.h"
39#include "ExRootAnalysis/ExRootFilter.h"
40#include "ExRootAnalysis/ExRootClassifier.h"
41
42#include "TMath.h"
43#include "TString.h"
44#include "TFormula.h"
45#include "TRandom3.h"
46#include "TObjArray.h"
47#include "TDatabasePDG.h"
48#include "TLorentzVector.h"
49
50#include <algorithm>
51#include <stdexcept>
52#include <iostream>
53#include <sstream>
54#include <vector>
55
56#include "fastjet/PseudoJet.hh"
57#include "fastjet/JetDefinition.hh"
58#include "fastjet/ClusterSequence.hh"
59#include "fastjet/Selector.hh"
60#include "fastjet/ClusterSequenceArea.hh"
61#include "fastjet/tools/JetMedianBackgroundEstimator.hh"
62
63#include "fastjet/plugins/SISCone/fastjet/SISConePlugin.hh"
64#include "fastjet/plugins/CDFCones/fastjet/CDFMidPointPlugin.hh"
65#include "fastjet/plugins/CDFCones/fastjet/CDFJetCluPlugin.hh"
66
67#include "fastjet/contribs/Nsubjettiness/Nsubjettiness.hh"
68#include "fastjet/contribs/Nsubjettiness/Njettiness.hh"
69#include "fastjet/contribs/Nsubjettiness/NjettinessPlugin.hh"
70#include "fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.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 NjettinessPlugin *njetPlugin = 0;
99
100 // read eta ranges
101
102 ExRootConfParam param = GetParam("RhoEtaRange");
103 Long_t i, size;
104
105 fEtaRangeMap.clear();
106 size = param.GetSize();
107 for(i = 0; i < size/2; ++i)
108 {
109 fEtaRangeMap[param[i*2].GetDouble()] = param[i*2 + 1].GetDouble();
110 }
111
112 // define algorithm
113
114 fJetAlgorithm = GetInt("JetAlgorithm", 6);
115 fParameterR = GetDouble("ParameterR", 0.5);
116
117 fConeRadius = GetDouble("ConeRadius", 0.5);
118 fSeedThreshold = GetDouble("SeedThreshold", 1.0);
119 fConeAreaFraction = GetDouble("ConeAreaFraction", 1.0);
120 fMaxIterations = GetInt("MaxIterations", 100);
121 fMaxPairSize = GetInt("MaxPairSize", 2);
122 fIratch = GetInt("Iratch", 1);
123 fAdjacencyCut = GetInt("AdjacencyCut", 2);
124 fOverlapThreshold = GetDouble("OverlapThreshold", 0.75);
125
126 fJetPTMin = GetDouble("JetPTMin", 10.0);
127
128 //-- N(sub)jettiness parameters --
129
130 fComputeNsubjettiness = GetBool("ComputeNsubjettiness", false);
131 fBeta = GetDouble("Beta", 1.0);
132 fAxisMode = GetInt("AxisMode", 1);
133 fRcutOff = GetDouble("RcutOff", 0.8); // used only if Njettiness is used as jet clustering algo (case 8)
134 fN = GetInt("N", 2); // used only if Njettiness is used as jet clustering algo (case 8)
135
136 // --- Jet Area Parameters ---
137 fAreaAlgorithm = GetInt("AreaAlgorithm", 0);
138 fComputeRho = GetBool("ComputeRho", false);
139
140 // - ghost based areas -
141 fGhostEtaMax = GetDouble("GhostEtaMax", 5.0);
142 fRepeat = GetInt("Repeat", 1);
143 fGhostArea = GetDouble("GhostArea", 0.01);
144 fGridScatter = GetDouble("GridScatter", 1.0);
145 fPtScatter = GetDouble("PtScatter", 0.1);
146 fMeanGhostPt = GetDouble("MeanGhostPt", 1.0E-100);
147
148 // - voronoi based areas -
149 fEffectiveRfact = GetDouble("EffectiveRfact", 1.0);
150
151 switch(fAreaAlgorithm)
152 {
153 case 1:
154 fAreaDefinition = new AreaDefinition(active_area_explicit_ghosts, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
155 break;
156 case 2:
157 fAreaDefinition = new AreaDefinition(one_ghost_passive_area, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
158 break;
159 case 3:
160 fAreaDefinition = new AreaDefinition(passive_area, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
161 break;
162 case 4:
163 fAreaDefinition = new AreaDefinition(VoronoiAreaSpec(fEffectiveRfact));
164 break;
165 case 5:
166 fAreaDefinition = new AreaDefinition(active_area, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
167 break;
168 default:
169 case 0:
170 fAreaDefinition = 0;
171 break;
172 }
173
174 switch(fJetAlgorithm)
175 {
176 case 1:
177 plugin = new CDFJetCluPlugin(fSeedThreshold, fConeRadius, fAdjacencyCut, fMaxIterations, fIratch, fOverlapThreshold);
178 fDefinition = new JetDefinition(plugin);
179 break;
180 case 2:
181 plugin = new CDFMidPointPlugin(fSeedThreshold, fConeRadius, fConeAreaFraction, fMaxPairSize, fMaxIterations, fOverlapThreshold);
182 fDefinition = new JetDefinition(plugin);
183 break;
184 case 3:
185 plugin = new SISConePlugin(fConeRadius, fOverlapThreshold, fMaxIterations, fJetPTMin);
186 fDefinition = new JetDefinition(plugin);
187 break;
188 case 4:
189 fDefinition = new JetDefinition(kt_algorithm, fParameterR);
190 break;
191 case 5:
192 fDefinition = new JetDefinition(cambridge_algorithm, fParameterR);
193 break;
194 default:
195 case 6:
196 fDefinition = new JetDefinition(antikt_algorithm, fParameterR);
197 break;
198 case 7:
199 recomb = new WinnerTakeAllRecombiner();
200 fDefinition = new JetDefinition(antikt_algorithm, fParameterR, recomb, Best);
201 break;
202 case 8:
203 njetPlugin = new NjettinessPlugin(fN, Njettiness::wta_kt_axes, Njettiness::unnormalized_cutoff_measure, fBeta, fRcutOff);
204 fDefinition = new JetDefinition(njetPlugin);
205 break;
206 }
207
208 fPlugin = plugin;
209 fRecomb = recomb;
210 fNjettinessPlugin = njetPlugin;
211
212 ClusterSequence::print_banner();
213
214 // import input array
215
216 fInputArray = ImportArray(GetString("InputArray", "Calorimeter/towers"));
217 fItInputArray = fInputArray->MakeIterator();
218
219 // create output arrays
220
221 fOutputArray = ExportArray(GetString("OutputArray", "jets"));
222 fRhoOutputArray = ExportArray(GetString("RhoOutputArray", "rho"));
223}
224
225//------------------------------------------------------------------------------
226
227void FastJetFinder::Finish()
228{
229 if(fItInputArray) delete fItInputArray;
230 if(fDefinition) delete fDefinition;
231 if(fAreaDefinition) delete fAreaDefinition;
232 if(fPlugin) delete static_cast<JetDefinition::Plugin*>(fPlugin);
233 if(fRecomb) delete static_cast<JetDefinition::Recombiner*>(fRecomb);
234 if(fNjettinessPlugin) delete static_cast<JetDefinition::Plugin*>(fNjettinessPlugin);
235}
236
237//------------------------------------------------------------------------------
238
239void FastJetFinder::Process()
240{
241 Candidate *candidate, *constituent;
242 TLorentzVector momentum;
243
244 TLorentzVector constmomentum;
245
246 Double_t deta, dphi, detaMax, dphiMax;
247 Double_t time, weightTime, avTime;
248 Int_t number;
249 Double_t rho = 0;
250 PseudoJet jet, area;
251 vector<PseudoJet> inputList, outputList;
252 ClusterSequence *sequence;
253 map< Double_t, Double_t >::iterator itEtaRangeMap;
254
255 DelphesFactory *factory = GetFactory();
256
257 inputList.clear();
258
259 // loop over input objects
260 fItInputArray->Reset();
261 number = 0;
262 while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
263 {
264 momentum = candidate->Momentum;
265 jet = PseudoJet(momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E());
266 jet.set_user_index(number);
267 inputList.push_back(jet);
268 ++number;
269 }
270
271 // construct jets
272 if(fAreaDefinition)
273 {
274 sequence = new ClusterSequenceArea(inputList, *fDefinition, *fAreaDefinition);
275 }
276 else
277 {
278 sequence = new ClusterSequence(inputList, *fDefinition);
279 }
280
281 // compute rho and store it
282 if(fComputeRho && fAreaDefinition)
283 {
284 for(itEtaRangeMap = fEtaRangeMap.begin(); itEtaRangeMap != fEtaRangeMap.end(); ++itEtaRangeMap)
285 {
286 Selector select_rapidity = SelectorAbsRapRange(itEtaRangeMap->first, itEtaRangeMap->second);
287 JetMedianBackgroundEstimator estimator(select_rapidity, *fDefinition, *fAreaDefinition);
288 estimator.set_particles(inputList);
289 rho = estimator.rho();
290
291 candidate = factory->NewCandidate();
292 candidate->Momentum.SetPtEtaPhiE(rho, 0.0, 0.0, rho);
293 candidate->Edges[0] = itEtaRangeMap->first;
294 candidate->Edges[1] = itEtaRangeMap->second;
295 fRhoOutputArray->Add(candidate);
296 }
297 }
298
299 outputList.clear();
300 outputList = sorted_by_pt(sequence->inclusive_jets(fJetPTMin));
301
302
303 // loop over all jets and export them
304 detaMax = 0.0;
305 dphiMax = 0.0;
306 vector<PseudoJet>::iterator itInputList, itOutputList;
307 for(itOutputList = outputList.begin(); itOutputList != outputList.end(); ++itOutputList)
308 {
309 jet = *itOutputList;
310 if(fJetAlgorithm == 7) jet = join(jet.constituents());
311
312 momentum.SetPxPyPzE(jet.px(), jet.py(), jet.pz(), jet.E());
313
314 area.reset(0.0, 0.0, 0.0, 0.0);
315 if(fAreaDefinition) area = itOutputList->area_4vector();
316
317 candidate = factory->NewCandidate();
318
319 time=0;
320 weightTime=0;
321
322 inputList.clear();
323 inputList = sequence->constituents(*itOutputList);
324
325 constmomentum.SetPxPyPzE(0.0,0.0,0.0,0.0);;
326
327 for(itInputList = inputList.begin(); itInputList != inputList.end(); ++itInputList)
328 {
329 constituent = static_cast<Candidate*>(fInputArray->At(itInputList->user_index()));
330
331 constmomentum += constituent->Momentum;
332
333 deta = TMath::Abs(momentum.Eta() - constituent->Momentum.Eta());
334 dphi = TMath::Abs(momentum.DeltaPhi(constituent->Momentum));
335 if(deta > detaMax) detaMax = deta;
336 if(dphi > dphiMax) dphiMax = dphi;
337
338 time += TMath::Sqrt(constituent->Momentum.E())*(constituent->Position.T());
339 weightTime += TMath::Sqrt(constituent->Momentum.E());
340
341 candidate->AddCandidate(constituent);
342 }
343
344 avTime = time/weightTime;
345
346 candidate->Momentum = momentum;
347 candidate->Position.SetT(avTime);
348 candidate->Area.SetPxPyPzE(area.px(), area.py(), area.pz(), area.E());
349
350 candidate->DeltaEta = detaMax;
351 candidate->DeltaPhi = dphiMax;
352
353 // --- compute N-subjettiness with N = 1,2,3,4,5 ----
354
355 if(fComputeNsubjettiness)
356 {
357 Njettiness::AxesMode axisMode;
358
359 switch(fAxisMode)
360 {
361 default:
362 case 1:
363 axisMode = Njettiness::wta_kt_axes;
364 break;
365 case 2:
366 axisMode = Njettiness::onepass_wta_kt_axes;
367 break;
368 case 3:
369 axisMode = Njettiness::kt_axes;
370 break;
371 case 4:
372 axisMode = Njettiness::onepass_kt_axes;
373 break;
374 }
375
376 Njettiness::MeasureMode measureMode = Njettiness::unnormalized_measure;
377
378 Nsubjettiness nSub1(1, axisMode, measureMode, fBeta);
379 Nsubjettiness nSub2(2, axisMode, measureMode, fBeta);
380 Nsubjettiness nSub3(3, axisMode, measureMode, fBeta);
381 Nsubjettiness nSub4(4, axisMode, measureMode, fBeta);
382 Nsubjettiness nSub5(5, axisMode, measureMode, fBeta);
383
384 candidate->Tau[0] = nSub1(*itOutputList);
385 candidate->Tau[1] = nSub2(*itOutputList);
386 candidate->Tau[2] = nSub3(*itOutputList);
387 candidate->Tau[3] = nSub4(*itOutputList);
388 candidate->Tau[4] = nSub5(*itOutputList);
389 }
390
391
392 fOutputArray->Add(candidate);
393 }
394 delete sequence;
395}
Note: See TracBrowser for help on using the repository browser.