Fork me on GitHub

source: svn/trunk/modules/FastJetFinder.cc@ 1373

Last change on this file since 1373 was 1372, checked in by Pavel Demin, 11 years ago

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