Fork me on GitHub

source: svn/trunk/modules/FastJetFinder.cc

Last change on this file was 1383, checked in by Pavel Demin, 11 years ago

fix namespace usage

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