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