[Rivet-svn] r2589 - in trunk: . include/Rivet/Projections src/Projections

[blackhole at projects.hepforge.org]

Author: buckley
Date: Wed Jul 14 09:51:19 2010
New Revision: 2589

Log:
Adding jet area functionality to FastJets -- a bit manky in terms of pointers/object lifetime/etc. but it should work. Feedback from users who requested it would be useful!

Modified:
   trunk/ChangeLog
   trunk/include/Rivet/Projections/FastJets.hh
   trunk/src/Projections/FastJets.cc

Modified: trunk/ChangeLog
==============================================================================
--- trunk/ChangeLog	Wed Jul 14 09:50:24 2010	(r2588)
+++ trunk/ChangeLog	Wed Jul 14 09:51:19 2010	(r2589)
@@ -1,3 +1,11 @@
+2010-07-14  Andy Buckley  <andy at insectnation.org>
+
+	* FastJets projection now allows setting of a jet area definition,
+	plus a hacky mapping for getting the area-enabled cluster
+	sequence. Requested by Pavel Starovoitov & Paolo Francavilla.
+
+	* Lots of script updates in last two weeks!
+
 2010-06-30  Andy Buckley  <andy at insectnation.org>
 
 	* Minimising amount of Log class mapped into SWIG.

Modified: trunk/include/Rivet/Projections/FastJets.hh
==============================================================================
--- trunk/include/Rivet/Projections/FastJets.hh	Wed Jul 14 09:50:24 2010	(r2588)
+++ trunk/include/Rivet/Projections/FastJets.hh	Wed Jul 14 09:51:19 2010	(r2589)
@@ -8,7 +8,10 @@
 #include "Rivet/Particle.hh"
 #include "Rivet/Jet.hh"
 
+#include "fastjet/JetDefinition.hh"
+#include "fastjet/AreaDefinition.hh"
 #include "fastjet/ClusterSequence.hh"
+#include "fastjet/ClusterSequenceArea.hh"
 #include "fastjet/PseudoJet.hh"
 
 namespace Rivet {
@@ -74,9 +77,16 @@
 
   public:
 
-    /// Reset the projection. Jet def etc are unchanged.
+    /// Reset the projection. Jet def, etc. are unchanged.
     void reset();
 
+    /// @brief Use provided jet area definition
+    /// @warning adef is NOT copied, the user must ensure that it remains valid!
+    /// Provide an adef null pointer to re-disable jet area calculation
+    void useJetArea(fastjet::AreaDefinition* adef) {
+      _adef = adef;
+    }
+
     /// Number of jets above the \f$ p_\perp \f$ cut.
     size_t numJets(double ptmin = 0.0) const;
 
@@ -120,11 +130,23 @@
       return _cseq.get();
     }
 
+    /// Return the cluster sequence (FastJet-specific).
+    const fastjet::ClusterSequenceArea* clusterSeqArea() const {
+      /// @todo Throw error if no area def? Or just blindly call dynamic_cast?
+      if (_adef == 0) return (fastjet::ClusterSequenceArea*) 0;
+      return dynamic_cast<fastjet::ClusterSequenceArea*>(_cseq.get());
+    }
+
     /// Return the jet definition (FastJet-specific).
     const fastjet::JetDefinition& jetDef() const {
       return _jdef;
     }
 
+    /// Return the area definition (FastJet-specific). May be null.
+    const fastjet::AreaDefinition* areaDef() const {
+      return _adef;
+    }
+
     /// Get the subjet splitting variables for the given jet.
     vector<double> ySubJet(const fastjet::PseudoJet& jet) const;
 
@@ -159,6 +181,9 @@
     /// Jet definition
     fastjet::JetDefinition _jdef;
 
+    /// Pointer to user-handled area definition
+    fastjet::AreaDefinition* _adef;
+
     /// Cluster sequence
     shared_ptr<fastjet::ClusterSequence> _cseq;
 

Modified: trunk/src/Projections/FastJets.cc
==============================================================================
--- trunk/src/Projections/FastJets.cc	Wed Jul 14 09:50:24 2010	(r2588)
+++ trunk/src/Projections/FastJets.cc	Wed Jul 14 09:51:19 2010	(r2589)
@@ -3,8 +3,6 @@
 #include "Rivet/Tools/Logging.hh"
 #include "Rivet/Projections/FastJets.hh"
 
-#include "fastjet/JetDefinition.hh"
-#include "fastjet/ClusterSequence.hh"
 #include "fastjet/SISConePlugin.hh"
 #include "fastjet/CDFJetCluPlugin.hh"
 #include "fastjet/CDFMidPointPlugin.hh"
@@ -17,7 +15,7 @@
 
 
   FastJets::FastJets(const FinalState& fsp, JetAlgName alg, double rparameter, double seed_threshold)
-    : JetAlg(fsp)
+    : JetAlg(fsp), _adef(0)
   {
     setName("FastJets");
     getLog() << Log::DEBUG << "R parameter = " << rparameter << endl;
@@ -59,16 +57,16 @@
 
   FastJets::FastJets(const FinalState& fsp, fastjet::JetAlgorithm type,
                      fastjet::RecombinationScheme recom, double rparameter)
-    : JetAlg(fsp)
+    : JetAlg(fsp), _adef(0)
   {
     setName("FastJets");
     _jdef = fastjet::JetDefinition(type, rparameter, recom);
   }
 
 
-  FastJets::FastJets(const FinalState& fsp, 
+  FastJets::FastJets(const FinalState& fsp,
                      fastjet::JetDefinition::Plugin& plugin)
-    : JetAlg(fsp)
+    : JetAlg(fsp), _adef(0)
   {
     setName("FastJets");
     /// @todo Should we be copying the plugin?
@@ -77,17 +75,6 @@
   }
 
 
-//   FastJets::FastJets(const FastJets& other)
-//     : JetAlg
-// //_cseq(other._cseq),
-//     _jdef(other._jdef),
-//     _plugin(other._plugin),
-//     _yscales(other._yscales)
-//   {
-//     setName("FastJets");
-//   }
-
-
   int FastJets::compare(const Projection& p) const {
     const FastJets& other = dynamic_cast<const FastJets&>(p);
     return \
@@ -95,7 +82,8 @@
       cmp(_jdef.jet_algorithm(), other._jdef.jet_algorithm()) ||
       cmp(_jdef.recombination_scheme(), other._jdef.recombination_scheme()) ||
       cmp(_jdef.plugin(), other._jdef.plugin()) ||
-      cmp(_jdef.R(), other._jdef.R());
+      cmp(_jdef.R(), other._jdef.R()) ||
+      cmp(_adef, other._adef);
   }
 
 
@@ -120,7 +108,13 @@
       ++counter;
     }
     getLog() << Log::DEBUG << "Running FastJet ClusterSequence construction" << endl;
-    _cseq.reset(new fastjet::ClusterSequence(vecs, _jdef));
+
+    // Choose CSeq as basic or area-calculating depending on whether _adef pointer is non-null.
+    if (_adef == 0) {
+      _cseq.reset(new fastjet::ClusterSequence(vecs, _jdef));
+    } else {
+      _cseq.reset(new fastjet::ClusterSequenceArea(vecs, _jdef, *_adef));
+    }
   }