source: trunk/Mars/mhist/MHCamEvent.cc@ 19345

Last change on this file since 19345 was 19345, checked in by tbretz, 3 years ago
Improves the behaviour with RecursiveRemove. Strictly speaking this change might only be necessary if a class contains more than one member which is bound to recursive remove. Onth other hand setting all members to NULL which might be affected by RecursiveRemove is not wrong either.
File size: 12.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): Thomas Bretz, 12/2002 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2003
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MHCamEvent
28//
29// A histogram to store the sum of camera events. This can be photons,
30// currents or enything else derived from MCamEvent
31//
32// Setup
33// =====
34//
35// To count how often a certain pixel is above or below a threshold do:
36// MHCamEvent::SetThreshold(5.5); // Default=LowerBound
37// MHCamEvent::SetThreshold(5.5, MHCamEvent::kIsUpperBound);
38//
39//
40// Axis titles
41// ===========
42//
43// 1) If no other title is given 'a.u.' is used.
44// 2) If the title of MHCamEvent is different from the default,
45// it is used as histogram title. You can use this to set the
46// axis title, too. For more information see TH1::SetTitle, eg.
47// SetTitle("MyHist;;y[cm];Counts");
48// Make sure to set an empty x-axis title.
49//
50// For example:
51// MHCamEvent myhist("Titele;;y [cm]");
52//
53// Class Version 3:
54// ----------------
55// + fNameSub
56// + fIsRatio
57//
58/////////////////////////////////////////////////////////////////////////////
59#include "MHCamEvent.h"
60
61#include <TF1.h>
62#include <TPaveStats.h>
63
64#include "MLog.h"
65#include "MLogManip.h"
66
67#include "MString.h"
68#include "MParList.h"
69#include "MCamEvent.h"
70#include "MHCamera.h"
71
72#include "MGeomCam.h"
73
74ClassImp(MHCamEvent);
75
76using namespace std;
77
78const TString MHCamEvent::gsDefName = "MHCamEvent";
79const TString MHCamEvent::gsDefTitle = "Average of MCamEvents";
80
81// --------------------------------------------------------------------------
82//
83// Initialize the name and title of the task
84//
85void MHCamEvent::Init(const char *name, const char *title)
86{
87 //
88 // set the name and title of this object
89 //
90 fName = name ? name : gsDefName.Data();
91 fTitle = title ? title : gsDefTitle.Data();
92
93 fNameGeom = "MGeomCam";
94}
95
96// --------------------------------------------------------------------------
97//
98// Initialize the name and title of the task. Set fType to 0
99//
100MHCamEvent::MHCamEvent(const char *name, const char *title)
101: fSum(NULL), fErr(NULL), fEvt(NULL), fSub(0), fType(0), fErrorSpread(kTRUE), fErrorRelative(kFALSE),
102fThreshold(0), fUseThreshold(0), fIsRatio(false)
103{
104 Init(name, title);
105
106 fSum = new MHCamera;
107}
108
109// --------------------------------------------------------------------------
110//
111// Initialize the name and title of the task. Set fType to type
112//
113MHCamEvent::MHCamEvent(Int_t type, const char *name, const char *title)
114: fSum(NULL), fErr(NULL), fEvt(NULL), fSub(0), fType(type), fErrorSpread(kTRUE), fErrorRelative(kFALSE),
115fThreshold(0), fUseThreshold(0), fIsRatio(false)
116{
117 Init(name, title);
118
119 fSum = new MHCamera;
120}
121
122// --------------------------------------------------------------------------
123//
124// Delete the corresponding camera display if available
125//
126MHCamEvent::~MHCamEvent()
127{
128 if (fSum)
129 {
130 delete fSum;
131 fSum = 0;
132 }
133}
134
135Bool_t MHCamEvent::InitGeom(const MParList &plist)
136{
137 MGeomCam *cam = (MGeomCam*)plist.FindObject(fNameGeom, "MGeomCam");
138 if (!cam)
139 {
140 *fLog << err << fNameGeom << " [MGeomCam] not found... abort." << endl;
141 return kFALSE;
142 }
143
144 // combine name
145 const TString name = fNameEvt.IsNull() ? fName : fNameEvt;
146
147 fSum->SetGeometry(*cam, name+";avg");
148 if (fTitle!=gsDefTitle)
149 fSum->SetTitle(fTitle);
150 if (!fTitle.Contains(";"))
151 fSum->SetYTitle("a.u.");
152
153 if (fErr)
154 fErr->SetGeometry(*fSum->GetGeometry(), fErr->GetName(), fErr->GetTitle());
155
156 return kTRUE;
157}
158
159// --------------------------------------------------------------------------
160//
161// Get the event (MCamEvent) the histogram might be filled with. If
162// it is not given, it is assumed, that it is filled with the argument
163// of the Fill function.
164// Looks for the camera geometry MGeomCam and resets the sum histogram.
165//
166Bool_t MHCamEvent::SetupFill(const MParList *plist)
167{
168 fEvt = dynamic_cast<MCamEvent*>(plist->FindObject(fNameEvt, "MCamEvent"));
169 if (!fEvt)
170 {
171 if (!fNameEvt.IsNull())
172 {
173 *fLog << err << GetDescriptor() << ": No " << fNameEvt <<" [MCamEvent] available..." << endl;
174 return kFALSE;
175 }
176 *fLog << inf << GetDescriptor() << ": Assuming to get MCamEvent from Fill()..." << endl;
177 }
178
179 fSub = 0;
180 if (!fNameSub.IsNull())
181 {
182 fSub = dynamic_cast<MCamEvent*>(plist->FindObject(fNameSub, "MCamEvent"));
183 if (!fSub)
184 {
185 *fLog << err << GetDescriptor() << ": No " << fNameSub <<" [MCamEvent] available..." << endl;
186 return kFALSE;
187 }
188 *fLog << inf << GetDescriptor() << ": Will subtract '" << fNameSub << "'" << endl;
189 }
190
191 fSum->Reset();
192
193 fNumReInit = 0;
194
195 if (!InitGeom(*plist))
196 return kFALSE;
197
198 if (fUseThreshold!=kCollectMin && fUseThreshold!=kCollectMax)
199 fSum->SetBit(MHCamera::kProfile);
200 if (!fErrorSpread)
201 fSum->SetBit(MHCamera::kErrorMean);
202
203 return kTRUE;
204}
205
206// --------------------------------------------------------------------------
207//
208// The geometry read from the RunHeaders might have changed. This does not
209// effect anything in PreProcess. So we set a new geometry. We don't move
210// this away from PreProcess to support also loops without calling ReInit.
211//
212Bool_t MHCamEvent::ReInit(MParList *plist)
213{
214 return fNumReInit++==0 ? InitGeom(*plist) : kTRUE;
215}
216
217
218// --------------------------------------------------------------------------
219//
220// Fill the histograms with data from a MCamEvent-Container.
221//
222Int_t MHCamEvent::Fill(const MParContainer *par, const Stat_t w)
223{
224 const MCamEvent *evt = par ? dynamic_cast<const MCamEvent*>(par) : fEvt;
225 if (!evt)
226 {
227 *fLog << err << dbginf << "Got no MCamEvent as argument of Fill()..." << endl;
228 return kERROR;
229 }
230
231 switch (fUseThreshold)
232 {
233 case kNoBound:
234 if (fSub)
235 {
236 if (fIsRatio)
237 fSum->AddCamRatio(*evt, *fSub, fType, w);
238 else
239 fSum->AddCamDifference(*evt, *fSub, fType, w);
240 }
241 else
242 fSum->AddCamContent(*evt, fType, w);
243 break;
244
245 case kIsLowerBound:
246 case kIsUpperBound:
247 fSum->CntCamContent(*evt, fThreshold, fType, fUseThreshold>0);
248 break;
249
250 case kCollectMin:
251 fSum->SetMinCamContent(*evt, /*fThreshold,*/ fType);
252 break;
253
254 case kCollectMax:
255 fSum->SetMaxCamContent(*evt, /*fThreshold,*/ fType);
256 break;
257
258 case kMeanShift:
259 fSum->AddMeanShift(*evt, fType, w);
260 break;
261
262 case kMedianShift:
263 fSum->AddMedianShift(*evt, fType, w);
264 break;
265
266 default:
267 *fLog << err << "ERROR - MHCamEvent::Fill: Unknown type." << endl;
268 return kERROR;
269 }
270
271 return kTRUE;
272}
273
274// --------------------------------------------------------------------------
275//
276// Take the mean of the sum histogram and print all pixel indices
277// which are above sum+s*rms
278//
279void MHCamEvent::PrintOutliers(Float_t s) const
280{
281 const Double_t mean = fSum->GetMean();
282 const Double_t rms = fSum->GetRMS();
283
284 *fLog << all << underline << GetDescriptor() << ": Mean=" << mean << ", Rms=" << rms << endl;
285
286 for (UInt_t i=0; i<fSum->GetNumPixels(); i++)
287 {
288 if (!fSum->IsUsed(i))
289 continue;
290
291 if ((*fSum)[i+1]>mean+s*rms)
292 *fLog << "Contents of Pixel-Index #" << i << ": " << (*fSum)[i+1] << " > " << s << "*rms" << endl;
293 }
294}
295
296// --------------------------------------------------------------------------
297//
298// Return fSum for "sum" and fRms for "rms"
299//
300TH1 *MHCamEvent::GetHistByName(const TString name) const
301{
302 return fSum;
303}
304
305// --------------------------------------------------------------------------
306//
307// Set the camera histogram to a clone of cam
308//
309void MHCamEvent::SetHist(const MHCamera &cam)
310{
311 if (fSum)
312 delete fSum;
313
314 fSum = static_cast<MHCamera*>(cam.Clone());
315}
316
317TString MHCamEvent::Format(const char *ext) const
318{
319 TString n = fSum->GetName();
320
321 const Ssiz_t pos = n.Last(';');
322 if (pos<0)
323 return "";
324
325 n = n(0, pos);
326 n += ";";
327 n += ext;
328 return n;
329}
330
331void MHCamEvent::Paint(Option_t *opt)
332{
333 TVirtualPad *pad = gPad;
334
335 pad->cd(2);
336 if (gPad->FindObject("proj"))
337 {
338 TH1 *h=fSum->Projection("proj");
339 if (h->GetMaximum()>0)
340 gPad->SetLogy();
341
342 if (TString(opt).Contains("gaus", TString::kIgnoreCase))
343 {
344 // Create a copy of the function if it is not yet there
345 if (!h->GetFunction("gaus"))
346 h->Fit("gaus", "q");
347
348 TF1 *f = h->GetFunction("gaus");
349 if (f)
350 {
351 f->SetParameters(h->GetMean(), h->GetRMS(), h->GetMaximum());
352 f->SetLineWidth(1);
353 h->Fit("gaus", "q");
354 }
355
356 TPaveStats *s1 = dynamic_cast<TPaveStats*>(h->GetListOfFunctions()->FindObject("stats"));
357 if (s1)
358 {
359 if (s1->GetOptStat()==1111)
360 s1->SetOptStat(1000000001);
361 s1->SetOptFit(1);
362 }
363 }
364 }
365
366 pad->cd(5);
367 if (gPad->FindObject("rad"))
368 fSum->RadialProfile("rad");
369
370 pad->cd(6);
371 if (gPad->FindObject("az"))
372 fSum->AzimuthProfile("az");
373
374 pad->cd(4);
375 gPad->cd(1);
376 MHCamera *cam = (MHCamera*)gPad->FindObject("err");
377 if (cam)
378 cam->SetCamContent(*fSum, fErrorRelative ? 2 : 1);
379}
380
381void MHCamEvent::Draw(Option_t *opt)
382{
383 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
384 const Int_t col = pad->GetFillColor();
385 pad->SetBorderMode(0);
386
387 AppendPad(opt);
388
389 TString name = MString::Format("%s_1", pad->GetName());
390 TPad *p = new TPad(name,name,0.005, 0.5, 0.66, 0.995,col,0,0);
391 p->SetFrameBorderMode(0);
392 p->SetNumber(1);
393 p->Draw();
394 p->cd();
395 fSum->Draw("EPhist");
396
397 pad->cd();
398 name = MString::Format("%s_2", pad->GetName());
399 p = new TPad(name,name,0.66, 0.5, 0.995, 0.995,col,0,0);
400 p->SetFrameBorderMode(0);
401 p->SetNumber(2);
402 p->Draw();
403 p->cd();
404 TH1 *h = fSum->Projection("proj", 50);
405 h->SetTitle("Projection");
406 h->SetBit(kCanDelete);
407 h->Draw();
408
409 pad->cd();
410 name = MString::Format("%s_3", pad->GetName());
411 p = new TPad(name,name,0.005, 0.005, 3./8, 0.5,col,0,0);
412 p->SetFrameBorderMode(0);
413 p->SetNumber(3);
414 p->Draw();
415 p->cd();
416 fSum->Draw();
417
418 pad->cd();
419 name = MString::Format("%s_4", pad->GetName());
420 p = new TPad(name,name,3./8, 0.005, 6./8-0.005, 0.5,col,0,0);
421 p->SetFrameBorderMode(0);
422 p->SetNumber(4);
423 p->Draw();
424 p->cd();
425
426 TString e = "Sqrt(Variance)";
427 if (fSum->TestBit(MHCamera::kErrorMean))
428 e += "/Sqrt(n_{i})";
429 if (fErrorRelative)
430 e += "/v_{i}";
431
432 fErr = new MHCamera(*fSum->GetGeometry());
433 fErr->SetName("err");
434 fErr->SetTitle(e);
435 fErr->SetYTitle(fSum->GetYaxis()->GetTitle());
436 fErr->SetCamContent(*fSum, fErrorRelative ? 2 : 1);
437 fErr->SetBit(kMustCleanup);
438 fErr->SetBit(kCanDelete);
439 fErr->Draw();
440
441 pad->cd();
442 name = MString::Format("%s_5", pad->GetName());
443 p = new TPad(name,name,6./8,0.25,0.995,0.5,col,0,0);
444 p->SetFrameBorderMode(0);
445 p->SetNumber(5);
446 p->Draw();
447 p->cd();
448 h = (TH1*)fSum->RadialProfile("rad", 20);
449 h->SetTitle("Radial Profile");
450 h->SetBit(kCanDelete|TH1::kNoStats);
451 h->Draw();
452
453 pad->cd();
454 name = MString::Format("%s_6", pad->GetName());
455 p = new TPad(name,name,6./8,0.005,0.995,0.25,col,0,0);
456 p->SetFrameBorderMode(0);
457 p->SetNumber(6);
458 p->Draw();
459 p->cd();
460 h = (TH1*)fSum->AzimuthProfile("az", 30);
461 h->SetTitle("Azimuth Profile");
462 h->SetBit(kCanDelete|TH1::kNoStats);
463 h->Draw();
464}
465
466
467void MHCamEvent::RecursiveRemove(TObject *obj)
468{
469 if (obj==fErr)
470 fErr = 0;
471 if (obj==fSum)
472 fSum = 0;
473
474 MH::RecursiveRemove(obj);
475}
Note: See TracBrowser for help on using the repository browser.