source: trunk/MagicSoft/Mars/mcalib/MCalibrationChargeCalc.cc@ 3281

Last change on this file since 3281 was 3264, checked in by gaug, 22 years ago
*** empty log message ***
File size: 17.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//
27// MCalibrationChargeCalc
28//
29// Task to calculate the calibration conversion factors from the FADC
30// time slices. The integrated time slices have to be delivered by an
31// MExtractedSignalCam. The pedestals by an MPedestalCam.
32//
33// The output container MCalibrationCam holds one entry of type MCalibrationChargePix
34// for every pixel. It is filled in the following way:
35//
36// ProProcess: Search for MPedestalCam, MExtractedSignalCam and MExtractedSignalBlindPixel
37// Initialize MCalibrationCam
38// Initialize pulser light wavelength
39//
40// ReInit: MCalibrationCam::InitSize(NumPixels) is called which allocates
41// memory in a TClonesArray of type MCalibrationChargePix
42// Initialize number of used FADC slices
43// Optionally exclude pixels from calibration
44//
45// Process: Every MCalibrationChargePix holds a histogram class,
46// MHCalibrationPixel which itself hold histograms of type:
47// HCharge(npix) (distribution of summed FADC time slice
48// entries)
49// HTime(npix) (distribution of position of maximum)
50// HChargevsN(npix) (distribution of charges vs. event number.
51//
52// PostProcess: All histograms HCharge(npix) are fitted to a Gaussian
53// All histograms HTime(npix) are fitted to a Gaussian
54// The histogram HBlindPixelCharge (blind pixel) is fitted to
55// a single PhE fit
56//
57// The histograms of the PIN Diode are fitted to Gaussians
58//
59// Fits can be excluded via the commands:
60// MalibrationCam::SkipBlindPixelFits() (skip all blind
61// pixel fits)
62//
63// Hi-Gain vs. Lo-Gain Calibration (very memory-intensive)
64// can be skipped with the command:
65// MalibrationCam::SkipHiLoGainCalibration()
66//
67// Input Containers:
68// MRawEvtData
69// MPedestalCam
70// MExtractedSignalCam
71// MExtractedSignalBlindPixel
72//
73// Output Containers:
74// MCalibrationCam
75//
76//////////////////////////////////////////////////////////////////////////////
77#include "MCalibrationChargeCalc.h"
78
79// FXIME: Usage of fstream is a preliminary workaround!
80#include <fstream>
81
82#include <TSystem.h>
83#include <TH1.h>
84
85#include "MLog.h"
86#include "MLogManip.h"
87
88#include "MParList.h"
89
90#include "MGeomCam.h"
91#include "MRawRunHeader.h"
92#include "MRawEvtPixelIter.h"
93
94#include "MPedestalCam.h"
95#include "MPedestalPix.h"
96
97#include "MCalibrationChargeCam.h"
98#include "MCalibrationChargePINDiode.h"
99#include "MCalibrationChargePix.h"
100
101#include "MExtractedSignalCam.h"
102#include "MExtractedSignalPix.h"
103#include "MExtractedSignalBlindPixel.h"
104
105#include "MCalibrationBlindPix.h"
106
107ClassImp(MCalibrationChargeCalc);
108
109using namespace std;
110
111const UInt_t MCalibrationChargeCalc::fgBlindPixelIdx = 559;
112const UInt_t MCalibrationChargeCalc::fgPINDiodeIdx = 9999;
113const UInt_t MCalibrationChargeCalc::fgBlindPixelSinglePheCut = 400;
114
115// --------------------------------------------------------------------------
116//
117// Default constructor.
118//
119MCalibrationChargeCalc::MCalibrationChargeCalc(const char *name, const char *title)
120 : fPedestals(NULL), fCam(NULL),
121 fRawEvt(NULL), fRunHeader(NULL), fEvtTime(NULL),
122 fSignals(NULL), fBlindPixel(NULL), fPINDiode(NULL)
123{
124
125 fName = name ? name : "MCalibrationChargeCalc";
126 fTitle = title ? title : "Task to calculate the calibration constants and MCalibrationCam ";
127
128 AddToBranchList("MRawEvtData.fHiGainPixId");
129 AddToBranchList("MRawEvtData.fLoGainPixId");
130 AddToBranchList("MRawEvtData.fHiGainFadcSamples");
131 AddToBranchList("MRawEvtData.fLoGainFadcSamples");
132
133 Clear();
134 SetBlindPixelIdx();
135 SetPINDiodeIdx();
136}
137
138void MCalibrationChargeCalc::Clear(const Option_t *o)
139{
140
141 SETBIT(fFlags, kUseBlindPixelFit);
142 SETBIT(fFlags, kUseQualityChecks);
143 SETBIT(fFlags, kHiLoGainCalibration);
144
145 fNumBlindPixelSinglePhe = 0;
146 fNumBlindPixelPedestal = 0;
147
148 fNumHiGainSamples = 0;
149 fNumLoGainSamples = 0;
150 fConversionHiLo = 0;
151 fNumExcludedPixels = 0;
152
153}
154
155
156// --------------------------------------------------------------------------
157//
158// The PreProcess searches for the following input containers:
159// - MRawEvtData
160// - MPedestalCam
161//
162// The following output containers are also searched and created if
163// they were not found:
164//
165// - MHCalibrationBlindPixel
166// - MCalibrationCam
167//
168// The following output containers are only searched, but not created
169//
170// - MTime
171//
172Int_t MCalibrationChargeCalc::PreProcess(MParList *pList)
173{
174
175 fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
176 if (!fRawEvt)
177 {
178 *fLog << err << dbginf << "MRawEvtData not found... aborting." << endl;
179 return kFALSE;
180 }
181
182 const MRawRunHeader *runheader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
183 if (!runheader)
184 *fLog << warn << dbginf << "Warning - cannot check file type, MRawRunHeader not found." << endl;
185 else
186 if (runheader->GetRunType() == kRTMonteCarlo)
187 {
188 return kTRUE;
189 }
190
191 fCam = (MCalibrationChargeCam*)pList->FindCreateObj("MCalibrationChargeCam");
192 if (!fCam)
193 {
194 *fLog << err << dbginf << "MCalibrationChargeCam could not be created ... aborting." << endl;
195 return kFALSE;
196 }
197
198 fPINDiode = (MCalibrationChargePINDiode*)pList->FindCreateObj("MCalibrationChargePINDiode");
199 if (!fPINDiode)
200 {
201 *fLog << err << dbginf << "MCalibrationChargePINDiode could not be created ... aborting." << endl;
202 return kFALSE;
203 }
204
205 fCam->SetPINDiode(fPINDiode);
206
207 fEvtTime = (MTime*)pList->FindObject("MTime");
208
209 fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam");
210 if (!fPedestals)
211 {
212 *fLog << err << dbginf << "Cannot find MPedestalCam ... aborting" << endl;
213 return kFALSE;
214 }
215
216
217 fSignals = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam");
218 if (!fSignals)
219 {
220 *fLog << err << dbginf << "Cannot find MExtractedSignalCam ... aborting" << endl;
221 return kFALSE;
222 }
223
224 fBlindPixel = (MExtractedSignalBlindPixel*)pList->FindObject("MExtractedSignalBlindPixel");
225 if (!fBlindPixel)
226 {
227 *fLog << err << dbginf << "Cannot find MExtractedSignalBlindPixel ... aborting" << endl;
228 return kFALSE;
229 }
230
231 return kTRUE;
232}
233
234
235// --------------------------------------------------------------------------
236//
237// The ReInit searches for the following input containers:
238// - MRawRunHeader
239//
240Bool_t MCalibrationChargeCalc::ReInit(MParList *pList )
241{
242
243 fRunHeader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
244 if (!fRunHeader)
245 {
246 *fLog << err << dbginf << ": MRawRunHeader not found... aborting." << endl;
247 return kFALSE;
248 }
249
250
251 MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam");
252 if (!cam)
253 {
254 *fLog << err << GetDescriptor() << ": No MGeomCam found... aborting." << endl;
255 return kFALSE;
256 }
257
258 fCam->SetGeomCam(cam);
259
260 fNumHiGainSamples = fSignals->GetNumUsedHiGainFADCSlices();
261 fNumLoGainSamples = fSignals->GetNumUsedLoGainFADCSlices();
262 fSqrtHiGainSamples = TMath::Sqrt((Float_t)fNumHiGainSamples);
263
264 UInt_t npixels = cam->GetNumPixels();
265
266 for (UInt_t i=0; i<npixels; i++)
267 {
268
269 MCalibrationChargePix &pix = (*fCam)[i];
270 pix.DefinePixId(i);
271
272 pix.SetAbsTimeBordersHiGain(fSignals->GetFirstUsedSliceHiGain(),
273 fSignals->GetLastUsedSliceHiGain());
274 pix.SetAbsTimeBordersLoGain(fSignals->GetFirstUsedSliceLoGain(),
275 fSignals->GetLastUsedSliceLoGain());
276
277
278 // Exclude the blind pixel and the PIN Diode from normal pixel calibration:
279 if (i == fBlindPixelIdx)
280 pix.SetExcluded();
281
282 if (i == fPINDiodeIdx)
283 pix.SetExcluded();
284
285 }
286
287 //
288 // Look for file to exclude pixels from analysis
289 //
290 if (!fExcludedPixelsFile.IsNull())
291 {
292
293 fExcludedPixelsFile = gSystem->ExpandPathName(fExcludedPixelsFile.Data());
294
295 //
296 // Initialize reading the file
297 //
298 ifstream in(fExcludedPixelsFile.Data(),ios::in);
299
300 if (in)
301 {
302 *fLog << inf << "Use excluded pixels from file: '" << fExcludedPixelsFile.Data() << "'" << endl;
303 //
304 // Read the file and count the number of entries
305 //
306 UInt_t pixel = 0;
307
308 while (++fNumExcludedPixels)
309 {
310
311 in >> pixel;
312
313 if (!in.good())
314 break;
315 //
316 // Check for out of range
317 //
318 if (pixel > npixels)
319 {
320 *fLog << warn << "WARNING: To be excluded pixel: " << pixel
321 << " is out of range " << endl;
322 continue;
323 }
324 //
325 // Exclude pixel
326 //
327 MCalibrationChargePix &pix = (*fCam)[pixel];
328 pix.SetExcluded();
329
330 *fLog << GetDescriptor() << inf << ": Exclude Pixel: " << pixel << endl;
331
332 }
333
334 if (--fNumExcludedPixels == 0)
335 *fLog << warn << "WARNING: File '" << fExcludedPixelsFile.Data()
336 << "'" << " is empty " << endl;
337 else
338 fCam->SetNumPixelsExcluded(fNumExcludedPixels);
339
340 }
341 else
342 *fLog << warn << dbginf << "Cannot open file '" << fExcludedPixelsFile.Data() << "'" << endl;
343 }
344
345 return kTRUE;
346}
347
348
349// --------------------------------------------------------------------------
350//
351// Calculate the integral of the FADC time slices and store them as a new
352// pixel in the MCerPhotEvt container.
353//
354Int_t MCalibrationChargeCalc::Process()
355{
356
357
358 MRawEvtPixelIter pixel(fRawEvt);
359
360 //
361 // Create a loop to fill the calibration histograms
362 // Search for: a signal in MExtractedSignalCam
363 // Search for: a signal in MExtractedSignalBlindPixel
364 // Fill histograms with:
365 // charge
366 // charge vs. event nr.
367 //
368 //
369 // Initialize pointers to blind pixel and individual pixels
370 //
371 // FIXME: filling the bind pixel histograms in this class is preliminary
372 // and will be replaced soon to fill them with MFillH
373 //
374 MCalibrationBlindPix &blindpixel = *(fCam->GetBlindPixel());
375 MExtractedSignalBlindPixel &blindsig = (*fBlindPixel);
376
377 const UInt_t signal = blindsig.GetExtractedSignal();
378
379 if (!blindpixel.FillCharge(signal))
380 *fLog << warn <<
381 "Overflow or Underflow occurred filling Blind Pixel sum = " << signal << endl;
382
383 blindpixel.FillGraphs(signal,0);
384
385 TH1I *hist;
386
387 if (signal > fBlindPixelSinglePheCut)
388 {
389 hist = (blindpixel.GetHist())->GetHSinglePheFADCSlices();
390 fNumBlindPixelSinglePhe++;
391 }
392 else
393 {
394 hist = (blindpixel.GetHist())->GetHPedestalFADCSlices();
395 fNumBlindPixelPedestal++;
396 }
397
398 pixel.Jump(fBlindPixelIdx);
399
400 const Byte_t *ptr = pixel.GetHiGainSamples();
401
402 for (Int_t i=1;i<16;i++)
403 hist->Fill(i,*ptr++);
404
405 ptr = pixel.GetLoGainSamples();
406 for (Int_t i=16;i<31;i++)
407 hist->Fill(i,*ptr++);
408
409 return kTRUE;
410}
411
412Int_t MCalibrationChargeCalc::PostProcess()
413{
414 *fLog << inf << GetDescriptor() << ": Cut Histogram Edges" << endl;
415
416 //
417 // Fit the blind pixel
418 //
419 if (TESTBIT(fFlags,kUseBlindPixelFit))
420 {
421 //
422 // Get pointer to blind pixel
423 //
424 MCalibrationBlindPix &blindpixel = *(fCam->GetBlindPixel());
425
426 *fLog << inf << GetDescriptor() << ": Fitting the Blind Pixel" << endl;
427
428 //
429 // retrieve the histogram containers
430 //
431 MHCalibrationBlindPixel *hist = blindpixel.GetHist();
432
433 //
434 // retrieve mean and sigma of the blind pixel pedestal,
435 // so that we can use it for the fit
436 //
437 // FIXME: This part has to be migrated into the future MHCalibrationChargeBlindPix class
438 //
439 const UInt_t nentries = fPedestals->GetTotalEntries();
440 const UInt_t nslices = 12;
441 const Float_t sqrslice = TMath::Sqrt((Float_t)nslices);
442
443 MPedestalPix &pedpix = (*fPedestals)[fBlindPixelIdx];
444 if (&pedpix)
445 {
446 const Float_t pedestal = pedpix.GetPedestal()*nslices;
447 const Float_t pederr = pedpix.GetPedestalRms()*nslices/nentries;
448 const Float_t pedsigma = pedpix.GetPedestalRms()*sqrslice;
449 const Float_t pedsigmaerr = pederr/2.;
450
451 hist->SetMeanPedestal(pedestal);
452 hist->SetMeanPedestalErr(pederr);
453 hist->SetSigmaPedestal(pedsigma);
454 hist->SetSigmaPedestalErr(pedsigmaerr);
455 }
456
457 if (!blindpixel.FitCharge())
458 {
459 *fLog << warn << "Could not fit the blind pixel! " << endl;
460 *fLog << warn << "Setting bit kBlindPixelMethodValid to FALSE in MCalibrationCam" << endl;
461 fCam->SetBlindPixelMethodValid(kFALSE);
462 }
463 else
464 fCam->SetBlindPixelMethodValid(kTRUE);
465
466 if (blindpixel.CheckOscillations())
467 fCam->SetBlindPixelMethodValid(kFALSE);
468
469 TH1I *sphehist = hist->GetHSinglePheFADCSlices();
470 TH1I *pedhist = hist->GetHPedestalFADCSlices();
471
472 if (fNumBlindPixelSinglePhe > 1)
473 sphehist->Scale(1./fNumBlindPixelSinglePhe);
474 if (fNumBlindPixelPedestal > 1)
475 pedhist->Scale(1./fNumBlindPixelPedestal);
476
477 blindpixel.DrawClone();
478 }
479 else
480 *fLog << inf << GetDescriptor() << ": Skipping Blind Pixel Fit " << endl;
481
482 *fLog << inf << "total: " << GetNumExecutions() << " sphe: " << fNumBlindPixelSinglePhe << " ped: " << fNumBlindPixelPedestal << endl;
483
484
485 *fLog << inf << GetDescriptor() << ": Fitting the Normal Pixels" << endl;
486
487 //
488 // loop over the pedestal events and check if we have calibration
489 //
490 for (Int_t pixid=0; pixid<fPedestals->GetSize(); pixid++)
491 {
492
493 MCalibrationChargePix &pix = (*fCam)[pixid];
494
495 //
496 // Check if the pixel has been excluded from the fits
497 //
498 if (pix.IsExcluded())
499 continue;
500
501 //
502 // get the pedestals
503 //
504 const Float_t ped = (*fPedestals)[pixid].GetPedestal();
505 const Float_t prms = (*fPedestals)[pixid].GetPedestalRms();
506 const Int_t num = fPedestals->GetTotalEntries();
507
508 //
509 // set them in the calibration camera
510 //
511 if (pix.IsHiGainSaturation())
512 {
513 pix.SetPedestal(ped * (Float_t)fNumLoGainSamples,
514 prms * TMath::Sqrt((Float_t)fNumLoGainSamples),
515 prms * (Float_t)fNumLoGainSamples / num);
516 pix.SetNumLoGainSamples((Float_t)fNumLoGainSamples);
517 pix.ApplyLoGainConversion();
518 }
519 else
520 {
521 pix.SetPedestal(ped * (Float_t)fNumHiGainSamples,
522 prms * TMath::Sqrt((Float_t)fNumHiGainSamples),
523 prms * (Float_t)fNumHiGainSamples / num);
524 }
525
526 if (!pix.CheckChargeValidity() || !pix.CheckTimeValidity())
527 continue;
528
529 if (!pix.CalcReducedSigma())
530 continue;
531
532 pix.CalcFFactorMethod();
533
534 }
535
536 if (TESTBIT(fFlags,kUseBlindPixelFit) && fCam->IsBlindPixelMethodValid())
537 {
538 if (!fCam->CalcFluxInsidePlexiglass())
539 {
540 *fLog << warn << "Could not calculate the number of photons from the blind pixel " << endl;
541 *fLog << "You can try to calibrate using the MCalibrationChargeCalc::SkipBlindPixelFit()" << endl;
542 fCam->SetBlindPixelMethodValid(kFALSE);
543 }
544 }
545 else
546 *fLog << inf << GetDescriptor() << ": Skipping Blind Pixel Calibration! " << endl;
547
548 if (!fPINDiode->CheckChargeFitValidity() || !fPINDiode->CheckTimeFitValidity())
549 {
550 *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, charge fit not valid " << endl;
551 fCam->SetPINDiodeMethodValid(kFALSE);
552 }
553 else
554 {
555 if (!fPINDiode->CalcFluxOutsidePlexiglass())
556 {
557 *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, will skip PIN Diode Calibration " << endl;
558 fCam->SetPINDiodeMethodValid(kFALSE);
559 }
560 else
561 {
562 fCam->SetPINDiodeMethodValid(kTRUE);
563 }
564 }
565 fCam->SetReadyToSave();
566
567 return kTRUE;
568}
569
570
571
572
Note: See TracBrowser for help on using the repository browser.