source: JEM-EUSO/esaf_cc_at_lal/packages/reconstruction/modules/shower/fitting/src/HoughFit.cc @ 114

Last change on this file since 114 was 114, checked in by moretto, 11 years ago

actual version of ESAF at CCin2p3
File size: 25.6 KB
Line 
1// ESAF : Euso Simulation and Analysis Framework
2// $Id: HoughFit.cc 2949 2011-06-26 18:39:34Z mabl $
3// Elena Taddei created Mar, 29 2005
4         
5//////////////////////////////////////////////////////////////////////////////
6//                                                                          //
7//  HoughFit                                                                //
8//                                                                          //
9//                                                                          //
10//////////////////////////////////////////////////////////////////////////////
11
12#include <iostream>
13#include "MedianFit.hh"
14#include "HoughFit.hh"
15#include "TMath.h"
16#include <TH2.h>
17#include <TCanvas.h>
18#include <TGraph.h>
19#include "Config.hh"
20#include <TMinuit.h>
21
22using namespace TMath;
23
24Double_t kHUGE = 1.e20;
25
26ClassImp(HoughFit)
27
28void HoughEstimator(Int_t &npar, Double_t *deriv, Double_t &f, Double_t *par, Int_t flag);
29
30//_________________________________________________________________________________________________________________________
31HoughFit::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) {
32    //
33    // Constructor
34    // 
35   
36    string fKindHough = Conf()->GetStr("HoughFit.fKindHough");
37    fDoGraph = Conf()->GetBool("HoughFit.fDoGraph");
38    fWidthFractionSinusoidal = Conf()->GetNum("HoughFit.fWidthFractionSinusoidal");
39    fWidthFractionLinear = Conf()->GetNum("HoughFit.fWidthFractionLinear");
40   
41    fListObj = new TList;   
42    Clear();
43   
44    fEpsX = ex;
45    fEpsY = ey;
46    fId = id;
47    fPointsX = x;
48    fPointsY = y;
49    fCounts = c;
50    fNumPoints = fPointsX.size();
51   
52    if ( fNumPoints < 5 ) return;
53   
54    if ( fKindHough == "sinusoidal" ) {
55        DoSinusoidal(1);
56        fCMaxToCAverage = fCMaxToCAverageFluttTmp;
57        Float_t minfrac = Conf()->GetNum("HoughFit.fMinFractionSinusoidal");
58        //cout<<"fFluttBackgHSpace="<<fFluttBackgHSpace<<"\tfCMaxToCAverageFluttTmp="<<fCMaxToCAverageFluttTmp<<" fNumPointsSelected="<< fNumPointsSelected<<endl;
59        if( fCMaxToCAverage > minfrac*fFluttBackgHSpace && fDoSelection == 1 && fNumPointsSelected > 5 ){
60                //cout<<"There could be a signal, try to do selection..";
61                fPointsX = fXSel;
62                fPointsY = fYSel;
63                fCounts = fCSel;
64                fNumPoints = fPointsX.size();
65                Clear();
66                DoSinusoidal(2);
67                fCMaxToCAverageSel = fCMaxToCAverageFluttTmp;
68                //cout<<" .. now fFluttBackgHSpace="<<fFluttBackgHSpace<<"\tfCMaxToCAverageFluttTmp="<<fCMaxToCAverageFluttTmp<<endl;
69                //cout<<" .. fCMaxToCAverageSel/fFluttBackgHSpace="<<fCMaxToCAverageSel/fFluttBackgHSpace<<endl;
70                if( fCMaxToCAverageSel > (minfrac+1)*fFluttBackgHSpace ) fThereIsSignal = kTRUE;
71        }
72    } else if( fKindHough == "linear" ) {
73        DoLinear2(1);
74        fCMaxToCAverage = fCMaxToCAverageFluttTmp;
75        Float_t minfrac = Conf()->GetNum("HoughFit.fMinFractionLinear");
76        //cout<<"fFluttBackgHSpace="<<fFluttBackgHSpace<<"\tfCMaxToCAverageFluttTmp="<<fCMaxToCAverageFluttTmp<<" fNumPointsSelected="<< fNumPointsSelected<<endl;
77        if( fCMaxToCAverage > minfrac*fFluttBackgHSpace &&  fDoSelection == 1 && fNumPointsSelected > 5 ){
78                //cout<<"There could be a signal, try to do selection..";
79                fPointsX = fXSel;
80                fPointsY = fYSel;
81                fCounts = fCSel;
82                fNumPoints = fPointsX.size();
83                Clear();
84                DoLinear2(2);
85                fCMaxToCAverageSel = fCMaxToCAverageFluttTmp;
86                //cout<<" .. now fFluttBackgHSpace="<<fFluttBackgHSpace<<"\tfCMaxToCAverageFluttTmp="<<fCMaxToCAverageFluttTmp<<endl;
87                //if( fCMaxToCAverageSel > fCMaxToCAverage )       fThereIsSignal=1;
88                //cout<<" .. fCMaxToCAverageSel/fFluttBackgHSpace="<<fCMaxToCAverageSel/fFluttBackgHSpace<<endl;
89                if( fCMaxToCAverageSel > (minfrac+1)*fFluttBackgHSpace ) fThereIsSignal = kTRUE;
90        }
91    } else if( fKindHough == "numerical" ){
92        //Msg(EsafMsg::Panic) << "DoNumerical() not implemented yet.." << MsgDispatch;
93        //return;
94        DoNumerical();
95        fThereIsSignal = kTRUE;
96    } else {
97        Msg(EsafMsg::Panic) << "Wrong fKindHough value in config file." << MsgDispatch;
98    }
99
100}
101
102//_______________________________________________________________________________
103HoughFit::~HoughFit() {
104    //
105    //Destructor
106    //
107   
108    /*if (fListObj && fListObj->GetEntries()) {
109        fListObj->Delete();
110        fListObj->Clear();
111    }*/
112    fId.clear();
113    fPointsX.clear();
114    fPointsY.clear();
115    fPointsXprime.clear();
116    fPointsYprime.clear();
117    fCounts.clear();
118    fNumPointsSelected=0;
119    fNumPoints=0;
120    fIdSel.clear();
121}
122
123//________________________________________________________________________________
124void HoughFit::Clear() {
125    //
126    // Clear method
127    //
128       
129    fThereIsSignal = 0;
130    fXSel.clear();
131    fYSel.clear();
132    fCSel.clear();
133    fSlopeToPass = 0;
134    fOffsetToPass = 0;
135    fSlope = 0;
136    fOffset = 0;
137    fWidth = 0;
138    fAbsDev = 0;       
139}
140
141//______________________________________________________________________________________
142void HoughFit::PreProcess() {
143    //
144    // Preprocess for all kind of hough fit
145    //
146
147    fMinX = kHUGE;
148    fMinY = kHUGE;
149    fMaxX = -kHUGE;
150    fMaxY = -kHUGE;
151   
152    fPointsXprime.clear();
153    fPointsYprime.clear();
154       
155    for( Int_t i=0; i<fNumPoints; i++ ) {
156        Double_t x = fPointsX[i];
157        Double_t y = fPointsY[i];
158        if ( x > fMaxX ) fMaxX=x;
159        if ( y > fMaxY ) fMaxY=y;
160        if ( x < fMinX ) fMinX=x;
161        if ( y < fMinY ) fMinY=y;   
162    }
163    fX0 = (fMaxX+fMinX)/2;
164    fY0 = (fMaxY+fMinY)/2;
165   
166    Float_t tmpfMinX_pre = kHUGE;
167    Float_t tmpfMinY_pre = kHUGE;
168    Float_t tmpfMaxX_pre = -kHUGE;
169    Float_t tmpfMaxY_pre = -kHUGE;
170    vector<Double_t> fPointsXprime_pre;
171    vector<Double_t> fPointsYprime_pre;
172    for( Int_t i=0; i<fNumPoints; i++ ) {
173        Double_t x = fPointsX[i]-fX0;
174        Double_t y = fPointsY[i]-fY0;
175        fPointsXprime_pre.push_back(x);
176        fPointsYprime_pre.push_back(y);
177        if ( x > tmpfMaxX_pre ) tmpfMaxX_pre=x;
178        if ( y > tmpfMaxY_pre ) tmpfMaxY_pre=y;
179        if ( x < tmpfMinX_pre ) tmpfMinX_pre=x;
180        if ( y < tmpfMinY_pre ) tmpfMinY_pre=y;   
181    }
182   
183    tmpfMinX = kHUGE;
184    tmpfMinY = kHUGE;
185    tmpfMaxX = -kHUGE;
186    tmpfMaxY = -kHUGE;
187    fScale = tmpfMaxX_pre/tmpfMaxY_pre;
188    for( Int_t i=0; i<fNumPoints; i++ ) {
189        Double_t x = fPointsXprime_pre[i];
190        Double_t y = fPointsYprime_pre[i]*fScale;
191        fPointsXprime.push_back(x);
192        fPointsYprime.push_back(y);
193        if ( x > tmpfMaxX ) tmpfMaxX=x;
194        if ( y > tmpfMaxY ) tmpfMaxY=y;
195        if ( x < tmpfMinX ) tmpfMinX=x;
196        if ( y < tmpfMinY ) tmpfMinY=y;   
197    }
198}
199
200//______________________________________________________________________________________
201void HoughFit::DoSinusoidal(Int_t firsttimedo) {
202    //
203    // Recognition of linear tracks with sinusoidal Hough Transform in the space of
204    // parameters (ro,theta) P_i: (x_i,y_i) --> S_i: ro=x_i*cos(theta)+y_i*sin(theta)   
205    //
206 
207    PreProcess();
208   
209    Int_t binmax(0), nx(0), binymax(0), binxmax(0);
210    Double_t thetamaxsingle(0), rhomaxsingle(0);
211    Float_t rhomax = Sqrt(tmpfMaxX*tmpfMaxX + tmpfMaxY*tmpfMaxY);
212    Float_t rhomin = -rhomax;
213    Float_t epsilon = Sqrt(fEpsX*fEpsX + fEpsY*fEpsY);
214    Float_t deltaro = 2*epsilon;
215    Float_t thetasinglemin = 0;
216    Float_t thetasinglemax = Pi();
217    Int_t binx = (Int_t)((tmpfMaxX - tmpfMinX)*1000.);
218    Int_t biny = (Int_t)((tmpfMaxY - tmpfMinY)*1000.);
219    Int_t numbinro = (Int_t)((rhomax - rhomin)/deltaro);
220    Float_t deltatheta = 2*deltaro/(rhomax - rhomin);
221    Int_t numbintheta = (Int_t)((thetasinglemax - thetasinglemin)/deltatheta);
222    Int_t numsteptheta = numbintheta;
223    Float_t steptheta=(thetasinglemax - thetasinglemin)/numsteptheta;
224    TH2F *single = new TH2F("single","single;theta;ro",numbintheta,thetasinglemin,thetasinglemax,numbinro,rhomin,rhomax);
225    TH2F *single2 = new TH2F("single2","single2",binx,tmpfMinX,tmpfMaxX,biny,tmpfMinY,tmpfMaxY);
226    fListObj->Add(single);  fListObj->Add(single2);
227    for( Int_t i=0; i<fNumPoints; i++ ) {
228        Double_t x = fPointsXprime[i];
229        Double_t y = fPointsYprime[i];
230        Int_t z = fCounts[i];
231        single2->Fill(x,y,z);
232        for(Int_t j=0; j<numsteptheta; j++) {
233                Double_t thetasingle = j*steptheta;
234                Double_t rosingle = x*cos(thetasingle) + y*sin(thetasingle);
235                single->Fill(thetasingle,rosingle,z);
236        } 
237    }
238 
239    binmax = single->GetMaximumBin();
240    nx = (single->GetXaxis())->GetNbins()+2;
241   
242    Int_t plus1=0;
243    Int_t plus2=0;
244    while( (binmax % nx < 2 || binmax % nx > nx-2) && plus1<10 && plus2<10) {
245        if (binmax % nx < 2)    plus1++;
246        if (binmax % nx > nx-2) plus2++;
247        single->GetXaxis()->SetRange(plus1,numbintheta-plus2);
248        binmax = single->GetMaximumBin();
249        nx = (single->GetXaxis())->GetNbins()+2;
250    }
251   
252    binxmax = binmax % nx;
253    binymax = binmax / nx;
254    thetamaxsingle = (single->GetXaxis())->GetBinCenter(binxmax);
255    rhomaxsingle = (single->GetYaxis())->GetBinCenter(binymax);
256    Float_t countmax=single->GetBinContent(binxmax,binymax);
257   
258    /*TH2F *singlelt = new TH2F("singlelt","singlelt;theta;ro",numbintheta,thetasinglemin,thetasinglemax,numbinro,rhomin,rhomax);
259    TH2F *singlelt1 = new TH2F("singlelt1","singlelt1;theta;ro",numbintheta,thetasinglemin,thetasinglemax,numbinro,rhomin,rhomax);
260    TH2F *singlelt2 = new TH2F("singlelt2","singlelt2;theta;ro",numbintheta,thetasinglemin,thetasinglemax,numbinro,rhomin,rhomax);*/
261       
262    Float_t consideredbin = 0;
263    Float_t counttotal = 0;
264    for(Int_t ith=1; ith < numbintheta; ith++) {
265        for(Int_t iro=1; iro<numbinro; iro++) {
266                Int_t count = (Int_t)single->GetBinContent(ith,iro);
267                /*if(count>0.25*countmax )singlelt->Fill((single->GetXaxis())->GetBinCenter(ith),(single->GetYaxis())->GetBinCenter(iro),count);
268                if(count>0.5*countmax )singlelt1->Fill((single->GetXaxis())->GetBinCenter(ith),(single->GetYaxis())->GetBinCenter(iro),count);
269                if(count>0.75*countmax )singlelt2->Fill((single->GetXaxis())->GetBinCenter(ith),(single->GetYaxis())->GetBinCenter(iro),count);*/
270                if(count< 0.9*countmax ){
271                        counttotal = counttotal + count;
272                        consideredbin++;
273                }
274        }
275    }
276   
277    /*TCanvas *cHTlt=new TCanvas("cHTlt","cHTlt",1);cHTlt->Divide(3,1);
278    cHTlt->cd(1);   singlelt->Draw();
279    cHTlt->cd(2);   singlelt1->Draw();
280    cHTlt->cd(3);   singlelt2->Draw();
281    cHTlt->Write();
282    delete cHTlt;       delete singlelt;    delete singlelt1;    delete singlelt2;*/
283   
284    Float_t countaverage = counttotal/(consideredbin);
285    fFluttBackgHSpace = Sqrt(countaverage);
286
287    fCMaxToCAverageTmp = (Float_t)(countmax/countaverage);
288    fCMaxToCAverageFluttTmp = (Float_t)(countmax/fFluttBackgHSpace);
289   
290    TH1D *singleproj; //=new TH1D("singleproj","singleproj",numbinro,rhomin,rhomax);
291    singleproj=single->ProjectionY("singleproj",binxmax-5,binxmax+5);
292    fListObj->Add(singleproj);
293       
294    fWidth = singleproj->GetRMS();
295    fSlope = -1/Tan(thetamaxsingle);
296    fOffset = rhomaxsingle/Sin(thetamaxsingle);
297    fSlopeToPass = fSlope/fScale;
298    fOffsetToPass = fY0+((fOffset-fSlope*fX0)/fScale);
299       
300    Double_t xg[fNumPoints],yg[fNumPoints];
301    fAbsDev = 0;
302    Double_t markercenterx[1] = {thetamaxsingle};
303    Double_t markercentery[1] = {rhomaxsingle};
304    TGraph *gmarkercenter = new TGraph(1,markercenterx,markercentery);
305    fListObj->Add(gmarkercenter);
306    gmarkercenter->SetMarkerColor(2); gmarkercenter->SetMarkerStyle(29); gmarkercenter->SetMarkerSize(1.6);
307       
308    for( Int_t i=0; i<fNumPoints; i++ ) {
309        xg[i] = fPointsXprime[i];
310        yg[i] = fSlope * fPointsXprime[i] + fOffset;
311        Float_t singlefAbsDev = Abs(yg[i] - fPointsYprime[i]);
312        fAbsDev += singlefAbsDev;
313    }
314        fAbsDev /= fNumPoints;
315       
316    if( fDoGraph ) {
317        TGraph *g = new TGraph(fNumPoints,xg,yg); g->SetLineColor(8);
318        fListObj->Add(g);
319        string namecHT;
320        if( firsttimedo==2 ) {
321                namecHT="cHTsel";
322        }else{
323                namecHT="cHT";
324        }
325        TCanvas *cHT = new TCanvas(namecHT.c_str(),namecHT.c_str(),1);cHT->Divide(3,1);
326        fListObj->Add(cHT);
327        cHT->cd(1);    single2->SetMarkerStyle(6); single2->Draw();g->Draw("PLsame");
328        cHT->cd(2);    single->Draw();gmarkercenter->Draw("Psame");
329        cHT->cd(3);    single->Draw("LEGO");
330            cHT->Write();
331        SafeDelete(cHT); SafeDelete(g);
332    }
333    SafeDelete(single2); SafeDelete(single); SafeDelete(gmarkercenter); SafeDelete(singleproj);
334   
335    fNumPointsSelected=0;
336    if( firsttimedo==1 && fWidth>1e-4 ){
337        for( Int_t i=0; i<fNumPoints; i++ ) {
338            Double_t x = fPointsXprime[i];
339            Double_t y = fPointsYprime[i];
340            Double_t rosingleprova = x*Cos(thetamaxsingle) + y*Sin(thetamaxsingle);
341            if(Abs(rosingleprova-rhomaxsingle)<fWidth*fWidthFractionSinusoidal) {
342                fXSel.push_back(x+fX0);
343                fYSel.push_back((y/fScale)+fY0);
344                fCSel.push_back(fCounts[i]);
345                if( fId.size()>0 ) fIdSel.push_back(fId[i]);
346                fNumPointsSelected++;
347            } 
348        }
349    }
350       
351}
352
353//_______________________________________________________________________________
354void HoughFit::DoLinear(Int_t firsttimedo) {
355    //
356    // Recognition of linear tracks with linear Hough Transform in the space
357    // of parameters (q,m) P_i: (x_i,y_i) --> L_i: q=-m*x_i+y_i         
358    //
359
360    PreProcess();
361
362    Int_t binmax(0),nx(0),binymax(0),binxmax(0);
363       
364    Double_t mhoughmax = tmpfMaxX;
365    Double_t qhoughmax = tmpfMaxY;
366    Double_t mmax = 2.;//to establish physically from EAS velocity, still to improve
367    Double_t mmin = -mmax;
368    Double_t qmax = mmax*mhoughmax + qhoughmax;
369    Double_t qmin = -qmax;
370    Float_t deltaq = Sqrt(mmax*mmax * fEpsX*fEpsX + fEpsY*fEpsY);
371    Int_t numbinq = (Int_t)((qmax-qmin)/deltaq);
372    Float_t deltam = deltaq/mhoughmax;
373    Int_t numbinm = (Int_t)((mmax-mmin)/deltam);
374    Int_t numstepm = numbinm;
375    Float_t stepm = deltam;
376       
377    Int_t binx=(Int_t)((tmpfMaxX-tmpfMinX)*1000.);//not important binning, is just to draw
378    Int_t biny=(Int_t)((tmpfMaxY-tmpfMinY)*1000.);//not important binning, is just to draw
379       
380    TH2F *single = new TH2F("single","single;m;q",numbinm,mmin,mmax,numbinq,qmin,qmax);
381    TH2F *single2 = new TH2F("single2","single2",binx,tmpfMinX,tmpfMaxX,biny,tmpfMinY,tmpfMaxY);
382    fListObj->Add(single);
383    fListObj->Add(single2);
384    for( Int_t i=0; i<fNumPoints; i++ ) {
385        Double_t x = fPointsXprime[i];
386        Double_t y = fPointsYprime[i];
387        Int_t z = fCounts[i];
388        single2->Fill(x,y,z);
389        for(Int_t j=0;j<numstepm;j++){
390                Double_t msingle = mmin+j*stepm;
391                Double_t qsingle = -msingle*x+y;
392                single->Fill(msingle,qsingle,z);
393        } 
394    }
395 
396    binmax = single->GetMaximumBin();
397    nx = (single->GetXaxis())->GetNbins()+2;
398   
399    Int_t plus1=0;
400    Int_t plus2=0;
401    while( (binmax % nx < 2 || binmax % nx > nx-2) && plus1<10 && plus2<10) {
402        if (binmax % nx < 2)    plus1++;
403        if (binmax % nx > nx-2) plus2++;
404        single->GetXaxis()->SetRange(plus1,numbinm-plus2);
405        binmax = single->GetMaximumBin();
406        nx = (single->GetXaxis())->GetNbins()+2;
407    }
408    binxmax = binmax % nx;
409    binymax = binmax / nx;
410    Float_t mmaxsingle = (single->GetXaxis())->GetBinCenter(binxmax);
411    Float_t qmaxsingle = (single->GetYaxis())->GetBinCenter(binymax);
412    Float_t countmax = single->GetBinContent(binxmax,binymax);
413   
414    fSlope = mmaxsingle;
415    fOffset = qmaxsingle;
416    fSlopeToPass = fSlope/fScale;
417    fOffsetToPass = fY0+((fOffset-fSlope*fX0)/fScale);
418       
419    Double_t markercenterx[1] = {mmaxsingle};
420    Double_t markercentery[1] = {qmaxsingle};
421    TGraph *gmarkercenter = new TGraph(1,markercenterx,markercentery);
422    fListObj->Add(gmarkercenter);
423    gmarkercenter->SetMarkerColor(2); gmarkercenter->SetMarkerStyle(29); gmarkercenter->SetMarkerSize(1.6);
424       
425    Double_t xg[fNumPoints],yg[fNumPoints];
426       
427    fAbsDev = 0;
428    for( Int_t i=0; i<fNumPoints; i++ ) {
429        xg[i] = fPointsXprime[i];
430        yg[i] = fSlope*fPointsXprime[i] + fOffset;
431        Float_t singlefAbsDev = fabs(yg[i] - fPointsYprime[i]);
432        fAbsDev += singlefAbsDev;
433    }
434        fAbsDev /= fNumPoints;
435       
436       
437    /*TH2F *singlelt = new TH2F("singlelt","singlelt;m;q",numbinm,mmin,mmax,numbinq,qmin,qmax);
438    TH2F *singlelt1 = new TH2F("singlelt1","singlelt1;m;q",numbinm,mmin,mmax,numbinq,qmin,qmax);
439    TH2F *singlelt2 = new TH2F("singlelt2","singlelt2;m;q",numbinm,mmin,mmax,numbinq,qmin,qmax);*/
440       
441    Float_t consideredbin = 0;
442    Float_t counttotal = 0;
443    for(Int_t im=1;im<numbinm; im++) {
444        for(Int_t iq=1; iq<numbinq; iq++) {
445            Int_t count=(Int_t)(single->GetBinContent(im,iq));
446            /*if(count>0.25*countmax )singlelt->Fill((single->GetXaxis())->GetBinCenter(im),(single->GetYaxis())->GetBinCenter(iq),count);
447            if(count>0.5*countmax )singlelt1->Fill((single->GetXaxis())->GetBinCenter(im),(single->GetYaxis())->GetBinCenter(iq),count);
448            if(count>0.75*countmax )singlelt2->Fill((single->GetXaxis())->GetBinCenter(im),(single->GetYaxis())->GetBinCenter(iq),count);*/
449            if(count<0.9*countmax ){
450                counttotal += count;
451                consideredbin++;
452            }
453        }
454    }
455    Float_t countaverage = counttotal/(consideredbin);
456    fFluttBackgHSpace = Sqrt(countaverage);
457
458    fCMaxToCAverageTmp = (Float_t)(countmax/countaverage);
459    fCMaxToCAverageFluttTmp = (Float_t)(countmax/fFluttBackgHSpace);
460   
461    TH1D *singleproj;//=new TH1D("singleproj","singleproj",numbinq,qmin,qmax);
462    singleproj=single->ProjectionY("singleproj",binxmax-5,binxmax+5);
463    fWidth = singleproj->GetRMS();
464   
465   
466   /* TCanvas *cHTlt=new TCanvas("cHTlt","cHTlt",1);cHTlt->Divide(2,2);
467    cHTlt->cd(1);   singlelt->Draw();
468    cHTlt->cd(2);   singlelt1->Draw();
469    cHTlt->cd(3);   singlelt2->Draw();
470    cHTlt->cd(4);   singleproj->Draw();
471    cHTlt->Write();
472    delete cHTlt; delete singlelt; delete singlelt1; delete singlelt2; delete singleproj;*/
473   
474    if(fDoGraph){
475        TGraph *g = new TGraph(fNumPoints,xg,yg); g->SetLineColor(8);
476        fListObj->Add(g);
477            string namecHT;
478        if (firsttimedo==2){
479                namecHT = "cHTsecond";
480        }else{
481                namecHT = "cHTfirst";
482        }
483        TCanvas *cHT=new TCanvas(namecHT.c_str(),namecHT.c_str(),1);cHT->Divide(3,1);
484        fListObj->Add(cHT);
485        cHT->cd(1);    single2->SetMarkerStyle(6); single2->Draw();g->Draw("PLsame");
486        cHT->cd(2);    single->Draw();gmarkercenter->Draw("Psame");
487        cHT->cd(3);    single->Draw("LEGO");
488            cHT->Write();
489            SafeDelete(cHT);
490        SafeDelete(g);
491    }
492    SafeDelete(single2); SafeDelete(single); SafeDelete(gmarkercenter); SafeDelete(singleproj); 
493   
494    fNumPointsSelected=0;
495    if(firsttimedo==1 && fWidth>1e-4){
496        for( Int_t i=0; i<fNumPoints; i++ ) {
497            Double_t x = fPointsXprime[i];
498            Double_t y = fPointsYprime[i];
499            Double_t qsingleprova = -mmaxsingle*x+y;
500            if(Abs(qsingleprova - qmaxsingle) < fWidth*fWidthFractionLinear) {
501                fXSel.push_back(x+fX0);
502                fYSel.push_back((y/fScale)+fY0);
503                fCSel.push_back(fCounts[i]);
504                if( fId.size()>0 ) fIdSel.push_back(fId[i]);
505                    fNumPointsSelected++;
506            } 
507        }
508        //cout<<"\nfNumPoints="<<fNumPoints<<"\tfNumPointsSelected="<<fNumPointsSelected<<"\tfWidth="<<fWidth<<endl;
509    }
510}
511
512//_______________________________________________________________________________
513void HoughFit::DoLinear2(Int_t firsttimedo) {
514    //
515    // Recognition of linear tracks with linear Hough Transform in the space
516    // of parameters (q,m) P_i: (x_i,y_i) --> L_i: q=-m*x_i+y_i         
517    //
518
519    PreProcess();
520
521    Int_t binmax(0),nx(0),binymax(0),binxmax(0);
522       
523    Double_t mhoughmax = tmpfMaxX;
524    Double_t qhoughmax = tmpfMaxY;
525    Double_t mmax = 2.;//to establish phisically from EAS velocity, still to improve
526    Double_t mmin = -mmax;
527    Double_t qmax = mmax*mhoughmax + qhoughmax;
528    Double_t qmin = -qmax;
529    Float_t deltaq = Sqrt(mmax*mmax * fEpsX*fEpsX + fEpsY*fEpsY);
530    Int_t numbinq = (Int_t)((qmax-qmin)/deltaq);
531    Float_t deltam = deltaq/mhoughmax;
532    Int_t numbinm = (Int_t)((mmax-mmin)/deltam);
533    Int_t numstepm = numbinm;
534    Float_t stepm = deltam;
535       
536    Int_t binx=(Int_t)((tmpfMaxX-tmpfMinX)*1000.);//not important binning, is just to draw
537    Int_t biny=(Int_t)((tmpfMaxY-tmpfMinY)*1000.);//not important binning, is just to draw
538       
539    TH2F *single = new TH2F("single","single;m;q",numbinm,mmin,mmax,numbinq,qmin,qmax);
540    TH2F *single2 = new TH2F("single2","single2",binx,tmpfMinX,tmpfMaxX,biny,tmpfMinY,tmpfMaxY);
541    fListObj->Add(single);
542    fListObj->Add(single2);
543    for( Int_t i=0; i<fNumPoints; i++ ) {
544        Double_t x = fPointsXprime[i];
545        Double_t y = fPointsYprime[i];
546        Int_t z = fCounts[i];
547        single2->Fill(x,y,z);
548        for(Int_t j=0;j<numstepm;j++){
549                Double_t msingle = mmin+j*stepm;
550                Double_t qsingle = -msingle*x+y;
551                single->Fill(msingle,qsingle,z);
552        } 
553    }
554 
555    binmax = single->GetMaximumBin();
556    nx = (single->GetXaxis())->GetNbins()+2;
557   
558    Int_t plus1=0;
559    Int_t plus2=0;
560    while( (binmax % nx < 2 || binmax % nx > nx-2) && plus1<10 && plus2<10) {
561        if (binmax % nx < 2)    plus1++;
562        if (binmax % nx > nx-2) plus2++;
563        single->GetXaxis()->SetRange(plus1,numbinm-plus2);
564        binmax = single->GetMaximumBin();
565        nx = (single->GetXaxis())->GetNbins()+2;
566    }
567    binxmax = binmax % nx;
568    binymax = binmax / nx;
569    Float_t mmaxsingle = (single->GetXaxis())->GetBinCenter(binxmax);
570    Float_t qmaxsingle = (single->GetYaxis())->GetBinCenter(binymax);
571    Float_t countmax = single->GetBinContent(binxmax,binymax);
572   
573    fSlope = mmaxsingle;
574    fOffset = qmaxsingle;
575    fSlopeToPass = fSlope/fScale;
576    fOffsetToPass = fY0+((fOffset-fSlope*fX0)/fScale);
577       
578    Double_t markercenterx[1] = {mmaxsingle};
579    Double_t markercentery[1] = {qmaxsingle};
580    TGraph *gmarkercenter = new TGraph(1,markercenterx,markercentery);
581    fListObj->Add(gmarkercenter);
582    gmarkercenter->SetMarkerColor(2); gmarkercenter->SetMarkerStyle(29); gmarkercenter->SetMarkerSize(1.6);
583       
584    Double_t xg[fNumPoints],yg[fNumPoints];
585       
586    fAbsDev = 0;
587    for( Int_t i=0; i<fNumPoints; i++ ) {
588        xg[i] = fPointsXprime[i];
589        yg[i] = fSlope*fPointsXprime[i] + fOffset;
590        Float_t singlefAbsDev = fabs(yg[i] - fPointsYprime[i]);
591        fAbsDev += singlefAbsDev;
592    }
593        fAbsDev /= fNumPoints;
594       
595       
596    /*TH2F *singlelt = new TH2F("singlelt","singlelt;m;q",numbinm,mmin,mmax,numbinq,qmin,qmax);
597    TH2F *singlelt1 = new TH2F("singlelt1","singlelt1;m;q",numbinm,mmin,mmax,numbinq,qmin,qmax);
598    TH2F *singlelt2 = new TH2F("singlelt2","singlelt2;m;q",numbinm,mmin,mmax,numbinq,qmin,qmax);*/
599       
600    Float_t consideredbin = 0;
601    Float_t counttotal = 0;
602    for(Int_t im=1;im<numbinm; im++) {
603        for(Int_t iq=1; iq<numbinq; iq++) {
604            Int_t count=(Int_t)(single->GetBinContent(im,iq));
605            /*if(count>0.25*countmax )singlelt->Fill((single->GetXaxis())->GetBinCenter(im),(single->GetYaxis())->GetBinCenter(iq),count);
606            if(count>0.5*countmax )singlelt1->Fill((single->GetXaxis())->GetBinCenter(im),(single->GetYaxis())->GetBinCenter(iq),count);
607            if(count>0.75*countmax )singlelt2->Fill((single->GetXaxis())->GetBinCenter(im),(single->GetYaxis())->GetBinCenter(iq),count);*/
608            if(count<0.9*countmax ){
609                counttotal += count;
610                consideredbin++;
611            }
612        }
613    }
614    Float_t countaverage = counttotal/(consideredbin);
615    fFluttBackgHSpace = Sqrt(countaverage);
616
617    fCMaxToCAverageTmp = (Float_t)(countmax/countaverage);
618    fCMaxToCAverageFluttTmp = (Float_t)(countmax/fFluttBackgHSpace);
619   
620    TH1D *singleproj;//=new TH1D("singleproj","singleproj",numbinq,qmin,qmax);
621    singleproj=single->ProjectionY("singleproj",binxmax-5,binxmax+5);
622    fWidth = singleproj->GetRMS();
623   
624   
625   /* TCanvas *cHTlt=new TCanvas("cHTlt","cHTlt",1);cHTlt->Divide(2,2);
626    cHTlt->cd(1);   singlelt->Draw();
627    cHTlt->cd(2);   singlelt1->Draw();
628    cHTlt->cd(3);   singlelt2->Draw();
629    cHTlt->cd(4);   singleproj->Draw();
630    cHTlt->Write();
631    delete cHTlt; delete singlelt; delete singlelt1; delete singlelt2; delete singleproj;*/
632   
633    if(fDoGraph){
634        TGraph *g = new TGraph(fNumPoints,xg,yg); g->SetLineColor(8);
635        fListObj->Add(g);
636            string namecHT;
637        if (firsttimedo==2){
638                namecHT = "cHTsecond";
639        }else{
640                namecHT = "cHTfirst";
641        }
642        TCanvas *cHT=new TCanvas(namecHT.c_str(),namecHT.c_str(),1);cHT->Divide(3,1);
643        fListObj->Add(cHT);
644        cHT->cd(1);    single2->SetMarkerStyle(6); single2->Draw();g->Draw("PLsame");
645        cHT->cd(2);    single->Draw();gmarkercenter->Draw("Psame");
646        cHT->cd(3);    single->Draw("LEGO");
647            cHT->Write();
648            SafeDelete(cHT);
649        SafeDelete(g);
650    }
651    SafeDelete(single2); SafeDelete(single); SafeDelete(gmarkercenter); SafeDelete(singleproj); 
652   
653    fNumPointsSelected=0;
654    if(firsttimedo==1 && fWidth>1e-4){
655        // calculate the RMS of the distance from points to the rect
656        Double_t *distance = new Double_t[(const Int_t)fNumPoints];
657            for( Int_t i=0; i<fNumPoints; i++ ) {
658            distance[i] = Abs(fPointsXprime[i]*mmaxsingle-fPointsYprime[i]+qmaxsingle);
659        }
660        Double_t distrms = RMS(fNumPoints,distance);
661        // select the points within 1 RMS from rect
662        for( Int_t i(0); i<fNumPoints; i++ ) {
663            Double_t x = fPointsXprime[i];
664            Double_t y = fPointsYprime[i];
665            //Double_t qsingleprova = -mmaxsingle*x+y;
666            if( distance[i] < fWidthFractionLinear*distrms) {
667                fXSel.push_back(x+fX0);
668                fYSel.push_back((y/fScale)+fY0);
669                fCSel.push_back(fCounts[i]);
670                if( fId.size()>0 ) fIdSel.push_back(fId[i]);
671                    fNumPointsSelected++;
672            }
673
674       
675        //cout<<"\nfNumPoints="<<fNumPoints<<"\tfNumPointsSelected="<<fNumPointsSelected<<"\tfWidth="<<fWidth<<endl;
676        }
677    }
678}
679
680//______________________________________________________________________________________
681void HoughFit::DoNumerical() {
682    //
683    // Numerical Hough Fit
684    // to be implemented..     
685    //
686}
687
688
Note: See TracBrowser for help on using the repository browser.