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

Last change on this file since 2987 was 2972, checked in by gaug, 21 years ago
*** empty log message ***
File size: 21.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// 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 and possibly
32// arrival times from MArrivalTime
33//
34// The output container MCalibrationCam holds one entry of type MCalibrationPix
35// for every pixel. It is filled in the following way:
36//
37// ProProcess: Search for MPedestalCam, MExtractedSignalCam
38// Initialize MCalibrationCam
39// Initialize MArrivalTime if exists
40// Initialize pulser light wavelength
41//
42// ReInit: MCalibrationCam::InitSize(NumPixels) is called which allocates
43// memory in a TClonesArray of type MCalibrationPix
44// Initialize number of used FADC slices
45// Optionally exclude pixels from calibration
46//
47// Process: Optionally, a cut on cosmics can be performed
48//
49// Every MCalibrationPix holds a histogram class,
50// MHCalibrationPixel which itself hold histograms of type:
51// HCharge(npix) (distribution of summed FADC time slice
52// entries)
53// HTime(npix) (distribution of position of maximum)
54// HChargevsN(npix) (distribution of charges vs. event number.
55//
56// PostProcess: All histograms HCharge(npix) are fitted to a Gaussian
57// All histograms HTime(npix) are fitted to a Gaussian
58// The histogram HBlindPixelCharge (blind pixel) is fitted to
59// a single PhE fit
60//
61// The histograms of the PIN Diode are fitted to Gaussians
62//
63// Fits can be excluded via the commands:
64// MalibrationCam::SkipTimeFits() (skip all time fits)
65// MalibrationCam::SkipBlindPixelFits() (skip all blind
66// pixel fits)
67// MalibrationCam::SkipPinDiodeFits() (skip all PIN Diode
68// fits)
69//
70// Input Containers:
71// MRawEvtData
72// MPedestalCam
73// MExtractedSignalCam
74//
75// Output Containers:
76// MCalibrationCam
77//
78//////////////////////////////////////////////////////////////////////////////
79#include "MCalibrationCalc.h"
80
81// FXIME: Usage of fstream is a preliminary workaround!
82#include <fstream>
83
84// FXIME: This has to be removed!!!! (YES, WHEN WE HAVE ACCESS TO THE DATABASE!!!!!)
85#include "MCalibrationConfig.h"
86
87#include <TSystem.h>
88
89#include "MLog.h"
90#include "MLogManip.h"
91
92#include "MParList.h"
93
94#include "MGeomCam.h"
95#include "MRawRunHeader.h"
96#include "MRawEvtPixelIter.h"
97
98#include "MPedestalCam.h"
99#include "MPedestalPix.h"
100
101#include "MCalibrationCam.h"
102#include "MCalibrationPix.h"
103
104#include "MExtractedSignalCam.h"
105#include "MExtractedSignalPix.h"
106
107#include "MCalibrationBlindPix.h"
108#include "MCalibrationPINDiode.h"
109
110#include "MArrivalTime.h"
111
112ClassImp(MCalibrationCalc);
113
114using namespace std;
115
116const Int_t MCalibrationCalc::fBlindPixelId = 559;
117const Int_t MCalibrationCalc::fPINDiodeId = 9999;
118
119// --------------------------------------------------------------------------
120//
121// Default constructor.
122//
123MCalibrationCalc::MCalibrationCalc(const char *name, const char *title)
124 : fPedestals(NULL), fCalibrations(NULL), fSignals(NULL),
125 fRawEvt(NULL), fRunHeader(NULL), fArrivalTime(NULL), fEvtTime(NULL)
126{
127
128 fName = name ? name : "MCalibrationCalc";
129 fTitle = title ? title : "Task to calculate the calibration constants and MCalibrationCam ";
130
131 AddToBranchList("MRawEvtData.fHiGainPixId");
132 AddToBranchList("MRawEvtData.fLoGainPixId");
133 AddToBranchList("MRawEvtData.fHiGainFadcSamples");
134 AddToBranchList("MRawEvtData.fLoGainFadcSamples");
135
136 Clear();
137}
138
139void MCalibrationCalc::Clear(const Option_t *o)
140{
141
142 SETBIT(fFlags, kUseTimes);
143 SETBIT(fFlags, kUseBlindPixelFit);
144 SETBIT(fFlags, kUseCosmicsRejection);
145 SETBIT(fFlags, kUseQualityChecks);
146
147 // As long as we don't have the PIN Diode:
148 CLRBIT(fFlags, kUsePinDiodeFit);
149
150 fEvents = 0;
151 fCosmics = 0;
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
166MCalibrationPINDiode *MCalibrationCalc::GetPINDiode() const
167{
168 return fCalibrations->GetPINDiode();
169}
170
171// --------------------------------------------------------------------------
172//
173// The PreProcess searches for the following input containers:
174// - MRawEvtData
175// - MPedestalCam
176//
177// The following output containers are also searched and created if
178// they were not found:
179//
180// - MHCalibrationBlindPixel
181// - MCalibrationCam
182//
183// The following output containers are only searched, but not created
184//
185// - MArrivaltime
186// - MTime
187//
188Int_t MCalibrationCalc::PreProcess(MParList *pList)
189{
190
191 fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
192 if (!fRawEvt)
193 {
194 *fLog << err << dbginf << "MRawEvtData not found... aborting." << endl;
195 return kFALSE;
196 }
197
198 const MRawRunHeader *runheader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
199 if (!runheader)
200 *fLog << warn << dbginf << "Warning - cannot check file type, MRawRunHeader not found." << endl;
201 else
202 if (runheader->GetRunType() == kRTMonteCarlo)
203 {
204 return kTRUE;
205 }
206
207 fCalibrations = (MCalibrationCam*)pList->FindCreateObj("MCalibrationCam");
208 if (!fCalibrations)
209 {
210 *fLog << err << dbginf << "MCalibrationCam could not be created ... aborting." << endl;
211 return kFALSE;
212 }
213
214 fArrivalTime = (MArrivalTime*)pList->FindObject("MArrivalTime");
215
216 if (!fArrivalTime)
217 CLRBIT(fFlags,kUseTimes);
218
219 fEvtTime = (MTime*)pList->FindObject("MTime");
220
221 switch (fColor)
222 {
223 case kEBlue:
224 fCalibrations->SetColor(MCalibrationCam::kECBlue);
225 break;
226 case kEGreen:
227 fCalibrations->SetColor(MCalibrationCam::kECGreen);
228 break;
229 case kEUV:
230 fCalibrations->SetColor(MCalibrationCam::kECUV);
231 break;
232 case kECT1:
233 fCalibrations->SetColor(MCalibrationCam::kECCT1);
234 break;
235 default:
236 fCalibrations->SetColor(MCalibrationCam::kECCT1);
237 }
238
239 fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam");
240 if (!fPedestals)
241 {
242 *fLog << err << dbginf << "Cannot find MPedestalCam ... aborting" << endl;
243 return kFALSE;
244 }
245
246
247 fSignals = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam");
248 if (!fSignals)
249 {
250 *fLog << err << dbginf << "Cannot find MExtractedSignalCam ... aborting" << endl;
251 return kFALSE;
252 }
253
254 return kTRUE;
255}
256
257
258// --------------------------------------------------------------------------
259//
260// The ReInit searches for the following input containers:
261// - MRawRunHeader
262//
263Bool_t MCalibrationCalc::ReInit(MParList *pList )
264{
265
266 fRunHeader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
267 if (!fRunHeader)
268 {
269 *fLog << err << dbginf << ": MRawRunHeader not found... aborting." << endl;
270 return kFALSE;
271 }
272
273
274 MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam");
275 if (!cam)
276 {
277 *fLog << err << GetDescriptor() << ": No MGeomCam found... aborting." << endl;
278 return kFALSE;
279 }
280
281 fNumHiGainSamples = fSignals->GetNumUsedHiGainFADCSlices();
282 fNumLoGainSamples = fSignals->GetNumUsedLoGainFADCSlices();
283 fSqrtHiGainSamples = TMath::Sqrt((Float_t)fNumHiGainSamples);
284
285 UInt_t npixels = cam->GetNumPixels();
286
287 for (UInt_t i=0; i<npixels; i++)
288 {
289
290 MCalibrationPix &pix = (*fCalibrations)[i];
291 pix.DefinePixId(i);
292 MHCalibrationPixel *hist = pix.GetHist();
293
294 hist->SetTimeFitRangesHiGain(fSignals->GetFirstUsedSliceHiGain(),
295 fSignals->GetLastUsedSliceHiGain());
296 hist->SetTimeFitRangesLoGain(fSignals->GetFirstUsedSliceLoGain(),
297 fSignals->GetLastUsedSliceLoGain());
298
299 if (!TESTBIT(fFlags,kUseQualityChecks))
300 pix.SetExcludeQualityCheck();
301
302 }
303
304 //
305 // Look for file to exclude pixels from analysis
306 //
307 if (!fExcludedPixelsFile.IsNull())
308 {
309
310 fExcludedPixelsFile = gSystem->ExpandPathName(fExcludedPixelsFile.Data());
311
312 //
313 // Initialize reading the file
314 //
315 ifstream in(fExcludedPixelsFile.Data(),ios::in);
316
317 if (in)
318 {
319 *fLog << inf << "Use excluded pixels from file: '" << fExcludedPixelsFile.Data() << "'" << endl;
320 //
321 // Read the file and count the number of entries
322 //
323 UInt_t pixel = 0;
324
325 while (++fNumExcludedPixels)
326 {
327
328 in >> pixel;
329
330 if (!in.good())
331 break;
332 //
333 // Check for out of range
334 //
335 if (pixel > npixels)
336 {
337 *fLog << warn << "WARNING: To be excluded pixel: " << pixel
338 << " is out of range " << endl;
339 continue;
340 }
341 //
342 // Exclude pixel
343 //
344 MCalibrationPix &pix = (*fCalibrations)[pixel];
345 pix.SetExcluded();
346
347 *fLog << GetDescriptor() << inf << ": Exclude Pixel: " << pixel << endl;
348
349 }
350
351 if (--fNumExcludedPixels == 0)
352 *fLog << warn << "WARNING: File '" << fExcludedPixelsFile.Data()
353 << "'" << " is empty " << endl;
354 else
355 fCalibrations->SetNumPixelsExcluded(fNumExcludedPixels);
356
357 }
358 else
359 *fLog << warn << dbginf << "Cannot open file '" << fExcludedPixelsFile.Data() << "'" << endl;
360 }
361
362 return kTRUE;
363}
364
365
366// --------------------------------------------------------------------------
367//
368// Calculate the integral of the FADC time slices and store them as a new
369// pixel in the MCerPhotEvt container.
370//
371Int_t MCalibrationCalc::Process()
372{
373
374 //
375 // Initialize pointers to blind pixel, PIN Diode and individual pixels
376 //
377 MCalibrationBlindPix &blindpixel = *(fCalibrations->GetBlindPixel());
378 MCalibrationPINDiode &pindiode = *(fCalibrations->GetPINDiode());
379
380 MRawEvtPixelIter pixel(fRawEvt);
381
382 //
383 // Perform cosmics cut
384 //
385 if (TESTBIT(fFlags,kUseCosmicsRejection))
386 {
387
388 Int_t cosmicpix = 0;
389
390 //
391 // Create a first loop to sort out the cosmics ...
392 //
393 // This is a very primitive check for the number of cosmicpixs
394 // The cut will be applied in the fit, but for the blind pixel,
395 // we need to remove this event
396 //
397 // FIXME: In the future need a much more sophisticated one!!!
398 //
399
400 while (pixel.Next())
401 {
402
403 const UInt_t pixid = pixel.GetPixelId();
404
405 MExtractedSignalPix &sig = (*fSignals)[pixid];
406 MPedestalPix &ped = (*fPedestals)[pixid];
407 Float_t pedrms = ped.GetPedestalRms()*fSqrtHiGainSamples;
408 Float_t sumhi = sig.GetExtractedSignalHiGain();
409
410 //
411 // We consider a pixel as presumably due to cosmics
412 // if its sum of FADC slices is lower than 3 pedestal RMS
413 //
414 if (sumhi < 3.*pedrms )
415 cosmicpix++;
416 }
417
418 //
419 // If the camera contains more than 230
420 // (this is the number of outer pixels plus about 50 inner ones)
421 // presumed pixels due to cosmics, then the event is discarted.
422 // This procedure is more or less equivalent to keeping only events
423 // with at least 350 pixels with high signals.
424 //
425 if (cosmicpix > 230.)
426 {
427 fCosmics++;
428 return kCONTINUE;
429 }
430
431 pixel.Reset();
432 }
433
434 fEvents++;
435
436 Int_t overflow = 0;
437
438 //
439 // Create a (second) loop to do fill the calibration histograms
440 //
441
442 while (pixel.Next())
443 {
444
445 const UInt_t pixid = pixel.GetPixelId();
446
447 MCalibrationPix &pix = (*fCalibrations)[pixid];
448
449 if (pix.IsExcluded())
450 continue;
451
452 MExtractedSignalPix &sig = (*fSignals)[pixid];
453
454 const Float_t sumhi = sig.GetExtractedSignalHiGain();
455 const Float_t sumlo = sig.GetExtractedSignalLoGain();
456
457 Double_t mtime = 0.;
458
459 if (TESTBIT(fFlags,kUseTimes))
460 mtime = (*fArrivalTime)[pixid];
461
462 switch(pixid)
463 {
464
465 case fBlindPixelId:
466
467 if (!blindpixel.FillCharge(sumhi))
468 *fLog << warn <<
469 "Overflow or Underflow occurred filling Blind Pixel sum = " << sumhi << endl;
470
471 if (TESTBIT(fFlags,kUseTimes))
472 {
473 if (!blindpixel.FillTime(mtime))
474 *fLog << warn <<
475 "Overflow or Underflow occurred filling Blind Pixel time = " << mtime << endl;
476 }
477
478 if (!blindpixel.FillRChargevsTime(sumhi,fEvents))
479 *fLog << warn <<
480 "Overflow or Underflow occurred filling Blind Pixel eventnr = " << fEvents << endl;
481 break;
482
483 case fPINDiodeId:
484
485 if (!pindiode.FillCharge(sumhi))
486 *fLog << warn <<
487 "Overflow or Underflow occurred filling PINDiode: sum = " << sumhi << endl;
488
489 if (TESTBIT(fFlags,kUseTimes))
490 {
491 if (!pindiode.FillTime(mtime))
492 *fLog << warn <<
493 "Overflow or Underflow occurred filling PINDiode: time = " << mtime << endl;
494 }
495
496 if (!pindiode.FillRChargevsTime(sumhi,fEvents))
497 *fLog << warn <<
498 "Overflow or Underflow occurred filling PINDiode: eventnr = " << fEvents << endl;
499 break;
500
501 default:
502
503 pix.FillChargesInGraph(sumhi,sumlo);
504
505 if (!pix.FillRChargevsTimeLoGain(sumlo,fEvents))
506 overflow++;
507
508 if (!pix.FillRChargevsTimeHiGain(sumhi,fEvents))
509 overflow++;
510
511 if (sig.IsLoGainUsed())
512 {
513
514 if (!pix.FillChargeLoGain(sumlo))
515 *fLog << warn << "Could not fill Lo Gain Charge of pixel: " << pixid
516 << " signal = " << sumlo << endl;
517
518 if (TESTBIT(fFlags,kUseTimes))
519 {
520 if (!pix.FillTimeLoGain(mtime))
521 *fLog << warn << "Could not fill Lo Gain Time of pixel: "
522 << pixid << " time = " << mtime << endl;
523 }
524 } /* if (sig.IsLoGainUsed()) */
525 else
526 {
527 if (!pix.FillChargeHiGain(sumhi))
528 *fLog << warn << "Could not fill Hi Gain Charge of pixel: " << pixid
529 << " signal = " << sumhi << endl;
530
531 if (TESTBIT(fFlags,kUseTimes))
532 {
533 if (!pix.FillTimeHiGain(mtime))
534 *fLog << warn << "Could not fill Hi Gain Time of pixel: "
535 << pixid << " time = " << mtime << endl;
536 }
537 } /* else (sig.IsLoGainUsed()) */
538 break;
539
540 } /* switch(pixid) */
541
542 } /* while (pixel.Next()) */
543
544 if (overflow)
545 *fLog << warn << "Overflow occurred filling Charges vs. EvtNr " << overflow << " times" << endl;
546
547 return kTRUE;
548}
549
550Int_t MCalibrationCalc::PostProcess()
551{
552
553 *fLog << inf << endl;
554
555 if (fEvents == 0)
556 {
557 *fLog << err << GetDescriptor()
558 << ": This run contains only cosmics or pedestals, "
559 << "cannot find events with more than 350 illuminated pixels. " << endl;
560 return kFALSE;
561 }
562
563 if (fEvents < fCosmics)
564 *fLog << warn << GetDescriptor()
565 << ": WARNING: Run contains more cosmics or pedestals than calibration events " << endl;
566
567
568 *fLog << inf << GetDescriptor() << ": Cut Histogram Edges" << endl;
569
570 //
571 // Cut edges to make fits and viewing of the hists easier
572 //
573 fCalibrations->CutEdges();
574
575 //
576 // Fit the blind pixel
577 //
578 if (TESTBIT(fFlags,kUseBlindPixelFit))
579 {
580 //
581 // Get pointer to blind pixel
582 //
583 MCalibrationBlindPix &blindpixel = *(fCalibrations->GetBlindPixel());
584
585 *fLog << inf << GetDescriptor() << ": Fitting the Blind Pixel" << endl;
586
587 //
588 // retrieve mean and sigma of the blind pixel pedestal,
589 // so that we can use it for the fit
590 //
591 if (fPedestals->IsUseHists())
592 {
593 //
594 // retrieve the pedestal pix of the blind pixel
595 //
596 MPedestalPix &ped = (*fPedestals)[fBlindPixelId];
597 //
598 // retrieve the blind pixel histogram container
599 //
600 MHCalibrationBlindPixel *hist = blindpixel.GetHist();
601 //
602 // Set the corresponding values
603 //
604 const Float_t peddiff = ped.GetMean()
605 - ped.GetPedestal()*fSignals->GetNumUsedFADCSlices();
606 const Float_t pederr = ped.GetMeanErr();
607 const Float_t pedsigma = ped.GetSigma();
608 const Float_t pedsigmaerr = ped.GetSigmaErr();
609
610 hist->SetMeanPedestal(peddiff);
611 hist->SetMeanPedestalErr(pederr);
612 hist->SetSigmaPedestal(pedsigma);
613 hist->SetSigmaPedestalErr(pedsigmaerr);
614 }
615
616 if (!blindpixel.FitCharge())
617 {
618 *fLog << warn << "Could not fit the blind pixel! " << endl;
619 *fLog << warn << "Setting bit kBlindPixelMethodValid to FALSE in MCalibrationCam" << endl;
620 fCalibrations->SetBlindPixelMethodValid(kFALSE);
621 }
622
623 fCalibrations->SetBlindPixelMethodValid(kTRUE);
624 blindpixel.DrawClone();
625 }
626 else
627 *fLog << inf << GetDescriptor() << ": Skipping Blind Pixel Fit " << endl;
628
629
630 *fLog << inf << GetDescriptor() << ": Fitting the Normal Pixels" << endl;
631
632 //
633 // loop over the pedestal events and check if we have calibration
634 //
635 for (Int_t pixid=0; pixid<fPedestals->GetSize(); pixid++)
636 {
637
638 MCalibrationPix &pix = (*fCalibrations)[pixid];
639
640 //
641 // get the pedestals
642 //
643 const Float_t ped = (*fPedestals)[pixid].GetPedestal() * fNumHiGainSamples;
644 const Float_t prms = (*fPedestals)[pixid].GetPedestalRms() * fSqrtHiGainSamples;
645
646 //
647 // set them in the calibration camera
648 //
649 pix.SetPedestal(ped,prms);
650
651
652 //
653 // Check if the pixel has been excluded from the fits
654 //
655 if (pix.IsExcluded())
656 continue;
657
658 //
659 // perform the Gauss fits to the charges
660 //
661 pix.FitCharge();
662
663 //
664 // Perform the Gauss fits to the arrival times
665 //
666 if (TESTBIT(fFlags,kUseTimes))
667 pix.FitTime();
668
669 }
670
671 if (TESTBIT(fFlags,kUseBlindPixelFit) && fCalibrations->IsBlindPixelMethodValid())
672 {
673
674 if (!fCalibrations->CalcNumPhotInsidePlexiglass())
675 {
676 *fLog << err
677 << "Could not calculate the number of photons from the blind pixel " << endl;
678 *fLog << err
679 << "You can try to calibrate using the MCalibrationCalc::SkipBlindPixelFit()" << endl;
680 fCalibrations->SetBlindPixelMethodValid(kFALSE);
681 }
682 }
683 else
684 *fLog << inf << GetDescriptor() << ": Skipping Blind Pixel Calibration! " << endl;
685
686
687 if (TESTBIT(fFlags,kUsePinDiodeFit) && fCalibrations->IsPINDiodeMethodValid())
688 {
689
690 if (!fCalibrations->CalcNumPhotInsidePlexiglass())
691 {
692 *fLog << err
693 << "Could not calculate the number of photons from the blind pixel " << endl;
694 *fLog << err
695 << "You can try to calibrate using the MCalibrationCalc::SkipPINDiodeFit()" << endl;
696 fCalibrations->SetPINDiodeMethodValid(kFALSE);
697 }
698 }
699 else
700 *fLog << inf << GetDescriptor() << ": Skipping PIN Diode Calibration! " << endl;
701
702 fCalibrations->SetReadyToSave();
703
704 if (GetNumExecutions()==0)
705 return kTRUE;
706
707 *fLog << inf << endl;
708 *fLog << dec << setfill(' ') << fCosmics << " Events presumably cosmics" << endl;
709
710 return kTRUE;
711}
712
Note: See TracBrowser for help on using the repository browser.