JetStability.cc

#include "JetStability.h"
#include "ClusterAnalysis.h"
#include "ClusteringComparison.h"

#include "DiffusionDistFastJet.h"
#include "fastjet/ClusterSequence.hh"

namespace stab {

void JetStability::run(AbsSmearingAgent &smearAgent)
{

 // initialize the smearing agent
  smearAgent.initialize(m_parts,m_info);
 
  // cycle over  distance cutoff
  const double dDistCut = (m_maxDcut-m_minDcut)/m_nScan;
  for ( unsigned d=0; d != m_nScan; d++ ) {
    Event smearedEvent;
    smearAgent.smear(smearedEvent);

    const double dCut = m_minDcut+dDistCut*d; 
    //fastjet::ClusterSequence cs(m_parts,*m_jetDef);
    stab::DiffusionDistFastJet fjDiffDist(0.5);
    fastjet::JetDefinition jet_def(&fjDiffDist);
    fastjet::ClusterSequence cs(smearedEvent,jet_def);
    //const Event & cluster = fastjet::sorted_by_pt(cs.exclusive_jets(dCut));
    //const std::vector<fastjet::ClusterSequence::history_element> & history = cs.history();
    //const History & history = cs.history();

    //ClusterAnalysis & cAnalysis = m_clusterAnalysis[d];
    ClusterAnalysis cAnalysis(dCut);
    assert( dCut == cAnalysis.dCut() );
    cAnalysis.initialize(cs);

 
    // cycle for nCycles smearing and reclustering
    for ( unsigned c=0; c != m_nCycles; c++ ) {
      smearedEvent.clear();
      // smear the event
      smearAgent.smear(smearedEvent);
      // re-cluster the event
      fastjet::ClusterSequence cst(smearedEvent,jet_def);
      cAnalysis.fill(cst);
    }

    // fill results vectors
    const stab::TransverseRandIndex2 randIndex;
    m_pairwiseRandStability[d].push_back(cAnalysis.doStability(randIndex));
    const stab::TransverseJaccardIndex jaccardIndex; 
    m_pairwiseJaccardStability[d].push_back(cAnalysis.doStability(jaccardIndex));
    m_pairwisePtMinMatchStability[d].push_back(cAnalysis.doStability(m_ptMinMatchIndex));

    m_fuzzinessNJet[d].push_back(cAnalysis.doFuzziness(m_fuzzyNjetCalc));
    m_fuzzinessPt[d].push_back(cAnalysis.doFuzziness(m_fuzzyPtSumCalc));
    m_fuzzinessNJetQ[d].push_back(cAnalysis.doFuzziness(m_fuzzyNjetQuadCalc));
    m_fuzzinessPtQ[d].push_back(cAnalysis.doFuzziness(m_fuzzyPtSumQuadCalc));

  }
   

}

} // end stab namespace