source: trunk/Mars/mcalib/MCalibrationRelTimeCalc.cc@ 11491

Last change on this file since 11491 was 9374, checked in by tbretz, 17 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): Markus Gaug 04/2004 <mailto:markus@ifae.es>
19! Author(s): Thomas Bretz <mailto:tbretz@astro.uni-wuerzburg.de>
20!
21! Copyright: MAGIC Software Development, 2000-2007
22!
23!
24\* ======================================================================== */
25
26//////////////////////////////////////////////////////////////////////////////
27//
28// MCalibrationRelTimeCalc
29//
30// Task to finalize the relative time calibration obtained
31// from the fit results to the summed FADC slice distributions delivered by
32// MCalibrationRelTimeCam, calculated and filled by MHCalibrationRelTimeCam,
33//
34// PreProcess(): Initialize pointers to MCalibrationRelTimeCam,
35//
36// ReInit(): Initializes pointer to MBadPixelsCam
37//
38// Process(): Nothing to be done, histograms getting filled by MHCalibrationChargeCam
39//
40// PostProcess(): - FinalizeRelTimes()
41// - FinalizeBadPixels()
42//
43// Class Version 2:
44// + Byte_t fCheckFlags; // Bit-field to hold the possible check flags
45//
46// ClassVersionb 3:
47// - TString fOutputPath; // Path to the output file
48// - TString fOutputFile; // Name of the output file
49//
50//
51// Input Containers:
52// MCalibrationRelTimeCam
53// MBadPixelsCam
54// MGeomCam
55//
56// Output Containers:
57// MCalibrationRelTimeCam
58// MBadPixelsCam
59//
60//
61//////////////////////////////////////////////////////////////////////////////
62#include "MCalibrationRelTimeCalc.h"
63
64#include "MLog.h"
65#include "MLogManip.h"
66
67#include "MMath.h"
68
69#include "MParList.h"
70
71#include "MStatusDisplay.h"
72
73#include "MGeomCam.h"
74#include "MGeom.h"
75
76#include "MCalibrationRelTimeCam.h"
77#include "MCalibrationRelTimePix.h"
78
79#include "MBadPixelsCam.h"
80#include "MBadPixelsPix.h"
81
82ClassImp(MCalibrationRelTimeCalc);
83
84using namespace std;
85
86const Float_t MCalibrationRelTimeCalc::fgRelTimeResolutionLimit = 5.0;
87
88// --------------------------------------------------------------------------
89//
90// Default constructor.
91//
92// Sets all pointers to NULL
93//
94// Initializes:
95// - fRelTimeResolutionLimit to fgRelTimeResolutionimit
96//
97// Calls:
98// - Clear()
99//
100MCalibrationRelTimeCalc::MCalibrationRelTimeCalc(const char *name, const char *title)
101 : fGeom(NULL), fFlags(0)
102{
103
104 fName = name ? name : "MCalibrationRelTimeCalc";
105 fTitle = title ? title : "Task to finalize the relative time calibration";
106
107 SetCheckFitResults ( kFALSE );
108 SetCheckDeviatingBehavior( kTRUE );
109 SetCheckHistOverflow ( kFALSE );
110 SetCheckOscillations ( kFALSE );
111
112 SetRelTimeResolutionLimit();
113
114 Clear();
115
116}
117
118// --------------------------------------------------------------------------
119//
120// Sets:
121// - all flags to kFALSE
122// - all pointers to NULL
123//
124void MCalibrationRelTimeCalc::Clear(const Option_t *o)
125{
126 fBadPixels = NULL;
127 fCam = NULL;
128}
129
130// --------------------------------------------------------------------------
131//
132// Search for the following input containers and abort if not existing:
133// - MGeomCam
134// - MCalibrationRelTimeCam
135// - MBadPixelsCam
136//
137// It defines the PixId of every pixel in:
138//
139// - MCalibrationRelTimeCam
140//
141// It sets all pixels in excluded which have the flag fBadBixelsPix::IsBad() set in:
142//
143// - MCalibrationRelTimePix
144//
145Bool_t MCalibrationRelTimeCalc::ReInit(MParList *pList )
146{
147
148 fGeom = (MGeomCam*)pList->FindObject("MGeomCam");
149 if (!fGeom)
150 {
151 *fLog << err << "No MGeomCam found... aborting." << endl;
152 return kFALSE;
153 }
154
155 fBadPixels = (MBadPixelsCam*)pList->FindObject(AddSerialNumber("MBadPixelsCam"));
156 if (!fBadPixels)
157 {
158 *fLog << err << "Cannot find MBadPixelsCam ... abort." << endl;
159 return kFALSE;
160 }
161
162 fCam = (MCalibrationRelTimeCam*)pList->FindObject(AddSerialNumber("MCalibrationRelTimeCam"));
163 if (!fCam)
164 {
165 *fLog << err << "Cannot find MCalibrationRelTimeCam ... abort." << endl;
166 return kFALSE;
167 }
168
169 if (IsDebug())
170 {
171 const UInt_t npixels = fGeom->GetNumPixels();
172 for (UInt_t i=0; i<npixels; i++)
173 (*fCam)[i].SetDebug();
174 }
175
176 return kTRUE;
177}
178
179// -----------------------------------------------------------------------
180//
181// Return if number of executions is null.
182//
183Int_t MCalibrationRelTimeCalc::PostProcess()
184{
185
186 if (GetNumExecutions()==0)
187 return kTRUE;
188
189 return Finalize();
190}
191
192// -----------------------------------------------------------------------
193//
194// First loop over pixels, average areas and sectors, call:
195// - FinalizeRelTimes()
196// for every entry. Count number of valid pixels in loop and return kFALSE
197// if there are none (the "Michele check").
198//
199// Call FinalizeBadPixels()
200//
201// Call MParContainer::SetReadyToSave() for fCam
202//
203// Print out some statistics
204//
205Int_t MCalibrationRelTimeCalc::Finalize()
206{
207
208 //
209 // First loop over pixels, call FinalizePedestals and FinalizeRelTimes
210 //
211 FinalizeRelTimes();
212
213 //
214 // Finalize Bad Pixels
215 //
216 FinalizeBadPixels();
217
218 //
219 // Finalize calibration statistics
220 //
221 FinalizeUnsuitablePixels();
222
223 fCam->SetReadyToSave();
224 fBadPixels->SetReadyToSave();
225
226 *fLog << inf << endl;
227 *fLog << GetDescriptor() << ": Errors statistics:" << endl;
228
229 PrintUncalibrated(MBadPixelsPix::kDeviatingRelTimeResolution,
230 Form("%s%2.1f%s","Rel.time rms more than ", fRelTimeResolutionLimit, " dev from median:"));
231 PrintUncalibrated(MBadPixelsPix::kRelTimeOscillating,
232 "Pixels with changing Rel. Times over time:");
233 PrintUncalibrated(MBadPixelsPix::kRelTimeNotFitted,
234 "Pixels with unsuccesful Gauss fit to the times:");
235
236 return kTRUE;
237}
238
239
240// ----------------------------------------------------------------------------------------------------
241//
242// Check for outliers. They are marked with
243// MBadPixelsPix::kDeviatingTimeResolution
244//
245// see also MCalibrationChargeCalc::FinalizeAbsTimes
246//
247void MCalibrationRelTimeCalc::FinalizeRelTimes()
248{
249 const Int_t npixels = fGeom->GetNumPixels();
250 const Int_t nareas = fGeom->GetNumAreas();
251
252 // Create an array capable of holding all pixels
253 TArrayF arr(npixels);
254
255 for (Int_t aidx=0; aidx<nareas; aidx++)
256 {
257 Int_t n = 0;
258 for (Int_t i=0; i<npixels; i++)
259 {
260 // Check for this aidx only
261 if ((*fGeom)[i].GetAidx()!=aidx)
262 continue;
263
264 // check if pixel may not contain a valid value
265 if ((*fBadPixels)[i].IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
266 continue;
267
268 // check if it was excluded for some reason
269 const MCalibrationRelTimePix &pix = (MCalibrationRelTimePix&)(*fCam)[i];
270 if (pix.IsExcluded())
271 continue;
272
273 // if TimePrecision is valid fill it into array
274 if (pix.GetTimePrecision()>0)
275 arr[n++] = pix.GetTimePrecision();
276 }
277
278 // Check the ratio of valid entries to the ratio of pixels
279 const Float_t ratio = 100*n/fGeom->GetNumPixWithAidx(aidx);
280 if (3*ratio<2)
281 *fLog << warn << "Area " << setw(4) << aidx << ": Only " << ratio << "% pixels with valid time resolution found." << endl;
282
283 // Calculate median and median deviation
284 Double_t med;
285 const Double_t dev = MMath::MedianDev(n, arr.GetArray(), med);
286
287 // Calculate upper and lower limit
288 const Float_t lolim = TMath::Max(med-fRelTimeResolutionLimit*dev, 0.);
289 const Float_t hilim = TMath::Max(med+fRelTimeResolutionLimit*dev, 0.);
290
291 // Now find the outliers
292 for (Int_t i=0; i<npixels; i++)
293 {
294 // Search only within this aidx
295 if ((*fGeom)[i].GetAidx()!=aidx)
296 continue;
297
298 // skip pixels already known to be unsuitable
299 if ((*fBadPixels)[i].IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
300 continue;
301
302 // check if a pixel has been excluded. This
303 const MCalibrationRelTimePix &pix = (MCalibrationRelTimePix&)(*fCam)[i];
304
305 // Check if time precision is valid (might be invalid
306 // for example in cae of empty histograms)
307 const Float_t res = pix.GetTimePrecision();
308 if (res<0) //FIXME!!! How does this happen?
309 {
310 *fLog << warn << "Pixel " << setw(4) << i << ": Rel-time rms could not be calculated." << endl;
311 (*fBadPixels)[i].SetUncalibrated(MBadPixelsPix::kDeviatingRelTimeResolution);
312 continue;
313 }
314
315 // Now compare to a lower and upper limit
316 if (res<=lolim || res>=hilim)
317 {
318 *fLog << warn << "Pixel " << setw(4) << i << ": Deviation from rel-time rms: "
319 << Form("%5.2f", res) << " out of range "
320 << Form("[%4.2f,%4.2f]", lolim, hilim) << endl;
321
322 (*fBadPixels)[i].SetUncalibrated(MBadPixelsPix::kDeviatingRelTimeResolution);
323 }
324 }
325 }
326}
327
328
329// -----------------------------------------------------------------------------------
330//
331// Sets pixel to MBadPixelsPix::kUnsuitableRun, if one of the following flags is set:
332// - MBadPixelsPix::kRelTimeIsPedestal
333// - MBadPixelsPix::kRelTimeErrNotValid
334// - MBadPixelsPix::kRelTimeRelErrNotValid
335//
336void MCalibrationRelTimeCalc::FinalizeBadPixels()
337{
338
339 for (Int_t i=0; i<fBadPixels->GetSize(); i++)
340 {
341 MBadPixelsPix &bad = (*fBadPixels)[i];
342
343 if (IsCheckDeviatingBehavior())
344 if (bad.IsUncalibrated(MBadPixelsPix::kDeviatingRelTimeResolution))
345 bad.SetUnsuitable(MBadPixelsPix::kUnsuitableRun);
346
347 if (IsCheckFitResults())
348 if (bad.IsUncalibrated(MBadPixelsPix::kRelTimeNotFitted))
349 bad.SetUnsuitable(MBadPixelsPix::kUnreliableRun);
350
351 if (IsCheckOscillations())
352 if (bad.IsUncalibrated(MBadPixelsPix::kRelTimeOscillating))
353 bad.SetUnsuitable(MBadPixelsPix::kUnreliableRun);
354 }
355
356}
357
358
359// -----------------------------------------------------------------------------------------------
360//
361// - Print out statistics about BadPixels of type UnsuitableType_t
362// - store numbers of bad pixels of each type in fIntensCam or fCam
363//
364void MCalibrationRelTimeCalc::FinalizeUnsuitablePixels()
365{
366 *fLog << inf << endl;
367 *fLog << GetDescriptor() << ": Rel. Times Calibration status:" << endl;
368 *fLog << dec;
369
370 const Int_t nareas = fGeom->GetNumAreas();
371
372 TArrayI unsuit(nareas);
373 TArrayI unrel(nareas);
374
375 for (int aidx=0; aidx<nareas; aidx++)
376 {
377 unsuit[aidx] += fBadPixels->GetNumUnsuitable(MBadPixelsPix::kUnsuitableRun, fGeom, aidx);
378 unrel[aidx] += fBadPixels->GetNumUnsuitable(MBadPixelsPix::kUnreliableRun, fGeom, aidx);
379 fCam->SetNumUnsuitable(unsuit[aidx], aidx);
380 fCam->SetNumUnreliable(unrel[aidx], aidx);
381 }
382
383 if (fGeom->InheritsFrom("MGeomCamMagic"))
384 {
385 PrintUncalibrated("Uncalibrated Pixels:", unsuit[0], unsuit[1]);
386 PrintUncalibrated("Unreliable Pixels:", unrel[0], unrel[1]);
387 }
388}
389
390// -----------------------------------------------------------------------------------------------
391//
392// Print out statistics about BadPixels of type UncalibratedType_t
393//
394void MCalibrationRelTimeCalc::PrintUncalibrated(MBadPixelsPix::UncalibratedType_t typ, const char *text) const
395{
396 UInt_t countinner = 0;
397 UInt_t countouter = 0;
398 for (Int_t i=0; i<fBadPixels->GetSize(); i++)
399 {
400 if ((*fBadPixels)[i].IsUncalibrated(typ))
401 {
402 if (fGeom->GetPixRatio(i) == 1.)
403 countinner++;
404 else
405 countouter++;
406 }
407 }
408
409 PrintUncalibrated(text, countinner, countouter);
410}
411
412void MCalibrationRelTimeCalc::PrintUncalibrated(const char *text, Int_t in, Int_t out) const
413{
414 *fLog << " " << setfill(' ') << setw(48) << setiosflags(ios::left) << text;
415 *fLog << " Inner: " << Form("%3i", in);
416 *fLog << " Outer: " << Form("%3i", out);
417 *fLog << resetiosflags(ios::left) << endl;
418
419}
420
421// --------------------------------------------------------------------------
422//
423// MCalibrationRelTimeCam.CheckFitResults: Yes
424// MCalibrationRelTimeCam.CheckDeviatingBehavior: Yes
425// MCalibrationRelTimeCam.CheckHistOverflow: Yes
426// MCalibrationRelTimeCam.CheckOscillations: Yes
427//
428Int_t MCalibrationRelTimeCalc::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
429{
430 Bool_t rc = kFALSE;
431
432 if (IsEnvDefined(env, prefix, "CheckFitResults", print))
433 {
434 SetCheckFitResults(GetEnvValue(env, prefix, "CheckFitResults", IsCheckFitResults()));
435 rc = kTRUE;
436 }
437 if (IsEnvDefined(env, prefix, "CheckDeviatingBehavior", print))
438 {
439 SetCheckDeviatingBehavior(GetEnvValue(env, prefix, "CheckDeviatingBehavior", IsCheckDeviatingBehavior()));
440 rc = kTRUE;
441 }
442 if (IsEnvDefined(env, prefix, "CheckHistOverflow", print))
443 {
444 SetCheckHistOverflow(GetEnvValue(env, prefix, "CheckHistOverflow", IsCheckHistOverflow()));
445 rc = kTRUE;
446 }
447
448 if (IsEnvDefined(env, prefix, "CheckOscillations", print))
449 {
450 SetCheckOscillations(GetEnvValue(env, prefix, "CheckOscillations", IsCheckOscillations()));
451 rc = kTRUE;
452 }
453 if (IsEnvDefined(env, prefix, "RelTimeResolutionLimit", print))
454 {
455 SetRelTimeResolutionLimit(GetEnvValue(env, prefix, "RelTimeResolutionLimit", fRelTimeResolutionLimit));
456 rc = kTRUE;
457 }
458
459 return rc;
460}
Note: See TracBrowser for help on using the repository browser.