Fork me on GitHub

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

Last change on this file since 1345 was 1345, checked in by Michele Selvaggi, 11 years ago

timing implemented

  • Property svn:keywords set to Id Revision Date
File size: 8.7 KB
Line 
1
2/** \class FastJetFinder
3 *
4 * Finds jets using FastJet library.
5 *
6 * $Date: 2013-12-21 14:00:11 +0000 (Sat, 21 Dec 2013) $
7 * $Revision: 1345 $
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
49using namespace std;
50using namespace fastjet;
51
52//------------------------------------------------------------------------------
53
54FastJetFinder::FastJetFinder() :
55 fPlugin(0), fDefinition(0), fAreaDefinition(0), fItInputArray(0)
56{
57
58}
59
60//------------------------------------------------------------------------------
61
62FastJetFinder::~FastJetFinder()
63{
64
65}
66
67//------------------------------------------------------------------------------
68
69void FastJetFinder::Init()
70{
71 JetDefinition::Plugin *plugin = NULL;
72
73 // read eta ranges
74
75 ExRootConfParam param = GetParam("RhoEtaRange");
76 Long_t i, size;
77
78 fEtaRangeMap.clear();
79 size = param.GetSize();
80 for(i = 0; i < size/2; ++i)
81 {
82 fEtaRangeMap[param[i*2].GetDouble()] = param[i*2 + 1].GetDouble();
83 }
84
85 // define algorithm
86
87 fJetAlgorithm = GetInt("JetAlgorithm", 6);
88 fParameterR = GetDouble("ParameterR", 0.5);
89
90 fConeRadius = GetDouble("ConeRadius", 0.5);
91 fSeedThreshold = GetDouble("SeedThreshold", 1.0);
92 fConeAreaFraction = GetDouble("ConeAreaFraction", 1.0);
93 fMaxIterations = GetInt("MaxIterations", 100);
94 fMaxPairSize = GetInt("MaxPairSize", 2);
95 fIratch = GetInt("Iratch", 1);
96 fAdjacencyCut = GetInt("AdjacencyCut", 2);
97 fOverlapThreshold = GetDouble("OverlapThreshold", 0.75);
98
99 fJetPTMin = GetDouble("JetPTMin", 10.0);
100
101 // --- Jet Area Parameters ---
102 fAreaAlgorithm = GetInt("AreaAlgorithm", 0);
103 fComputeRho = GetBool("ComputeRho", false);
104 // - ghost based areas -
105 fGhostEtaMax = GetDouble("GhostEtaMax", 5.0);
106 fRepeat = GetInt("Repeat", 1);
107 fGhostArea = GetDouble("GhostArea", 0.01);
108 fGridScatter = GetDouble("GridScatter", 1.0);
109 fPtScatter = GetDouble("PtScatter", 0.1);
110 fMeanGhostPt = GetDouble("MeanGhostPt", 1.0E-100);
111 // - voronoi based areas -
112 fEffectiveRfact = GetDouble("EffectiveRfact", 1.0);
113
114 switch(fAreaAlgorithm)
115 {
116 case 1:
117 fAreaDefinition = new fastjet::AreaDefinition(active_area_explicit_ghosts, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
118 break;
119 case 2:
120 fAreaDefinition = new fastjet::AreaDefinition(one_ghost_passive_area, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
121 break;
122 case 3:
123 fAreaDefinition = new fastjet::AreaDefinition(passive_area, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
124 break;
125 case 4:
126 fAreaDefinition = new fastjet::AreaDefinition(VoronoiAreaSpec(fEffectiveRfact));
127 break;
128 case 5:
129 fAreaDefinition = new fastjet::AreaDefinition(active_area, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
130 break;
131 default:
132 case 0:
133 fAreaDefinition = 0;
134 break;
135 }
136
137 switch(fJetAlgorithm)
138 {
139 case 1:
140 plugin = new fastjet::CDFJetCluPlugin(fSeedThreshold, fConeRadius, fAdjacencyCut, fMaxIterations, fIratch, fOverlapThreshold);
141 fDefinition = new fastjet::JetDefinition(plugin);
142 break;
143 case 2:
144 plugin = new fastjet::CDFMidPointPlugin(fSeedThreshold, fConeRadius, fConeAreaFraction, fMaxPairSize, fMaxIterations, fOverlapThreshold);
145 fDefinition = new fastjet::JetDefinition(plugin);
146 break;
147 case 3:
148 plugin = new fastjet::SISConePlugin(fConeRadius, fOverlapThreshold, fMaxIterations, fJetPTMin);
149 fDefinition = new fastjet::JetDefinition(plugin);
150 break;
151 case 4:
152 fDefinition = new fastjet::JetDefinition(fastjet::kt_algorithm, fParameterR);
153 break;
154 case 5:
155 fDefinition = new fastjet::JetDefinition(fastjet::cambridge_algorithm, fParameterR);
156 break;
157 default:
158 case 6:
159 fDefinition = new fastjet::JetDefinition(fastjet::antikt_algorithm, fParameterR);
160 break;
161 }
162
163 fPlugin = plugin;
164
165 ClusterSequence::print_banner();
166
167 // import input array
168
169 fInputArray = ImportArray(GetString("InputArray", "Calorimeter/towers"));
170 fItInputArray = fInputArray->MakeIterator();
171
172 // create output arrays
173
174 fOutputArray = ExportArray(GetString("OutputArray", "jets"));
175 fRhoOutputArray = ExportArray(GetString("RhoOutputArray", "rho"));
176}
177
178//------------------------------------------------------------------------------
179
180void FastJetFinder::Finish()
181{
182 if(fItInputArray) delete fItInputArray;
183 if(fDefinition) delete fDefinition;
184 if(fAreaDefinition) delete fAreaDefinition;
185 if(fPlugin) delete static_cast<JetDefinition::Plugin*>(fPlugin);
186}
187
188//------------------------------------------------------------------------------
189
190void FastJetFinder::Process()
191{
192 Candidate *candidate, *constituent;
193 TLorentzVector momentum;
194 Double_t deta, dphi, detaMax, dphiMax;
195 Double_t time, weightTime, avTime;
196 Int_t number;
197 Double_t rho = 0;
198 PseudoJet jet, area;
199 vector<PseudoJet> inputList, outputList;
200 ClusterSequence *sequence;
201 map< Double_t, Double_t >::iterator itEtaRangeMap;
202
203 DelphesFactory *factory = GetFactory();
204
205 inputList.clear();
206
207 // loop over input objects
208 fItInputArray->Reset();
209 number = 0;
210 while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
211 {
212 momentum = candidate->Momentum;
213 jet = PseudoJet(momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E());
214 jet.set_user_index(number);
215 inputList.push_back(jet);
216 ++number;
217 }
218
219 // construct jets
220 if(fAreaDefinition)
221 {
222 sequence = new ClusterSequenceArea(inputList, *fDefinition, *fAreaDefinition);
223 }
224 else
225 {
226 sequence = new ClusterSequence(inputList, *fDefinition);
227 }
228
229 // compute rho and store it
230 if(fComputeRho && fAreaDefinition)
231 {
232 for(itEtaRangeMap = fEtaRangeMap.begin(); itEtaRangeMap != fEtaRangeMap.end(); ++itEtaRangeMap)
233 {
234 Selector select_rapidity = SelectorAbsRapRange(itEtaRangeMap->first, itEtaRangeMap->second);
235 JetMedianBackgroundEstimator estimator(select_rapidity, *fDefinition, *fAreaDefinition);
236 estimator.set_particles(inputList);
237 rho = estimator.rho();
238
239 candidate = factory->NewCandidate();
240 candidate->Momentum.SetPtEtaPhiE(rho, 0.0, 0.0, rho);
241 candidate->Edges[0] = itEtaRangeMap->first;
242 candidate->Edges[1] = itEtaRangeMap->second;
243 fRhoOutputArray->Add(candidate);
244 }
245 }
246
247 outputList.clear();
248 outputList = sorted_by_pt(sequence->inclusive_jets(fJetPTMin));
249
250 // loop over all jets and export them
251 detaMax = 0.0;
252 dphiMax = 0.0;
253 vector<PseudoJet>::iterator itInputList, itOutputList;
254 for(itOutputList = outputList.begin(); itOutputList != outputList.end(); ++itOutputList)
255 {
256 momentum.SetPxPyPzE(itOutputList->px(), itOutputList->py(), itOutputList->pz(), itOutputList->E());
257 area.reset(0.0, 0.0, 0.0, 0.0);
258 if(fAreaDefinition) area = itOutputList->area_4vector();
259
260 candidate = factory->NewCandidate();
261
262 time=0;
263 weightTime=0;
264
265 inputList.clear();
266 inputList = sequence->constituents(*itOutputList);
267 for(itInputList = inputList.begin(); itInputList != inputList.end(); ++itInputList)
268 {
269 constituent = static_cast<Candidate*>(fInputArray->At(itInputList->user_index()));
270
271 deta = TMath::Abs(momentum.Eta() - constituent->Momentum.Eta());
272 dphi = TMath::Abs(momentum.DeltaPhi(constituent->Momentum));
273 if(deta > detaMax) detaMax = deta;
274 if(dphi > dphiMax) dphiMax = dphi;
275
276 time += TMath::Sqrt(constituent->Momentum.E())*(constituent->Position.T());
277 weightTime += TMath::Sqrt(constituent->Momentum.E());
278
279 candidate->AddCandidate(constituent);
280 }
281
282 avTime = time/weightTime;
283
284 candidate->Momentum = momentum;
285 candidate->Position.SetT(avTime);
286 candidate->Area.SetPxPyPzE(area.px(), area.py(), area.pz(), area.E());
287
288 candidate->DeltaEta = detaMax;
289 candidate->DeltaPhi = dphiMax;
290
291 fOutputArray->Add(candidate);
292 }
293 delete sequence;
294}
Note: See TracBrowser for help on using the repository browser.