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

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