Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

BaseClusteringModule - source file

// ESAF : Euso Simulation and Analysis Framework
// $Id: BaseClusteringModule.cc,v 1.34 2005/06/06 14:07:27 pesce Exp $
// R. Pesce created Jan, 19 2004

#include "BaseClusteringModule.hh"
#include "RecoEvent.hh"
#include "RecoRootEvent.hh"
#include "RecoCellInfo.hh"
#include "TDirectory.h"
#include "TMath.h"
#include "Config.hh"
#include "RecoPixelData.hh"
#include "EusoCluster.hh"
#include "MedianFit.hh"

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

ClassImp(BaseClusteringModule)

// ctor
 BaseClusteringModule::BaseClusteringModule() : 
    RecoModule(string("BaseClustering")) {
}

// dtor
 BaseClusteringModule::~BaseClusteringModule() {
}

// init
 Bool_t BaseClusteringModule::Init() {
    fEv = (RecoEvent*)0;
    Msg(EsafMsg::Info) << "Initializing " << MsgDispatch;
    fAlpha = 0.55;
    fNumHitsMinimum = (Int_t)Conf()->GetNum("BaseClusteringModule.NumHitsMinimum");
    if ( fNumHitsMinimum <= 0 )
        Msg(EsafMsg::Panic) << "You must specify a positive number in BaseClusteringModule.NumHitsMinimum" << MsgDispatch;  
    fSignificanceLevel = (Int_t)Conf()->GetNum("BaseClusteringModule.SignificanceLevel");
    if ( fSignificanceLevel <= 0 )
        Msg(EsafMsg::Panic) << "You must specify a positive number in BaseClusteringModule.SignificanceLevel" << MsgDispatch;  
    
    return kTRUE;
}

// pre-process
 Bool_t BaseClusteringModule::PreProcess() {
    fTotalNumPoints = 0;
    fNumPoints = 0;
    fNumClusters = 0;
    fNumSigClusters = 0;
    fNumNodes = 0;
    fNumPointsCluster = 0;
    fBookmark = 0;
    fHitted.clear();
    fId.clear();
    fFlag1.clear();
    fFlag2.clear();
    fClusters.clear();
        
    return kTRUE;
}

// process
 Bool_t BaseClusteringModule::Process( RecoEvent* ev ) {
    fEv = ev;
    fTotalNumPoints = fEv->GetHeader().GetNumActiveFee();
    

    if ( fTotalNumPoints <= 0 ) {       
        Msg(EsafMsg::Info) << "No points in event! Exit from this event." << MsgDispatch;
        return kFALSE;
    }
    
    // get id of most populous macrocell (main macrocell)
    fMainMacroCellId = fEv->GetMainMacroCell(1);
    if ( fEv->GetRecoCellInfo(fMainMacroCellId) == NULL )
        Msg(EsafMsg::Panic) << "Main macrocell does not correspond to any macrocell simulated in the event." << MsgDispatch;

    // search pixels with num hits > minimum
    for( Int_t i=0; i<fTotalNumPoints; i++ ) {
        if ( fEv->GetRecoPixelData(i)->GetCounts() >= fNumHitsMinimum ) {
            fHitted.push_back(i);
            fNumPoints++;
        }
    }
    Msg(EsafMsg::Debug) << "Num. of pixels with at least " << fNumHitsMinimum << " hits = " << fNumPoints << MsgDispatch;
    if (fNumPoints<10) {
        Msg(EsafMsg::Warning) << "Number of points "<< fNumPoints << " is less then 10 "<< " Exit from this event." << MsgDispatch;
        fHitted.clear();
        fMainMacroCellId = 0;
        fNumPoints = 0;
        return kFALSE;
    }

    // fill id and flag vectors
    for( Int_t i=0; i<fNumPoints; i++ ) {
        fId.push_back(0);
        fFlag1.push_back(kFALSE);
        fFlag2.push_back(kFALSE);
    }
    
    if (Conf()->GetStr("BaseClusteringModule.UseNaturalValues")=="yes") {
        CalculateDensity();
	CalculateNaturalThreshold();
	Msg(EsafMsg::Debug) <<"Density: " << fDensity << "  Natural Distance Threshold: " << fThreshold << MsgDispatch;
	CalculateNumClustersUniform();
	CalculateNumPointsMinimum();
	Msg(EsafMsg::Debug) <<"Number of clusters in uniform distribution: " << fNumClustersUniform << MsgDispatch;
	Msg(EsafMsg::Debug) << "Number of points minimum in a significant cluster: " << fNumPointsMinimum << MsgDispatch;
    } else if (Conf()->GetStr("BaseClusteringModule.UseNaturalValues")=="no"){
	fThreshold = Conf()->GetNum("BaseClusteringModule.Threshold");
	fNumPointsMinimum = (Int_t)Conf()->GetNum("BaseClusteringModule.NumPointsMinimum");
	Msg(EsafMsg::Debug) << "Distance Threshold: " << fThreshold << MsgDispatch;
	Msg(EsafMsg::Debug) << "Number of points minimum in a significant cluster: " << fNumPointsMinimum << MsgDispatch;
    } else {
        Msg(EsafMsg::Panic) << "You must specify yes or no in BaseClusteringModule.UseNaturalValues."<< MsgDispatch;
    }
    SearchClusters();
    Msg(EsafMsg::Debug) << "Number of significant clusters: " << fNumSigClusters << " in " << fNumClusters << " clusters" << MsgDispatch;
    
    if ( fNumSigClusters != 0 ) {
        for( Int_t i=0; i<fNumSigClusters; i++ ) {
	    MsgForm(EsafMsg::Debug, "Cluster n.%d with %d points",i,fClusters[i]->GetNumPoints());
            if (fClusters[i]->GetFitted()) {
	        MsgForm(EsafMsg::Debug,"Slope: %f Offset: %f Absolute deviation %f ",fClusters[i]->GetSlope(),fClusters[i]->GetOffset(),fClusters[i]->GetAbsoluteDeviation());
            } 
        }

        // add data to event
        vector<EusoCluster*> *pclu = &fClusters;
        MyData()->Add("Clusters", pclu);
        MyData()->Add("NumHitsMinimum",fNumHitsMinimum);
        MyData()->Add("SignificanceLevel",fSignificanceLevel);
        MyData()->Add("NumPointsMinimum",fNumPointsMinimum);
        MyData()->Add("DistanceThreshold",fThreshold);
        return kTRUE;
    }
    
    return kFALSE;
}

// post-process
 Bool_t BaseClusteringModule::PostProcess() {
    if ( fTotalNumPoints == 0 ) return kTRUE;
    if ( fNumPoints < 10 ) {
        fEv = (RecoEvent*)0; 
        return kTRUE;
    }
    if ( Conf()->GetStr("BaseClusteringModule.Histogram")=="yes" ) {
        if ( fClusters.size() != 0 ) HistoCluster();
        fEv = (RecoEvent*)0;
        return kTRUE;
    } else if ( Conf()->GetStr("BaseClusteringModule.Histogram")=="no" ) {
        fEv = (RecoEvent*)0;
        return kTRUE;
    } else {
        Msg(EsafMsg::Panic) << "Wrong config value BaseClusteringModule.Histogram" << MsgDispatch;
        return kFALSE;
    }
}

 Bool_t BaseClusteringModule::SaveRootData(RecoRootEvent *fRecoRootEvent) {
     return kTRUE;
}

// done
 Bool_t BaseClusteringModule::Done() {
    fHitted.clear();
    fId.clear();
    fFlag1.clear();
    fFlag2.clear();
    fClusters.clear();
    Msg(EsafMsg::Info) << "Completed" << MsgDispatch;
    return kTRUE;
}

// memory clean
 void BaseClusteringModule::UserMemoryClean() {
    if ( fTotalNumPoints == 0 || fNumPoints < 10 || fNumSigClusters == 0 ) return;

    vector<EusoCluster*> *pclu = (vector<EusoCluster*>*)MyData()->GetObj("Clusters");
    EusoCluster *dummy(NULL);
    for(size_t i(0); i<pclu->size(); i++) {
        dummy = (*pclu)[i];
        delete dummy;
    }
    (*pclu).clear();
    MyData()->RemoveObj("Clusters");
}

// calculate density of points
 void BaseClusteringModule::CalculateDensity() {
    if ( Conf()->GetStr("BaseClusteringModule.DensityMethod")=="mainmc" ) {
        DensityMainMC();
    } else if ( Conf()->GetStr("BaseClusteringModule.DensityMethod")=="event" ) {
        DensityEvent();
    } else {
        Msg(EsafMsg::Panic) <<"Wrong config value BaseClusteringModule.DensityMethod" << MsgDispatch;  
    }
}

// Calculate density from main macrocell
 void BaseClusteringModule::DensityMainMC() {
    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;
    if (fDensity==0) DensityEvent();
}

// Calculate density from complete event
 void BaseClusteringModule::DensityEvent() {
    Float_t thmin = fEv->GetRecoPixelData(fHitted[0])->GetTheta();
    Float_t thmax = fEv->GetRecoPixelData(fHitted[0])->GetTheta();
    Float_t phimin = fEv->GetRecoPixelData(fHitted[0])->GetPhi();
    Float_t phimax = fEv->GetRecoPixelData(fHitted[0])->GetPhi();
    
    for(Int_t i=1; i<fNumPoints; i++) {
        Int_t j = fHitted[i];
        if ( fEv->GetRecoPixelData(j)->GetTheta() < thmin )
            thmin = fEv->GetRecoPixelData(j)->GetTheta();
        if ( fEv->GetRecoPixelData(j)->GetTheta() > thmax )
            thmax = fEv->GetRecoPixelData(j)->GetTheta();
        if ( fEv->GetRecoPixelData(j)->GetPhi() < phimin )
            phimin = fEv->GetRecoPixelData(i)->GetPhi();
        if ( fEv->GetRecoPixelData(j)->GetPhi() > phimax )
            phimax = fEv->GetRecoPixelData(i)->GetPhi();
    }
    Double_t area = (thmax-thmin)*(phimax-phimin);
    fDensity = (Double_t) fTotalNumPoints / area;
}

// calculate natural threshold
 void BaseClusteringModule::CalculateNaturalThreshold() {
    fThreshold = TMath::Sqrt(TMath::Log(2)/
        (fDensity*fAlpha*TMath::Pi()));
}

// calculate distance between two points
 Double_t BaseClusteringModule::CalculateDistance(Int_t n1, Int_t n2) {
    RecoPixelData *px1=fEv->GetRecoPixelData(n1);
    RecoPixelData *px2=fEv->GetRecoPixelData(n2);
    return TMath::ACos( TMath::Sin(px1->GetTheta())*TMath::Sin(px2->GetTheta())*
           TMath::Cos(px1->GetPhi())*TMath::Cos(px2->GetPhi()) +
           TMath::Sin(px1->GetTheta())*TMath::Sin(px2->GetTheta())*
           TMath::Sin(px1->GetPhi())*TMath::Sin(px2->GetPhi()) +
           TMath::Cos(px1->GetTheta())*TMath::Cos(px2->GetTheta()) );
}

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

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

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

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

// fill a cluster
 void BaseClusteringModule::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 ( CalculateDistance(fHitted[i], fHitted[k]) <= fThreshold ) {
                        fFlag2[i] = kTRUE;
                        fNumPointsCluster++;
                        fId[fNumPointsCluster-1] = i;
                        fNumNodes++;
                    }
                    fBookmark = i;
                }
            }
        }
        fCounter++;
    }
}

// close and write the current cluster if is significant
 void BaseClusteringModule::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 );
    FitCluster( cl );
    fClusters.push_back(cl);
}


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

// fit the points of a cluster; use median fit
 void BaseClusteringModule::FitCluster( EusoCluster *clu ) {
    if ( Conf()->GetStr("BaseClusteringModule.FitMethod")=="median" ) {
        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() );
    } else if ( Conf()->GetStr("BaseClusteringModule.FitMethod")=="none" ) {
        clu->SetFitted( kFALSE );
        clu->SetSlope( 0. );
        clu->SetOffset( 0. );
        clu->SetAbsoluteDeviation( 0. );
        return;
    } else {
        Msg(EsafMsg::Panic) << "Wrong config value BaseClusteringModule.FitMethod" << MsgDispatch;  
    }
}

// make an histogram of the cluster
// FIXME
 void BaseClusteringModule::HistoCluster() {

    gDirectory->cd("/");
    string pathname = Form("cluster%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<fTotalNumPoints; 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("BaseClusteringModule.HistogramBinning");
    TH3F *clhisto = new TH3F("ClHisto","Clustered Points",bin,phimin,phimax,bin,thmin,thmax,bin,tmin,tmax);
    TH3F *hithisto = new TH3F("HitHisto","Points with minimum hits",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 *clhisto2 = new TH2F("ClHisto2","Clustered Points",bin,phimin,phimax,bin,thmin,thmax);
    TH2F *hithisto2 = new TH2F("HitHisto2","Points with minimum hits",bin,phimin,phimax,bin,thmin,thmax);
    TH2F *histo2 = new TH2F("Histo2","Points",bin,phimin,phimax,bin,thmin,thmax);
 
 
    for( Int_t i=0; i<fTotalNumPoints; i++ ) {
        histo->Fill( fEv->GetRecoPixelData(i)->GetPhi(),fEv->GetRecoPixelData(i)->GetTheta(), 
                     fEv->GetRecoPixelData(i)->GetGtu() );
        histo2->Fill( fEv->GetRecoPixelData(i)->GetPhi(),fEv->GetRecoPixelData(i)->GetTheta());
    }
    histo->GetXaxis()->SetTitle("#varphi (rad)");
    histo->GetYaxis()->SetTitle("#theta (rad)");
    histo->GetZaxis()->SetTitle("GTU");
    histo2->GetXaxis()->SetTitle("#varphi (rad)");
    histo2->GetYaxis()->SetTitle("#theta (rad)");
    
    for( Int_t i=0; i<fNumPoints; i++ ) {
        hithisto->Fill(fEv->GetRecoPixelData(fHitted[i])->GetPhi(),
            fEv->GetRecoPixelData(fHitted[i])->GetTheta(), fEv->GetRecoPixelData(fHitted[i])->GetGtu() );
        hithisto2->Fill(fEv->GetRecoPixelData(fHitted[i])->GetPhi(),
            fEv->GetRecoPixelData(fHitted[i])->GetTheta());
    }
    hithisto->GetXaxis()->SetTitle("#varphi (rad)");
    hithisto->GetYaxis()->SetTitle("#theta (rad)");
    hithisto->GetZaxis()->SetTitle("GTU");
    hithisto2->GetXaxis()->SetTitle("#varphi (rad)");
    hithisto2->GetYaxis()->SetTitle("#theta (rad)");
    
    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));
            clhisto->Fill( pix->GetPhi(), pix->GetTheta(), pix->GetGtu() );
            clhisto2->Fill( pix->GetPhi(), pix->GetTheta());
        }
    }
    clhisto->GetXaxis()->SetTitle("#varphi (rad)");
    clhisto->GetYaxis()->SetTitle("#theta (rad)");
    clhisto->GetZaxis()->SetTitle("GTU");
    clhisto2->GetXaxis()->SetTitle("#varphi (rad)");
    clhisto2->GetYaxis()->SetTitle("#theta (rad)");
    
    clhisto->Write();
    hithisto->Write();
    histo->Write();
    delete clhisto;
    delete hithisto;
    delete histo;
    clhisto2->Write();
    hithisto2->Write();
    histo2->Write();
    delete clhisto2;
    delete hithisto2;
    delete histo2;
}
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