ClusterAnalysisHighPt.cc

#include <climits>

#include "ClusterAnalysisHighPt.h"

#include "matchOneToOne.hh"

namespace stab {

void ClusterAnalysisHighPt::initialize(const fastjet::ClusterSequence &cs)
{
 
  // clear the class
  clear();

  // how many particles do we have
  m_nEvent = cs.n_particles();

  // check if there is room in our maps and vectors 
  if ( m_nEvent > m_nEventSize ) {
    m_nEventSize = m_nEvent;
    m_event.resize(m_nEvent);
    m_partMap.resize(m_nEvent);
    m_cpStab.resize(m_nEvent);
  }  

  // zero the cp stability 
  for ( unsigned i=0; i != m_nEvent; i++ )
    m_cpStab[0] = 0.;

  // fill the initial event into m_event
  const Event & allJets = cs.jets();
  for ( unsigned jp=0; jp != m_nEvent; jp++ ) {
    const int jpHistIndex = allJets[jp].cluster_sequence_history_index();
    assert(jpHistIndex >= 0);
    m_event[jp] = allJets[jp]; 
    assert( m_event[jp].pt() > 0. );
  }

 
  // perform initial clustering
  const Event & cluJets = fastjet::sorted_by_pt(cs.exclusive_jets(m_dCut));
 
  // make sure there is room in m_jets and the like 
  m_nJets=cluJets.size();
  m_nJetsInCycle = m_nJets;
  if ( m_nJetSize < m_nJets ) {
    m_jets.resize(m_nJets);
    m_jetCycleMap.resize(m_nJets);
  } 
    
  const std::vector<int> & histMap = cs.particle_jet_indices(cluJets);
  assert(m_nEvent == histMap.size());
  for ( unsigned p=0; p != m_nEvent; p++ ) {
    const int index = histMap[p];
    if ( index >= 0 && index < 2 ) m_partMap[p] = index;
    else if ( m_useJunkJet && index > 1 ) m_partMap[p] = 2U;
    else {
    m_partMap[p] = UINT_MAX;
    }
  }

}


void ClusterAnalysisHighPt::fill(const fastjet::ClusterSequence &cs)
{

  // get the clustered jets
  const Event & cluJets = fastjet::sorted_by_pt(cs.exclusive_jets(m_dCut));
  m_nJetsInCycle = cluJets.size();

  // make sure we don't exceed the number of jets seen thus far
  if ( m_nJetsInCycle > m_nJetSize ) {
    m_nJetSize = m_nJetsInCycle;
    m_jets.resize(m_nJetSize);
    m_jetCycleMap.resize(m_nJetSize);
  }
  
  // map cycle's jets to m_jets indices
  for ( unsigned j=0; j != m_nJetsInCycle; j++ ) 
    m_jetCycleMap[j]=0U;

  // match one-to-one
  std::vector<int> matchIndex;
  deltaRCalc calc;
  matchOneToOne(cluJets,m_jets,m_nJetsInCycle,m_nJets,calc,&matchIndex);

  for ( unsigned j=0; j != m_nJetsInCycle; j++ ) {
    const int index = matchIndex[j];
    if ( index >=0 ) m_jetCycleMap[j] = index;
    else { // we have encountered a new jet
      if ( m_nJetSize == m_nJets ) m_jets.resize(++m_nJetSize);
      m_jets[m_nJets] = cluJets[j];
      m_jetCycleMap[j] = m_nJets++;
    }
  }

  // fill particle maps and histograms
  std::vector<unsigned> partMap(m_nEvent);
  const std::vector<int> & histMap = cs.particle_jet_indices(cluJets);  
  assert(m_nEvent == histMap.size());
  for ( unsigned p=0; p != m_nEvent; p++ ) {
    const int index = histMap[p];
    if ( index >= 0 && index < 2 ) partMap[p] = index;
    else if ( m_useJunkJet && index > 1 ) partMap[p] = 2U;
    else partMap[p] = UINT_MAX; 
  }
  m_partMapVec.push_back(partMap); 

  // increment m_nCycle
  m_nCycle++;

  assert( m_nJetsInCycle );
}

void ClusterAnalysisHighPt::clear()
{
  m_nCycle=0U;

  // clear event
  for ( unsigned p=0; p != m_nEventSize; p++ ) {
    m_event[p] = PseudoJet();
    m_partMap[p] = 0U;
  }
  m_nEvent = 0U;
  m_partMapVec.clear();

  for ( unsigned j=0; j != m_nJetSize; j++ ) {
    m_jets[j] = PseudoJet();
    m_jetCycleMap[j] = 0U;
  } 
  m_nJets = 0U;

}


double ClusterAnalysisHighPt::doStability(const ClusteringComparisonIndexCP& index, const StabilityCalcMethod method)
{ 
  double stab=0.;
  
  switch(method) {
    case PAIRWISE:
      stab = pairwiseStability(index);
      break;
    case EXHAUSTIVE:
      stab = exhaustiveStability(index);
      break;
    case UNSMEARED:
      stab = unsmearedStability(index);
      break;
    case RESAMPLED:
      stab = resampledStability(index);
      break;
  }

  return stab;
}


double ClusterAnalysisHighPt::pairwiseStability(const ClusteringComparisonIndexCP &index)
{

  long double stab=0.L;

  for ( unsigned c=0; c != m_nCycle-1U; c++ ) {
    std::vector<unsigned> & clustMap1 = m_partMapVec[c];
    std::vector<unsigned> & clustMap2 = m_partMapVec[c+1U];
    const double temp = index.index(&m_event[0],&clustMap1[0],&clustMap2[0],m_nEvent,&m_cpStab[0]);
    stab += temp;
  } 

  // normalize the constituent perspective stability measures
  for ( unsigned i=0; i != m_nEvent; i++ )
    m_cpStab[i] /= (m_nCycle-1U);

  stab /= (m_nCycle-1);
  //std::cout << "stab=" << stab << std::endl;

  return stab; 

}

double ClusterAnalysisHighPt::exhaustiveStability(const ClusteringComparisonIndexCP &index)
{

  double stab=0.;
  
  for ( unsigned c1=0; c1 != m_nCycle-1U; c1++ ) {
    std::vector<unsigned> & clustMap1 = m_partMapVec[c1];
    for ( unsigned c2=c1+1; c2 != m_nCycle; c2++ ) {
      std::vector<unsigned> & clustMap2 = m_partMapVec[c2];
      stab += index.index(&m_event[0],&clustMap1[0],&clustMap2[0],m_nEvent,&m_cpStab[0]);
    }
  }

  // normalize the constituent perspective stability measures
  for ( unsigned i=0; i != m_nEvent; i++ )
    m_cpStab[i] /= (double)m_nCycle*((double)m_nCycle-1.)/2.;

  stab /= (double)m_nCycle*((double)m_nCycle-1.)/2.;
  return stab;
}

double ClusterAnalysisHighPt::unsmearedStability(const ClusteringComparisonIndexCP &index)
{
  double stab=0.;

  for ( unsigned c=0; c != m_nCycle; c++ ) {
    std::vector<unsigned> & clustMap = m_partMapVec[c];
    stab += index.index(&m_event[0],&m_partMap[0],&clustMap[0],m_nEvent,&m_cpStab[0]);
  }

  // normalize the constituent perspective stability measures
  for ( unsigned i=0; i != m_nEvent; i++ )
    m_cpStab[i] /= m_nCycle;

  stab /= m_nCycle;
  return stab;
}

double ClusterAnalysisHighPt::resampledStability(const ClusteringComparisonIndexCP &index)
{
  return 0.;
}

double ClusterAnalysisHighPt::doFuzziness(const AbsFuzzinessCalculator &fuzzyCalc)
{

  // construct weights
  const unsigned weightSize = m_nJets*m_nEvent;
  double weights[weightSize];
  
  for ( unsigned i=0; i != weightSize; i++ ) 
    weights[i] = 0.0;

  // fill weights as probability of initial particle to go into a final jet
  // ** this could use some optimization **
  for ( unsigned i=0; i != m_nJets; i++ ) {
    for ( unsigned j=0; j != m_nEvent; j++ ) {
      unsigned count = 0;
      for ( unsigned k=0; k != m_nCycle; k++ ) 
    if ( m_partMapVec[k][j] == i ) count++;
      weights[i*m_nEvent+j] = (double)count/(double)m_nCycle;
      assert( weights[i*m_nEvent+j] >= 0. ); 
    }
  }

  bool lastFinalJetIsUnclustered=false;
  const unsigned lenJetFuzziness = m_nJets;
  double jetFuzziness[lenJetFuzziness];

  return fuzzyCalc.calculate(m_event, m_jets, weights, m_nJets, m_nEvent, lastFinalJetIsUnclustered,
    jetFuzziness, lenJetFuzziness);
   
  return 0.;
}

} // end stab namespace