Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

GTUClusteringModule - source file

// ESAF : Euso Simulation and Analysis Framework
// $Id: GTUClusteringModule.cc,v 1.21 2005/01/18 10:06:16 pesce Exp $
// R. Pesce created Mar, 11 2004

#include "GTUClusteringModule.hh"
#include "RecoEvent.hh"
#include "RecoPixelData.hh"
#include "EusoCluster.hh"
#include "RecoRootEvent.hh"


#include <TDirectory.h>
#include <map>
#include <vector>
#include <TMath.h>
#include <TH3F.h>

ClassImp(GTUClusteringModule)
//_____________________________________________________________________________
// ctor
 GTUClusteringModule::GTUClusteringModule() : RecoModule("GTUClustering"){
}

//_____________________________________________________________________________
// dtor
 GTUClusteringModule::~GTUClusteringModule() {
}

//_____________________________________________________________________________
// init
 Bool_t GTUClusteringModule::Init() {
    fEv = (RecoEvent*)0;
    Msg(EsafMsg::Info) << "Initializing " << MsgDispatch;
    Msg(EsafMsg::Debug) << "with GTUClusteringModule.NaturalValues = " << Conf()->GetNum("GTUClusteringModule.NaturalValues")<< MsgDispatch;

    fNumHitsMinimum = (Int_t)Conf()->GetNum("GTUClusteringModule.NumHitsMinimum");
    if ( fNumHitsMinimum <= 0 ) {
        Msg(EsafMsg::Panic) << "You must specify a positive number in GTUClusteringModule.NumHitsMinimum" << MsgDispatch;
    }
    if ((Int_t)Conf()->GetNum("GTUClusteringModule.NaturalValues")==1) {
        fSignificanceLevel = (Int_t)Conf()->GetNum("GTUClusteringModule.SignificanceLevel");
        if ( fSignificanceLevel <= 0 ) {
            Msg(EsafMsg::Panic) << "You must specify a positive number in GTUClusteringModule.SignificanceLevel" << MsgDispatch;
        }
        fNumPointsMinimum = 0;
        fUseNatural = kTRUE;
    } else if ((Int_t)Conf()->GetNum("GTUClusteringModule.NaturalValues")==0) {
        fNumPointsMinimum = (Int_t)Conf()->GetNum("GTUClusteringModule.NumPointsMinimum");
        if ( fNumPointsMinimum <= 0 ) {
            Msg(EsafMsg::Panic) << "You must specify a positive number in GTUClusteringModule.NumPointsMinimum" << MsgDispatch;
         }
        fSignificanceLevel = 0;
        fUseNatural = kFALSE;
        fThreshold = (Double_t)Conf()->GetNum("GTUClusteringModule.DistanceThreshold");
        if ( fThreshold <= 0 ) {
            Msg(EsafMsg::Panic) << "You must specify a positive number in GTUClusteringModule.DistanceThreshold" << MsgDispatch;
        }
    } else {
        Msg(EsafMsg::Panic) << "You must specify 0 or 1 in GTUClusteringModule.NaturalValues" << MsgDispatch;
    }
    
    return kTRUE;
}
//_____________________________________________________________________________
// pre-process
 Bool_t GTUClusteringModule::PreProcess() {
    fSpClusters.clear();
    fPixels.clear();
    fCluPixels.clear();
    fId.clear();
    fFlag1.clear();
    fFlag2.clear();
    fNumPoints = 0;
    if ( fUseNatural==kTRUE ) {
        fNumPointsMinimum = 0;
        fThreshold = 0;
        fDensity = 0;
    }
        
    return kTRUE;
}

//_____________________________________________________________________________
// process
 Bool_t GTUClusteringModule::Process( RecoEvent* ev ) {
    fEv = ev;
    const RecoModuleData *camd = fEv->GetModuleData("ClusterAnalysis");
    if ( camd == NULL ) {
        Msg(EsafMsg::Panic) << "ClusterAnalysisModule is not found." << MsgDispatch;
    } else {
        if ( camd->GetObj("Clusters") == NULL ) {
	    Msg(EsafMsg::Panic) << "No clusters to process." << MsgDispatch;
            return kFALSE;
        }
        fSpClusters = *(vector<EusoCluster*>*) camd->GetObj("Clusters");
        Msg(EsafMsg::Debug) << "ClusterAnalysisModule data found." << MsgDispatch;
    }
    Int_t hitsmin = camd->GetInt("NumHitsMinimum");
    // search pixels in angular range with num hits minimum
    if (hitsmin==1) {
        for( Int_t i=0; i<(Int_t)fSpClusters.size(); i++ ) {
            for( Int_t j=0; j<fSpClusters[i]->GetNumPoints(); j++ ) {
                fPixels.push_back(fSpClusters[i]->GetPixelId(j));
            }
        }
    } else {
        for( Int_t i=0; i<(Int_t)fSpClusters.size(); i++ ) {
            Double_t thmin = fSpClusters[i]->GetThetaMin();
            Double_t thmax = fSpClusters[i]->GetThetaMax();
            Double_t phimin = fSpClusters[i]->GetPhiMin();
            Double_t phimax = fSpClusters[i]->GetPhiMax();
            Double_t tmin = fSpClusters[i]->GetGtuMin();
            Double_t tmax = fSpClusters[i]->GetGtuMax();
            for( Int_t j=0; j<fEv->GetHeader().GetNumActiveFee(); j++ ) {
                RecoPixelData *pix = fEv->GetRecoPixelData(j);
                if ( pix->GetTheta() >= thmin && pix->GetTheta() <= thmax &&
                    pix->GetPhi() >= phimin && pix->GetPhi() <= phimax &&
                    pix->GetGtu() >= tmin && pix->GetGtu() <= tmax &&
                    pix->GetCounts() >= fNumHitsMinimum )
                    fPixels.push_back(j);
            }

        }
    }
    fNumPoints = fPixels.size();

    Msg(EsafMsg::Debug) << "Pixels found = " << fNumPoints << MsgDispatch;

    if ( fUseNatural == kTRUE ) {
        CalculateDensity();
        CalculateThreshold();
        CalculateNumClustersUniform();
        CalculateNumPointsMinimum();
        //cout << "##### " << fDensity << " " << fThreshold <<" " << fNumClustersUniform 
             //<< " " << fNumPointsMinimum << endl;
    }
    Prepare();
    SearchClusters();

    Msg(EsafMsg::Debug) << "found " << fTmClusters.size() << " clusters." << MsgDispatch;
    if ( fTmClusters.size() == 0 ) return kFALSE;

    vector<Int_t> medTime;
    vector<Int_t> numPC;
    for( size_t i(0); i<fTmClusters.size(); i++ ) {
        Int_t sum;
        for( size_t j(0); j<(fTmClusters[i]).size(); j++ ) {
            sum += fEv->GetRecoPixelData((fTmClusters[i])[j])->GetGtu();
        }
        medTime.push_back(sum/(fTmClusters[i]).size());
        numPC.push_back((fTmClusters[i]).size());
    }
    Int_t cluMostPop(0);
    for( size_t i(0); i<numPC.size(); i++ ) {
        if ( numPC[i] > numPC[cluMostPop] ) cluMostPop = i;
    }
    Int_t timeMP = medTime[cluMostPop];

    vector<Int_t> selTClu;
    selTClu.push_back(cluMostPop);

    for( size_t i(0); i<medTime.size(); i++ ) {
        if( TMath::Abs(timeMP-medTime[i]) < (size_t)2*fThreshold && i!=(size_t)cluMostPop ) {
            selTClu.push_back(i);
        }
    }
    Msg(EsafMsg::Debug) << "Selected " << selTClu.size() << " clusters." << MsgDispatch;        
    vector<Int_t> tclu;
    for( size_t i=0; i<selTClu.size(); i++ ) {
        tclu = fTmClusters[selTClu[i]];
        for(Int_t j=0; j<(Int_t)tclu.size(); j++) {
            fCluPixels.push_back(tclu[j]);
        }
    }

    return kTRUE;
}

//_____________________________________________________________________________
// post-process
 Bool_t GTUClusteringModule::PostProcess() {
    if ( fPixels.size() == 0 ) return kTRUE;
    if ((Int_t)Conf()->GetNum("GTUClusteringModule.DoHistogram") == 1) 
        Histo();

    vector< vector<Int_t> > *pclu = &fTmClusters;
    vector<Int_t> *pix = &fCluPixels;
    MyData()->Add("TmClusters", pclu);
    MyData()->Add("CluPixels", pix);
    
    fEv = (RecoEvent*)0;
    
    return kTRUE;
}

//_____________________________________________________________________________
 Bool_t GTUClusteringModule::SaveRootData(RecoRootEvent *fRecoRootEvent) {
     return kTRUE;
}

//_____________________________________________________________________________
// done
 Bool_t GTUClusteringModule::Done() {
    Msg(EsafMsg::Info) << "Completed" << MsgDispatch;
    return kTRUE;
}

//_____________________________________________________________________________
// memory clean
 void GTUClusteringModule::UserMemoryClean() {
    if ( fPixels.size() == 0 ) return;

    if ( fTmClusters.size() != 0 ) {
        vector< vector<Int_t> > *pclu = (vector< vector<Int_t> >*)MyData()->GetObj("TmClusters");
        if ( pclu != 0 && (*pclu).size() != 0 ) {
            (*pclu).clear();
            MyData()->RemoveObj("TmClusters");
        }
    }    
    
    if ( fCluPixels.size() !=0 ) {
        vector<Int_t> *pix = (vector<Int_t>*)MyData()->GetObj("CluPixels");
        if ( pix != NULL && (*pix).size() != 0 ) {
            (*pix).clear();
            MyData()->RemoveObj("CluPixels");
        }
    }
}

// density: is simply the num of pixel minimum in a gtu
 void GTUClusteringModule::CalculateDensity() {
    Int_t tmin = 0;
    Int_t tmax = 0;
    for( Int_t i=0; i<(Int_t)fPixels.size(); i++ ) {
        Int_t gtuid = fEv->GetRecoPixelData(fPixels[i])->GetGtu();
        if ( gtuid < tmin ) tmin = gtuid;
        if ( gtuid > tmax ) tmax = gtuid;
    }
    fDensity = fPixels.size()/(tmax-tmin+1);
}

// gtu distance
 Int_t GTUClusteringModule::GTUDistance(Int_t n1, Int_t n2) {
    return TMath::Abs( fEv->GetRecoPixelData(fPixels[n1])->GetGtu() - 
                       fEv->GetRecoPixelData(fPixels[n2])->GetGtu() );
}

// natural threshold
 void GTUClusteringModule::CalculateThreshold() {
    fThreshold = TMath::Log(2)/fDensity;
}

// calculate number of clusters expected in a distribution uniform
 void GTUClusteringModule::CalculateNumClustersUniform() {
    fNumClustersUniform = (Int_t)( 1 + ( fNumPoints - 1 ) *
        TMath::Exp( - fDensity * fThreshold ) );
}

// calculate the minimum number of points in a significant cluster
 void GTUClusteringModule::CalculateNumPointsMinimum() {
    fNumPointsMinimum = (Int_t) ( fNumPoints / fNumClustersUniform 
        + fSignificanceLevel * TMath::Sqrt( (Double_t)( fNumPoints / fNumClustersUniform ) ) );
}

// preparation for clustering
 void GTUClusteringModule::Prepare() {
    for( Int_t i=0; i<fNumPoints; i++ ) {
        fId.push_back(0);
        fFlag1.push_back(kFALSE);
        fFlag2.push_back(kFALSE);
    }
    fNumClusters = 0;
    fNumNodes = 0;
    fNumSigClusters = 0;
    fBookmark = 0;
    fTmClusters.clear();
}

// search for clusters
 void GTUClusteringModule::SearchClusters() {
    while( fNumNodes < fNumPoints ) {
        NewCluster();
        FillCluster();
        WriteCluster();
    }
}

// init a new cluster
 void GTUClusteringModule::NewCluster() {
    fNumPointsCluster = 1;
    fNumClusters++;
    fNumNodes++;
    fFlag2[fBookmark] = kTRUE;
    fId[0] = fBookmark;
    fCounter = 0;
}

// fill a cluster
 void GTUClusteringModule::FillCluster() {
    while ( fCounter < fNumPointsCluster ) {
    
        Int_t k = fId[fCounter];
        if ( !fFlag1[k] ) {
            fFlag1[k] = kTRUE;
            for (Int_t i=0; i<fNumPoints; i++) {
                if ( !fFlag2[i] ) {
                    if ( GTUDistance(i, k) <= fThreshold ) {
                        fFlag2[i] = kTRUE;
                        fNumPointsCluster++;
                        fId[fNumPointsCluster-1] = i;
                        fNumNodes++;
                    }
                    fBookmark = i;
                }
            }
        }
        fCounter++;
    }
}

// write a cluster
 void GTUClusteringModule::WriteCluster() {
    // check cluster significance
    if ( fNumPointsCluster < fNumPointsMinimum )
        return;
    
    vector<Int_t> clu;
    for( Int_t i=0; i<fNumPointsCluster; i++ ) {
        clu.push_back(fPixels[fId[i]]);
    }
    fTmClusters.push_back(clu);
}

// histogram
 void GTUClusteringModule::Histo() {

    gDirectory->cd("/");
    string pathname = Form("timecluster%d", fEv->GetHeader().GetNum());
    TDirectory *path = new TDirectory(pathname.c_str(),pathname.c_str());
    path->cd();

    Float_t thmin = fEv->GetRecoPixelData(0)->GetTheta();
    Float_t thmax = fEv->GetRecoPixelData(0)->GetTheta();
    Float_t phimin = fEv->GetRecoPixelData(0)->GetPhi();
    Float_t phimax = fEv->GetRecoPixelData(0)->GetPhi();
    Float_t tmin = fEv->GetRecoPixelData(0)->GetGtu();
    Float_t tmax = fEv->GetRecoPixelData(0)->GetGtu();

    for( Int_t i=1; i<(Int_t)fEv->GetHeader().GetNumActiveFee(); i++ ) {
        if ( fEv->GetRecoPixelData(i)->GetTheta() < thmin )
            thmin = fEv->GetRecoPixelData(i)->GetTheta();
        if ( fEv->GetRecoPixelData(i)->GetTheta() > thmax )
            thmax = fEv->GetRecoPixelData(i)->GetTheta();
        if ( fEv->GetRecoPixelData(i)->GetPhi() < phimin )
            phimin = fEv->GetRecoPixelData(i)->GetPhi();
        if ( fEv->GetRecoPixelData(i)->GetPhi() > phimax )
            phimax = fEv->GetRecoPixelData(i)->GetPhi();
        if ( fEv->GetRecoPixelData(i)->GetGtu() < tmin )
            tmin = fEv->GetRecoPixelData(i)->GetGtu();
        if ( fEv->GetRecoPixelData(i)->GetGtu() > tmax )
            tmax = fEv->GetRecoPixelData(i)->GetGtu();
}

    Int_t bin = (Int_t)Conf()->GetNum("GTUClusteringModule.HistogramBinning");;
    TH3F *histo = new TH3F("Histo","Points",bin,phimin,phimax,bin,thmin,thmax,bin,tmin,tmax);
    
    for( Int_t i=0; i<(Int_t)fEv->GetHeader().GetNumActiveFee(); i++ ) {
        histo->Fill( fEv->GetRecoPixelData(i)->GetPhi(),fEv->GetRecoPixelData(i)->GetTheta(),
                     fEv->GetRecoPixelData(i)->GetGtu() );
    }
    histo->GetXaxis()->SetTitle("#varphi (rad)");
    histo->GetYaxis()->SetTitle("#theta (rad)");
    histo->GetZaxis()->SetTitle("GTU");

    TH3F *tmhisto = new TH3F("TmHisto","Clustered Points",bin,phimin,phimax,bin,thmin,thmax,
                              bin,tmin,tmax);
    thmin = fEv->GetRecoPixelData(fCluPixels[0])->GetTheta();
    thmax = fEv->GetRecoPixelData(fCluPixels[0])->GetTheta();
    phimin = fEv->GetRecoPixelData(fCluPixels[0])->GetPhi();
    phimax = fEv->GetRecoPixelData(fCluPixels[0])->GetPhi();
    tmin = fEv->GetRecoPixelData(fCluPixels[0])->GetGtu();
    tmax = fEv->GetRecoPixelData(fCluPixels[0])->GetGtu();
    
    for(Int_t i=0; i<(Int_t)fCluPixels.size(); i++) {
        RecoPixelData *pix = fEv->GetRecoPixelData(fCluPixels[i]);
        if (pix->GetTheta() < thmin ) thmin = pix->GetTheta();
        if (pix->GetTheta() > thmax ) thmax = pix->GetTheta();
        if (pix->GetPhi() < phimin ) phimin = pix->GetPhi();
        if (pix->GetPhi() > phimax ) phimax = pix->GetPhi();
        if (pix->GetGtu() < tmin ) tmin = pix->GetGtu();
        if (pix->GetGtu() > tmax ) tmax = pix->GetGtu();
        tmhisto->Fill( pix->GetPhi(), pix->GetTheta(), pix->GetGtu() );
    }
    tmhisto->GetXaxis()->SetTitle("#varphi (rad)");
    tmhisto->GetYaxis()->SetTitle("#theta (rad)");
    tmhisto->GetZaxis()->SetTitle("GTU");

    TH3F *zoomtmhisto = new TH3F("ZTmHisto","Clustered points",bin,phimin,phimax,bin,thmin,thmax,
                                  bin,tmin,tmax);
    for(Int_t i=0; i<(Int_t)fCluPixels.size(); i++) {
        RecoPixelData *pix = fEv->GetRecoPixelData(fCluPixels[i]);
        zoomtmhisto->Fill( pix->GetPhi(), pix->GetTheta(), pix->GetGtu() );
    }
    zoomtmhisto->GetXaxis()->SetTitle("#varphi (rad)");
    zoomtmhisto->GetYaxis()->SetTitle("#theta (rad)");
    zoomtmhisto->GetZaxis()->SetTitle("GTU");
    
    histo->Write();
    tmhisto->Write();
    zoomtmhisto->Write();
    delete histo;
    delete tmhisto;
    delete zoomtmhisto;
}
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