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

Last change on this file since 8167 was 8142, checked in by tbretz, 19 years ago
*** empty log message ***
File size: 16.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): 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 "MParList.h"
67
68#include "MStatusDisplay.h"
69
70#include "MGeomCam.h"
71#include "MGeomPix.h"
72
73#include "MCalibrationIntensityRelTimeCam.h"
74#include "MCalibrationRelTimeCam.h"
75#include "MCalibrationRelTimePix.h"
76
77#include "MBadPixelsIntensityCam.h"
78#include "MBadPixelsCam.h"
79#include "MBadPixelsPix.h"
80
81ClassImp(MCalibrationRelTimeCalc);
82
83using namespace std;
84
85const Float_t MCalibrationRelTimeCalc::fgRelTimeResolutionLimit = 1.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( kFALSE );
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
126 fIntensBad = NULL;
127 fBadPixels = NULL;
128 fCam = NULL;
129 fIntensCam = NULL;
130
131}
132
133// --------------------------------------------------------------------------
134//
135// Search for the following input containers and abort if not existing:
136// - MGeomCam
137// - MCalibrationIntensityRelTimeCam or MCalibrationRelTimeCam
138// - MBadPixelsIntensityCam or MBadPixelsCam
139//
140// It defines the PixId of every pixel in:
141//
142// - MCalibrationRelTimeCam
143//
144// It sets all pixels in excluded which have the flag fBadBixelsPix::IsBad() set in:
145//
146// - MCalibrationRelTimePix
147//
148Bool_t MCalibrationRelTimeCalc::ReInit(MParList *pList )
149{
150
151 fGeom = (MGeomCam*)pList->FindObject("MGeomCam");
152 if (!fGeom)
153 {
154 *fLog << err << "No MGeomCam found... aborting." << endl;
155 return kFALSE;
156 }
157
158 fIntensBad = (MBadPixelsIntensityCam*)pList->FindObject(AddSerialNumber("MBadPixelsIntensityCam"));
159 if (fIntensBad)
160 *fLog << inf << "Found MBadPixelsIntensityCam ... " << endl;
161 else
162 {
163 fBadPixels = (MBadPixelsCam*)pList->FindObject(AddSerialNumber("MBadPixelsCam"));
164 if (!fBadPixels)
165 {
166 *fLog << err << "Cannot find MBadPixelsCam ... abort." << endl;
167 return kFALSE;
168 }
169 }
170
171 fIntensCam = (MCalibrationIntensityRelTimeCam*)pList->FindObject(AddSerialNumber("MCalibrationIntensityRelTimeCam"));
172 if (fIntensCam)
173 *fLog << inf << "Found MCalibrationIntensityRelTimeCam ... " << endl;
174 else
175 {
176 fCam = (MCalibrationRelTimeCam*)pList->FindObject(AddSerialNumber("MCalibrationRelTimeCam"));
177 if (!fCam)
178 {
179 *fLog << err << "Cannot find MCalibrationRelTimeCam ... abort." << endl;
180 *fLog << err << "Maybe you forget to call an MFillH for the MHCalibrationRelTimeCam before..." << endl;
181 return kFALSE;
182 }
183 }
184
185 UInt_t npixels = fGeom->GetNumPixels();
186
187 MCalibrationRelTimeCam *relcam = fIntensCam
188 ? (MCalibrationRelTimeCam*)fIntensCam->GetCam() : fCam;
189 MBadPixelsCam *badcam = fIntensBad
190 ? (MBadPixelsCam*) fIntensBad->GetCam() : fBadPixels;
191
192 for (UInt_t i=0; i<npixels; i++)
193 {
194
195 MCalibrationRelTimePix &pix = (MCalibrationRelTimePix&)(*relcam)[i];
196 MBadPixelsPix &bad = (*badcam)[i];
197
198 if (bad.IsBad())
199 {
200 pix.SetExcluded();
201 continue;
202 }
203
204 if (IsDebug())
205 pix.SetDebug();
206 }
207
208 return kTRUE;
209}
210
211// -----------------------------------------------------------------------
212//
213// Return if number of executions is null.
214//
215Int_t MCalibrationRelTimeCalc::PostProcess()
216{
217
218 if (GetNumExecutions()==0)
219 return kTRUE;
220
221 return Finalize();
222}
223
224// -----------------------------------------------------------------------
225//
226// First loop over pixels, average areas and sectors, call:
227// - FinalizeRelTimes()
228// for every entry. Count number of valid pixels in loop and return kFALSE
229// if there are none (the "Michele check").
230//
231// Call FinalizeBadPixels()
232//
233// Call MParContainer::SetReadyToSave() for fCam
234//
235// Print out some statistics
236//
237Int_t MCalibrationRelTimeCalc::Finalize()
238{
239
240 //
241 // First loop over pixels, call FinalizePedestals and FinalizeRelTimes
242 //
243 FinalizeRelTimes();
244
245 //
246 // Finalize Bad Pixels
247 //
248 FinalizeBadPixels();
249
250 //
251 // Finalize calibration statistics
252 //
253 FinalizeUnsuitablePixels();
254
255 if (fIntensCam)
256 fIntensCam->SetReadyToSave();
257 else
258 fCam ->SetReadyToSave();
259
260 if (fIntensBad)
261 fIntensBad->SetReadyToSave();
262 else
263 fBadPixels->SetReadyToSave();
264
265 *fLog << inf << endl;
266 *fLog << GetDescriptor() << ": Errors statistics:" << endl;
267
268 PrintUncalibrated(MBadPixelsPix::kDeviatingTimeResolution,
269 Form("%s%2.1f%s","Time resol. less than ", fRelTimeResolutionLimit, " FADC sl. from Mean:"));
270 PrintUncalibrated(MBadPixelsPix::kRelTimeOscillating,
271 "Pixels with changing Rel. Times over time:");
272 PrintUncalibrated(MBadPixelsPix::kRelTimeNotFitted,
273 "Pixels with unsuccesful Gauss fit to the times:");
274
275 return kTRUE;
276}
277
278
279// ----------------------------------------------------------------------------------------------------
280//
281//
282// First loop: Calculate a mean and mean RMS of time resolution per area index
283// Include only pixels which are not MBadPixelsPix::kUnsuitableRun or
284// MBadPixelsPix::kUnreliableRun (see FinalizeBadPixels())
285//
286// Second loop: Exclude those deviating by more than fRelTimeResolutionLimit FADC slices
287// from the mean (obtained in first loop). Set
288// MBadPixelsPix::kDeviatingTimeResolution if excluded.
289//
290void MCalibrationRelTimeCalc::FinalizeRelTimes()
291{
292
293 MCalibrationRelTimeCam *relcam = fIntensCam
294 ? (MCalibrationRelTimeCam*)fIntensCam->GetCam() : fCam;
295 MBadPixelsCam *badcam = fIntensBad
296 ? (MBadPixelsCam*) fIntensBad->GetCam() : fBadPixels;
297
298 const UInt_t npixels = fGeom->GetNumPixels();
299 const UInt_t nareas = fGeom->GetNumAreas();
300
301 TArrayF lowlim (nareas);
302 TArrayF upplim (nareas);
303 TArrayF areasum (nareas);
304 // Float_t areasum2 [nareas];
305 TArrayI numareavalid (nareas);
306 TArrayI useunreliable(nareas);
307
308 //
309 // Apero loop: Count number of unreliable pixels:
310 //
311 for (UInt_t i=0; i<npixels; i++)
312 {
313 MBadPixelsPix &bad = (*badcam)[i];
314 const Int_t aidx = (*fGeom)[i].GetAidx();
315
316 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
317 continue;
318
319 if (bad.IsUnsuitable(MBadPixelsPix::kUnreliableRun))
320 continue;
321
322 numareavalid[aidx] ++;
323 }
324
325 for (UInt_t aidx=0; aidx<nareas; aidx++)
326 if (numareavalid[aidx] < 100)
327 useunreliable[aidx] = 1;
328
329 numareavalid.Reset();
330 //
331 // First loop: Get mean time resolution the RMS
332 // The loop is only to recognize later pixels with very deviating numbers
333 //
334 for (UInt_t i=0; i<npixels; i++)
335 {
336
337 MCalibrationRelTimePix &pix = (MCalibrationRelTimePix&)(*relcam)[i];
338 MBadPixelsPix &bad = (*badcam)[i];
339
340 if (pix.IsExcluded())
341 continue;
342
343 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
344 continue;
345
346 const Int_t aidx = (*fGeom)[i].GetAidx();
347
348 if (!useunreliable[aidx])
349 if (bad.IsUnsuitable(MBadPixelsPix::kUnreliableRun))
350 continue;
351
352 const Float_t res = pix.GetTimePrecision();
353
354 areasum [aidx] += res;
355 // areasum2 [aidx] += res*res;
356 numareavalid[aidx] ++;
357 }
358
359
360 for (UInt_t aidx=0; aidx<nareas; aidx++)
361 {
362 if (numareavalid[aidx] < 20)
363 {
364 *fLog << warn << GetDescriptor() << ": Less than 20 pixels with valid time resolution found "
365 << "in area index: " << aidx << endl;
366 continue;
367 }
368
369 // Calculate the rms out of sum2:
370 /*
371 areasum2[aidx] = (areasum2[aidx] - areasum[aidx]*areasum[aidx]/numareavalid[aidx]);
372 areasum2[aidx] /= (numareavalid[aidx]-1.);
373 */
374 areasum [aidx] /= numareavalid[aidx];
375 lowlim [aidx] = 0.;
376 upplim [aidx] = areasum [aidx] + fRelTimeResolutionLimit;
377
378 }
379 *fLog << endl;
380
381
382 for (UInt_t i=0; i<npixels; i++)
383 {
384
385 MCalibrationRelTimePix &pix = (MCalibrationRelTimePix&)(*relcam)[i];
386 MBadPixelsPix &bad = (*badcam)[i];
387
388 if (pix.IsExcluded())
389 continue;
390
391 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
392 continue;
393
394 const Float_t res = pix.GetTimePrecision();
395 const Int_t aidx = (*fGeom)[i].GetAidx();
396
397 if ( res < lowlim[aidx] || res > upplim[aidx] )
398 {
399 *fLog << warn << "Deviating time resolution: "
400 << Form("%4.2f",res) << " out of range ["
401 << Form("%4.2f,%4.2f",lowlim[aidx],upplim[aidx]) << "] in pixel " << i << endl;
402 bad.SetUncalibrated( MBadPixelsPix::kDeviatingTimeResolution);
403 pix.SetExcluded();
404 }
405 }
406}
407
408
409// -----------------------------------------------------------------------------------
410//
411// Sets pixel to MBadPixelsPix::kUnsuitableRun, if one of the following flags is set:
412// - MBadPixelsPix::kRelTimeIsPedestal
413// - MBadPixelsPix::kRelTimeErrNotValid
414// - MBadPixelsPix::kRelTimeRelErrNotValid
415//
416// - Call MCalibrationPix::SetExcluded() for the bad pixels
417//
418void MCalibrationRelTimeCalc::FinalizeBadPixels()
419{
420
421 MCalibrationRelTimeCam *relcam = fIntensCam
422 ? (MCalibrationRelTimeCam*)fIntensCam->GetCam() : fCam;
423 MBadPixelsCam *badcam = fIntensBad
424 ? (MBadPixelsCam*) fIntensBad->GetCam() : fBadPixels;
425
426 for (Int_t i=0; i<badcam->GetSize(); i++)
427 {
428
429 MBadPixelsPix &bad = (*badcam)[i];
430 MCalibrationRelTimePix &pix = (MCalibrationRelTimePix&)(*relcam)[i];
431
432 if (IsCheckDeviatingBehavior())
433 if (bad.IsUncalibrated(MBadPixelsPix::kDeviatingTimeResolution))
434 bad.SetUnsuitable(MBadPixelsPix::kUnreliableRun);
435
436 if (IsCheckFitResults())
437 if (bad.IsUncalibrated(MBadPixelsPix::kRelTimeNotFitted))
438 bad.SetUnsuitable(MBadPixelsPix::kUnreliableRun);
439
440 if (IsCheckOscillations())
441 if (bad.IsUncalibrated(MBadPixelsPix::kRelTimeOscillating))
442 bad.SetUnsuitable(MBadPixelsPix::kUnreliableRun);
443
444 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
445 pix.SetExcluded();
446 }
447
448}
449
450
451// -----------------------------------------------------------------------------------------------
452//
453// - Print out statistics about BadPixels of type UnsuitableType_t
454// - store numbers of bad pixels of each type in fIntensCam or fCam
455//
456void MCalibrationRelTimeCalc::FinalizeUnsuitablePixels()
457{
458 *fLog << inf << endl;
459 *fLog << GetDescriptor() << ": Rel. Times Calibration status:" << endl;
460 *fLog << dec;
461
462 MCalibrationRelTimeCam *relcam = fIntensCam ? (MCalibrationRelTimeCam*)fIntensCam->GetCam() : fCam;
463 const MBadPixelsCam *badcam = fIntensBad ? (MBadPixelsCam*)fIntensBad->GetCam() : fBadPixels;
464
465 const Int_t nareas = fGeom->GetNumAreas();
466
467 TArrayI unsuit(nareas);
468 TArrayI unrel(nareas);
469
470 for (int aidx=0; aidx<nareas; aidx++)
471 {
472 unsuit[aidx] += badcam->GetNumUnsuitable(MBadPixelsPix::kUnsuitableRun, fGeom, aidx);
473 unrel[aidx] += badcam->GetNumUnsuitable(MBadPixelsPix::kUnreliableRun, fGeom, aidx);
474 relcam->SetNumUnsuitable(unsuit[aidx], aidx);
475 relcam->SetNumUnreliable(unrel[aidx], aidx);
476 }
477
478 if (fGeom->InheritsFrom("MGeomCamMagic"))
479 {
480 PrintUncalibrated("Uncalibrated Pixels:", unsuit[0], unsuit[1]);
481 PrintUncalibrated("Unreliable Pixels:", unrel[0], unrel[1]);
482 }
483}
484
485// -----------------------------------------------------------------------------------------------
486//
487// Print out statistics about BadPixels of type UncalibratedType_t
488//
489void MCalibrationRelTimeCalc::PrintUncalibrated(MBadPixelsPix::UncalibratedType_t typ, const char *text) const
490{
491 const MBadPixelsCam *badcam = fIntensBad ? (MBadPixelsCam*)fIntensBad->GetCam() : fBadPixels;
492
493 UInt_t countinner = 0;
494 UInt_t countouter = 0;
495 for (Int_t i=0; i<badcam->GetSize(); i++)
496 {
497 if ((*badcam)[i].IsUncalibrated(typ))
498 {
499 if (fGeom->GetPixRatio(i) == 1.)
500 countinner++;
501 else
502 countouter++;
503 }
504 }
505
506 PrintUncalibrated(text, countinner, countouter);
507}
508
509void MCalibrationRelTimeCalc::PrintUncalibrated(const char *text, Int_t in, Int_t out) const
510{
511 *fLog << " " << setfill(' ') << setw(56) << setiosflags(ios::left) << text;
512 *fLog << " Inner: " << Form("%3i", in);
513 *fLog << " Outer: " << Form("%3i", out);
514 *fLog << resetiosflags(ios::left) << endl;
515
516}
517
518// --------------------------------------------------------------------------
519//
520// MCalibrationRelTimeCam.CheckFitResults: Yes
521// MCalibrationRelTimeCam.CheckDeviatingBehavior: Yes
522// MCalibrationRelTimeCam.CheckHistOverflow: Yes
523// MCalibrationRelTimeCam.CheckOscillations: Yes
524//
525Int_t MCalibrationRelTimeCalc::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
526{
527 Bool_t rc = kFALSE;
528
529 if (IsEnvDefined(env, prefix, "CheckFitResults", print))
530 {
531 SetCheckFitResults(GetEnvValue(env, prefix, "CheckFitResults", IsCheckFitResults()));
532 rc = kTRUE;
533 }
534 if (IsEnvDefined(env, prefix, "CheckDeviatingBehavior", print))
535 {
536 SetCheckDeviatingBehavior(GetEnvValue(env, prefix, "CheckDeviatingBehavior", IsCheckDeviatingBehavior()));
537 rc = kTRUE;
538 }
539 if (IsEnvDefined(env, prefix, "CheckHistOverflow", print))
540 {
541 SetCheckHistOverflow(GetEnvValue(env, prefix, "CheckHistOverflow", IsCheckHistOverflow()));
542 rc = kTRUE;
543 }
544
545 if (IsEnvDefined(env, prefix, "CheckOscillations", print))
546 {
547 SetCheckOscillations(GetEnvValue(env, prefix, "CheckOscillations", IsCheckOscillations()));
548 rc = kTRUE;
549 }
550 if (IsEnvDefined(env, prefix, "RelTimeResolutionLimit", print))
551 {
552 SetRelTimeResolutionLimit(GetEnvValue(env, prefix, "RelTimeResolutionLimit", fRelTimeResolutionLimit));
553 rc = kTRUE;
554 }
555
556 return rc;
557}
Note: See TracBrowser for help on using the repository browser.