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

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