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

Last change on this file since 4682 was 4560, checked in by wittek, 21 years ago
*** empty log message ***
File size: 15.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!
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 "MStarLocalCam.h"
51#include "MStarLocalPos.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 fChi2Prob = new TH1D("Chi2-Prob", "Chi2 probability", 40, 0.0, 1.0);
81 fChi2Prob->SetDirectory(NULL);
82 fChi2Prob->SetXTitle("Chi2 probability");
83 fChi2Prob->SetYTitle("Counts");
84
85
86 fNumIter = new TH1D("NumIter", "Number of iterations", 50, -0.5, 49.5);
87 fNumIter->SetDirectory(NULL);
88 fNumIter->SetXTitle("Number of iterations");
89 fNumIter->SetYTitle("Counts");
90
91
92 fLambda = new TH1D("Lambda", "Scale factor lambda", 80, 0.90, 1.10);
93 fLambda->SetDirectory(NULL);
94 fLambda->SetXTitle("Scale factor lambda");
95 fLambda->SetYTitle("Counts");
96
97
98 fAlfa = new TH1D("Alfa", "Rotation angle alfa", 100, -2.5, 2.5);
99 fAlfa->SetDirectory(NULL);
100 fAlfa->SetXTitle("Rotation angle alfa [\\circ]");
101 fAlfa->SetYTitle("Counts");
102
103
104 fShift = new TH2D("Shift", "Sky-Cam transformnation : (x,y) shift",
105 51, -445, 445, 51, -445, 445);
106 fShift->SetDirectory(NULL);
107 fShift->SetZTitle("Counts");
108 fShift->SetXTitle("x-shift [\\circ]");
109 fShift->SetYTitle("y-shift [\\circ]");
110
111
112 fEstPos1 = new TH2D("EstPos1", "Estimated position1",
113 51, -445, 445, 51, -445, 445);
114 fEstPos1->SetDirectory(NULL);
115 fEstPos1->SetZTitle("Counts");
116 fEstPos1->SetXTitle("Estimated position1 x [\\circ]");
117 fEstPos1->SetYTitle("Estimated position1 y [\\circ]");
118
119 fEstPos2 = new TH2D("EstPos2", "Estimated position2",
120 51, -445, 445, 51, -445, 445);
121 fEstPos2->SetDirectory(NULL);
122 fEstPos2->SetZTitle("Counts");
123 fEstPos2->SetXTitle("Estimated position2 x [\\circ]");
124 fEstPos2->SetYTitle("Estimated position2 y [\\circ]");
125
126 // default input type : results from correlated Gauss fit
127 fInputType = 2;
128}
129
130// --------------------------------------------------------------------------
131//
132// Set the type of the input
133//
134// type = 0 calculated star positions (by averaging)
135// type = 1 fitted star positions (by uncorrelated Gauss fit)
136// type = 2 fitted star positions (by correlated Gauss fit)
137//
138void MHTelAxisFromStars::SetInputType(Int_t type)
139{
140 *fLog << all << "MHTelAxisFromStars::SetInputType; input type is set equal to : "
141 << type << endl;
142
143 fInputType = type;
144}
145
146
147// --------------------------------------------------------------------------
148//
149// Delete the histograms
150//
151MHTelAxisFromStars::~MHTelAxisFromStars()
152{
153 delete fNStars;
154 delete fNdoF;
155 delete fChi2Prob;
156 delete fNumIter;
157 delete fLambda;
158 delete fAlfa;
159
160 delete fShift;
161 delete fEstPos1;
162 delete fEstPos2;
163}
164
165// --------------------------------------------------------------------------
166//
167// Setup the Binning for the histograms.
168//
169// Find the pointers to the containers.
170//
171// Use this function if you want to set the conversion factor which
172// is used to convert the mm-scale in the camera plane into the deg-scale.
173// The conversion factor is part of the camera geometry. Please create a
174// corresponding MGeomCam container.
175//
176Bool_t MHTelAxisFromStars::SetupFill(const MParList *plist)
177{
178
179 fStarLocalCam = (MStarLocalCam*)plist->FindObject("MStarLocalCam", "MStarLocalCam");
180 if (!fStarLocalCam)
181 {
182 *fLog << err << "MHTelAxisFromStars::SetupFill; container 'MStarLocalCam' not found... aborting." << endl;
183 return kFALSE;
184 }
185
186
187 fSourceLocalCam = (MStarLocalCam*)plist->FindObject("MSourceLocalCam", "MStarLocalCam");
188 if (!fSourceLocalCam)
189 {
190 *fLog << err << "MHTelAxisFromStars::SetupFill; container 'MSourceLocalCam' not found... aborting." << endl;
191 return kFALSE;
192 }
193
194
195
196 fSrcPos = (MSrcPosCam*)plist->FindObject(AddSerialNumber("MSrcPosCam"));
197 if (!fSrcPos)
198 {
199 *fLog << err << "MHTelAxisFromStars::SetupFill; MSrcPosCam not found... aborting" << endl;
200 return kFALSE;
201 }
202
203 fSkyCamTrans = (MSkyCamTrans*)plist->FindObject(AddSerialNumber("MSkyCamTrans"));
204 if (!fSkyCamTrans)
205 {
206 *fLog << err << "MHTelAxisFromStars::SetupFill; MSkyCamTrans not found... aborting" << endl;
207 return kFALSE;
208 }
209
210
211 //------------------------------------------------------------------
212 const MGeomCam *geom = (MGeomCam*)plist->FindObject("MGeomCam");
213 if (!geom)
214 {
215 *fLog << warn << GetDescriptor()
216 << ": No Camera Geometry available. Using mm-scale for histograms."
217 << endl;
218 SetMmScale(kTRUE);
219 }
220 else
221 {
222 fMm2Deg = geom->GetConvMm2Deg();
223 SetMmScale(kFALSE);
224 }
225
226 ApplyBinning(*plist, "NStars", fNStars);
227 ApplyBinning(*plist, "NdoF", fNdoF);
228 ApplyBinning(*plist, "Chi2Prob", fChi2Prob);
229 ApplyBinning(*plist, "NumIter", fNumIter);
230 ApplyBinning(*plist, "Lambda", fLambda);
231 ApplyBinning(*plist, "Alfa", fAlfa);
232
233 const MBinning *bins = (MBinning*)plist->FindObject("BinningCamera");
234 if (!bins)
235 {
236 float r = geom ? geom->GetMaxRadius() : 600;
237 r *= 0.9;
238 if (!fUseMmScale)
239 r *= fMm2Deg;
240
241 MBinning b;
242 b.SetEdges(61, -r, r);
243 SetBinning(fShift, &b, &b);
244 SetBinning(fEstPos1, &b, &b);
245 SetBinning(fEstPos2, &b, &b);
246 }
247 else
248 {
249 SetBinning(fShift, bins, bins);
250 SetBinning(fEstPos1, bins, bins);
251 SetBinning(fEstPos2, bins, bins);
252 }
253
254 //-----------------------------------------------
255 // copy names from MStarLocalCam to the histograms
256 MStarLocalPos *starSource = 0;
257 TIter nextSource(fSourceLocalCam->GetList());
258
259 starSource = (MStarLocalPos*)nextSource();
260 if (starSource)
261 {
262 fEstPos1->SetName(starSource->GetName());
263
264 starSource = (MStarLocalPos*)nextSource();
265 if (starSource)
266 fEstPos2->SetName(starSource->GetName());
267 }
268
269 return kTRUE;
270}
271
272// --------------------------------------------------------------------------
273//
274// Use this function to setup your own conversion factor between degrees
275// and millimeters. The conversion factor should be the one calculated in
276// MGeomCam. Use this function with Caution: You could create wrong values
277// by setting up your own scale factor.
278//
279void MHTelAxisFromStars::SetMm2Deg(Float_t mmdeg)
280{
281 if (mmdeg<0)
282 {
283 *fLog << warn << dbginf
284 << "Warning - Conversion factor < 0 - nonsense. Ignored."
285 << endl;
286 return;
287 }
288
289 if (fMm2Deg>=0)
290 *fLog << warn << dbginf
291 << "Warning - Conversion factor already set. Overwriting"
292 << endl;
293
294 fMm2Deg = mmdeg;
295}
296
297// --------------------------------------------------------------------------
298//
299// With this function you can convert the histogram ('on the fly') between
300// degrees and millimeters.
301//
302void MHTelAxisFromStars::SetMmScale(Bool_t mmscale)
303{
304 if (fUseMmScale == mmscale)
305 return;
306
307 if (fMm2Deg<0)
308 {
309 *fLog << warn << dbginf
310 << "Warning - Sorry, no conversion factor for conversion available."
311 << endl;
312 return;
313 }
314
315 const Double_t scale = mmscale ? 1./fMm2Deg : fMm2Deg;
316 MH::ScaleAxis(fShift, scale, scale);
317 MH::ScaleAxis(fEstPos1, scale, scale);
318 MH::ScaleAxis(fEstPos2, scale, scale);
319
320 if (mmscale)
321 {
322 fShift->SetXTitle("x-shift [mm]");
323 fShift->SetYTitle("y-shift [mm]");
324
325 fEstPos1->SetXTitle("Estimated position1 x [mm]");
326 fEstPos1->SetYTitle("Estimated position1 y [mm]");
327
328 fEstPos2->SetXTitle("Estimated position2 x [mm]");
329 fEstPos2->SetYTitle("Estimated position2 y [mm]");
330 }
331 else
332 {
333 fShift->SetXTitle("x-shift [\\circ]");
334 fShift->SetYTitle("y-shift [\\circ]");
335
336 fEstPos1->SetXTitle("Estimated position1 x [\\circ]");
337 fEstPos1->SetYTitle("Estimated position1 y [\\circ]");
338
339 fEstPos2->SetXTitle("Estimated position2 x [\\circ]");
340 fEstPos2->SetYTitle("Estimated position2 y [\\circ]");
341 }
342
343 fUseMmScale = mmscale;
344}
345
346// --------------------------------------------------------------------------
347//
348// Fill the histograms
349//
350//
351Bool_t MHTelAxisFromStars::Fill(const MParContainer *par, const Stat_t w)
352{
353 Int_t Ndof = fSkyCamTrans->GetNdof();
354 if (Ndof < 0)
355 return kTRUE;
356
357 const Double_t scale = fUseMmScale ? 1 : fMm2Deg;
358
359
360 fNStars ->Fill(fSkyCamTrans->GetNStars(), w);
361 fNdoF ->Fill(fSkyCamTrans->GetNdof(), w);
362 fChi2Prob ->Fill(fSkyCamTrans->GetChiSquareProb(), w);
363 fNumIter ->Fill(fSkyCamTrans->GetNumIter(), w);
364 fLambda ->Fill(fSkyCamTrans->GetLambda(), w);
365 fAlfa ->Fill(fSkyCamTrans->GetAlfa(), w);
366
367
368 Double_t lowx;
369 Double_t higx;
370 Double_t lowy;
371 Double_t higy;
372 Double_t x;
373 Double_t y;
374 Int_t nbix;
375 Int_t nbiy;
376
377 nbix = fShift->GetXaxis()->GetNbins();
378 lowx = fShift->GetXaxis()->GetBinLowEdge(1);
379 higx = fShift->GetXaxis()->GetBinLowEdge(nbix+1);
380
381 nbiy = fShift->GetYaxis()->GetNbins();
382 lowy = fShift->GetYaxis()->GetBinLowEdge(1);
383 higy = fShift->GetYaxis()->GetBinLowEdge(nbiy+1);
384
385 x = scale * (fSkyCamTrans->GetShiftD())[0];
386 y = scale * (fSkyCamTrans->GetShiftD())[1];
387 if (x>lowx && x<higx && y>lowy && y<higy)
388 {
389 fShift ->Fill(x, y, w);
390 }
391
392
393 MStarLocalPos *starSource = 0;
394 TIter nextSource(fSourceLocalCam->GetList());
395
396 starSource = (MStarLocalPos*)nextSource();
397 if (starSource)
398 {
399 if (fInputType == 2)
400 {
401 x = scale * starSource->GetMeanXCGFit();
402 y = scale * starSource->GetMeanYCGFit();
403 }
404 else if (fInputType == 1)
405 {
406 x = scale * starSource->GetMeanXFit();
407 y = scale * starSource->GetMeanYFit();
408 }
409 else
410 {
411 x = scale * starSource->GetMeanXCalc();
412 y = scale * starSource->GetMeanYCalc();
413 }
414
415 if (x>lowx && x<higx && y>lowy && y<higy)
416 {
417 fEstPos1 ->Fill(x, y, w);
418 }
419
420
421 starSource = (MStarLocalPos*)nextSource();
422 if (starSource)
423 {
424 if (fInputType == 2)
425 {
426 x = scale * starSource->GetMeanXCGFit();
427 y = scale * starSource->GetMeanYCGFit();
428 }
429 else if (fInputType == 1)
430 {
431 x = scale * starSource->GetMeanXFit();
432 y = scale * starSource->GetMeanYFit();
433 }
434 else
435 {
436 x = scale * starSource->GetMeanXCalc();
437 y = scale * starSource->GetMeanYCalc();
438 }
439
440 if (x>lowx && x<higx && y>lowy && y<higy)
441 {
442 fEstPos2 ->Fill(x, y, w);
443 }
444 }
445 }
446
447 return kTRUE;
448}
449
450// --------------------------------------------------------------------------
451//
452// Finalize :
453//
454// it is called in the postprocessing
455// reset pointers in order to allow cloning of the object
456//
457Bool_t MHTelAxisFromStars::Finalize()
458{
459
460 return kTRUE;
461}
462
463// --------------------------------------------------------------------------
464//
465// Creates a new canvas and draws the histograms.
466// Be careful: The histograms belongs to this object and won't get deleted
467// together with the canvas.
468//
469
470void MHTelAxisFromStars::Draw(Option_t *opt)
471{
472 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
473 pad->SetBorderMode(0);
474
475 AppendPad("");
476
477 //TCanvas *pad = new TCanvas("TelAxisFromStars", "TelAxis plots", 900, 900);
478 //gROOT->SetSelectedPad(NULL);
479
480 pad->Divide(3,3);
481
482 pad->cd(1);
483 gPad->SetBorderMode(0);
484 fNStars->Draw(opt);
485
486 pad->cd(2);
487 gPad->SetBorderMode(0);
488 fNdoF->Draw(opt);
489
490 pad->cd(3);
491 gPad->SetBorderMode(0);
492 fChi2Prob->Draw(opt);
493
494 pad->cd(4);
495 gPad->SetBorderMode(0);
496 fNumIter->Draw(opt);
497
498 pad->cd(5);
499 gPad->SetBorderMode(0);
500 fLambda->Draw(opt);
501
502 pad->cd(6);
503 gPad->SetBorderMode(0);
504 fAlfa->Draw(opt);
505
506 pad->cd(7);
507 gPad->SetBorderMode(0);
508 //SetColors();
509 //fShift->Draw("colz");
510 fShift->Draw("");
511
512 pad->cd(8);
513 gPad->SetBorderMode(0);
514 //SetColors();
515 fEstPos1->Draw("");
516
517 pad->cd(9);
518 gPad->SetBorderMode(0);
519 //SetColors();
520 fEstPos2->Draw("");
521
522 pad->Modified();
523 pad->Update();
524}
525
526//---------------------------------------------------------------------
527//
528
529TH1 *MHTelAxisFromStars::GetHistByName(const TString name)
530{
531 if (name.Contains("NStars", TString::kIgnoreCase))
532 return fNStars;
533 if (name.Contains("NdoF", TString::kIgnoreCase))
534 return fNdoF;
535 if (name.Contains("Chi2Prob", TString::kIgnoreCase))
536 return fChi2Prob;
537 if (name.Contains("NumIter", TString::kIgnoreCase))
538 return fNumIter;
539 if (name.Contains("Lambda", TString::kIgnoreCase))
540 return fLambda;
541 if (name.Contains("Alfa", TString::kIgnoreCase))
542 return fAlfa;
543
544 if (name.Contains("Shift", TString::kIgnoreCase))
545 return fShift;
546 if (name.Contains("EstPos1", TString::kIgnoreCase))
547 return fEstPos1;
548 if (name.Contains("EstPos2", TString::kIgnoreCase))
549 return fEstPos2;
550
551 return NULL;
552}
553
554// --------------------------------------------------------------------------
555//
556// Setup an inversed deep blue sea palette.
557//
558
559void MHTelAxisFromStars::SetColors() const
560{
561 gStyle->SetPalette(51, NULL);
562 Int_t c[50];
563 for (int i=0; i<50; i++)
564 c[49-i] = gStyle->GetColorPalette(i);
565 gStyle->SetPalette(50, c);
566}
567
568
569//---------------------------------------------------------------------
570//
571
572void MHTelAxisFromStars::Paint(Option_t *opt)
573{
574 SetColors();
575 MH::Paint();
576}
577
578//==========================================================================
579
580
581
582
583
584
585
586
Note: See TracBrowser for help on using the repository browser.