Fork me on GitHub

source: git/modules/FastJetFinder.cc@ fb21bc8

ImprovedOutputFile Timing dual_readout llp
Last change on this file since fb21bc8 was 8336b6e, checked in by pavel <pavel@…>, 11 years ago

add eta ranges for rho

  • Property mode set to 100644
File size: 8.4 KB
Line 
1
2/** \class FastJetFinder
3 *
4 * Finds jets using FastJet library.
5 *
6 * $Date$
7 * $Revision$
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 = GetDouble("AdjacencyCut", 2.0);
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 Int_t number;
196 Double_t rho = 0;
197 PseudoJet jet, area;
198 vector<PseudoJet> inputList, outputList;
199 ClusterSequence *sequence;
200 map< Double_t, Double_t >::iterator itEtaRangeMap;
201
202 DelphesFactory *factory = GetFactory();
203
204 inputList.clear();
205
206 // loop over input objects
207 fItInputArray->Reset();
208 number = 0;
209 while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
210 {
211 momentum = candidate->Momentum;
212 jet = PseudoJet(momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E());
213 jet.set_user_index(number);
214 inputList.push_back(jet);
215 ++number;
216 }
217
218 // construct jets
219 if(fAreaDefinition)
220 {
221 sequence = new ClusterSequenceArea(inputList, *fDefinition, *fAreaDefinition);
222 }
223 else
224 {
225 sequence = new ClusterSequence(inputList, *fDefinition);
226 }
227
228 // compute rho and store it
229 if(fComputeRho && fAreaDefinition)
230 {
231 for(itEtaRangeMap = fEtaRangeMap.begin(); itEtaRangeMap != fEtaRangeMap.end(); ++itEtaRangeMap)
232 {
233 Selector select_rapidity = SelectorAbsRapRange(itEtaRangeMap->first, itEtaRangeMap->second);
234 JetMedianBackgroundEstimator estimator(select_rapidity, *fDefinition, *fAreaDefinition);
235 estimator.set_particles(inputList);
236 rho = estimator.rho();
237
238 candidate = factory->NewCandidate();
239 candidate->Momentum.SetPtEtaPhiE(rho, 0.0, 0.0, rho);
240 candidate->Edges[0] = itEtaRangeMap->first;
241 candidate->Edges[1] = itEtaRangeMap->second;
242 fRhoOutputArray->Add(candidate);
243 }
244 }
245
246 outputList.clear();
247 outputList = sorted_by_pt(sequence->inclusive_jets(fJetPTMin));
248
249 // loop over all jets and export them
250 detaMax = 0.0;
251 dphiMax = 0.0;
252 vector<PseudoJet>::iterator itInputList, itOutputList;
253 for(itOutputList = outputList.begin(); itOutputList != outputList.end(); ++itOutputList)
254 {
255 momentum.SetPxPyPzE(itOutputList->px(), itOutputList->py(), itOutputList->pz(), itOutputList->E());
256 area.reset(0.0, 0.0, 0.0, 0.0);
257 if(fAreaDefinition) area = itOutputList->area_4vector();
258
259 candidate = factory->NewCandidate();
260
261 inputList.clear();
262 inputList = sequence->constituents(*itOutputList);
263 for(itInputList = inputList.begin(); itInputList != inputList.end(); ++itInputList)
264 {
265 constituent = static_cast<Candidate*>(fInputArray->At(itInputList->user_index()));
266
267 deta = TMath::Abs(momentum.Eta() - constituent->Momentum.Eta());
268 dphi = TMath::Abs(momentum.DeltaPhi(constituent->Momentum));
269 if(deta > detaMax) detaMax = deta;
270 if(dphi > dphiMax) dphiMax = dphi;
271
272 candidate->AddCandidate(constituent);
273 }
274
275 candidate->Momentum = momentum;
276 candidate->Area.SetPxPyPzE(area.px(), area.py(), area.pz(), area.E());
277
278 candidate->DeltaEta = detaMax;
279 candidate->DeltaPhi = dphiMax;
280
281 fOutputArray->Add(candidate);
282 }
283 delete sequence;
284}
Note: See TracBrowser for help on using the repository browser.