optimize and clean up EventShapeVariables code, add Fox-Wolfram moments

PhysicsTools/CandAlgos/plugins/EventShapeVarsProducer.cc

@@ -4,10 +4,14 @@
 
 #include "PhysicsTools/CandUtils/interface/Thrust.h"
 
+#include <vector>
+#include <memory>
+
 EventShapeVarsProducer::EventShapeVarsProducer(const edm::ParameterSet& cfg)
 {
   srcToken_ = consumes<edm::View<reco::Candidate> >(cfg.getParameter<edm::InputTag>("src"));
   r_ = cfg.exists("r") ? cfg.getParameter<double>("r") : 2.;
+  fwmax_ = cfg.exists("fwmax") ? cfg.getParameter<unsigned>("fwmax") : 0;
 
   produces<double>("thrust");
   //produces<double>("oblateness");
@@ -17,7 +21,7 @@ EventShapeVarsProducer::EventShapeVarsProducer(const edm::ParameterSet& cfg)
   produces<double>("aplanarity");
   produces<double>("C");
   produces<double>("D");
-
+  if(fwmax_ > 0) produces<std::vector<double>>("FWmoments");
 }
 
 void put(edm::Event& evt, double value, const char* instanceName)
@@ -27,7 +31,6 @@ void put(edm::Event& evt, double value, const char* instanceName)
 
 void EventShapeVarsProducer::produce(edm::Event& evt, const edm::EventSetup&)
 {
-  //std::cout << "<EventShapeVarsProducer::produce>:" << std::endl;
 
   edm::Handle<edm::View<reco::Candidate> > objects;
   evt.getByToken(srcToken_, objects);
@@ -37,12 +40,18 @@ void EventShapeVarsProducer::produce(edm::Event& evt, const edm::EventSetup&)
   //put(evt, thrustAlgo.oblateness(), "oblateness");
 
   EventShapeVariables eventShapeVarsAlgo(*objects);
+  eventShapeVarsAlgo.set_r(r_);
   put(evt, eventShapeVarsAlgo.isotropy(), "isotropy");
   put(evt, eventShapeVarsAlgo.circularity(), "circularity");
-  put(evt, eventShapeVarsAlgo.sphericity(r_), "sphericity");
-  put(evt, eventShapeVarsAlgo.aplanarity(r_), "aplanarity");
-  put(evt, eventShapeVarsAlgo.C(r_), "C");
-  put(evt, eventShapeVarsAlgo.D(r_), "D");
+  put(evt, eventShapeVarsAlgo.sphericity(), "sphericity");
+  put(evt, eventShapeVarsAlgo.aplanarity(), "aplanarity");
+  put(evt, eventShapeVarsAlgo.C(), "C");
+  put(evt, eventShapeVarsAlgo.D(), "D");
+  if(fwmax_ > 0){
+    eventShapeVarsAlgo.setFWmax(fwmax_);
+    auto vfw = std::make_unique<std::vector<double>>(eventShapeVarsAlgo.getFWmoments());
+    evt.put(std::move(vfw),"FWmoments");
+  }
 }
 
 #include "FWCore/Framework/interface/MakerMacros.h"

PhysicsTools/CandAlgos/plugins/EventShapeVarsProducer.h

@@ -17,10 +17,10 @@
  *  - aplanarity
  *  - C
  *  - D
+ *  - Fox-Wolfram moments
  *
- *  See http://cepa.fnal.gov/psm/simulation/mcgen/lund/pythia_manual/pythia6.3/pythia6301/node213.html
- *  ( http://cepa.fnal.gov/psm/simulation/mcgen/lund/pythia_manual/pythia6.3/pythia6301/node214.html )
- *  for an explanation of sphericity, aplanarity and the quantities C and D (thrust and oblateness).
+ *  See https://arxiv.org/pdf/hep-ph/0603175v2.pdf#page=524
+ *  for an explanation of sphericity, aplanarity, the quantities C and D, thrust, oblateness, Fox-Wolfram moments.
  *
  * \author Christian Veelken, UC Davis
  *
@@ -47,6 +47,7 @@ class EventShapeVarsProducer : public edm::EDProducer
 
   edm::EDGetTokenT<edm::View<reco::Candidate> > srcToken_;
   double r_;
+  unsigned fwmax_;
 
   void beginJob() override {}
   void produce(edm::Event&, const edm::EventSetup&) override;

PhysicsTools/CandAlgos/python/EventShapeVars_cff.py

@@ -8,11 +8,15 @@
 
     # momentum dependence of sphericity and aplanarity variables and of C and D quantities 
     # (r = 2. corresponds to the conventionally used default, but raises issues of infrared safety in QCD calculations;
-    #  see http://cepa.fnal.gov/psm/simulation/mcgen/lund/pythia_manual/pythia6.3/pythia6301/node213.html for more details)
-    r = cms.double(2.)
+    #  see https://arxiv.org/pdf/hep-ph/0603175v2.pdf#page=524 for more details)
+    r = cms.double(2.),
+
+    # number of Fox-Wolfram moments to compute
+    fwmax = cms.uint32(0),
 )
 
-pfEventShapeVars = copy.deepcopy(caloEventShapeVars)
-pfEventShapeVars.src = cms.InputTag("pfNoPileUp")
+pfEventShapeVars = caloEventShapeVars.clone(
+    src = cms.InputTag("pfNoPileUp")
+)
 
 produceEventShapeVars = cms.Sequence( caloEventShapeVars * pfEventShapeVars )

PhysicsTools/CandUtils/interface/EventShapeVariables.h

@@ -12,7 +12,7 @@
    cartesian, cylindrical or polar coordinates. It exploits the ROOT::TMatrixDSym 
    for the calculation of the sphericity and aplanarity.
 
-   See http://cepa.fnal.gov/psm/simulation/mcgen/lund/pythia_manual/pythia6.3/pythia6301/node213.html
+   See https://arxiv.org/pdf/hep-ph/0603175v2.pdf#page=524
    for an explanation of sphericity, aplanarity and the quantities C and D.
 
    Author: Sebastian Naumann-Emme, University of Hamburg
@@ -52,31 +52,59 @@ class EventShapeVariables {
   /// number of steps to determine how fine the granularity of the algorithm in phi should be
   double circularity(const unsigned int& numberOfSteps = 1000) const;
 
-  /// 1.5*(q1+q2) where 0<=q1<=q2<=q3 are the eigenvalues of the momemtum tensor 
+  /// set exponent for computation of momentum tensor and related products
+  void set_r(double r);
+
+  /// set number of Fox-Wolfram moments to compute
+  void setFWmax(unsigned m);
+
+  /// 1.5*(q1+q2) where q0>=q1>=q2>=0 are the eigenvalues of the momentum tensor 
   /// sum{p_j[a]*p_j[b]}/sum{p_j**2} normalized to 1. Return values are 1 for spherical, 3/4 for 
   /// plane and 0 for linear events
-  double sphericity(double = 2.)  const;
-  /// 1.5*q1 where 0<=q1<=q2<=q3 are the eigenvalues of the momemtum tensor 
+  double sphericity();
+  /// 1.5*q2 where q0>=q1>=q2>=0 are the eigenvalues of the momentum tensor 
   /// sum{p_j[a]*p_j[b]}/sum{p_j**2} normalized to 1. Return values are 0.5 for spherical and 0 
   /// for plane and linear events
-  double aplanarity(double = 2.)  const;
-  /// 3.*(q1*q2+q1*q3+q2*q3) where 0<=q1<=q2<=q3 are the eigenvalues of the momemtum tensor 
+  double aplanarity();
+  /// 3.*(q0*q1+q0*q2+q1*q2) where q0>=q1>=q2>=0 are the eigenvalues of the momentum tensor 
   /// sum{p_j[a]*p_j[b]}/sum{p_j**2} normalized to 1. Return value is between 0 and 1 
   /// and measures the 3-jet structure of the event (C vanishes for a "perfect" 2-jet event)
-  double C(double = 2.) const;
-  /// 27.*(q1*q2*q3) where 0<=q1<=q2<=q3 are the eigenvalues of the momemtum tensor 
+  double C();
+  /// 27.*(q0*q1*q2) where q0>=q1>=q2>=0 are the eigenvalues of the momentum tensor 
   /// sum{p_j[a]*p_j[b]}/sum{p_j**2} normalized to 1. Return value is between 0 and 1 
   /// and measures the 4-jet structure of the event (D vanishes for a planar event)
-  double D(double = 2.) const;
-
+  double D();
+
+  const std::vector<double>& getEigenValues() { if(!tensors_computed_) compTensorsAndVectors(); return eigenValues_; }
+  const std::vector<double>& getEigenValuesNoNorm() { if(!tensors_computed_) compTensorsAndVectors(); return eigenValuesNoNorm_; }
+  const TMatrixD& getEigenVectors() { if(!tensors_computed_) compTensorsAndVectors(); return eigenVectors_; }
+
+  double getFWmoment( unsigned l ) ;
+  const std::vector<double>& getFWmoments();
+
  private:
-  /// helper function to fill the 3 dimensional momentum tensor from the inputVectors where 
-  /// needed
-  TMatrixDSym compMomentumTensor(double = 2.) const;
-  TVectorD compEigenValues(double = 2.) const;
+  /// helper function to fill the 3 dimensional momentum tensor from the inputVectors where needed
+  /// also fill the 3 dimensional vectors of eigen-values and eigen-vectors;
+  /// the largest (smallest) eigen-value is stored at index position 0 (2)
+  void compTensorsAndVectors() ;
+
+  void computeFWmoments() ;
 
-  /// cashing of input vectors
+  /// caching of input vectors
   std::vector<math::XYZVector> inputVectors_;
+
+  /// caching of output
+  double r_;
+  bool tensors_computed_;
+  TMatrixD eigenVectors_;
+  TVectorD eigenValuesTmp_, eigenValuesNoNormTmp_;
+  std::vector<double> eigenValues_, eigenValuesNoNorm_;
+
+  /// Owen ; save computed Fox-Wolfram moments
+  unsigned fwmom_maxl_;
+  std::vector<double> fwmom_;
+  bool   fwmom_computed_ ;
+
 };
 
 #endif