source: trunk/MagicSoft/Mars/mcalib/MCalibrationRelTimeCalc.cc@ 8462

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