Fork me on GitHub

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

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

add more AreaDefinition parameters for active_area

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