source: JEM-EUSO/esaf_cc_at_lal/macros/selftest/CheckHoughTransform.C @ 114

/*
To run this macro you should set the kDevelop mode in rootlogon.C and do:

root [] .x CheckHoughTransform.C+               

then answer to questions,
to choose the kind of HoughTransform to use look in 
file config/Reco/HoughModule.cfg and change 
HoughModule.KindHough=linear or sinusoidal


*/
#include <vector>
#include "EsafRandom.hh"
#include "HoughFit.hh"

#include <iostream>

#include "TBrowser.h"
#include "TH1.h"
#include "TH2F.h"
#include "TVector3.h"
#include "TGraph.h"
#include "TCanvas.h"

using namespace std;
void DoHoughTransform(Int_t dohough);
 
//______________________________________________________________________________
void CheckHoughTransform() {
        Int_t dohough;
        cout<<"\nHoughFit of signal+background(1),signal(2),background(3)"<<endl;
        cout<<"Which do you want?  ";
        cin>>dohough;
        if( dohough!=1 && dohough!=2 && dohough!=3 ){
                cout<<"Not valid choice, return"<<endl;
                return;
        }
        DoHoughTransform(dohough);
}

//______________________________________________________________________________
void DoHoughTransform(Int_t dohough) {
        
        TRandom* rndm = EsafRandom::Get();
        
        
        Double_t m=-1+2*rndm->Rndm();
        Double_t frac=0;
        string nameroot;
        if(dohough==1 || dohough==2){
                cout<<"Fraction of signal (between 0 and 1)= ";
                cin>>frac;
                if(dohough==1)  nameroot="hfbkgsign.root";
                if(dohough==2)  nameroot="hfsign.root";
        }else{
                nameroot="hfbkg.root";
        }
        TFile *f = TFile::Open(nameroot.c_str(), "RECREATE" );
        
        Float_t ex=0.025;//error=1gtu
        Float_t ey=0.004;//error=4mm
        Double_t x1=0;//tempo se metto gtu*2.5*0.01
        Double_t x2=2;//tempo se metto gtu*2.5*0.01
        
        /*
        Float_t ex=0.004;//error=4mm
        Float_t ey=0.004;//error=4mm
        Double_t x1=-1;//m
        Double_t x2=1;//m
        */
        Double_t y1=-1;//m
        Double_t y2=1;//m

        Double_t xs1=rndm->Rndm()+x1;
        Double_t xs2=xs1+1;
        Double_t ys1=y1+(rndm->Rndm())*(y2-y1);
        Double_t ys2=ys1+m*(xs2-xs1);
        Double_t q=ys1-m*xs1;
        
        Int_t bin=100;
        TH2F *hrndmbkgsign=new TH2F("hrndmbkgsign","hrndmbkgsign",bin,x1,x2,bin,y1,y2); 
        
        vector<Double_t> xbkg;
        vector<Double_t> ybkg;
        vector<Int_t> cbkg;
        vector<Double_t> xsign;
        vector<Double_t> ysign;
        vector<Int_t> csign;
        vector<Double_t> xbkgsign;
        vector<Double_t> ybkgsign;
        vector<Int_t> cbkgsign;
        Double_t meangaussr=0.005;
        Double_t sigmagaussr=0.002;
        Double_t mu=0.3*2.5;//rate per pixel per gtu
        Double_t timelenght=10;//gtu
        TVector2 v1(xs1,ys1);
        TVector2 v2(xs2,ys2);
        TVector2 vmax(xs1+(xs2-xs1)*0.5,ys1+(ys2-ys1)*0.5);
        Double_t sigmadistv1v2=(v1-v2).Mod()*0.5;
        for(Int_t i=0;i<1000;i++){
                Double_t xb=x1+(rndm->Rndm())*(x2-x1);
                Double_t yb=y1+(rndm->Rndm())*(y2-y1);
                Int_t cb=(Int_t)(rndm->Poisson(mu)*timelenght);
                xbkg.push_back(xb);
                ybkg.push_back(yb);
                cbkg.push_back(cb);
                xbkgsign.push_back(xb);
                ybkgsign.push_back(yb);
                cbkgsign.push_back(cb);
                if(rndm->Rndm()<frac){
                        Double_t xsexact=xs1+(rndm->Rndm())*(xs2-xs1);
                        Double_t ysexact=(m*xsexact+q);
                        TVector2 vcurr(xsexact,ysexact);
                        Double_t distexact=(vcurr-vmax).Mod();
                        Double_t cs=10*TMath::Gaus(distexact,0,sigmadistv1v2,kFALSE);
                        for(Int_t ci=0;ci<cs;ci++){
                                Double_t r=rndm->Gaus(meangaussr,sigmagaussr);
                                Double_t phi=2*TMath::Pi()*rndm->Rndm();
                                Double_t xstep=r*cos(phi);
                                Double_t ystep=r*sin(phi);
                                TVector2 vstep(xstep,ystep);
                                vcurr=vcurr+vstep;
                                Double_t xs=vcurr.X();
                                Double_t ys=vcurr.Y();
                                xsign.push_back(xs);
                                ysign.push_back(ys);
                                csign.push_back(1);
                                xbkgsign.push_back(xs);
                                ybkgsign.push_back(ys);
                                cbkgsign.push_back(1);
                                hrndmbkgsign->Fill(xs,ys,1);
                        }
                }
                hrndmbkgsign->Fill(xb,yb,cb);
        }
        
        TCanvas *b=new TCanvas("b","b",1);
        b->Divide(2,1);
        b->cd(1);hrndmbkgsign->Draw();
        b->cd(2);hrndmbkgsign->Draw("LEGO");
        b->Write();
        delete b;
        delete hrndmbkgsign;
        
        vector<Double_t> x;
        vector<Double_t> y;
        vector<Int_t> c;
        if(dohough==1){
                x=xbkgsign;     y=ybkgsign;     c=cbkgsign;
        }else{
                if(dohough==2){
                        x=xsign;        y=ysign;        c=csign;
                }else{
                        x=xbkg; y=ybkg; c=cbkg;
                }
        }
        
        HoughFit *hf = new HoughFit(x, y, c, 1, ex, ey);
        Double_t mh = hf->GetSlope();
        Double_t qh = hf->GetOffset();
        Float_t cmaxtocmedio = hf->GetMaxToMedium();
        cout<<"\nfrom HoughFit:  mh="<<mh<<" (m_sim="<<m<<")"<<endl;
        cout<<"                qh="<<qh<<" (q_sim="<<q<<")"<<endl;
        cout<<"                cmaxtocmedio="<<cmaxtocmedio<<endl;
        f->Close();
        delete hf;
        
        
        cout<<"\nlook at the results in file --> "<<nameroot.c_str()<<"\n"<<endl;
        TBrowser* bro = new TBrowser();
        
}