Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

HoughModule - source file


#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;
}


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