Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

HoughFit - source file

// 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;
}
About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update Wed Nov 16 16:57:39 2005 Wed Nov 16 16:29:22 2005