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

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