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

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