10 #include "Utilities/func/MathUtil.h" 50 this->
SeedA(minClust);
53 while (
fA.size() <=maxClust){
67 if (
fA.size() > maxClust)
fA.pop_back();
76 for(
unsigned int j=0;
j<
fA.size();
j++){
79 for (
unsigned int j=0;
j<
fX.size(); ++
j) {
80 for (
unsigned int k=0; k<
fA.size(); ++k) {
100 for (i=0; i<
n; ++
i) {
101 fD[
i].resize(
fX.size());
102 fU[
i].resize(
fX.size());
107 for(
unsigned int i=0; i<
n; i++){
137 for(
unsigned int i=0;
i<
fHL.size();
i++){
139 if (
fHL[
i].fView == view){
147 if (dist == 0.0) dist = 0.0001;
148 fDist.push_back(dist);
171 for (
unsigned int i=0;
i<
n; ++
i) {
181 for (
unsigned int k=0; k<
w.size(); k++){
183 for (
unsigned int j=0;
j<
fX.size();
j++){
193 for (
unsigned int k=0; k<
w.size(); k++) {
195 if ((
w[k] >=
w.back()) && (
w[k] >
w[k+1])) {
199 else if (k==
w.size()-1) {
200 if ((
w[k] >=
w[k-1]) && (
w[k] >
w.front())) {
205 if ((
w[k] >=
w[k-1]) && (
w[k] >
w[k+1])) {
215 for (
unsigned int k=0; k<
n; k++){
240 for (i=0; i<
fA.size(); ++
i) {
241 for (j=0; j<
fX.size(); ++
j) {
245 std::cout <<
"D[" << i <<
"," << j <<
"]=(" 246 <<
fX[
j] <<
"-" <<
fA[
i] <<
")/" 260 for (j=0; j<
fX.size(); ++
j) {
261 for (i=0; i<
fA.size(); ++
i) {
263 if (
fU[i][j] !=
fU[i][j])
fU[
i][
j] = 1
e-50;
264 if (
fU[i][j] < 1
e-50)
fU[
i][
j] = 1
e-50;
266 std::cout <<
"U[" << i <<
"," << j <<
"]=" 267 <<
fU[
i][
j] <<
" | ";
286 for (i=0; i<
fA.size(); ++
i) {
291 for (j=0; j<
fX.size(); ++
j) {
292 double dtheta =
fX[
j]-musave;
293 while (dtheta<-
M_PI) dtheta += 2*
M_PI;
294 while (dtheta>
M_PI) dtheta -= 2*
M_PI;
299 double thetaNew = musave + s/
w;
300 while (thetaNew<-
M_PI) thetaNew += 2*
M_PI;
301 while (thetaNew>
M_PI) thetaNew -= 2*
M_PI;
316 for (i=0; i<
fA.size(); ++
i) {
317 for (j=0; j<
fX.size(); ++
j) {
332 bool fContinue =
false;
339 for (
unsigned int j=0;
j<
fX.size();
j++){
341 for (
unsigned int i=0;
i<
fA.size();
i++){
358 for (
unsigned int k=0; k<
w.size(); k++){
360 for (
unsigned int j=0;
j<
fUnused.size();
j++){
362 while (diff < -
M_PI) diff += 2*
M_PI;
370 for (
unsigned int k=0; k<
w.size(); k++) {
372 if ((
w[k] >=
w.back()) && (
w[k] >
w[k+1])) {
376 else if (k==
w.size()-1) {
377 if ((
w[k] >=
w[k-1]) && (
w[k] >
w.front())) {
382 if ((
w[k] >=
w[k-1]) && (
w[k] >
w[k+1])) {
396 for (
unsigned int i=0;
i<
fLMax.size();
i++){
399 for (
unsigned int j=0;
j<
fSeed.size();
j++){
433 bool recluster =
false;
434 for (
unsigned int i=0;
i<
fA.size();
i++){
435 for (
unsigned int j=
i+1;
j<
fA.size();
j++){
442 for (
unsigned int k=0; k<
fLMax.size(); k++){
444 for(
unsigned int m=0;
m<
fSeed.size();
m++){
459 fA.erase(
fA.begin()+
j);
476 unsigned int size =
fA.size();
480 for (
unsigned int j=i+1;
j<
fA.size();
j++){
486 std::vector<int> iClustCells;
487 std::vector<int> jClustCells;
489 for (
unsigned int k=0; k<
fX.size(); k++){
497 if (iClustCells.size() <= jClustCells.size()){
498 frac = (double)overlap/iClustCells.size();
503 for (
unsigned int k=0; k<
fX.size(); k++){
509 for (
unsigned int k=0; k<
fX.size(); ++k) {
515 fA.erase(
fA.begin()+
i);
516 fU.erase(
fU.begin()+
i);
517 fD.erase(
fD.begin()+
i);
524 frac = (double)overlap/jClustCells.size();
529 for (
unsigned int k=0; k<
fX.size(); k++){
535 for (
unsigned int k=0; k<
fX.size(); ++k) {
541 fA.erase(
fA.begin()+
j);
542 fU.erase(
fU.begin()+
j);
543 fD.erase(
fD.begin()+
j);
552 if (size !=
fA.size()){
double Jpfuz()
Calculate and return the score evaluating the clustering.
std::vector< std::vector< double > > fD
Dist clust i to sub. j.
double Distance(double mu, double x, double dx)
Calculate and return the distance from a hit to a cluster center.
bool ComputeA()
Update the mean cluster center positions.
fvar< T > fabs(const fvar< T > &x)
std::vector< double > fSeed
List of used seeds.
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
std::vector< double > fDist
Distance from vertex.
std::vector< DAQHit > HitList
std::vector< double > fdX
Uncertainty in angle.
A collection of associated CellHits.
bool FindUnclustered()
A function to collect the hits that have not made it to a cluster yet and assign new cluster seeds ba...
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
std::vector< double > fX
Angle to vertex.
void MakeAngles(double vXY, double vZ, geo::View_t view)
Function to turn the 2D hit coordinates into angles (radians) with respect to the vertex for clusteri...
Atom< unsigned int > MaxItr
void SeedA(unsigned int n)
A function to calculate the seeds used for the first round of clustering.
std::vector< std::pair< double, double > > fLMax
List of max densities.
void ComputeDij()
Compute all distances for j hits to i cluster centers.
std::vector< double > fA
Cluster centers.
std::vector< std::vector< double > > fU
Assign cluster i to sub. j.
void DoClustering(unsigned int minClust, unsigned int maxClust, double xy, double z, geo::View_t view)
A function to control all steps in clustering in 2D on a group of hits.
Utility class for making prongs from vertex.
fvar< T > exp(const fvar< T > &x)
rb::HitList fHL
List of hits in slice.
Fuzzy k-Means prong-finding algorithm.
double frac(double x)
Fractional part.
std::vector< double > fW
Hit weight, currently PE.
double pythag(double x, double y)
2D Euclidean distance
std::vector< std::pair< double, double > > fUnused
Unused cell pos., uncert.
bool CheckUnique()
Check if each of the cluster centers (after the algorith has converged) is unique, remove duplicates to reseed.
FuzzyKMeanAlg(const FuzzyKMeanParams ¶ms)
void ComputeUij()
Compute the membership values for j hits to i cluster centers.
void MergeClusters()
Function to merge cluster centers together that have significant overlap.
void LoadCluster(rb::Cluster slice)
Function to take in a cluster and build the hitlist needed to do the clustering.
void MakeClusters(unsigned int n)
A function to calculate membership weights for a specified number of clusters.