// ESAF : Euso Simulation and Analysis Framework
// $Id: AutomatedClustering.cc,v 1.5 2004/11/25 18:49:30 naumov Exp $
// R. Pesce created Mar, 4 2004
#include "TMath.h"
#include "AutomatedClustering.hh"
#include "MedianFit.hh"
#include "RecoEvent.hh"
#include "RecoPixelData.hh"
#include "EusoCluster.hh"
ClassImp(AutomatedClustering)
// ctor with given values
AutomatedClustering::AutomatedClustering( RecoEvent *ev, Int_t numhitsmin, Double_t thres, Int_t numpmin,
Bool_t fit ) {
fEv = ev;
fNumHitsMinimum = numhitsmin;
fThreshold = thres;
fNumPointsMinimum = numpmin;
gDoFit = fit;
Clustering();
}
// ctor with natural values
AutomatedClustering::AutomatedClustering( RecoEvent *ev, Int_t numhitsmin, Int_t siglev, Bool_t fit ) {
fEv = ev;
fNumHitsMinimum = numhitsmin;
fSignificanceLevel = siglev;
fThreshold = -1.;
gDoFit = fit;
Clustering();
}
// dtor
AutomatedClustering::~AutomatedClustering() {
}
// process
void AutomatedClustering::Clustering() {
fTotalNumPoints = fEv->GetHeader().GetNumActiveFee();
fNumPoints = 0;
fNumClusters = 0;
fNumSigClusters = 0;
fNumNodes = 0;
fNumPointsCluster = 0;
fBookmark = 0;
fHitted.clear();
fId.clear();
fFlag1.clear();
fFlag2.clear();
for(Int_t i=0; i<fTotalNumPoints; i++) {
fHitted.push_back(0);
fId.push_back(0);
fFlag1.push_back(kFALSE);
fFlag2.push_back(kFALSE);
}
// get id of most populous macrocell (main macrocell)
fMainMacroCellId = fEv->GetMainMacroCell(1);
// search pixels with num hits > minimum
for( Int_t i=0; i<fTotalNumPoints; i++ ) {
if ( fEv->GetRecoPixelData(i)->GetCounts() >= fNumHitsMinimum ) {
fHitted[fNumPoints] = i;
fNumPoints++;
}
}
if ( fThreshold == -1. ) {
CalculateDensity();
CalculateNaturalThreshold();
CalculateNumClustersUniform();
CalculateNumPointsMinimum();
}
SearchClusters();
}
// write cluster
void AutomatedClustering::WriteCluster() {
// check cluster significance
if ( fNumPointsCluster < fNumPointsMinimum )
return;
// write the cluster
fNumSigClusters++;
EusoCluster *cl = new EusoCluster();
Double_t thmin, thmax, phimin, phimax;
Int_t gtumin, gtumax;
thmin = fEv->GetRecoPixelData(fHitted[fId[0]])->GetTheta();
thmax = fEv->GetRecoPixelData(fHitted[fId[0]])->GetTheta();
phimin = fEv->GetRecoPixelData(fHitted[fId[0]])->GetPhi();
phimax = fEv->GetRecoPixelData(fHitted[fId[0]])->GetPhi();
gtumin = fEv->GetRecoPixelData(fHitted[fId[0]])->GetGtu();
gtumax = fEv->GetRecoPixelData(fHitted[fId[0]])->GetGtu();
Double_t theta(0),phi(0);
Int_t gtu(0);
for( Int_t i=0; i<fNumPointsCluster; i++ ) {
theta = fEv->GetRecoPixelData(fHitted[fId[i]])->GetTheta();
phi = fEv->GetRecoPixelData(fHitted[fId[i]])->GetPhi();
gtu = fEv->GetRecoPixelData(fHitted[fId[i]])->GetGtu();
if ( theta < thmin ) thmin = theta;
if ( theta > thmax ) thmax = theta;
if ( phi < phimin ) phimin = phi;
if ( phi > phimax ) phimax = phi;
if ( gtu < gtumin ) gtumin = gtu;
if ( gtu > gtumax ) gtumax = gtu;
cl->AddPoint( fHitted[fId[i]] );
}
cl->SetThetaMin( thmin );
cl->SetThetaMax( thmax );
cl->SetPhiMin( phimin );
cl->SetPhiMax( phimax );
cl->SetGtuMin( gtumin );
cl->SetGtuMax( gtumax );
if ( gDoFit ) FitCluster( cl );
fClusters.push_back(cl);
}
// fit
void AutomatedClustering::FitCluster( EusoCluster *clu ) {
vector<Double_t> thetavector;
vector<Double_t> phivector;
for( Int_t i=0; i<clu->GetNumPoints(); i++ ) {
RecoPixelData *pix = fEv->GetRecoPixelData( clu->GetPixelId(i) );
thetavector.push_back(pix->GetTheta());
phivector.push_back(pix->GetPhi());
}
MedianFit *mf = new MedianFit( clu->GetNumPoints(), phivector, thetavector );
clu->SetFitted( kTRUE );
clu->SetSlope( mf->GetSlope() );
clu->SetOffset( mf->GetOffset() );
clu->SetAbsoluteDeviation( mf->GetAbsoluteDeviation() );
}
// calculate density of points by the main macrocell
void AutomatedClustering::CalculateDensity() {
Double_t nph = (Double_t)fEv->GetNumPx(fMainMacroCellId, fNumHitsMinimum);
Double_t area = (Double_t)fEv->GetThetaRange(fMainMacroCellId, fNumHitsMinimum);
area *= (Double_t)fEv->GetPhiRange(fMainMacroCellId, fNumHitsMinimum);
fDensity = nph / area;
}