CRichLikeCluster.C
//-----------------------------------------------------------------------------
// $Header: /asis/offline/ceres/cool/project/RCS/CRichLikeCluster.C,v 3.0 1996/10/02 09:40:39 voigt Exp $
//
// COOL Program Library
// Copyright (C) CERES collaboration, 1996
//
// Implementation of class CRichLikeCluster.
//
//-----------------------------------------------------------------------------
#include <iostream.h>
#include <math.h>
#include <limits.h>
#include <iomanip.h>
#include "CRichLikeCluster.h"
#include "cool.h"
#include "CCollection.h"
#include "CPad.h"
#include "rw/tstack.h"
#include "rw/tpordvec.h"
#include "rw/tvordvec.h"
CRichLikeCluster::CRichLikeCluster()
{
npads = 0;
firstPad = -1;
nlocMax = 0;
minAmpLocMax = 1000;
ampSum = 0;
amp2Sum = 0;
rms2 = 0;
}
int CRichLikeCluster::split( RWTPtrOrderedVector<CRichLikeCluster>& cluster,
CCollection<CPad>& pads, int id)
{
enum padstatus { useThis=10000, alreadyUsed, doNotUse, done, locMax};
RWTPtrOrderedVector<CRichLikeCluster> newcluster; // container for subcluster
RWTPtrOrderedVector<CPad> highPads; // container for amp>=minAmp
RWTPtrOrderedVector<CPad> lowPads; // container for amp<minAmp
RWTStack<CPad*, RWTValOrderedVector<CPad*> > padStack; // stack for amp>=minAmp
if ( cluster.length()==0 || pads.length()==0 ) return 0;
//
// First get the reference amplitude
//
int ioff = firstPad;
int iend = ioff+npads;
const float startAmp = 300;
const float ampDiff = 0.001;
float maxAmp = 0;
float minAmp = startAmp;
int j;
for ( j=ioff; j<iend; j++ ) {
if ( pads(j)->getAmp() > maxAmp ) maxAmp = pads(j)->getAmp();
if ( pads(j)->getAmp() < minAmp ) minAmp = pads(j)->getAmp();
}
if (maxAmp-minAmp < ampDiff || cluster(id)->getNLocMax() == 0 ||
maxAmp-minAmpLocMax < ampDiff) return 1; // nothing to do
float amp=0, avrAmp=0, refAmp=0;
const int maxtries=3;
short clusterId;
int ix, iy, lx, ly;
for ( int tries=0; tries<=maxtries; tries++) { // main loop: try to split
lowPads.clear();
highPads.clear();
newcluster.clearAndDestroy();
padStack.clear();
//
// For each try get a new average amplitude
//
if ( tries==0 && minAmpLocMax>minAmp )
avrAmp = minAmpLocMax;
else {
refAmp = startAmp;
for (j=ioff; j<iend; j++ ) {
if (pads(j)->getIsLocalMax() && pads(j)->getAmp()>avrAmp &&
pads(j)->getAmp()<refAmp) refAmp = pads(j)->getAmp();
}
if ( refAmp-startAmp < ampDiff)
avrAmp = minAmp+0.8*(maxAmp-minAmp); // investigate these cases later
else
avrAmp = refAmp;
}
if ( avrAmp<minAmp || avrAmp>maxAmp ) {
cerr << "CRichLikeCluster::split()" << endl;
cerr << "\tERROR" << endl;
cerr << "\tWrong reference amplitude for cluster "
<< id << " after " << tries << " tries " << endl;
cerr << "\tminAmp\t" << minAmp << "\tavrAmp\t" << avrAmp
<< "\tmaxAmp\t" << maxAmp << endl;
return 1;
}
//
// Divide pads into two classes: above and below this average amplitude
//
for ( j=ioff; j<iend; j++ ) { // copy lowPads, mark others for use
if ( pads(j)->getAmp() < avrAmp ) {
lowPads.insert( pads(j) );
pads(j)->setTmp(doNotUse);
}
else
pads(j)->setTmp(useThis);
}
//
// Start splitting of high pads here
//
for ( j=ioff; j<iend; j++ ) {
if ( pads(j)->getTmp() == useThis ) {
padStack.push(pads(j)); // put pad onto the stack
pads(j)->setTmp(alreadyUsed); // avoid second usage
newcluster.append(new CRichLikeCluster); // get the subcluster
newcluster.last()->setFirstPad(highPads.length()); // get entry to new list
newcluster.last()->setMinAmpLocMax(startAmp);
while ( !padStack.isEmpty() ) {
highPads.append(padStack.pop());
highPads.last()->setIsLocalMax(True); // default
//
// the first subcluster will replace the parent in the cluster list
// all other subclusters are appended so :
clusterId = (newcluster.length()==1 ? id : cluster.length()+newcluster.length()-2);
highPads.last()->setClusterNumber(clusterId);
ix = highPads.last()->getX();
iy = highPads.last()->getY();
amp = highPads.last()->getAmp();
ly = iy;
for ( lx=ix-1; lx<ix+2; lx+=2 ) { // loop neighbours in x
if ( pads(lx,ly) == 0 ) continue; // pad must exist !!!!
if ( (pads(lx,ly)->getAmp() > amp) ||
(pads(lx,ly)->getAmp() > avrAmp &&
pads(lx,ly+1) != 0 &&
pads(lx,ly+1)->getAmp()>amp && // check diagonal
(pads(lx,ly+1)->getTmp() == alreadyUsed || // connected neighbours
pads(lx,ly+1)->getTmp() == useThis)) ||
(pads(lx,ly-1) != 0 &&
pads(lx,ly-1)->getAmp()>amp &&
(pads(lx,ly-1)->getTmp() == alreadyUsed ||
pads(lx,ly-1)->getTmp() == useThis)) ) {
highPads.last()->setIsLocalMax(False);
}
if ( pads(lx,ly)->getTmp() == useThis ) {
padStack.push(pads(lx,ly));
pads(lx,ly)->setTmp(alreadyUsed);
}
} // end loop neighbours in x
lx = ix;
for ( ly=iy-1; ly<iy+2; ly+=2 ) { // loop neighbours in y
if ( pads(lx,ly) == 0 ) continue; // pad must exist !!!!
if ( (pads(lx,ly)->getAmp() > amp) ||
(pads(lx,ly)->getAmp() > avrAmp &&
pads(lx+1,ly) != 0 &&
pads(lx+1,ly)->getAmp()>amp && // check diagonal
(pads(lx+1,ly)->getTmp() == alreadyUsed || // connected neighbours
pads(lx+1,ly)->getTmp() == useThis)) ||
(pads(lx-1,ly) != 0 &&
pads(lx-1,ly)->getAmp()>amp &&
(pads(lx-1,ly)->getTmp() == alreadyUsed ||
pads(lx-1,ly)->getTmp() == useThis)) ) {
highPads.last()->setIsLocalMax(False);
}
if ( pads(lx,ly)->getTmp() == useThis ) {
padStack.push(pads(lx,ly));
pads(lx,ly)->setTmp(alreadyUsed);
}
} // end loop neighbours in y
//
// now update all cluster properties
//
newcluster.last()->incrNPads();
newcluster.last()->updateAmp((int)highPads.last()->getAmp());
if ( highPads.last()->getIsLocalMax() ) newcluster.last()->incrNLocMax();
if ( highPads.last()->getIsLocalMax() &&
highPads.last()->getAmp()<newcluster.last()->getMinAmpLocMax() ) {
newcluster.last()->setMinAmpLocMax((int)highPads.last()->getAmp());
}
} // end while !padstack.isEmpty()
} // end if pads(j)->getTmp..
} // end split highPads
if ( newcluster.length()>1 ) break; // cluster did split, so leave split loop
} // end try to split
for ( j=ioff; j<iend; j++ )
pads(j)->setTmp(done); // avoid second usage
if ( newcluster.length()==1 ) { // cluster did not split
lowPads.clear();
highPads.clear();
newcluster.clearAndDestroy();
padStack.clear();
return 1;
}
//
// distribute the lowPads to the cluster of the closest local maximum
//
float smalldist, oldlocmax, dist;
int k,l, istart,istop, idx, idy;
for (j=0; j<lowPads.length(); j++) {
smalldist = FLT_MAX;
oldlocmax = 0;
ix = lowPads(j)->getX();
iy = lowPads(j)->getY();
for ( k=0;k<newcluster.length(); k++) { // loop over all subclusters
istart = newcluster(k)->getFirstPad();
istop = istart+newcluster(k)->getNPads();
for ( l=istart; l<istop; l++) { // loop over all pads of that subcluster
if (newcluster(k)->getNPads() > 1)
if ( !highPads(l)->getIsLocalMax()) continue;
idx = ix-highPads(l)->getX();
idy = iy-highPads(l)->getY();
dist = idx*idx+idy*idy;
// use the clusternumber at this point for the recognition of
// the subcluster, the real clusternumber will be updated when the
// pad is reinserted in the original padlist
if ( dist<smalldist ) { // this lowpad belongs to subcluster k now
smalldist = dist;
//
// here the inverse number is used to be able to set the clusterid
// for the lowPads later directly when they are reinserted
//
lowPads(j)->setClusterNumber(-k);
}
else if ( dist==smalldist ) { // in the middle, so put it to
if ( oldlocmax<highPads(l)->getAmp()) { // the higher local maximum
oldlocmax = highPads(l)->getAmp();
lowPads(j)->setClusterNumber(-k);
}
}
} // end loop over all pads of subcluster
} // end loop over all subclusters
} // end connect low pads to local maximum
//
// re-insert the lowPads into the padlist
//
int entry=firstPad, firstHigh, lastHigh;
for ( k=0; k<newcluster.length(); k++) { // copy_back: loop over all subclusters
firstHigh = newcluster(k)->getFirstPad();
lastHigh = firstHigh+newcluster(k)->getNPads();
newcluster(k)->setFirstPad(entry);
for ( l=firstHigh;l<lastHigh;l++) { // copy high pad pointer
pads(entry) = highPads(l);
entry++;
}
//
// copy low pad pointer back into the padlist, since all pads are distributed
// update the cluster id for each pad now.
//
for (l=0;l<lowPads.length();l++) {
if (lowPads(l)->getClusterNumber()==-k) {
pads(entry) = lowPads(l);
pads(entry)->setTmp(done); // avoid second usage !!!!!
//
// again: the first subcluster replaces its parent, the others are
// appended to the clusterlist
//
clusterId = (k==0 ? id : cluster.length()+k-1);
pads(entry)->setClusterNumber(clusterId);
entry++;
if (lowPads(l)->getIsLocalMax()) {
newcluster(k)->incrNLocMax();
}
newcluster(k)->incrNPads();
newcluster(k)->updateAmp(int(lowPads(l)->getAmp())); // change class cluster
}
}
newcluster(k)->rms2Calc();
} // end copy_back
//
// insert the first cluster here, avoid to resize the cluster list
// so delete the old cluster, if it did not split, newcluster(0) contains
// a deep copy
//
delete cluster(id);
cluster(id) = newcluster(0);
int nsubcluster = newcluster.length(); // needed for the return value
for (int i=1; i<nsubcluster; i++) // update cluster List
cluster.append(newcluster(i));
//
// clean up programming utilities
//
newcluster.clear();
highPads.clear();
lowPads.clear();
return nsubcluster;
} // end of splitcluster
void CRichLikeCluster::dump(const CCollection<CPad>& pads,ostream& ost)
{
const ListingNameWidth = 25;
const ListingNameWidth1 = 8;
ost.setf(ios::left);
ost << setw(ListingNameWidth) << "number of pads:" << npads << endl;
ost << setw(ListingNameWidth) << "first pad index:" << firstPad << endl;
ost << setw(ListingNameWidth) << "last pad index:" << firstPad+npads-1 << endl;
ost << setw(ListingNameWidth) << "number of local maxima:" << nlocMax << endl;
ost << setw(ListingNameWidth) << "minimal amp of loc max:" << minAmpLocMax << endl;
ost << setw(ListingNameWidth) << "sum of amplitudes:" << ampSum << endl;
ost << setw(ListingNameWidth) << "sum of amplitudes2:" << amp2Sum << endl;
ost << setw(ListingNameWidth) << "rms2:" << rms2 << endl;
for ( int j=firstPad; j<firstPad+npads; j++ ) {
ost << setw(ListingNameWidth1) << "index\t" << j
<< "\tx\t " << pads(j)->getX()
<< "\ty\t" << pads(j)->getY()
<< "\tamp\t" << pads(j)->getAmp()
<< "\tid\t" << pads(j)->getClusterNumber() << endl;
}
ost.unsetf(ios::left);
}