Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

AutomatedClustering - source file

// 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;
}

About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update Wed Nov 16 16:57:39 2005 Wed Nov 16 16:29:22 2005