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 | |
---|
22 | using namespace TMath; |
---|
23 | |
---|
24 | Double_t kHUGE = 1.e20; |
---|
25 | |
---|
26 | ClassImp(HoughFit) |
---|
27 | |
---|
28 | void HoughEstimator(Int_t &npar, Double_t *deriv, Double_t &f, Double_t *par, Int_t flag); |
---|
29 | |
---|
30 | //_________________________________________________________________________________________________________________________ |
---|
31 | HoughFit::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 | //_______________________________________________________________________________ |
---|
103 | HoughFit::~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 | //________________________________________________________________________________ |
---|
124 | void 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 | //______________________________________________________________________________________ |
---|
142 | void 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 | //______________________________________________________________________________________ |
---|
201 | void 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 | //_______________________________________________________________________________ |
---|
354 | void 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 | //_______________________________________________________________________________ |
---|
513 | void 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 | //______________________________________________________________________________________ |
---|
681 | void HoughFit::DoNumerical() { |
---|
682 | // |
---|
683 | // Numerical Hough Fit |
---|
684 | // to be implemented.. |
---|
685 | // |
---|
686 | } |
---|
687 | |
---|
688 | |
---|