Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

HoughModule - source file


// $Id: HoughModule.cc
// Author: Elena Taddei   

/*****************************************************************************
 * ESAF: Euso Simulation and Analysis Framework                              *
 *                                                                           *
 *  Id: HoughModule                                                          *
 *                                                                           *
 *****************************************************************************/

//_____________________________________________________________________________
//
// HoughModule
//
// <extensive class description>
//
//   Config file parameters
//   ======================
//
//   <parameter name>: <parameter description>
//   -Valid options: <available options>
//

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

using namespace TMath;

ClassImp(HoughModule)

// FIXME memory leaks to be fixed asap. D.Naumov

//_____________________________________________________________________________
 HoughModule::HoughModule():RecoModule("HoughTransform") {
    //
    // ctor
    //
}

//_____________________________________________________________________________
 HoughModule::~HoughModule() {
    //
    // dtor
    //
}

//_____________________________________________________________________________
 Bool_t HoughModule::Init() {
    //
    // init
    //
    
    fEv = (RecoEvent*)0;
    
    fNumPointsMinimum = (Int_t)Conf()->GetNum("HoughModule.NumPointsMinimum");
    fProcedure = Conf()->GetStr("HoughModule.Procedure");
    Msg(EsafMsg::Info) << "Initializing... Procedure: ""<< fProcedure.c_str() << """ << MsgDispatch;
    
    return kTRUE;
}

//_____________________________________________________________________________
 Bool_t HoughModule::PreProcess() {
    //
    // pre-process
    //
	
    fTotalNumPoints = 0;
    fNumPointsAll = 0;
   
    fHittedFinal.clear();
    fHittedFinalMacrocell.clear();
    fHittedSelHough.clear();
    fHittedSel2Hough.clear();
    fHittedAll.clear();
    fIdMacrocellsSelected.clear();
       
    fClusters.clear();
        
    return kTRUE;
}

// process
//_____________________________________________________________________________
 Bool_t HoughModule::Process( RecoEvent* ev ) {
	
    fEv = ev;
    fTotalNumPoints = fEv->GetHeader().GetNumActiveFee();
    fIdEv = fEv->GetHeader().GetNum();
	
    
    if ( fTotalNumPoints <= fNumPointsMinimum ) {       
        Msg(EsafMsg::Info) << "Not enough points in event! Exit from this event." << MsgDispatch;
        return kFALSE;
    }else{
        Msg(EsafMsg::Debug) << "fTotalNumPoints = " << fTotalNumPoints << MsgDispatch;
    }
    
    fMainMacroCellId = fEv->GetMainMacroCell(1);
    if ( fEv->GetRecoCellInfo(fMainMacroCellId) == NULL )
        Msg(EsafMsg::Info) << "Main macrocell does not correspond to any macrocell hit." << MsgDispatch;
      
    Bool_t ProcessOnlySignal = Config::Get()->GetCF("Reco","RootInputModule")->GetBool("RootInputModule.ProcessOnlySignal");
    if (ProcessOnlySignal)  {
        for( Int_t i=0; i<fTotalNumPoints; i++ ) {
             RecoPixelData *pix = (RecoPixelData*)fEv->GetRecoPixelData(i);
             Int_t signal_counts = pix->GetCountsSignal();
             if ( signal_counts > 0 ) fHittedFinal.push_back(i);
        }
        fFinalNumHitsMinimum=1;
        WriteSingleCluster();
        return kTRUE;
    }
    gDirectory->cd("/");
    string pathnameevent = Form("HoughModule%d", fIdEv);
    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){
        Msg(EsafMsg::Info) << "Problems ret=0.. " << MsgDispatch;
    	return kFALSE;
    }
    if(ret==1){
    	if(fHittedFinal.size()<10 || fHittedFinal.size()>100000)	return kFALSE;
	WriteSingleCluster();
    	return kTRUE;
    } 
    return kFALSE; 
    
}

//_____________________________________________________________________________
 Int_t HoughModule::ProcessHoughSingleMacrocell(){
    //
    // Do Hough Transform for each macrocell hitted and trying to select 
    // which are hitted by real signal
    //
	
    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> empty;
    		pair< Int_t,vector<Int_t> > p(i,empty);
    		mapIdfHittedSingleCell.insert(p);
    	}
    }
    Msg(EsafMsg::Info) << "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);
	}
    }
   
    Msg(EsafMsg::Info) << "selected pixel with counts > " << numcountsmin << MsgDispatch;
    for( it = mapIdfHittedSingleCell.begin(); it != mapIdfHittedSingleCell.end(); it++ ){
    	
    	Msg(EsafMsg::Info)<<"n\tt---->macrocellid = "<<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;
		
	TH1F *h1 = new TH1F("h1","h1",10,0,10);
    	for(Int_t i=0;i<(Int_t)fHittedmc.size();i++){
    		h1->Fill(fEv->GetRecoPixelData(fHittedmc[i])->GetCounts());
    		Int_t pixid=fEv->GetRecoPixelData(fHittedmc[i])->GetPixelId();
		TVector3 dummy2 = runpars->PixelCenter(pixid);
    		xmc.push_back(0.001*dummy2.x());//m
    		ymc.push_back(0.001*dummy2.y());//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/", fIdEv);
    	gDirectory->cd(pathname.c_str());
	string pathnames = Form("SingleHoughMacrocell%d", idmacrocell);
    	TDirectory *path = new TDirectory(pathnames.c_str(),pathnames.c_str());
    	path->cd();
    	
    	h1->Draw();
	h1->Write();
	delete h1;
    	
    	Bool_t possiblesignal = HoughTransform(xmc,ymc,tmc,cmc);
    	
    	if( possiblesignal ) {
    		Msg(EsafMsg::Info) << "Macrocell " << idmacrocell << " probably HAS signal!!!" << MsgDispatch;
    		fIdMacrocellsSelected.push_back(idmacrocell);
    		fHittedFinalMacrocell.insert(fHittedFinalMacrocell.end(),fHittedmc.begin(),fHittedmc.end());
    	}else{
    		Msg(EsafMsg::Info) << "Macrocell " << idmacrocell << " probably DOESN'T HAVE a signal,";
    		if ( idmacrocell==fMainMacroCellId ){
    			Msg(EsafMsg::Info) << " but is added because is the main!" << MsgDispatch;
    			fHittedFinalMacrocell.insert(fHittedFinalMacrocell.end(),fHittedmc.begin(),fHittedmc.end());
    			fIdMacrocellsSelected.push_back(idmacrocell);
    		} else {
    			Msg(EsafMsg::Info) << " discarded!" << MsgDispatch;
    		}
    		
    	}
    	
    	xmc.clear();
    	ymc.clear();
    	tmc.clear();
    	cmc.clear();
    }
    
    if (fIdMacrocellsSelected.size()<1 || fHittedFinalMacrocell.size()<10)	return 0;
    Msg(EsafMsg::Info) << "Re-process " << fIdMacrocellsSelected.size() << " selected macrocells by ProcessHoughSingleMacrocell()" << MsgDispatch;  
    Int_t ret = ProcessHoughAll();

    return ret;    
}

//_____________________________________________________________________________
 Int_t HoughModule::ProcessVelocity(){
    //
    // Do histograms of velocities for each macrocell hitted and trying to select 
    // which are hitted by real signal
    //
	
    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> empty;
    		pair< Int_t,vector<Int_t> > p(i,empty);
    		mapIdfHittedSingleCell.insert(p);
    	}
    }
    Msg(EsafMsg::Info) << "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 = 20;
    Msg(EsafMsg::Info) << "selected pixel with counts > " << numcountsmin << MsgDispatch;
    for( it = mapIdfHittedSingleCell.begin(); it != mapIdfHittedSingleCell.end(); it++ ) {
    	Int_t entriesTH2(0);    	
    	Msg(EsafMsg::Info) << "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); 
    	TH2F *sxy = new TH2F("sxy","velocita",binv,-2,2,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 pointstocouple = fHittedmc.size();
    	    for(Int_t j=i+1;j<pointstocouple;j++){
    	        Int_t pixtocouple = fEv->GetRecoPixelData(fHittedmc[j])->GetPixelId();
		TVector3 dummycouple = runpars->PixelCenter(pixtocouple);
		Int_t tmccouple = fEv->GetRecoPixelData(fHittedmc[j])->GetGtu();
		if (Abs(dummy2.x()-dummycouple.x())>0 && Abs(dummy2.y()-dummycouple.y())>0 && tmccouple-tmcgtu>0 && tmccouple-tmcgtu<20 && Abs(dummy2.x()-dummycouple.x())<100 && Abs(dummy2.y()-dummycouple.y())<100 ) {
		    Int_t countproduct = (Int_t)Sqrt(1.e0*counts*fEv->GetRecoPixelData(fHittedmc[j])->GetCounts());
		    Double_t m = (dummy2.y()-dummycouple.y())/(dummy2.x()-dummycouple.x());
	            Double_t v = Sqrt((dummy2.y()-dummycouple.y())*(dummy2.y()-dummycouple.y())+(dummy2.x()-dummycouple.x())*(dummy2.x()-dummycouple.x()))/((tmccouple-tmcgtu)*2.5);// mm/us
			Double_t vx = (dummy2.x()-dummycouple.x())/((tmcgtu-tmccouple)*2.5);
			Double_t vy = (dummy2.y()-dummycouple.y())/((tmcgtu-tmccouple)*2.5);
			entriesTH2++;
			s->Fill(v,m,countproduct);
			sx->Fill(vx,countproduct);
			sy->Fill(vy,countproduct);
			sxy->Fill(vx,vy,countproduct);
		}
    	    }
    	}    
    
        Int_t binmaxv = sxy->GetMaximumBin();
        Int_t nxv = (sxy->GetXaxis())->GetNbins()+2;
        Int_t binxmaxv = binmaxv % nxv;
        Int_t binymaxv = binmaxv / nxv;
        Float_t vxmax[1], vymax[1];
        vxmax[0] = (sxy->GetXaxis())->GetBinCenter(binxmaxv);
        vymax[0] = (sxy->GetYaxis())->GetBinCenter(binymaxv);
        TGraph *gvmax=new TGraph(1,vxmax,vymax);
        gvmax->SetMarkerColor(2); gvmax->SetMarkerStyle(29); gvmax->SetMarkerSize(1.6);
        //Float_t contmaxv=sxy->GetBinContent(binxmaxv,binymaxv);
        Stat_t countx(0),county(0);
        for(Int_t bin=1;bin<binv;bin++){
    	    countx += sx->GetBinContent(bin);
        	county += sy->GetBinContent(bin);
        }
        Float_t countmaxx = sx->GetMaximum();
        Float_t countmaxy = sy->GetMaximum();
        Float_t countmeanx = (Float_t)((Float_t)countx/binv);
        Float_t countmeany = (Float_t)((Float_t)county/binv);
        Float_t ms = 20.;
        /*cout<<"countmeanx="<<countmeanx<<"tcontmaxx="<<countmaxx<<"tcontmaxx-countmeanx="<<countmaxx-countmeanx<<"tsig="<<ms*TMath::Sqrt(countmeanx)<<endl;
    	cout<<"countmeany="<<countmeany<<"tcontmaxy="<<countmaxy<<"tcontmaxy-countmeany="<<countmaxy-countmeany<<"tsig="<<ms*TMath::Sqrt(countmeany)<<endl;*/
        string a = "a"; string b = "b";
        string pathname = Form("/HoughModule%d/", fIdEv);
        gDirectory->cd(pathname.c_str());
        string pathnames;
        if( ((countmaxx-countmeanx) > ms*Sqrt(countmeanx) || (countmaxy-countmeany) > ms*Sqrt(countmeany)) && sx->GetEntries()>10 ) {
            fHittedFinalMacrocell.insert(fHittedFinalMacrocell.end(),fHittedmc.begin(),fHittedmc.end());
    	    fIdMacrocellsSelected.push_back(idmacrocell);
            Msg(EsafMsg::Info) << "Macrocell " << idmacrocell << " probably HAS signal!!!" << MsgDispatch;
            if (idmacrocell==fMainMacroCellId) {
    	        pathnames = Form("VelocityMacrocellmain%d", idmacrocell);
            } else {
                pathnames = Form("VelocityMacrocell%d", idmacrocell);
    	    }
        } else{
            Msg(EsafMsg::Info) << "Macrocell " << idmacrocell << " probably DOESN'T HAVE a signal,";
	    if (idmacrocell==fMainMacroCellId) {
	        pathnames = Form("VelocityMacrocellmain%d", idmacrocell);
    	        Msg(EsafMsg::Info) <<" but is added because is the main!"<< MsgDispatch;
    	        fHittedFinalMacrocell.insert(fHittedFinalMacrocell.end(),fHittedmc.begin(),fHittedmc.end());
    	        fIdMacrocellsSelected.push_back(idmacrocell);
            } else {
    	        pathnames = Form("VelocityMacrocelldiscarded%d", idmacrocell);
       	        Msg(EsafMsg::Info) <<" 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,2);
        bc->cd(1);hsxt->Draw(); 
        bc->cd(2);hsyt->Draw(); 
        bc->cd(3);sxy->Draw();gvmax->Draw("same"); 
        bc->cd(4);sxy->Draw("LEGO"); 
        bc->Write();
        delete bc;	delete hsxt;	delete hsyt;
        delete sxy;delete gvmax;
    }
    
    if (fIdMacrocellsSelected.size()<1 || fHittedFinalMacrocell.size()<10)	return 0;
	
    Msg(EsafMsg::Info) << "Re-process " << fIdMacrocellsSelected.size() << " selected macrocells by ProcessVelocity()" << MsgDispatch;  
    Int_t ret = ProcessHoughAll();
    
    return ret;

}

//_____________________________________________________________________________
 Int_t HoughModule::ProcessHoughAll(){
    //
    // Doing Hough Transform to all macrocells selected by previous 
    // methods or by euso trigger
    //
	
    Msg(EsafMsg::Info) << " ...ProcessHoughAll()" << MsgDispatch;
	
    string pathname = Form("/HoughModule%d/", fIdEv);
    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::Info) << "Let's try beginning with MinCountsHoughAll = " << fNumHitsMinimum << MsgDispatch;
    }
    
    Int_t regr=0;
    
    Float_t gtumin = 1000.;
    Float_t gtumax = 0.;
    Float_t deltagtu = 0.;
    
    fNumPointsAll = 0;
    fHittedAll.clear();
    if(!(fHittedFinalMacrocell.size())){
	Msg(EsafMsg::Info) << "processing all macrocells!" << MsgDispatch;
	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);
        }
    } else {
        Msg(EsafMsg::Info)<<"processing fHittedFinalMacrocell.size()="<<fHittedFinalMacrocell.size()<<MsgDispatch;
        for( Int_t i=0; i<(Int_t)(fHittedFinalMacrocell.size()); i++ ){
    	    Int_t i2=fHittedFinalMacrocell[i];
    	    if (fEv->GetRecoPixelData(i2)->GetGtu()<gtumin) gtumin=fEv->GetRecoPixelData(i2)->GetGtu();
    	    if (fEv->GetRecoPixelData(i2)->GetGtu()>gtumax) gtumax=fEv->GetRecoPixelData(i2)->GetGtu();
    	    if (fEv->GetRecoPixelData(i2)->GetCounts() >= fNumHitsMinimum) fHittedAll.push_back(i2);
        }
    }
        
    fNumPointsAll = fHittedAll.size();
        
    fTmin = gtumin*2.5*0.01; //FIXME
    fTmax = gtumax*2.5*0.01; //FIXME
    deltagtu=gtumax-gtumin;
   
    Msg(EsafMsg::Info) << "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,fTmin,fTmax,bin,-1,1);
	TH2F *single2yt = new TH2F("single2yt","single2yt;t;y(m)",bin,fTmin,fTmax,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(0.001*dummy2.x());
    	y.push_back(0.001*dummy2.y());
    	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]);
    }
    
    
    TCanvas *HMinit=new TCanvas("HMinit","HMinit",1);
    HMinit->Divide(3,1);
    HMinit->cd(1);  single2->SetMarkerStyle(6);  single2->Draw();
    HMinit->cd(2);  single2xt->SetMarkerStyle(6);  single2xt->Draw();
    HMinit->cd(3);  single2yt->SetMarkerStyle(6);  single2yt->Draw();
    HMinit->Write();
    delete HMinit;
    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,fXmin,fXmax,bin,fYmin,fYmax);
    TH2F *single2xtsel = new TH2F("single2xtsel","single2xtsel;t;x(m)",bin,fTmin,fTmax,bin,fXmin,fXmax);
    TH2F *single2ytsel = new TH2F("single2ytsel","single2ytsel;t;y(m)",bin,fTmin,fTmax,bin,fYmin,fYmax);
    for( Int_t i=0; i<fTotalNumPoints; i++ ) {
    	Int_t pixid=fEv->GetRecoPixelData(i)->GetPixelId();
	TVector3 dummy2 = runpars->PixelCenter(pixid);
    	Double_t x=0.001*dummy2.x();
    	Double_t y=0.001*dummy2.y();
    	Double_t t=(Double_t)fEv->GetRecoPixelData(i)->GetGtu()*2.5*0.01;
    	Int_t c=fEv->GetRecoPixelData(i)->GetCounts();
        if ( x>=fXmin && x<=fXmax && y>=fYmin && y<=fYmax && 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::Info) << "first selected zone: pixels with at least " << fNumHitsMinimum-regr << " hits = " << fHittedSelHough.size() << MsgDispatch;
    if(fHittedSelHough.size()<10){
    	Msg(EsafMsg::Info) << "too few! exit from this event" <<MsgDispatch;
    	return kFALSE;
    }
    
    
    TCanvas *HMsel=new TCanvas("HMsel","HMsel",1);
    HMsel->Divide(3,1);
    HMsel->cd(1);  single2sel->SetMarkerStyle(6);  single2sel->Draw();
    HMsel->cd(2);  single2xtsel->SetMarkerStyle(6);  single2xtsel->Draw();
    HMsel->cd(3);  single2ytsel->SetMarkerStyle(6);  single2ytsel->Draw();
    HMsel->Write();
    delete HMsel;
    delete single2xtsel;
    delete single2ytsel;
    delete single2sel;
    
    HoughTransform(xsel,ysel,tsel,csel); 
    
    xsel.clear();  ysel.clear();  tsel.clear();  csel.clear();
    if(fNumHitsMinimum-regr > 2)	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,fXmin,fXmax,bin,fYmin,fYmax);
    TH2F *single2xtselsel = new TH2F("single2xtselsel","single2xtselsel;t;x(m)",bin,fTmin,fTmax,bin,fXmin,fXmax);
    TH2F *single2ytselsel = new TH2F("single2ytselsel","single2ytselsel;t;y(m)",bin,fTmin,fTmax,bin,fYmin,fYmax);
    for( Int_t i=0; i<fTotalNumPoints; i++ ) {
    	Int_t pixid=fEv->GetRecoPixelData(i)->GetPixelId();
	TVector3 dummy2 = runpars->PixelCenter(pixid);
    	Double_t x=0.001*dummy2.x();
    	Double_t y=0.001*dummy2.y();
    	Double_t t=(Double_t)fEv->GetRecoPixelData(i)->GetGtu()*2.5*0.01;
    	Int_t c=fEv->GetRecoPixelData(i)->GetCounts();
        if ( x>=fXmin && x<=fXmax && y>=fYmin && y<=fYmax && 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);
        	fHittedSel2Hough.push_back(i);
        }
    }
    Msg(EsafMsg::Info) << "second selected zone: pixels with at least " << fNumHitsMinimum-regr << " hits = " << fHittedSel2Hough.size()<< MsgDispatch;
    if(fHittedSel2Hough.size()<10){
    	Msg(EsafMsg::Info) << "too few! exit from this event" <<MsgDispatch;
    	return kFALSE;
    }
    
    TCanvas *HMselsel=new TCanvas("HMselsel","HMselsel",1);
    HMselsel->Divide(3,1);
    HMselsel->cd(1);  single2selsel->SetMarkerStyle(6);  single2selsel->Draw();
    HMselsel->cd(2);  single2xtselsel->SetMarkerStyle(6);  single2xtselsel->Draw();
    HMselsel->cd(3);  single2ytselsel->SetMarkerStyle(6);  single2ytselsel->Draw();
    HMselsel->Write();
    delete HMselsel;
    delete single2xtselsel;
    delete single2ytselsel;
    delete single2selsel;
    
    
    HoughTransform(xselsel,yselsel,tselsel,cselsel,fHittedSel2Hough); 
    fHittedFinal=fHittedHoughMerged;
    xselsel.clear();  yselsel.clear();  tselsel.clear();  cselsel.clear();
    
    vector<Double_t> xend;
    vector<Double_t> yend;
    vector<Double_t> tend;
    vector<Int_t> cend;
    
    TH2F *single2end = new TH2F("single2end","single2end;x(m);y(m)",bin,fXmin,fXmax,bin,fYmin,fYmax);
    TH2F *single2xtend = new TH2F("single2xtend","single2xtend;t;x(m)",bin,fTmin,fTmax,bin,fXmin,fXmax);
	TH2F *single2ytend = new TH2F("single2ytend","single2ytend;t;y(m)",bin,fTmin,fTmax,bin,fYmin,fYmax);
	for( Int_t i=0; i<(Int_t)fHittedHoughMerged.size(); i++ ) {
    	Int_t pixid=fEv->GetRecoPixelData(fHittedHoughMerged[i])->GetPixelId();
		TVector3 dummy2 = runpars->PixelCenter(pixid);
    	Double_t x=0.001*dummy2.x();
    	Double_t y=0.001*dummy2.y();
    	Double_t t=(Double_t)fEv->GetRecoPixelData(fHittedHoughMerged[i])->GetGtu()*2.5*0.01;
    	Int_t c=fEv->GetRecoPixelData(fHittedHoughMerged[i])->GetCounts();
        single2end->Fill(x,y,c);
      	single2xtend->Fill(t,x,c);
       	single2ytend->Fill(t,y,c);
       	xend.push_back(x);
       	yend.push_back(y);
       	tend.push_back(t);
       	cend.push_back(c);
    }
    Msg(EsafMsg::Info) << "final selected zone: pixels with at least " << fNumHitsMinimum-regr << " hits = " << fHittedHoughMerged.size()<< MsgDispatch;
    if(fHittedHoughMerged.size()<10){
    	Msg(EsafMsg::Info) << "too few! exit from this event" <<MsgDispatch;
    	return kFALSE;
    }
    
    TCanvas *HMfinal=new TCanvas("HMfinal","HMfinal",1);
    HMfinal->Divide(3,1);
    HMfinal->cd(1);  single2end->SetMarkerStyle(6);  single2end->Draw();
    HMfinal->cd(2);  single2xtend->SetMarkerStyle(6);  single2xtend->Draw();
    HMfinal->cd(3);  single2ytend->SetMarkerStyle(6);  single2ytend->Draw();
    HMfinal->Write();
    delete HMfinal;
    delete single2xtend;
    delete single2ytend;
    delete single2end;
    xend.clear();yend.clear();tend.clear();cend.clear();
    fFinalNumHitsMinimum=fNumHitsMinimum-regr;
    return 1;
    
}


//_____________________________________________________________________________
 Bool_t HoughModule::HoughTransform(vector<Double_t> x,vector<Double_t> y,vector<Double_t> t,vector<Int_t> c,vector<Int_t> id){
    //
    // Implements the Hough Transform
    // Return true if there is a possible signal
    //

    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,id);
    //Double_t vx = xtFit->GetSlope();//Double_t wx = xtFit->GetOffset();
    fXmin = xtFit->GetMinY();
    fXmax = xtFit->GetMaxY();
    Double_t tmin1 = xtFit->GetMinX();
    Double_t tmax1 = xtFit->GetMaxX();
    fAbsDevX = xtFit->GetAbsDev();
    Bool_t possiblesignalxt = xtFit->IsSignal();
    fCmaxToCmeanX = xtFit->GetMaxToMedium();
    vector<Int_t> idselx = xtFit->GetIdSel();
    delete xtFit;
	
    HoughFit *ytFit = new HoughFit(t, y, c, optionselection,et,ey,id);
    //Double_t vy = ytFit->GetSlope();//Double_t wy = ytFit->GetOffset();
    fYmin = ytFit->GetMinY();
    fYmax = ytFit->GetMaxY();
    Double_t tmin2 = ytFit->GetMinX();
    Double_t tmax2 = ytFit->GetMaxX();
    fAbsDevY = ytFit->GetAbsDev();
    Bool_t possiblesignalyt = ytFit->IsSignal();
    fCmaxToCmeanY = ytFit->GetMaxToMedium();
    vector<Int_t> idsely = ytFit->GetIdSel();
    delete ytFit;
	
    if(id.size()>0) {
	Int_t sumsel=idsely.size()+idselx.size();
	//cout<<"id.size()="<<id.size()<<"tidselx.size()="<<idselx.size()<<"tidsely.size()="<<idsely.size()<<endl;
        vector<Int_t> idselectedmerged(sumsel);
	merge(idselx.begin(), idselx.end(),idsely.begin(), idsely.end(),idselectedmerged.begin());
	vector<Int_t>::iterator diversi=unique(idselectedmerged.begin(),idselectedmerged.end());
	idselectedmerged.erase(diversi,idselectedmerged.end());
	//cout<<"idselectedmerged.size()="<<idselectedmerged.size()<<endl;
	fHittedHoughMerged=idselectedmerged;
    }
	
    //Msg(EsafMsg::Info) <<"possiblesignalxt="<<possiblesignalxt<<"tpossiblesignalyt="<<possiblesignalyt<<MsgDispatch;
    //Msg(EsafMsg::Info) <<"fCmaxToCmeanX="<<fCmaxToCmeanX<<"tcmaxtocmedioy="<<fCmaxToCmeanY<<MsgDispatch;
    Bool_t possible = kFALSE;
    if( possiblesignalxt || possiblesignalyt ) possible = kTRUE;
	
    fTmin = (tmin1<tmin2) ? tmin1 : tmin2;
    fTmax = (tmax1<tmax2) ? tmax2 : tmax1;
	
    x.clear();  y.clear();  t.clear();  c.clear();
    
    return possible;
    
}

//_____________________________________________________________________________
 Bool_t HoughModule::PostProcess() {
    //
    // post-process
    //
    
    if ( !fTotalNumPoints ) return kTRUE;
    
    fEv = (RecoEvent*)0;
    return kTRUE;
    	
}

//____________________________________________________________________________
 Bool_t HoughModule::SaveRootData(RecoRootEvent *fRecoRootEvent) {
    //
    // Save reco root data
    //
    
    return kTRUE;
}

//____________________________________________________________________________
 Bool_t HoughModule::Done() {
    //
    // done
    //
    
    fHittedAll.clear();
    fHittedFinal.clear();
    fHittedSelHough.clear();
    fHittedSel2Hough.clear();
    fIdMacrocellsSelected.clear();
   
    fClusters.clear();
    Msg(EsafMsg::Info) << "Completed" << MsgDispatch;
    return kTRUE;
}

//_____________________________________________________________________________
 void HoughModule::UserMemoryClean() {
    //
    // memory clean
    //
	
    Msg(EsafMsg::Debug) << "UserMemoryClean()" << MsgDispatch;
	
    Int_t fHittedFinalsize = fHittedFinal.size();

    if(fHittedFinalsize<10 || !fTotalNumPoints || fHittedFinalsize>100000)		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");
    
    vector<Int_t> *pix = (vector<Int_t>*)MyData()->GetObj("CluPixels");
    pix->clear();
    MyData()->RemoveObj("CluPixels");
   
}

//_____________________________________________________________________________
 void HoughModule::WriteSingleCluster() {
    //
    // close and write the current cluster if is significant
    //
    
    // 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;
    
    // add data to event
    vector<EusoCluster*> *pclu = &fClusters;
    MyData()->Add("Clusters", pclu);
    MyData()->Add("NumHitsMinimum",fFinalNumHitsMinimum);
    vector<Int_t> *pix = &fHittedFinal;
    MyData()->Add("CluPixels", pix);
}

//_____________________________________________________________________________
 EusoCluster* HoughModule::GetCluster( Int_t i ) {
    //
    // get a specified cluster
    //

    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 Wed Nov 16 16:57:39 2005 Wed Nov 16 16:29:22 2005