Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

ClusterAnalysisModule - source file

// ESAF : Euso Simulation and Analysis Framework
// $Id: ClusterAnalysisModule.cc,v 1.20 2005/01/18 09:07:31 pesce Exp $
// R. Pesce created Mar,  4 2004

#include "TMath.h"
#include "RecoFramework.hh"
#include "ClusterAnalysisModule.hh"
#include "EusoCluster.hh"
#include "AutomatedClustering.hh"
#include "RecoEvent.hh"
#include "RecoPixelData.hh"
#include "RecoRootEvent.hh"

#include <TH3F.h>
#include <TH2F.h>
#include <TDirectory.h>

ClassImp(ClusterAnalysisModule)

//_____________________________________________________________________________
// ctor
 ClusterAnalysisModule::ClusterAnalysisModule() : RecoModule("ClusterAnalysis") {
}

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

//_____________________________________________________________________________
// get a specified cluster
 EusoCluster* ClusterAnalysisModule::GetCluster( Int_t i ) {
    if ( i >= 0 && i < (Int_t)fClusters.size() )
        return fClusters[i];
    else
        return NULL;
}

//_____________________________________________________________________________
// init
 Bool_t ClusterAnalysisModule::Init() {
    fEv = (RecoEvent*)0;
    Msg(EsafMsg::Info) << "Initializing " << MsgDispatch;
    fNumThreshold = (Double_t)Conf()->GetNum("ClusterAnalysisModule.NumThreshold");
    if ( fNumThreshold <= 0 ) {
        Msg(EsafMsg::Panic) << "You must specify a positive number in ClusterAnalysisModule.NumThreshold" << MsgDispatch;
    }

    fRelNumPointsMin = (Double_t)Conf()->GetNum("ClusterAnalysisModule.RelNumPointsMin");
    if ( fRelNumPointsMin <= 0 ) {
        Msg(EsafMsg::Panic) << "You must specify a positive number in ClusterAnalysisModule.RelNumPointsMin" << MsgDispatch;
    }
    
    fRelNumHits = (Double_t)Conf()->GetNum("ClusterAnalysisModule.RelNumHits");
    if ( fRelNumHits <= 0 ) {
        Msg(EsafMsg::Panic) << "You must specify a positive number in ClusterAnalysisModule.RelNumHits" << MsgDispatch;
    }

    gDoFit = (Bool_t)Conf()->GetNum("ClusterAnalysisModule.DoMedianFit");
    return true;
}

//_____________________________________________________________________________
// pre-process
 Bool_t ClusterAnalysisModule::PreProcess() {
    fMainCluster = 0;
    fThetaMean.clear();
    fPhiMean.clear();
    fClusters.clear();
    fSelectedClusters.clear();
    fId.clear();
    fFlag1.clear();
    fFlag2.clear();
    return true;
}

//_____________________________________________________________________________
// process
 Bool_t ClusterAnalysisModule::Process( RecoEvent *ev ) {
    fEv = ev;
    Int_t siglevel(0), hitsmin(0), pointsmin(0);
    const RecoModuleData *bcmd = fEv->GetModuleData("BaseClustering");
    if ( bcmd == NULL ) {
        Msg(EsafMsg::Panic) << "BaseClusteringModule is not found" << MsgDispatch;
    } else {
        if ( bcmd->GetObj("Clusters") == NULL ) return kFALSE;
        siglevel = bcmd->GetInt("SignificanceLevel");
        hitsmin = bcmd->GetInt("NumHitsMinimum");
        fThreshold = bcmd->GetDouble("DistanceThreshold");
        pointsmin = bcmd->GetInt("NumPointsMinimum");
        fClusters = *(vector<EusoCluster*>*) bcmd->GetObj("Clusters");
	Msg(EsafMsg::Debug) << "BaseClusteringModule data found." << MsgDispatch;
    }
    
    if ( fClusters.size() == 0 ) {
        Msg(EsafMsg::Panic) << "There are no significant clusters. " << fSelectedClusters.size() << MsgDispatch;
        return kFALSE;
    } else if ( fClusters.size() == 1 ) {
        Msg(EsafMsg::Debug) << "There is only one significant cluster." << MsgDispatch;
        fSelectedClusters.push_back(0);
        MyData()->Add("NumHitsMinimum", hitsmin);
       	return kTRUE;
    } else {
        Msg(EsafMsg::Debug)<< "There are " << fClusters.size() << " significant clusters." << MsgDispatch;
    }
    

    if(Conf()->GetStr("ClusterAnalysisModule.ReClustering") == "yes") {
        fClusters.clear();
        pointsmin = (Int_t) (pointsmin*fRelNumPointsMin);
        hitsmin = (Int_t) (hitsmin*fRelNumHits);
        AutomatedClustering *ac = 
            new AutomatedClustering( fEv, hitsmin, fThreshold*fNumThreshold, pointsmin, gDoFit );
        fClusters = ac->GetClusters();
        Msg(EsafMsg::Debug)<< "ReClustering: found " << fClusters.size() << " significant clusters." << MsgDispatch;
        for( Int_t i=0; i<(Int_t)fClusters.size(); i++ ) {
	    Msg(EsafMsg::Debug) << "Cluster n. " << i << ": Number of points = " << fClusters[i]->GetNumPoints() << MsgDispatch;
            if (gDoFit) {
                 Msg(EsafMsg::Debug)<< " Slope = " << fClusters[i]->GetSlope() << " Offset = "
                     << fClusters[i]->GetOffset() << " AbsoluteDeviation = "
                     << fClusters[i]->GetAbsoluteDeviation() << MsgDispatch;
            }
        }
    } else if(Conf()->GetStr("ClusterAnalysisModule.ReClustering") == "no") {
        Msg(EsafMsg::Debug) << "No ReClustering." << MsgDispatch;
    } else {
        Msg(EsafMsg::Panic) << "You must specify yes or no in ClusterAnalysisModule.RiClustering" << MsgDispatch; 
    }

    SearchMainCluster();
    
    for(Int_t i=0; i<(Int_t)fClusters.size(); i++) {
        Double_t thmin = fClusters[i]->GetThetaMin();
        Double_t thmax = fClusters[i]->GetThetaMax();
        Double_t phimin = fClusters[i]->GetPhiMin();
        Double_t phimax = fClusters[i]->GetPhiMax();
        fThetaMean.push_back( (thmin + thmax) / 2 );
        fPhiMean.push_back( (phimin + phimax) / 2 );
        Double_t sum = AngularDistance(fThetaMean[i], fPhiMean[i], thmin, phimin);
        sum += AngularDistance(fThetaMean[i], fPhiMean[i], thmin, phimax); 
        sum += AngularDistance(fThetaMean[i], fPhiMean[i], thmax, phimin); 
        sum += AngularDistance(fThetaMean[i], fPhiMean[i], thmax, phimax);
        fDim.push_back(sum / 4.);
    }

    InitSuperClustering();
    SearchSuperClusters();

    if ( fSelectedClusters.size() == 0 ) {
        fSelectedClusters.push_back(fMainCluster);
	Msg(EsafMsg::Debug) << "Selected main cluster." << MsgDispatch;
    } else {
        Msg(EsafMsg::Debug) << "Selected " << fSelectedClusters.size() << " clusters" << MsgDispatch;
        Int_t pointsum = 0;
        Bool_t mainselect = kFALSE;
        for(Int_t i=0; i<(Int_t)fSelectedClusters.size(); i++) {
            pointsum += fClusters[fSelectedClusters[i]]->GetNumPoints();
            if (fSelectedClusters[i] == fMainCluster) mainselect = kTRUE;
        }
        Msg(EsafMsg::Debug) << "Number of points clusterized = " << pointsum << MsgDispatch;
        if (mainselect==kFALSE) {
	    Msg(EsafMsg::Warning) << "Main Cluster is not selected." << MsgDispatch;
            Double_t thmean(0), phimean(0), dim(0);
            for(Int_t i=0; i<(Int_t)fSelectedClusters.size(); i++) {
                thmean += fThetaMean[fSelectedClusters[i]];
                phimean += fPhiMean[fSelectedClusters[i]];
                dim += fDim[fSelectedClusters[i]]+fThreshold;
            }
            Double_t dist = AngularDistance(thmean, phimean, fThetaMean[fMainCluster],
                            fPhiMean[fMainCluster]);
            dist -= (fDim[fMainCluster]+dim);
            if (dist <= 2.*fThreshold) {
                Msg(EsafMsg::Debug) << "Checking: Main Cluster selected" << MsgDispatch;
                pointsum += fClusters[fMainCluster]->GetNumPoints();
                Msg(EsafMsg::Debug) << "Total points clusterized = " << pointsum << MsgDispatch;
                fSelectedClusters.push_back(fMainCluster);
            } else {
                if (fClusters[fMainCluster]->GetNumPoints() > 2*pointsum) {
                    fSelectedClusters.clear();
                    fSelectedClusters.push_back(fMainCluster);
                    Msg(EsafMsg::Debug) << "Keeped Main Cluster." << MsgDispatch;
                } else if (pointsum > 2*fClusters[fMainCluster]->GetNumPoints()) {
                    Msg(EsafMsg::Debug) << "Keeped selected clusters." << MsgDispatch;
                } else {
                    fSelectedClusters.push_back(fMainCluster);
                    Msg(EsafMsg::Debug) << "Doubdt case: keeped selected clusters and main cluster" << MsgDispatch;
                }
            }
        }
    }
    MyData()->Add("NumHitsMinimum", hitsmin);

    return kTRUE;
}

//_____________________________________________________________________________
// post-process
 Bool_t ClusterAnalysisModule::PostProcess() {
    if (fSelectedClusters.size() == 0) return true;

    if ((Bool_t)Conf()->GetNum("ClusterAnalysisModule.DoHistogram")) HistoClusters();
    
    // keep only selected clusters
    vector<EusoCluster*> temp;
    EusoCluster *dummy;
    for(Int_t i=0; i<(Int_t)fSelectedClusters.size(); i++) {
        dummy = new EusoCluster();
        (*dummy) = *(fClusters[fSelectedClusters[i]]);
        temp.push_back(dummy);
    }
    fClusters.clear();
    fClusters = temp;
    temp.clear();

 
    // add data to event
    vector<EusoCluster*> *pclu = &fClusters;
    MyData()->Add("Clusters", pclu);

    fEv = (RecoEvent*)0;

    return true;
}

 Bool_t ClusterAnalysisModule::SaveRootData(RecoRootEvent *fRecoRootEvent) {
     return kTRUE;
}
// done
 Bool_t ClusterAnalysisModule::Done() {
   Msg(EsafMsg::Info) << "Completed" << MsgDispatch;
   return true;
}

// user memory clean
 void ClusterAnalysisModule::UserMemoryClean() {
    if ( fSelectedClusters.size() == 0 ) return;
    
    if (fClusters.size() != 0) {
        vector<EusoCluster*> *pclu = (vector<EusoCluster*>*)MyData()->GetObj("Clusters");
        EusoCluster *dummy;
        for(size_t i(0); i<pclu->size(); i++) {
            dummy = (*pclu)[i];
            delete dummy;
        }
        (*pclu).clear();
    MyData()->RemoveObj("Clusters");
    }
}

// get main cluster id
 void ClusterAnalysisModule::SearchMainCluster() {
    fMainCluster = 0;
    Int_t pointsmain = fClusters[0]->GetNumPoints();
    for( Int_t i=0; i<(Int_t)fClusters.size(); i++ ) {
        if ( fClusters[i]->GetNumPoints() > pointsmain ) {
            fMainCluster = i;
            pointsmain = fClusters[i]->GetNumPoints();
        }
    }
    Msg(EsafMsg::Debug) << "Cluster most populous: n." << fMainCluster << " with " << pointsmain << " points" << MsgDispatch;
}

// angular distance between two points
 Double_t ClusterAnalysisModule::AngularDistance(Double_t th1, Double_t phi1, Double_t th2, Double_t phi2) {
    return TMath::ACos( TMath::Sin(th1)*TMath::Sin(th2)*TMath::Cos(phi1)*TMath::Cos(phi2) +
           TMath::Sin(th1)*TMath::Sin(th2)*TMath::Sin(phi1)*TMath::Sin(phi2) +
           TMath::Cos(th1)*TMath::Cos(th2) );
}

// init superclustering process
 void ClusterAnalysisModule::InitSuperClustering() {
    fNumSuperClusters = 0;
    fNumNodes = 0;
    fNumPointsCluster = 0;
    fBookmark = 0;
    for(Int_t i=0; i<(Int_t)fClusters.size(); i++) {
        fId.push_back(0);
        fFlag1.push_back(kFALSE);
        fFlag2.push_back(kFALSE);
    }
}

// search for superclusters
 void ClusterAnalysisModule::SearchSuperClusters() {
    while( fNumNodes < (Int_t)fClusters.size() ) {
        NewSuperCluster();
        FillSuperCluster();
        WriteSuperCluster();
    }
}

// init a new supercluster
 void ClusterAnalysisModule::NewSuperCluster() {
    fNumPointsCluster = 1;
    fNumSuperClusters++;
    fNumNodes++;
    fFlag2[fBookmark] = kTRUE;
    fId[0] = fBookmark;
    fCounter = 0;
}

// fill a supercluster
 void ClusterAnalysisModule::FillSuperCluster() {
    while ( fCounter < fNumPointsCluster ) {
    
        Int_t k = fId[fCounter];
        if ( !fFlag1[k] ) {
            fFlag1[k] = kTRUE;
            for (Int_t i=0; i<(Int_t)fClusters.size(); i++) {
                if ( !fFlag2[i] ) {
                    Double_t dist = AngularDistance(fThetaMean[i], fPhiMean[i], fThetaMean[k], fPhiMean[k]);
                    dist -= (fDim[i]+fDim[k]);
                    if ( dist <= 2.*fThreshold ) {
                        fFlag2[i] = kTRUE;
                        fNumPointsCluster++;
                        fId[fNumPointsCluster-1] = i;
                        fNumNodes++;
                    }
                    fBookmark = i;
                }
            }
        }
        fCounter++;
    }
}

// write a supercluster
 void ClusterAnalysisModule::WriteSuperCluster() {
    // check cluster significance
    if ( fNumPointsCluster < 2 )
        return;
    
    for( Int_t i=0; i<fNumPointsCluster; i++ ) {
        fSelectedClusters.push_back(fId[i]);
    }
}

// histogram of clustering
 void ClusterAnalysisModule::HistoClusters() {

    gDirectory->cd("/");
    string pathname = Form("selectcluster%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("ClusterAnalysisModule.HistogramBinning");
    TH3F *selhisto = new TH3F("SelHisto","Selected Clustered Points",bin,phimin,phimax,bin,thmin,thmax,
                              bin,tmin,tmax);
    TH3F *reclhisto = new TH3F("RiClHisto","RiClusteredPoints",bin,phimin,phimax,bin,thmin,thmax,
                               bin,tmin,tmax);
    TH3F *histo = new TH3F("Histo","Points",bin,phimin,phimax,bin,thmin,thmax,bin,tmin,tmax);
    TH2F *selhisto2 = new TH2F("SelHisto2","Selected Clustered Points",bin,phimin,phimax,bin,thmin,thmax);
    TH2F *reclhisto2 = new TH2F("RiClHisto2","RiClusteredPoints",bin,phimin,phimax,bin,thmin,thmax);


    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");

    for(Int_t i=0; i<(Int_t)fClusters.size(); i++) {
        EusoCluster *clu = fClusters[i];
        for( Int_t j=0; j<clu->GetNumPoints(); j++) {
            RecoPixelData *pix = fEv->GetRecoPixelData(clu->GetPixelId(j));
            reclhisto->Fill( pix->GetPhi(), pix->GetTheta(), pix->GetGtu() );
            reclhisto2->Fill( pix->GetPhi(), pix->GetTheta());
        }
    }
    reclhisto->GetXaxis()->SetTitle("#varphi (rad)");
    reclhisto->GetYaxis()->SetTitle("#theta (rad)");
    reclhisto->GetZaxis()->SetTitle("GTU");
    reclhisto2->GetXaxis()->SetTitle("#varphi (rad)");
    reclhisto2->GetYaxis()->SetTitle("#theta (rad)");

    thmin = fEv->GetRecoPixelData(fClusters[fSelectedClusters[0]]->GetPixelId(0))->GetTheta();
    thmax = fEv->GetRecoPixelData(fClusters[fSelectedClusters[0]]->GetPixelId(0))->GetTheta();
    phimin = fEv->GetRecoPixelData(fClusters[fSelectedClusters[0]]->GetPixelId(0))->GetPhi();
    phimax = fEv->GetRecoPixelData(fClusters[fSelectedClusters[0]]->GetPixelId(0))->GetPhi();
    tmin = fEv->GetRecoPixelData(fClusters[fSelectedClusters[0]]->GetPixelId(0))->GetGtu();
    tmax = fEv->GetRecoPixelData(fClusters[fSelectedClusters[0]]->GetPixelId(0))->GetGtu();
    
    for(Int_t i=0; i<(Int_t)fSelectedClusters.size(); i++) {
        EusoCluster *clu = fClusters[fSelectedClusters[i]];
        for( Int_t j=0; j<clu->GetNumPoints(); j++) {
            RecoPixelData *pix = fEv->GetRecoPixelData(clu->GetPixelId(j));
            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();
            selhisto->Fill( pix->GetPhi(), pix->GetTheta(), pix->GetGtu() );
            selhisto2->Fill( pix->GetPhi(), pix->GetTheta());
        }
    }
    selhisto->GetXaxis()->SetTitle("#varphi (rad)");
    selhisto->GetYaxis()->SetTitle("#theta (rad)");
    selhisto->GetZaxis()->SetTitle("GTU");
    selhisto->GetXaxis()->SetTitle("#varphi (rad)");
    selhisto->GetYaxis()->SetTitle("#theta (rad)");


    TH3F *zoomselhisto = new TH3F("ZSelHisto","Selected Clustered points",bin,phimin,phimax,bin,thmin,thmax,
                                  bin,tmin,tmax);
    TH2F *zoomselhisto2 = new TH2F("ZSelHisto2","Selected Clustered points",bin,phimin,phimax,bin,thmin,thmax);
    for(Int_t i=0; i<(Int_t)fSelectedClusters.size(); i++) {
        EusoCluster *clu = fClusters[fSelectedClusters[i]];
        for( Int_t j=0; j<clu->GetNumPoints(); j++) {
            RecoPixelData *pix = fEv->GetRecoPixelData(clu->GetPixelId(j));
            zoomselhisto->Fill( pix->GetPhi(), pix->GetTheta(), pix->GetGtu() );
            zoomselhisto2->Fill( pix->GetPhi(), pix->GetTheta() );
        }
    }
    zoomselhisto->GetXaxis()->SetTitle("#varphi (rad)");
    zoomselhisto->GetYaxis()->SetTitle("#theta (rad)");
    zoomselhisto->GetZaxis()->SetTitle("GTU");
    zoomselhisto2->GetXaxis()->SetTitle("#varphi (rad)");
    zoomselhisto2->GetYaxis()->SetTitle("#theta (rad)");
    
    histo->Write();
    selhisto->Write();
    reclhisto->Write();
    zoomselhisto->Write();
    delete histo;
    delete selhisto;
    delete zoomselhisto;
    delete reclhisto;
    selhisto2->Write();
    reclhisto2->Write();
    zoomselhisto2->Write();
    delete selhisto2;
    delete zoomselhisto2;
    delete reclhisto2;
}
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