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

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