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

Last change on this file since 13398 was 13398, checked in by tbretz, 10 years ago
Update the statistics box with fit information if a gauss fit is done.
File size: 12.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): 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 delete fSum;
130}
131
132Bool_t MHCamEvent::InitGeom(const MParList &plist)
133{
134 MGeomCam *cam = (MGeomCam*)plist.FindObject(fNameGeom, "MGeomCam");
135 if (!cam)
136 {
137 *fLog << err << fNameGeom << " [MGeomCam] not found... abort." << endl;
138 return kFALSE;
139 }
140
141 // combine name
142 const TString name = fNameEvt.IsNull() ? fName : fNameEvt;
143
144 fSum->SetGeometry(*cam, name+";avg");
145 if (fTitle!=gsDefTitle)
146 fSum->SetTitle(fTitle);
147 if (!fTitle.Contains(";"))
148 fSum->SetYTitle("a.u.");
149
150 if (fErr)
151 fErr->SetGeometry(*fSum->GetGeometry(), fErr->GetName(), fErr->GetTitle());
152
153 return kTRUE;
154}
155
156// --------------------------------------------------------------------------
157//
158// Get the event (MCamEvent) the histogram might be filled with. If
159// it is not given, it is assumed, that it is filled with the argument
160// of the Fill function.
161// Looks for the camera geometry MGeomCam and resets the sum histogram.
162//
163Bool_t MHCamEvent::SetupFill(const MParList *plist)
164{
165 fEvt = dynamic_cast<MCamEvent*>(plist->FindObject(fNameEvt, "MCamEvent"));
166 if (!fEvt)
167 {
168 if (!fNameEvt.IsNull())
169 {
170 *fLog << err << GetDescriptor() << ": No " << fNameEvt <<" [MCamEvent] available..." << endl;
171 return kFALSE;
172 }
173 *fLog << inf << GetDescriptor() << ": Assuming to get MCamEvent from Fill()..." << endl;
174 }
175
176 fSub = 0;
177 if (!fNameSub.IsNull())
178 {
179 fSub = dynamic_cast<MCamEvent*>(plist->FindObject(fNameSub, "MCamEvent"));
180 if (!fSub)
181 {
182 *fLog << err << GetDescriptor() << ": No " << fNameSub <<" [MCamEvent] available..." << endl;
183 return kFALSE;
184 }
185 *fLog << inf << GetDescriptor() << ": Will subtract '" << fNameSub << "'" << endl;
186 }
187
188 fSum->Reset();
189
190 fNumReInit = 0;
191
192 if (!InitGeom(*plist))
193 return kFALSE;
194
195 if (fUseThreshold!=kCollectMin && fUseThreshold!=kCollectMax)
196 fSum->SetBit(MHCamera::kProfile);
197 if (!fErrorSpread)
198 fSum->SetBit(MHCamera::kErrorMean);
199
200 return kTRUE;
201}
202
203// --------------------------------------------------------------------------
204//
205// The geometry read from the RunHeaders might have changed. This does not
206// effect anything in PreProcess. So we set a new geometry. We don't move
207// this away from PreProcess to support also loops without calling ReInit.
208//
209Bool_t MHCamEvent::ReInit(MParList *plist)
210{
211 return fNumReInit++==0 ? InitGeom(*plist) : kTRUE;
212}
213
214
215// --------------------------------------------------------------------------
216//
217// Fill the histograms with data from a MCamEvent-Container.
218//
219Int_t MHCamEvent::Fill(const MParContainer *par, const Stat_t w)
220{
221 const MCamEvent *evt = par ? dynamic_cast<const MCamEvent*>(par) : fEvt;
222 if (!evt)
223 {
224 *fLog << err << dbginf << "Got no MCamEvent as argument of Fill()..." << endl;
225 return kERROR;
226 }
227
228 switch (fUseThreshold)
229 {
230 case kNoBound:
231 if (fSub)
232 {
233 if (fIsRatio)
234 fSum->AddCamRatio(*evt, *fSub, fType, w);
235 else
236 fSum->AddCamDifference(*evt, *fSub, fType, w);
237 }
238 else
239 fSum->AddCamContent(*evt, fType, w);
240 break;
241
242 case kIsLowerBound:
243 case kIsUpperBound:
244 fSum->CntCamContent(*evt, fThreshold, fType, fUseThreshold>0);
245 break;
246
247 case kCollectMin:
248 fSum->SetMinCamContent(*evt, /*fThreshold,*/ fType);
249 break;
250
251 case kCollectMax:
252 fSum->SetMaxCamContent(*evt, /*fThreshold,*/ fType);
253 break;
254
255 case kMeanShift:
256 fSum->AddMeanShift(*evt, fType, w);
257 break;
258
259 case kMedianShift:
260 fSum->AddMedianShift(*evt, fType, w);
261 break;
262
263 default:
264 *fLog << err << "ERROR - MHCamEvent::Fill: Unknown type." << endl;
265 return kERROR;
266 }
267
268 return kTRUE;
269}
270
271// --------------------------------------------------------------------------
272//
273// Take the mean of the sum histogram and print all pixel indices
274// which are above sum+s*rms
275//
276void MHCamEvent::PrintOutliers(Float_t s) const
277{
278 const Double_t mean = fSum->GetMean();
279 const Double_t rms = fSum->GetRMS();
280
281 *fLog << all << underline << GetDescriptor() << ": Mean=" << mean << ", Rms=" << rms << endl;
282
283 for (UInt_t i=0; i<fSum->GetNumPixels(); i++)
284 {
285 if (!fSum->IsUsed(i))
286 continue;
287
288 if ((*fSum)[i+1]>mean+s*rms)
289 *fLog << "Contents of Pixel-Index #" << i << ": " << (*fSum)[i+1] << " > " << s << "*rms" << endl;
290 }
291}
292
293// --------------------------------------------------------------------------
294//
295// Return fSum for "sum" and fRms for "rms"
296//
297TH1 *MHCamEvent::GetHistByName(const TString name) const
298{
299 return fSum;
300}
301
302// --------------------------------------------------------------------------
303//
304// Set the camera histogram to a clone of cam
305//
306void MHCamEvent::SetHist(const MHCamera &cam)
307{
308 if (fSum)
309 delete fSum;
310
311 fSum = static_cast<MHCamera*>(cam.Clone());
312}
313
314TString MHCamEvent::Format(const char *ext) const
315{
316 TString n = fSum->GetName();
317
318 const Ssiz_t pos = n.Last(';');
319 if (pos<0)
320 return "";
321
322 n = n(0, pos);
323 n += ";";
324 n += ext;
325 return n;
326}
327
328void MHCamEvent::Paint(Option_t *opt)
329{
330 TVirtualPad *pad = gPad;
331
332 pad->cd(2);
333 if (gPad->FindObject("proj"))
334 {
335 TH1 *h=fSum->Projection("proj");
336 if (h->GetMaximum()>0)
337 gPad->SetLogy();
338
339 if (TString(opt).Contains("gaus", TString::kIgnoreCase))
340 {
341 // Create a copy of the function if it is not yet there
342 if (!h->GetFunction("gaus"))
343 h->Fit("gaus", "q");
344
345 TF1 *f = h->GetFunction("gaus");
346 if (f)
347 {
348 f->SetParameters(h->GetMean(), h->GetRMS(), h->GetMaximum());
349 f->SetLineWidth(1);
350 h->Fit("gaus", "q");
351 }
352
353 TPaveStats *s1 = dynamic_cast<TPaveStats*>(h->GetListOfFunctions()->FindObject("stats"));
354 if (s1)
355 {
356 if (s1->GetOptStat()==1111)
357 s1->SetOptStat(1000000001);
358 s1->SetOptFit(1);
359 }
360 }
361 }
362
363 pad->cd(5);
364 if (gPad->FindObject("rad"))
365 fSum->RadialProfile("rad");
366
367 pad->cd(6);
368 if (gPad->FindObject("az"))
369 fSum->AzimuthProfile("az");
370
371 pad->cd(4);
372 gPad->cd(1);
373 MHCamera *cam = (MHCamera*)gPad->FindObject("err");
374 if (cam)
375 cam->SetCamContent(*fSum, fErrorRelative ? 2 : 1);
376}
377
378void MHCamEvent::Draw(Option_t *opt)
379{
380 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
381 const Int_t col = pad->GetFillColor();
382 pad->SetBorderMode(0);
383
384 AppendPad(opt);
385
386 TString name = MString::Format("%s_1", pad->GetName());
387 TPad *p = new TPad(name,name,0.005, 0.5, 0.66, 0.995,col,0,0);
388 p->SetFrameBorderMode(0);
389 p->SetNumber(1);
390 p->Draw();
391 p->cd();
392 fSum->Draw("EPhist");
393
394 pad->cd();
395 name = MString::Format("%s_2", pad->GetName());
396 p = new TPad(name,name,0.66, 0.5, 0.995, 0.995,col,0,0);
397 p->SetFrameBorderMode(0);
398 p->SetNumber(2);
399 p->Draw();
400 p->cd();
401 TH1 *h = fSum->Projection("proj", 50);
402 h->SetTitle("Projection");
403 h->SetBit(kCanDelete);
404 h->Draw();
405
406 pad->cd();
407 name = MString::Format("%s_3", pad->GetName());
408 p = new TPad(name,name,0.005, 0.005, 3./8, 0.5,col,0,0);
409 p->SetFrameBorderMode(0);
410 p->SetNumber(3);
411 p->Draw();
412 p->cd();
413 fSum->Draw();
414
415 pad->cd();
416 name = MString::Format("%s_4", pad->GetName());
417 p = new TPad(name,name,3./8, 0.005, 6./8-0.005, 0.5,col,0,0);
418 p->SetFrameBorderMode(0);
419 p->SetNumber(4);
420 p->Draw();
421 p->cd();
422
423 TString e = "Sqrt(Variance)";
424 if (fSum->TestBit(MHCamera::kErrorMean))
425 e += "/Sqrt(n_{i})";
426 if (fErrorRelative)
427 e += "/v_{i}";
428
429 fErr = new MHCamera(*fSum->GetGeometry());
430 fErr->SetName("err");
431 fErr->SetTitle(e);
432 fErr->SetYTitle(fSum->GetYaxis()->GetTitle());
433 fErr->SetCamContent(*fSum, fErrorRelative ? 2 : 1);
434 fErr->SetBit(kMustCleanup);
435 fErr->SetBit(kCanDelete);
436 fErr->Draw();
437
438 pad->cd();
439 name = MString::Format("%s_5", pad->GetName());
440 p = new TPad(name,name,6./8,0.25,0.995,0.5,col,0,0);
441 p->SetFrameBorderMode(0);
442 p->SetNumber(5);
443 p->Draw();
444 p->cd();
445 h = (TH1*)fSum->RadialProfile("rad", 20);
446 h->SetTitle("Radial Profile");
447 h->SetBit(kCanDelete|TH1::kNoStats);
448 h->Draw();
449
450 pad->cd();
451 name = MString::Format("%s_6", pad->GetName());
452 p = new TPad(name,name,6./8,0.005,0.995,0.25,col,0,0);
453 p->SetFrameBorderMode(0);
454 p->SetNumber(6);
455 p->Draw();
456 p->cd();
457 h = (TH1*)fSum->AzimuthProfile("az", 30);
458 h->SetTitle("Azimuth Profile");
459 h->SetBit(kCanDelete|TH1::kNoStats);
460 h->Draw();
461}
462
463
464void MHCamEvent::RecursiveRemove(TObject *obj)
465{
466 if (obj==fErr)
467 fErr = 0;
468 if (obj==fSum)
469 fSum = 0;
470
471 MH::RecursiveRemove(obj);
472}
Note: See TracBrowser for help on using the repository browser.