source: trunk/MagicSoft/Mars/mhcalib/MHCalibrationHiLoCam.cc@ 5762

Last change on this file since 5762 was 5749, checked in by gaug, 20 years ago
*** empty log message ***
File size: 16.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 02/2004 <mailto:markus@ifae.es>
19!
20! Copyright: MAGIC Software Development, 2000-2004
21!
22!
23\* ======================================================================== */
24/////////////////////////////////////////////////////////////////////////////
25//
26// MHCalibrationHiLoCam
27//
28// Fills the extracted high-gain low-gain charge ratios of MArrivalTimeCam into
29// the MHCalibrationPix-classes MHCalibrationPix for every:
30//
31// - Pixel, stored in the TObjArray's MHCalibrationCam::fHiGainArray
32// or MHCalibrationCam::fHiGainArray, respectively, depending if
33// MArrivalTimePix::IsLoGainUsed() is set.
34//
35// - Average pixel per AREA index (e.g. inner and outer for the MAGIC camera),
36// stored in the TObjArray's MHCalibrationCam::fAverageHiGainAreas and
37// MHCalibrationCam::fAverageHiGainAreas
38//
39// - Average pixel per camera SECTOR (e.g. sectors 1-6 for the MAGIC camera),
40// stored in the TObjArray's MHCalibrationCam::fAverageHiGainSectors
41// and MHCalibrationCam::fAverageHiGainSectors
42//
43// The histograms are fitted to a Gaussian, mean and sigma with its errors
44// and the fit probability are extracted. If none of these values are NaN's and
45// if the probability is bigger than MHGausEvents::fProbLimit (default: 0.5%),
46// the fit is declared valid.
47// Otherwise, the fit is repeated within ranges of the previous mean
48// +- MHCalibrationPix::fPickupLimit (default: 5) sigma (see MHCalibrationPix::RepeatFit())
49// In case this does not make the fit valid, the histogram means and RMS's are
50// taken directly (see MHCalibrationPix::BypassFit()) and the following flags are set:
51// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiLoNotFitted ) and
52// - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun )
53//
54// Outliers of more than MHCalibrationPix::fPickupLimit (default: 5) sigmas
55// from the mean are counted as Pickup events (stored in MHCalibrationPix::fPickup)
56//
57// The class also fills arrays with the signal vs. event number, creates a fourier
58// spectrum (see MHGausEvents::CreateFourierSpectrum()) and investigates if the
59// projected fourier components follow an exponential distribution.
60// In case that the probability of the exponential fit is less than
61// MHGausEvents::fProbLimit (default: 0.5%), the following flags are set:
62// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiLoOscillating ) and
63// - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun )
64//
65// This same procedure is performed for the average pixels.
66//
67// The following results are written into MCalibrationHiLoCam:
68//
69// - MCalibrationPix::SetMean()
70// - MCalibrationPix::SetMeanErr()
71// - MCalibrationPix::SetSigma()
72// - MCalibrationPix::SetSigmaErr()
73// - MCalibrationPix::SetProb()
74// - MCalibrationPix::SetNumPickup()
75//
76// For all averaged areas, the fitted sigma is multiplied with the square root of
77// the number involved pixels in order to be able to compare it to the average of
78// sigmas in the camera.
79//
80/////////////////////////////////////////////////////////////////////////////
81#include "MHCalibrationHiLoCam.h"
82#include "MHCalibrationPix.h"
83
84#include "MLog.h"
85#include "MLogManip.h"
86
87#include "MParList.h"
88
89#include "MCalibrationIntensityHiLoCam.h"
90
91#include "MCalibrationHiLoCam.h"
92#include "MCalibrationCam.h"
93#include "MCalibrationPix.h"
94
95#include "MExtractedSignalCam.h"
96#include "MExtractedSignalPix.h"
97
98#include "MRawEvtPixelIter.h"
99
100#include "MGeomCam.h"
101#include "MGeomPix.h"
102
103#include "MBadPixelsIntensityCam.h"
104#include "MBadPixelsCam.h"
105#include "MBadPixelsPix.h"
106
107#include <TOrdCollection.h>
108#include <TPad.h>
109#include <TVirtualPad.h>
110#include <TCanvas.h>
111#include <TStyle.h>
112#include <TF1.h>
113#include <TLine.h>
114#include <TLatex.h>
115#include <TLegend.h>
116#include <TGraph.h>
117
118ClassImp(MHCalibrationHiLoCam);
119
120using namespace std;
121
122const Int_t MHCalibrationHiLoCam::fgNbins = 200;
123const Axis_t MHCalibrationHiLoCam::fgFirst = 0.;
124const Axis_t MHCalibrationHiLoCam::fgLast = 20.;
125const Float_t MHCalibrationHiLoCam::fgProbLimit = 0.;
126const TString MHCalibrationHiLoCam::gsHistName = "HiLo";
127const TString MHCalibrationHiLoCam::gsHistTitle = "HiGain vs. LoGain";
128const TString MHCalibrationHiLoCam::gsHistXTitle = "Amplification Ratio [1]";
129const TString MHCalibrationHiLoCam::gsHistYTitle = "Nr. events";
130const Byte_t MHCalibrationHiLoCam::fgLowerLim = 200;
131const Byte_t MHCalibrationHiLoCam::fgUpperLim = 252;
132// --------------------------------------------------------------------------
133//
134// Default Constructor.
135//
136// Sets:
137// - fNbins to fgNbins
138// - fFirst to fgFirst
139// - fLast to fgLast
140//
141// - fHistName to gsHistName
142// - fHistTitle to gsHistTitle
143// - fHistXTitle to gsHistXTitle
144// - fHistYTitle to gsHistYTitle
145//
146// - fLowerLimt to fgLowerLim
147// - fUpperLimt to fgUpperLim
148//
149MHCalibrationHiLoCam::MHCalibrationHiLoCam(const char *name, const char *title)
150{
151
152 fName = name ? name : "MHCalibrationHiLoCam";
153 fTitle = title ? title : "Histogram class for the high-gain vs. low-gain amplification ratio calibration";
154
155 SetNbins(fgNbins);
156 SetFirst(fgFirst);
157 SetLast (fgLast );
158
159 SetProbLimit(fgProbLimit);
160
161 SetHistName (gsHistName .Data());
162 SetHistTitle (gsHistTitle .Data());
163 SetHistXTitle(gsHistXTitle.Data());
164 SetHistYTitle(gsHistYTitle.Data());
165
166 SetLowerLim();
167 SetUpperLim();
168}
169
170// --------------------------------------------------------------------------
171//
172// Creates new MHCalibrationHiLoCam only with the averaged areas:
173// the rest has to be retrieved directly, e.g. via:
174// MHCalibrationHiLoCam *cam = MParList::FindObject("MHCalibrationHiLoCam");
175// - cam->GetAverageSector(5).DrawClone();
176// - (*cam)[100].DrawClone()
177//
178TObject *MHCalibrationHiLoCam::Clone(const char *) const
179{
180
181 MHCalibrationHiLoCam *cam = new MHCalibrationHiLoCam();
182
183 //
184 // Copy the data members
185 //
186 cam->fColor = fColor;
187 cam->fRunNumbers = fRunNumbers;
188 cam->fPulserFrequency = fPulserFrequency;
189 cam->fFlags = fFlags;
190 cam->fNbins = fNbins;
191 cam->fFirst = fFirst;
192 cam->fLast = fLast;
193
194 //
195 // Copy the MArrays
196 //
197 cam->fAverageAreaRelSigma = fAverageAreaRelSigma;
198 cam->fAverageAreaRelSigmaVar = fAverageAreaRelSigmaVar;
199 cam->fAverageAreaSat = fAverageAreaSat;
200 cam->fAverageAreaSigma = fAverageAreaSigma;
201 cam->fAverageAreaSigmaVar = fAverageAreaSigmaVar;
202 cam->fAverageAreaNum = fAverageAreaNum;
203 cam->fAverageSectorNum = fAverageSectorNum;
204
205 if (!IsAverageing())
206 return cam;
207
208 const Int_t navhi = fAverageHiGainAreas->GetSize();
209
210 for (int i=0; i<navhi; i++)
211 cam->fAverageHiGainAreas->AddAt(GetAverageHiGainArea(i).Clone(),i);
212
213 return cam;
214}
215
216// --------------------------------------------------------------------------
217//
218// Gets the pointers to:
219// - MRawEvtData
220//
221Bool_t MHCalibrationHiLoCam::SetupHists(const MParList *pList)
222{
223
224 fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
225 if (!fRawEvt)
226 {
227 *fLog << err << dbginf << "MRawEvtData not found... aborting." << endl;
228 return kFALSE;
229 }
230
231 return kTRUE;
232}
233
234// --------------------------------------------------------------------------
235//
236// Gets or creates the pointers to:
237// - MCalibrationHiLoCam
238//
239// Searches pointer to:
240// - MExtractedSignalCam
241//
242// Calls:
243// - MHCalibrationCam::InitHiGainArrays()
244//
245// Sets:
246// - SetLoGain(kFALSE);
247// - fSumarea to nareas
248// - fSumsector to nareas
249// - fNumarea to nareas
250// - fNumsector to nareas
251//
252Bool_t MHCalibrationHiLoCam::ReInitHists(MParList *pList)
253{
254
255 fIntensCam = (MCalibrationIntensityCam*)pList->FindObject(AddSerialNumber("MCalibrationIntensityHiLoCam"));
256 if (fIntensCam)
257 *fLog << inf << "Found MCalibrationIntensityHiLoCam ... " << endl;
258 else
259 {
260 fCam = (MCalibrationCam*)pList->FindObject(AddSerialNumber("MCalibrationHiLoCam"));
261 if (!fCam)
262 {
263 fCam = (MCalibrationCam*)pList->FindCreateObj(AddSerialNumber("MCalibrationHiLoCam"));
264 if (!fCam)
265 {
266 *fLog << err << "Cannot find nor create MCalibrationHiLoCam ... abort." << endl;
267 return kFALSE;
268 }
269 fCam->Init(*fGeom);
270 }
271 }
272
273 MExtractedSignalCam *signal = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam");
274 if (!signal)
275 {
276 *fLog << err << "MExtractedSignalCam not found... abort." << endl;
277 return kFALSE;
278 }
279
280 const Int_t npixels = fGeom->GetNumPixels();
281 const Int_t nsectors = fGeom->GetNumSectors();
282 const Int_t nareas = fGeom->GetNumAreas();
283
284 InitHiGainArrays(npixels,nareas,nsectors);
285
286 fSumarea .Set(nareas);
287 fSumsector.Set(nsectors);
288 fNumarea .Set(nareas);
289 fNumsector.Set(nsectors);
290
291 SetLoGain(kFALSE);
292
293 return kTRUE;
294}
295
296// -------------------------------------------------------------------------------
297//
298// Retrieves pointer to MExtractedSignalCam:
299//
300// Retrieves from MGeomCam:
301// - number of pixels
302// - number of pixel areas
303// - number of sectors
304//
305// Fills histograms (MHGausEvents::FillHistAndArray()) with:
306// - MExtractedSignalPix::GetExtractedSignalHiGain(pixid) / MExtractedSignalPix::GetExtractedSignalLoGain;
307// if the high-gain signal lies in between the limits: fLowerLim and fUpperLim
308//
309Bool_t MHCalibrationHiLoCam::FillHists(const MParContainer *par, const Stat_t w)
310{
311
312 MExtractedSignalCam *signal = (MExtractedSignalCam*)par;
313 if (!signal)
314 {
315 gLog << err << "No argument in MExtractedSignal::Fill... abort." << endl;
316 return kFALSE;
317 }
318
319 const Int_t nareas = fGeom->GetNumAreas();
320 const Int_t nsectors = fGeom->GetNumSectors();
321
322 fSumarea .Reset();
323 fSumsector.Reset();
324 fNumarea .Reset();
325 fNumsector.Reset();
326
327 MRawEvtPixelIter pixel(fRawEvt);
328
329 while (pixel.Next())
330 {
331
332 const Byte_t max = pixel.GetMaxHiGainSample();
333
334 if (max < fLowerLim || max > fUpperLim)
335 continue;
336
337 const UInt_t pixidx = pixel.GetPixelId();
338
339 MHCalibrationPix &hist = (*this)[pixidx];
340
341 if (hist.IsExcluded())
342 continue;
343
344 const MExtractedSignalPix &pix = (*signal)[pixidx];
345 const Int_t aidx = (*fGeom)[pixidx].GetAidx();
346 const Int_t sector = (*fGeom)[pixidx].GetSector();
347
348 const Float_t ratio = pix.GetExtractedSignalHiGain() / pix.GetExtractedSignalLoGain();
349
350 hist.FillHistAndArray(ratio) ;
351 fSumarea [aidx] += ratio;
352 fNumarea [aidx] ++;
353 fSumsector[sector] += ratio;
354 fNumsector[sector] ++;
355 }
356
357 for (Int_t j=0; j<nareas; j++)
358 {
359 MHCalibrationPix &hist = GetAverageHiGainArea(j);
360 hist.FillHistAndArray(fNumarea[j] == 0 ? 0. : fSumarea[j]/fNumarea[j]);
361 }
362
363 for (Int_t j=0; j<nsectors; j++)
364 {
365 MHCalibrationPix &hist = GetAverageHiGainSector(j);
366 hist.FillHistAndArray(fNumsector[j] == 0 ? 0. : fSumsector[j]/fNumsector[j]);
367 }
368
369 return kTRUE;
370}
371
372// --------------------------------------------------------------------------
373//
374// Calls:
375// - MHCalibrationCam::FitHiGainArrays() with flags:
376// MBadPixelsPix::kHiLoNotFitted and MBadPixelsPix::kHiLoOscillating
377// - MHCalibrationCam::FitLoGainArrays() with flags:
378// MBadPixelsPix::kHiLoNotFitted and MBadPixelsPix::kHiLoOscillating
379//
380Bool_t MHCalibrationHiLoCam::FinalizeHists()
381{
382
383 *fLog << endl;
384
385 MCalibrationCam *hilocam = fIntensCam ? fIntensCam->GetCam() : fCam;
386 MBadPixelsCam *badcam = fIntensBad ? fIntensBad->GetCam() : fBadPixels;
387
388 const Int_t nareas = fAverageHiGainAreas->GetSize();
389 const Int_t nsectors = fAverageHiGainSectors->GetSize();
390
391 TArrayI satarea(nareas);
392 TArrayI satsect(nsectors);
393 fNumarea .Reset();
394 fNumsector.Reset();
395
396 for (Int_t i=0; i<fHiGainArray->GetSize(); i++)
397 {
398
399 MHCalibrationPix &hist = (*this)[i];
400
401 if (hist.IsExcluded())
402 continue;
403
404 const Int_t aidx = (*fGeom)[i].GetAidx();
405 const Int_t sector = (*fGeom)[i].GetSector();
406
407 fNumarea[aidx]++;
408 fNumsector[sector]++;
409 //
410 // Check histogram overflow
411 //
412 CheckOverflow(hist);
413 }
414
415 for (Int_t j=0; j<nareas; j++)
416 {
417
418 MHCalibrationPix &hist = GetAverageHiGainArea(j);
419 //
420 // Check histogram overflow
421 //
422 CheckOverflow(hist);
423 }
424
425 for (Int_t j=0; j<fAverageHiGainSectors->GetSize(); j++)
426 {
427
428 MHCalibrationPix &hist = GetAverageHiGainSector(j);
429 //
430 // Check histogram overflow
431 //
432 CheckOverflow(hist);
433 }
434
435 FitHiGainArrays(*hilocam,*badcam,
436 MBadPixelsPix::kHiLoNotFitted,
437 MBadPixelsPix::kHiLoOscillating);
438
439 return kTRUE;
440}
441
442// --------------------------------------------------------------------------
443//
444// Sets all pixels to MBadPixelsPix::kUnreliableRun, if following flags are set:
445// - MBadPixelsPix::kHiLoNotFitted
446// - MBadPixelsPix::kHiLoOscillating
447//
448void MHCalibrationHiLoCam::FinalizeBadPixels()
449{
450
451 MBadPixelsCam *badcam = fIntensBad ? fIntensBad->GetCam() : fBadPixels;
452
453 for (Int_t i=0; i<badcam->GetSize(); i++)
454 {
455 MBadPixelsPix &bad = (*badcam)[i];
456
457 if (bad.IsUncalibrated( MBadPixelsPix::kHiLoNotFitted ))
458 bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun );
459
460 if (bad.IsUncalibrated( MBadPixelsPix::kHiLoOscillating))
461 bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun );
462
463 }
464}
465
466// --------------------------------------------------------------------------
467//
468// The types are as follows:
469//
470// Fitted values:
471// ==============
472//
473// 0: Fitted Mean High-Gain Low-Gain Charge Ratio in FADC slices (MHGausEvents::GetMean()
474// 1: Error Mean High-Gain Low-Gain Charge Ratio in FADC slices (MHGausEvents::GetMeanErr()
475// 2: Sigma fitted High-Gain Low-Gain Charge Ratio in FADC slices (MHGausEvents::GetSigma()
476// 3: Error Sigma High-Gain Low-Gain Charge Ratio in FADC slices (MHGausEvents::GetSigmaErr()
477//
478// Useful variables derived from the fit results:
479// =============================================
480//
481// 4: Returned probability of Gauss fit (calls: MHGausEvents::GetProb())
482//
483// Localized defects:
484// ==================
485//
486// 5: Gaus fit not OK (calls: MHGausEvents::IsGausFitOK())
487// 6: Fourier spectrum not OK (calls: MHGausEvents::IsFourierSpectrumOK())
488//
489Bool_t MHCalibrationHiLoCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const
490{
491
492 if (fHiGainArray->GetSize() <= idx)
493 return kFALSE;
494
495 const MHCalibrationPix &pix = (*this)[idx];
496
497 switch (type)
498 {
499 case 0:
500 val = pix.GetMean();
501 break;
502 case 1:
503 val = pix.GetMeanErr();
504 break;
505 case 2:
506 val = pix.GetSigma();
507 break;
508 case 3:
509 val = pix.GetSigmaErr();
510 break;
511 case 4:
512 val = pix.GetProb();
513 break;
514 case 5:
515 if (!pix.IsGausFitOK())
516 val = 1.;
517 break;
518 case 6:
519 if (!pix.IsFourierSpectrumOK())
520 val = 1.;
521 break;
522 default:
523 return kFALSE;
524 }
525 return kTRUE;
526}
527
528// --------------------------------------------------------------------------
529//
530// Calls MHCalibrationPix::DrawClone() for pixel idx
531//
532void MHCalibrationHiLoCam::DrawPixelContent(Int_t idx) const
533{
534 (*this)[idx].DrawClone();
535}
536
537void MHCalibrationHiLoCam::CheckOverflow( MHCalibrationPix &pix )
538{
539
540 if (pix.IsExcluded())
541 return;
542
543 TH1F *hist = pix.GetHGausHist();
544
545 Stat_t overflow = hist->GetBinContent(hist->GetNbinsX()+1);
546 if (overflow > 0.0005*hist->GetEntries())
547 {
548 *fLog << warn << "Hist-overflow " << overflow
549 << " times in " << pix.GetName() << " (w/o saturation!) " << endl;
550 }
551
552 overflow = hist->GetBinContent(0);
553 if (overflow > 0.0005*hist->GetEntries())
554 {
555 *fLog << warn << "Hist-underflow " << overflow
556 << " times in " << pix.GetName() << " (w/o saturation!) " << endl;
557 }
558}
559
Note: See TracBrowser for help on using the repository browser.