// $Id: HoughModule.cc
// Author: Elena Taddei
/*****************************************************************************
* ESAF: Euso Simulation and Analysis Framework *
* *
* Id: HoughModule *
* *
*****************************************************************************/
//_____________________________________________________________________________
//
// HoughModule
//
// <extensive class description>
//
// Config file parameters
// ======================
//
// <parameter name>: <parameter description>
// -Valid options: <available options>
//
#include "HoughModule.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 "HoughFit.hh"
#include "ERunParameters.hh"
#include "EDetector.hh"
#include <TH3F.h>
#include <TH2F.h>
#include <TH1F.h>
#include <TGraph.h>
#include <TCanvas.h>
using namespace TMath;
ClassImp(HoughModule)
// FIXME memory leaks to be fixed asap. D.Naumov
//_____________________________________________________________________________
HoughModule::HoughModule():RecoModule("HoughTransform") {
//
// ctor
//
}
//_____________________________________________________________________________
HoughModule::~HoughModule() {
//
// dtor
//
}
//_____________________________________________________________________________
Bool_t HoughModule::Init() {
//
// init
//
fEv = (RecoEvent*)0;
fNumPointsMinimum = (Int_t)Conf()->GetNum("HoughModule.NumPointsMinimum");
fProcedure = Conf()->GetStr("HoughModule.Procedure");
Msg(EsafMsg::Info) << "Initializing... Procedure: ""<< fProcedure.c_str() << """ << MsgDispatch;
return kTRUE;
}
//_____________________________________________________________________________
Bool_t HoughModule::PreProcess() {
//
// pre-process
//
fTotalNumPoints = 0;
fNumPointsAll = 0;
fHittedFinal.clear();
fHittedFinalMacrocell.clear();
fHittedSelHough.clear();
fHittedSel2Hough.clear();
fHittedAll.clear();
fIdMacrocellsSelected.clear();
fClusters.clear();
return kTRUE;
}
// process
//_____________________________________________________________________________
Bool_t HoughModule::Process( RecoEvent* ev ) {
fEv = ev;
fTotalNumPoints = fEv->GetHeader().GetNumActiveFee();
fIdEv = fEv->GetHeader().GetNum();
if ( fTotalNumPoints <= fNumPointsMinimum ) {
Msg(EsafMsg::Info) << "Not enough points in event! Exit from this event." << MsgDispatch;
return kFALSE;
}else{
Msg(EsafMsg::Debug) << "fTotalNumPoints = " << fTotalNumPoints << MsgDispatch;
}
fMainMacroCellId = fEv->GetMainMacroCell(1);
if ( fEv->GetRecoCellInfo(fMainMacroCellId) == NULL )
Msg(EsafMsg::Info) << "Main macrocell does not correspond to any macrocell hit." << MsgDispatch;
Bool_t ProcessOnlySignal = Config::Get()->GetCF("Reco","RootInputModule")->GetBool("RootInputModule.ProcessOnlySignal");
if (ProcessOnlySignal) {
for( Int_t i=0; i<fTotalNumPoints; i++ ) {
RecoPixelData *pix = (RecoPixelData*)fEv->GetRecoPixelData(i);
Int_t signal_counts = pix->GetCountsSignal();
if ( signal_counts > 0 ) fHittedFinal.push_back(i);
}
fFinalNumHitsMinimum=1;
WriteSingleCluster();
return kTRUE;
}
gDirectory->cd("/");
string pathnameevent = Form("HoughModule%d", fIdEv);
TDirectory *path = new TDirectory(pathnameevent.c_str(),pathnameevent.c_str());
path->cd();
Int_t ret;
if( fProcedure == "single" ) ret = ProcessHoughSingleMacrocell();
if( fProcedure == "velocity" ) ret = ProcessVelocity();
if( fProcedure == "all" ) ret = ProcessHoughAll();
if(!ret){
Msg(EsafMsg::Info) << "Problems ret=0.. " << MsgDispatch;
return kFALSE;
}
if(ret==1){
if(fHittedFinal.size()<10 || fHittedFinal.size()>100000) return kFALSE;
WriteSingleCluster();
return kTRUE;
}
return kFALSE;
}
//_____________________________________________________________________________
Int_t HoughModule::ProcessHoughSingleMacrocell(){
//
// Do Hough Transform for each macrocell hitted and trying to select
// which are hitted by real signal
//
Msg(EsafMsg::Info) << " ...ProcessHoughSingleMacrocell()" << MsgDispatch;
Int_t numcountsmin = (Int_t)Conf()->GetNum("HoughModule.MinCountsHoughSingleMacrocell");
RecoEventHeader header = fEv->GetHeader();
ERunParameters *runpars = header.GetRunPars();
map< Int_t,vector<Int_t> > mapIdfHittedSingleCell;
for(Int_t i=0;i<200;i++){
RecoCellInfo *info = fEv->GetRecoCellInfo(i);
if( info != NULL){
vector<Int_t> empty;
pair< Int_t,vector<Int_t> > p(i,empty);
mapIdfHittedSingleCell.insert(p);
}
}
Msg(EsafMsg::Info) << "Number of macrocells hitted = " << mapIdfHittedSingleCell.size() << MsgDispatch;
map< Int_t,vector<Int_t> >::iterator it;
for( Int_t i=0; i<fTotalNumPoints; i++ ){
if( fEv->GetRecoPixelData(i)->GetCounts() > numcountsmin ){
Int_t cellid=fEv->GetRecoPixelData(i)->GetMacroCellId();
it = mapIdfHittedSingleCell.find(cellid);
it->second.push_back(i);
}
}
Msg(EsafMsg::Info) << "selected pixel with counts > " << numcountsmin << MsgDispatch;
for( it = mapIdfHittedSingleCell.begin(); it != mapIdfHittedSingleCell.end(); it++ ){
Msg(EsafMsg::Info)<<"n\tt---->macrocellid = "<<it->first<<"t with "<<(it->second).size()<<" points "<<MsgDispatch;
if((it->second).size()<10) continue;
Int_t idmacrocell=it->first;
vector<Int_t> fHittedmc=it->second;
vector<Double_t> xmc;
vector<Double_t> ymc;
vector<Double_t> tmc;
vector<Int_t> cmc;
TH1F *h1 = new TH1F("h1","h1",10,0,10);
for(Int_t i=0;i<(Int_t)fHittedmc.size();i++){
h1->Fill(fEv->GetRecoPixelData(fHittedmc[i])->GetCounts());
Int_t pixid=fEv->GetRecoPixelData(fHittedmc[i])->GetPixelId();
TVector3 dummy2 = runpars->PixelCenter(pixid);
xmc.push_back(0.001*dummy2.x());//m
ymc.push_back(0.001*dummy2.y());//m
Int_t tmcgtu=fEv->GetRecoPixelData(fHittedmc[i])->GetGtu();
tmc.push_back(tmcgtu*2.5*0.01);
cmc.push_back(fEv->GetRecoPixelData(fHittedmc[i])->GetCounts());
}
string pathname = Form("/HoughModule%d/", fIdEv);
gDirectory->cd(pathname.c_str());
string pathnames = Form("SingleHoughMacrocell%d", idmacrocell);
TDirectory *path = new TDirectory(pathnames.c_str(),pathnames.c_str());
path->cd();
h1->Draw();
h1->Write();
delete h1;
Bool_t possiblesignal = HoughTransform(xmc,ymc,tmc,cmc);
if( possiblesignal ) {
Msg(EsafMsg::Info) << "Macrocell " << idmacrocell << " probably HAS signal!!!" << MsgDispatch;
fIdMacrocellsSelected.push_back(idmacrocell);
fHittedFinalMacrocell.insert(fHittedFinalMacrocell.end(),fHittedmc.begin(),fHittedmc.end());
}else{
Msg(EsafMsg::Info) << "Macrocell " << idmacrocell << " probably DOESN'T HAVE a signal,";
if ( idmacrocell==fMainMacroCellId ){
Msg(EsafMsg::Info) << " but is added because is the main!" << MsgDispatch;
fHittedFinalMacrocell.insert(fHittedFinalMacrocell.end(),fHittedmc.begin(),fHittedmc.end());
fIdMacrocellsSelected.push_back(idmacrocell);
} else {
Msg(EsafMsg::Info) << " discarded!" << MsgDispatch;
}
}
xmc.clear();
ymc.clear();
tmc.clear();
cmc.clear();
}
if (fIdMacrocellsSelected.size()<1 || fHittedFinalMacrocell.size()<10) return 0;
Msg(EsafMsg::Info) << "Re-process " << fIdMacrocellsSelected.size() << " selected macrocells by ProcessHoughSingleMacrocell()" << MsgDispatch;
Int_t ret = ProcessHoughAll();
return ret;
}
//_____________________________________________________________________________
Int_t HoughModule::ProcessVelocity(){
//
// Do histograms of velocities for each macrocell hitted and trying to select
// which are hitted by real signal
//
Msg(EsafMsg::Info) << " ...ProcessVelocity()" << MsgDispatch;
Int_t numcountsmin = (Int_t)Conf()->GetNum("HoughModule.MinCountsVelocity");
RecoEventHeader header = fEv->GetHeader();
ERunParameters *runpars = header.GetRunPars();
map< Int_t,vector<Int_t> > mapIdfHittedSingleCell;
for(Int_t i=0;i<200;i++){
RecoCellInfo *info = fEv->GetRecoCellInfo(i);
if( info != NULL){
vector<Int_t> empty;
pair< Int_t,vector<Int_t> > p(i,empty);
mapIdfHittedSingleCell.insert(p);
}
}
Msg(EsafMsg::Info) << "Number of macrocells hitted = " << mapIdfHittedSingleCell.size() << MsgDispatch;
map< Int_t,vector<Int_t> >::iterator it;
for( Int_t i=0; i<fTotalNumPoints; i++ ){
if( fEv->GetRecoPixelData(i)->GetCounts() > numcountsmin ){
Int_t cellid=fEv->GetRecoPixelData(i)->GetMacroCellId();
it = mapIdfHittedSingleCell.find(cellid);
it->second.push_back(i);
}
}
Int_t binv = 20;
Msg(EsafMsg::Info) << "selected pixel with counts > " << numcountsmin << MsgDispatch;
for( it = mapIdfHittedSingleCell.begin(); it != mapIdfHittedSingleCell.end(); it++ ) {
Int_t entriesTH2(0);
Msg(EsafMsg::Info) << "macrocellid = " << it->first << "t with " << (it->second).size() << " points" << MsgDispatch;
if((it->second).size()<10) continue;
TH2F *hsxt = new TH2F("hsxt","hsxt;t(mus);x(mm)",1000,0,200,1000,-1000,1000);
TH2F *hsyt = new TH2F("hsyt","hsyt;t(mus);y(mm)",1000,0,200,1000,-1000,1000);
TH2F *s = new TH2F("s","s;velocita;coeffang",200,0,2,200,-4,4);
TH1F *sx = new TH1F("sx","sx;velocita",binv,-2,2);
TH1F *sy = new TH1F("sy","sy;velocita",binv,-2,2);
TH2F *sxy = new TH2F("sxy","velocita",binv,-2,2,binv,-2,2);
Int_t idmacrocell = it->first;
vector<Int_t> fHittedmc = it->second;
for(Int_t i=0;i<(Int_t)fHittedmc.size();i++){
Int_t pixid = fEv->GetRecoPixelData(fHittedmc[i])->GetPixelId();
TVector3 dummy2 = runpars->PixelCenter(pixid);//mm
Int_t tmcgtu = fEv->GetRecoPixelData(fHittedmc[i])->GetGtu();
Int_t counts = fEv->GetRecoPixelData(fHittedmc[i])->GetCounts();
hsxt->Fill(tmcgtu*2.5,dummy2.x(),counts);
hsyt->Fill(tmcgtu*2.5,dummy2.y(),counts);
Int_t pointstocouple = fHittedmc.size();
for(Int_t j=i+1;j<pointstocouple;j++){
Int_t pixtocouple = fEv->GetRecoPixelData(fHittedmc[j])->GetPixelId();
TVector3 dummycouple = runpars->PixelCenter(pixtocouple);
Int_t tmccouple = fEv->GetRecoPixelData(fHittedmc[j])->GetGtu();
if (Abs(dummy2.x()-dummycouple.x())>0 && Abs(dummy2.y()-dummycouple.y())>0 && tmccouple-tmcgtu>0 && tmccouple-tmcgtu<20 && Abs(dummy2.x()-dummycouple.x())<100 && Abs(dummy2.y()-dummycouple.y())<100 ) {
Int_t countproduct = (Int_t)Sqrt(1.e0*counts*fEv->GetRecoPixelData(fHittedmc[j])->GetCounts());
Double_t m = (dummy2.y()-dummycouple.y())/(dummy2.x()-dummycouple.x());
Double_t v = Sqrt((dummy2.y()-dummycouple.y())*(dummy2.y()-dummycouple.y())+(dummy2.x()-dummycouple.x())*(dummy2.x()-dummycouple.x()))/((tmccouple-tmcgtu)*2.5);// mm/us
Double_t vx = (dummy2.x()-dummycouple.x())/((tmcgtu-tmccouple)*2.5);
Double_t vy = (dummy2.y()-dummycouple.y())/((tmcgtu-tmccouple)*2.5);
entriesTH2++;
s->Fill(v,m,countproduct);
sx->Fill(vx,countproduct);
sy->Fill(vy,countproduct);
sxy->Fill(vx,vy,countproduct);
}
}
}
Int_t binmaxv = sxy->GetMaximumBin();
Int_t nxv = (sxy->GetXaxis())->GetNbins()+2;
Int_t binxmaxv = binmaxv % nxv;
Int_t binymaxv = binmaxv / nxv;
Float_t vxmax[1], vymax[1];
vxmax[0] = (sxy->GetXaxis())->GetBinCenter(binxmaxv);
vymax[0] = (sxy->GetYaxis())->GetBinCenter(binymaxv);
TGraph *gvmax=new TGraph(1,vxmax,vymax);
gvmax->SetMarkerColor(2); gvmax->SetMarkerStyle(29); gvmax->SetMarkerSize(1.6);
//Float_t contmaxv=sxy->GetBinContent(binxmaxv,binymaxv);
Stat_t countx(0),county(0);
for(Int_t bin=1;bin<binv;bin++){
countx += sx->GetBinContent(bin);
county += sy->GetBinContent(bin);
}
Float_t countmaxx = sx->GetMaximum();
Float_t countmaxy = sy->GetMaximum();
Float_t countmeanx = (Float_t)((Float_t)countx/binv);
Float_t countmeany = (Float_t)((Float_t)county/binv);
Float_t ms = 20.;
/*cout<<"countmeanx="<<countmeanx<<"tcontmaxx="<<countmaxx<<"tcontmaxx-countmeanx="<<countmaxx-countmeanx<<"tsig="<<ms*TMath::Sqrt(countmeanx)<<endl;
cout<<"countmeany="<<countmeany<<"tcontmaxy="<<countmaxy<<"tcontmaxy-countmeany="<<countmaxy-countmeany<<"tsig="<<ms*TMath::Sqrt(countmeany)<<endl;*/
string a = "a"; string b = "b";
string pathname = Form("/HoughModule%d/", fIdEv);
gDirectory->cd(pathname.c_str());
string pathnames;
if( ((countmaxx-countmeanx) > ms*Sqrt(countmeanx) || (countmaxy-countmeany) > ms*Sqrt(countmeany)) && sx->GetEntries()>10 ) {
fHittedFinalMacrocell.insert(fHittedFinalMacrocell.end(),fHittedmc.begin(),fHittedmc.end());
fIdMacrocellsSelected.push_back(idmacrocell);
Msg(EsafMsg::Info) << "Macrocell " << idmacrocell << " probably HAS signal!!!" << MsgDispatch;
if (idmacrocell==fMainMacroCellId) {
pathnames = Form("VelocityMacrocellmain%d", idmacrocell);
} else {
pathnames = Form("VelocityMacrocell%d", idmacrocell);
}
} else{
Msg(EsafMsg::Info) << "Macrocell " << idmacrocell << " probably DOESN'T HAVE a signal,";
if (idmacrocell==fMainMacroCellId) {
pathnames = Form("VelocityMacrocellmain%d", idmacrocell);
Msg(EsafMsg::Info) <<" but is added because is the main!"<< MsgDispatch;
fHittedFinalMacrocell.insert(fHittedFinalMacrocell.end(),fHittedmc.begin(),fHittedmc.end());
fIdMacrocellsSelected.push_back(idmacrocell);
} else {
pathnames = Form("VelocityMacrocelldiscarded%d", idmacrocell);
Msg(EsafMsg::Info) <<" discarded!"<< MsgDispatch;
}
}
s->SetTitle(pathnames.c_str());
TDirectory *path = new TDirectory(pathnames.c_str(),pathnames.c_str());
path->cd();
TCanvas *ac=new TCanvas(a.c_str(),a.c_str(),1);ac->Divide(2,2);
ac->cd(1);s->Draw(); ac->cd(2);s->Draw("LEGO");
ac->cd(3);sx->Draw(); ac->cd(4);sy->Draw();
ac->Write();
delete ac; delete s; delete sx; delete sy;
TCanvas *bc=new TCanvas(b.c_str(),b.c_str(),1);bc->Divide(2,2);
bc->cd(1);hsxt->Draw();
bc->cd(2);hsyt->Draw();
bc->cd(3);sxy->Draw();gvmax->Draw("same");
bc->cd(4);sxy->Draw("LEGO");
bc->Write();
delete bc; delete hsxt; delete hsyt;
delete sxy;delete gvmax;
}
if (fIdMacrocellsSelected.size()<1 || fHittedFinalMacrocell.size()<10) return 0;
Msg(EsafMsg::Info) << "Re-process " << fIdMacrocellsSelected.size() << " selected macrocells by ProcessVelocity()" << MsgDispatch;
Int_t ret = ProcessHoughAll();
return ret;
}
//_____________________________________________________________________________
Int_t HoughModule::ProcessHoughAll(){
//
// Doing Hough Transform to all macrocells selected by previous
// methods or by euso trigger
//
Msg(EsafMsg::Info) << " ...ProcessHoughAll()" << MsgDispatch;
string pathname = Form("/HoughModule%d/", fIdEv);
gDirectory->cd(pathname.c_str());
string pathnametotal = Form("HoughAll");
TDirectory *pathtotal = new TDirectory(pathnametotal.c_str(),pathnametotal.c_str());
pathtotal->cd();
RecoEventHeader header = fEv->GetHeader();
ERunParameters *runpars = header.GetRunPars();
fNumHitsMinimum = (Int_t)Conf()->GetNum("HoughModule.MinCountsHoughAll");
if ( fNumHitsMinimum <= 0 ){
Msg(EsafMsg::Panic) << "You must specify a positive number in HoughModule.MinCountsHoughAll" << MsgDispatch;
}else{
Msg(EsafMsg::Info) << "Let's try beginning with MinCountsHoughAll = " << fNumHitsMinimum << MsgDispatch;
}
Int_t regr=0;
Float_t gtumin = 1000.;
Float_t gtumax = 0.;
Float_t deltagtu = 0.;
fNumPointsAll = 0;
fHittedAll.clear();
if(!(fHittedFinalMacrocell.size())){
Msg(EsafMsg::Info) << "processing all macrocells!" << MsgDispatch;
for( Int_t i=0; i<fTotalNumPoints; i++ ) {
if (fEv->GetRecoPixelData(i)->GetGtu()<gtumin) gtumin=fEv->GetRecoPixelData(i)->GetGtu();
if (fEv->GetRecoPixelData(i)->GetGtu()>gtumax) gtumax=fEv->GetRecoPixelData(i)->GetGtu();
if (fEv->GetRecoPixelData(i)->GetCounts() >= fNumHitsMinimum) fHittedAll.push_back(i);
}
} else {
Msg(EsafMsg::Info)<<"processing fHittedFinalMacrocell.size()="<<fHittedFinalMacrocell.size()<<MsgDispatch;
for( Int_t i=0; i<(Int_t)(fHittedFinalMacrocell.size()); i++ ){
Int_t i2=fHittedFinalMacrocell[i];
if (fEv->GetRecoPixelData(i2)->GetGtu()<gtumin) gtumin=fEv->GetRecoPixelData(i2)->GetGtu();
if (fEv->GetRecoPixelData(i2)->GetGtu()>gtumax) gtumax=fEv->GetRecoPixelData(i2)->GetGtu();
if (fEv->GetRecoPixelData(i2)->GetCounts() >= fNumHitsMinimum) fHittedAll.push_back(i2);
}
}
fNumPointsAll = fHittedAll.size();
fTmin = gtumin*2.5*0.01; //FIXME
fTmax = gtumax*2.5*0.01; //FIXME
deltagtu=gtumax-gtumin;
Msg(EsafMsg::Info) << "Num. of pixels with at least " << fNumHitsMinimum << " hits = " << fNumPointsAll << " over fTotalNumPoints=" << fTotalNumPoints << MsgDispatch;
if (fNumPointsAll<10 || deltagtu<5) {
if( fNumPointsAll<10 ) Msg(EsafMsg::Warning) << "Number of points "<< fNumPointsAll << " is less then 10, Exit from this event." << MsgDispatch;
if ( deltagtu<5 ) Msg(EsafMsg::Warning) << " deltagtu = "<< deltagtu << " is less then 5, Exit from this event." << MsgDispatch;
fHittedAll.clear();
fMainMacroCellId = 0;
fNumPointsAll = 0;
return 0;
}
vector<Double_t> x;
vector<Double_t> y;
vector<Double_t> t;
vector<Int_t> c;
Int_t bin=2000;
TH2F *single2 = new TH2F("single2","single2;x(m);y(m)",bin,-1,1,bin,-1,1);
TH2F *single2xt = new TH2F("single2xt","single2xt;t;x(m)",bin,fTmin,fTmax,bin,-1,1);
TH2F *single2yt = new TH2F("single2yt","single2yt;t;y(m)",bin,fTmin,fTmax,bin,-1,1);
for( Int_t i=0; i < fNumPointsAll; i++ ) {
Int_t pixid=fEv->GetRecoPixelData(fHittedAll[i])->GetPixelId();
TVector3 dummy2 = runpars->PixelCenter(pixid);
x.push_back(0.001*dummy2.x());
y.push_back(0.001*dummy2.y());
t.push_back((Double_t)fEv->GetRecoPixelData(fHittedAll[i])->GetGtu()*2.5*0.01);
c.push_back(fEv->GetRecoPixelData(fHittedAll[i])->GetCounts());
single2->Fill(x[i],y[i],c[i]);
single2xt->Fill(t[i],x[i],c[i]);
single2yt->Fill(t[i],y[i],c[i]);
}
TCanvas *HMinit=new TCanvas("HMinit","HMinit",1);
HMinit->Divide(3,1);
HMinit->cd(1); single2->SetMarkerStyle(6); single2->Draw();
HMinit->cd(2); single2xt->SetMarkerStyle(6); single2xt->Draw();
HMinit->cd(3); single2yt->SetMarkerStyle(6); single2yt->Draw();
HMinit->Write();
delete HMinit;
delete single2xt;
delete single2yt;
delete single2;
HoughTransform(x,y,t,c);
x.clear(); y.clear(); t.clear(); c.clear();
if(fNumHitsMinimum-regr > 3) regr++;
vector<Double_t> xsel;
vector<Double_t> ysel;
vector<Double_t> tsel;
vector<Int_t> csel;
TH2F *single2sel = new TH2F("single2sel","single2sel;x(m);y(m)",bin,fXmin,fXmax,bin,fYmin,fYmax);
TH2F *single2xtsel = new TH2F("single2xtsel","single2xtsel;t;x(m)",bin,fTmin,fTmax,bin,fXmin,fXmax);
TH2F *single2ytsel = new TH2F("single2ytsel","single2ytsel;t;y(m)",bin,fTmin,fTmax,bin,fYmin,fYmax);
for( Int_t i=0; i<fTotalNumPoints; i++ ) {
Int_t pixid=fEv->GetRecoPixelData(i)->GetPixelId();
TVector3 dummy2 = runpars->PixelCenter(pixid);
Double_t x=0.001*dummy2.x();
Double_t y=0.001*dummy2.y();
Double_t t=(Double_t)fEv->GetRecoPixelData(i)->GetGtu()*2.5*0.01;
Int_t c=fEv->GetRecoPixelData(i)->GetCounts();
if ( x>=fXmin && x<=fXmax && y>=fYmin && y<=fYmax && c>=(fNumHitsMinimum-regr) ) {
single2sel->Fill(x,y,c);
single2xtsel->Fill(t,x,c);
single2ytsel->Fill(t,y,c);
xsel.push_back(x);
ysel.push_back(y);
tsel.push_back(t);
csel.push_back(c);
fHittedSelHough.push_back(i);
}
}
Msg(EsafMsg::Info) << "first selected zone: pixels with at least " << fNumHitsMinimum-regr << " hits = " << fHittedSelHough.size() << MsgDispatch;
if(fHittedSelHough.size()<10){
Msg(EsafMsg::Info) << "too few! exit from this event" <<MsgDispatch;
return kFALSE;
}
TCanvas *HMsel=new TCanvas("HMsel","HMsel",1);
HMsel->Divide(3,1);
HMsel->cd(1); single2sel->SetMarkerStyle(6); single2sel->Draw();
HMsel->cd(2); single2xtsel->SetMarkerStyle(6); single2xtsel->Draw();
HMsel->cd(3); single2ytsel->SetMarkerStyle(6); single2ytsel->Draw();
HMsel->Write();
delete HMsel;
delete single2xtsel;
delete single2ytsel;
delete single2sel;
HoughTransform(xsel,ysel,tsel,csel);
xsel.clear(); ysel.clear(); tsel.clear(); csel.clear();
if(fNumHitsMinimum-regr > 2) regr++;
vector<Double_t> xselsel;
vector<Double_t> yselsel;
vector<Double_t> tselsel;
vector<Int_t> cselsel;
TH2F *single2selsel = new TH2F("single2selsel","single2selsel;x(m);y(m)",bin,fXmin,fXmax,bin,fYmin,fYmax);
TH2F *single2xtselsel = new TH2F("single2xtselsel","single2xtselsel;t;x(m)",bin,fTmin,fTmax,bin,fXmin,fXmax);
TH2F *single2ytselsel = new TH2F("single2ytselsel","single2ytselsel;t;y(m)",bin,fTmin,fTmax,bin,fYmin,fYmax);
for( Int_t i=0; i<fTotalNumPoints; i++ ) {
Int_t pixid=fEv->GetRecoPixelData(i)->GetPixelId();
TVector3 dummy2 = runpars->PixelCenter(pixid);
Double_t x=0.001*dummy2.x();
Double_t y=0.001*dummy2.y();
Double_t t=(Double_t)fEv->GetRecoPixelData(i)->GetGtu()*2.5*0.01;
Int_t c=fEv->GetRecoPixelData(i)->GetCounts();
if ( x>=fXmin && x<=fXmax && y>=fYmin && y<=fYmax && c>=(fNumHitsMinimum-regr) ) {
single2selsel->Fill(x,y,c);
single2xtselsel->Fill(t,x,c);
single2ytselsel->Fill(t,y,c);
xselsel.push_back(x);
yselsel.push_back(y);
tselsel.push_back(t);
cselsel.push_back(c);
fHittedSel2Hough.push_back(i);
}
}
Msg(EsafMsg::Info) << "second selected zone: pixels with at least " << fNumHitsMinimum-regr << " hits = " << fHittedSel2Hough.size()<< MsgDispatch;
if(fHittedSel2Hough.size()<10){
Msg(EsafMsg::Info) << "too few! exit from this event" <<MsgDispatch;
return kFALSE;
}
TCanvas *HMselsel=new TCanvas("HMselsel","HMselsel",1);
HMselsel->Divide(3,1);
HMselsel->cd(1); single2selsel->SetMarkerStyle(6); single2selsel->Draw();
HMselsel->cd(2); single2xtselsel->SetMarkerStyle(6); single2xtselsel->Draw();
HMselsel->cd(3); single2ytselsel->SetMarkerStyle(6); single2ytselsel->Draw();
HMselsel->Write();
delete HMselsel;
delete single2xtselsel;
delete single2ytselsel;
delete single2selsel;
HoughTransform(xselsel,yselsel,tselsel,cselsel,fHittedSel2Hough);
fHittedFinal=fHittedHoughMerged;
xselsel.clear(); yselsel.clear(); tselsel.clear(); cselsel.clear();
vector<Double_t> xend;
vector<Double_t> yend;
vector<Double_t> tend;
vector<Int_t> cend;
TH2F *single2end = new TH2F("single2end","single2end;x(m);y(m)",bin,fXmin,fXmax,bin,fYmin,fYmax);
TH2F *single2xtend = new TH2F("single2xtend","single2xtend;t;x(m)",bin,fTmin,fTmax,bin,fXmin,fXmax);
TH2F *single2ytend = new TH2F("single2ytend","single2ytend;t;y(m)",bin,fTmin,fTmax,bin,fYmin,fYmax);
for( Int_t i=0; i<(Int_t)fHittedHoughMerged.size(); i++ ) {
Int_t pixid=fEv->GetRecoPixelData(fHittedHoughMerged[i])->GetPixelId();
TVector3 dummy2 = runpars->PixelCenter(pixid);
Double_t x=0.001*dummy2.x();
Double_t y=0.001*dummy2.y();
Double_t t=(Double_t)fEv->GetRecoPixelData(fHittedHoughMerged[i])->GetGtu()*2.5*0.01;
Int_t c=fEv->GetRecoPixelData(fHittedHoughMerged[i])->GetCounts();
single2end->Fill(x,y,c);
single2xtend->Fill(t,x,c);
single2ytend->Fill(t,y,c);
xend.push_back(x);
yend.push_back(y);
tend.push_back(t);
cend.push_back(c);
}
Msg(EsafMsg::Info) << "final selected zone: pixels with at least " << fNumHitsMinimum-regr << " hits = " << fHittedHoughMerged.size()<< MsgDispatch;
if(fHittedHoughMerged.size()<10){
Msg(EsafMsg::Info) << "too few! exit from this event" <<MsgDispatch;
return kFALSE;
}
TCanvas *HMfinal=new TCanvas("HMfinal","HMfinal",1);
HMfinal->Divide(3,1);
HMfinal->cd(1); single2end->SetMarkerStyle(6); single2end->Draw();
HMfinal->cd(2); single2xtend->SetMarkerStyle(6); single2xtend->Draw();
HMfinal->cd(3); single2ytend->SetMarkerStyle(6); single2ytend->Draw();
HMfinal->Write();
delete HMfinal;
delete single2xtend;
delete single2ytend;
delete single2end;
xend.clear();yend.clear();tend.clear();cend.clear();
fFinalNumHitsMinimum=fNumHitsMinimum-regr;
return 1;
}
//_____________________________________________________________________________
Bool_t HoughModule::HoughTransform(vector<Double_t> x,vector<Double_t> y,vector<Double_t> t,vector<Int_t> c,vector<Int_t> id){
//
// Implements the Hough Transform
// Return true if there is a possible signal
//
if( x.size()<10 ) return 0;
Int_t optionselection=1;
Float_t ex=0.004;
Float_t ey=0.004;
Float_t et=2.5*0.01;
HoughFit *xtFit = new HoughFit(t, x, c, optionselection,et,ex,id);
//Double_t vx = xtFit->GetSlope();//Double_t wx = xtFit->GetOffset();
fXmin = xtFit->GetMinY();
fXmax = xtFit->GetMaxY();
Double_t tmin1 = xtFit->GetMinX();
Double_t tmax1 = xtFit->GetMaxX();
fAbsDevX = xtFit->GetAbsDev();
Bool_t possiblesignalxt = xtFit->IsSignal();
fCmaxToCmeanX = xtFit->GetMaxToMedium();
vector<Int_t> idselx = xtFit->GetIdSel();
delete xtFit;
HoughFit *ytFit = new HoughFit(t, y, c, optionselection,et,ey,id);
//Double_t vy = ytFit->GetSlope();//Double_t wy = ytFit->GetOffset();
fYmin = ytFit->GetMinY();
fYmax = ytFit->GetMaxY();
Double_t tmin2 = ytFit->GetMinX();
Double_t tmax2 = ytFit->GetMaxX();
fAbsDevY = ytFit->GetAbsDev();
Bool_t possiblesignalyt = ytFit->IsSignal();
fCmaxToCmeanY = ytFit->GetMaxToMedium();
vector<Int_t> idsely = ytFit->GetIdSel();
delete ytFit;
if(id.size()>0) {
Int_t sumsel=idsely.size()+idselx.size();
//cout<<"id.size()="<<id.size()<<"tidselx.size()="<<idselx.size()<<"tidsely.size()="<<idsely.size()<<endl;
vector<Int_t> idselectedmerged(sumsel);
merge(idselx.begin(), idselx.end(),idsely.begin(), idsely.end(),idselectedmerged.begin());
vector<Int_t>::iterator diversi=unique(idselectedmerged.begin(),idselectedmerged.end());
idselectedmerged.erase(diversi,idselectedmerged.end());
//cout<<"idselectedmerged.size()="<<idselectedmerged.size()<<endl;
fHittedHoughMerged=idselectedmerged;
}
//Msg(EsafMsg::Info) <<"possiblesignalxt="<<possiblesignalxt<<"tpossiblesignalyt="<<possiblesignalyt<<MsgDispatch;
//Msg(EsafMsg::Info) <<"fCmaxToCmeanX="<<fCmaxToCmeanX<<"tcmaxtocmedioy="<<fCmaxToCmeanY<<MsgDispatch;
Bool_t possible = kFALSE;
if( possiblesignalxt || possiblesignalyt ) possible = kTRUE;
fTmin = (tmin1<tmin2) ? tmin1 : tmin2;
fTmax = (tmax1<tmax2) ? tmax2 : tmax1;
x.clear(); y.clear(); t.clear(); c.clear();
return possible;
}
//_____________________________________________________________________________
Bool_t HoughModule::PostProcess() {
//
// post-process
//
if ( !fTotalNumPoints ) return kTRUE;
fEv = (RecoEvent*)0;
return kTRUE;
}
//____________________________________________________________________________
Bool_t HoughModule::SaveRootData(RecoRootEvent *fRecoRootEvent) {
//
// Save reco root data
//
return kTRUE;
}
//____________________________________________________________________________
Bool_t HoughModule::Done() {
//
// done
//
fHittedAll.clear();
fHittedFinal.clear();
fHittedSelHough.clear();
fHittedSel2Hough.clear();
fIdMacrocellsSelected.clear();
fClusters.clear();
Msg(EsafMsg::Info) << "Completed" << MsgDispatch;
return kTRUE;
}
//_____________________________________________________________________________
void HoughModule::UserMemoryClean() {
//
// memory clean
//
Msg(EsafMsg::Debug) << "UserMemoryClean()" << MsgDispatch;
Int_t fHittedFinalsize = fHittedFinal.size();
if(fHittedFinalsize<10 || !fTotalNumPoints || fHittedFinalsize>100000) 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");
vector<Int_t> *pix = (vector<Int_t>*)MyData()->GetObj("CluPixels");
pix->clear();
MyData()->RemoveObj("CluPixels");
}
//_____________________________________________________________________________
void HoughModule::WriteSingleCluster() {
//
// close and write the current cluster if is significant
//
// write the cluster
EusoCluster *cl = new EusoCluster();
Int_t fHittedFinalsize = fHittedFinal.size();
for( Int_t i=0; i<fHittedFinalsize; i++ ) {
cl->AddPoint( fHittedFinal[i] );
}
fClusters.push_back(cl);
Msg(EsafMsg::Info) << "SingleCluster saved with "<<fHittedFinalsize<< " points" << MsgDispatch;
// add data to event
vector<EusoCluster*> *pclu = &fClusters;
MyData()->Add("Clusters", pclu);
MyData()->Add("NumHitsMinimum",fFinalNumHitsMinimum);
vector<Int_t> *pix = &fHittedFinal;
MyData()->Add("CluPixels", pix);
}
//_____________________________________________________________________________
EusoCluster* HoughModule::GetCluster( Int_t i ) {
//
// get a specified cluster
//
if ( i >= 0 && i < (Int_t)fClusters.size() )
return fClusters[i];
else
return NULL;
}