// ESAF : Euso Simulation and Analysis Framework
// $Id: HoughFit.cc,v 1.8 2005/09/25 13:26:46 naumov 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>
using namespace TMath;
Double_t kHUGE = 1.e20;
ClassImp(HoughFit)
void HoughEstimator(Int_t &npar, Double_t *deriv, Double_t &f, Double_t *par, Int_t flag);
//_________________________________________________________________________________________________________________________
HoughFit::HoughFit( vector<Double_t> x, vector<Double_t> y, vector<Int_t> c, Int_t fDoSelection,Float_t ex,Float_t ey,vector<Int_t> id) {
//
// Constructor
//
string fKindHough = Conf()->GetStr("HoughFit.fKindHough");
fDoGraph = Conf()->GetBool("HoughFit.fDoGraph");
Clear();
fEpsX = ex;
fEpsY = ey;
fId = id;
fPointsX = x;
fPointsY = y;
fCounts = c;
fNumPoints = fPointsX.size();
if ( fNumPoints < 5 ) return;
if ( fKindHough == "sinusoidal" ) {
DoSinusoidal(1);
fCMaxToCAverage = fCMaxToCAverageFluttTmp;
Float_t minfrac = Conf()->GetNum( "HoughFit.fMinFractionSinusoidal" );
//cout<<"fFluttBackgHSpace="<<fFluttBackgHSpace<<"tfCMaxToCAverageFluttTmp="<<fCMaxToCAverageFluttTmp<<" fNumPointsSelected="<< fNumPointsSelected<<endl;
if( fCMaxToCAverage > minfrac*fFluttBackgHSpace && fDoSelection == 1 && fNumPointsSelected > 5 ){
//cout<<"There could be a signal, try to do selection..";
fPointsX = fXSel;
fPointsY = fYSel;
fCounts = fCSel;
fNumPoints = fPointsX.size();
Clear();
DoSinusoidal(2);
fCMaxToCAverageSel = fCMaxToCAverageFluttTmp;
//cout<<" .. now fFluttBackgHSpace="<<fFluttBackgHSpace<<"tfCMaxToCAverageFluttTmp="<<fCMaxToCAverageFluttTmp<<endl;
//cout<<" .. fCMaxToCAverageSel/fFluttBackgHSpace="<<fCMaxToCAverageSel/fFluttBackgHSpace<<endl;
if( fCMaxToCAverageSel > (minfrac+1)*fFluttBackgHSpace ) fThereIsSignal = kTRUE;
}
} else if( fKindHough == "linear" ) {
DoLinear(1);
fCMaxToCAverage = fCMaxToCAverageFluttTmp;
Float_t minfrac = Conf()->GetNum("HoughFit.fMinFractionLinear");
//cout<<"fFluttBackgHSpace="<<fFluttBackgHSpace<<"tfCMaxToCAverageFluttTmp="<<fCMaxToCAverageFluttTmp<<" fNumPointsSelected="<< fNumPointsSelected<<endl;
if( fCMaxToCAverage > minfrac*fFluttBackgHSpace && fDoSelection == 1 && fNumPointsSelected > 5 ){
//cout<<"There could be a signal, try to do selection..";
fPointsX = fXSel;
fPointsY = fYSel;
fCounts = fCSel;
fNumPoints = fPointsX.size();
Clear();
DoLinear(2);
fCMaxToCAverageSel = fCMaxToCAverageFluttTmp;
//cout<<" .. now fFluttBackgHSpace="<<fFluttBackgHSpace<<"tfCMaxToCAverageFluttTmp="<<fCMaxToCAverageFluttTmp<<endl;
//if( fCMaxToCAverageSel > fCMaxToCAverage ) fThereIsSignal=1;
//cout<<" .. fCMaxToCAverageSel/fFluttBackgHSpace="<<fCMaxToCAverageSel/fFluttBackgHSpace<<endl;
if( fCMaxToCAverageSel > (minfrac+1)*fFluttBackgHSpace ) fThereIsSignal = kTRUE;
}
} else if( fKindHough == "numerical" ){
//Msg(EsafMsg::Panic) << "DoNumerical() not implemented yet.." << MsgDispatch;
//return;
DoNumerical();
fThereIsSignal = kTRUE;
} else {
Msg(EsafMsg::Panic) << "Wrong fKindHough value in config file." << MsgDispatch;
}
}
//_______________________________________________________________________________
HoughFit::~HoughFit() {
//
//Destructor
//
fId.clear();
fPointsX.clear();
fPointsY.clear();
fPointsXprime.clear();
fPointsYprime.clear();
fCounts.clear();
fNumPointsSelected=0;
fNumPoints=0;
fIdSel.clear();
}
//________________________________________________________________________________
void HoughFit::Clear() {
//
// Clear method
//
fThereIsSignal = 0;
fXSel.clear();
fYSel.clear();
fCSel.clear();
fSlopeToPass = 0;
fOffsetToPass = 0;
fSlope = 0;
fOffset = 0;
fWidth = 0;
fAbsDev = 0;
}
//______________________________________________________________________________________
void HoughFit::PreProcess() {
//
// Preprocess for all kind of hough fit
//
fMinX = kHUGE;
fMinY = kHUGE;
fMaxX = -kHUGE;
fMaxY = -kHUGE;
fPointsXprime.clear();
fPointsYprime.clear();
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;
}
fX0 = (fMaxX+fMinX)/2;
fY0 = (fMaxY+fMinY)/2;
Float_t tmpfMinX_pre = kHUGE;
Float_t tmpfMinY_pre = kHUGE;
Float_t tmpfMaxX_pre = -kHUGE;
Float_t tmpfMaxY_pre = -kHUGE;
vector<Double_t> fPointsXprime_pre;
vector<Double_t> fPointsYprime_pre;
for( Int_t i=0; i<fNumPoints; i++ ) {
Double_t x = fPointsX[i]-fX0;
Double_t y = fPointsY[i]-fY0;
fPointsXprime_pre.push_back(x);
fPointsYprime_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;
}
tmpfMinX = kHUGE;
tmpfMinY = kHUGE;
tmpfMaxX = -kHUGE;
tmpfMaxY = -kHUGE;
fScale = tmpfMaxX_pre/tmpfMaxY_pre;
for( Int_t i=0; i<fNumPoints; i++ ) {
Double_t x = fPointsXprime_pre[i];
Double_t y = fPointsYprime_pre[i]*fScale;
fPointsXprime.push_back(x);
fPointsYprime.push_back(y);
if ( x > tmpfMaxX ) tmpfMaxX=x;
if ( y > tmpfMaxY ) tmpfMaxY=y;
if ( x < tmpfMinX ) tmpfMinX=x;
if ( y < tmpfMinY ) tmpfMinY=y;
}
}
//______________________________________________________________________________________
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)
//
PreProcess();
Int_t binmax(0), nx(0), binymax(0), binxmax(0);
Double_t thetamaxsingle(0), rhomaxsingle(0);
Float_t rhomax = Sqrt(tmpfMaxX*tmpfMaxX + tmpfMaxY*tmpfMaxY);
Float_t rhomin = -rhomax;
Float_t epsilon = Sqrt(fEpsX*fEpsX + fEpsY*fEpsY);
Float_t deltaro = 2*epsilon;
Float_t thetasinglemin = 0;
Float_t thetasinglemax = Pi();
Int_t binx = (Int_t)((tmpfMaxX - tmpfMinX)*1000.);
Int_t biny = (Int_t)((tmpfMaxY - tmpfMinY)*1000.);
Int_t numbinro = (Int_t)((rhomax - rhomin)/deltaro);
Float_t deltatheta = 2*deltaro/(rhomax - rhomin);
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,rhomin,rhomax);
TH2F *single2 = new TH2F("single2","single2",binx,tmpfMinX,tmpfMaxX,biny,tmpfMinY,tmpfMaxY);
for( Int_t i=0; i<fNumPoints; i++ ) {
Double_t x = fPointsXprime[i];
Double_t y = fPointsYprime[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);
rhomaxsingle = (single->GetYaxis())->GetBinCenter(binymax);
Float_t countmax=single->GetBinContent(binxmax,binymax);
/*TH2F *singlelt = new TH2F("singlelt","singlelt;theta;ro",numbintheta,thetasinglemin,thetasinglemax,numbinro,rhomin,rhomax);
TH2F *singlelt1 = new TH2F("singlelt1","singlelt1;theta;ro",numbintheta,thetasinglemin,thetasinglemax,numbinro,rhomin,rhomax);
TH2F *singlelt2 = new TH2F("singlelt2","singlelt2;theta;ro",numbintheta,thetasinglemin,thetasinglemax,numbinro,rhomin,rhomax);*/
Float_t consideredbin = 0;
Float_t counttotal = 0;
for(Int_t ith=1; ith < numbintheta; ith++) {
for(Int_t iro=1; iro<numbinro; iro++) {
Int_t count = (Int_t)single->GetBinContent(ith,iro);
/*if(count>0.25*countmax )singlelt->Fill((single->GetXaxis())->GetBinCenter(ith),(single->GetYaxis())->GetBinCenter(iro),count);
if(count>0.5*countmax )singlelt1->Fill((single->GetXaxis())->GetBinCenter(ith),(single->GetYaxis())->GetBinCenter(iro),count);
if(count>0.75*countmax )singlelt2->Fill((single->GetXaxis())->GetBinCenter(ith),(single->GetYaxis())->GetBinCenter(iro),count);*/
if(count< 0.9*countmax ){
counttotal = counttotal + count;
consideredbin++;
}
}
}
/*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;*/
Float_t countaverage = counttotal/(consideredbin);
fFluttBackgHSpace = Sqrt(countaverage);
fCMaxToCAverageTmp = (Float_t)(countmax/countaverage);
fCMaxToCAverageFluttTmp = (Float_t)(countmax/fFluttBackgHSpace);
TH1D *singleproj=new TH1D("singleproj","singleproj",numbinro,rhomin,rhomax);
singleproj=single->ProjectionY("singleproj",binxmax-5,binxmax+5);
fWidth = singleproj->GetRMS();
fSlope = -1/Tan(thetamaxsingle);
fOffset = rhomaxsingle/Sin(thetamaxsingle);
fSlopeToPass = fSlope/fScale;
fOffsetToPass = fY0+((fOffset-fSlope*fX0)/fScale);
Double_t xg[fNumPoints],yg[fNumPoints];
fAbsDev = 0;
Double_t markercenterx[1] = {thetamaxsingle};
Double_t markercentery[1] = {rhomaxsingle};
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] = fPointsXprime[i];
yg[i] = fSlope * fPointsXprime[i] + fOffset;
Float_t singlefAbsDev = Abs(yg[i] - fPointsYprime[i]);
fAbsDev += singlefAbsDev;
}
fAbsDev /= fNumPoints;
if( fDoGraph ) {
TGraph *g = new TGraph(fNumPoints,xg,yg); g->SetLineColor(8);
string namecHT;
if( firsttimedo==2 ) {
namecHT="cHTsel";
}else{
namecHT="cHT";
}
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 g;
}
delete single2; delete single; delete gmarkercenter; delete singleproj;
fNumPointsSelected=0;
if( firsttimedo==1 && fWidth>1e-4 ){
for( Int_t i=0; i<fNumPoints; i++ ) {
Double_t x = fPointsXprime[i];
Double_t y = fPointsYprime[i];
Int_t z = fCounts[i];
Double_t rosingleprova = x*cos(thetamaxsingle) + y*sin(thetamaxsingle);
if(fabs(rosingleprova-rhomaxsingle)<fWidth*0.3) {
fXSel.push_back(x+fX0);
fYSel.push_back((y/fScale)+fY0);
fCSel.push_back(z);
if( fId.size()>0 ) fIdSel.push_back(fId[i]);
fNumPointsSelected++;
}
}
}
}
//_______________________________________________________________________________
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
//
PreProcess();
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 = Sqrt(mmax*mmax * fEpsX*fEpsX + fEpsY*fEpsY);
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)*1000.);//not important binning, is just to draw
Int_t biny=(Int_t)((tmpfMaxY-tmpfMinY)*1000.);//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 = fPointsXprime[i];
Double_t y = fPointsYprime[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 countmax = single->GetBinContent(binxmax,binymax);
fSlope = mmaxsingle;
fOffset = qmaxsingle;
fSlopeToPass = fSlope/fScale;
fOffsetToPass = fY0+((fOffset-fSlope*fX0)/fScale);
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];
fAbsDev = 0;
for( Int_t i=0; i<fNumPoints; i++ ) {
xg[i] = fPointsXprime[i];
yg[i] = fSlope*fPointsXprime[i] + fOffset;
Float_t singlefAbsDev = fabs(yg[i] - fPointsYprime[i]);
fAbsDev += singlefAbsDev;
}
fAbsDev /= 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 counttotal = 0;
for(Int_t im=1;im<numbinm; im++) {
for(Int_t iq=1; iq<numbinq; iq++) {
Int_t count=(Int_t)(single->GetBinContent(im,iq));
/*if(count>0.25*countmax )singlelt->Fill((single->GetXaxis())->GetBinCenter(im),(single->GetYaxis())->GetBinCenter(iq),count);
if(count>0.5*countmax )singlelt1->Fill((single->GetXaxis())->GetBinCenter(im),(single->GetYaxis())->GetBinCenter(iq),count);
if(count>0.75*countmax )singlelt2->Fill((single->GetXaxis())->GetBinCenter(im),(single->GetYaxis())->GetBinCenter(iq),count);*/
if(count<0.9*countmax ){
counttotal += count;
consideredbin++;
}
}
}
Float_t countaverage = counttotal/(consideredbin);
fFluttBackgHSpace = Sqrt(countaverage);
fCMaxToCAverageTmp = (Float_t)(countmax/countaverage);
fCMaxToCAverageFluttTmp = (Float_t)(countmax/fFluttBackgHSpace);
TH1D *singleproj=new TH1D("singleproj","singleproj",numbinq,qmin,qmax);
singleproj=single->ProjectionY("singleproj",binxmax-5,binxmax+5);
fWidth = singleproj->GetRMS();
/* 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;*/
if(fDoGraph==1){
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 singleproj; delete g;
fNumPointsSelected=0;
if(firsttimedo==1 && fWidth>1e-4){
for( Int_t i=0; i<fNumPoints; i++ ) {
Double_t x = fPointsXprime[i];
Double_t y = fPointsYprime[i];
Int_t z = fCounts[i];
Double_t qsingleprova = -mmaxsingle*x+y;
if(Abs(qsingleprova - qmaxsingle) < fWidth*0.3) {
fXSel.push_back(x+fX0);
fYSel.push_back((y/fScale)+fY0);
fCSel.push_back(z);
if( fId.size()>0 ) fIdSel.push_back(fId[i]);
fNumPointsSelected++;
}
}
//cout<<"nfNumPoints="<<fNumPoints<<"tfNumPointsSelected="<<fNumPointsSelected<<"tfWidth="<<fWidth<<endl;
}
}
//______________________________________________________________________________________
void HoughFit::DoNumerical() {
//
// Numerical Hough Fit
// to be implemented..
//
PreProcess();
// fill the container for data
HoughData fData;
fData.fNumPoints = fNumPoints;
fData.fX = fPointsXprime;
fData.fY = fPointsYprime;
fData.fCounts = fCounts;
fData.fRadius = 0.8; //FIXME FIXME FIXME
// initialize the fit
TMinuit *fHoughFit = new TMinuit(2);
fHoughFit->SetPrintLevel(-1);
fHoughFit->SetObjectFit(&fData);
fHoughFit->SetFCN(HoughEstimator);
fHoughFit->SetErrorDef(1.);
// initialize parameters
Double_t par[2];
Double_t XSum = 0;
Double_t XSqSum = 0;
Double_t YSum = 0;
Double_t XYSum = 0;
for( Int_t i=0; i<fNumPoints; i++ ) {
XSum += fPointsXprime[i];
XSqSum += fPointsXprime[i] * fPointsXprime[i];
YSum += fPointsYprime[i];
XYSum += fPointsXprime[i] * fPointsYprime[i];
}
// method of least squares for a first guess for slope and offset
par[0] = (fNumPoints * XYSum - XSum * YSum) / (fNumPoints * XSqSum - XSum * XSum); // slope
par[1] = (XSqSum * YSum - XSum * XYSum) / (fNumPoints * XSqSum - XSum * XSum); // offset
par[2] = 0.8; //FIXME
Double_t min[3] = {0.0, 0.0, 0.0}; //FIXME
Double_t max[3] = {0.0, 0.0, 0.0}; //FIXME
string cpar[3] = {"m","q","r"};
Double_t step[3];
/*Double_t step[3] = {0.000, 0.000, 0.001}; //FIXME
for( Int_t i=0; i<3; i++ ) {
fHoughFit->DefineParameter(i,cpar[i].c_str(),par[i],step[i],min[i],max[i]);
}
fHoughFit->Migrad();
Double_t err_r;
fHoughFit->GetParameter(2,par[2],err_r);*/
step[0] = 0.001; step[1] = 0.001; step[2] = 0.000; //FIXME
for( Int_t i=0; i<3; i++ ) {
fHoughFit->DefineParameter(i,cpar[i].c_str(),par[i],step[i],min[i],max[i]);
}
// retrieve results
fHoughFit->Migrad();
Double_t outpar[2],err[2];
for( Int_t i=0; i<2; i++ ) {
fHoughFit->GetParameter(i,outpar[i],err[i]);
}
fSlope = outpar[0];
fOffset = outpar[1];
fSlopeToPass = fSlope/fScale;
fOffsetToPass = fY0+((fOffset-fSlope*fX0)/fScale);
fWidth = Abs(0.3*fSlopeToPass); // FIXME
// selection of points
fNumPointsSelected=0;
if(fWidth>1e-4){
for( Int_t i=0; i<fNumPoints; i++ ) {
Double_t x = fPointsXprime[i];
Double_t y = fPointsYprime[i];
Int_t z = fCounts[i];
Double_t qsingleprova = -fSlope*x+y;
if(Abs(qsingleprova - fOffset) < 3.*fWidth) { //FIXME
fXSel.push_back(x+fX0);
fYSel.push_back((y/fScale)+fY0);
fCSel.push_back(z);
if( fId.size()>0 ) fIdSel.push_back(fId[i]);
fNumPointsSelected++;
}
}
}
fData.Clear();
}
//______________________________________________________________________________________
void HoughData::Clear() {
//
// Clear container of data
//
fX.clear();
fY.clear();
fCounts.clear();
fRadius = -1.;
fNumPoints = 0;
}
//______________________________________________________________________________
void HoughEstimator(Int_t &npar, Double_t *deriv,Double_t &f, Double_t *par, Int_t flag) {
f = 0;
HoughData *data = (HoughData*)gMinuit->GetObjectFit();
Double_t sum = 0;
Int_t tot_counts = 0;
for( Int_t i(0); i<data->fNumPoints; i++ ) {
if ( Power( Abs(data->fX[i]*par[0] + par[1] - data->fY[i]),2)
<= Power(par[2],2) * ( Power(data->fX[i],2) + 1 ) ) sum += data->fCounts[i];
tot_counts += data->fCounts[i];
}
f = tot_counts / sum;
}