#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>
ClassImp(HoughModule)
using namespace TMath;
// ctor
HoughModule::HoughModule() :
RecoModule(string("HoughTransform")) {
}
// dtor
HoughModule::~HoughModule() {
}
// init
Bool_t HoughModule::Init() {
fEv = (RecoEvent*)0;
fNumPointsMinimum = (Int_t)Conf()->GetNum("HoughModule.NumPointsMinimum");
fProcedure = Conf()->GetStr("HoughModule.Procedure");
fKindHough = Conf()->GetStr("HoughModule.KindHough");
Msg(EsafMsg::Info) << "Initializing... Procedure: ""<< fProcedure.c_str() << "", KindHough: "" <<fKindHough.c_str() << """<<MsgDispatch;
return kTRUE;
}
// pre-process
Bool_t HoughModule::PreProcess() {
fTotalNumPoints = 0;
fNumPointsall = 0;
fHittedFinal.clear();
fHittedFinalMacrocell.clear();
fHittedselhough.clear();
fHittedselselhough.clear();
fHittedall.clear();
idmacrocellselected.clear();
fClusters.clear();
return kTRUE;
}
// process
Bool_t HoughModule::Process( RecoEvent* ev ) {
fEv = ev;
fTotalNumPoints = fEv->GetHeader().GetNumActiveFee();
Idev = fEv->GetHeader().GetNum();
if ( fTotalNumPoints <= fNumPointsMinimum ) {
Msg(EsafMsg::Info) << "Not enough points in event! Exit from this event." << MsgDispatch;
return kFALSE;
}
fMainMacroCellId = fEv->GetMainMacroCell(1);
if ( fEv->GetRecoCellInfo(fMainMacroCellId) == NULL )
Msg(EsafMsg::Debug) << "main macrocell does not correspond to any macrocell hit." << MsgDispatch;
gDirectory->cd("/");
string pathnameevent = Form("HoughModule%d", Idev);
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==0) return kFALSE;
if(ret==1){
if(fHittedFinal.size()<10 || fHittedFinal.size()>10000) return kFALSE;
WriteSingleCluster();
// add data to event
vector<EusoCluster*> *pclu = &fClusters;
MyData()->Add("Clusters", pclu);
MyData()->Add("NumHitsMinimum",fFinalNumHitsMinimum);
return kTRUE;
}
return kFALSE;
}
////////////////////////////////////////////////////////
//Doing Hough Transform for each macrocell hitted and trying to select which are hitted by real signal
Int_t HoughModule::ProcessHoughSingleMacrocell(){
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> vuoto;
pair< Int_t,vector<Int_t> > p(i,vuoto);
mapIdfHittedSingleCell.insert(p);
}
}
Msg(EsafMsg::Debug) << "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 mincountsmacrocell=100000;
Int_t idmincountsmacrocell(0);
vector<Int_t> fHittedmcmin;
for( it = mapIdfHittedSingleCell.begin(); it != mapIdfHittedSingleCell.end(); it++ ){
if((Int_t)(it->second).size()<mincountsmacrocell){
mincountsmacrocell=(it->second).size();
idmincountsmacrocell=it->first;
fHittedmcmin=it->second;
}
}
Msg(EsafMsg::Debug)<<" n mincountsmacrocell="<<mincountsmacrocell<<"tidmincountsmacrocell="<<idmincountsmacrocell<<MsgDispatch;
vector<Double_t> xmcmin;
vector<Double_t> ymcmin;
vector<Double_t> tmcmin;
vector<Int_t> cmcmin;
for(Int_t i=0;i<(Int_t)fHittedmcmin.size();i++){
Int_t pixid=fEv->GetRecoPixelData(fHittedmcmin[i])->GetPixelId();
TVector3 dummy2 = runpars->PixelCenter(pixid);
xmcmin.push_back(dummy2.x()/1000);//m
ymcmin.push_back(dummy2.y()/1000);//m
tmcmin.push_back((Double_t)fEv->GetRecoPixelData(fHittedmcmin[i])->GetGtu()*2.5*0.01);
cmcmin.push_back(fEv->GetRecoPixelData(fHittedmcmin[i])->GetCounts());
}
HoughTransform(xmcmin,ymcmin,tmcmin,cmcmin);
cmaxtocmediominx=cmaxtocmediox;
cmaxtocmediominy=cmaxtocmedioy;
/////////////////////////////////////////////////////////////////////////////
Msg(EsafMsg::Debug)<<"selected pixel with counts > "<<numcountsmin<<MsgDispatch;
for( it = mapIdfHittedSingleCell.begin(); it != mapIdfHittedSingleCell.end(); it++ ){
Msg(EsafMsg::Debug)<<"nmacrocellid = "<<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;
for(Int_t i=0;i<(Int_t)fHittedmc.size();i++){
Int_t pixid=fEv->GetRecoPixelData(fHittedmc[i])->GetPixelId();
TVector3 dummy2 = runpars->PixelCenter(pixid);
xmc.push_back(dummy2.x()/1000);//m
ymc.push_back(dummy2.y()/1000);//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/", Idev);
gDirectory->cd(pathname.c_str());
string pathnames = Form("SingleHoughMacrocell%d", idmacrocell);
TDirectory *path = new TDirectory(pathnames.c_str(),pathnames.c_str());
path->cd();
Int_t possiblesignal=HoughTransform(xmc,ymc,tmc,cmc);
if( possiblesignal ) {
Msg(EsafMsg::Debug) <<"Macrocell "<<idmacrocell<<" probably HAS signal!!!"<< MsgDispatch;
idmacrocellselected.push_back(idmacrocell);
fHittedFinalMacrocell.insert(fHittedFinalMacrocell.end(),fHittedmc.begin(),fHittedmc.end());
}else{
Msg(EsafMsg::Debug) <<"Macrocell "<<idmacrocell<<" probably DOESN'T HAVE a signal,";
if(idmacrocell==fMainMacroCellId){
Msg(EsafMsg::Debug) <<" but is added because is the main!"<< MsgDispatch;
fHittedFinalMacrocell.insert(fHittedFinalMacrocell.end(),fHittedmc.begin(),fHittedmc.end());
idmacrocellselected.push_back(idmacrocell);
}else{
Msg(EsafMsg::Debug) <<" discarded!"<< MsgDispatch;
}
}
xmc.clear();
ymc.clear();
tmc.clear();
cmc.clear();
}
if (idmacrocellselected.size()<1 || fHittedFinalMacrocell.size()<10) return 0;
Msg(EsafMsg::Info)<<"Re-process "<<idmacrocellselected.size()<<" selected macrocells by ProcessHoughSingleMacrocell()"<< MsgDispatch;
Int_t ret=ProcessHoughAll();
return ret;
}
////////////////////////////////////////////////////////
//Doing histograms of velocities for each macrocell hitted and trying to select which are hitted by real signal
Int_t HoughModule::ProcessVelocity(){
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> vuoto;
pair< Int_t,vector<Int_t> > p(i,vuoto);
mapIdfHittedSingleCell.insert(p);
}
}
Msg(EsafMsg::Debug) << "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=50;
Msg(EsafMsg::Debug)<<"selected pixel with counts > "<<numcountsmin<<MsgDispatch;
for( it = mapIdfHittedSingleCell.begin(); it != mapIdfHittedSingleCell.end(); it++ ){
Int_t entriesTH2(0);
Msg(EsafMsg::Debug)<<"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);
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 puntidaaccoppiare=fHittedmc.size();
for(Int_t j=i+1;j<puntidaaccoppiare;j++){
Int_t pixidaccop=fEv->GetRecoPixelData(fHittedmc[j])->GetPixelId();
TVector3 dummyaccop = runpars->PixelCenter(pixidaccop);
Int_t tmcgtuaccop=fEv->GetRecoPixelData(fHittedmc[j])->GetGtu();
if(fabs(dummy2.x()-dummyaccop.x())>0 && fabs(dummy2.y()-dummyaccop.y())>0 && tmcgtuaccop-tmcgtu>0 && tmcgtuaccop-tmcgtu<20 && fabs(dummy2.x()-dummyaccop.x())<100 && fabs(dummy2.y()-dummyaccop.y())<100 ){
Int_t countprodotto=(Int_t)sqrt(1.e0*counts*fEv->GetRecoPixelData(fHittedmc[j])->GetCounts());
Double_t m=(dummy2.y()-dummyaccop.y())/(dummy2.x()-dummyaccop.x());
Double_t v=sqrt((dummy2.y()-dummyaccop.y())*(dummy2.y()-dummyaccop.y())+(dummy2.x()-dummyaccop.x())*(dummy2.x()-dummyaccop.x()))/((tmcgtuaccop-tmcgtu)*2.5);//millimetri/microsecondo
Double_t vx=(dummy2.x()-dummyaccop.x())/((tmcgtu-tmcgtuaccop)*2.5);
Double_t vy=(dummy2.y()-dummyaccop.y())/((tmcgtu-tmcgtuaccop)*2.5);
entriesTH2++;
s->Fill(v,m,countprodotto);
sx->Fill(vx,countprodotto);
sy->Fill(vy,countprodotto);
}
}
}
Stat_t contx(0),conty(0);
for(Int_t bin=1;bin<binv;bin++){
contx=contx+sx->GetBinContent(bin);
conty=conty+sy->GetBinContent(bin);
}
Float_t contmaxx=sx->GetMaximum();
Float_t contmaxy=sy->GetMaximum();
Float_t contmediox=(Float_t)((Float_t)contx/binv);
Float_t contmedioy=(Float_t)((Float_t)conty/binv);
Float_t ms=10;
/*cout<<"contmediox="<<contmediox<<"tcontmaxx="<<contmaxx<<"tcontmaxx-contmediox="<<contmaxx-contmediox<<"tsig="<<ms*sqrt(contmediox)<<endl;
cout<<"contmedioy="<<contmedioy<<"tcontmaxy="<<contmaxy<<"tcontmaxy-contmedioy="<<contmaxy-contmedioy<<"tsig="<<ms*sqrt(contmedioy)<<endl;*/
string a="a";
string b="b";
string pathname = Form("/HoughModule%d/", Idev);
gDirectory->cd(pathname.c_str());
string pathnames;
if( ((contmaxx-contmediox) > ms*sqrt(contmediox) || (contmaxy-contmedioy) > ms*sqrt(contmedioy)) && sx->GetEntries()>10 ){
fHittedFinalMacrocell.insert(fHittedFinalMacrocell.end(),fHittedmc.begin(),fHittedmc.end());
idmacrocellselected.push_back(idmacrocell);
Msg(EsafMsg::Debug) <<"Macrocell "<<idmacrocell<<" probably HAS signal!!!"<< MsgDispatch;
if (idmacrocell==fMainMacroCellId){
pathnames = Form("VelocityMacrocellmain%d", idmacrocell);
}else{
pathnames = Form("VelocityMacrocell%d", idmacrocell);
}
}else{
Msg(EsafMsg::Debug) <<"Macrocell "<<idmacrocell<<" probably DOESN'T HAVE a signal,";
if(idmacrocell==fMainMacroCellId){
pathnames = Form("VelocityMacrocellmain%d", idmacrocell);
Msg(EsafMsg::Debug) <<" but is added because is the main!"<< MsgDispatch;
fHittedFinalMacrocell.insert(fHittedFinalMacrocell.end(),fHittedmc.begin(),fHittedmc.end());
idmacrocellselected.push_back(idmacrocell);
}else{
pathnames = Form("VelocityMacrocelldiscarded%d", idmacrocell);
Msg(EsafMsg::Debug) <<" 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,1);
bc->cd(1);hsxt->Draw();
bc->cd(2);hsyt->Draw();
bc->Write();
delete bc; delete hsxt; delete hsyt;
}
if (idmacrocellselected.size()<1 || fHittedFinalMacrocell.size()<10) return 0;
Msg(EsafMsg::Info)<<"Re-process "<<idmacrocellselected.size()<<" selected macrocells by ProcessVelocity()"<< MsgDispatch;
Int_t ret=ProcessHoughAll();
return ret;
}
////////////////////////////////////////////////////////
//Doing Hough Transform to all macrocells selected by previous methods or by euso trigger
Int_t HoughModule::ProcessHoughAll(){
Msg(EsafMsg::Info) << " ...ProcessHoughAll()" << MsgDispatch;
string pathname = Form("/HoughModule%d/", Idev);
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::Debug) << "Let's try beginning with MinCountsHoughAll = " << fNumHitsMinimum << MsgDispatch;
}
Int_t regr=0;
Float_t gtumin=1000;
Float_t gtumax=0;
Float_t deltagtu=0;
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);
}
}
if(fHittedFinalMacrocell.size()>0){
fNumPointsall=0;
Msg(EsafMsg::Debug)<<"processing fHittedFinalMacrocell.size()="<<fHittedFinalMacrocell.size()<<MsgDispatch;
fHittedall.clear();
for( Int_t i=0; i<(Int_t)(fHittedFinalMacrocell.size()); i++ ){
if ( fEv->GetRecoPixelData(fHittedFinalMacrocell[i])->GetCounts() >= fNumHitsMinimum ) {
fHittedall.push_back(fHittedFinalMacrocell[i]);
}
}
}
fNumPointsall=fHittedall.size();
tmin=gtumin*2.5*0.01;
tmax=gtumax*2.5*0.01;
deltagtu=gtumax-gtumin;
Msg(EsafMsg::Debug) << "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,tmin,tmax,bin,-1,1);
TH2F *single2yt = new TH2F("single2yt","single2yt;t;y(m)",bin,tmin,tmax,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(dummy2.x()/1000);
y.push_back(dummy2.y()/1000);
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]);
}
Msg(EsafMsg::Debug) << "Num. of pixels with at least " << fNumHitsMinimum << " hits = " << fNumPointsall << MsgDispatch;
TCanvas *cHTsingle2hough=new TCanvas("cHTsingle2hough","cHTsingle2hough",1);
cHTsingle2hough->Divide(3,1);
cHTsingle2hough->cd(1); single2->SetMarkerStyle(6); single2->Draw();
cHTsingle2hough->cd(2); single2xt->SetMarkerStyle(6); single2xt->Draw();
cHTsingle2hough->cd(3); single2yt->SetMarkerStyle(6); single2yt->Draw();
cHTsingle2hough->Write();
delete cHTsingle2hough;
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,xmin,xmax,bin,ymin,ymax);
TH2F *single2xtsel = new TH2F("single2xtsel","single2xtsel;t;x(m)",bin,tmin,tmax,bin,xmin,xmax);
TH2F *single2ytsel = new TH2F("single2ytsel","single2ytsel;t;y(m)",bin,tmin,tmax,bin,ymin,ymax);
for( Int_t i=0; i<fTotalNumPoints; i++ ) {
Int_t pixid=fEv->GetRecoPixelData(i)->GetPixelId();
TVector3 dummy2 = runpars->PixelCenter(pixid);
Double_t x=dummy2.x()/1000;
Double_t y=dummy2.y()/1000;
Double_t t=(Double_t)fEv->GetRecoPixelData(i)->GetGtu()*2.5*0.01;
Int_t c=fEv->GetRecoPixelData(i)->GetCounts();
if ( x>=xmin && x<=xmax && y>=ymin && y<=ymax && 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::Debug) << "first selected zone: pixels with at least " << fNumHitsMinimum-regr << " hits = " << fHittedselhough.size() << MsgDispatch;
if(fHittedselhough.size()<10){
Msg(EsafMsg::Debug) << "too few! exit from this event" <<MsgDispatch;
return kFALSE;
}
TCanvas *cHTsingle2selhough=new TCanvas("cHTsingle2selhough","cHTsingle2selhough",1);
cHTsingle2selhough->Divide(3,1);
cHTsingle2selhough->cd(1); single2sel->SetMarkerStyle(6); single2sel->Draw();
cHTsingle2selhough->cd(2); single2xtsel->SetMarkerStyle(6); single2xtsel->Draw();
cHTsingle2selhough->cd(3); single2ytsel->SetMarkerStyle(6); single2ytsel->Draw();
cHTsingle2selhough->Write();
delete cHTsingle2selhough;
delete single2xtsel;
delete single2ytsel;
delete single2sel;
HoughTransform(xsel,ysel,tsel,csel);
xsel.clear();
ysel.clear();
tsel.clear();
csel.clear();
if(fNumHitsMinimum-regr>3) 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,xmin,xmax,bin,ymin,ymax);
TH2F *single2xtselsel = new TH2F("single2xtselsel","single2xtselsel;t;x(m)",bin,tmin,tmax,bin,xmin,xmax);
TH2F *single2ytselsel = new TH2F("single2ytselsel","single2ytselsel;t;y(m)",bin,tmin,tmax,bin,ymin,ymax);
for( Int_t i=0; i<fTotalNumPoints; i++ ) {
Int_t pixid=fEv->GetRecoPixelData(i)->GetPixelId();
TVector3 dummy2 = runpars->PixelCenter(pixid);
Double_t x=dummy2.x()/1000;
Double_t y=dummy2.y()/1000;
Double_t t=(Double_t)fEv->GetRecoPixelData(i)->GetGtu()*2.5*0.01;
Int_t c=fEv->GetRecoPixelData(i)->GetCounts();
if ( x>=xmin && x<=xmax && y>=ymin && y<=ymax && 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);
fHittedselselhough.push_back(i);
}
}
Msg(EsafMsg::Debug) << "second selected zone: pixels with at least " << fNumHitsMinimum-regr << " hits = " << fHittedselselhough.size()<< MsgDispatch;
if(fHittedselselhough.size()<10){
Msg(EsafMsg::Debug) << "too few! exit from this event" <<MsgDispatch;
return kFALSE;
}
TCanvas *cHTsingle2selselhough=new TCanvas("cHTsingle2selselhough","cHTsingle2selselhough",1);
cHTsingle2selselhough->Divide(3,1);
cHTsingle2selselhough->cd(1); single2selsel->SetMarkerStyle(6); single2selsel->Draw();
cHTsingle2selselhough->cd(2); single2xtselsel->SetMarkerStyle(6); single2xtselsel->Draw();
cHTsingle2selselhough->cd(3); single2ytselsel->SetMarkerStyle(6); single2ytselsel->Draw();
cHTsingle2selselhough->Write();
delete cHTsingle2selselhough;
delete single2xtselsel;
delete single2ytselsel;
delete single2selsel;
HoughTransform(xselsel,yselsel,tselsel,cselsel);
xselsel.clear();
yselsel.clear();
tselsel.clear();
cselsel.clear();
if( fNumHitsMinimum-regr > 3 ) regr++;
for( Int_t i=0; i<fTotalNumPoints; i++ ) {
Int_t pixid=fEv->GetRecoPixelData(i)->GetPixelId();
TVector3 dummy2 = runpars->PixelCenter(pixid);
Double_t x=dummy2.x()/1000;
Double_t y=dummy2.y()/1000;
Int_t c=fEv->GetRecoPixelData(i)->GetCounts();
if ( x>=xmin && x<=xmax && y>=ymin && y<=ymax && c>=(fNumHitsMinimum-regr)) {
fHittedFinal.push_back(i);
}
}
//fHittedFinal=fHittedselhough;
Msg(EsafMsg::Debug) << "last selected zone: pixels with at least " << fNumHitsMinimum-regr << " hits = " << fHittedFinal.size() <<MsgDispatch;
fFinalNumHitsMinimum=fNumHitsMinimum-regr;
return 1;
}
/////////////////////////////////////////////////////////
//
Int_t HoughModule::HoughTransform(vector<Double_t> x,vector<Double_t> y,vector<Double_t> t,vector<Int_t> c){
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);
//Double_t vx = xtFit->GetSlope();//Double_t wx = xtFit->GetOffset();
xmin=xtFit->GetMiny();
xmax=xtFit->GetMaxy();
Double_t tmin1=xtFit->GetMinx();
Double_t tmax1=xtFit->GetMaxx();
absdevx = xtFit->GetAbsDev();
Int_t possiblesignalxt=xtFit->fthereissignal;
cmaxtocmediox=xtFit->GetMaxToMedium();
delete xtFit;
HoughFit *ytFit = new HoughFit(t, y, c, optionselection,et,ey);
//Double_t vy = ytFit->GetSlope();//Double_t wy = ytFit->GetOffset();
ymin=ytFit->GetMiny();
ymax=ytFit->GetMaxy();
Double_t tmin2=ytFit->GetMinx();
Double_t tmax2=ytFit->GetMaxx();
absdevy = ytFit->GetAbsDev();
Int_t possiblesignalyt=ytFit->fthereissignal;
cmaxtocmedioy=ytFit->GetMaxToMedium();
delete ytFit;
Msg(EsafMsg::Debug) <<"possiblesignalxt="<<possiblesignalxt<<"tpossiblesignalyt="<<possiblesignalyt<<MsgDispatch;
Msg(EsafMsg::Debug) <<"cmaxtocmediox="<<cmaxtocmediox<<"tcmaxtocmedioy="<<cmaxtocmedioy<<MsgDispatch;
Int_t possible(0);
if( (possiblesignalxt==1 || possiblesignalyt==1) && (cmaxtocmediox>1.5*cmaxtocmediominx || cmaxtocmedioy>1.5*cmaxtocmediominy)){
possible=1;
}else{possible=0;}
tmin = ( tmin1<tmin2 ? tmin1 : tmin2);
tmax = ( tmax1<tmax2 ? tmax2 : tmax1);
x.clear();
y.clear();
t.clear();
c.clear();
return possible;
}
//////////////////////////////////////////////////////////////////////////////////////////////////
// post-process
Bool_t HoughModule::PostProcess() {
if ( fTotalNumPoints == 0 ) return kTRUE;
fEv = (RecoEvent*)0;
return kTRUE;
}
Bool_t HoughModule::SaveRootData(RecoRootEvent *fRecoRootEvent) {
return kTRUE;
}
// done
Bool_t HoughModule::Done() {
fHittedall.clear();
fHittedFinal.clear();
fHittedselhough.clear();
fHittedselselhough.clear();
idmacrocellselected.clear();
fClusters.clear();
Msg(EsafMsg::Info) << "Completed" << MsgDispatch;
return kTRUE;
}
// memory clean
void HoughModule::UserMemoryClean() {
Int_t fHittedFinalsize = fHittedFinal.size();
if(fHittedFinalsize<10 || fTotalNumPoints == 0 || fHittedFinalsize>10000) 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");
}
// close and write the current cluster if is significant
void HoughModule::WriteSingleCluster() {
// 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;
}
// get a specified cluster
EusoCluster* HoughModule::GetCluster( Int_t i ) {
if ( i >= 0 && i < (Int_t)fClusters.size() )
return fClusters[i];
else
return NULL;
}