source: trunk/MagicSoft/Mars/mcalib/MCalibrationTestCalc.cc@ 5710

Last change on this file since 5710 was 5346, checked in by gaug, 20 years ago
*** empty log message ***
File size: 20.3 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 08/2004 <mailto:markus@ifae.es>
19!
20! Copyright: MAGIC Software Development, 2000-2004
21!
22!
23\* ======================================================================== */
24
25//////////////////////////////////////////////////////////////////////////////
26//
27// MCalibrationTestCalc
28//
29// PreProcess(): Initialize pointers to MHCalibrationTestCam
30//
31// ReInit(): MCalibrationCam::InitSize(NumPixels) is called from MGeomApply (which allocates
32// memory in a TClonesArray of type MCalibrationChargePix)
33// Initializes pointer to MBadPixelsCam
34//
35// Process(): Nothing to be done, histograms getting filled by MHCalibrationTestCam
36//
37// PostProcess(): Print out interpolation results to file
38//
39// Input Containers:
40// MHCalibrationTestCam
41// MBadPixelsCam
42// MGeomCam
43//
44// Output Containers:
45// none
46//
47//////////////////////////////////////////////////////////////////////////////
48#include "MCalibrationTestCalc.h"
49
50#include <TSystem.h>
51#include <TH1.h>
52#include <TF1.h>
53
54#include "MLog.h"
55#include "MLogManip.h"
56
57#include "MParList.h"
58
59#include "MGeomCam.h"
60#include "MGeomPix.h"
61#include "MHCamera.h"
62
63#include "MHCalibrationTestCam.h"
64#include "MHCalibrationPix.h"
65
66#include "MCalibrationIntensityTestCam.h"
67#include "MCalibrationTestCam.h"
68#include "MCalibrationTestPix.h"
69
70#include "MBadPixelsIntensityCam.h"
71#include "MBadPixelsCam.h"
72#include "MBadPixelsPix.h"
73
74ClassImp(MCalibrationTestCalc);
75
76using namespace std;
77
78const Float_t MCalibrationTestCalc::fgPhotErrLimit = 4.5;
79// --------------------------------------------------------------------------
80//
81// Default constructor.
82//
83// Sets the pointer to fTestCam and fGeom to NULL
84// Sets outputpath to "."
85// Sets outputfile to "TestCalibStat.txt"
86// Sets fPhotErrLimit to fgPhotErrLimit
87//
88// Calls:
89// - Clear()
90//
91MCalibrationTestCalc::MCalibrationTestCalc(const char *name, const char *title)
92 : fIntensBad(NULL), fBadPixels(NULL),
93 fTestCam(NULL), fIntensCam(NULL), fCam(NULL),
94 fGeom(NULL)
95{
96
97 fName = name ? name : "MCalibrationTestCalc";
98 fTitle = title ? title : "Task to output the results of MHCalibrationTestCam ";
99
100 SetPhotErrLimit();
101
102 SetOutputPath();
103 SetOutputFile();
104
105}
106
107
108// --------------------------------------------------------------------------
109//
110// Search for the following input containers and abort if not existing:
111// - MGeomCam
112// - MHCalibrationTestCam
113// - MCalibrationIntensityTestCam or MCalibrationTestCam
114// - MBadPixelsIntensityCam or MBadPixelsCam
115//
116Bool_t MCalibrationTestCalc::ReInit(MParList *pList )
117{
118
119 fGeom = (MGeomCam*)pList->FindObject("MGeomCam");
120 if (!fGeom)
121 {
122 *fLog << err << "No MGeomCam found... aborting." << endl;
123 return kFALSE;
124 }
125
126 fTestCam = (MHCalibrationTestCam*)pList->FindObject("MHCalibrationTestCam");
127 if (!fTestCam)
128 {
129 *fLog << err << "Cannot find MHCalibrationTestCam... aborting" << endl;
130 *fLog << err << "Maybe you forget to call an MFillH for the MHCalibrationTestCam before..." << endl;
131 return kFALSE;
132 }
133
134 fIntensCam = (MCalibrationIntensityTestCam*)pList->FindObject(AddSerialNumber("MCalibrationIntensityTestCam"));
135 if (fIntensCam)
136 *fLog << inf << "Found MCalibrationIntensityTestCam ... " << endl;
137 else
138 {
139 fCam = (MCalibrationTestCam*)pList->FindObject(AddSerialNumber("MCalibrationTestCam"));
140 if (!fCam)
141 {
142 *fLog << err << "Cannot find MCalibrationTestCam ... abort." << endl;
143 *fLog << err << "Maybe you forget to call an MFillH for the MHCalibrationTestCam before..." << endl;
144 return kFALSE;
145 }
146 }
147
148 fIntensBad = (MBadPixelsIntensityCam*)pList->FindObject(AddSerialNumber("MBadPixelsIntensityCam"));
149 if (fIntensBad)
150 *fLog << inf << "Found MBadPixelsIntensityCam ... " << endl;
151 else
152 {
153 fBadPixels = (MBadPixelsCam*)pList->FindObject(AddSerialNumber("MBadPixelsCam"));
154 if (!fBadPixels)
155 {
156 *fLog << err << "Cannot find MBadPixelsCam ... abort." << endl;
157 return kFALSE;
158 }
159 }
160
161
162 return kTRUE;
163}
164
165// ----------------------------------------------------------------------------------
166//
167// Nothing to be done in Process, but have a look at MHCalibrationTestCam, instead
168//
169Int_t MCalibrationTestCalc::Process()
170{
171 return kTRUE;
172}
173
174// -----------------------------------------------------------------------
175//
176// Return if number of executions is null.
177//
178// Print out some statistics
179//
180Int_t MCalibrationTestCalc::PostProcess()
181{
182
183 if (GetNumExecutions()==0)
184 return kFALSE;
185
186 MCalibrationTestCam *testcam = fIntensCam
187 ? (MCalibrationTestCam*)fIntensCam->GetCam() : fCam;
188
189 //
190 // Re-direct the output to an ascii-file from now on:
191 //
192 MLog asciilog;
193 asciilog.SetOutputFile(GetOutputFile(),kTRUE);
194 SetLogStream(&asciilog);
195 //
196 // Finalize calibration statistics
197 //
198 FinalizeCalibratedPhotons();
199 FinalizeNotInterpolated();
200 const Int_t maxbad = CalcMaxNumBadPixelsCluster();
201
202
203 *fLog << inf << endl;
204 *fLog << GetDescriptor() << ": Pixel Interpolation status:" << endl;
205
206 if (fGeom->InheritsFrom("MGeomCamMagic"))
207 {
208 *fLog << " " << setw(7) << "Not interpolateable Pixels: "
209 << Form("%s%3i%s%3i","Inner: ",testcam->GetNumUninterpolated(0),
210 " Outer: ",testcam->GetNumUninterpolated(1)) << endl;
211 *fLog << " " << setw(7) << "Biggest not-interpolateable cluster: "
212 << maxbad << endl;
213 }
214
215 testcam->SetNumUninterpolatedInMaxCluster(maxbad);
216
217 *fLog << endl;
218 SetLogStream(&gLog);
219
220 return kTRUE;
221}
222
223
224// ------------------------------------------------------------------------
225//
226//
227// First loop: Calculate a mean and mean RMS of calibrated photons per area index
228// Include only MHCalibrationTestPix's which are not empty (not interpolated)
229//
230// Second loop: Get weighted mean number of calibrated photons and its RMS
231// excluding those deviating by more than fPhotErrLimit
232// sigmas from the mean (obtained in first loop). Set
233// MBadPixelsPix::kDeviatingNumPhots if excluded.
234//
235void MCalibrationTestCalc::FinalizeCalibratedPhotons() const
236{
237
238 const UInt_t npixels = fGeom->GetNumPixels();
239 const UInt_t nareas = fGeom->GetNumAreas();
240 const UInt_t nsectors = fGeom->GetNumSectors();
241
242 MCalibrationTestCam *testcam = fIntensCam
243 ? (MCalibrationTestCam*)fIntensCam->GetCam() : fCam;
244 MBadPixelsCam *badcam = fIntensBad
245 ? (MBadPixelsCam*) fIntensBad->GetCam() : fBadPixels;
246
247 TArrayD lowlim (nareas);
248 TArrayD upplim (nareas);
249 TArrayD areaphotons (nareas);
250 TArrayD sectorphotons(nsectors);
251 TArrayD areavars (nareas);
252 TArrayD sectorvars (nsectors);
253 TArrayD fittedmean (nareas);
254 TArrayD fittedsigma (nareas);
255 TArrayI numareavalid(nareas);
256 TArrayI numsectorvalid(nsectors);
257
258 //
259 // First loop: Get mean number of calibrated photons and the RMS
260 // The loop is only to recognize later pixels with very deviating numbers
261 //
262 MHCamera camphotons(*fGeom,"Camphotons","Photons in Camera");
263
264 for (UInt_t i=0; i<npixels; i++)
265 {
266
267 MHCalibrationPix &hist = (*fTestCam)[i];
268 MCalibrationTestPix &pix = (MCalibrationTestPix&)(*testcam)[i];
269 //
270 // We assume that the pixels have been interpolated so far.
271 // The MBadPixelsCam does not give any more information
272 //
273 if (hist.IsEmpty())
274 {
275 pix.SetExcluded();
276 continue;
277 }
278
279 const Double_t nphot = hist.GetMean();
280 const Double_t nphoterr = hist.GetMeanErr();
281 const Int_t aidx = (*fGeom)[i].GetAidx();
282
283 camphotons.Fill(i,nphot);
284 camphotons.SetUsed(i);
285
286 pix.SetNumPhotons ( nphot );
287 pix.SetNumPhotonsErr( nphoterr );
288
289 areaphotons [aidx] += nphot;
290 areavars [aidx] += nphot*nphot;
291 numareavalid[aidx] ++;
292 }
293
294 for (UInt_t aidx=0; aidx<nareas; aidx++)
295 {
296 if (numareavalid[aidx] == 0)
297 {
298 *fLog << warn << GetDescriptor() << ": No pixels with valid number of calibrated photons found "
299 << "in area index: " << aidx << endl;
300 continue;
301 }
302
303 if (numareavalid[aidx] == 1)
304 areavars[aidx] = 0.;
305 else if (numareavalid[aidx] == 0)
306 {
307 areaphotons[aidx] = -1.;
308 areavars[aidx] = -1.;
309 }
310 else
311 {
312 areavars[aidx] = (areavars[aidx] - areaphotons[aidx]*areaphotons[aidx]/numareavalid[aidx])
313 / (numareavalid[aidx]-1.);
314 areaphotons[aidx] = areaphotons[aidx] / numareavalid[aidx];
315 }
316
317 if (areavars[aidx] < 0.)
318 {
319 *fLog << warn << GetDescriptor() << ": No pixels with valid variance of calibrated photons found "
320 << "in area index: " << aidx << endl;
321 continue;
322 }
323
324 const Float_t areamean = areaphotons[aidx];
325 const Float_t areaerr = TMath::Sqrt(areavars[aidx]);
326
327 MCalibrationTestPix &avpix = (MCalibrationTestPix&)testcam->GetAverageArea(aidx);
328 avpix.SetNumPhotons (areamean);
329 avpix.SetNumPhotonsErr(areaerr );
330
331 lowlim [aidx] = areamean - fPhotErrLimit*areaerr;
332 upplim [aidx] = areamean + fPhotErrLimit*areaerr;
333
334 TArrayI area(1);
335 area[0] = aidx;
336
337 TH1D *hist = camphotons.ProjectionS(TArrayI(),area,"_py",100);
338 hist->Fit("gaus","Q");
339 const Double_t mean = hist->GetFunction("gaus")->GetParameter(1);
340 const Double_t sigma = hist->GetFunction("gaus")->GetParameter(2);
341 const Int_t ndf = hist->GetFunction("gaus")->GetNDF();
342
343 if (ndf < 2)
344 {
345 *fLog << warn << GetDescriptor() << ": Cannot use a Gauss fit to the number of calibrated photons "
346 << "in the camera with area index: " << aidx << endl;
347 *fLog << warn << GetDescriptor() << ": Number of dof.: " << ndf << " is smaller than 2 " << endl;
348 *fLog << warn << GetDescriptor() << ": Will use the simple mean and rms " << endl;
349 delete hist;
350 continue;
351 }
352
353 const Double_t prob = hist->GetFunction("gaus")->GetProb();
354
355 if (prob < 0.001)
356 {
357 *fLog << warn << GetDescriptor() << ": Cannot use a Gauss fit to the number of calibrated photons "
358 << "in the camera with area index: " << aidx << endl;
359 *fLog << warn << GetDescriptor() << ": Fit probability " << prob
360 << " is smaller than 0.001 " << endl;
361 *fLog << warn << GetDescriptor() << ": Will use the simple mean and rms " << endl;
362 delete hist;
363 continue;
364 }
365
366 fittedmean [aidx] = mean;
367 fittedsigma[aidx] = sigma;
368
369 avpix.SetNumPhotons (mean );
370 avpix.SetNumPhotonsErr(sigma);
371
372 lowlim [aidx] = mean - fPhotErrLimit*sigma;
373 upplim [aidx] = mean + fPhotErrLimit*sigma;
374
375 *fLog << inf << GetDescriptor()
376 << ": Fitted number of calib. equiv. Cher. photons in area index " << aidx
377 << ": " << Form("%7.2f +- %6.2f",fittedmean[aidx],fittedsigma[aidx]) << endl;
378
379 delete hist;
380 }
381
382 *fLog << endl;
383
384 numareavalid.Reset();
385 areaphotons .Reset();
386 areavars .Reset();
387
388 //
389 // Second loop: Get mean number of calibrated photons and its RMS excluding
390 // pixels deviating by more than fPhotErrLimit sigma.
391 //
392 for (UInt_t i=0; i<npixels; i++)
393 {
394
395 MHCalibrationPix &hist = (*fTestCam)[i];
396 MCalibrationTestPix &pix = (MCalibrationTestPix&) (*testcam)[i];
397
398 const Int_t aidx = (*fGeom)[i].GetAidx();
399 const Int_t sector = (*fGeom)[i].GetSector();
400 const Double_t nphot = hist.GetMean();
401 const Double_t nphotpera = nphot / (*fGeom)[i].GetA();
402 const Double_t nphotperaerr = hist.GetMeanErr()/ (*fGeom)[i].GetA();
403
404 pix.SetNumPhotonsPerArea ( nphotpera );
405 pix.SetNumPhotonsPerAreaErr( nphotperaerr );
406
407 if ( nphot < lowlim[aidx] || nphot > upplim[aidx] )
408 {
409 *fLog << warn << GetDescriptor() << ": Number of photons: "
410 << Form("%8.2f out of %3.1f sigma limit: ",nphot,fPhotErrLimit)
411 << Form("[%8.2f,%8.2f] pixel%4i",lowlim[aidx],upplim[aidx],i) << endl;
412 MBadPixelsPix &bad = (*badcam)[i];
413 bad.SetUncalibrated( MBadPixelsPix::kDeviatingNumPhots );
414 bad.SetUnsuitable ( MBadPixelsPix::kUnsuitableRun );
415 continue;
416 }
417
418 areavars [aidx] += nphotpera*nphotpera;
419 areaphotons [aidx] += nphotpera;
420 numareavalid [aidx] ++;
421
422 sectorvars [sector] += nphotpera*nphotpera;
423 sectorphotons [sector] += nphotpera;
424 numsectorvalid[sector] ++;
425 }
426
427 *fLog << endl;
428
429 for (UInt_t aidx=0; aidx<nareas; aidx++)
430 {
431
432 if (numareavalid[aidx] == 1)
433 areavars[aidx] = 0.;
434 else if (numareavalid[aidx] == 0)
435 {
436 areaphotons[aidx] = -1.;
437 areavars[aidx] = -1.;
438 }
439 else
440 {
441 areavars[aidx] = (areavars[aidx] - areaphotons[aidx]*areaphotons[aidx]/numareavalid[aidx])
442 / (numareavalid[aidx]-1.);
443 areaphotons[aidx] /= numareavalid[aidx];
444 }
445
446
447 MCalibrationTestPix &avpix = (MCalibrationTestPix&)testcam->GetAverageArea(aidx);
448
449 if (areavars[aidx] < 0. || areaphotons[aidx] <= 0.)
450 {
451 *fLog << warn << GetDescriptor()
452 << ": Mean number of photons per area in area index "
453 << aidx << " could not be calculated! Mean: " << areaphotons[aidx]
454 << " Variance: " << areavars[aidx] << endl;
455 avpix.SetExcluded();
456 continue;
457 }
458
459 const Float_t areaerr = TMath::Sqrt(areavars[aidx]);
460
461 avpix.SetNumPhotonsPerArea (areaphotons[aidx]);
462 avpix.SetNumPhotonsPerAreaErr(areaerr );
463
464 *fLog << inf << GetDescriptor() << ": Mean number of equiv. Cher. photons "
465 << "per area in area idx " << aidx << ": "
466 << Form("%5.3f+-%5.4f [ph/mm^2]",areaphotons[aidx],areaerr) << endl;
467 }
468
469 *fLog << endl;
470
471 for (UInt_t sector=0; sector<nsectors; sector++)
472 {
473
474 if (numsectorvalid[sector] == 1)
475 sectorvars[sector] = 0.;
476 else if (numsectorvalid[sector] == 0)
477 {
478 sectorphotons[sector] = -1.;
479 sectorvars[sector] = -1.;
480 }
481 else
482 {
483 sectorvars[sector] = (sectorvars[sector]
484 - sectorphotons[sector]*sectorphotons[sector]/numsectorvalid[sector]
485 )
486 / (numsectorvalid[sector]-1.);
487 sectorphotons[sector] /= numsectorvalid[sector];
488 }
489
490 MCalibrationTestPix &avpix = (MCalibrationTestPix&)testcam->GetAverageSector(sector);
491
492 if (sectorvars[sector] < 0. || sectorphotons[sector] <= 0.)
493 {
494 *fLog << warn << GetDescriptor()
495 << ": Mean number of calibrated photons per area in sector "
496 << sector << " could not be calculated! Mean: " << sectorphotons[sector]
497 << " Variance: " << sectorvars[sector] << endl;
498 avpix.SetExcluded();
499 continue;
500 }
501
502
503 const Float_t sectorerr = TMath::Sqrt(sectorvars[sector]);
504
505 avpix.SetNumPhotonsPerArea (sectorphotons[sector]);
506 avpix.SetNumPhotonsPerAreaErr(sectorerr );
507
508 *fLog << inf << GetDescriptor() << ": Mean number of equiv. Cher. photons "
509 << "per area in sector " << sector << ": "
510 << Form("%5.3f+-%5.4f [ph/mm^2]",sectorphotons[sector],sectorerr) << endl;
511 }
512
513 return;
514}
515
516
517// -----------------------------------------------------------------------------------------------
518//
519// Print out statistics about not interpolated pixels
520//
521void MCalibrationTestCalc::FinalizeNotInterpolated()
522{
523
524 MCalibrationTestCam *testcam = fIntensCam
525 ? (MCalibrationTestCam*)fIntensCam->GetCam() : fCam;
526
527 const Int_t areas = fGeom->GetNumAreas();
528 TArrayI *newarr[areas];
529
530 for (Int_t aidx=0; aidx<areas; aidx++)
531 newarr[aidx] = new TArrayI(0);
532
533 for (Int_t i=0; i<testcam->GetSize(); i++)
534 {
535 const Int_t aidx = (*fGeom)[i].GetAidx();
536 MCalibrationTestPix &pix = (MCalibrationTestPix&)(*testcam)[i];
537 if (pix.IsExcluded())
538 {
539 const Int_t size = newarr[aidx]->GetSize();
540 newarr[aidx]->Set(size+1);
541 newarr[aidx]->AddAt(i,size);
542 }
543 }
544
545 Int_t num = 0;
546
547 for (Int_t aidx = 0; aidx<areas; aidx++)
548 {
549 *fLog << endl;
550 *fLog << " " << setw(7)
551 << "Not interpolated pixels by in area idx " << aidx << ": ";
552 for (Int_t i=0; i<newarr[aidx]->GetSize(); i++)
553 {
554 *fLog << newarr[aidx]->At(i) << " ";
555 num++;
556 }
557 testcam->SetNumUninterpolated(newarr[aidx]->GetSize(),aidx);
558 *fLog << endl;
559 }
560
561 *fLog << " " << setw(7) << num << " total not interpolateable pixels " << endl;
562
563}
564
565Int_t MCalibrationTestCalc::CalcMaxNumBadPixelsCluster()
566{
567
568 MCalibrationTestCam *testcam = fIntensCam
569 ? (MCalibrationTestCam*)fIntensCam->GetCam() : fCam;
570
571 TArrayI arr(0);
572
573 for (Int_t i=0; i<testcam->GetSize(); i++)
574 {
575 MCalibrationTestPix &pix = (MCalibrationTestPix&)(*testcam)[i];
576 if (pix.IsExcluded())
577 {
578 const Int_t s = arr.GetSize();
579 arr.Set(s+1);
580 arr[s] = i;
581 }
582 }
583
584 const Int_t size = arr.GetSize();
585
586 if (size == 0)
587 return 0;
588
589 if (size == 1)
590 return 1;
591
592 TArrayI knownpixels(0);
593 Int_t clustersize = 1;
594 Int_t oldclustersize = 0;
595 //
596 // Loop over the not-interpolateable pixels:
597 //
598 for (Int_t i=0; i<size; i++)
599 {
600
601 const Int_t id = arr[i];
602 const Int_t knownsize = knownpixels.GetSize();
603 knownpixels.Set(knownsize+1);
604 knownpixels[knownsize] = id;
605 LoopNeighbours(arr, knownpixels, clustersize, id);
606 if (clustersize > oldclustersize)
607 oldclustersize = clustersize;
608 clustersize = 1;
609 }
610
611 return oldclustersize;
612
613}
614
615
616void MCalibrationTestCalc::LoopNeighbours( const TArrayI &arr, TArrayI &knownpixels, Int_t &clustersize, const Int_t idx )
617{
618
619 const MGeomPix &pix = (*fGeom)[idx];
620 const Byte_t neighbours = pix.GetNumNeighbors();
621
622 //
623 // Loop over the next neighbours:
624 // - Check if they are already in the list of known pixels
625 // - If not, call loopneighbours for the rest
626 // - grow clustersize for those
627 //
628 for (Int_t i=0;i<neighbours;i++)
629 {
630 const Int_t newid = pix.GetNeighbor(i);
631 Bool_t known = kFALSE;
632
633 for (Int_t j=knownpixels.GetSize()-1;j>=0;j--)
634 if (newid == knownpixels.At(j))
635 {
636 known = kTRUE;
637 break;
638 }
639 if (known)
640 continue;
641
642 for (Int_t k=0;k<arr.GetSize();k++)
643 if (newid == arr.At(k))
644 {
645 // This is an unknown, new pixel in the cluster!!
646 clustersize++;
647 const Int_t knownsize = knownpixels.GetSize();
648 knownpixels.Set(knownsize+1);
649 knownpixels[knownsize] = newid;
650 LoopNeighbours(arr, knownpixels, clustersize, newid);
651 }
652 }
653}
654
655// --------------------------------------------------------------------------
656//
657// Set the path for output file
658//
659void MCalibrationTestCalc::SetOutputPath(TString path)
660{
661 fOutputPath = path;
662 if (fOutputPath.EndsWith("/"))
663 fOutputPath = fOutputPath(0, fOutputPath.Length()-1);
664}
665
666void MCalibrationTestCalc::SetOutputFile(TString file)
667{
668 fOutputFile = file;
669}
670
671// --------------------------------------------------------------------------
672//
673// Get the output file
674//
675const char* MCalibrationTestCalc::GetOutputFile()
676{
677 return Form("%s/%s", (const char*)fOutputPath, (const char*)fOutputFile);
678}
679
Note: See TracBrowser for help on using the repository browser.