source: trunk/MagicSoft/Mars/mcalib/MCalibrationCalc.cc@ 3243

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