// ESAF : Euso Simulation and Analysis Framework
// $Id: HoughFit.cc,v 1.3 2005/03/30 10:35:07 etaddei Exp $
// Elena Taddei created Mar, 29 2005
//////////////////////////////////////////////////////////////////////////////
// //
// HoughFit //
// //
// //
//////////////////////////////////////////////////////////////////////////////
#include <iostream>
#include "MedianFit.hh"
#include "HoughFit.hh"
#include "TMath.h"
#include <TH2.h>
#include <TCanvas.h>
#include <TGraph.h>
#include "Config.hh"
#include <TMinuit.h>
ClassImp(HoughFit)
//_________________________________________________________________________________________________________________________
HoughFit::HoughFit( vector<Double_t> x, vector<Double_t> y, vector<Int_t> c, Int_t fDoSelection, Float_t ex, Float_t ey) {
//
//Constructor
//
ConfigFileParser *pConfigHoughFit = Config::Get()->GetCF("Reco","HoughModule");
string fKindHough = pConfigHoughFit->GetStr("HoughModule.KindHough");
Clear();
epsx=ex;
epsy=ey;
fPointsX = x;
fPointsY = y;
fCounts = c;
fNumPoints=fPointsX.size();
if(fNumPoints<5) return;
if( fKindHough=="sinusoidal" ){
DoSinusoidal(1);
cmaxtocmedio=cmaxtocmediotmp;
Float_t minfrac=pConfigHoughFit->GetNum("HoughModule.fMinFractionSinusoidal");
if( cmaxtocmedio > minfrac && fDoSelection==1 && fNumPointsselected>5 ){
//cout<<"There could be a signal, try to do selection.."<<endl;
fPointsX = xsel;
fPointsY = ysel;
fCounts = csel;
fNumPoints=fPointsX.size();
Clear();
DoSinusoidal(2);
if( cmaxtocmediotmp>cmaxtocmedio || cmaxtocmediotmp>minfrac ) fthereissignal=1;
}
}
if( fKindHough=="linear" ){
DoLinear(1);
cmaxtocmedio=cmaxtocmediotmp;
Float_t minfrac=pConfigHoughFit->GetNum("HoughModule.fMinFractionLinear");
if( cmaxtocmedio > minfrac && fDoSelection==1 && fNumPointsselected>5 ){
//cout<<"There could be a signal, try to do selection.."<<endl;
fPointsX = xsel;
fPointsY = ysel;
fCounts = csel;
fNumPoints=fPointsX.size();
Clear();
DoLinear(2);
if( cmaxtocmediotmp>cmaxtocmedio || cmaxtocmediotmp>minfrac ) fthereissignal=1;
}
}
if( fKindHough=="numerical" ){
cout<<"DoNumerical() not implemented yet.."<<endl;//DoNumerical();
return;
}
}
//____________________________________
HoughFit::~HoughFit() {
//
//Destructor
//
fPointsX.clear();
fPointsY.clear();
fCounts.clear();
fNumPointsselected=0;
fNumPoints=0;
}
//____________________________________________________________
void HoughFit::Clear() {
fthereissignal=0;
xsel.clear();
ysel.clear();
csel.clear();
fSlopetopass = 0;
fOffsettopass = 0;
fSlope = 0;
fOffset = 0;
fwidth = 0;
absdev = 0;
}
//______________________________________________________________________________________
void HoughFit::DoSinusoidal(Int_t firsttimedo) {
//
//recognition of linear tracks with sinusoidal Hough Transform in the space of parameters (ro,theta)
// P_i: (x_i,y_i) --> S_i: ro=x_i*cos(theta)+y_i*sin(theta)
//
fminx=100;
fminy=100;
fmaxx=-100;
fmaxy=-100;
for( Int_t i=0; i<fNumPoints; i++ ) {
Double_t x=fPointsX[i];
Double_t y=fPointsY[i];
if(x>fmaxx) fmaxx=x;
if(y>fmaxy) fmaxy=y;
if(x<fminx) fminx=x;
if(y<fminy) fminy=y;
}
Double_t x0=(fmaxx+fminx)/2;
Double_t y0=(fmaxy+fminy)/2;
Float_t tmpfminx_pre=100;
Float_t tmpfminy_pre=100;
Float_t tmpfmaxx_pre=-100;
Float_t tmpfmaxy_pre=-100;
vector<Double_t> fPointsXprimo_pre;
vector<Double_t> fPointsYprimo_pre;
for( Int_t i=0; i<fNumPoints; i++ ) {
Double_t x=fPointsX[i]-x0;
Double_t y=fPointsY[i]-y0;
fPointsXprimo_pre.push_back(x);
fPointsYprimo_pre.push_back(y);
if(x>tmpfmaxx_pre) tmpfmaxx_pre=x;
if(y>tmpfmaxy_pre) tmpfmaxy_pre=y;
if(x<tmpfminx_pre) tmpfminx_pre=x;
if(y<tmpfminy_pre) tmpfminy_pre=y;
}
Float_t tmpfminx=100;
Float_t tmpfminy=100;
Float_t tmpfmaxx=-100;
Float_t tmpfmaxy=-100;
Float_t scale=tmpfmaxx_pre/tmpfmaxy_pre;
vector<Double_t> fPointsXprimo;
vector<Double_t> fPointsYprimo;
for( Int_t i=0; i<fNumPoints; i++ ) {
Double_t x=fPointsXprimo_pre[i];
Double_t y=fPointsYprimo_pre[i]*scale;
fPointsXprimo.push_back(x);
fPointsYprimo.push_back(y);
if(x>tmpfmaxx) tmpfmaxx=x;
if(y>tmpfmaxy) tmpfmaxy=y;
if(x<tmpfminx) tmpfminx=x;
if(y<tmpfminy) tmpfminy=y;
}
Int_t binmax(0),nx(0),binymax(0),binxmax(0);
Double_t thetamaxsingle,romaxsingle;
Float_t romax=TMath::Sqrt(tmpfmaxx*tmpfmaxx+tmpfmaxy*tmpfmaxy);
Float_t romin=-romax;
Float_t epsilon=TMath::Sqrt(epsx*epsx+epsy*epsy);
Float_t deltaro=2*epsilon;
Float_t thetasinglemin=0;
Float_t thetasinglemax=TMath::Pi();
Int_t binx=(Int_t)((tmpfmaxx-tmpfminx)/0.001);
Int_t biny=(Int_t)((tmpfmaxy-tmpfminy)/0.001);
Int_t numbinro=(Int_t)((romax-romin)/deltaro);
Float_t deltatheta=2*deltaro/(romax-romin);
Int_t numbintheta=(Int_t)((thetasinglemax-thetasinglemin)/deltatheta);
Int_t numsteptheta=numbintheta;
Float_t steptheta=(thetasinglemax-thetasinglemin)/numsteptheta;
TH2F *single = new TH2F("single","single;theta;ro",numbintheta,thetasinglemin,thetasinglemax,numbinro,romin,romax);
TH2F *single2 = new TH2F("single2","single2",binx,tmpfminx,tmpfmaxx,biny,tmpfminy,tmpfmaxy);
for( Int_t i=0; i<fNumPoints; i++ ) {
Double_t x=fPointsXprimo[i];
Double_t y=fPointsYprimo[i];
Int_t z=fCounts[i];
single2->Fill(x,y,z);
for(Int_t j=0;j<numsteptheta;j++){
Double_t thetasingle=j*steptheta;
Double_t rosingle=x*cos(thetasingle)+y*sin(thetasingle);
single->Fill(thetasingle,rosingle,z);
}
}
binmax = single->GetMaximumBin();
nx = (single->GetXaxis())->GetNbins()+2;
Int_t plus1=0;
Int_t plus2=0;
while( (binmax % nx < 2 || binmax % nx > nx-2) && plus1<10 && plus2<10){
if(binmax % nx < 2) plus1++;
if(binmax % nx > nx-2) plus2++;
single->GetXaxis()->SetRange(plus1,numbintheta-plus2);
binmax = single->GetMaximumBin();
nx = (single->GetXaxis())->GetNbins()+2;
}
binxmax = binmax % nx;
binymax = binmax / nx;
thetamaxsingle = (single->GetXaxis())->GetBinCenter(binxmax);
romaxsingle = (single->GetYaxis())->GetBinCenter(binymax);
Float_t contmax=single->GetBinContent(binxmax,binymax);
TH2F *singlelt = new TH2F("singlelt","singlelt;theta;ro",numbintheta,thetasinglemin,thetasinglemax,numbinro,romin,romax);
TH2F *singlelt1 = new TH2F("singlelt1","singlelt1;theta;ro",numbintheta,thetasinglemin,thetasinglemax,numbinro,romin,romax);
TH2F *singlelt2 = new TH2F("singlelt2","singlelt2;theta;ro",numbintheta,thetasinglemin,thetasinglemax,numbinro,romin,romax);
Float_t consideredbin=0;
Float_t conttotale=0;
for(Int_t ith=1;ith<numbintheta;ith++){
for(Int_t iro=1;iro<numbinro;iro++){
Int_t cont=(Int_t)single->GetBinContent(ith,iro);
if(cont>0.25*contmax )singlelt->Fill((single->GetXaxis())->GetBinCenter(ith),(single->GetYaxis())->GetBinCenter(iro),cont);
if(cont>0.5*contmax )singlelt1->Fill((single->GetXaxis())->GetBinCenter(ith),(single->GetYaxis())->GetBinCenter(iro),cont);
if(cont>0.75*contmax )singlelt2->Fill((single->GetXaxis())->GetBinCenter(ith),(single->GetYaxis())->GetBinCenter(iro),cont);
if(cont<0.9*contmax ){
conttotale=conttotale+cont;
consideredbin++;
}
}
}
Float_t contmedio=conttotale/(consideredbin);
//Float_t flutt=TMath::Sqrt(contmedio);
cmaxtocmediotmp=(Float_t)(contmax/contmedio);
TCanvas *cHTlt=new TCanvas("cHTlt","cHTlt",1);cHTlt->Divide(3,1);
cHTlt->cd(1); singlelt->Draw();
cHTlt->cd(2); singlelt1->Draw();
cHTlt->cd(3); singlelt2->Draw();
cHTlt->Write();
delete cHTlt;
delete singlelt;
delete singlelt1;
delete singlelt2;
TH1D *singleproj=new TH1D("singleproj","singleproj",numbinro,romin,romax);
singleproj=single->ProjectionY("singleproj",binxmax-5,binxmax+5);
fwidth = singleproj->GetRMS();
fSlope = -1/tan(thetamaxsingle);
fOffset = romaxsingle/sin(thetamaxsingle);
fSlopetopass=fSlope/scale;
fOffsettopass=fOffset/scale;
Double_t xg[fNumPoints],yg[fNumPoints];
absdev=0;
Double_t markercenterx[1]={thetamaxsingle};
Double_t markercentery[1]={romaxsingle};
TGraph *gmarkercenter = new TGraph(1,markercenterx,markercentery);
gmarkercenter->SetMarkerColor(2); gmarkercenter->SetMarkerStyle(29); gmarkercenter->SetMarkerSize(1.6);
for( Int_t i=0; i<fNumPoints; i++ ) {
xg[i]=fPointsXprimo[i];
yg[i]=fSlope*fPointsXprimo[i]+fOffset;
Float_t singleabsdev=fabs(yg[i]-fPointsYprimo[i]);
absdev=absdev+singleabsdev;
}
absdev=absdev/fNumPoints;
TGraph *g = new TGraph(fNumPoints,xg,yg);g->SetLineColor(8);
string namecHT;
if(firsttimedo==2){
namecHT="cHTsecond";
}else{
namecHT="cHTfirst";
}
TCanvas *cHT=new TCanvas(namecHT.c_str(),namecHT.c_str(),1);cHT->Divide(3,1);
cHT->cd(1);
single2->SetMarkerStyle(6);
single2->Draw();
g->Draw("PLsame");
cHT->cd(2); single->Draw();gmarkercenter->Draw("Psame");
cHT->cd(3); single->Draw("LEGO");
cHT->Write();
delete cHT;delete single2;delete single;delete gmarkercenter;delete g;
delete singleproj;
fNumPointsselected=0;
if(firsttimedo==1 && fwidth>1e-4){
for( Int_t i=0; i<fNumPoints; i++ ) {
Double_t x=fPointsXprimo[i];
Double_t y=fPointsYprimo[i];
Int_t z=fCounts[i];
Double_t rosingleprova = x*cos(thetamaxsingle)+y*sin(thetamaxsingle);
if(fabs(rosingleprova-romaxsingle)<fwidth*0.3){
xsel.push_back(x+x0);
ysel.push_back((y/scale)+y0);
csel.push_back(z);
fNumPointsselected++;
}
}
//cout<<"nfNumPoints="<<fNumPoints<<"tfNumPointsselected="<<fNumPointsselected<<"tfwidth="<<fwidth<<endl;
}
}
//_______________________________________________________________________________
void HoughFit::DoLinear(Int_t firsttimedo) {
//
//recognition of linear tracks with linear Hough Transform in the space of parameters (q,m)
// P_i: (x_i,y_i) --> L_i: q=-m*x_i+y_i
//
fminx=100;
fminy=100;
fmaxx=-100;
fmaxy=-100;
for( Int_t i=0; i<fNumPoints; i++ ) {
Double_t x=fPointsX[i];
Double_t y=fPointsY[i];
if(x>fmaxx) fmaxx=x;
if(y>fmaxy) fmaxy=y;
if(x<fminx) fminx=x;
if(y<fminy) fminy=y;
}
Double_t x0=(fmaxx+fminx)/2;
Double_t y0=(fmaxy+fminy)/2;
Float_t tmpfminx_pre=100;
Float_t tmpfminy_pre=100;
Float_t tmpfmaxx_pre=-100;
Float_t tmpfmaxy_pre=-100;
vector<Double_t> fPointsXprimo_pre;
vector<Double_t> fPointsYprimo_pre;
for( Int_t i=0; i<fNumPoints; i++ ) {
Double_t x=fPointsX[i]-x0;
Double_t y=fPointsY[i]-y0;
fPointsXprimo_pre.push_back(x);
fPointsYprimo_pre.push_back(y);
if(x>tmpfmaxx_pre) tmpfmaxx_pre=x;
if(y>tmpfmaxy_pre) tmpfmaxy_pre=y;
if(x<tmpfminx_pre) tmpfminx_pre=x;
if(y<tmpfminy_pre) tmpfminy_pre=y;
}
Float_t tmpfminx=100;
Float_t tmpfminy=100;
Float_t tmpfmaxx=-100;
Float_t tmpfmaxy=-100;
Float_t scale=tmpfmaxx_pre/tmpfmaxy_pre;
vector<Double_t> fPointsXprimo;
vector<Double_t> fPointsYprimo;
for( Int_t i=0; i<fNumPoints; i++ ) {
Double_t x=fPointsXprimo_pre[i];
Double_t y=fPointsYprimo_pre[i]*scale;
fPointsXprimo.push_back(x);
fPointsYprimo.push_back(y);
if(x>tmpfmaxx) tmpfmaxx=x;
if(y>tmpfmaxy) tmpfmaxy=y;
if(x<tmpfminx) tmpfminx=x;
if(y<tmpfminy) tmpfminy=y;
}
Int_t binmax(0),nx(0),binymax(0),binxmax(0);
Double_t mhoughmax=tmpfmaxx;
Double_t qhoughmax=tmpfmaxy;
Double_t mmax=2;//to establish phisically from EAS velocity, still to improve
Double_t mmin=-mmax;
Double_t qmax=mmax*mhoughmax+qhoughmax;
Double_t qmin=-qmax;
Float_t deltaq=TMath::Sqrt(mmax*mmax*epsx*epsx+epsy*epsy);
Int_t numbinq=(Int_t)((qmax-qmin)/deltaq);
Float_t deltam=deltaq/mhoughmax;
Int_t numbinm=(Int_t)((mmax-mmin)/deltam);
Int_t numstepm=numbinm;
Float_t stepm=deltam;
Int_t binx=(Int_t)((tmpfmaxx-tmpfminx)/0.001);//not important binning, is just to draw
Int_t biny=(Int_t)((tmpfmaxy-tmpfminy)/0.001);//not important binning, is just to draw
TH2F *single = new TH2F("single","single;m;q",numbinm,mmin,mmax,numbinq,qmin,qmax);
TH2F *single2 = new TH2F("single2","single2",binx,tmpfminx,tmpfmaxx,biny,tmpfminy,tmpfmaxy);
for( Int_t i=0; i<fNumPoints; i++ ) {
Double_t x=fPointsXprimo[i];
Double_t y=fPointsYprimo[i];
Int_t z=fCounts[i];
single2->Fill(x,y,z);
for(Int_t j=0;j<numstepm;j++){
Double_t msingle=mmin+j*stepm;
Double_t qsingle=-msingle*x+y;
single->Fill(msingle,qsingle,z);
}
}
binmax = single->GetMaximumBin();
nx = (single->GetXaxis())->GetNbins()+2;
Int_t plus1=0;
Int_t plus2=0;
while( (binmax % nx < 2 || binmax % nx > nx-2) && plus1<10 && plus2<10){
if(binmax % nx < 2) plus1++;
if(binmax % nx > nx-2) plus2++;
single->GetXaxis()->SetRange(plus1,numbinm-plus2);
binmax = single->GetMaximumBin();
nx = (single->GetXaxis())->GetNbins()+2;
}
binxmax = binmax % nx;
binymax = binmax / nx;
Float_t mmaxsingle = (single->GetXaxis())->GetBinCenter(binxmax);
Float_t qmaxsingle = (single->GetYaxis())->GetBinCenter(binymax);
Float_t contmax=single->GetBinContent(binxmax,binymax);
fSlope=mmaxsingle;
fOffset=qmaxsingle;
fSlopetopass=fSlope/scale;
fOffsettopass=fOffset/scale;
Double_t markercenterx[1]={mmaxsingle};
Double_t markercentery[1]={qmaxsingle};
TGraph *gmarkercenter = new TGraph(1,markercenterx,markercentery);
gmarkercenter->SetMarkerColor(2); gmarkercenter->SetMarkerStyle(29); gmarkercenter->SetMarkerSize(1.6);
Double_t xg[fNumPoints],yg[fNumPoints];
absdev=0;
for( Int_t i=0; i<fNumPoints; i++ ) {
xg[i]=fPointsXprimo[i];
yg[i]=fSlope*fPointsXprimo[i]+fOffset;
Float_t singleabsdev=fabs(yg[i]-fPointsYprimo[i]);
absdev+=singleabsdev;
}
absdev=absdev/fNumPoints;
TGraph *g = new TGraph(fNumPoints,xg,yg);g->SetLineColor(8);
TH2F *singlelt = new TH2F("singlelt","singlelt;m;q",numbinm,mmin,mmax,numbinq,qmin,qmax);
TH2F *singlelt1 = new TH2F("singlelt1","singlelt1;m;q",numbinm,mmin,mmax,numbinq,qmin,qmax);
TH2F *singlelt2 = new TH2F("singlelt2","singlelt2;m;q",numbinm,mmin,mmax,numbinq,qmin,qmax);
Float_t consideredbin=0;
Float_t conttotale=0;
for(Int_t im=1;im<numbinm;im++){
for(Int_t iq=1;iq<numbinq;iq++){
Int_t cont=(Int_t)(single->GetBinContent(im,iq));
if(cont>0.25*contmax )singlelt->Fill((single->GetXaxis())->GetBinCenter(im),(single->GetYaxis())->GetBinCenter(iq),cont);
if(cont>0.5*contmax )singlelt1->Fill((single->GetXaxis())->GetBinCenter(im),(single->GetYaxis())->GetBinCenter(iq),cont);
if(cont>0.75*contmax )singlelt2->Fill((single->GetXaxis())->GetBinCenter(im),(single->GetYaxis())->GetBinCenter(iq),cont);
if(cont<0.9*contmax ){
conttotale=conttotale+cont;
consideredbin++;
}
}
}
Float_t contmedio=conttotale/(consideredbin);
//Float_t flutt=TMath::Sqrt(contmedio);
cmaxtocmediotmp=(Float_t)(contmax/contmedio);
TH1D *singleproj=new TH1D("singleproj","singleproj",numbinq,qmin,qmax);
singleproj=single->ProjectionY("singleproj",binxmax-5,binxmax+5);
fwidth = singleproj->GetRMS();
singleproj->Draw();
TCanvas *cHTlt=new TCanvas("cHTlt","cHTlt",1);cHTlt->Divide(2,2);
cHTlt->cd(1); singlelt->Draw();
cHTlt->cd(2); singlelt1->Draw();
cHTlt->cd(3); singlelt2->Draw();
cHTlt->cd(4); singleproj->Draw();
cHTlt->Write();
delete cHTlt; delete singlelt; delete singlelt1; delete singlelt2; delete singleproj;
string namecHT;
if(firsttimedo==2){
namecHT="cHTlinearsecond";
}else{
namecHT="cHTlinearfirst";
}
TCanvas *cHT=new TCanvas(namecHT.c_str(),namecHT.c_str(),1);cHT->Divide(3,1);
cHT->cd(1);
single2->SetMarkerStyle(6);
single2->Draw();g->Draw("PLsame");
cHT->cd(2); single->Draw();gmarkercenter->Draw("Psame");
cHT->cd(3); single->Draw("LEGO");
cHT->Write();
delete cHT;delete single2;delete single;delete gmarkercenter;delete g;
fNumPointsselected=0;
if(firsttimedo==1 && fwidth>1e-4){
for( Int_t i=0; i<fNumPoints; i++ ) {
Double_t x=fPointsXprimo[i];
Double_t y=fPointsYprimo[i];
Int_t z=fCounts[i];
Double_t qsingleprova = -mmaxsingle*x+y;
if(fabs(qsingleprova-qmaxsingle)<fwidth*0.3){
xsel.push_back(x+x0);
ysel.push_back((y/scale)+y0);
csel.push_back(z);
fNumPointsselected++;
}
}
//cout<<"nfNumPoints="<<fNumPoints<<"tfNumPointsselected="<<fNumPointsselected<<"tfwidth="<<fwidth<<endl;
}
}
//_________________________________________
void HoughFit::DoNumerical() {
//
//to be implemented..
//
}