JetStability.h

#ifndef FastJetStability_JetStability_h
#define FastJetStability_JetStability_h

#include "JetStabilityDefs.h"
#include "StabilityUserInfo.h"
#include "AbsSmearingAgent.h"
#include "ClusterAnalysis.h"
#include "PtPartitioningIndex.h"
#include "TransverseFuzzinessCalculator.h"

#include "fastjet/ClusterSequence.hh"
#include "fastjet/JetDefinition.hh"

namespace stab { 

class JetStability {

public:

  inline JetStability ( const fastjet::JetDefinition *jetDef)
    :m_jetDef(jetDef)
    ,m_nCycles(100U)
    ,m_fuzzyNjetCalc(1,0,0)
    ,m_fuzzyPtSumCalc(0,1,0)
    ,m_fuzzyNjetQuadCalc(1,0,1)
    ,m_fuzzyPtSumQuadCalc(0,1,1)
   { }

  inline void reset( const Event &parts, const StabilityEventInfo &info) {
    m_parts = parts;
    m_partsSize = parts.size();
    m_info = info;
    m_infoSize = info.size();
  
    assert( m_partsSize == m_infoSize );
  
    for ( unsigned d=0; d != m_nScan; d++ ) {
      m_clusterAnalysis.clear();
    }
    
  }

  void run(AbsSmearingAgent &);  

  inline void setScan(const unsigned nScan, const double min, const double max)
  {
    m_nScan = nScan;
    m_minDcut = min;
    m_maxDcut = max;

    m_clusterAnalysis.resize(m_nScan);
    m_pairwiseRandStability.resize(m_nScan);
    m_pairwiseJaccardStability.resize(m_nScan);
    m_pairwisePtMinMatchStability.resize(m_nScan);

    m_fuzzinessNJet.resize(m_nScan);
    m_fuzzinessPt.resize(m_nScan);
    m_fuzzinessNJetQ.resize(m_nScan);
    m_fuzzinessPtQ.resize(m_nScan);

    const double dDistCut = (m_maxDcut-m_minDcut)/m_nScan;
    for ( unsigned d=0; d != m_nScan; d++ ) {
      const double dCut = m_minDcut+dDistCut*d;
      m_clusterAnalysis[d] = ClusterAnalysis(dCut);
    }
  }

  inline void setNCycles(const unsigned nCycles) { 
    m_nCycles = nCycles; 
  }

  //--------------- accessor methods
  inline const std::vector<std::vector<double> > & pairwiseRandStability() const
    { return m_pairwiseRandStability; }

  inline const std::vector<std::vector<double> > & fuzzinessNJet() const
    { return m_fuzzinessNJet; }

  inline const double scanMin() const { return m_minDcut; }
  inline const double scanMax() const { return m_maxDcut; }
  inline const unsigned nScan() const { return m_nScan; }
  inline const unsigned nCycles() const { return m_nCycles; }


private:

  // --------------------  private member data ------------------------
  // initial un-clustered event (no smearing)
  Event m_parts;
  unsigned m_partsSize;

  StabilityEventInfo m_info;
  unsigned m_infoSize;

  double m_minDcut;
  double m_maxDcut;
  unsigned m_nScan;

  // FastJet Jet definition
  const fastjet::JetDefinition *m_jetDef;

  // number of smearing cycles
  unsigned m_nCycles;

  // ClusterAnalysis handles the stability calculations
  // it has a length of m_nScan
  std::vector<ClusterAnalysis> m_clusterAnalysis;

  // stability and fuzziness calculators
  stab::BestPtMatchIndex m_ptMinMatchIndex;

  stab::TransverseFuzzinessCalculator m_fuzzyNjetCalc;
  stab::TransverseFuzzinessCalculator m_fuzzyPtSumCalc;
  stab::TransverseFuzzinessCalculator m_fuzzyNjetQuadCalc;
  stab::TransverseFuzzinessCalculator m_fuzzyPtSumQuadCalc;


  // accumulate the results of the cluser analysis 
  // each of these vectors also has legnth=m_nScan
  std::vector<std::vector<double> > m_pairwiseRandStability;
  std::vector<std::vector<double> > m_pairwiseJaccardStability;
  std::vector<std::vector<double> > m_pairwisePtMinMatchStability;

  std::vector<std::vector<double> > m_fuzzinessNJet;
  std::vector<std::vector<double> > m_fuzzinessPt;
  std::vector<std::vector<double> > m_fuzzinessNJetQ;  // added in quadrature
  std::vector<std::vector<double> > m_fuzzinessPtQ;    // added in quadrature

};


} // end stab namespace 


#endif