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