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

Last change on this file since 3253 was 3247, checked in by gaug, 21 years ago
*** empty log message ***
File size: 18.2 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 09/2003 <mailto:markus@ifae.es>
19!
20! Copyright: MAGIC Software Development, 2000-2001
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 MCalibrationPix
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 MCalibrationPix
42// Initialize number of used FADC slices
43// Optionally exclude pixels from calibration
44//
45// Process: Every MCalibrationPix 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 "MCalibrationPix.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 CLRBIT(fFlags, kHiGainOverFlow);
146 CLRBIT(fFlags, kLoGainOverFlow);
147
148 fNumBlindPixelSinglePhe = 0;
149 fNumBlindPixelPedestal = 0;
150
151 fNumHiGainSamples = 0;
152 fNumLoGainSamples = 0;
153 fConversionHiLo = 0;
154 fNumExcludedPixels = 0;
155
156}
157
158
159// --------------------------------------------------------------------------
160//
161// The PreProcess searches for the following input containers:
162// - MRawEvtData
163// - MPedestalCam
164//
165// The following output containers are also searched and created if
166// they were not found:
167//
168// - MHCalibrationBlindPixel
169// - MCalibrationCam
170//
171// The following output containers are only searched, but not created
172//
173// - MTime
174//
175Int_t MCalibrationChargeCalc::PreProcess(MParList *pList)
176{
177
178 fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
179 if (!fRawEvt)
180 {
181 *fLog << err << dbginf << "MRawEvtData not found... aborting." << endl;
182 return kFALSE;
183 }
184
185 const MRawRunHeader *runheader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
186 if (!runheader)
187 *fLog << warn << dbginf << "Warning - cannot check file type, MRawRunHeader not found." << endl;
188 else
189 if (runheader->GetRunType() == kRTMonteCarlo)
190 {
191 return kTRUE;
192 }
193
194 fCam = (MCalibrationChargeCam*)pList->FindCreateObj("MCalibrationChargeCam");
195 if (!fCam)
196 {
197 *fLog << err << dbginf << "MCalibrationChargeCam could not be created ... aborting." << endl;
198 return kFALSE;
199 }
200
201 fPINDiode = (MCalibrationChargePINDiode*)pList->FindCreateObj("MCalibrationChargePINDiode");
202 if (!fPINDiode)
203 {
204 *fLog << err << dbginf << "MCalibrationChargePINDiode could not be created ... aborting." << endl;
205 return kFALSE;
206 }
207
208 fCam->SetPINDiode(fPINDiode);
209
210 fEvtTime = (MTime*)pList->FindObject("MTime");
211
212 fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam");
213 if (!fPedestals)
214 {
215 *fLog << err << dbginf << "Cannot find MPedestalCam ... aborting" << endl;
216 return kFALSE;
217 }
218
219
220 fSignals = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam");
221 if (!fSignals)
222 {
223 *fLog << err << dbginf << "Cannot find MExtractedSignalCam ... aborting" << endl;
224 return kFALSE;
225 }
226
227 fBlindPixel = (MExtractedSignalBlindPixel*)pList->FindObject("MExtractedSignalBlindPixel");
228 if (!fBlindPixel)
229 {
230 *fLog << err << dbginf << "Cannot find MExtractedSignalBlindPixel ... aborting" << endl;
231 return kFALSE;
232 }
233
234 return kTRUE;
235}
236
237
238// --------------------------------------------------------------------------
239//
240// The ReInit searches for the following input containers:
241// - MRawRunHeader
242//
243Bool_t MCalibrationChargeCalc::ReInit(MParList *pList )
244{
245
246 fRunHeader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
247 if (!fRunHeader)
248 {
249 *fLog << err << dbginf << ": MRawRunHeader not found... aborting." << endl;
250 return kFALSE;
251 }
252
253
254 MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam");
255 if (!cam)
256 {
257 *fLog << err << GetDescriptor() << ": No MGeomCam found... aborting." << endl;
258 return kFALSE;
259 }
260
261 fCam->SetGeomCam(cam);
262
263 fNumHiGainSamples = fSignals->GetNumUsedHiGainFADCSlices();
264 fNumLoGainSamples = fSignals->GetNumUsedLoGainFADCSlices();
265 fSqrtHiGainSamples = TMath::Sqrt((Float_t)fNumHiGainSamples);
266
267 UInt_t npixels = cam->GetNumPixels();
268
269 for (UInt_t i=0; i<npixels; i++)
270 {
271
272 MCalibrationPix &pix = (*fCam)[i];
273 pix.DefinePixId(i);
274
275 pix.SetAbsTimeBordersHiGain(fSignals->GetFirstUsedSliceHiGain(),
276 fSignals->GetLastUsedSliceHiGain());
277 pix.SetAbsTimeBordersLoGain(fSignals->GetFirstUsedSliceLoGain(),
278 fSignals->GetLastUsedSliceLoGain());
279
280 if (!TESTBIT(fFlags,kUseQualityChecks))
281 pix.SetExcludeQualityCheck();
282
283 // Exclude the blind pixel and the PIN Diode from normal pixel calibration:
284 if (i == fBlindPixelIdx)
285 pix.SetExcluded();
286
287 if (i == fPINDiodeIdx)
288 pix.SetExcluded();
289
290 }
291
292 //
293 // Look for file to exclude pixels from analysis
294 //
295 if (!fExcludedPixelsFile.IsNull())
296 {
297
298 fExcludedPixelsFile = gSystem->ExpandPathName(fExcludedPixelsFile.Data());
299
300 //
301 // Initialize reading the file
302 //
303 ifstream in(fExcludedPixelsFile.Data(),ios::in);
304
305 if (in)
306 {
307 *fLog << inf << "Use excluded pixels from file: '" << fExcludedPixelsFile.Data() << "'" << endl;
308 //
309 // Read the file and count the number of entries
310 //
311 UInt_t pixel = 0;
312
313 while (++fNumExcludedPixels)
314 {
315
316 in >> pixel;
317
318 if (!in.good())
319 break;
320 //
321 // Check for out of range
322 //
323 if (pixel > npixels)
324 {
325 *fLog << warn << "WARNING: To be excluded pixel: " << pixel
326 << " is out of range " << endl;
327 continue;
328 }
329 //
330 // Exclude pixel
331 //
332 MCalibrationPix &pix = (*fCam)[pixel];
333 pix.SetExcluded();
334
335 *fLog << GetDescriptor() << inf << ": Exclude Pixel: " << pixel << endl;
336
337 }
338
339 if (--fNumExcludedPixels == 0)
340 *fLog << warn << "WARNING: File '" << fExcludedPixelsFile.Data()
341 << "'" << " is empty " << endl;
342 else
343 fCam->SetNumPixelsExcluded(fNumExcludedPixels);
344
345 }
346 else
347 *fLog << warn << dbginf << "Cannot open file '" << fExcludedPixelsFile.Data() << "'" << endl;
348 }
349
350 return kTRUE;
351}
352
353
354// --------------------------------------------------------------------------
355//
356// Calculate the integral of the FADC time slices and store them as a new
357// pixel in the MCerPhotEvt container.
358//
359Int_t MCalibrationChargeCalc::Process()
360{
361
362
363 MRawEvtPixelIter pixel(fRawEvt);
364
365 //
366 // Create a loop to fill the calibration histograms
367 // Search for: a signal in MExtractedSignalCam
368 // Search for: a signal in MExtractedSignalBlindPixel
369 // Fill histograms with:
370 // charge
371 // charge vs. event nr.
372 //
373 //
374 // Initialize pointers to blind pixel and individual pixels
375 //
376 // FIXME: filling the bind pixel histograms in this class is preliminary
377 // and will be replaced soon to fill them with MFillH
378 //
379 MCalibrationBlindPix &blindpixel = *(fCam->GetBlindPixel());
380 MExtractedSignalBlindPixel &blindsig = (*fBlindPixel);
381
382 const UInt_t signal = blindsig.GetExtractedSignal();
383
384 if (!blindpixel.FillCharge(signal))
385 *fLog << warn <<
386 "Overflow or Underflow occurred filling Blind Pixel sum = " << signal << endl;
387
388 blindpixel.FillGraphs(signal,0);
389
390 TH1I *hist;
391
392 if (signal > fBlindPixelSinglePheCut)
393 {
394 hist = (blindpixel.GetHist())->GetHSinglePheFADCSlices();
395 fNumBlindPixelSinglePhe++;
396 }
397 else
398 {
399 hist = (blindpixel.GetHist())->GetHPedestalFADCSlices();
400 fNumBlindPixelPedestal++;
401 }
402
403 pixel.Jump(fBlindPixelIdx);
404
405 const Byte_t *ptr = pixel.GetHiGainSamples();
406
407 for (Int_t i=1;i<16;i++)
408 hist->Fill(i,*ptr++);
409
410 ptr = pixel.GetLoGainSamples();
411 for (Int_t i=16;i<31;i++)
412 hist->Fill(i,*ptr++);
413
414 pixel.Reset();
415
416 while (pixel.Next())
417 {
418
419 const UInt_t pixid = pixel.GetPixelId();
420
421 MCalibrationPix &pix = (*fCam)[pixid];
422 MExtractedSignalPix &sig = (*fSignals) [pixid];
423
424 const Float_t sumhi = sig.GetExtractedSignalHiGain();
425 const Float_t sumlo = sig.GetExtractedSignalLoGain();
426
427 Float_t abstime = 0.;
428
429 if (sig.IsLoGainUsed())
430 abstime = (Float_t)pixel.GetIdxMaxLoGainSample();
431 else
432 abstime = (Float_t)pixel.GetIdxMaxHiGainSample();
433
434 if (pix.IsExcluded())
435 continue;
436
437 pix.FillGraphs(sumhi,sumlo);
438
439 if (sig.IsLoGainUsed())
440 {
441
442 if (!pix.FillChargeLoGain(sumlo))
443 *fLog << warn << "Could not fill Lo Gain Charge of pixel: " << pixid
444 << " signal = " << sumlo << endl;
445
446 if (!pix.FillAbsTimeLoGain(abstime))
447 *fLog << warn << "Could not fill Lo Gain Abs. Time of pixel: "
448 << pixid << " time = " << abstime << endl;
449 }
450 else
451 {
452 if (!pix.FillChargeHiGain(sumhi))
453 *fLog << warn << "Could not fill Hi Gain Charge of pixel: " << pixid
454 << " signal = " << sumhi << endl;
455
456 if (!pix.FillAbsTimeHiGain(abstime))
457 *fLog << warn << "Could not fill Hi Gain Abs. Time of pixel: "
458 << pixid << " time = " << abstime << endl;
459 }
460
461 } /* while (pixel.Next()) */
462
463 return kTRUE;
464}
465
466Int_t MCalibrationChargeCalc::PostProcess()
467{
468 *fLog << inf << GetDescriptor() << ": Cut Histogram Edges" << endl;
469
470 //
471 // Cut edges to make fits and viewing of the hists easier
472 //
473 fCam->CutEdges();
474
475 //
476 // Fit the blind pixel
477 //
478 if (TESTBIT(fFlags,kUseBlindPixelFit))
479 {
480 //
481 // Get pointer to blind pixel
482 //
483 MCalibrationBlindPix &blindpixel = *(fCam->GetBlindPixel());
484
485 *fLog << inf << GetDescriptor() << ": Fitting the Blind Pixel" << endl;
486
487 //
488 // retrieve the histogram containers
489 //
490 MHCalibrationBlindPixel *hist = blindpixel.GetHist();
491
492 //
493 // retrieve mean and sigma of the blind pixel pedestal,
494 // so that we can use it for the fit
495 //
496 // FIXME: This part has to be migrated into the future MHCalibrationChargeBlindPix class
497 //
498 const UInt_t nentries = fPedestals->GetTotalEntries();
499 const UInt_t nslices = 12;
500 const Float_t sqrslice = TMath::Sqrt((Float_t)nslices);
501
502 MPedestalPix &pedpix = (*fPedestals)[fBlindPixelIdx];
503 if (&pedpix)
504 {
505 const Float_t pedestal = pedpix.GetPedestal()*nslices;
506 const Float_t pederr = pedpix.GetPedestalRms()*nslices/nentries;
507 const Float_t pedsigma = pedpix.GetPedestalRms()*sqrslice;
508 const Float_t pedsigmaerr = pederr/2.;
509
510 hist->SetMeanPedestal(pedestal);
511 hist->SetMeanPedestalErr(pederr);
512 hist->SetSigmaPedestal(pedsigma);
513 hist->SetSigmaPedestalErr(pedsigmaerr);
514 }
515
516 if (!blindpixel.FitCharge())
517 {
518 *fLog << warn << "Could not fit the blind pixel! " << endl;
519 *fLog << warn << "Setting bit kBlindPixelMethodValid to FALSE in MCalibrationCam" << endl;
520 fCam->SetBlindPixelMethodValid(kFALSE);
521 }
522 else
523 fCam->SetBlindPixelMethodValid(kTRUE);
524
525 if (blindpixel.CheckOscillations())
526 fCam->SetBlindPixelMethodValid(kFALSE);
527
528 TH1I *sphehist = hist->GetHSinglePheFADCSlices();
529 TH1I *pedhist = hist->GetHPedestalFADCSlices();
530
531 if (fNumBlindPixelSinglePhe > 1)
532 sphehist->Scale(1./fNumBlindPixelSinglePhe);
533 if (fNumBlindPixelPedestal > 1)
534 pedhist->Scale(1./fNumBlindPixelPedestal);
535
536 blindpixel.DrawClone();
537 }
538 else
539 *fLog << inf << GetDescriptor() << ": Skipping Blind Pixel Fit " << endl;
540
541 *fLog << inf << "total: " << GetNumExecutions() << " sphe: " << fNumBlindPixelSinglePhe << " ped: " << fNumBlindPixelPedestal << endl;
542
543
544 *fLog << inf << GetDescriptor() << ": Fitting the Normal Pixels" << endl;
545
546 //
547 // loop over the pedestal events and check if we have calibration
548 //
549 for (Int_t pixid=0; pixid<fPedestals->GetSize(); pixid++)
550 {
551
552 MCalibrationPix &pix = (*fCam)[pixid];
553
554 //
555 // Check if the pixel has been excluded from the fits
556 //
557 if (pix.IsExcluded())
558 continue;
559
560 //
561 // get the pedestals
562 //
563 const Float_t ped = (*fPedestals)[pixid].GetPedestal();
564 const Float_t prms = (*fPedestals)[pixid].GetPedestalRms();
565
566 //
567 // set them in the calibration camera
568 //
569 pix.SetPedestal(ped,prms,(Float_t)fNumHiGainSamples,(Float_t)fNumLoGainSamples);
570
571 //
572 // perform the Gauss fits to the charges
573 //
574 pix.FitCharge();
575
576 //
577 // check also for oscillations
578 //
579 pix.CheckOscillations();
580
581 //
582 // calculate the F-Factor method
583 //
584 pix.CalcFFactorMethod();
585 }
586
587 if (TESTBIT(fFlags,kUseBlindPixelFit) && fCam->IsBlindPixelMethodValid())
588 {
589 if (!fCam->CalcFluxInsidePlexiglass())
590 {
591 *fLog << warn << "Could not calculate the number of photons from the blind pixel " << endl;
592 *fLog << "You can try to calibrate using the MCalibrationChargeCalc::SkipBlindPixelFit()" << endl;
593 fCam->SetBlindPixelMethodValid(kFALSE);
594 }
595 }
596 else
597 *fLog << inf << GetDescriptor() << ": Skipping Blind Pixel Calibration! " << endl;
598
599 if (!fPINDiode->CalcFluxOutsidePlexiglass())
600 {
601 *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, will skip PIN Diode Calibration " << endl;
602 fCam->SetPINDiodeMethodValid(kFALSE);
603 }
604 else
605 {
606 fCam->SetPINDiodeMethodValid(kTRUE);
607 }
608
609 fCam->SetReadyToSave();
610
611 return kTRUE;
612}
613
614
Note: See TracBrowser for help on using the repository browser.