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

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