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

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