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

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