Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

RobustModule.cc @ 117

Last change on this file since 117 was 117, checked in by moretto, 11 years ago
ESAF version compilable on mac OS
File size: 20.6 KB

Line
1	// $Id$
2	// Author: Svetlana Biktemerova Okt, 20 2011
3
4	/*****************************************************************************
5	* ESAF: Euso Simulation and Analysis Framework *
6	* *
7	* Id: RobustModule *
8	* Package: <packagename> *
9	* Coordinator: <coordinator> *
10	* *
11	*****************************************************************************/
12
13	//_____________________________________________________________________________
14	//
15	// RobustModule
16	//
17	// <extensive class description>
18	//
19	// Config file parameters
20	// ======================
21	//
22	// <parameter name>: <parameter description>
23	// -Valid options: <available options>
24	//
25	#include <TVirtualFitter.h>
26	#include "RobustModule.hh"
27	#include "RecoEvent.hh"
28	#include "EEvent.hh"
29	#include "EGeometry.hh"
30	#include "EShower.hh"
31	#include "ETruth.hh"
32	#include "RecoPixelData.hh"
33	#include "ERunParameters.hh"
34	#include "EusoCluster.hh"
35	#include "RecoGlobalData.hh"
36	#include <TVector3.h>
37	#include <TGraph.h>
38	#include <TH2F.h>
39	#include <TF1.h>
40	#include <TFile.h>
41	#include <TStyle.h>
42	#include <TFitResultPtr.h>
43	#include <TFitResult.h>
44	#include <TH1D.h>
45	#include <TH2D.h>
46	#include <TGraph2DErrors.h>
47	#include <TCanvas.h>
48	#include <TH2F.h>
49	#include <TStyle.h>
50	#include <TGraphErrors.h>
51	#include <TTree.h>
52
53	#include "EConst.hh"
54
55	#include <vector>
56	using std::vector;
57
58	// #include <globals.hh>
59	using namespace TMath;
60	using namespace std;
61	ClassImp(RobustModule);
62	//ClassImp(RobustModule::TrackInformation);
63
64	//#define RobustDebug
65	#define minTOLERANCE 1.e-9
66	RobustModule *RobustModule::fMe = NULL;
67
68	//_____________________________________________________________________________
69	RobustModule::RobustModule() : RecoModule("RobustModule"), fWeights(0) {
70	// ctor
71
72	if ( fMe == NULL )
73	fMe = this;
74
75	}
76
77	//_____________________________________________________________________________
78	RobustModule::~RobustModule() {
79	// dtor
80	if(fWeights) delete [] fWeights;
81
82	}
83
84	//_____________________________________________________________________________
85	RobustModule* RobustModule::Get() {
86	return fMe;
87	}
88
89	//_____________________________________________________________________________
90	Bool_t RobustModule::Init() {
91	return kTRUE;
92	}
93
94	//_____________________________________________________________________________
95	Bool_t RobustModule::PreProcess() {
96	fx1=fy1=1e+20;
97	fx2=fy2=-1e+20;
98	if ( fWeights ) delete fWeights;
99	fWeights=0;
100
101	Msg(EsafMsg::Info) << "RobustModule" << MsgDispatch;
102	return kTRUE;
103	}
104
105	//_____________________________________________________________________________
106	Bool_t RobustModule::Process(RecoEvent *ev) {
107
108	fEvent = ev;
109	fRunpars = fEvent->GetHeader().GetRunPars();
110	if(!FindTrackOnFocalPlane())return kFALSE;
111	SaveGlobalData();
112	return kTRUE;
113
114	}
115
116	//_____________________________________________________________________________
117	Bool_t RobustModule::PostProcess() {
118	fEvent = (RecoEvent*)0;
119	return kTRUE;
120	}
121
122	//_____________________________________________________________________________
123	Bool_t RobustModule::Done() {
124	Msg(EsafMsg::Info) << "RobustModule done. " << MsgDispatch;
125	return kTRUE;
126	}
127
128	//_____________________________________________________________________________
129	void RobustModule::UserMemoryClean() {
130	//
131	// memory clean
132	//
133	}
134
135	//_____________________________________________________________________________
136	Bool_t RobustModule::SaveRootData(RecoRootEvent *fRecoRootEvent) {
137	return kTRUE;
138	}
139
140	//_____________________________________________________________________________
141	Bool_t RobustModule::SaveGlobalData() {
142	RecoGlobalData* data = fEvent->FindCreateGlobalData("PatternRecognition");
143	data->CleanMaps();
144	data->SetIssuer(GetName());
145	data->Add("SelectedPixels",&sigpix);
146	data->Add("TrackInformation",&ftrack_info);
147	return kTRUE;
148	}
149	//_____________________________________________________________________________
150
151	Bool_t RobustModule::FindTrackOnFocalPlane(double *coeff) {
152
153	//get triggered macrocells (pdms)
154	map< int, map< int, vector<RecoPixelData*> > > cell = fEvent->GetCellGtuRecoPixelData();
155	map< int, map< int, vector<RecoPixelData*> > >::iterator itcell;
156
157	if(cell.size()==0){
158	Msg(EsafMsg::Info) << "There are not triggered pdms. " << MsgDispatch;
159	return kFALSE;
160	}
161	// itpixels->first - pixels id;
162	// itpixels->second - TVector3 v, where v.X() and v.Y() - coordinates of pixel center, v.Z()- the number of counts in pixel;
163
164	map <int, TVector3 > pixels;
165	map <int, TVector3 >::iterator itpixels;
166	int average_counts=0;
167	// fill pixels
168	int fgtu1, fgtu2;
169	fgtu1=fgtu2=-10;
170	for(itcell=cell.begin(); itcell!=cell.end(); itcell++){
171
172	map< int, vector<RecoPixelData*> >& recopixels = itcell->second;
173	map< int, vector<RecoPixelData*> >::iterator it;
174
175	for (it= recopixels.begin();it!=recopixels.end();it++) {
176	vector<RecoPixelData*>& vpix = it->second;
177	for(unsigned int i=0;i<vpix.size();i++){
178	RecoPixelData* pix = vpix[i];
179	TVector3 center = fRunpars->PixelCenter(pix->GetPixelId());
180	average_counts+=pix->GetCounts();
181	int gtu = pix->GetGtu();
182	if(gtu<fgtu1) fgtu1=gtu;
183	if(gtu>fgtu2) fgtu2=gtu;
184	if(center.X()<fx1) fx1=center.X();
185	if(center.X()>fx2) fx2=center.X();
186	if(center.Y()<fy1) fy1=center.Y();
187	if(center.Y()>fy2) fy2=center.Y();
188
189	itpixels=pixels.find(pix->GetPixelId());
190	if(itpixels != pixels.end()){
191	TVector3& v = itpixels->second;
192	v.SetZ(v.Z()+pix->GetCounts());
193	pixels[itpixels->first]=v;
194
195	}
196	else {
197	TVector3 v;
198	v.SetXYZ(center.X(),center.Y(),pix->GetCounts());
199	pixels[pix->GetPixelId()]=v;
200	}
201	}
202	}
203	}
204
205	//cpixels=pixels;
206
207	//temporary method of background cut
208	//average_counts - average number of counts in pixel
209	int npix = pixels.size();
210	average_counts/=npix;
211	for ( itpixels=pixels.begin(); itpixels!=pixels.end(); itpixels++) {
212	TVector3& pixel_info= itpixels->second;
213	//pixel_info.SetZ(pixel_info.Z()-average_counts);
214	if(int(pixel_info.Z())-average_counts<=0) pixels.erase (itpixels);
215
216	}
217
218
219
220	#ifdef RobustDebug
221	//number of counts with background
222	TH1D* hcounts1 = new TH1D("counts1","Number of counts with bg",300,0.,300.);
223	TH1D* weighted_counts1 = new TH1D("weighted_counts1","Number of counts with bg",300,0.,300.);
224	//number of counts minus average number of counts
225	TH1D* hcounts2 = new TH1D("counts2","Number of counts - average_counts",300,0.,300.);
226	TH1D* weighted_counts2 = new TH1D("weighted_counts2","Number of counts - average_counts",300,0.,300.);
227	hcounts1->GetXaxis()->SetTitle("N counts");
228	weighted_counts1->GetXaxis()->SetTitle("N counts");
229	weighted_counts1->GetYaxis()->SetTitle("N^{2} counts");
230	hcounts2->GetXaxis()->SetTitle("(N - N_{average})");
231	weighted_counts2->GetXaxis()->SetTitle("(N - N_{average})");
232	weighted_counts2->GetYaxis()->SetTitle("(N - N_{average})^{2}");
233
234
235	for ( itpixels=pixels.begin(); itpixels!=pixels.end(); itpixels++) {
236	TVector3& pixel_info= itpixels->second;
237	hcounts1->Fill(pixel_info.Z());
238	weighted_counts1->Fill(pixel_info.Z(),pow(pixel_info.Z(),2.));
239	double nc = pixel_info.Z()-average_counts;
240	if(nc>0){
241	hcounts2->Fill(nc);
242	weighted_counts2->Fill(nc,pow(nc,2.));
243	}
244	}
245	TCanvas* c = new TCanvas("c2");
246	c->Print("counts.eps[");
247	hcounts1->Draw();
248	c->Print("counts.eps");
249	weighted_counts1->Draw();
250	c->Print("counts.eps");
251	hcounts2->Draw();
252	c->Print("counts.eps");
253	weighted_counts2->Draw();
254	c->Print("counts.eps");
255	c->Print("counts.eps]");
256
257	delete hcounts1;
258	delete weighted_counts1;
259	delete hcounts2;
260	delete weighted_counts2;
261	#endif
262
263	if(fWeights) delete [] fWeights;
264	fWeights=new double[npix];
265
266	for(int i(0);i<npix;i++)fWeights[i]=1.;
267
268	double ax(0.), bx(0.),a,b;
269	double min_disp(100.), disp, sumw;
270	int counts=0;
271	int i=0;
272
273	// Ax and Bx parameters of line
274	Coefficients(ax,bx,pixels);
275
276	do {
277	disp = 0.;
278	a = ax; b = bx;
279	sumw=0.;
280
281	// Calculate dispersion
282
283	i=0;
284	for ( itpixels=pixels.begin(); itpixels!=pixels.end(); itpixels++) {
285	TVector3& pixel_info= itpixels->second;
286	double dist = (pixel_info.Y() - apixel_info.X() - b)/Sqrt(1.+aa);
287	disp += distdistfWeights[i];
288	sumw += fWeights[i];
289	i++;
290	}
291
292	disp/=sumw;
293
294	if (disp < min_disp) disp = min_disp;
295
296	// Calculate fWeights
297	i=0;
298	for ( itpixels=pixels.begin(); itpixels!=pixels.end(); itpixels++) {
299	TVector3& pixel_info = itpixels->second;
300	double dist = (pixel_info.Y() - apixel_info.X() - b)/Sqrt(1+aa);
301	double n = (pixel_info.Z()-average_counts);
302	fWeights[i] = pow(n<0.?0.:n,2.)Exp(-distdist/disp);
303	i++;
304	}
305
306	// Ax and Bx parameters
307	Coefficients(ax,bx,pixels);
308	counts++;
309
310	}
311	while (((a-ax)(a-ax) + (b-bx)(b-bx)) > 0.00005 && counts < 1000);
312
313	TF1* f = &ftrack_info.fLine;
314	f->SetRange(fx1, fx2);
315	f->FixParameter(0, ax);
316	f->FixParameter(1, bx);
317
318	#ifdef RobustDebug
319	int csize = cell.size();
320	TH2D* hcell_c = new TH2D(Form("Cell %d, counts",csize),Form("Number of pdm = %d, counts",csize),100,fx1,fx2,100,fy1,fy2);
321	TH2D* hcell_w = new TH2D(Form("Cell %d, weights",csize),Form("Number of pdm = %d, weights",csize),100,fx1,fx2,100,fy1,fy2);
322	hcell_c->SetStats(0);
323	hcell_w->SetStats(0);
324
325	i=0;
326	for ( itpixels=pixels.begin(); itpixels!=pixels.end(); itpixels++) {
327	TVector3& pixel_info = itpixels->second;
328	double n = pixel_info.Z();
329	hcell_c->SetBinContent(hcell_c->GetXaxis()->FindBin(pixel_info.X()),hcell_c->GetYaxis()->FindBin(pixel_info.Y()),n);
330	hcell_w->SetBinContent(hcell_c->GetXaxis()->FindBin(pixel_info.X()),hcell_c->GetYaxis()->FindBin(pixel_info.Y()),fWeights[i]);
331	i++;
332	}
333
334	gStyle->SetPalette(1);
335	c->Print("outline.eps[");
336	hcell_c->Draw("colorz");
337	f->Draw("same");
338	c->Print("outline.eps");
339	hcell_w->Draw("colorz");
340	f->Draw("same");
341	c->Print("outline.eps");
342	c->Print("outline.eps]");
343
344	delete hcell_c;
345	delete hcell_w;
346
347	delete c;
348	#endif
349
350	if(!GetSignalPixels(a,b,pixels)) return kFALSE;
351	return kTRUE;
352	}
353
354	//_____________________________________________________________________________
355	Bool_t RobustModule::GetSignalPixels(double a, double b, map <int, TVector3 >& pixels){
356
357	double px[2];
358	double py[2];
359	double x[4];
360	x[0]=fx1;
361	x[1]=ftrack_info.fLine.Eval(fx1);
362	x[2]=fx2;
363	x[3]=ftrack_info.fLine.Eval(fx2);
364	//find intersection between line and border of triggered pdms
365	PointsIntersection(px, py, x);
366	//set begin and end of line
367	ftrack_info.Clear();
368	ftrack_info.fLine_point1.Set(px[0],py[0]);
369	ftrack_info.fLine_point2.Set(px[1],py[1]);
370
371	//d - length of line
372	double d = sqrt(pow(px[0]-px[1],2.)+pow(py[0]-py[1],2.));
373	double pix_size = 3.375; //pixel size
374	int av_counts=0;
375	int n = (d/pix_size); //bins number along l
376	int m = n*2.; //bins number across l
377	double k=d*0.5;
378	TH2D* h = new TH2D(0,"hm",n,0.,d,m,-k,k);
379	map <int, TVector3 >::iterator itpixels;
380	for ( itpixels=pixels.begin(); itpixels!=pixels.end(); itpixels++) {
381	TVector3& pixel_info= itpixels->second;
382	double dist = (pixel_info.Y() - apixel_info.X() - b)/Sqrt(1.+aa);
383	double y = pixel_info.Y()-distSqrt(1.+aa);
384	double x = ftrack_info.fLine.GetX(y);
385	double l =sqrt(pow(x-px[0],2.)+pow(py[0]-y,2.));
386	h->Fill(l,dist,pixel_info.Z());
387	pixel_info.SetXYZ(l,dist,pixel_info.Z());
388	av_counts+=pixel_info.Z();
389
390	}
391	av_counts/=pixels.size();
392
393	TH1D* hl = new TH1D(0,"Maximum counts along l",n,0.,d);
394	GetProjectionX(h,hl);
395
396	TH1D* hd = new TH1D(0,"Maximum counts across l",m,-k,k);
397	GetProjectionY(h,hd);
398
399
400	#ifdef RobustDebug
401	TCanvas* c = new TCanvas("c");
402	c->Print("max_counts.eps[");
403	h->SetTitle("color - tne number of counts");
404	h->GetXaxis()->SetTitle("l- length along line (bin = 1 size of pixel (3.375 mm))");
405	h->GetYaxis()->SetTitle("d - length across line (bin = 0.5 size of pixel");
406	h->Draw("colz");
407	c->Print("max_counts.eps");
408	h->Draw("LEGO");
409	c->Print("max_counts.eps");
410	hl->SetTitle("Maximum counts along l ");
411	hl->GetXaxis()->SetTitle("l- length along line (mm)");
412	hl->GetYaxis()->SetTitle("N_{counts}");
413	hl->Draw();
414	c->Print("max_counts.eps");
415	hd->SetTitle("Maximum counts across l ");
416	hd->GetXaxis()->SetTitle("d- length across line (mm)");
417	hd->GetYaxis()->SetTitle("N_{counts}");
418	hd->Draw();
419	c->Print("max_counts.eps");
420	#endif
421
422	TF1 *fgaus = new TF1("fgaus","gaus(0)+[3]");
423	int mbin=hl->GetMaximumBin();
424	fgaus->SetParameter(0,hl->GetBinContent(mbin));
425	fgaus->SetParameter(1,hl->GetBinCenter(mbin));
426	fgaus->SetParameter(2,40.);//half pdm81.
427	fgaus->SetParameter(3,av_counts);
428
429	TFitResultPtr fit = hl->Fit("fgaus","MSWQ0");
430	const double *par=fit->GetParams();
431	#ifdef RobustDebug
432	hl->Draw();
433	c->Print("max_counts.eps");
434	#endif
435
436	double bg=par[3];
437	SubtractBg(hl,bg,2.);
438	SubtractBg(hd,bg,1.);
439
440	#ifdef RobustDebug
441	hl->GetFunction("fgaus")->Delete();
442	hl->SetTitle("Maximum counts along l - BG ");
443	hd->SetTitle("Maximum counts across l - BG");
444	hl->Draw();
445
446	c->Print("max_counts.eps");
447	hd->Draw();
448	c->Print("max_counts.eps");
449	#endif
450	double l1, l2;
451	GetLimits(hl,par[1],0.98,l1,l2);
452
453
454	mbin=hd->GetMaximumBin();
455	fgaus->SetParameter(0,hd->GetBinContent(mbin));
456	fgaus->SetParameter(1,hd->GetBinCenter(mbin));
457	fgaus->SetParameter(2,10.); // pmt size
458	fgaus->SetParameter(3,0.);
459	fit=hd->Fit("fgaus","MSWQ0");
460	par=fit->GetParams();
461	#ifdef RobustDebug
462	hd->Draw();
463	c->Print("max_counts.eps");
464	TH1D* hcopy = new TH1D("hlcopy","Maximum counts along l",n,0.,d);
465
466	for(int i(hl->FindBin(l1)); i<=hl->FindBin(l2);i++){
467	hcopy->SetBinContent(i,hl->GetBinContent(i));
468
469	}
470	hcopy->SetFillColor(2);
471	hl->Draw();
472	hcopy->Draw("same");
473	c->Print("max_counts.eps");
474
475	int s1 = 2*(fabs(fx1-fx2)/pix_size);
476	int s2 = 2*(fabs(fy1-fy2)/pix_size);
477	TH2D* hcell = new TH2D("cell","",s1,fx1,fx2,s2,fy1,fy2);
478	TH2D* hcell_c = new TH2D("cell copy","",s1,fx1,fx2,s2,fy1,fy2);
479	#endif
480
481
482	map< int, map< int, vector<RecoPixelData*> > > cell = fEvent->GetCellGtuRecoPixelData();
483	map< int, map< int, vector<RecoPixelData*> > >::iterator itcell;
484	for(itcell= cell.begin();itcell!=cell.end();itcell++){
485	map< int, vector<RecoPixelData*> >& pixel = itcell->second;
486	map< int, vector<RecoPixelData*> >::iterator itpixel;
487
488	for (itpixel= pixel.begin();itpixel!=pixel.end();itpixel++) {
489	vector<RecoPixelData*>& vpix = itpixel->second;
490	vector<RecoPixelData*> new_pix;
491	for(unsigned int i=0;i<vpix.size();i++){
492	RecoPixelData* pix = vpix[i];
493	TVector3 pixel_info = fRunpars->PixelCenter(pix->GetPixelId());
494	#ifdef RobustDebug
495	hcell->Fill(pixel_info.X(),pixel_info.Y(),pix->GetCounts());
496	#endif
497	double dist = (pixel_info.Y() - apixel_info.X() - b)/Sqrt(1.+aa);
498	double y = pixel_info.Y()-distSqrt(1.+aa);
499	double x = ftrack_info.fLine.GetX(y);
500	double l =sqrt(pow(x-px[0],2.)+pow(py[0]-y,2.));
501	if(l>=(l1+minTOLERANCE) && l<=(l2-minTOLERANCE) && abs(dist)<=fabs(par[2]*3.) ){
502	new_pix.push_back(pix);
503
504	#ifdef RobustDebug
505	hcell_c->SetBinContent(hcell_c->GetXaxis()->FindBin(pixel_info.X()),hcell_c->GetYaxis()->FindBin(pixel_info.Y()),1);
506	#endif
507	}
508	}
509
510	ftrack_info.fSignalPerGtu[itpixel->first]=new_pix;
511
512	}
513	}
514
515	sigpix.clear();
516	map< int, vector<RecoPixelData*> >::iterator itsignal_pixels;
517	for (itsignal_pixels=ftrack_info.fSignalPerGtu.begin();itsignal_pixels!=ftrack_info.fSignalPerGtu.end();itsignal_pixels++) {
518	vector<RecoPixelData*>& vpix = itsignal_pixels->second;
519	for (unsigned int i=0; i<vpix.size();i++ ) {
520	RecoPixelData* pix = vpix[i];
521
522	int id = pix->GetPixelId();
523	bool signal = true;
524	for (unsigned int i=0; i<sigpix.size();i++ ) {
525	if(sigpix[i]==id){
526	signal=false;
527	break;
528	}
529	}
530	if(signal)sigpix.push_back(id);
531
532	}
533	}
534	ftrack_info.fSig_point1=l1;
535	ftrack_info.fSig_point2=l2;
536
537	#ifdef RobustDebug
538	hcell->SetStats(0);
539	hcell->Draw("colz");
540	c->Print("max_counts.eps");
541	hcell_c->SetMarkerStyle(2);
542	hcell_c->SetMarkerSize(0.5);
543	hcell->Draw("colz");
544	hcell_c->Draw("same");
545	c->Print("max_counts.eps");
546	c->Print("max_counts.eps]");
547	delete c;
548	delete hcell;
549	delete hcell_c;
550	delete hcopy;
551	#endif
552
553	delete h;
554	delete hl;
555	delete hd;
556	delete fgaus;
557
558	if ( !sigpix.size() ) {
559	Msg(EsafMsg::Info) << "None of the signal pixels was found" << MsgDispatch;
560	return kFALSE;
561	}
562	Msg(EsafMsg::Info) << "Number of the signal pixels:"<< sigpix.size() << MsgDispatch;
563	return kTRUE;
564	}
565
566	//_____________________________________________________________________________
567	void RobustModule::GetProjectionX(TH2D* h,TH1D* hh){
568	TH1D* h1 = new TH1D(0,"h1",h->GetNbinsY(),h->GetYaxis()->GetXmin(),h->GetYaxis()->GetXmax());
569	for(int i(1); i<=h->GetNbinsX();i++){
570	for(int j(1); j<=h->GetNbinsY();j++){
571	h1->Fill(h->GetYaxis()->GetBinCenter(j),h->GetBinContent(i,j));
572	}
573	hh->Fill(h->GetXaxis()->GetBinCenter(i),h1->GetBinContent(h1->GetMaximumBin()));
574	h1->Reset();
575	}
576
577	delete h1;
578	}
579
580	//_____________________________________________________________________________
581	void RobustModule::GetProjectionY(TH2D* h,TH1D* hh){
582	TH1D* h1 = new TH1D(0,"h1",h->GetNbinsX(),h->GetXaxis()->GetXmin(),h->GetXaxis()->GetXmax());
583	for(int i(1); i<=h->GetNbinsY();i++){
584	for(int j(1); j<=h->GetNbinsX();j++){
585	h1->Fill(h->GetXaxis()->GetBinCenter(j),h->GetBinContent(j,i));
586
587	}
588	hh->Fill(h->GetYaxis()->GetBinCenter(i),h1->GetBinContent(h1->GetMaximumBin()));
589	h1->Reset();
590	}
591	delete h1;
592	}
593
594	//_____________________________________________________________________________
595	void RobustModule::GetLimits(TH1D* h,double par, double opt, double& l1, double& l2){
596	int maximum = h->GetXaxis()->FindBin(par);
597	double integral=h->Integral();
598	double r=integral*opt; //half 90% of integral
599	double limit(0.);
600	int j=maximum;
601	int i=maximum-1;
602	do{
603	limit+=h->GetBinContent(j)+h->GetBinContent(i);
604	j++;
605	i--;
606	}while(limit<r && j<h->GetNbinsX() && i>0);
607
608	l1 = h->GetBinCenter(i+1);
609	l2 = h->GetBinCenter(j-1);
610	}
611
612	//_____________________________________________________________________________
613	void RobustModule::SubtractBg(TH1D* h, double bg, double power){
614	for(int i(1); i<=h->GetNbinsX();i++){
615	double c = h->GetBinContent(i)-bg;
616	h->SetBinContent(i,c<0.?0.:pow(c,power));
617	}
618	}
619
620	//_____________________________________________________________________________
621	void RobustModule::Coefficients(double &a, double &b,const map <int, TVector3 >& pixels){
622	// minimization of (dist)**2
623	// dist = (y - ax - b)/Sqrt(1.+aa)
624
625	double y, x, N, x_y,x_2,y_2;
626	y=x=N=x_y=x_2=y_2=0.;
627
628	int j=0;
629	map <int, TVector3 >::const_iterator itpixels;
630	for ( itpixels=pixels.begin(); itpixels!=pixels.end(); itpixels++) {
631	const TVector3& pixel_info = itpixels->second;
632	y += fWeights[j]*pixel_info.Y();
633	x += fWeights[j]*pixel_info.X();
634	x_2 += fWeights[j]pixel_info.X()pixel_info.X();
635	y_2 += fWeights[j]pixel_info.Y()pixel_info.Y();
636	x_y += fWeights[j]pixel_info.Y()pixel_info.X();
637	N += fWeights[j];
638	j++;
639	}
640
641	double det = sqrt((xx-yy+Ny_2-Nx_2)(xx-yy+Ny_2-Nx_2)-4.(xy-Nx_y)(Nx_y-x*y));
642	double a1(0.),a2(0.), b1(0.),b2(0.);
643	double den=2.(xy-N*x_y);
644
645	if(det==0.){
646	a1=-(xx-yy+Ny_2-Nx_2)/den;
647	b1=(y-a*x)/N;
648	a=a1;
649	b=b1;
650	}
651	else{
652	a1 = (-xx+yy-Ny_2+Nx_2+det)/den;
653	a2 = (-xx+yy-Ny_2+Nx_2-det)/den;
654
655	b1=(y-a1*x)/N;
656	b2=(y-a2*x)/N;
657
658	double chi1 = chi(pixels, a1, b1);
659	double chi2 = chi(pixels, a2, b2);
660
661	if(chi2>chi1){a=a1;b=b1;}
662	else{a=a2;b=b2;}
663	}
664	}
665
666	//_____________________________________________________________________________
667	double RobustModule::chi(const map <int, TVector3 >& pixels,double a,double b){
668	int j=0;
669	double chi=0.;
670	map <int, TVector3 >::const_iterator itpixels;
671	for ( itpixels=pixels.begin(); itpixels!=pixels.end(); itpixels++) {
672	const TVector3& pixel_info = itpixels->second;
673	double dist = (pixel_info.Y() - apixel_info.X() - b)/Sqrt(1.+aa);
674	chi+=fWeights[j]distdist ;
675	j++;
676	}
677	return chi;
678
679	}
680
681	//_____________________________________________________________________________
682	void RobustModule::PointsIntersection(double px, double py, double *x){
683	double x1, x2, y1, y2;
684	x1=x[0]; y1=x[1];
685	x2=x[2]; y2=x[3];
686
687	double x3[4]={fx1,fx2,fx1,fx1};
688	double x4[4]={fx1,fx2,fx2,fx2};
689	double y3[4]={fy1,fy1,fy1,fy2};
690	double y4[4]={fy2,fy2,fy1,fy2};
691
692	int j=0;
693	for(int i=0;i<4;i++){
694	double x=((x1y2-y1x2)(x3[i]-x4[i])-(x1-x2)(x3[i]y4[i]-y3[i]x4[i]))/((x1-x2)(y3[i]-y4[i])-(y1-y2)(x3[i]-x4[i]));
695	double y=((x1y2-y1x2)(y3[i]-y4[i])-(y1-y2)(x3[i]y4[i]-y3[i]x4[i]))/((x1-x2)(y3[i]-y4[i])-(y1-y2)(x3[i]-x4[i]));
696	if((x<=fx2+minTOLERANCE && x>=fx1-minTOLERANCE) && (y<=fy2+minTOLERANCE && y>=fy1-minTOLERANCE)){
697	px[j]=x;
698	py[j]=y;
699	j++;
700	}
701	}
702	}
703
704	//_____________________________________________________________________________

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: JEM-EUSO/esaf_lal/tags/v1_r0/esaf/packages/reconstruction/modules/base/clustering/src/RobustModule.cc @ 117

Download in other formats: