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

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