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

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