ClusteringComparisonCP.cc

#include <cassert>
#include "ClusteringComparisonCP.h"

namespace stab {

double TransversePairwiseIndexCP::index(const PseudoJet * inputParticles,
                    const unsigned * clustering1,
                    const unsigned * clustering2,
                    const unsigned len,
                    double * constitStab ) const
{

  assert(inputParticles);
  assert(clustering1);
  assert(clustering2);
  assert(len);
  assert(constitStab);

  long double c1 = 0.0L, c2 = 0.0L, c3 = 0.0L, c4 = 0.0L;

  for ( unsigned i=0; i<len; i++ ) {
    const double pt_i = inputParticles[i].pt();
    const unsigned cl_i_1 = clustering1[i];
    const unsigned cl_i_2 = clustering2[i];

    // weights associated to only pairs involving a particular candidate.
    long double ci1 = 0.L, ci2 = 0.L, ci3 = 0.L, ci4 = 0.L;

    for ( unsigned j=0; j<len; ++j ) {
      if ( j == i ) continue; // skip  ith particle
 
      const double weight = inputParticles[j].pt()*pt_i;
      const bool same_in_1 = (cl_i_1 == clustering1[j]);
      const bool same_in_2 = (cl_i_2 == clustering2[j]);

      if ( same_in_1 && same_in_2)
    ci1 += weight;
      else if ( same_in_1 && ! same_in_2)
    ci3 += weight;
      else if ( !same_in_1 && same_in_2)
    ci4 += weight;
      else
    ci2 += weight;
    }

    assert( ci1 + ci2 + ci3 + ci4 > 0.L );
    constitStab[i] += combineCounts(ci1,ci2,ci3,ci4);

    c1 += ci1;
    c2 += ci2;
    c3 += ci3;
    c4 += ci4;
  }

  assert(c1 + c2 + c3 + c4 > 0.L);

  // weights are divided by two because of the double counting
  // in the double cycle of i and j.  Pairs are counted twice since
  // i and j both run over the entire collection input particles.
  return combineCounts(c1/2.,c2/2.,c3/2.,c4/2.);

}  // end index

} // end stab namespace