source: trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargeBlindCam.cc@ 8688

Last change on this file since 8688 was 8434, checked in by tbretz, 19 years ago
*** empty log message ***
File size: 15.1 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// MHCalibrationChargeBlindCam
27//
28// Histogram class for blind pixels in the camera. Incorporates the TObjArray's:
29// - fBlindPixelsArray (for calibrated High Gains per pixel)
30//
31/////////////////////////////////////////////////////////////////////////////
32#include "MHCalibrationChargeBlindCam.h"
33
34#include <TVirtualPad.h>
35#include <TCanvas.h>
36#include <TPad.h>
37#include <TOrdCollection.h>
38
39#include "MLog.h"
40#include "MLogManip.h"
41
42#include "MRawEvtData.h"
43#include "MRawEvtPixelIter.h"
44
45#include "MExtractedSignalBlindPixel.h"
46
47#include "MCalibrationBlindPix.h"
48
49#include "MParList.h"
50
51#include "MRawRunHeader.h"
52
53ClassImp(MHCalibrationChargeBlindCam);
54
55using namespace std;
56const Int_t MHCalibrationChargeBlindCam::fgNbins = 128;
57const Axis_t MHCalibrationChargeBlindCam::fgFirst = -0.5;
58const Axis_t MHCalibrationChargeBlindCam::fgLast = 511.5;
59const Axis_t MHCalibrationChargeBlindCam::fgSPheCut = 20.;
60const TString MHCalibrationChargeBlindCam::gsHistName = "ChargeBlind";
61const TString MHCalibrationChargeBlindCam::gsHistTitle = "Signals Blind";
62const TString MHCalibrationChargeBlindCam::gsHistXTitle = "Signal [FADC counts]";
63const TString MHCalibrationChargeBlindCam::gsHistYTitle = "Nr. events";
64// --------------------------------------------------------------------------
65//
66// Default Constructor.
67//
68// Sets:
69// - all pointers to NULL
70//
71// - fFitFunc to kEPoisson4
72// - fNbins to fgNbins
73// - fFirst to fgFirst
74// - fLast to fgLast
75// - fSPheCut to fgSPheCut
76//
77// - fHistName to gsHistName
78// - fHistTitle to gsHistTitle
79// - fHistXTitle to gsHistXTitle
80// - fHistYTitle to gsHistYTitle
81//
82// - SetAverageing (kFALSE);
83// - SetLoGain (kFALSE);
84// - SetOscillations(kFALSE);
85// - SetSizeCheck (kFALSE);
86//
87MHCalibrationChargeBlindCam::MHCalibrationChargeBlindCam(const char *name, const char *title)
88 : fRawEvt(NULL)
89{
90
91 fName = name ? name : "MHCalibrationChargeBlindCam";
92 fTitle = title ? title : "Class to fille the blind pixel histograms";
93
94 SetSPheCut();
95
96 SetBinning(fgNbins, fgFirst, fgLast);
97
98 SetHistName (gsHistName .Data());
99 SetHistTitle (gsHistTitle .Data());
100 SetHistXTitle(gsHistXTitle.Data());
101 SetHistYTitle(gsHistYTitle.Data());
102
103 SetAverageing (kFALSE);
104 SetLoGain (kFALSE);
105 SetOscillations(kFALSE);
106 SetSizeCheck (kFALSE);
107
108 SetFitFunc(MHCalibrationChargeBlindPix::kEPoisson4);
109}
110
111// --------------------------------------------------------------------------
112//
113// Gets the pointers to:
114// - MRawEvtData
115//
116Bool_t MHCalibrationChargeBlindCam::SetupHists(const MParList *pList)
117{
118
119 fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
120 if (!fRawEvt)
121 {
122 *fLog << err << dbginf << "MRawEvtData not found... aborting." << endl;
123 return kFALSE;
124 }
125
126 return kTRUE;
127}
128
129// --------------------------------------------------------------------------
130//
131// Gets or creates the pointers to:
132// - MExtractedSignalBlindPixel
133// - MCalibrationChargeCam or MCalibrationBlindCam
134//
135// Initializes the number of used FADC slices from MExtractedSignalCam
136// into MCalibrationChargeCam and test for changes in that variable
137//
138// Retrieve:
139// - fRunHeader->GetNumSamplesHiGain();
140//
141// Initializes the High Gain Arrays:
142//
143// - Expand fHiGainArrays to nblindpixels
144//
145// - For every entry in the expanded arrays:
146// * Initialize an MHCalibrationPix
147// * Set Binning from fNbins, fFirst and fLast
148// * Set Histgram names and titles from fHistName and fHistTitle
149// * Set X-axis and Y-axis titles from fHistXTitle and fHistYTitle
150// * Call InitHists
151//
152Bool_t MHCalibrationChargeBlindCam::ReInitHists(MParList *pList)
153{
154
155 if (!InitCams(pList,"Blind"))
156 return kFALSE;
157
158 MExtractedSignalBlindPixel *signal =
159 (MExtractedSignalBlindPixel*)pList->FindObject(AddSerialNumber("MExtractedSignalBlindPixel"));
160 if (!signal)
161 {
162 *fLog << err << "MExtractedSignalBlindPixel not found... abort." << endl;
163 return kFALSE;
164 }
165
166 const Int_t nblindpixels = signal->GetNumBlindPixels();
167 const Int_t samples = signal->GetNumFADCSamples();
168 const Int_t integ = signal->IsExtractionType( MExtractBlindPixel::kIntegral );
169
170 if (fHiGainArray->GetSize()==0)
171 {
172 for (Int_t i=0; i<nblindpixels; i++)
173 {
174 fHiGainArray->AddAt(new MHCalibrationChargeBlindPix(Form("%s%s%d",fHistName.Data(),"Pix",i),
175 Form("%s%s%d",fHistTitle.Data()," Pixel ",i)),i);
176
177 MHCalibrationChargeBlindPix &pix = (MHCalibrationChargeBlindPix&)(*this)[i];
178
179 const Axis_t last = integ ? ((fLast+0.5)*samples)-0.5 : fLast;
180
181 pix.SetBinning(fNbins, fFirst, last);
182 pix.SetSinglePheCut ( integ ? fSPheCut * samples : fSPheCut );
183 pix.SetFitFunc ( integ ? MHCalibrationChargeBlindPix::kEPoisson5 : fFitFunc );
184
185 pix.InitBins();
186
187 }
188 }
189 return kTRUE;
190}
191
192// --------------------------------------------------------------------------
193//
194// Resets the histogram titles for each entry in:
195// - fHiGainArray
196//
197void MHCalibrationChargeBlindCam::ResetHistTitles()
198{
199
200 TH1F *h;
201
202 if (fHiGainArray)
203 for (Int_t i=0;i<fHiGainArray->GetSize(); i++)
204 {
205 MHCalibrationPix &pix = (*this)[i];
206 h = pix.GetHGausHist();
207 h->SetName (Form("%s%s%s%d","H",fHistName.Data(),"Pix",i));
208 h->SetTitle(Form("%s%s%d%s",fHistTitle.Data()," Pixel ",i," Runs: "));
209 h->SetXTitle(fHistXTitle.Data());
210 h->SetYTitle(fHistYTitle.Data());
211 }
212}
213
214// --------------------------------------------------------------------------
215//
216// Retrieves from MExtractedSignalBlindPixel:
217// - number of blind pixels
218//
219// Retrieves from MExtractedSignalBlindPixel:
220// - number of FADC samples
221// - extracted signal
222// - blind Pixel ID
223//
224// Resizes (if necessary):
225// - fASinglePheFADCSlices to sum of HiGain and LoGain samples
226// - fAPedestalFADCSlices to sum of HiGain and LoGain samples
227//
228// Fills the following histograms:
229// - MHGausEvents::FillHistAndArray(signal)
230//
231// Creates MRawEvtPixelIter, jumps to blind pixel ID,
232// fills the vectors fASinglePheFADCSlices and fAPedestalFADCSlices
233// with the full FADC slices, depending on the size of the signal w.r.t. fSinglePheCut
234//
235Bool_t MHCalibrationChargeBlindCam::FillHists(const MParContainer *par, const Stat_t w)
236{
237
238 MExtractedSignalBlindPixel *signal = (MExtractedSignalBlindPixel*)par;
239 if (!signal)
240 {
241 *fLog << err << "No argument in MExtractedSignalBlindCam::Fill... abort." << endl;
242 return kFALSE;
243 }
244
245 const Int_t nblindpixels = signal->GetNumBlindPixels();
246
247 if (GetSize() != nblindpixels)
248 {
249 gLog << err << "ERROR - Size mismatch... abort." << endl;
250 return kFALSE;
251 }
252
253 Float_t slices = (Float_t)signal->GetNumFADCSamples();
254
255 if (slices == 0.)
256 {
257 *fLog << err << dbginf
258 << "Number of used signal slices in MExtractedSignalBlindPix "
259 << "is zero ... abort."
260 << endl;
261 return kFALSE;
262 }
263
264 for (Int_t i=0; i<nblindpixels; i++)
265 {
266
267 //
268 // Signal extraction and histogram filling
269 // If filter has been applied, sig has been set to -1.
270 //
271 const Float_t sig = signal->GetExtractedSignal(i);
272
273 if (sig < -0.5)
274 continue;
275
276 MHCalibrationChargeBlindPix &hist = (MHCalibrationChargeBlindPix&)(*this)[i];
277
278 hist.FillHist(sig);
279 //
280 // In order to study the single-phe posistion, we extract the slices
281 //
282 const Int_t blindpixIdx = signal->GetBlindPixelIdx(i);
283
284 MRawEvtPixelIter pixel(fRawEvt);
285 pixel.Jump(blindpixIdx);
286/*
287 if (sig > fSPheCut)
288 hist.FillSinglePheFADCSlices(pixel);
289 else
290 hist.FillPedestalFADCSlices(pixel);
291 */
292 }
293
294 return kTRUE;
295}
296
297// --------------------------------------------------------------------------
298//
299// For all TObjArray's (including the averaged ones), the following steps are performed:
300//
301// 1) Returns if the pixel is excluded.
302// 2) Tests saturation. In case yes, set the flag: MCalibrationPix::SetHiGainSaturation()
303// or the flag: MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainSaturated )
304// 3) Store the absolute arrival times in the MCalibrationChargePix's. If flag
305// MCalibrationPix::IsHiGainSaturation() is set, the Low-Gain arrival times are stored,
306// otherwise the Hi-Gain ones.
307// 4) Calls to MHCalibrationCam::FitHiGainArrays() and MCalibrationCam::FitLoGainArrays()
308// with the flags:
309// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainNotFitted )
310// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainNotFitted )
311// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainOscillating )
312// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainOscillating )
313//
314Bool_t MHCalibrationChargeBlindCam::FinalizeHists()
315{
316
317 *fLog << endl;
318
319 for (Int_t i=0; i<fHiGainArray->GetSize(); i++)
320 {
321 MHCalibrationChargeBlindPix &hist = (MHCalibrationChargeBlindPix&)(*this)[i];
322
323 TH1F *h = hist.GetHGausHist();
324
325 switch (fColor)
326 {
327 case MCalibrationCam::kNONE:
328 break;
329 case MCalibrationCam::kBLUE:
330 h->SetTitle( Form("%s%s", h->GetTitle(),"BLUE "));
331 break;
332 case MCalibrationCam::kGREEN:
333 h->SetTitle( Form("%s%s", h->GetTitle(),"GREEN "));
334 break;
335 case MCalibrationCam::kUV:
336 h->SetTitle( Form("%s%s", h->GetTitle(),"UV "));
337 break;
338 case MCalibrationCam::kCT1:
339 h->SetTitle( Form("%s%s", h->GetTitle(),"CT1-Pulser "));
340 break;
341 }
342
343 Stat_t overflow = h->GetBinContent(h->GetNbinsX()+1);
344 if (overflow > 0.1)
345 {
346 *fLog << warn << GetDescriptor()
347 << ": Histogram Overflow occurred " << overflow
348 << " times in blind pixel: " << i << endl;
349 }
350
351 overflow = h->GetBinContent(0);
352 if (overflow > 0.1)
353 {
354 *fLog << warn << GetDescriptor()
355 << ": Histogram Underflow occurred " << overflow
356 << " times in blind pixel: " << i << endl;
357 }
358
359 FitBlindPixel(hist, (MCalibrationBlindPix&)(*fCam)[i]);
360 }
361
362 return kTRUE;
363}
364
365
366// --------------------------------------------------------------------------
367//
368// Returns kFALSE, if empty
369//
370// - Creates the fourier spectrum and sets bit MHGausEvents::IsFourierSpectrumOK()
371// - Retrieves the pedestals from MExtractedSignalBlindPixel
372// - Normalizes fASinglePheFADCSlices and fAPedestalFADCSlices
373// - Executes FitPedestal()
374// - Executes FitSinglePhe()
375// - Retrieves fit results and stores them in MCalibrationBlindPix
376//
377void MHCalibrationChargeBlindCam::FitBlindPixel(MHCalibrationChargeBlindPix &hist, MCalibrationBlindPix &pix)
378{
379
380 if (hist.IsEmpty())
381 {
382 *fLog << err << GetDescriptor() << " ID: " << hist.GetName()
383 << " My histogram has not been filled !! " << endl;
384 return;
385 }
386
387 //hist.FinalizeSinglePheSpectrum();
388
389 hist.FitPedestal();
390
391 pix.SetValid(kTRUE);
392
393 if (hist.FitSinglePhe())
394 pix.SetSinglePheFitOK();
395 else
396 pix.SetValid(hist.IsPedestalFitOK());
397
398 pix.SetLambda ( hist.GetLambda () );
399 pix.SetLambdaVar ( hist.GetLambdaErr()*hist.GetLambdaErr() );
400 pix.SetMu0 ( hist.GetMu0 () );
401 pix.SetMu0Err ( hist.GetMu0Err () );
402 pix.SetMu1 ( hist.GetMu1 () );
403 pix.SetMu1Err ( hist.GetMu1Err () );
404 pix.SetSigma0 ( hist.GetSigma0 () );
405 pix.SetSigma0Err ( hist.GetSigma0Err() );
406 pix.SetSigma1 ( hist.GetSigma1 () );
407 pix.SetSigma1Err ( hist.GetSigma1Err() );
408 pix.SetProb ( hist.GetProb () );
409
410 pix.SetLambdaCheck ( hist.GetLambdaCheck() );
411 pix.SetLambdaCheckErr ( hist.GetLambdaCheckErr() );
412}
413
414// --------------------------------------------------------------------------
415//
416// This Clone-function has to be different from the MHCalibrationCam
417// Clone function since it does not store and display the averaged values
418// (in fAverageAreas), but the blind pixels stored in fHiGainArray.
419//
420// Creates new MHCalibrationChargeBlindCam and
421// Clones MHCalibrationChargeBlindPix's individually
422//
423#if 0
424TObject *MHCalibrationChargeBlindCam::Clone(const char *name) const
425{
426
427 MHCalibrationChargeBlindCam *cam = new MHCalibrationChargeBlindCam();
428
429 //
430 // Copy the data members
431 //
432 cam->fRunNumbers = fRunNumbers;
433 cam->fPulserFrequency = fPulserFrequency;
434 cam->fFlags = fFlags;
435 cam->fNbins = fNbins;
436 cam->fFirst = fFirst;
437 cam->fLast = fLast;
438 cam->fFitFunc = fFitFunc;
439
440 const Int_t nhi = fHiGainArray->GetSize();
441
442 for (int i=0; i<nhi; i++)
443 cam->fHiGainArray->AddAt((*this)[i].Clone(),i);
444
445 return cam;
446}
447#endif
448// -----------------------------------------------------------------------------
449//
450// Default draw:
451//
452// Displays the averaged areas, both High Gain and Low Gain
453//
454// Calls the Draw of the fAverageHiGainAreas and fAverageLoGainAreas objects with options
455//
456void MHCalibrationChargeBlindCam::Draw(Option_t *opt)
457{
458
459 const Int_t size = fHiGainArray->GetSize();
460
461 if (size == 0)
462 return;
463
464 TString option(opt);
465 option.ToLower();
466
467 TVirtualPad *pad = gPad ? gPad : MH::MakeDefCanvas(this);
468 pad->SetBorderMode(0);
469
470 switch (size)
471 {
472 case 1:
473 break;
474 case 2:
475 pad->Divide(2,1);
476 break;
477 case 3:
478 case 4:
479 pad->Divide(2,2);
480 break;
481 default:
482 pad->Divide(size/2+1,size/2+1);
483 break;
484 }
485
486 for (Int_t i=0; i<size;i++)
487 {
488 pad->cd(i+1);
489 (*this)[i].Draw(option);
490 }
491
492 pad->Modified();
493 pad->Update();
494
495}
496
497Int_t MHCalibrationChargeBlindCam::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
498{
499
500 Bool_t rc = kFALSE;
501
502 if (MHCalibrationCam::ReadEnv(env,prefix,print))
503 rc = kTRUE;
504
505 if (IsEnvDefined(env, prefix, "SPheCut", print))
506 {
507 SetSPheCut(GetEnvValue(env, prefix, "SPheCut", fSPheCut));
508 rc = kTRUE;
509 }
510
511 // FIXME: GetEnvValue does not work with enums yet
512 /*
513 if (IsEnvDefined(env, prefix, "FitFunc", print))
514 {
515 SetFitFunc((Int_t)GetEnvValue(env, prefix, "FitFunc", fFitFunc));
516 rc = kTRUE;
517 }
518 */
519 return rc;
520}
Note: See TracBrowser for help on using the repository browser.