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

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