source: trunk/MagicSoft/Mars/mimage/MHNewImagePar.cc@ 8321

Last change on this file since 8321 was 7687, checked in by tbretz, 20 years ago
*** empty log message ***
File size: 13.8 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, 03/2003 <mailto:wittek@mppmu.mpg.de>
19! Author(s): Thomas Bretz, 04/2003 <mailto:tbretz@astro.uni-wuerzburg.de>
20!
21! Copyright: MAGIC Software Development, 2000-2006
22!
23!
24\* ======================================================================== */
25
26/////////////////////////////////////////////////////////////////////////////
27//
28// MHNewImagePar
29//
30//
31// Version 2:
32// ----------
33// + TH1F fHistConcCOG; // [ratio] concentration around the center of gravity (all pixels)
34// + TH1F fHistConcCore; // [ratio] concentration of signals inside ellipse (used pixels)
35//
36////////////////////////////////////////////////////////////////////////////
37#include "MHNewImagePar.h"
38
39#include <math.h>
40
41#include <TH1.h>
42#include <TPad.h>
43#include <TCanvas.h>
44
45#include "MLog.h"
46#include "MLogManip.h"
47
48#include "MGeomCam.h"
49#include "MBinning.h"
50#include "MParList.h"
51
52#include "MHillas.h"
53#include "MNewImagePar.h"
54
55ClassImp(MHNewImagePar);
56
57using namespace std;
58
59// --------------------------------------------------------------------------
60//
61// Setup histograms
62//
63MHNewImagePar::MHNewImagePar(const char *name, const char *title)
64 : fMm2Deg(1), fUseMmScale(kTRUE)
65{
66 fName = name ? name : "MHNewImagePar";
67 fTitle = title ? title : "Histograms of new image parameters";
68
69 fHistLeakage1.SetName("Leakage1");
70 fHistLeakage1.SetTitle("Leakage_{1}");
71 fHistLeakage1.SetXTitle("Leakage");
72 fHistLeakage1.SetYTitle("Counts");
73 fHistLeakage1.SetDirectory(NULL);
74 fHistLeakage1.UseCurrentStyle();
75 fHistLeakage1.SetFillStyle(4000);
76
77 fHistLeakage2.SetName("Leakage2");
78 fHistLeakage2.SetTitle("Leakage_{2}");
79 fHistLeakage2.SetXTitle("Leakage");
80 fHistLeakage2.SetYTitle("Counts");
81 fHistLeakage2.SetDirectory(NULL);
82 fHistLeakage2.UseCurrentStyle();
83 fHistLeakage2.SetLineColor(kBlue);
84 fHistLeakage2.SetFillStyle(4000);
85
86 fHistUsedPix.SetName("UsedPix");
87 fHistUsedPix.SetTitle("Number of used pixels");
88 fHistUsedPix.SetXTitle("Number of Pixels");
89 fHistUsedPix.SetYTitle("Counts");
90 fHistUsedPix.SetDirectory(NULL);
91 fHistUsedPix.UseCurrentStyle();
92 fHistUsedPix.SetLineColor(kBlue);
93 fHistUsedPix.SetFillStyle(4000);
94
95 fHistCorePix.SetName("CorePix");
96 fHistCorePix.SetTitle("Number of core pixels");
97 fHistCorePix.SetXTitle("Number of Pixels");
98 fHistCorePix.SetYTitle("Counts");
99 fHistCorePix.SetDirectory(NULL);
100 fHistCorePix.UseCurrentStyle();
101 fHistCorePix.SetLineColor(kBlack);
102 fHistCorePix.SetFillStyle(4000);
103
104 fHistUsedArea.SetName("UsedArea");
105 fHistUsedArea.SetTitle("Area of used pixels");
106 fHistUsedArea.SetXTitle("Area [m^2]");
107 fHistUsedArea.SetYTitle("Counts");
108 fHistUsedArea.SetDirectory(NULL);
109 fHistUsedArea.UseCurrentStyle();
110 fHistUsedArea.SetLineColor(kBlue);
111 fHistUsedArea.SetFillStyle(4000);
112
113 fHistCoreArea.SetName("CoreArea");
114 fHistCoreArea.SetTitle("Area of core pixels");
115 fHistCoreArea.SetXTitle("Area [m^2]");
116 fHistCoreArea.SetYTitle("Counts");
117 fHistCoreArea.SetDirectory(NULL);
118 fHistCoreArea.UseCurrentStyle();
119 fHistCoreArea.SetLineColor(kBlack);
120 fHistCoreArea.SetFillStyle(4000);
121
122 fHistConc.SetDirectory(NULL);
123 fHistConc1.SetDirectory(NULL);
124 fHistConcCOG.SetDirectory(NULL);
125 fHistConcCore.SetDirectory(NULL);
126 fHistConc.SetName("Conc2");
127 fHistConc1.SetName("Conc1");
128 fHistConcCOG.SetName("ConcCOG");
129 fHistConcCore.SetName("ConcCore");
130 fHistConc.SetTitle("Ratio: Conc");
131 fHistConc1.SetTitle("Ratio: Conc1");
132 fHistConcCOG.SetTitle("Ratio: ConcCOG");
133 fHistConcCore.SetTitle("Ratio: ConcCore");
134 fHistConc.SetXTitle("Ratio");
135 fHistConc1.SetXTitle("Ratio");
136 fHistConcCOG.SetXTitle("Ratio");
137 fHistConcCore.SetXTitle("Ratio");
138 fHistConc.SetYTitle("Counts");
139 fHistConc1.SetYTitle("Counts");
140 fHistConcCOG.SetYTitle("Counts");
141 fHistConcCore.SetYTitle("Counts");
142 fHistConc.UseCurrentStyle();
143 fHistConc1.UseCurrentStyle();
144 fHistConcCOG.UseCurrentStyle();
145 fHistConcCore.UseCurrentStyle();
146 fHistConc.SetFillStyle(4000);
147 fHistConc1.SetFillStyle(4000);
148 fHistConcCOG.SetFillStyle(4000);
149 fHistConcCore.SetFillStyle(4000);
150 fHistConc1.SetLineColor(kBlue);
151 fHistConcCOG.SetLineColor(kBlue);
152
153 MBinning bins;
154
155 bins.SetEdges(100, 0, 1);
156 bins.Apply(fHistLeakage1);
157 bins.Apply(fHistLeakage2);
158 bins.Apply(fHistConc);
159 bins.Apply(fHistConc1);
160 bins.Apply(fHistConcCOG);
161 bins.Apply(fHistConcCore);
162
163 bins.SetEdges(75, 0.5, 150.5);
164 bins.Apply(fHistUsedPix);
165 bins.Apply(fHistCorePix);
166
167 //bins.SetEdges(75, 0, 0.249);
168 //bins.Apply(fHistUsedArea);
169 //bins.Apply(fHistCoreArea);
170}
171
172// --------------------------------------------------------------------------
173//
174// Setup the Binning for the histograms automatically if the correct
175// instances of MBinning
176//
177Bool_t MHNewImagePar::SetupFill(const MParList *plist)
178{
179 MGeomCam *geom = (MGeomCam*)plist->FindObject("MGeomCam");
180 if (!geom)
181 *fLog << warn << GetDescriptor() << ": No Camera Geometry available. Using mm-scale for histograms." << endl;
182 else
183 {
184 fMm2Deg = geom->GetConvMm2Deg();
185 SetMmScale(kFALSE);
186 }
187
188 const MBinning *bins = (MBinning*)plist->FindObject("BinningArea");
189 if (!bins)
190 {
191 float r = geom ? 1.5 : 0.133;
192
193 MBinning b;
194 b.SetEdges(50, 0, r);
195 b.Apply(fHistUsedArea);
196 b.Apply(fHistCoreArea);
197 }
198 else
199 {
200 bins->Apply(fHistUsedArea);
201 bins->Apply(fHistCoreArea);
202 }
203
204 ApplyBinning(*plist, "Leakage", &fHistLeakage1);
205 ApplyBinning(*plist, "Leakage", &fHistLeakage2);
206
207 ApplyBinning(*plist, "Pixels", &fHistUsedPix);
208 ApplyBinning(*plist, "Pixels", &fHistCorePix);
209
210 //ApplyBinning(*plist, "Area", &fHistUsedArea);
211 //ApplyBinning(*plist, "Area", &fHistCoreArea);
212
213 ApplyBinning(*plist, "Conc", &fHistConc);
214 ApplyBinning(*plist, "Conc1", &fHistConc1);
215 ApplyBinning(*plist, "ConcCOG", &fHistConcCOG);
216 ApplyBinning(*plist, "ConcCore", &fHistConcCore);
217
218 return kTRUE;
219}
220
221
222// --------------------------------------------------------------------------
223//
224// Fill the histograms with data from a MNewImagePar container.
225//
226Bool_t MHNewImagePar::Fill(const MParContainer *par, const Stat_t w)
227{
228 if (!par)
229 {
230 *fLog << err << "MHNewImagePar::Fill: Pointer (!=NULL) expected." << endl;
231 return kFALSE;
232 }
233
234 const Double_t scale = fUseMmScale ? 1e-6 : fMm2Deg*fMm2Deg;
235
236 const MNewImagePar &h = *(MNewImagePar*)par;
237
238 fHistLeakage1.Fill(h.GetLeakage1(), w);
239 fHistLeakage2.Fill(h.GetLeakage2(), w);
240
241 fHistUsedPix.Fill(h.GetNumUsedPixels(), w);
242 fHistCorePix.Fill(h.GetNumCorePixels(), w);
243
244 fHistUsedArea.Fill(h.GetUsedArea()*scale, w);
245 fHistCoreArea.Fill(h.GetCoreArea()*scale, w);
246
247 fHistConc.Fill(h.GetConc(), w);
248 fHistConc1.Fill(h.GetConc1(), w);
249 fHistConcCOG.Fill(h.GetConcCOG(), w);
250 fHistConcCore.Fill(h.GetConcCore(), w);
251
252 return kTRUE;
253}
254
255// --------------------------------------------------------------------------
256//
257// With this function you can convert the histogram ('on the fly') between
258// degrees and millimeters.
259//
260void MHNewImagePar::SetMmScale(Bool_t mmscale)
261{
262 if (fUseMmScale == mmscale)
263 return;
264
265 if (fMm2Deg<0)
266 {
267 *fLog << warn << dbginf << "Warning - Sorry, no conversion factor for conversion available." << endl;
268 return;
269 }
270
271 const Double_t scale = mmscale ? 1./(fMm2Deg*fMm2Deg) : (fMm2Deg*fMm2Deg);
272 MH::ScaleAxis(&fHistUsedArea, scale);
273 MH::ScaleAxis(&fHistCoreArea, scale);
274
275 if (mmscale)
276 {
277 fHistUsedArea.SetXTitle("A [m^{2}]");
278 fHistCoreArea.SetXTitle("A [m^{2}]");
279 }
280 else
281 {
282 fHistUsedArea.SetXTitle("A [deg^{2}]");
283 fHistCoreArea.SetXTitle("A [deg^{2}]");
284 }
285
286 fUseMmScale = mmscale;
287}
288
289// --------------------------------------------------------------------------
290//
291// Use this function to setup your own conversion factor between degrees
292// and millimeters. The conversion factor should be the one calculated in
293// MGeomCam. Use this function with Caution: You could create wrong values
294// by setting up your own scale factor.
295//
296void MHNewImagePar::SetMm2Deg(Float_t mmdeg)
297{
298 if (mmdeg<0)
299 {
300 *fLog << warn << dbginf << "Warning - Conversion factor < 0 - nonsense. Ignored." << endl;
301 return;
302 }
303
304 if (fMm2Deg>=0)
305 *fLog << warn << dbginf << "Warning - Conversion factor already set. Overwriting" << endl;
306
307 fMm2Deg = mmdeg;
308}
309
310void MHNewImagePar::Paint(Option_t *o)
311{
312 if (fHistLeakage1.GetMaximum()>0 && gPad->GetPad(1) && gPad->GetPad(1)->GetPad(1))
313 gPad->GetPad(1)->GetPad(1)->SetLogy();
314}
315
316// --------------------------------------------------------------------------
317//
318// Creates a new canvas and draws the two histograms into it.
319// Be careful: The histograms belongs to this object and won't get deleted
320// together with the canvas.
321//
322void MHNewImagePar::Draw(Option_t *o)
323{
324 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
325 pad->SetBorderMode(0);
326
327 AppendPad("");
328
329 // FIXME: If same-option given make two independant y-axis!
330 const TString opt(o);
331 const Bool_t same = opt.Contains("same");
332
333 if (!same)
334 pad->Divide(2, 1);
335 else
336 {
337 fHistLeakage1.SetName("Leakage1Same");
338 fHistLeakage2.SetName("Leakage2Same");
339 fHistUsedPix.SetName("UsedPixSame");
340 fHistCorePix.SetName("CorePixSame");
341 fHistUsedArea.SetName("UsedAreaSame");
342 fHistCoreArea.SetName("CoreAreaSame");
343 fHistConcCOG.SetName("ConcCOGSame");
344 fHistConcCore.SetName("ConcCoreSame");
345 fHistConc1.SetName("Conc1Same");
346 fHistConc.SetName("Conc2Same");
347
348 fHistLeakage1.SetDirectory(0);
349 fHistLeakage2.SetDirectory(0);
350 fHistUsedPix.SetDirectory(0);
351 fHistCorePix.SetDirectory(0);
352 fHistUsedArea.SetDirectory(0);
353 fHistCoreArea.SetDirectory(0);
354 fHistConcCOG.SetDirectory(0);
355 fHistConcCore.SetDirectory(0);
356 fHistConc1.SetDirectory(0);
357 fHistConc.SetDirectory(0);
358
359 fHistLeakage1.SetLineColor(kMagenta);
360 fHistLeakage1.SetLineColor(kCyan);
361 fHistCorePix.SetLineColor(kMagenta);
362 fHistUsedPix.SetLineColor(kCyan);
363 fHistConc1.SetLineColor(kMagenta);
364 fHistConc.SetLineColor(kCyan);
365 fHistConcCOG.SetLineColor(kMagenta);
366 fHistConcCore.SetLineColor(kCyan);
367 fHistCoreArea.SetLineColor(kMagenta);
368 fHistUsedArea.SetLineColor(kCyan);
369 }
370
371 pad->cd(1);
372 TVirtualPad *pad1=gPad;
373 pad1->SetBorderMode(0);
374 if (!same)
375 pad1->Divide(1,3, 0.001, 0.001);
376
377 pad1->cd(1);
378 gPad->SetBorderMode(0);
379 gPad->SetGridx();
380 TAxis &x = *fHistLeakage1.GetXaxis();
381 x.SetRangeUser(0.0, x.GetXmax());
382 RemoveFromPad("Leakage1Same");
383 RemoveFromPad("Leakage2Same");
384 MH::DrawSame(fHistLeakage1, fHistLeakage2, "Leakage1 and Leakage2", same);
385 fHistLeakage1.SetMinimum();
386 fHistLeakage2.SetMinimum();
387 fHistLeakage2.SetMaximum(0.1); // dummy value to allow log-scale
388
389 pad1->cd(2);
390 gPad->SetBorderMode(0);
391 gPad->SetGridx();
392 RemoveFromPad("UsedPixSame");
393 RemoveFromPad("CorePixSame");
394 MH::DrawSame(fHistCorePix, fHistUsedPix, "Number of core/used Pixels", same);
395
396 pad1->cd(3);
397 gPad->SetBorderMode(0);
398 gPad->SetGridx();
399 RemoveFromPad("CoreAreaSame");
400 RemoveFromPad("UsedAreaSame");
401 MH::DrawSame(fHistCoreArea, fHistUsedArea, "Area of core/used Pixels", same);
402
403 pad->cd(2);
404 TVirtualPad *pad2=gPad;
405 pad2->SetBorderMode(0);
406 if (!same)
407 pad2->Divide(1, 2, 0.001, 0.001);
408
409 pad2->cd(1);
410 gPad->SetBorderMode(0);
411 gPad->SetGridx();
412 RemoveFromPad("Conc1Same");
413 RemoveFromPad("Conc2Same");
414 MH::DrawSame(fHistConc1, fHistConc, "Concentrations", same);
415
416 pad2->cd(2);
417 gPad->SetBorderMode(0);
418 gPad->SetGridx();
419 RemoveFromPad("ConcCOGSame");
420 RemoveFromPad("ConcCoreSame");
421 MH::DrawSame(fHistConcCore, fHistConcCOG, "Concentrations", same);
422}
423
424TH1 *MHNewImagePar::GetHistByName(const TString name) const
425{
426 if (name.Contains("Leakage1", TString::kIgnoreCase))
427 return const_cast<TH1F*>(&fHistLeakage1);
428 if (name.Contains("Leakage2", TString::kIgnoreCase))
429 return const_cast<TH1F*>(&fHistLeakage2);
430 if (name.Contains("ConcCOG", TString::kIgnoreCase))
431 return const_cast<TH1F*>(&fHistConcCOG);
432 if (name.Contains("ConcCore", TString::kIgnoreCase))
433 return const_cast<TH1F*>(&fHistConcCore);
434 if (name.Contains("Conc1", TString::kIgnoreCase)) // must be first!
435 return const_cast<TH1F*>(&fHistConc1);
436 if (name.Contains("Conc", TString::kIgnoreCase))
437 return const_cast<TH1F*>(&fHistConc);
438 if (name.Contains("UsedPix", TString::kIgnoreCase))
439 return const_cast<TH1F*>(&fHistUsedPix);
440 if (name.Contains("CorePix", TString::kIgnoreCase))
441 return const_cast<TH1F*>(&fHistCorePix);
442 if (name.Contains("UsedArea", TString::kIgnoreCase))
443 return const_cast<TH1F*>(&fHistUsedArea);
444 if (name.Contains("CoreArea", TString::kIgnoreCase))
445 return const_cast<TH1F*>(&fHistCoreArea);
446
447 return NULL;
448}
Note: See TracBrowser for help on using the repository browser.