Context Navigation

source: trunk/MagicSoft/Mars/macros/pointspreadfunction.C@ 3117

Visit:

Last change on this file since 3117 was 3114, checked in by jlopez, 21 years ago
* empty log message *
File size: 21.6 KB

Line
1	/* ======================================================================== *\
2	!
3	! *
4	! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5	! * Software. It is distributed to you in the hope that it can be a useful
6	! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7	! * It is distributed WITHOUT ANY WARRANTY.
8	! *
9	! * Permission to use, copy, modify and distribute this software and its
10	! * documentation for any purpose is hereby granted without fee,
11	! * provided that the above copyright notice appear in all copies and
12	! * that both that copyright notice and this permission notice appear
13	! * in supporting documentation. It is provided "as is" without express
14	! * or implied warranty.
15	! *
16	!
17	! Author(s): Javier Lopez, 12/2003 <mailto:jlopez@ifae.es>
18	! Author(s): Alex Armada, 1/2004 <mailto:armada@ifae.es>
19	!
20	! Copyright: MAGIC Software Development, 2000-2003
21	!
22	!
23	\* ======================================================================== */
24
25
26	//------------------------------------------------------------------------- //
27	// //
28	// This macro fits the dc signal of a star using a two dimension gaussian //
29	// for each dc measurement. Then the values of parameters of the fit are //
30	// stored in histograms and shown at the end of the macro. //
31	// //
32	// USAGE: //
33	// It has two arguments, //
34	// 1- The first one is the dc file with the tracked star //
35	// 2- The second one is a continuos light file used to intercalibrate //
36	// the gain of the photomultipliers. //
37	// (It's possible not to use the calibration file and then the gain //
38	// of the pmts are supouse to be the same for all of them. //
39	// 3- The third argument is just the number of dc measurements you want //
40	// analize. If you put 0 it just stops after each fit and show you //
41	// results. //
42	// //
43	//--------------------------------------------------------------------------//
44
45	const Int_t numPixels = 577;
46	Int_t nPixelsFitted;
47	Bool_t isPixelsFitted[numPixels];
48	Float_t z[numPixels],x[numPixels],y[numPixels],errorz[numPixels];
49	Float_t chisquare;
50
51	//______________________________________________________________________________
52	//
53	// Function used by Minuit to do the fit
54	//
55	void fcn(Int_t &npar, Double_t gin, Double_t &f, Double_t par, Int_t iflag)
56	{
57	Int_t i;
58
59	//calculate chisquare
60	Double_t chisq = 0;
61	Double_t delta;
62	nPixelsFitted=0;
63	for (i=1;i<numPixels; i++) {
64	if (isPixelsFitted[i])
65	{
66	if (errorz[i] != 0.0)
67	{
68	delta = (z[i]-func(x[i],y[i],par))/errorz[i];
69	chisq += delta*delta;
70	nPixelsFitted++;
71	}
72	else
73	cout << "This should never happen errorz[" << i << "] " << errorz[i] << endl;
74	}
75	}
76	f = chisq;
77	chisquare = chisq;
78	}
79
80	//______________________________________________________________________________
81	//
82	// The 2D gaussian fucntion used to fit the spot of the star
83	//
84	Double_t func(float x,float y,Double_t *par)
85	{
86	Double_t value=par[0]TMath::exp(-(x-par[1])(x-par[1])/(2par[2]par[2]))TMath::exp(-(y-par[3])(y-par[3])/(2par[4]par[4]));
87	return value;
88	}
89
90	Bool_t HandleInput()
91	{
92	TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
93
94	while (1)
95	{
96	//
97	// While reading the input process gui events asynchronously
98	//
99	timer.TurnOn();
100	TString input = Getline("Type 'q' to exit, <return> to go on: ");
101	timer.TurnOff();
102
103	if (input=="q\n")
104	return kFALSE;
105
106	if (input=="\n")
107	return kTRUE;
108	};
109
110	return kFALSE;
111	}
112
113
114	void pointspreadfunction(TString fname, TString clname = "NULL", Int_t userNumLines = 1000)
115	{
116
117
118	Float_t currmean[numPixels];
119	Float_t currsquaremean[numPixels];
120	Float_t currrms[numPixels];
121	Float_t meanofcurrmean = 0;
122
123	for (Int_t swpix=0; swpix<numPixels; swpix++)
124	{
125	currmean[swpix] = 0.;
126	currsquaremean[swpix] = 0.;
127	currrms[swpix] = 0.;
128	}
129
130	Int_t numLines=0;
131
132	//containers
133	MGeomCamMagic geomcam;
134	MCameraDC curr;
135	MCameraDC currbis;
136
137	const Float_t conv4mm2deg = geomcam.GetConvMm2Deg();
138
139
140	if (clname != "NULL")
141	{
142
143	//
144	// First run over continuos light files to have a DC calibration
145	//
146
147	MParList plist0;
148
149	MTaskList tlist0;
150	plist0.AddToList(&tlist0);
151
152	plist0.AddToList(&geomcam);
153	plist0.AddToList(&curr);
154
155	MReportFileRead read0(clname);
156	read0.SetHasNoHeader();
157	read0.AddToList("MReportCurrents");
158
159	tlist0.AddToList(&read0);
160
161	MEvtLoop evtloop0;
162	evtloop0.SetParList(&plist0);
163
164
165	if (!evtloop0.PreProcess())
166	return;
167
168	while (tlist0.Process())
169	{
170	for (Int_t swpix=0; swpix<numPixels; swpix++)
171	{
172	meanofcurrmean += curr[swpix];
173	currmean[swpix] += curr[swpix];
174	currsquaremean[swpix] += curr[swpix]*curr[swpix];
175	}
176	numLines++;
177	}
178
179	evtloop0.PostProcess();
180
181	meanofcurrmean /= (numLines*numPixels);
182	for (Int_t swpix=0; swpix<numPixels; swpix++)
183	{
184
185	currmean[swpix] /= numLines;
186	currsquaremean[swpix] /= numLines;
187	currrms[swpix] = sqrt(fabs(currsquaremean[swpix] - currmean[swpix]*currmean[swpix]));
188
189	curr[swpix] = currmean[swpix];
190	currbis[swpix] = currrms[swpix];
191
192	currmean[swpix] /= meanofcurrmean;
193	currrms[swpix] /= meanofcurrmean;
194
195	}
196
197
198	/* MHCamera display0(geomcam);
199	display0.SetPrettyPalette();
200	display0.Draw();
201
202	// curr.Print();
203	display0.SetCamContent(currbis);
204	cout << "PSF>> DC mean values drawn" << endl;
205	// Remove the comments if you want to go through the file
206	// event-by-event:
207	if (!HandleInput())
208	break;
209	*/
210	}
211	else
212	{
213	// If you don't use the continuous light this is the currrms[] array
214	// is the error associated to the currents TMinuit will use.
215	for (Int_t swpix=0; swpix<numPixels; swpix++)
216	{
217	currmean[swpix] = 1.;
218	currrms[swpix] = 0.2;
219	}
220
221	}
222
223	//
224	// Now we can run over the dc data to extract the psf.
225	//
226	const Int_t maxNumLines = 10000;
227
228	Double_t ux[maxNumLines];
229	Double_t uy[maxNumLines];
230	Double_t sx[maxNumLines];
231	Double_t sy[maxNumLines];
232	Double_t chisqu[maxNumLines];
233	Double_t time[maxNumLines];
234
235	MParList plist;
236
237	MGeomCamMagic geomcam;
238	MCameraDC curr;
239	MTaskList tlist;
240
241	plist.AddToList(&geomcam);
242	plist.AddToList(&curr);
243	plist.AddToList(&tlist);
244
245	MReportFileRead read(fname);
246	read.SetHasNoHeader();
247	read.AddToList("MReportCurrents");
248
249	tlist.AddToList(&read);
250
251	MEvtLoop evtloop;
252	evtloop.SetParList(&plist);
253
254	if (!evtloop.PreProcess())
255	return;
256
257	MHCamera display(geomcam);
258	display.SetPrettyPalette();
259	display.Draw();
260	gPad->SetLogy();
261	gPad->cd(1);
262
263	Double_t amin,edm,errdef;
264	Int_t nvpar,nparx,icstat;
265
266	Double_t max,maxerror;
267	Double_t xmean,xsigma,ymean,ysigma;
268	Double_t xmeanerror,xsigmaerror,ymeanerror,ysigmaerror;
269
270	TEllipse ellipse;
271	ellipse.SetFillStyle(4000);
272	ellipse.SetLineWidth(2);
273	ellipse.SetLineColor(2);
274
275	ellipse.Draw();
276
277	Int_t nbinsxy = 80;
278	Float_t minxy = -600*conv4mm2deg;
279	Float_t maxxy = 600*conv4mm2deg;
280	Float_t fromdegtobin = (maxxy-minxy)/nbinsxy;
281
282	TH2D psfhist("psfhist","",nbinsxy,minxy,maxxy,nbinsxy,minxy,maxxy);
283	psfhist->GetXaxis()->SetTitle("[deg]");
284	psfhist->GetYaxis()->SetTitle("[deg]");
285	psfhist->GetZaxis()->SetTitle("DC [uA]");
286
287	TCanvas *psfcanvas = new TCanvas("psfcanvas","Point Spread Funtion 2D",200,20,900,700);
288
289
290	//
291	// Using the first dc measurement we search the pixels which contains the star and define
292	// an area to be fitted by Minuit which is 3 rings of neightbords around the peak of the star.
293	//
294
295	tlist.Process();
296
297	Int_t numLines=0;
298	Float_t minDCStar = 6.0;
299
300	Int_t numPixelsInStar = 0;
301	Float_t maxDC = 0;
302	Int_t swpixelmaxDC;
303
304	Bool_t isPixelsFittedTmp[numPixels];
305
306	for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
307	isPixelsFittedTmp[swpixel] = kFALSE;
308
309	for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
310	{
311	if(curr[swpixel] > maxDC)
312	{
313	swpixelmaxDC = swpixel;
314	maxDC = curr[swpixelmaxDC];
315	}
316
317	if(curr[swpixel]>minDCStar)
318	{
319	numPixelsInStar++;
320	isPixelsFitted[swpixel] = kTRUE;
321	}
322	else
323	isPixelsFitted[swpixel] = kFALSE;
324	}
325
326	if (numPixelsInStar == 0)
327	{
328	cout << "PSF>> Warning: none pixel over minDCStar(" << minDCStar << ')' << endl;
329	return;
330	}
331
332	//1st neighboor ring
333	for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
334	if (isPixelsFitted[swpixel])
335	for(Int_t next=0; next<geomcam[swpixel].GetNumNeighbors(); next++)
336	isPixelsFittedTmp[geomcam[swpixel].GetNeighbor(next)] = kTRUE;
337
338	for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
339	if (isPixelsFittedTmp[swpixel])
340	isPixelsFitted[swpixel] = kTRUE;
341
342	//2on neighboor ring
343	for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
344	if (isPixelsFitted[swpixel])
345	for(Int_t next=0; next<geomcam[swpixel].GetNumNeighbors(); next++)
346	isPixelsFittedTmp[geomcam[swpixel].GetNeighbor(next)] = kTRUE;
347
348	for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
349	if (isPixelsFittedTmp[swpixel])
350	isPixelsFitted[swpixel] = kTRUE;
351
352
353	//3rt neighboor ring
354	for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
355	if (isPixelsFitted[swpixel])
356	for(Int_t next=0; next<geomcam[swpixel].GetNumNeighbors(); next++)
357	isPixelsFittedTmp[geomcam[swpixel].GetNeighbor(next)] = kTRUE;
358
359	for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
360	if (isPixelsFittedTmp[swpixel])
361	isPixelsFitted[swpixel] = kTRUE;
362
363
364	for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
365	curr[swpixel] = (Float_t)isPixelsFitted[swpixel];
366
367	/* MHCamera display0(geomcam);
368	display0.SetPrettyPalette();
369	display0.Draw();
370
371	display0.SetCamContent(curr);
372	cout << "PSF>> Fitted pixels drawn" << endl;
373	// Remove the comments if you want to go through the file
374	// event-by-event:
375	if (!HandleInput())
376	break;
377	*/
378
379	// Minuit initialization
380
381	TMinuit *gMinuit = new TMinuit(7); //initialize TMinuit with a maximum of 5 params
382	gMinuit->SetFCN(fcn);
383
384	Double_t arglist[10];
385	Int_t ierflg = 0;
386
387	arglist[0] = 1;
388	gMinuit->mnexcm("SET ERR", arglist ,1,ierflg);
389	// arglist[0] = -1;
390	arglist[0] = 0;
391	gMinuit->mnexcm("SET PRI", arglist ,1,ierflg);
392
393	// Set starting values and step sizes for parameters
394	Double_t vstart[5];
395	Double_t step[5];
396	Double_t lowlimit[5] = {minDCStar, -2., 0.05, -2, 0.05};
397	Double_t uplimit[5] = {30., 2., 1., 2., 1.};
398
399	vstart[0] = maxDC;
400	vstart[1] = geomcam[swpixelmaxDC].GetX()*conv4mm2deg;
401	vstart[2] = 30TMath::Sqrt(numPixelsInStar/2)conv4mm2deg;
402	vstart[3] = geomcam[swpixelmaxDC].GetY()*conv4mm2deg;
403	vstart[4] = 30TMath::Sqrt(numPixelsInStar/2)conv4mm2deg;
404
405	for(Int_t i=0; i<5; i++)
406	{
407	if (vstart[i] != 0)
408	step[i] = TMath::Abs(vstart[i]/sqrt(2));
409	else
410	step[i] = uplimit[i]/2;
411	}
412
413	gMinuit->mnparm(0, "max", vstart[0], step[0], lowlimit[0], uplimit[0],ierflg);
414	gMinuit->mnparm(1, "xmean", vstart[1], step[1], lowlimit[1], uplimit[1],ierflg);
415	gMinuit->mnparm(2, "xsigma", vstart[2], step[2], lowlimit[2], uplimit[2],ierflg);
416	gMinuit->mnparm(3, "ymean", vstart[3], step[3], lowlimit[3], uplimit[3],ierflg);
417	gMinuit->mnparm(4, "ysigma", vstart[4], step[4], lowlimit[4], uplimit[4],ierflg);
418
419	while (tlist.Process() && numLines < maxNumLines)
420	{
421
422	for (Int_t swpixel=1; swpixel<577; swpixel++)
423	{
424
425	x[swpixel] = geomcam[swpixel].GetX()*conv4mm2deg;
426	y[swpixel] = geomcam[swpixel].GetY()*conv4mm2deg;
427	z[swpixel] = curr[swpixel]/currmean[swpixel];
428	errorz[swpixel] = TMath::Sqrt((curr[swpixel]currrms[swpixel]/(currmean[swpixel]currmean[swpixel]))(curr[swpixel]currrms[swpixel]/(currmean[swpixel]currmean[swpixel]))+(0.1)/currmean[swpixel](0.1)/currmean[swpixel]);
429
430
431	psfhist->SetBinContent((Int_t)((x[swpixel]+600conv4mm2deg)/fromdegtobin),(Int_t)((y[swpixel]+600conv4mm2deg)/fromdegtobin),z[swpixel]);
432	}
433
434	psfcanvas->cd(1);
435	psfhist->Draw("lego2");
436
437	// Now ready for minimization step
438	arglist[0] = 500;
439	arglist[1] = 1.;
440	gMinuit->mnexcm("MIGRAD", arglist ,2,ierflg);
441
442	// Print results
443	/* gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
444	gMinuit->mnprin(3,amin);
445	*/
446	gMinuit->GetParameter(0,max,maxerror);
447	gMinuit->GetParameter(1,xmean,xmeanerror);
448	gMinuit->GetParameter(2,xsigma,xsigmaerror);
449	gMinuit->GetParameter(3,ymean,ymeanerror);
450	gMinuit->GetParameter(4,ysigma,ysigmaerror);
451
452	/* cout << endl;
453	cout << "numLine " << numLines << endl;
454	cout << "max \t" << max << " +- " << maxerror << endl;
455	cout << "xmean \t" << xmean << " +- " << xmeanerror << endl;
456	cout << "xsigma \t" << TMath::Abs(xsigma) << " +- " << xsigmaerror << endl;
457	cout << "ymean \t" << ymean << " +- " << ymeanerror << endl;
458	cout << "ysigma \t" << TMath::Abs(ysigma) << " +- " << ysigmaerror << endl;
459	cout << "chisquare/ndof \t" << chisquare/(nPixelsFitted-5) << endl;
460	*/
461
462	chisqu[numLines]=chisquare/(nPixelsFitted-5);
463
464	if(chisqu[numLines]<100.)
465	{
466	ux[numLines]=xmean;
467	uy[numLines]=ymean;
468	sx[numLines]=TMath::Abs(xsigma);
469	sy[numLines]=TMath::Abs(ysigma);
470	time[numLines]=numLines;
471
472	display.SetCamContent(curr);
473	gPad->cd(1);
474	ellipse.SetX1(xmean/conv4mm2deg);
475	ellipse.SetY1(ymean/conv4mm2deg);
476	ellipse.SetR1(TMath::Abs(xsigma)/conv4mm2deg);
477	ellipse.SetR2(TMath::Abs(ysigma)/conv4mm2deg);
478
479	gPad->Modified();
480	gPad->Update();
481
482	// Remove the comments if you want to go through the file
483	//event-by-event:
484	if (userNumLines>0)
485	{
486	if (numLines>userNumLines)
487	break;
488	}
489	else
490	{
491	if (!HandleInput())
492	break;
493	}
494	numLines++;
495	}
496	}
497
498
499	evtloop.PostProcess();
500
501	//
502	// Draw the ditributions of the sigmas the point spread function
503	//
504
505	cout<<"PSF>> Number of lines "<<numLines<<endl;
506
507	gROOT->Reset();
508	gStyle->SetCanvasColor(0);
509	gStyle->SetCanvasBorderMode(0);
510	gStyle->SetPadBorderMode(0);
511	gStyle->SetFrameBorderMode(0);
512	gStyle->SetOptStat(00000000);
513
514	c1 = new TCanvas("c1","Time evolution & distributions",0,0,1200,850);
515	c1->Divide(3,2);
516
517	c1->cd(1);
518
519	TMath math;
520
521	Double_t minmeanx, maxmeanx;
522	minmeanx = ux[math.LocMin(numLines,ux)];
523	maxmeanx = ux[math.LocMax(numLines,ux)];
524
525	Double_t minmeany, maxmeany;
526	minmeany = uy[math.LocMin(numLines,uy)];
527	maxmeany = uy[math.LocMax(numLines,uy)];
528
529	Double_t minmean, maxmean;
530	minmean = math.Min(minmeanx,minmeany);
531	maxmean = math.Max(maxmeanx,maxmeany);
532
533	Double_t diff;
534	diff = maxmean - minmean;
535	diff = 0.1*diff;
536	minmean = minmean - diff;
537	maxmean = maxmean + diff;
538
539	Double_t mintime, maxtime;
540	mintime = time[math.LocMin(numLines,time)];
541	maxtime = time[math.LocMax(numLines,time)];
542
543	TH2D *h1 = new TH2D("h1",fname,1,mintime-1,maxtime+1,1,minmean,maxmean);
544	h1->GetXaxis()->SetTitle("Event number");
545	h1->GetYaxis()->SetTitle("mean position (deg)");
546	h1->Draw();
547
548	TGraph *grtimeevolmeanx = new TGraph(numLines,time,ux);
549	grtimeevolmeanx->SetMarkerColor(3);
550	grtimeevolmeanx->SetMarkerStyle(20);
551	grtimeevolmeanx->SetMarkerSize (0.4);
552	grtimeevolmeanx->Draw("P");
553
554	TGraph *grtimeevolmeany = new TGraph(numLines,time,uy);
555	grtimeevolmeany->SetMarkerColor(6);
556	grtimeevolmeany->SetMarkerStyle(24);
557	grtimeevolmeany->SetMarkerSize (0.4);
558	grtimeevolmeany->Draw("P");
559
560	legmeanxy = new TLegend(0.8,0.85,0.95,0.95);
561	legmeanxy.SetTextSize(0.03);
562	legmeanxy.AddEntry(grtimeevolmeanx,"mean x","P");
563	legmeanxy.AddEntry(grtimeevolmeany,"mean y","P");
564	legmeanxy.Draw();
565
566	c1->cd(2);
567
568	TMath math;
569
570	Double_t minsigmax, maxsigmax;
571	minsigmax = sx[math.LocMin(numLines,sx)];
572	maxsigmax = sx[math.LocMax(numLines,sx)];
573
574	Double_t minsigmay, maxsigmay;
575	minsigmay = sy[math.LocMin(numLines,sy)];
576	maxsigmay = sy[math.LocMax(numLines,sy)];
577
578	Double_t minsigma, maxsigma;
579	minsigma = math.Min(minsigmax,minsigmay);
580	maxsigma = math.Max(maxsigmax,maxsigmay);
581
582	diff = maxsigma - minsigma;
583	diff = 0.1*diff;
584	minsigma = minsigma - diff;
585	maxsigma = maxsigma + diff;
586
587	TH2D *h2 = new TH2D("h2","",1,mintime-1,maxtime+1,1,minsigma,maxsigma);
588	h2->GetXaxis()->SetTitle("Event number");
589	h2->GetYaxis()->SetTitle("PSF Rms (deg)");
590	h2->Draw();
591
592	TGraph* grtimeevolsigmax= new TGraph(numLines,time,sx);
593	grtimeevolsigmax->SetMarkerColor(3);
594	grtimeevolsigmax->SetMarkerStyle(20);
595	grtimeevolsigmax->SetMarkerSize (0.4);
596	grtimeevolsigmax->Draw("P");
597
598	TGraph* grtimeevolsigmay= new TGraph(numLines,time,sy);
599	grtimeevolsigmay->SetMarkerColor(6);
600	grtimeevolsigmay->SetMarkerStyle(24);
601	grtimeevolsigmay->SetMarkerSize (0.4);
602	grtimeevolsigmay->Draw("P");
603
604	legsigmaxy = new TLegend(0.8,0.85,0.95,0.95);
605	legsigmaxy.SetTextSize(0.03);
606	legsigmaxy.AddEntry(grtimeevolsigmax,"sigma x","P");
607	legsigmaxy.AddEntry(grtimeevolsigmay,"sigma y","P");
608	legsigmaxy.Draw();
609
610	c1->cd(3);
611
612	Double_t minchisqu, maxchisqu;
613
614	minchisqu = chisqu[math.LocMin(numLines,chisqu)];
615	maxchisqu = chisqu[math.LocMax(numLines,chisqu)];
616
617	diff = maxchisqu - minchisqu;
618	diff = 0.1*diff;
619	minchisqu = minchisqu - diff;
620	maxchisqu = maxchisqu + diff;
621
622	TH2D *h3 = new TH2D("h3","",1,mintime-1,maxtime+1,1,minchisqu,maxchisqu);
623	h3->GetXaxis()->SetTitle("Event number");
624	h3->Draw();
625
626	TGraph * grtimeevolchisqu = new TGraph(numLines,time,chisqu);
627	grtimeevolchisqu->SetMarkerColor(215);
628	grtimeevolchisqu->SetMarkerStyle(20);
629	grtimeevolchisqu->SetMarkerSize(0.4);
630	grtimeevolchisqu->Draw("P");
631
632	legchisqu = new TLegend(0.55,0.90,0.95,0.95);
633	legchisqu.SetTextSize(0.03);
634	legchisqu.AddEntry(grtimeevolchisqu,"chi square / ndof","P");
635	legchisqu.Draw();
636
637	//***************************************
638
639	const Int_t nbins = 100;
640
641	TH1D *xsigmahist = new TH1D("xsigmahist","",nbins,minsigma,maxsigma);
642	TH1D *ysigmahist = new TH1D("ysigmahist","",nbins,minsigma,maxsigma);
643	TH1D *xmeanhist = new TH1D("xmeanhist","",nbins,minmean,maxmean);
644	TH1D *ymeanhist = new TH1D("ymeanhist","",nbins,minmean,maxmean);
645	TH1D *chisquhist = new TH1D("chisquhist","",nbins,minchisqu,maxchisqu);
646
647	for (Int_t i=0; i<numLines; i++)
648	{
649	xsigmahist->Fill(TMath::Abs(sx[i]));
650	ysigmahist->Fill(TMath::Abs(sy[i]));
651	xmeanhist->Fill(ux[i]);
652	ymeanhist->Fill(uy[i]);
653	chisquhist->Fill(chisqu[i]);
654
655	}
656
657	c1->cd(5);
658
659	TMath math;
660	Double_t maxsigma;
661	Int_t binmaxx, binmaxy;
662	xsigmahist->SetLineColor(3);
663	xsigmahist->SetLineWidth(2);
664	xsigmahist->SetXTitle("RMS [deg]");
665	binmaxx = xsigmahist->GetMaximumBin();
666	binmaxx = xsigmahist->GetBinContent(binmaxx);
667
668	ysigmahist->SetLineColor(6);
669	ysigmahist->SetLineWidth(2);
670	binmaxy = ysigmahist->GetMaximumBin();
671	binmaxy = ysigmahist->GetBinContent(binmaxy);
672
673	maxsigma = math.Max(binmaxx,binmaxy);
674	maxsigma += 1;
675
676	xsigmahist->SetMaximum(maxsigma);
677	ysigmahist->SetMaximum(maxsigma);
678	xsigmahist->DrawCopy();
679	ysigmahist->DrawCopy("Same");
680
681	TLegend *legendsigma = new TLegend(.3,.8,.95,.95);
682	legendsigma->SetTextSize(0.03);
683	char xsigmatitle[100];
684	char ysigmatitle[100];
685	sprintf(xsigmatitle,"PSF Rms on X axis -- %5.2f +/- %5.2f deg",xsigmahist->GetMean(),xsigmahist->GetRMS());
686	sprintf(ysigmatitle,"PSF Rms on Y axis -- %5.2f +/- %5.2f deg",ysigmahist->GetMean(),ysigmahist->GetRMS());
687	legendsigma->AddEntry(xsigmahist,xsigmatitle,"F");
688	legendsigma->AddEntry(ysigmahist,ysigmatitle,"F");
689	legendsigma->Draw();
690
691	c1->cd(4);
692
693	Double_t maxmean;
694
695	xmeanhist->SetLineColor(3);
696	xmeanhist->SetLineWidth(2);
697	xmeanhist->SetXTitle("mean [deg]");
698	binmaxx = xmeanhist->GetMaximumBin();
699	binmaxx = xmeanhist->GetBinContent(binmaxx);
700
701	ymeanhist->SetLineColor(6);
702	ymeanhist->SetLineWidth(2);
703	binmaxy = ymeanhist->GetMaximumBin();
704	binmaxy = ymeanhist->GetBinContent(binmaxy);
705
706	maxmean = math.Max(binmaxx,binmaxy);
707	maxmean += 1;
708
709	xmeanhist->SetMaximum(maxmean);
710	ymeanhist->SetMaximum(maxmean);
711	xmeanhist->DrawCopy();
712	ymeanhist->DrawCopy("Same");
713
714	TLegend *legendmean = new TLegend(.35,.8,.95,.95);
715	legendmean->SetTextSize(0.03);
716	char xmeantitle[100];
717	char ymeantitle[100];
718	sprintf(xmeantitle,"mean on X axis -- %5.2f +/- %5.2f deg",xmeanhist->GetMean(),xmeanhist->GetRMS());
719	sprintf(ymeantitle,"mean on Y axis -- %5.2f +/- %5.2f deg",ymeanhist->GetMean(),ymeanhist->GetRMS());
720	legendmean->AddEntry(xmeanhist,xmeantitle,"F");
721	legendmean->AddEntry(ymeanhist,ymeantitle,"F");
722	legendmean->Draw();
723
724	//meancanvas->Modified();
725	//meancanvas->Update();
726
727	c1->cd(6);
728
729	chisquhist->SetLineColor(215);
730	chisquhist->SetLineWidth(2);
731	chisquhist->SetXTitle("chi square / ndof");
732	TAxis * axis = chisquhist->GetXaxis();
733	axis->SetLabelSize(0.025);
734	chisquhist->DrawCopy();
735
736	TLegend *legendchisqu = new TLegend(.4,.85,.95,.95);
737	legendchisqu->SetTextSize(0.03);
738	char chisqutitle[100];
739	sprintf(chisqutitle,"chi square / ndof -- %5.2f +/- %5.2f ",chisquhist->GetMean(),chisquhist->GetRMS());
740	legendchisqu->AddEntry(chisquhist,chisqutitle,"F");
741	legendchisqu->Draw();
742
743
744	return;
745
746	}
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825

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

Download in other formats: