Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

HoughFit - source file

// 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..	
//
}



About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update 14-May-2005 21:31