source: tags/Mars-V0.8.6/mcalib/MCalibrationTestCalc.cc

Last change on this file was 5046, checked in by gaug, 20 years ago
*** empty log message ***
File size: 20.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 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 "MHCalibrationTestPix.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 MHCalibrationTestPix &hist = (MHCalibrationTestPix&)(*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
280 const Double_t nphot = hist.GetMean();
281 const Double_t nphoterr = hist.GetMeanErr();
282 const Int_t aidx = (*fGeom)[i].GetAidx();
283
284 camphotons.Fill(i,nphot);
285 camphotons.SetUsed(i);
286
287 pix.SetNumPhotons ( nphot );
288 pix.SetNumPhotonsErr( nphoterr );
289
290 areaphotons [aidx] += nphot;
291 areavars [aidx] += nphot*nphot;
292 numareavalid[aidx] ++;
293 }
294
295 for (UInt_t aidx=0; aidx<nareas; aidx++)
296 {
297 if (numareavalid[aidx] == 0)
298 {
299 *fLog << warn << GetDescriptor() << ": No pixels with valid number of calibrated photons found "
300 << "in area index: " << aidx << endl;
301 continue;
302 }
303
304 if (numareavalid[aidx] == 1)
305 areavars[aidx] = 0.;
306 else if (numareavalid[aidx] == 0)
307 {
308 areaphotons[aidx] = -1.;
309 areavars[aidx] = -1.;
310 }
311 else
312 {
313 areavars[aidx] = (areavars[aidx] - areaphotons[aidx]*areaphotons[aidx]/numareavalid[aidx])
314 / (numareavalid[aidx]-1.);
315 areaphotons[aidx] = areaphotons[aidx] / numareavalid[aidx];
316 }
317
318 if (areavars[aidx] < 0.)
319 {
320 *fLog << warn << GetDescriptor() << ": No pixels with valid variance of calibrated photons found "
321 << "in area index: " << aidx << endl;
322 continue;
323 }
324
325 const Float_t areamean = areaphotons[aidx];
326 const Float_t areaerr = TMath::Sqrt(areavars[aidx]);
327
328 MCalibrationTestPix &avpix = (MCalibrationTestPix&)testcam->GetAverageArea(aidx);
329 avpix.SetNumPhotons (areamean);
330 avpix.SetNumPhotonsErr(areaerr );
331
332 lowlim [aidx] = areamean - fPhotErrLimit*areaerr;
333 upplim [aidx] = areamean + fPhotErrLimit*areaerr;
334
335 TArrayI area(1);
336 area[0] = aidx;
337
338 TH1D *hist = camphotons.ProjectionS(TArrayI(),area,"_py",100);
339 hist->Fit("gaus","Q");
340 const Double_t mean = hist->GetFunction("gaus")->GetParameter(1);
341 const Double_t sigma = hist->GetFunction("gaus")->GetParameter(2);
342 const Int_t ndf = hist->GetFunction("gaus")->GetNDF();
343
344 if (ndf < 2)
345 {
346 *fLog << warn << GetDescriptor() << ": Cannot use a Gauss fit to the number of calibrated photons "
347 << "in the camera with area index: " << aidx << endl;
348 *fLog << warn << GetDescriptor() << ": Number of dof.: " << ndf << " is smaller than 2 " << endl;
349 *fLog << warn << GetDescriptor() << ": Will use the simple mean and rms " << endl;
350 delete hist;
351 continue;
352 }
353
354 const Double_t prob = hist->GetFunction("gaus")->GetProb();
355
356 if (prob < 0.001)
357 {
358 *fLog << warn << GetDescriptor() << ": Cannot use a Gauss fit to the number of calibrated photons "
359 << "in the camera with area index: " << aidx << endl;
360 *fLog << warn << GetDescriptor() << ": Fit probability " << prob
361 << " is smaller than 0.001 " << endl;
362 *fLog << warn << GetDescriptor() << ": Will use the simple mean and rms " << endl;
363 delete hist;
364 continue;
365 }
366
367 fittedmean [aidx] = mean;
368 fittedsigma[aidx] = sigma;
369
370 avpix.SetNumPhotons (mean );
371 avpix.SetNumPhotonsErr(sigma);
372
373 lowlim [aidx] = mean - fPhotErrLimit*sigma;
374 upplim [aidx] = mean + fPhotErrLimit*sigma;
375
376 *fLog << inf << GetDescriptor()
377 << ": Fitted number of calib. equiv. Cher. photons in area index " << aidx
378 << ": " << Form("%7.2f +- %6.2f",fittedmean[aidx],fittedsigma[aidx]) << endl;
379
380 delete hist;
381 }
382
383 *fLog << endl;
384
385 numareavalid.Reset();
386 areaphotons .Reset();
387 areavars .Reset();
388
389 //
390 // Second loop: Get mean number of calibrated photons and its RMS excluding
391 // pixels deviating by more than fPhotErrLimit sigma.
392 //
393 for (UInt_t i=0; i<npixels; i++)
394 {
395
396 MHCalibrationTestPix &hist = (MHCalibrationTestPix&)(*fTestCam)[i];
397 MCalibrationTestPix &pix = (MCalibrationTestPix&) (*testcam)[i];
398
399 const Int_t aidx = (*fGeom)[i].GetAidx();
400 const Int_t sector = (*fGeom)[i].GetSector();
401 const Double_t nphot = hist.GetMean();
402 const Double_t nphotpera = nphot / (*fGeom)[i].GetA();
403 const Double_t nphotperaerr = hist.GetMeanErr()/ (*fGeom)[i].GetA();
404
405 pix.SetNumPhotonsPerArea ( nphotpera );
406 pix.SetNumPhotonsPerAreaErr( nphotperaerr );
407
408 if ( nphot < lowlim[aidx] || nphot > upplim[aidx] )
409 {
410 *fLog << warn << GetDescriptor() << ": Number of photons: "
411 << Form("%8.2f out of %3.1f sigma limit: ",nphot,fPhotErrLimit)
412 << Form("[%8.2f,%8.2f] pixel%4i",lowlim[aidx],upplim[aidx],i) << endl;
413 MBadPixelsPix &bad = (*badcam)[i];
414 bad.SetUncalibrated( MBadPixelsPix::kDeviatingNumPhots );
415 bad.SetUnsuitable ( MBadPixelsPix::kUnsuitableRun );
416 continue;
417 }
418
419 areavars [aidx] += nphotpera*nphotpera;
420 areaphotons [aidx] += nphotpera;
421 numareavalid [aidx] ++;
422
423 sectorvars [sector] += nphotpera*nphotpera;
424 sectorphotons [sector] += nphotpera;
425 numsectorvalid[sector] ++;
426 }
427
428 *fLog << endl;
429
430 for (UInt_t aidx=0; aidx<nareas; aidx++)
431 {
432
433 if (numareavalid[aidx] == 1)
434 areavars[aidx] = 0.;
435 else if (numareavalid[aidx] == 0)
436 {
437 areaphotons[aidx] = -1.;
438 areavars[aidx] = -1.;
439 }
440 else
441 {
442 areavars[aidx] = (areavars[aidx] - areaphotons[aidx]*areaphotons[aidx]/numareavalid[aidx])
443 / (numareavalid[aidx]-1.);
444 areaphotons[aidx] /= numareavalid[aidx];
445 }
446
447
448 MCalibrationTestPix &avpix = (MCalibrationTestPix&)testcam->GetAverageArea(aidx);
449
450 if (areavars[aidx] < 0. || areaphotons[aidx] <= 0.)
451 {
452 *fLog << warn << GetDescriptor()
453 << ": Mean number of photons per area in area index "
454 << aidx << " could not be calculated! Mean: " << areaphotons[aidx]
455 << " Variance: " << areavars[aidx] << endl;
456 avpix.SetExcluded();
457 continue;
458 }
459
460 const Float_t areaerr = TMath::Sqrt(areavars[aidx]);
461
462 avpix.SetNumPhotonsPerArea (areaphotons[aidx]);
463 avpix.SetNumPhotonsPerAreaErr(areaerr );
464
465 *fLog << inf << GetDescriptor() << ": Mean number of equiv. Cher. photons "
466 << "per area in area idx " << aidx << ": "
467 << Form("%5.3f+-%5.4f [ph/mm^2]",areaphotons[aidx],areaerr) << endl;
468 }
469
470 *fLog << endl;
471
472 for (UInt_t sector=0; sector<nsectors; sector++)
473 {
474
475 if (numsectorvalid[sector] == 1)
476 sectorvars[sector] = 0.;
477 else if (numsectorvalid[sector] == 0)
478 {
479 sectorphotons[sector] = -1.;
480 sectorvars[sector] = -1.;
481 }
482 else
483 {
484 sectorvars[sector] = (sectorvars[sector]
485 - sectorphotons[sector]*sectorphotons[sector]/numsectorvalid[sector]
486 )
487 / (numsectorvalid[sector]-1.);
488 sectorphotons[sector] /= numsectorvalid[sector];
489 }
490
491 MCalibrationTestPix &avpix = (MCalibrationTestPix&)testcam->GetAverageSector(sector);
492
493 if (sectorvars[sector] < 0. || sectorphotons[sector] <= 0.)
494 {
495 *fLog << warn << GetDescriptor()
496 << ": Mean number of calibrated photons per area in sector "
497 << sector << " could not be calculated! Mean: " << sectorphotons[sector]
498 << " Variance: " << sectorvars[sector] << endl;
499 avpix.SetExcluded();
500 continue;
501 }
502
503
504 const Float_t sectorerr = TMath::Sqrt(sectorvars[sector]);
505
506 avpix.SetNumPhotonsPerArea (sectorphotons[sector]);
507 avpix.SetNumPhotonsPerAreaErr(sectorerr );
508
509 *fLog << inf << GetDescriptor() << ": Mean number of equiv. Cher. photons "
510 << "per area in sector " << sector << ": "
511 << Form("%5.3f+-%5.4f [ph/mm^2]",sectorphotons[sector],sectorerr) << endl;
512 }
513
514 return;
515}
516
517
518// -----------------------------------------------------------------------------------------------
519//
520// Print out statistics about not interpolated pixels
521//
522void MCalibrationTestCalc::FinalizeNotInterpolated()
523{
524
525 MCalibrationTestCam *testcam = fIntensCam
526 ? (MCalibrationTestCam*)fIntensCam->GetCam() : fCam;
527
528 const Int_t areas = fGeom->GetNumAreas();
529 TArrayI *newarr[areas];
530
531 for (Int_t aidx=0; aidx<areas; aidx++)
532 newarr[aidx] = new TArrayI(0);
533
534 for (Int_t i=0; i<testcam->GetSize(); i++)
535 {
536 const Int_t aidx = (*fGeom)[i].GetAidx();
537 MCalibrationTestPix &pix = (MCalibrationTestPix&)(*testcam)[i];
538 if (pix.IsExcluded())
539 {
540 const Int_t size = newarr[aidx]->GetSize();
541 newarr[aidx]->Set(size+1);
542 newarr[aidx]->AddAt(i,size);
543 }
544 }
545
546 Int_t num = 0;
547
548 for (Int_t aidx = 0; aidx<areas; aidx++)
549 {
550 *fLog << endl;
551 *fLog << " " << setw(7)
552 << "Not interpolated pixels by in area idx " << aidx << ": ";
553 for (Int_t i=0; i<newarr[aidx]->GetSize(); i++)
554 {
555 *fLog << newarr[aidx]->At(i) << " ";
556 num++;
557 }
558 testcam->SetNumUninterpolated(newarr[aidx]->GetSize(),aidx);
559 *fLog << endl;
560 }
561
562 *fLog << " " << setw(7) << num << " total not interpolateable pixels " << endl;
563
564}
565
566Int_t MCalibrationTestCalc::CalcMaxNumBadPixelsCluster()
567{
568
569 MCalibrationTestCam *testcam = fIntensCam
570 ? (MCalibrationTestCam*)fIntensCam->GetCam() : fCam;
571
572 TArrayI arr(0);
573
574 for (Int_t i=0; i<testcam->GetSize(); i++)
575 {
576 MCalibrationTestPix &pix = (MCalibrationTestPix&)(*testcam)[i];
577 if (pix.IsExcluded())
578 {
579 const Int_t s = arr.GetSize();
580 arr.Set(s+1);
581 arr[s] = i;
582 }
583 }
584
585 const Int_t size = arr.GetSize();
586
587 if (size == 0)
588 return 0;
589
590 if (size == 1)
591 return 1;
592
593 TArrayI knownpixels(0);
594 Int_t clustersize = 1;
595 Int_t oldclustersize = 0;
596 //
597 // Loop over the not-interpolateable pixels:
598 //
599 for (Int_t i=0; i<size; i++)
600 {
601
602 const Int_t id = arr[i];
603 const Int_t knownsize = knownpixels.GetSize();
604 knownpixels.Set(knownsize+1);
605 knownpixels[knownsize] = id;
606 LoopNeighbours(arr, knownpixels, clustersize, id);
607 if (clustersize > oldclustersize)
608 oldclustersize = clustersize;
609 clustersize = 1;
610 }
611
612 return oldclustersize;
613
614}
615
616
617void MCalibrationTestCalc::LoopNeighbours( const TArrayI &arr, TArrayI &knownpixels, Int_t &clustersize, const Int_t idx )
618{
619
620 const MGeomPix &pix = (*fGeom)[idx];
621 const Byte_t neighbours = pix.GetNumNeighbors();
622
623 //
624 // Loop over the next neighbours:
625 // - Check if they are already in the list of known pixels
626 // - If not, call loopneighbours for the rest
627 // - grow clustersize for those
628 //
629 for (Int_t i=0;i<neighbours;i++)
630 {
631 const Int_t newid = pix.GetNeighbor(i);
632 Bool_t known = kFALSE;
633
634 for (Int_t j=knownpixels.GetSize()-1;j>=0;j--)
635 if (newid == knownpixels.At(j))
636 {
637 known = kTRUE;
638 break;
639 }
640 if (known)
641 continue;
642
643 for (Int_t k=0;k<arr.GetSize();k++)
644 if (newid == arr.At(k))
645 {
646 // This is an unknown, new pixel in the cluster!!
647 clustersize++;
648 const Int_t knownsize = knownpixels.GetSize();
649 knownpixels.Set(knownsize+1);
650 knownpixels[knownsize] = newid;
651 LoopNeighbours(arr, knownpixels, clustersize, newid);
652 }
653 }
654}
655
656// --------------------------------------------------------------------------
657//
658// Set the path for output file
659//
660void MCalibrationTestCalc::SetOutputPath(TString path)
661{
662 fOutputPath = path;
663 if (fOutputPath.EndsWith("/"))
664 fOutputPath = fOutputPath(0, fOutputPath.Length()-1);
665}
666
667void MCalibrationTestCalc::SetOutputFile(TString file)
668{
669 fOutputFile = file;
670}
671
672// --------------------------------------------------------------------------
673//
674// Get the output file
675//
676const char* MCalibrationTestCalc::GetOutputFile()
677{
678 return Form("%s/%s", (const char*)fOutputPath, (const char*)fOutputFile);
679}
680
Note: See TracBrowser for help on using the repository browser.