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

Last change on this file since 5963 was 4799, checked in by wittek, 20 years ago
*** empty log message ***
File size: 17.5 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 {
449 if (fInputType == 1)
450 {
451 x = scale * starSource->GetMeanXFit();
452 y = scale * starSource->GetMeanYFit();
453 }
454 else
455 {
456 x = scale * starSource->GetMeanXCalc();
457 y = scale * starSource->GetMeanYCalc();
458 }
459
460 if (x>lowx && x<higx && y>lowy && y<higy)
461 {
462 if ( fNstarnames > 0 && starSource->GetName() == fStarnames[0] )
463 fEstPos1->Fill(x, y, w);
464
465 else if( fNstarnames > 1 && starSource->GetName() == fStarnames[1] )
466 fEstPos2->Fill(x, y, w);
467
468 else if( fNstarnames > 2 && starSource->GetName() == fStarnames[2] )
469 fEstPos3->Fill(x, y, w);
470 }
471 }
472 }
473
474 return kTRUE;
475}
476
477// --------------------------------------------------------------------------
478//
479// Finalize :
480//
481// it is called in the postprocessing
482// reset pointers in order to allow cloning of the object
483//
484Bool_t MHTelAxisFromStars::Finalize()
485{
486 //*fLog << "MHTelAxisFromStars::Finalize; fSourceCam = "
487 // << fSourceCam << endl;
488
489 return kTRUE;
490}
491
492// --------------------------------------------------------------------------
493//
494// Creates a new canvas and draws the histograms.
495// Be careful: The histograms belongs to this object and won't get deleted
496// together with the canvas.
497//
498
499void MHTelAxisFromStars::Draw(Option_t *opt)
500{
501 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
502 pad->SetBorderMode(0);
503
504 AppendPad("");
505
506 //TCanvas *pad = new TCanvas("TelAxisFromStars", "TelAxis plots", 1200, 900);
507 //gROOT->SetSelectedPad(NULL);
508
509 pad->Divide(4,3);
510
511 pad->cd(1);
512 gPad->SetBorderMode(0);
513 fNStars->Draw(opt);
514
515 pad->cd(2);
516 gPad->SetBorderMode(0);
517 fNdoF->Draw(opt);
518
519 pad->cd(3);
520 gPad->SetBorderMode(0);
521 fLog10Chi2->Draw(opt);
522
523 pad->cd(4);
524 gPad->SetBorderMode(0);
525 fChi2Prob->Draw(opt);
526
527 pad->cd(5);
528 gPad->SetBorderMode(0);
529 fNumIter->Draw(opt);
530
531 pad->cd(6);
532 gPad->SetBorderMode(0);
533 fLambda->Draw(opt);
534
535 pad->cd(7);
536 gPad->SetBorderMode(0);
537 fAlfa->Draw(opt);
538
539 pad->cd(8);
540 gPad->SetBorderMode(0);
541 //SetColors();
542 //fShift->Draw("colz");
543 fShift->Draw("");
544
545
546 //-----------------------------------------------
547 // plot the expected positions of some sources
548 //*fLog << "fSourcsCam = " << fSourceCam << endl;
549
550 if (fSourceCam)
551 {
552 *fLog << "MHTelAxisFromSrars::Draw; plotting" << endl;
553
554 pad->cd(9);
555 gPad->SetBorderMode(0);
556 //SetColors();
557 fEstPos1->Draw("");
558
559 pad->cd(10);
560 gPad->SetBorderMode(0);
561 //SetColors();
562 fEstPos2->Draw("");
563
564 pad->cd(11);
565 gPad->SetBorderMode(0);
566 //SetColors();
567 fEstPos3->Draw("");
568 }
569
570 pad->Modified();
571 pad->Update();
572}
573
574//---------------------------------------------------------------------
575//
576
577TH1 *MHTelAxisFromStars::GetHistByName(const TString name)
578{
579 if (name.Contains("NStars", TString::kIgnoreCase))
580 return fNStars;
581 if (name.Contains("NdoF", TString::kIgnoreCase))
582 return fNdoF;
583 if (name.Contains("Log10Chi2", TString::kIgnoreCase))
584 return fLog10Chi2;
585 if (name.Contains("Chi2Prob", TString::kIgnoreCase))
586 return fChi2Prob;
587 if (name.Contains("NumIter", TString::kIgnoreCase))
588 return fNumIter;
589 if (name.Contains("Lambda", TString::kIgnoreCase))
590 return fLambda;
591 if (name.Contains("Alfa", TString::kIgnoreCase))
592 return fAlfa;
593
594 if (name.Contains("Shift", TString::kIgnoreCase))
595 return fShift;
596 if (name.Contains("EstPos1", TString::kIgnoreCase))
597 return fEstPos1;
598 if (name.Contains("EstPos2", TString::kIgnoreCase))
599 return fEstPos2;
600 if (name.Contains("EstPos3", TString::kIgnoreCase))
601 return fEstPos3;
602
603 return NULL;
604}
605
606// --------------------------------------------------------------------------
607//
608// Setup an inversed deep blue sea palette.
609//
610
611void MHTelAxisFromStars::SetColors() const
612{
613 gStyle->SetPalette(51, NULL);
614 Int_t c[50];
615 for (int i=0; i<50; i++)
616 c[49-i] = gStyle->GetColorPalette(i);
617 gStyle->SetPalette(50, c);
618}
619
620
621//---------------------------------------------------------------------
622//
623
624void MHTelAxisFromStars::Paint(Option_t *opt)
625{
626 SetColors();
627 MH::Paint();
628}
629
630//==========================================================================
631
632
633
634
635
636
637
638
Note: See TracBrowser for help on using the repository browser.