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

Last change on this file since 3917 was 3917, checked in by gaug, 20 years ago
*** empty log message ***
File size: 14.4 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// Input Containers:
43// MCalibrationRelTimeCam
44// MBadPixelsCam
45// MGeomCam
46//
47// Output Containers:
48// MCalibrationRelTimeCam
49// MBadPixelsCam
50//
51//
52//////////////////////////////////////////////////////////////////////////////
53#include "MCalibrationRelTimeCalc.h"
54
55#include "MLog.h"
56#include "MLogManip.h"
57
58#include "MParList.h"
59
60#include "MGeomCam.h"
61#include "MGeomPix.h"
62
63#include "MCalibrationRelTimeCam.h"
64#include "MCalibrationRelTimePix.h"
65
66#include "MBadPixelsCam.h"
67#include "MBadPixelsPix.h"
68
69
70ClassImp(MCalibrationRelTimeCalc);
71
72using namespace std;
73
74const Float_t MCalibrationRelTimeCalc::fgRelTimeRelErrLimit = 3.;
75// --------------------------------------------------------------------------
76//
77// Default constructor.
78//
79// Sets all pointers to NULL
80//
81// Initializes:
82// - fRelTimeRelErrLimit to fgRelTimeRelErrLimit
83//
84// Calls:
85// - Clear()
86//
87MCalibrationRelTimeCalc::MCalibrationRelTimeCalc(const char *name, const char *title)
88 : fBadPixels(NULL), fCam(NULL), fGeom(NULL)
89{
90
91 fName = name ? name : "MCalibrationRelTimeCalc";
92 fTitle = title ? title : "Task to finalize the relative time calibration";
93
94 SetRelTimeRelErrLimit();
95
96 Clear();
97
98}
99
100// --------------------------------------------------------------------------
101//
102// Sets:
103// - all flags to kFALSE
104//
105void MCalibrationRelTimeCalc::Clear(const Option_t *o)
106{
107 SkipHiLoGainCalibration( kFALSE );
108}
109
110
111// -----------------------------------------------------------------------------------
112//
113// The following containers are searched and created if they were not found:
114//
115// - MBadPixelsCam
116//
117Int_t MCalibrationRelTimeCalc::PreProcess(MParList *pList)
118{
119
120
121 fBadPixels = (MBadPixelsCam*)pList->FindCreateObj("MBadPixelsCam");
122 if (!fBadPixels)
123 {
124 *fLog << err << "Could not find or create MBadPixelsCam ... aborting." << endl;
125 return kFALSE;
126 }
127
128
129
130 return kTRUE;
131}
132
133
134// --------------------------------------------------------------------------
135//
136// Search for the following input containers and abort if not existing:
137// - MGeomCam
138// - MCalibrationRelTimeCam
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 fCam = (MCalibrationRelTimeCam*)pList->FindObject("MCalibrationRelTimeCam");
159 if (!fCam)
160 {
161 *fLog << err << "Cannot find MCalibrationRelTimeCam... aborting" << endl;
162 *fLog << err << "Maybe you forget to call an MFillH for the MHCalibrationRelTimeCam before..." << endl;
163 return kFALSE;
164 }
165
166
167 UInt_t npixels = fGeom->GetNumPixels();
168
169 for (UInt_t i=0; i<npixels; i++)
170 {
171
172 MCalibrationRelTimePix &pix = (MCalibrationRelTimePix&)(*fCam) [i];
173 MBadPixelsPix &bad = (*fBadPixels)[i];
174
175 pix.SetPixId(i);
176
177 if (bad.IsBad())
178 {
179 pix.SetExcluded();
180 continue;
181 }
182
183 }
184
185 return kTRUE;
186}
187
188// ----------------------------------------------------------------------------------
189//
190// Nothing to be done in Process, but have a look at MHCalibrationRelTimeCam, instead
191//
192Int_t MCalibrationRelTimeCalc::Process()
193{
194 return kTRUE;
195}
196
197// -----------------------------------------------------------------------
198//
199// Return if number of executions is null.
200//
201// First loop over pixels, average areas and sectors, call:
202// - FinalizeRelTimes()
203// for every entry. Count number of valid pixels in loop and return kFALSE
204// if there are none (the "Michele check").
205//
206// Call FinalizeBadPixels()
207//
208// Call MParContainer::SetReadyToSave() for fCam
209//
210// Print out some statistics
211//
212Int_t MCalibrationRelTimeCalc::PostProcess()
213{
214
215 if (GetNumExecutions()==0)
216 return kFALSE;
217
218 //
219 // First loop over pixels, call FinalizePedestals and FinalizeRelTimes
220 //
221 Int_t nvalid = 0;
222
223 for (Int_t pixid=0; pixid<fCam->GetSize(); pixid++)
224 {
225
226 MCalibrationRelTimePix &pix = (MCalibrationRelTimePix&)(*fCam)[pixid];
227 //
228 // Check if the pixel has been excluded from the fits
229 //
230 if (pix.IsExcluded())
231 continue;
232
233 MBadPixelsPix &bad = (*fBadPixels)[pixid];
234
235 if (FinalizeRelTimes(pix,bad))
236 nvalid++;
237 }
238
239 //
240 // The Michele check ...
241 //
242 if (nvalid == 0)
243 {
244 *fLog << err << GetDescriptor() << ": All pixels have non-valid calibration. "
245 << "Did you forget to fill the histograms "
246 << "(filling MHCalibrationRelTimeCam from MArrivalTimeCam using MFillH) ? " << endl;
247 *fLog << err << GetDescriptor() << ": Or, maybe, you have used a pedestal run "
248 << "instead of a calibration run " << endl;
249 return kFALSE;
250 }
251
252 for (UInt_t aidx=0; aidx<fGeom->GetNumAreas(); aidx++)
253 {
254
255 MCalibrationRelTimePix &pix = (MCalibrationRelTimePix&)fCam->GetAverageArea(aidx);
256 FinalizeRelTimes(pix, fCam->GetAverageBadArea(aidx));
257 }
258
259 for (UInt_t sector=0; sector<fGeom->GetNumSectors(); sector++)
260 {
261
262 MCalibrationRelTimePix &pix = (MCalibrationRelTimePix&)fCam->GetAverageSector(sector);
263 FinalizeRelTimes(pix, fCam->GetAverageBadSector(sector));
264 }
265
266 //
267 // Finalize Bad Pixels
268 //
269 FinalizeBadPixels();
270
271 //
272 // Finalize calibration statistics
273 //
274 FinalizeUnsuitablePixels();
275
276 fCam ->SetReadyToSave();
277 fBadPixels->SetReadyToSave();
278
279 *fLog << inf << endl;
280 *fLog << GetDescriptor() << ": Errors statistics:" << endl;
281
282 PrintUncalibrated(MBadPixelsPix::kDeviatingTimeResolution,
283 Form("%s%2.1f%s","Time resolution less than ",fRelTimeRelErrLimit," sigma from Mean: "));
284 PrintUncalibrated(MBadPixelsPix::kRelTimeOscillating,
285 "Pixels with changing Rel. Times over time: ");
286 PrintUncalibrated(MBadPixelsPix::kRelTimeNotFitted,
287 "Pixels with unsuccesful Gauss fit to the times: ");
288 return kTRUE;
289}
290
291
292// ----------------------------------------------------------------------------------------------------
293//
294//
295// First loop: Calculate a mean and mean RMS of time resolution per area index
296// Include only pixels which are not MBadPixelsPix::kUnsuitableRun or
297// MBadPixelsPix::kUnreliableRun (see FinalizeBadPixels())
298//
299// Second loop: Exclude those deviating by more than fRelTimeRelErrLimit mean
300// sigmas from the mean (obtained in first loop). Set
301// MBadPixelsPix::kDeviatingTimeResolution if excluded.
302//
303Bool_t MCalibrationRelTimeCalc::FinalizeRelTimes(MCalibrationRelTimePix &cal, MBadPixelsPix &bad)
304{
305
306 const UInt_t npixels = fGeom->GetNumPixels();
307 const UInt_t nareas = fGeom->GetNumAreas();
308
309 Float_t lowlim [nareas];
310 Float_t upplim [nareas];
311 Float_t areaerrs [nareas];
312 Float_t areamean [nareas];
313 Int_t numareavalid[nareas];
314
315 memset(lowlim ,0, nareas * sizeof(Float_t));
316 memset(upplim ,0, nareas * sizeof(Float_t));
317 memset(areaerrs ,0, nareas * sizeof(Float_t));
318 memset(areamean ,0, nareas * sizeof(Float_t));
319 memset(numareavalid ,0, nareas * sizeof(Int_t ));
320
321 //
322 // First loop: Get mean time resolution the RMS
323 // The loop is only to recognize later pixels with very deviating numbers
324 //
325 for (UInt_t i=0; i<npixels; i++)
326 {
327
328 MCalibrationRelTimePix &pix = (MCalibrationRelTimePix&)(*fCam) [i];
329 MBadPixelsPix &bad = (*fBadPixels)[i];
330
331 if (!pix.IsExcluded())
332 continue;
333
334 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
335 continue;
336
337 if (bad.IsUnsuitable(MBadPixelsPix::kUnreliableRun))
338 continue;
339
340 const Float_t res = pix.GetTimePrecision();
341 const Float_t perr = pix.GetTimePrecisionErr();
342 const Int_t aidx = (*fGeom)[i].GetAidx();
343
344 areamean [aidx] += res;
345 areaerrs [aidx] += perr;
346 numareavalid[aidx] ++;
347 }
348
349
350
351 for (UInt_t i=0; i<nareas; i++)
352 {
353 if (numareavalid[i] == 0)
354 {
355 *fLog << warn << GetDescriptor() << ": No pixels with valid number of photo-electrons found "
356 << "in area index: " << i << endl;
357 continue;
358 }
359
360 areamean[i] = areamean[i] / numareavalid[i];
361 areaerrs[i] = areaerrs[i] / numareavalid[i];
362 lowlim [i] = areamean[i] - fRelTimeRelErrLimit*areaerrs[i];
363 upplim [i] = areamean[i] + fRelTimeRelErrLimit*areaerrs[i];
364 }
365
366
367
368 //
369 // Second loop: Exclude pixels deviating by more than fRelTimeErrLimit sigma.
370 //
371 for (UInt_t i=0; i<npixels; i++)
372 {
373
374 MCalibrationRelTimePix &pix = (MCalibrationRelTimePix&)(*fCam)[i];
375
376 if (!pix.IsExcluded())
377 continue;
378
379 const Float_t res = pix.GetTimePrecision();
380
381 MBadPixelsPix &bad = (*fBadPixels)[i];
382 const Int_t aidx = (*fGeom)[i].GetAidx();
383
384 if ( res < lowlim[aidx] || res > upplim[aidx] )
385 {
386 *fLog << warn << GetDescriptor() << ": Deviating time resolution: "
387 << Form("%4.2f",res) << " out of accepted limits: ["
388 << Form("%4.2f%s%4.2f",lowlim[aidx],",",upplim[aidx]) << "] in pixel " << i << endl;
389 bad.SetUncalibrated( MBadPixelsPix::kDeviatingTimeResolution);
390 bad.SetUnsuitable ( MBadPixelsPix::kUnsuitableRun );
391 pix.SetExcluded();
392 }
393 }
394
395 return kTRUE;
396}
397
398
399
400// -----------------------------------------------------------------------------------
401//
402// Sets pixel to MBadPixelsPix::kUnsuitableRun, if one of the following flags is set:
403// - MBadPixelsPix::kRelTimeIsPedestal
404// - MBadPixelsPix::kRelTimeErrNotValid
405// - MBadPixelsPix::kRelTimeRelErrNotValid
406//
407// - Call MCalibrationPix::SetExcluded() for the bad pixels
408//
409void MCalibrationRelTimeCalc::FinalizeBadPixels()
410{
411
412 for (Int_t i=0; i<fBadPixels->GetSize(); i++)
413 {
414
415 MBadPixelsPix &bad = (*fBadPixels)[i];
416 MCalibrationPix &pix = (*fCam)[i];
417
418 if (bad.IsUncalibrated( MBadPixelsPix::kDeviatingTimeResolution))
419 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );
420
421 if (bad.IsUncalibrated( MBadPixelsPix::kRelTimeNotFitted))
422 bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun );
423
424 if (bad.IsUncalibrated( MBadPixelsPix::kRelTimeOscillating))
425 bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun );
426
427 if (bad.IsUnsuitable( MBadPixelsPix::kUnsuitableRun ))
428 pix.SetExcluded();
429
430 }
431}
432
433
434// -----------------------------------------------------------------------------------------------
435//
436// - Print out statistics about BadPixels of type UnsuitableType_t
437// - store numbers of bad pixels of each type in fCam
438//
439void MCalibrationRelTimeCalc::FinalizeUnsuitablePixels()
440{
441
442 *fLog << inf << endl;
443 *fLog << GetDescriptor() << ": Calibration statistics:" << endl;
444 *fLog << dec << setfill(' ');
445
446 const Int_t nareas = fGeom->GetNumAreas();
447
448 Int_t counts[nareas];
449 memset(counts,0,nareas*sizeof(Int_t));
450
451 for (Int_t i=0; i<fBadPixels->GetSize(); i++)
452 {
453 MBadPixelsPix &bad = (*fBadPixels)[i];
454 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
455 {
456 const Int_t aidx = (*fGeom)[i].GetAidx();
457 counts[aidx]++;
458 }
459 }
460
461 for (Int_t aidx=0; aidx<nareas; aidx++)
462 fCam->SetNumUncalibrated(counts[aidx], aidx);
463
464 if (fGeom->InheritsFrom("MGeomCamMagic"))
465 *fLog << " " << setw(7) << "Uncalibrated Pixels: "
466 << Form("%s%3i%s%3i","Inner: ",counts[0]," Outer: ",counts[1]) << endl;
467
468 memset(counts,0,nareas*sizeof(Int_t));
469
470 for (Int_t i=0; i<fBadPixels->GetSize(); i++)
471 {
472 MBadPixelsPix &bad = (*fBadPixels)[i];
473 if (bad.IsUnsuitable(MBadPixelsPix::kUnreliableRun))
474 {
475 const Int_t aidx = (*fGeom)[i].GetAidx();
476 counts[aidx]++;
477 }
478 }
479
480 for (Int_t aidx=0; aidx<nareas; aidx++)
481 fCam->SetNumUnreliable(counts[aidx], aidx);
482
483 *fLog << " " << setw(7) << "Unreliable Pixels: "
484 << Form("%s%3i%s%3i","Inner: ",counts[0]," Outer: ",counts[1]) << endl;
485
486}
487
488// -----------------------------------------------------------------------------------------------
489//
490// Print out statistics about BadPixels of type UncalibratedType_t
491//
492void MCalibrationRelTimeCalc::PrintUncalibrated(MBadPixelsPix::UncalibratedType_t typ, const char *text) const
493{
494
495 UInt_t countinner = 0;
496 UInt_t countouter = 0;
497 for (Int_t i=0; i<fBadPixels->GetSize(); i++)
498 {
499 MBadPixelsPix &bad = (*fBadPixels)[i];
500 if (bad.IsUncalibrated(typ))
501 {
502 if (fGeom->GetPixRatio(i) == 1.)
503 countinner++;
504 else
505 countouter++;
506 }
507 }
508
509 *fLog << " " << setw(7) << text
510 << Form("%s%3i%s%3i","Inner: ",countinner," Outer: ",countouter) << endl;
511}
512
Note: See TracBrowser for help on using the repository browser.