2 #define FASTJETUTIL_H 1
28 #include "FastJetProcessor.h"
29 #include "EClusterMode.h"
31 #include "LCIOSTLTypes.h"
32 #include "marlin/Processor.h"
35 #include <fastjet/PseudoJet.hh>
36 #include <fastjet/SISConePlugin.hh>
37 #include <fastjet/SISConeSphericalPlugin.hh>
38 #include <fastjet/CDFMidPointPlugin.hh>
39 #include <fastjet/CDFJetCluPlugin.hh>
40 #include <fastjet/NestedDefsPlugin.hh>
41 #include <fastjet/EECambridgePlugin.hh>
42 #include <fastjet/JadePlugin.hh>
44 #include <fastjet/contrib/ValenciaPlugin.hh>
49 #define ITERATIVE_INCLUSIVE_MAX_ITERATIONS 20
50 typedef std::vector< fastjet::PseudoJet > PseudoJetList;
70 _jetAlgoNameAndParams( EVENT::StringVec() ),
74 _clusterModeNameAndParam( EVENT::StringVec() ),
77 _jetRecoSchemeName(
""),
81 _requestedNumberOfJets(0),
89 _cs(
new fastjet::ClusterSequence(*(rhs._cs))),
90 _jetAlgoNameAndParams( rhs._jetAlgoNameAndParams ),
91 _jetAlgoName(rhs._jetAlgoName),
92 _jetAlgo(
new fastjet::JetDefinition(*(rhs._jetAlgo))),
93 _jetAlgoType(rhs._jetAlgoType),
94 _clusterModeNameAndParam(rhs._clusterModeNameAndParam),
95 _clusterModeName(rhs._clusterModeName),
96 _clusterMode(rhs._clusterMode),
97 _jetRecoSchemeName(rhs._jetRecoSchemeName),
98 _jetRecoScheme(rhs._jetRecoScheme),
99 _strategyName(rhs._strategyName),
100 _strategy(rhs._strategy),
101 _requestedNumberOfJets(rhs._requestedNumberOfJets),
112 this->_cs =
new fastjet::ClusterSequence(*(rhs._cs));
113 delete this->_jetAlgo;
114 this->_jetAlgo =
new fastjet::JetDefinition(*rhs._jetAlgo);
126 fastjet::ClusterSequence *_cs;
129 EVENT::StringVec _jetAlgoNameAndParams;
130 std::string _jetAlgoName;
131 fastjet::JetDefinition* _jetAlgo;
132 fastjet::JetAlgorithm _jetAlgoType;
135 EVENT::StringVec _clusterModeNameAndParam;
136 std::string _clusterModeName;
137 EClusterMode _clusterMode;
140 std::string _jetRecoSchemeName;
141 fastjet::RecombinationScheme _jetRecoScheme;
144 std::string _strategyName;
145 fastjet::Strategy _strategy;
148 unsigned _requestedNumberOfJets;
162 inline PseudoJetList
clusterJets(PseudoJetList& pjList, LCCollection* reconstructedPars);
170 inline bool isJetAlgo(std::string algo,
int nrParams,
int supportedModes);
173 PseudoJetList doIterativeInclusiveClustering(PseudoJetList& pjList);
181 EVENT::StringVec defAlgoAndParam;
182 defAlgoAndParam.push_back(
"kt_algorithm");
183 defAlgoAndParam.push_back(
"0.7");
184 proc->registerProcessorParameter(
186 "Selects the algorithm and its parameters. E.g. 'kt_algorithm 0.7' or 'ee_kt_algorithm'. For a full list of supported algorithms, see the logfile after execution.",
187 _jetAlgoNameAndParams,
190 proc->registerProcessorParameter(
191 "recombinationScheme",
192 "The recombination scheme used when merging 2 particles. Usually there is no need to use anything else than 4-Vector addition: E_scheme",
194 std::string(
"E_scheme"));
196 EVENT::StringVec defClusterMode;
197 defClusterMode.push_back(
"Inclusive");
198 proc->registerProcessorParameter(
200 "One of 'Inclusive <minPt>', 'InclusiveIterativeNJets <nrJets> <minE>', 'ExclusiveNJets <nrJets>', 'ExclusiveYCut <yCut>'. Note: not all modes are available for all algorithms.",
201 _clusterModeNameAndParam,
216 void FastJetUtil::initJetAlgo() {
221 if (_jetAlgoNameAndParams.size() < 1)
222 throw Exception(
"No Jet algorithm provided!");
225 this->_jetAlgoName = _jetAlgoNameAndParams[0];
230 streamlog_out(MESSAGE) <<
"Algorithms: ";
233 if (isJetAlgo(
"kt_algorithm", 1, FJ_inclusive | FJ_exclusive_nJets | FJ_exclusive_yCut | OWN_inclusiveIteration)) {
234 _jetAlgoType = fastjet::kt_algorithm;
235 _jetAlgo =
new fastjet::JetDefinition(
236 _jetAlgoType, atof(_jetAlgoNameAndParams[1].c_str()), _jetRecoScheme, _strategy);
240 if (isJetAlgo(
"cambridge_algorithm", 1, FJ_inclusive | FJ_exclusive_nJets | FJ_exclusive_yCut | OWN_inclusiveIteration)) {
241 _jetAlgoType = fastjet::cambridge_algorithm;
242 _jetAlgo =
new fastjet::JetDefinition(
243 _jetAlgoType, atof(_jetAlgoNameAndParams[1].c_str()), _jetRecoScheme, _strategy);
246 if (isJetAlgo(
"antikt_algorithm", 1, FJ_inclusive | OWN_inclusiveIteration)) {
247 _jetAlgoType = fastjet::antikt_algorithm;
248 _jetAlgo =
new fastjet::JetDefinition(
249 _jetAlgoType, atof(_jetAlgoNameAndParams[1].c_str()), _jetRecoScheme, _strategy);
252 if (isJetAlgo(
"genkt_algorithm", 2, FJ_inclusive | OWN_inclusiveIteration | FJ_exclusive_nJets | FJ_exclusive_yCut)) {
253 _jetAlgoType = fastjet::genkt_algorithm;
254 _jetAlgo =
new fastjet::JetDefinition(
255 _jetAlgoType, atof(_jetAlgoNameAndParams[1].c_str()), atof(_jetAlgoNameAndParams[2].c_str()), _jetRecoScheme, _strategy);
258 if (isJetAlgo(
"cambridge_for_passive_algorithm", 1, FJ_inclusive | OWN_inclusiveIteration | FJ_exclusive_nJets | FJ_exclusive_yCut)) {
259 _jetAlgoType = fastjet::cambridge_for_passive_algorithm;
260 _jetAlgo =
new fastjet::JetDefinition(
261 _jetAlgoType, atof(_jetAlgoNameAndParams[1].c_str()), _jetRecoScheme, _strategy);
264 if (isJetAlgo(
"genkt_for_passive_algorithm", 1, FJ_inclusive | OWN_inclusiveIteration)) {
265 _jetAlgoType = fastjet::genkt_for_passive_algorithm;
266 _jetAlgo =
new fastjet::JetDefinition(
267 _jetAlgoType, atof(_jetAlgoNameAndParams[1].c_str()), _jetRecoScheme, _strategy);
270 if (isJetAlgo(
"ee_kt_algorithm", 0, FJ_exclusive_nJets | FJ_exclusive_yCut)) {
271 _jetAlgoType = fastjet::ee_kt_algorithm;
272 _jetAlgo =
new fastjet::JetDefinition(
273 _jetAlgoType, _jetRecoScheme, _strategy);
276 if (isJetAlgo(
"ee_genkt_algorithm", 1, FJ_exclusive_nJets | FJ_exclusive_yCut)) {
277 _jetAlgoType = fastjet::ee_genkt_algorithm;
278 _jetAlgo =
new fastjet::JetDefinition(
279 _jetAlgoType, atof(_jetAlgoNameAndParams[1].c_str()), _jetRecoScheme, _strategy);
282 if (isJetAlgo(
"SISConePlugin", 2, FJ_inclusive | OWN_inclusiveIteration)) {
283 fastjet::SISConePlugin* pl;
285 pl =
new fastjet::SISConePlugin(
286 atof(_jetAlgoNameAndParams[1].c_str()),
287 atof(_jetAlgoNameAndParams[2].c_str())
289 _jetAlgo =
new fastjet::JetDefinition(pl);
292 if (isJetAlgo(
"SISConeSphericalPlugin", 2, FJ_inclusive | OWN_inclusiveIteration)) {
294 fastjet::SISConeSphericalPlugin* pl;
295 pl =
new fastjet::SISConeSphericalPlugin(
296 atof(_jetAlgoNameAndParams[1].c_str()),
297 atof(_jetAlgoNameAndParams[2].c_str())
300 _jetAlgo =
new fastjet::JetDefinition(pl);
303 if (isJetAlgo(
"ValenciaPlugin", 3, FJ_exclusive_nJets | FJ_exclusive_yCut)) {
305 fastjet::contrib::ValenciaPlugin* pl;
306 pl =
new fastjet::contrib::ValenciaPlugin(
307 atof(_jetAlgoNameAndParams[1].c_str()),
308 atof(_jetAlgoNameAndParams[2].c_str()),
309 atof(_jetAlgoNameAndParams[3].c_str())
312 _jetAlgo =
new fastjet::JetDefinition(pl);
336 streamlog_out(MESSAGE) << std::endl;
339 streamlog_out(ERROR) <<
"The given algorithm \"" << _jetAlgoName
340 <<
"\" is unknown to me!" << std::endl;
341 throw Exception(
"Unknown FastJet algorithm.");
344 streamlog_out(MESSAGE) <<
"jet algorithm: " << _jetAlgo->description() << std::endl;
348 bool FastJetUtil::isJetAlgo(std::string algo,
int nrParams,
int supportedModes)
350 streamlog_out(MESSAGE) <<
" " << algo;
353 if (_jetAlgoName.compare(algo) != 0) {
356 streamlog_out(MESSAGE) <<
"*";
361 if ((
int)_jetAlgoNameAndParams.size() - 1 != nrParams) {
362 streamlog_out(ERROR) << std::endl
363 <<
"Wrong numbers of parameters for algorithm: " << algo << std::endl
364 <<
"We need " << nrParams <<
" params, but we got " << _jetAlgoNameAndParams.size() - 1 << std::endl;
365 throw Exception(
"You have insufficient number of parameters for this algorithm! See log for more details.");
369 if ((supportedModes & _clusterMode) != _clusterMode) {
370 streamlog_out(ERROR) << std::endl
371 <<
"This algorithm is not capable of running in this clustering mode ("
372 << _clusterMode <<
"). Sorry!" << std::endl;
373 throw Exception(
"This algorithm is not capable of running in this mode");
382 if (_jetRecoSchemeName.compare(
"E_scheme") == 0)
383 _jetRecoScheme = fastjet::E_scheme;
384 else if (_jetRecoSchemeName.compare(
"pt_scheme") == 0)
385 _jetRecoScheme = fastjet::pt_scheme;
386 else if (_jetRecoSchemeName.compare(
"pt2_scheme") == 0)
387 _jetRecoScheme = fastjet::pt2_scheme;
388 else if (_jetRecoSchemeName.compare(
"Et_scheme") == 0)
389 _jetRecoScheme = fastjet::Et_scheme;
390 else if (_jetRecoSchemeName.compare(
"Et2_scheme") == 0)
391 _jetRecoScheme = fastjet::Et2_scheme;
392 else if (_jetRecoSchemeName.compare(
"BIpt_scheme") == 0)
393 _jetRecoScheme = fastjet::BIpt_scheme;
394 else if (_jetRecoSchemeName.compare(
"BIpt2_scheme") == 0)
395 _jetRecoScheme = fastjet::BIpt2_scheme;
398 <<
"Unknown recombination scheme: " << _jetRecoSchemeName << std::endl;
399 throw Exception(
"Unknown FastJet recombination scheme! See log for more details.");
402 streamlog_out(MESSAGE) <<
"recombination scheme: " << _jetRecoSchemeName << std::endl;
405 void FastJetUtil::initStrategy() {
409 _strategy = fastjet::Best;
410 _strategyName =
"Best";
411 streamlog_out(MESSAGE) <<
"Strategy: " << _strategyName << std::endl;
418 if (_clusterModeNameAndParam.size() == 0)
419 throw Exception(
"Cluster mode not specified");
422 _clusterModeName = _clusterModeNameAndParam[0];
427 if (_clusterModeName.compare(
"Inclusive") == 0) {
429 if (_clusterModeNameAndParam.size() != 2) {
430 throw Exception(
"Wrong Parameter(s) for Clustering Mode. Expected:\n <parameter name=\"clusteringMode\" type=\"StringVec\"> Inclusive <minPt> </parameter>");
433 _minPt = atof(_clusterModeNameAndParam[1].c_str());
434 _clusterMode = FJ_inclusive;
436 }
else if (_clusterModeName.compare(
"InclusiveIterativeNJets") == 0) {
438 if (_clusterModeNameAndParam.size() != 3) {
439 throw Exception(
"Wrong Parameter(s) for Clustering Mode. Expected:\n <parameter name=\"clusteringMode\" type=\"StringVec\"> InclusiveIterativeNJets <NJets> <minE> </parameter>");
442 _requestedNumberOfJets = atoi(_clusterModeNameAndParam[1].c_str());
443 _minE = atoi(_clusterModeNameAndParam[2].c_str());
444 _clusterMode = OWN_inclusiveIteration;
446 }
else if (_clusterModeName.compare(
"ExclusiveNJets") == 0) {
448 if (_clusterModeNameAndParam.size() != 2) {
449 throw Exception(
"Wrong Parameter(s) for Clustering Mode. Expected:\n <parameter name=\"clusteringMode\" type=\"StringVec\"> ExclusiveNJets <NJets> </parameter>");
452 _requestedNumberOfJets = atoi(_clusterModeNameAndParam[1].c_str());
453 _clusterMode = FJ_exclusive_nJets;
455 }
else if (_clusterModeName.compare(
"ExclusiveYCut") == 0) {
457 if (_clusterModeNameAndParam.size() != 2) {
458 throw Exception(
"Wrong Parameter(s) for Clustering Mode. Expected:\n <parameter name=\"clusteringMode\" type=\"StringVec\"> ExclusiveYCut <yCut> </parameter>");
461 _yCut = atof(_clusterModeNameAndParam[1].c_str());
462 _clusterMode = FJ_exclusive_yCut;
465 throw Exception(
"Unknown cluster mode.");
468 streamlog_out(MESSAGE) <<
"cluster mode: " << _clusterMode << std::endl;
474 _cs =
new fastjet::ClusterSequence(pjList, *_jetAlgo);
478 if (_clusterMode == FJ_inclusive) {
480 jets = _cs->inclusive_jets(_minPt);
482 }
else if (_clusterMode == FJ_exclusive_yCut) {
484 jets = _cs->exclusive_jets_ycut(_yCut);
486 }
else if (_clusterMode == FJ_exclusive_nJets) {
489 if (reconstructedPars->getNumberOfElements() < (int)_requestedNumberOfJets) {
491 streamlog_out(WARNING) <<
"Not enough elements in the input collection to create " << _requestedNumberOfJets <<
" jets." << std::endl;
496 jets = _cs->exclusive_jets((
int)(_requestedNumberOfJets));
500 }
else if (_clusterMode == OWN_inclusiveIteration) {
503 if (reconstructedPars->getNumberOfElements() < (int)_requestedNumberOfJets) {
505 streamlog_out(WARNING) <<
"Not enough elements in the input collection to create " << _requestedNumberOfJets <<
" jets." << std::endl;
510 jets = doIterativeInclusiveClustering(pjList);
523 PseudoJetList pjList;
524 for (
int i = 0; i < recCol->getNumberOfElements(); ++i) {
525 ReconstructedParticle* par =
static_cast<ReconstructedParticle*
> (recCol->getElementAt(i));
526 pjList.push_back( fastjet::PseudoJet( par->getMomentum()[0],
527 par->getMomentum()[1],
528 par->getMomentum()[2],
529 par->getEnergy() ) );
530 pjList.back().set_user_index(i);
536 PseudoJetList FastJetUtil::doIterativeInclusiveClustering( PseudoJetList& pjList) {
542 double RDiff = R / 2;
544 PseudoJetList jetsReturn;
553 fastjet::SISConePlugin* pluginSisCone = NULL;
554 fastjet::SISConeSphericalPlugin* pluginSisConeSph = NULL;
556 bool useSisCone = _jetAlgoName.compare(
"SISConePlugin") == 0;
557 bool useSisConeSph = _jetAlgoName.compare(
"SISConeSphericalPlugin") == 0;
559 double sisConeOverlapThreshold = 0;
560 if (useSisCone || useSisConeSph)
561 sisConeOverlapThreshold = atof(_jetAlgoNameAndParams[2].c_str());
572 for (iIter=0; iIter<ITERATIVE_INCLUSIVE_MAX_ITERATIONS; iIter++) {
575 fastjet::JetDefinition* jetDefinition = NULL;
579 pluginSisCone =
new fastjet::SISConePlugin(R, sisConeOverlapThreshold);
580 jetDefinition =
new fastjet::JetDefinition(pluginSisCone);
581 }
else if (useSisConeSph) {
582 pluginSisConeSph =
new fastjet::SISConeSphericalPlugin(R, sisConeOverlapThreshold);
583 jetDefinition =
new fastjet::JetDefinition(pluginSisConeSph);
585 jetDefinition =
new fastjet::JetDefinition(_jetAlgoType, R, _jetRecoScheme, _strategy);
589 fastjet::ClusterSequence cs(pjList, *jetDefinition);
591 jets = cs.inclusive_jets(0);
596 for (
unsigned j=0; j<jets.size(); j++)
597 if (jets[j].E() > _minE)
598 jetsReturn.push_back(jets[j]);
599 nJets = jetsReturn.size();
601 streamlog_out(DEBUG) << iIter <<
" " << R <<
" " << jets.size() <<
" " << nJets << std::endl;
603 if (nJets == _requestedNumberOfJets) {
604 delete pluginSisCone; pluginSisCone = NULL;
605 delete pluginSisConeSph; pluginSisConeSph = NULL;
606 delete jetDefinition;
609 }
else if (nJets < _requestedNumberOfJets) {
614 }
else if (nJets > _requestedNumberOfJets) {
622 delete pluginSisCone; pluginSisCone = NULL;
623 delete pluginSisConeSph; pluginSisConeSph = NULL;
624 delete jetDefinition;
628 if (iIter == ITERATIVE_INCLUSIVE_MAX_ITERATIONS) {
629 streamlog_out(WARNING) <<
"Maximum number of iterations reached. Canceling" << std::endl;
639 #endif // FastJetUtil_h
void registerFastJetParameters(T *proc)
call in processor constructor (c'tor) to register parameters
Definition: FastJetUtil.h:179
void init()
call in processor init
Definition: FastJetUtil.h:207
PseudoJetList convertFromRecParticle(LCCollection *recCol)
convert reconstructed particles to pseudo jets
Definition: FastJetUtil.h:521
Definition: FastJetUtil.h:59
void initRecoScheme()
parse the steering parameter for the recombination scheme
Definition: FastJetUtil.h:380
Definition: FastJetUtil.h:66
Definition: FastJetProcessor.h:33
Definition: FastJetUtil.h:54
PseudoJetList clusterJets(PseudoJetList &pjList, LCCollection *reconstructedPars)
does the actual clustering
Definition: FastJetUtil.h:471
void initClusterMode()
parse the clustermode string
Definition: FastJetUtil.h:415