// ESAF : Euso Simulation and Analysis Framework
// $Id: BaseClusteringModule.cc,v 1.33 2005/01/18 09:07:31 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::Info) << "main macrocell does not correspond to any macrocell hit." << 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;
}