MarlinFastJet  0.3.0
FastJetUtil.h
1 #ifndef FASTJETUTIL_H
2 #define FASTJETUTIL_H 1
3 /*
4  * In order to use in a processor
5  * do forward declaration to FastJetUtil in your MyJetProcessor.h
6  * class FastJetUtil;
7  *
8  * make FastJetUtil* a member object
9  * make FastJetUtil a friend class of your Processor (registerProcessorParameter is protected so we need to be a friend)
10  * friend class FastJetUtil;
11  *
12  * in the constructor call
13  * registerFastJetParameters( this )
14  *
15  * in init call
16  * init()
17 
18  * in procesEvent
19  * call:
20  * PseudoJetList convertFromRecParticle(LCCollection* recCol);
21  * inline PseudoJetList clusterJets(PseudoJetList& pjList, LCCollection* reconstructedPars);
22  *
23  */
24 
25 
26 class FastJetProcessor;
27 
28 #include "FastJetProcessor.h"
29 #include "EClusterMode.h"
30 
31 #include "LCIOSTLTypes.h"
32 #include "marlin/Processor.h"
33 
34 //FastJet
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>
43 
44 #include <fastjet/contrib/ValenciaPlugin.hh>
45 
46 #include <stdexcept>
47 #include <string>
48 
49 #define ITERATIVE_INCLUSIVE_MAX_ITERATIONS 20
50 typedef std::vector< fastjet::PseudoJet > PseudoJetList;
51 
52 
53 
54 class SkippedFixedNrJetException: public std::runtime_error {
55 public:
56  SkippedFixedNrJetException():std::runtime_error("") {}
57 };
58 
59 class SkippedMaxIterationException: public std::runtime_error {
60 public:
61  SkippedMaxIterationException( PseudoJetList& jets) :std::runtime_error(""), _jets(jets) {}
62  PseudoJetList& _jets;
63 };
64 
65 
66 class FastJetUtil {
67 
68 public:
69  FastJetUtil(): _cs(NULL),
70  _jetAlgoNameAndParams( EVENT::StringVec() ),
71  _jetAlgoName(""),
72  _jetAlgo(NULL),
73  _jetAlgoType(),
74  _clusterModeNameAndParam( EVENT::StringVec() ),
75  _clusterModeName(""),
76  _clusterMode( NONE ),
77  _jetRecoSchemeName(""),
78  _jetRecoScheme(),
79  _strategyName(""),
80  _strategy(),
81  _requestedNumberOfJets(0),
82  _yCut(0.0),
83  _minPt(0.0),
84  _minE(0.0)
85  {}
86 
87 
88  FastJetUtil(const FastJetUtil& rhs):
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),
102  _yCut(rhs._yCut),
103  _minPt(rhs._minPt),
104  _minE(rhs._minE)
105  {}
106 
107  FastJetUtil& operator=(const FastJetUtil& rhs) {
108  if( this == &rhs ){
109  return *this;
110  }
111  delete this->_cs;
112  this->_cs = new fastjet::ClusterSequence(*(rhs._cs));
113  delete this->_jetAlgo;
114  this->_jetAlgo = new fastjet::JetDefinition(*rhs._jetAlgo);
115  return *this;
116  }
117 
118 
119  ~FastJetUtil() {
120  delete _cs;
121  delete _jetAlgo;
122  }
123 
124 public:
125 
126  fastjet::ClusterSequence *_cs;
127 
128  // jet algorithm
129  EVENT::StringVec _jetAlgoNameAndParams;
130  std::string _jetAlgoName;
131  fastjet::JetDefinition* _jetAlgo;
132  fastjet::JetAlgorithm _jetAlgoType;
133 
134  // clustering mode
135  EVENT::StringVec _clusterModeNameAndParam;
136  std::string _clusterModeName;
137  EClusterMode _clusterMode;
138 
139  // jet reco scheme
140  std::string _jetRecoSchemeName;
141  fastjet::RecombinationScheme _jetRecoScheme;
142 
143  // jet strategy
144  std::string _strategyName;
145  fastjet::Strategy _strategy;
146 
147  // parameters
148  unsigned _requestedNumberOfJets;
149  double _yCut;
150  double _minPt;
151  double _minE;
152 
153 public:
155  template< class T>
156  inline void registerFastJetParameters(T* proc);
158  inline void init();
160  inline PseudoJetList convertFromRecParticle(LCCollection* recCol);
162  inline PseudoJetList clusterJets(PseudoJetList& pjList, LCCollection* reconstructedPars);
163 
164 protected:
165  // helper functions to init the jet algorithms in general
166  void initJetAlgo();
167  void initRecoScheme();
168  void initStrategy();
169  void initClusterMode();
170  inline bool isJetAlgo(std::string algo, int nrParams, int supportedModes);
171 
172  // special clustering function, called from clusterJets
173  PseudoJetList doIterativeInclusiveClustering(PseudoJetList& pjList);
174 
175 }; //end class FastJetUtil
176 
177 
178 template< class T >
180 
181  EVENT::StringVec defAlgoAndParam;
182  defAlgoAndParam.push_back("kt_algorithm");
183  defAlgoAndParam.push_back("0.7");
184  proc->registerProcessorParameter(
185  "algorithm",
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,
188  defAlgoAndParam);
189 
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",
193  _jetRecoSchemeName,
194  std::string("E_scheme"));
195 
196  EVENT::StringVec defClusterMode;
197  defClusterMode.push_back("Inclusive");
198  proc->registerProcessorParameter(
199  "clusteringMode",
200  "One of 'Inclusive <minPt>', 'InclusiveIterativeNJets <nrJets> <minE>', 'ExclusiveNJets <nrJets>', 'ExclusiveYCut <yCut>'. Note: not all modes are available for all algorithms.",
201  _clusterModeNameAndParam,
202  defClusterMode);
203 
204 
205 }
206 
208 
209  initStrategy();
210  initRecoScheme();
211  initClusterMode();
212  initJetAlgo();
213 
214 }
215 
216 void FastJetUtil::initJetAlgo() {
217 
218  // parse the given jet algorithm string
219 
220  // sanity check
221  if (_jetAlgoNameAndParams.size() < 1)
222  throw Exception("No Jet algorithm provided!");
223 
224  // save the name
225  this->_jetAlgoName = _jetAlgoNameAndParams[0];
226 
227  // check all supported algorithms and create the appropriate FJ instance
228 
229  _jetAlgo = NULL;
230  streamlog_out(MESSAGE) << "Algorithms: "; // the isJetAlgo function will write to streamlog_out(MESSAGE), so that we get a list of available algorithms in the log
231 
232  // example: kt_algorithm, needs 1 parameter, supports inclusive, inclusiveIterative, exlusiveNJets and exlusiveYCut clustering
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);
237 
238  }
239 
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);
244  }
245 
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);
250  }
251 
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);
256  }
257 
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);
262  }
263 
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);
268  }
269 
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);
274  }
275 
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);
280  }
281 
282  if (isJetAlgo("SISConePlugin", 2, FJ_inclusive | OWN_inclusiveIteration)) {
283  fastjet::SISConePlugin* pl;
284 
285  pl = new fastjet::SISConePlugin(
286  atof(_jetAlgoNameAndParams[1].c_str()),
287  atof(_jetAlgoNameAndParams[2].c_str())
288  );
289  _jetAlgo = new fastjet::JetDefinition(pl);
290  }
291 
292  if (isJetAlgo("SISConeSphericalPlugin", 2, FJ_inclusive | OWN_inclusiveIteration)) {
293 
294  fastjet::SISConeSphericalPlugin* pl;
295  pl = new fastjet::SISConeSphericalPlugin(
296  atof(_jetAlgoNameAndParams[1].c_str()),
297  atof(_jetAlgoNameAndParams[2].c_str())
298  );
299 
300  _jetAlgo = new fastjet::JetDefinition(pl);
301  }
302 
303  if (isJetAlgo("ValenciaPlugin", 3, FJ_exclusive_nJets | FJ_exclusive_yCut)) {
304 
305  fastjet::contrib::ValenciaPlugin* pl;
306  pl = new fastjet::contrib::ValenciaPlugin(
307  atof(_jetAlgoNameAndParams[1].c_str()), // R value
308  atof(_jetAlgoNameAndParams[2].c_str()), // beta value
309  atof(_jetAlgoNameAndParams[3].c_str()) // gamma value
310  );
311 
312  _jetAlgo = new fastjet::JetDefinition(pl);
313  }
314 
315  // TODO: Maybe we should complete the user defined (ATLAS, CMS, ..) algorithms
316  // if (isJetAlgo("ATLASConePlugin", 1))
317  // {
318  // fastjet::ATLASConePlugin* pl;
319  // pl = new fastjet::ATLASConePlugin(
320  // atof(_jetAlgoAndParamsString[1].c_str())
321  // );
322  //
323  // _jetAlgo = new fastjet::JetDefinition(pl);
324  // }
325  //
326  // if (isJetAlgo("CMSIterativeConePlugin", 1))
327  // {
328  // fastjet::CMSIterativeConePlugin* pl;
329  // pl = new fastjet::CMSIterativeConePlugin(
330  // atof(_jetAlgoAndParamsString[1].c_str())
331  // );
332  //
333  // _jetAlgo = new fastjet::JetDefinition(pl);
334  // }
335 
336  streamlog_out(MESSAGE) << std::endl; // end of list of available algorithms
337 
338  if (!_jetAlgo) {
339  streamlog_out(ERROR) << "The given algorithm \"" << _jetAlgoName
340  << "\" is unknown to me!" << std::endl;
341  throw Exception("Unknown FastJet algorithm.");
342  }
343 
344  streamlog_out(MESSAGE) << "jet algorithm: " << _jetAlgo->description() << std::endl;
345 }
346 
347 
348 bool FastJetUtil::isJetAlgo(std::string algo, int nrParams, int supportedModes)
349 {
350  streamlog_out(MESSAGE) << " " << algo;
351 
352  // check if the chosen algorithm is the same as it was passed
353  if (_jetAlgoName.compare(algo) != 0) {
354  return false;
355  }
356  streamlog_out(MESSAGE) << "*"; // mark the current algorithm as the selected
357  // one, even before we did our checks on nr of
358  // parameters etc.
359 
360  // check if we have enough number of parameters
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.");
366  }
367 
368  // check if the mode is supported via a binary AND
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");
374  }
375 
376  return true;
377 }
378 
381 
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;
396  else {
397  streamlog_out(ERROR)
398  << "Unknown recombination scheme: " << _jetRecoSchemeName << std::endl;
399  throw Exception("Unknown FastJet recombination scheme! See log for more details.");
400  }
401 
402  streamlog_out(MESSAGE) << "recombination scheme: " << _jetRecoSchemeName << std::endl;
403 }
404 
405 void FastJetUtil::initStrategy() {
406  // we could provide a steering parameter for this and it would be parsed here.
407  // however, when using the 'Best' clustering strategy FJ will chose automatically from a big list
408  // changing this would (most likely) only change the speed of calculation, not the outcome
409  _strategy = fastjet::Best;
410  _strategyName = "Best";
411  streamlog_out(MESSAGE) << "Strategy: " << _strategyName << std::endl;
412 }
413 
416 
417  // at least a name has to be given
418  if (_clusterModeNameAndParam.size() == 0)
419  throw Exception("Cluster mode not specified");
420 
421  // save the name of the cluster mode
422  _clusterModeName = _clusterModeNameAndParam[0];
423 
424  _clusterMode = NONE;
425 
426  // check the different cluster mode possibilities, and check if the number of parameters are correct
427  if (_clusterModeName.compare("Inclusive") == 0) {
428 
429  if (_clusterModeNameAndParam.size() != 2) {
430  throw Exception("Wrong Parameter(s) for Clustering Mode. Expected:\n <parameter name=\"clusteringMode\" type=\"StringVec\"> Inclusive <minPt> </parameter>");
431  }
432 
433  _minPt = atof(_clusterModeNameAndParam[1].c_str());
434  _clusterMode = FJ_inclusive;
435 
436  } else if (_clusterModeName.compare("InclusiveIterativeNJets") == 0) {
437 
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>");
440  }
441 
442  _requestedNumberOfJets = atoi(_clusterModeNameAndParam[1].c_str());
443  _minE = atoi(_clusterModeNameAndParam[2].c_str());
444  _clusterMode = OWN_inclusiveIteration;
445 
446  } else if (_clusterModeName.compare("ExclusiveNJets") == 0) {
447 
448  if (_clusterModeNameAndParam.size() != 2) {
449  throw Exception("Wrong Parameter(s) for Clustering Mode. Expected:\n <parameter name=\"clusteringMode\" type=\"StringVec\"> ExclusiveNJets <NJets> </parameter>");
450  }
451 
452  _requestedNumberOfJets = atoi(_clusterModeNameAndParam[1].c_str());
453  _clusterMode = FJ_exclusive_nJets;
454 
455  } else if (_clusterModeName.compare("ExclusiveYCut") == 0) {
456 
457  if (_clusterModeNameAndParam.size() != 2) {
458  throw Exception("Wrong Parameter(s) for Clustering Mode. Expected:\n <parameter name=\"clusteringMode\" type=\"StringVec\"> ExclusiveYCut <yCut> </parameter>");
459  }
460 
461  _yCut = atof(_clusterModeNameAndParam[1].c_str());
462  _clusterMode = FJ_exclusive_yCut;
463 
464  } else {
465  throw Exception("Unknown cluster mode.");
466  }
467 
468  streamlog_out(MESSAGE) << "cluster mode: " << _clusterMode << std::endl;
469 }
470 
471 PseudoJetList FastJetUtil::clusterJets( PseudoJetList& pjList, LCCollection* reconstructedPars) {
473  // do the jet finding for the user defined parameter jet finder
474  _cs = new fastjet::ClusterSequence(pjList, *_jetAlgo);
475 
476  PseudoJetList jets; // will contain the found jets
477 
478  if (_clusterMode == FJ_inclusive) {
479 
480  jets = _cs->inclusive_jets(_minPt);
481 
482  } else if (_clusterMode == FJ_exclusive_yCut) {
483 
484  jets = _cs->exclusive_jets_ycut(_yCut);
485 
486  } else if (_clusterMode == FJ_exclusive_nJets) {
487 
488  // sanity check: if we have not enough particles, FJ will cause an assert
489  if (reconstructedPars->getNumberOfElements() < (int)_requestedNumberOfJets) {
490 
491  streamlog_out(WARNING) << "Not enough elements in the input collection to create " << _requestedNumberOfJets << " jets." << std::endl;
493 
494  } else {
495 
496  jets = _cs->exclusive_jets((int)(_requestedNumberOfJets));
497 
498  }
499 
500  } else if (_clusterMode == OWN_inclusiveIteration) {
501 
502  // sanity check: if we have not enough particles, FJ will cause an assert
503  if (reconstructedPars->getNumberOfElements() < (int)_requestedNumberOfJets) {
504 
505  streamlog_out(WARNING) << "Not enough elements in the input collection to create " << _requestedNumberOfJets << " jets." << std::endl;
507 
508  } else {
509 
510  jets = doIterativeInclusiveClustering(pjList);
511 
512  }
513 
514  }
515 
516  return jets;
517 
518 }
519 
520 
521 PseudoJetList FastJetUtil::convertFromRecParticle(LCCollection* recCol) {
522  // foreach RecoParticle in the LCCollection: convert it into a PseudoJet and and save it in our list
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); // save the id of this recParticle
531  }
532 
533  return pjList;
534 }
535 
536 PseudoJetList FastJetUtil::doIterativeInclusiveClustering( PseudoJetList& pjList) {
537  // lets do a iterative procedure until we found the correct number of jets
538  // for that we will do inclusive clustering, modifying the R parameter in some kind of minimization
539  // this is based on Marco Battaglia's FastJetClustering
540 
541  double R = M_PI_4; // maximum of R is Pi/2, minimum is 0. So we start hat Pi/4
542  double RDiff = R / 2; // the step size we modify the R parameter at each iteration. Its size for the n-th step is R/(2n), i.e. starts with R/2
543  PseudoJetList jets;
544  PseudoJetList jetsReturn;
545  unsigned nJets;
546  int iIter = 0; // nr of current iteration
547 
548  // these variables are only used if the SisCone(Spherical)Plugin is selected
549  // This is necessary, as these are plugins and hence use a different constructor than
550  // the built in fastjet algorithms
551 
552  // here we save pointer to the plugins, so that if they are created, we can delete them again after usage
553  fastjet::SISConePlugin* pluginSisCone = NULL;
554  fastjet::SISConeSphericalPlugin* pluginSisConeSph = NULL;
555  // check if we use the siscones
556  bool useSisCone = _jetAlgoName.compare("SISConePlugin") == 0;
557  bool useSisConeSph = _jetAlgoName.compare("SISConeSphericalPlugin") == 0;
558  // save the 2nd parameter of the SisCones
559  double sisConeOverlapThreshold = 0;
560  if (useSisCone || useSisConeSph)
561  sisConeOverlapThreshold = atof(_jetAlgoNameAndParams[2].c_str());
562 
563  // do a maximum of N iterations
564  // For each iteration we will modify the Radius R by RDiff. RDiff will be reduced by 2 for each iteration
565  // i.e. RDiff = Pi/8, Pi/16, Pi/32, Pi/64, ...
566  // e.g.
567  // R = Pi/4
568  // R = Pi/4 + Pi/8
569  // R = Pi/4 + Pi/8 - Pi/16
570  // R = Pi/4 + Pi/8 - Pi/16 + Pi/32
571  // ...
572  for (iIter=0; iIter<ITERATIVE_INCLUSIVE_MAX_ITERATIONS; iIter++) {
573 
574  // do the clustering for this value of R. For this we need to re-initialize the JetDefinition, as it takes the R parameter
575  fastjet::JetDefinition* jetDefinition = NULL;
576 
577  // unfortunately SisCone(spherical) are being initialized differently, so we have to check for this
578  if (useSisCone) {
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);
584  } else {
585  jetDefinition = new fastjet::JetDefinition(_jetAlgoType, R, _jetRecoScheme, _strategy);
586  }
587 
588  // now we can finally create the cluster sequence
589  fastjet::ClusterSequence cs(pjList, *jetDefinition);
590 
591  jets = cs.inclusive_jets(0); // no pt cut, we will do an energy cut
592  jetsReturn.clear();
593 
594  // count the number of jets above threshold
595  nJets = 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();
600 
601  streamlog_out(DEBUG) << iIter << " " << R << " " << jets.size() << " " << nJets << std::endl;
602 
603  if (nJets == _requestedNumberOfJets) { // if the number of jets is correct: success!
604  delete pluginSisCone; pluginSisCone = NULL;
605  delete pluginSisConeSph; pluginSisConeSph = NULL;
606  delete jetDefinition;
607  break;
608 
609  } else if (nJets < _requestedNumberOfJets) {
610  // if number of jets is too small: we need a smaller Radius per jet (so
611  // that we get more jets)
612  R -= RDiff;
613 
614  } else if (nJets > _requestedNumberOfJets) { // if the number of jets is too
615  // high: increase the Radius
616  R += RDiff;
617  }
618 
619  RDiff /= 2;
620 
621  // clean up
622  delete pluginSisCone; pluginSisCone = NULL;
623  delete pluginSisConeSph; pluginSisConeSph = NULL;
624  delete jetDefinition;
625 
626  }
627 
628  if (iIter == ITERATIVE_INCLUSIVE_MAX_ITERATIONS) {
629  streamlog_out(WARNING) << "Maximum number of iterations reached. Canceling" << std::endl;
630  throw SkippedMaxIterationException( jets );
631  // Currently we will return the latest results, independent if the number is actually matched
632  // jetsReturn.clear();
633  }
634 
635  return jetsReturn;
636 }
637 
638 
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