source: trunk/MagicSoft/Mars/mtemp/MHTelAxisFromStars.cc@ 4707

Last change on this file since 4707 was 4705, checked in by wittek, 21 years ago
*** empty log message ***
File size: 19.2 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!
18! Author(s): Wolfgang Wittek 07/2004 <mailto:wittek@mppmu.mpg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2004
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MHTelAxisFromStars
28//
29// This class contains histograms for MTelAxisFromStars
30//
31/////////////////////////////////////////////////////////////////////////////
32#include "MHTelAxisFromStars.h"
33
34#include <math.h>
35
36#include <TH1.h>
37#include <TH2.h>
38#include <TPad.h>
39#include <TStyle.h>
40#include <TCanvas.h>
41
42#include "MLog.h"
43#include "MLogManip.h"
44
45#include "MParList.h"
46
47#include "MGeomCam.h"
48#include "MBinning.h"
49
50#include "MStarCam.h"
51#include "MStarPos.h"
52#include "MSkyCamTrans.h"
53#include "MSrcPosCam.h"
54
55
56ClassImp(MHTelAxisFromStars);
57
58using namespace std;
59
60// --------------------------------------------------------------------------
61//
62// Setup the histograms
63//
64MHTelAxisFromStars::MHTelAxisFromStars(const char *name, const char *title)
65 : fMm2Deg(1), fUseMmScale(kFALSE)
66{
67 fName = name ? name : "MHTelAxisFromStars";
68 fTitle = title ? title : "Plots for MTelAxisFromStars";
69
70 fNStars = new TH1D("NStars", "No. of stars", 10, -0.5, 9.5);
71 fNStars->SetDirectory(NULL);
72 fNStars->SetXTitle("Numder of stars");
73 fNStars->SetYTitle("Counts");
74
75 fNdoF = new TH1D("NdoF", "No. of degrees of freedom", 20, -0.5, 19.5);
76 fNdoF->SetDirectory(NULL);
77 fNdoF->SetXTitle("Numder of degrees of freedom");
78 fNdoF->SetYTitle("Counts");
79
80 fLog10Chi2 = new TH1D("log10Chi2", "log10Chi2", 60, -10, 5);
81 fLog10Chi2->SetDirectory(NULL);
82 fLog10Chi2->SetXTitle("log10(Chi2)");
83 fLog10Chi2->SetYTitle("Counts");
84
85 fChi2Prob = new TH1D("Chi2-Prob", "Chi2 probability", 40, 0.0, 1.0);
86 fChi2Prob->SetDirectory(NULL);
87 fChi2Prob->SetXTitle("Chi2 probability");
88 fChi2Prob->SetYTitle("Counts");
89
90
91 fNumIter = new TH1D("NumIter", "Number of iterations", 50, -0.5, 49.5);
92 fNumIter->SetDirectory(NULL);
93 fNumIter->SetXTitle("Number of iterations");
94 fNumIter->SetYTitle("Counts");
95
96
97 fLambda = new TH1D("Lambda", "Scale factor lambda", 80, 0.90, 1.10);
98 fLambda->SetDirectory(NULL);
99 fLambda->SetXTitle("Scale factor lambda");
100 fLambda->SetYTitle("Counts");
101
102
103 fAlfa = new TH1D("Alfa", "Rotation angle alfa", 100, -2.5, 2.5);
104 fAlfa->SetDirectory(NULL);
105 fAlfa->SetXTitle("Rotation angle alfa [\\circ]");
106 fAlfa->SetYTitle("Counts");
107
108
109 fShift = new TH2D("Shift", "Sky-Cam transformnation : (x,y) shift",
110 72, -534.0707, 534.0707, 72, -534.0707, 534.0707);
111 fShift->SetDirectory(NULL);
112 fShift->SetZTitle("Counts");
113 fShift->SetXTitle("x-shift [\\circ]");
114 fShift->SetYTitle("y-shift [\\circ]");
115
116
117 fEstPos1 = new TH2D("EstPos1", "Estimated position1",
118 72, -534.0707, 534.0707, 72, -534.0707, 534.0707);
119 fEstPos1->SetDirectory(NULL);
120 fEstPos1->SetZTitle("Counts");
121 fEstPos1->SetXTitle("Estimated position1 x [\\circ]");
122 fEstPos1->SetYTitle("Estimated position1 y [\\circ]");
123
124 fEstPos2 = new TH2D("EstPos2", "Estimated position2",
125 72, -534.0707, 534.0707, 72, -534.0707, 534.0707);
126 fEstPos2->SetDirectory(NULL);
127 fEstPos2->SetZTitle("Counts");
128 fEstPos2->SetXTitle("Estimated position2 x [\\circ]");
129 fEstPos2->SetYTitle("Estimated position2 y [\\circ]");
130
131 fEstPos3 = new TH2D("EstPos3", "Estimated position3",
132 72, -534.0707, 534.0707, 72, -534.0707, 534.0707);
133 fEstPos3->SetDirectory(NULL);
134 fEstPos3->SetZTitle("Counts");
135 fEstPos3->SetXTitle("Estimated position3 x [\\circ]");
136 fEstPos3->SetYTitle("Estimated position3 y [\\circ]");
137
138
139 // default input type : results from Gauss fit
140 fInputType = 1;
141}
142
143// --------------------------------------------------------------------------
144//
145// Set the type of the input
146//
147// type = 0 calculated star positions (by averaging)
148// type = 1 fitted star positions (by Gauss fit)
149//
150void MHTelAxisFromStars::SetInputType(Int_t type)
151{
152 *fLog << all << "MHTelAxisFromStars::SetInputType; input type is set equal to : "
153 << type << endl;
154
155 fInputType = type;
156}
157
158
159// --------------------------------------------------------------------------
160//
161// Delete the histograms
162//
163MHTelAxisFromStars::~MHTelAxisFromStars()
164{
165 //*fLog << "MHTelAxisFromStars::Destructor" << endl;
166
167 delete fNStars;
168 delete fNdoF;
169 delete fLog10Chi2;
170 delete fChi2Prob;
171 delete fNumIter;
172 delete fLambda;
173 delete fAlfa;
174
175 delete fShift;
176 delete fEstPos1;
177 delete fEstPos2;
178 delete fEstPos3;
179}
180
181// --------------------------------------------------------------------------
182//
183// Setup the Binning for the histograms.
184//
185// Find the pointers to the containers.
186//
187// Use this function if you want to set the conversion factor which
188// is used to convert the mm-scale in the camera plane into the deg-scale.
189// The conversion factor is part of the camera geometry. Please create a
190// corresponding MGeomCam container.
191//
192Bool_t MHTelAxisFromStars::SetupFill(const MParList *plist)
193{
194
195 fStarCam = (MStarCam*)plist->FindObject("MStarCam", "MStarCam");
196 if (!fStarCam)
197 {
198 *fLog << err << "MHTelAxisFromStars::SetupFill; container 'MStarCam' not found... aborting." << endl;
199 return kFALSE;
200 }
201
202
203 fSourceCam = (MStarCam*)plist->FindObject("MSourceCam", "MStarCam");
204 if (!fSourceCam)
205 {
206 *fLog << err << "MHTelAxisFromStars::SetupFill; container 'MSourceCam' not found... " << endl;
207 }
208
209
210
211 fSrcPos = (MSrcPosCam*)plist->FindObject(AddSerialNumber("MSrcPosCam"));
212 if (!fSrcPos)
213 {
214 *fLog << err << "MHTelAxisFromStars::SetupFill; MSrcPosCam not found... aborting" << endl;
215 return kFALSE;
216 }
217
218 fSkyCamTrans = (MSkyCamTrans*)plist->FindObject(AddSerialNumber("MSkyCamTrans"));
219 if (!fSkyCamTrans)
220 {
221 *fLog << err << "MHTelAxisFromStars::SetupFill; MSkyCamTrans not found... aborting" << endl;
222 return kFALSE;
223 }
224
225
226 //------------------------------------------------------------------
227 const MGeomCam *geom = (MGeomCam*)plist->FindObject("MGeomCam");
228 if (!geom)
229 {
230 *fLog << warn << GetDescriptor()
231 << ": No Camera Geometry available. Using mm-scale for histograms."
232 << endl;
233 SetMmScale(kTRUE);
234 }
235 else
236 {
237 fMm2Deg = geom->GetConvMm2Deg();
238 SetMmScale(kFALSE);
239 }
240
241 ApplyBinning(*plist, "NStars", fNStars);
242 ApplyBinning(*plist, "NdoF", fNdoF);
243 ApplyBinning(*plist, "Log10Chi2",fLog10Chi2);
244 ApplyBinning(*plist, "Chi2Prob", fChi2Prob);
245 ApplyBinning(*plist, "NumIter", fNumIter);
246 ApplyBinning(*plist, "Lambda", fLambda);
247 ApplyBinning(*plist, "Alfa", fAlfa);
248
249 const MBinning *bins = (MBinning*)plist->FindObject("BinningCamera");
250 if (!bins)
251 {
252 float r = geom ? geom->GetMaxRadius() : 600;
253 r *= 0.9;
254 if (!fUseMmScale)
255 r *= fMm2Deg;
256
257 MBinning b;
258 b.SetEdges(61, -r, r);
259 SetBinning(fShift, &b, &b);
260 SetBinning(fEstPos1, &b, &b);
261 SetBinning(fEstPos2, &b, &b);
262 SetBinning(fEstPos3, &b, &b);
263 }
264 else
265 {
266 SetBinning(fShift, bins, bins);
267 SetBinning(fEstPos1, bins, bins);
268 SetBinning(fEstPos2, bins, bins);
269 SetBinning(fEstPos3, bins, bins);
270 }
271
272
273 // copy names from MStarCam to the histograms
274 if (fSourceCam)
275 {
276 *fLog << "Source container exists " << endl;
277
278 MStarPos *starSource = 0;
279 TIter nextSource(fSourceCam->GetList());
280
281 *fLog << "&nextSource = " << &nextSource << endl;
282
283 starSource = (MStarPos*)nextSource();
284 if (starSource)
285 {
286 *fLog << "name of source container1 : " << starSource->GetName()
287 << endl;
288
289 fEstPos1->SetName(starSource->GetName());
290
291 starSource = (MStarPos*)nextSource();
292 if (starSource)
293 {
294 *fLog << "name of source container2 : " << starSource->GetName()
295 << endl;
296
297
298 fEstPos2->SetName(starSource->GetName());
299
300 starSource = (MStarPos*)nextSource();
301 if (starSource)
302 {
303
304 *fLog << "name of source container3 : " << starSource->GetName()
305 << endl;
306
307 fEstPos3->SetName(starSource->GetName());
308 }
309 }
310 }
311 }
312
313
314 return kTRUE;
315}
316
317// --------------------------------------------------------------------------
318//
319// Use this function to setup your own conversion factor between degrees
320// and millimeters. The conversion factor should be the one calculated in
321// MGeomCam. Use this function with Caution: You could create wrong values
322// by setting up your own scale factor.
323//
324void MHTelAxisFromStars::SetMm2Deg(Float_t mmdeg)
325{
326 if (mmdeg<0)
327 {
328 *fLog << warn << dbginf
329 << "Warning - Conversion factor < 0 - nonsense. Ignored."
330 << endl;
331 return;
332 }
333
334 if (fMm2Deg>=0)
335 *fLog << warn << dbginf
336 << "Warning - Conversion factor already set. Overwriting"
337 << endl;
338
339 fMm2Deg = mmdeg;
340}
341
342// --------------------------------------------------------------------------
343//
344// With this function you can convert the histogram ('on the fly') between
345// degrees and millimeters.
346//
347void MHTelAxisFromStars::SetMmScale(Bool_t mmscale)
348{
349 if (fUseMmScale == mmscale)
350 return;
351
352 if (fMm2Deg<0)
353 {
354 *fLog << warn << dbginf
355 << "Warning - Sorry, no conversion factor for conversion available."
356 << endl;
357 return;
358 }
359
360 const Double_t scale = mmscale ? 1./fMm2Deg : fMm2Deg;
361 MH::ScaleAxis(fShift, scale, scale);
362 MH::ScaleAxis(fEstPos1, scale, scale);
363 MH::ScaleAxis(fEstPos2, scale, scale);
364 MH::ScaleAxis(fEstPos3, scale, scale);
365
366 if (mmscale)
367 {
368 fShift->SetXTitle("x-shift [mm]");
369 fShift->SetYTitle("y-shift [mm]");
370
371 fEstPos1->SetXTitle("Estimated position1 x [mm]");
372 fEstPos1->SetYTitle("Estimated position1 y [mm]");
373
374 fEstPos2->SetXTitle("Estimated position2 x [mm]");
375 fEstPos2->SetYTitle("Estimated position2 y [mm]");
376
377 fEstPos3->SetXTitle("Estimated position3 x [mm]");
378 fEstPos3->SetYTitle("Estimated position3 y [mm]");
379 }
380 else
381 {
382 fShift->SetXTitle("x-shift [\\circ]");
383 fShift->SetYTitle("y-shift [\\circ]");
384
385 fEstPos1->SetXTitle("Estimated position1 x [\\circ]");
386 fEstPos1->SetYTitle("Estimated position1 y [\\circ]");
387
388 fEstPos2->SetXTitle("Estimated position2 x [\\circ]");
389 fEstPos2->SetYTitle("Estimated position2 y [\\circ]");
390
391 fEstPos3->SetXTitle("Estimated position3 x [\\circ]");
392 fEstPos3->SetYTitle("Estimated position3 y [\\circ]");
393 }
394
395 fUseMmScale = mmscale;
396}
397
398// --------------------------------------------------------------------------
399//
400// Fill the histograms
401//
402//
403Bool_t MHTelAxisFromStars::Fill(const MParContainer *par, const Stat_t w)
404{
405 Int_t Ndof = fSkyCamTrans->GetNdof();
406 if (Ndof < 0)
407 return kTRUE;
408
409 const Double_t scale = fUseMmScale ? 1 : fMm2Deg;
410
411
412 fNStars ->Fill(fSkyCamTrans->GetNStars(), w);
413 fNdoF ->Fill(fSkyCamTrans->GetNdof(), w);
414 if (fSkyCamTrans->GetChiSquare() > 0.0)
415 fLog10Chi2->Fill( log10(fSkyCamTrans->GetChiSquare() ), w);
416 fChi2Prob ->Fill(fSkyCamTrans->GetChiSquareProb(), w);
417 fNumIter ->Fill(fSkyCamTrans->GetNumIter(), w);
418 fLambda ->Fill(fSkyCamTrans->GetLambda(), w);
419 fAlfa ->Fill(fSkyCamTrans->GetAlfa(), w);
420
421
422 Double_t lowx;
423 Double_t higx;
424 Double_t lowy;
425 Double_t higy;
426 Double_t x;
427 Double_t y;
428 Int_t nbix;
429 Int_t nbiy;
430
431 nbix = fShift->GetXaxis()->GetNbins();
432 lowx = fShift->GetXaxis()->GetBinLowEdge(1);
433 higx = fShift->GetXaxis()->GetBinLowEdge(nbix+1);
434
435 nbiy = fShift->GetYaxis()->GetNbins();
436 lowy = fShift->GetYaxis()->GetBinLowEdge(1);
437 higy = fShift->GetYaxis()->GetBinLowEdge(nbiy+1);
438
439 x = scale * (fSkyCamTrans->GetShiftD())[0];
440 y = scale * (fSkyCamTrans->GetShiftD())[1];
441 if (x>lowx && x<higx && y>lowy && y<higy)
442 {
443 fShift ->Fill(x, y, w);
444 }
445
446 if (fSourceCam)
447 {
448 MStarPos *starSource = 0;
449 TIter nextSource(fSourceCam->GetList());
450
451 starSource = (MStarPos*)nextSource();
452 if (starSource)
453 {
454 if (fInputType == 1)
455 {
456 x = scale * starSource->GetMeanXFit();
457 y = scale * starSource->GetMeanYFit();
458 }
459 else
460 {
461 x = scale * starSource->GetMeanXCalc();
462 y = scale * starSource->GetMeanYCalc();
463 }
464
465 if (x>lowx && x<higx && y>lowy && y<higy)
466 {
467 fEstPos1 ->Fill(x, y, w);
468 //fEstPos1->SetName(starSource->GetName());
469 }
470
471
472 starSource = (MStarPos*)nextSource();
473 if (starSource)
474 {
475 if (fInputType == 1)
476 {
477 x = scale * starSource->GetMeanXFit();
478 y = scale * starSource->GetMeanYFit();
479 }
480 else
481 {
482 x = scale * starSource->GetMeanXCalc();
483 y = scale * starSource->GetMeanYCalc();
484 }
485
486 if (x>lowx && x<higx && y>lowy && y<higy)
487 {
488 fEstPos2 ->Fill(x, y, w);
489 //fEstPos2->SetName(starSource->GetName());
490 }
491
492 starSource = (MStarPos*)nextSource();
493 if (starSource)
494 {
495 if (fInputType == 1)
496 {
497 x = scale * starSource->GetMeanXFit();
498 y = scale * starSource->GetMeanYFit();
499 }
500 else
501 {
502 x = scale * starSource->GetMeanXCalc();
503 y = scale * starSource->GetMeanYCalc();
504 }
505
506 if (x>lowx && x<higx && y>lowy && y<higy)
507 {
508 fEstPos3 ->Fill(x, y, w);
509 //fEstPos3->SetName(starSource->GetName());
510 }
511 }
512 }
513 }
514 }
515
516 return kTRUE;
517}
518
519// --------------------------------------------------------------------------
520//
521// Finalize :
522//
523// it is called in the postprocessing
524// reset pointers in order to allow cloning of the object
525//
526Bool_t MHTelAxisFromStars::Finalize()
527{
528 //*fLog << "MHTelAxisFromStars::Finalize; fSourceCam = "
529 // << fSourceCam << endl;
530
531 return kTRUE;
532}
533
534// --------------------------------------------------------------------------
535//
536// Creates a new canvas and draws the histograms.
537// Be careful: The histograms belongs to this object and won't get deleted
538// together with the canvas.
539//
540
541void MHTelAxisFromStars::Draw(Option_t *opt)
542{
543 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
544 pad->SetBorderMode(0);
545
546 AppendPad("");
547
548 //TCanvas *pad = new TCanvas("TelAxisFromStars", "TelAxis plots", 1200, 900);
549 //gROOT->SetSelectedPad(NULL);
550
551 pad->Divide(4,3);
552
553 pad->cd(1);
554 gPad->SetBorderMode(0);
555 fNStars->Draw(opt);
556
557 pad->cd(2);
558 gPad->SetBorderMode(0);
559 fNdoF->Draw(opt);
560
561 pad->cd(3);
562 gPad->SetBorderMode(0);
563 fLog10Chi2->Draw(opt);
564
565 pad->cd(4);
566 gPad->SetBorderMode(0);
567 fChi2Prob->Draw(opt);
568
569 pad->cd(5);
570 gPad->SetBorderMode(0);
571 fNumIter->Draw(opt);
572
573 pad->cd(6);
574 gPad->SetBorderMode(0);
575 fLambda->Draw(opt);
576
577 pad->cd(7);
578 gPad->SetBorderMode(0);
579 fAlfa->Draw(opt);
580
581 pad->cd(8);
582 gPad->SetBorderMode(0);
583 //SetColors();
584 //fShift->Draw("colz");
585 fShift->Draw("");
586
587
588 //-----------------------------------------------
589 // plot the expected positions of some sources
590 *fLog << "fSourcsCam = " << fSourceCam << endl;
591
592 /*
593 // copy names from MStarCam to the histograms
594 if (fSourceCam)
595 {
596 *fLog << "Source container exists " << endl;
597
598 MStarPos *starSource = 0;
599 TIter nextSource(fSourceCam->GetList());
600
601 *fLog << "&nextSource = " << &nextSource << endl;
602
603 starSource = (MStarPos*)nextSource();
604 if (starSource)
605 {
606 *fLog << "name of source container1 : " << starSource->GetName()
607 << endl;
608
609 fEstPos1->SetName(starSource->GetName());
610
611 starSource = (MStarPos*)nextSource();
612 if (starSource)
613 {
614 *fLog << "name of source container2 : " << starSource->GetName()
615 << endl;
616
617
618 fEstPos2->SetName(starSource->GetName());
619
620 starSource = (MStarPos*)nextSource();
621 if (starSource)
622 {
623
624 *fLog << "name of source container3 : " << starSource->GetName()
625 << endl;
626
627 fEstPos3->SetName(starSource->GetName());
628 }
629 }
630 }
631 }
632 */
633
634
635 if (fSourceCam)
636 {
637 *fLog << "MHTelAxisFromSrars::Draw; plotting" << endl;
638
639 pad->cd(9);
640 gPad->SetBorderMode(0);
641 //SetColors();
642 fEstPos1->Draw("");
643
644 pad->cd(10);
645 gPad->SetBorderMode(0);
646 //SetColors();
647 fEstPos2->Draw("");
648
649 pad->cd(11);
650 gPad->SetBorderMode(0);
651 //SetColors();
652 fEstPos3->Draw("");
653 }
654
655 pad->Modified();
656 pad->Update();
657}
658
659//---------------------------------------------------------------------
660//
661
662TH1 *MHTelAxisFromStars::GetHistByName(const TString name)
663{
664 if (name.Contains("NStars", TString::kIgnoreCase))
665 return fNStars;
666 if (name.Contains("NdoF", TString::kIgnoreCase))
667 return fNdoF;
668 if (name.Contains("Log10Chi2", TString::kIgnoreCase))
669 return fLog10Chi2;
670 if (name.Contains("Chi2Prob", TString::kIgnoreCase))
671 return fChi2Prob;
672 if (name.Contains("NumIter", TString::kIgnoreCase))
673 return fNumIter;
674 if (name.Contains("Lambda", TString::kIgnoreCase))
675 return fLambda;
676 if (name.Contains("Alfa", TString::kIgnoreCase))
677 return fAlfa;
678
679 if (name.Contains("Shift", TString::kIgnoreCase))
680 return fShift;
681 if (name.Contains("EstPos1", TString::kIgnoreCase))
682 return fEstPos1;
683 if (name.Contains("EstPos2", TString::kIgnoreCase))
684 return fEstPos2;
685 if (name.Contains("EstPos3", TString::kIgnoreCase))
686 return fEstPos3;
687
688 return NULL;
689}
690
691// --------------------------------------------------------------------------
692//
693// Setup an inversed deep blue sea palette.
694//
695
696void MHTelAxisFromStars::SetColors() const
697{
698 gStyle->SetPalette(51, NULL);
699 Int_t c[50];
700 for (int i=0; i<50; i++)
701 c[49-i] = gStyle->GetColorPalette(i);
702 gStyle->SetPalette(50, c);
703}
704
705
706//---------------------------------------------------------------------
707//
708
709void MHTelAxisFromStars::Paint(Option_t *opt)
710{
711 SetColors();
712 MH::Paint();
713}
714
715//==========================================================================
716
717
718
719
720
721
722
723
Note: See TracBrowser for help on using the repository browser.