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

Last change on this file since 8199 was 8192, checked in by tbretz, 18 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 << "Area " << setw(4) << aidx << ": Less than 20 pixels with valid time resolution found." << endl;
365 continue;
366 }
367
368 // Calculate the rms out of sum2:
369 /*
370 areasum2[aidx] = (areasum2[aidx] - areasum[aidx]*areasum[aidx]/numareavalid[aidx]);
371 areasum2[aidx] /= (numareavalid[aidx]-1.);
372 */
373 areasum [aidx] /= numareavalid[aidx];
374 lowlim [aidx] = 0.;
375 upplim [aidx] = areasum [aidx] + fRelTimeResolutionLimit;
376
377 }
378 *fLog << endl;
379
380
381 for (UInt_t i=0; i<npixels; i++)
382 {
383
384 MCalibrationRelTimePix &pix = (MCalibrationRelTimePix&)(*relcam)[i];
385 MBadPixelsPix &bad = (*badcam)[i];
386
387 if (pix.IsExcluded())
388 continue;
389
390 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
391 continue;
392
393 const Float_t res = pix.GetTimePrecision();
394 const Int_t aidx = (*fGeom)[i].GetAidx();
395
396 if ( res < lowlim[aidx] || res > upplim[aidx] )
397 {
398 *fLog << warn << "Pixel " << setw(4) << i << ": Deviating time resolution: "
399 << Form("%4.2f",res) << " out of range ["
400 << Form("%4.2f,%4.2f",lowlim[aidx],upplim[aidx]) << "]" << endl;
401 bad.SetUncalibrated( MBadPixelsPix::kDeviatingTimeResolution);
402 pix.SetExcluded();
403 }
404 }
405}
406
407
408// -----------------------------------------------------------------------------------
409//
410// Sets pixel to MBadPixelsPix::kUnsuitableRun, if one of the following flags is set:
411// - MBadPixelsPix::kRelTimeIsPedestal
412// - MBadPixelsPix::kRelTimeErrNotValid
413// - MBadPixelsPix::kRelTimeRelErrNotValid
414//
415// - Call MCalibrationPix::SetExcluded() for the bad pixels
416//
417void MCalibrationRelTimeCalc::FinalizeBadPixels()
418{
419
420 MCalibrationRelTimeCam *relcam = fIntensCam
421 ? (MCalibrationRelTimeCam*)fIntensCam->GetCam() : fCam;
422 MBadPixelsCam *badcam = fIntensBad
423 ? (MBadPixelsCam*) fIntensBad->GetCam() : fBadPixels;
424
425 for (Int_t i=0; i<badcam->GetSize(); i++)
426 {
427
428 MBadPixelsPix &bad = (*badcam)[i];
429 MCalibrationRelTimePix &pix = (MCalibrationRelTimePix&)(*relcam)[i];
430
431 if (IsCheckDeviatingBehavior())
432 if (bad.IsUncalibrated(MBadPixelsPix::kDeviatingTimeResolution))
433 bad.SetUnsuitable(MBadPixelsPix::kUnreliableRun);
434
435 if (IsCheckFitResults())
436 if (bad.IsUncalibrated(MBadPixelsPix::kRelTimeNotFitted))
437 bad.SetUnsuitable(MBadPixelsPix::kUnreliableRun);
438
439 if (IsCheckOscillations())
440 if (bad.IsUncalibrated(MBadPixelsPix::kRelTimeOscillating))
441 bad.SetUnsuitable(MBadPixelsPix::kUnreliableRun);
442
443 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
444 pix.SetExcluded();
445 }
446
447}
448
449
450// -----------------------------------------------------------------------------------------------
451//
452// - Print out statistics about BadPixels of type UnsuitableType_t
453// - store numbers of bad pixels of each type in fIntensCam or fCam
454//
455void MCalibrationRelTimeCalc::FinalizeUnsuitablePixels()
456{
457 *fLog << inf << endl;
458 *fLog << GetDescriptor() << ": Rel. Times Calibration status:" << endl;
459 *fLog << dec;
460
461 MCalibrationRelTimeCam *relcam = fIntensCam ? (MCalibrationRelTimeCam*)fIntensCam->GetCam() : fCam;
462 const MBadPixelsCam *badcam = fIntensBad ? (MBadPixelsCam*)fIntensBad->GetCam() : fBadPixels;
463
464 const Int_t nareas = fGeom->GetNumAreas();
465
466 TArrayI unsuit(nareas);
467 TArrayI unrel(nareas);
468
469 for (int aidx=0; aidx<nareas; aidx++)
470 {
471 unsuit[aidx] += badcam->GetNumUnsuitable(MBadPixelsPix::kUnsuitableRun, fGeom, aidx);
472 unrel[aidx] += badcam->GetNumUnsuitable(MBadPixelsPix::kUnreliableRun, fGeom, aidx);
473 relcam->SetNumUnsuitable(unsuit[aidx], aidx);
474 relcam->SetNumUnreliable(unrel[aidx], aidx);
475 }
476
477 if (fGeom->InheritsFrom("MGeomCamMagic"))
478 {
479 PrintUncalibrated("Uncalibrated Pixels:", unsuit[0], unsuit[1]);
480 PrintUncalibrated("Unreliable Pixels:", unrel[0], unrel[1]);
481 }
482}
483
484// -----------------------------------------------------------------------------------------------
485//
486// Print out statistics about BadPixels of type UncalibratedType_t
487//
488void MCalibrationRelTimeCalc::PrintUncalibrated(MBadPixelsPix::UncalibratedType_t typ, const char *text) const
489{
490 const MBadPixelsCam *badcam = fIntensBad ? (MBadPixelsCam*)fIntensBad->GetCam() : fBadPixels;
491
492 UInt_t countinner = 0;
493 UInt_t countouter = 0;
494 for (Int_t i=0; i<badcam->GetSize(); i++)
495 {
496 if ((*badcam)[i].IsUncalibrated(typ))
497 {
498 if (fGeom->GetPixRatio(i) == 1.)
499 countinner++;
500 else
501 countouter++;
502 }
503 }
504
505 PrintUncalibrated(text, countinner, countouter);
506}
507
508void MCalibrationRelTimeCalc::PrintUncalibrated(const char *text, Int_t in, Int_t out) const
509{
510 *fLog << " " << setfill(' ') << setw(56) << setiosflags(ios::left) << text;
511 *fLog << " Inner: " << Form("%3i", in);
512 *fLog << " Outer: " << Form("%3i", out);
513 *fLog << resetiosflags(ios::left) << endl;
514
515}
516
517// --------------------------------------------------------------------------
518//
519// MCalibrationRelTimeCam.CheckFitResults: Yes
520// MCalibrationRelTimeCam.CheckDeviatingBehavior: Yes
521// MCalibrationRelTimeCam.CheckHistOverflow: Yes
522// MCalibrationRelTimeCam.CheckOscillations: Yes
523//
524Int_t MCalibrationRelTimeCalc::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
525{
526 Bool_t rc = kFALSE;
527
528 if (IsEnvDefined(env, prefix, "CheckFitResults", print))
529 {
530 SetCheckFitResults(GetEnvValue(env, prefix, "CheckFitResults", IsCheckFitResults()));
531 rc = kTRUE;
532 }
533 if (IsEnvDefined(env, prefix, "CheckDeviatingBehavior", print))
534 {
535 SetCheckDeviatingBehavior(GetEnvValue(env, prefix, "CheckDeviatingBehavior", IsCheckDeviatingBehavior()));
536 rc = kTRUE;
537 }
538 if (IsEnvDefined(env, prefix, "CheckHistOverflow", print))
539 {
540 SetCheckHistOverflow(GetEnvValue(env, prefix, "CheckHistOverflow", IsCheckHistOverflow()));
541 rc = kTRUE;
542 }
543
544 if (IsEnvDefined(env, prefix, "CheckOscillations", print))
545 {
546 SetCheckOscillations(GetEnvValue(env, prefix, "CheckOscillations", IsCheckOscillations()));
547 rc = kTRUE;
548 }
549 if (IsEnvDefined(env, prefix, "RelTimeResolutionLimit", print))
550 {
551 SetRelTimeResolutionLimit(GetEnvValue(env, prefix, "RelTimeResolutionLimit", fRelTimeResolutionLimit));
552 rc = kTRUE;
553 }
554
555 return rc;
556}
Note: See TracBrowser for help on using the repository browser.