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

Last change on this file since 2994 was 2994, checked in by gaug, 21 years ago
*** empty log message ***
File size: 21.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 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 fEvtTime = (MTime*)pList->FindObject("MTime");
217
218 switch (fColor)
219 {
220 case kEBlue:
221 fCalibrations->SetColor(MCalibrationCam::kECBlue);
222 break;
223 case kEGreen:
224 fCalibrations->SetColor(MCalibrationCam::kECGreen);
225 break;
226 case kEUV:
227 fCalibrations->SetColor(MCalibrationCam::kECUV);
228 break;
229 case kECT1:
230 fCalibrations->SetColor(MCalibrationCam::kECCT1);
231 break;
232 default:
233 fCalibrations->SetColor(MCalibrationCam::kECCT1);
234 }
235
236 fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam");
237 if (!fPedestals)
238 {
239 *fLog << err << dbginf << "Cannot find MPedestalCam ... aborting" << endl;
240 return kFALSE;
241 }
242
243
244 fSignals = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam");
245 if (!fSignals)
246 {
247 *fLog << err << dbginf << "Cannot find MExtractedSignalCam ... aborting" << endl;
248 return kFALSE;
249 }
250
251 return kTRUE;
252}
253
254
255// --------------------------------------------------------------------------
256//
257// The ReInit searches for the following input containers:
258// - MRawRunHeader
259//
260Bool_t MCalibrationCalc::ReInit(MParList *pList )
261{
262
263 fRunHeader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
264 if (!fRunHeader)
265 {
266 *fLog << err << dbginf << ": MRawRunHeader not found... aborting." << endl;
267 return kFALSE;
268 }
269
270
271 MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam");
272 if (!cam)
273 {
274 *fLog << err << GetDescriptor() << ": No MGeomCam found... aborting." << endl;
275 return kFALSE;
276 }
277
278 fNumHiGainSamples = fSignals->GetNumUsedHiGainFADCSlices();
279 fNumLoGainSamples = fSignals->GetNumUsedLoGainFADCSlices();
280 fSqrtHiGainSamples = TMath::Sqrt((Float_t)fNumHiGainSamples);
281
282 UInt_t npixels = cam->GetNumPixels();
283
284 for (UInt_t i=0; i<npixels; i++)
285 {
286
287 MCalibrationPix &pix = (*fCalibrations)[i];
288 pix.DefinePixId(i);
289 MHCalibrationPixel *hist = pix.GetHist();
290
291 hist->SetTimeFitRangesHiGain(fSignals->GetFirstUsedSliceHiGain(),
292 fSignals->GetLastUsedSliceHiGain());
293 hist->SetTimeFitRangesLoGain(fSignals->GetFirstUsedSliceLoGain(),
294 fSignals->GetLastUsedSliceLoGain());
295
296 if (!TESTBIT(fFlags,kUseQualityChecks))
297 pix.SetExcludeQualityCheck();
298
299 }
300
301 //
302 // Look for file to exclude pixels from analysis
303 //
304 if (!fExcludedPixelsFile.IsNull())
305 {
306
307 fExcludedPixelsFile = gSystem->ExpandPathName(fExcludedPixelsFile.Data());
308
309 //
310 // Initialize reading the file
311 //
312 ifstream in(fExcludedPixelsFile.Data(),ios::in);
313
314 if (in)
315 {
316 *fLog << inf << "Use excluded pixels from file: '" << fExcludedPixelsFile.Data() << "'" << endl;
317 //
318 // Read the file and count the number of entries
319 //
320 UInt_t pixel = 0;
321
322 while (++fNumExcludedPixels)
323 {
324
325 in >> pixel;
326
327 if (!in.good())
328 break;
329 //
330 // Check for out of range
331 //
332 if (pixel > npixels)
333 {
334 *fLog << warn << "WARNING: To be excluded pixel: " << pixel
335 << " is out of range " << endl;
336 continue;
337 }
338 //
339 // Exclude pixel
340 //
341 MCalibrationPix &pix = (*fCalibrations)[pixel];
342 pix.SetExcluded();
343
344 *fLog << GetDescriptor() << inf << ": Exclude Pixel: " << pixel << endl;
345
346 }
347
348 if (--fNumExcludedPixels == 0)
349 *fLog << warn << "WARNING: File '" << fExcludedPixelsFile.Data()
350 << "'" << " is empty " << endl;
351 else
352 fCalibrations->SetNumPixelsExcluded(fNumExcludedPixels);
353
354 }
355 else
356 *fLog << warn << dbginf << "Cannot open file '" << fExcludedPixelsFile.Data() << "'" << endl;
357 }
358
359 return kTRUE;
360}
361
362
363// --------------------------------------------------------------------------
364//
365// Calculate the integral of the FADC time slices and store them as a new
366// pixel in the MCerPhotEvt container.
367//
368Int_t MCalibrationCalc::Process()
369{
370
371 //
372 // Initialize pointers to blind pixel, PIN Diode and individual pixels
373 //
374 MCalibrationBlindPix &blindpixel = *(fCalibrations->GetBlindPixel());
375 MCalibrationPINDiode &pindiode = *(fCalibrations->GetPINDiode());
376
377 MRawEvtPixelIter pixel(fRawEvt);
378
379 //
380 // Perform cosmics cut
381 //
382 if (TESTBIT(fFlags,kUseCosmicsRejection))
383 {
384
385 Int_t cosmicpix = 0;
386
387 //
388 // Create a first loop to sort out the cosmics ...
389 //
390 // This is a very primitive check for the number of cosmicpixs
391 // The cut will be applied in the fit, but for the blind pixel,
392 // we need to remove this event
393 //
394 // FIXME: In the future need a much more sophisticated one!!!
395 //
396
397 while (pixel.Next())
398 {
399
400 const UInt_t pixid = pixel.GetPixelId();
401
402 MExtractedSignalPix &sig = (*fSignals)[pixid];
403 MPedestalPix &ped = (*fPedestals)[pixid];
404 Float_t pedrms = ped.GetPedestalRms()*fSqrtHiGainSamples;
405 Float_t sumhi = sig.GetExtractedSignalHiGain();
406
407 //
408 // We consider a pixel as presumably due to cosmics
409 // if its sum of FADC slices is lower than 3 pedestal RMS
410 //
411 if (sumhi < 3.*pedrms )
412 cosmicpix++;
413 }
414
415 //
416 // If the camera contains more than 230
417 // (this is the number of outer pixels plus about 50 inner ones)
418 // presumed pixels due to cosmics, then the event is discarted.
419 // This procedure is more or less equivalent to keeping only events
420 // with at least 350 pixels with high signals.
421 //
422 if (cosmicpix > 230.)
423 {
424 fCosmics++;
425 return kCONTINUE;
426 }
427
428 pixel.Reset();
429 }
430
431 fEvents++;
432
433 Int_t overflow = 0;
434
435 //
436 // Create a (second) loop to do fill the calibration histograms
437 //
438
439 while (pixel.Next())
440 {
441
442 const UInt_t pixid = pixel.GetPixelId();
443
444 MCalibrationPix &pix = (*fCalibrations)[pixid];
445
446 if (pix.IsExcluded())
447 continue;
448
449 MExtractedSignalPix &sig = (*fSignals)[pixid];
450
451 const Float_t sumhi = sig.GetExtractedSignalHiGain();
452 const Float_t sumlo = sig.GetExtractedSignalLoGain();
453
454 Double_t mtime = 0.;
455
456 if (TESTBIT(fFlags,kUseTimes))
457 {
458 //
459 // first, have a look in MArrivalTime,
460 // otherwise search the position of maximum bin
461 // in MRawEvtData
462 //
463 if (fArrivalTime)
464 mtime = (*fArrivalTime)[pixid];
465 else
466 if (sig.IsLoGainUsed())
467 fRawEvt->GetPixelContent(mtime, pixid, NULL, 4);
468 else
469 fRawEvt->GetPixelContent(mtime, pixid, NULL, 3);
470 }
471
472 switch(pixid)
473 {
474
475 case fBlindPixelId:
476
477 if (!blindpixel.FillCharge(sumhi))
478 *fLog << warn <<
479 "Overflow or Underflow occurred filling Blind Pixel sum = " << sumhi << endl;
480
481 if (TESTBIT(fFlags,kUseTimes))
482 {
483 if (!blindpixel.FillTime(mtime))
484 *fLog << warn <<
485 "Overflow or Underflow occurred filling Blind Pixel time = " << mtime << endl;
486 }
487
488 if (!blindpixel.FillRChargevsTime(sumhi,fEvents))
489 *fLog << warn <<
490 "Overflow or Underflow occurred filling Blind Pixel eventnr = " << fEvents << endl;
491 break;
492
493 case fPINDiodeId:
494
495 if (!pindiode.FillCharge(sumhi))
496 *fLog << warn <<
497 "Overflow or Underflow occurred filling PINDiode: sum = " << sumhi << endl;
498
499 if (TESTBIT(fFlags,kUseTimes))
500 {
501 if (!pindiode.FillTime(mtime))
502 *fLog << warn <<
503 "Overflow or Underflow occurred filling PINDiode: time = " << mtime << endl;
504 }
505
506 if (!pindiode.FillRChargevsTime(sumhi,fEvents))
507 *fLog << warn <<
508 "Overflow or Underflow occurred filling PINDiode: eventnr = " << fEvents << endl;
509 break;
510
511 default:
512
513 pix.FillChargesInGraph(sumhi,sumlo);
514
515 if (!pix.FillRChargevsTimeLoGain(sumlo,fEvents))
516 overflow++;
517
518 if (!pix.FillRChargevsTimeHiGain(sumhi,fEvents))
519 overflow++;
520
521 if (sig.IsLoGainUsed())
522 {
523
524 if (!pix.FillChargeLoGain(sumlo))
525 *fLog << warn << "Could not fill Lo Gain Charge of pixel: " << pixid
526 << " signal = " << sumlo << endl;
527
528 if (TESTBIT(fFlags,kUseTimes))
529 {
530 if (!pix.FillTimeLoGain(mtime))
531 *fLog << warn << "Could not fill Lo Gain Time of pixel: "
532 << pixid << " time = " << mtime << endl;
533 }
534 } /* if (sig.IsLoGainUsed()) */
535 else
536 {
537 if (!pix.FillChargeHiGain(sumhi))
538 *fLog << warn << "Could not fill Hi Gain Charge of pixel: " << pixid
539 << " signal = " << sumhi << endl;
540
541 if (TESTBIT(fFlags,kUseTimes))
542 {
543 if (!pix.FillTimeHiGain(mtime))
544 *fLog << warn << "Could not fill Hi Gain Time of pixel: "
545 << pixid << " time = " << mtime << endl;
546 }
547 } /* else (sig.IsLoGainUsed()) */
548 break;
549
550 } /* switch(pixid) */
551
552 } /* while (pixel.Next()) */
553
554 if (overflow)
555 *fLog << warn << "Overflow occurred filling Charges vs. EvtNr " << overflow << " times" << endl;
556
557 return kTRUE;
558}
559
560Int_t MCalibrationCalc::PostProcess()
561{
562
563 *fLog << inf << endl;
564
565 if (fEvents == 0)
566 {
567 *fLog << err << GetDescriptor()
568 << ": This run contains only cosmics or pedestals, "
569 << "cannot find events with more than 350 illuminated pixels. " << endl;
570 return kFALSE;
571 }
572
573 if (fEvents < fCosmics)
574 *fLog << warn << GetDescriptor()
575 << ": WARNING: Run contains more cosmics or pedestals than calibration events " << endl;
576
577
578 *fLog << inf << GetDescriptor() << ": Cut Histogram Edges" << endl;
579
580 //
581 // Cut edges to make fits and viewing of the hists easier
582 //
583 fCalibrations->CutEdges();
584
585 //
586 // Fit the blind pixel
587 //
588 if (TESTBIT(fFlags,kUseBlindPixelFit))
589 {
590 //
591 // Get pointer to blind pixel
592 //
593 MCalibrationBlindPix &blindpixel = *(fCalibrations->GetBlindPixel());
594
595 *fLog << inf << GetDescriptor() << ": Fitting the Blind Pixel" << endl;
596
597 //
598 // retrieve mean and sigma of the blind pixel pedestal,
599 // so that we can use it for the fit
600 //
601 if (fPedestals->IsUseHists())
602 {
603 //
604 // retrieve the pedestal pix of the blind pixel
605 //
606 MPedestalPix &ped = (*fPedestals)[fBlindPixelId];
607 //
608 // retrieve the blind pixel histogram container
609 //
610 MHCalibrationBlindPixel *hist = blindpixel.GetHist();
611 //
612 // Set the corresponding values
613 //
614 const Float_t peddiff = ped.GetMean()
615 - ped.GetPedestal()*fSignals->GetNumUsedFADCSlices();
616 const Float_t pederr = ped.GetMeanErr();
617 const Float_t pedsigma = ped.GetSigma();
618 const Float_t pedsigmaerr = ped.GetSigmaErr();
619
620 hist->SetMeanPedestal(peddiff);
621 hist->SetMeanPedestalErr(pederr);
622 hist->SetSigmaPedestal(pedsigma);
623 hist->SetSigmaPedestalErr(pedsigmaerr);
624 }
625
626 if (!blindpixel.FitCharge())
627 {
628 *fLog << warn << "Could not fit the blind pixel! " << endl;
629 *fLog << warn << "Setting bit kBlindPixelMethodValid to FALSE in MCalibrationCam" << endl;
630 fCalibrations->SetBlindPixelMethodValid(kFALSE);
631 }
632
633 fCalibrations->SetBlindPixelMethodValid(kTRUE);
634 blindpixel.DrawClone();
635 }
636 else
637 *fLog << inf << GetDescriptor() << ": Skipping Blind Pixel Fit " << endl;
638
639
640 *fLog << inf << GetDescriptor() << ": Fitting the Normal Pixels" << endl;
641
642 //
643 // loop over the pedestal events and check if we have calibration
644 //
645 for (Int_t pixid=0; pixid<fPedestals->GetSize(); pixid++)
646 {
647
648 MCalibrationPix &pix = (*fCalibrations)[pixid];
649
650 //
651 // get the pedestals
652 //
653 const Float_t ped = (*fPedestals)[pixid].GetPedestal() * fNumHiGainSamples;
654 const Float_t prms = (*fPedestals)[pixid].GetPedestalRms() * fSqrtHiGainSamples;
655
656 //
657 // set them in the calibration camera
658 //
659 pix.SetPedestal(ped,prms);
660
661
662 //
663 // Check if the pixel has been excluded from the fits
664 //
665 if (pix.IsExcluded())
666 continue;
667
668 //
669 // perform the Gauss fits to the charges
670 //
671 pix.FitCharge();
672
673 //
674 // Perform the Gauss fits to the arrival times
675 //
676 if (TESTBIT(fFlags,kUseTimes))
677 pix.FitTime();
678
679 }
680
681 if (TESTBIT(fFlags,kUseBlindPixelFit) && fCalibrations->IsBlindPixelMethodValid())
682 {
683
684 if (!fCalibrations->CalcNumPhotInsidePlexiglass())
685 {
686 *fLog << err
687 << "Could not calculate the number of photons from the blind pixel " << endl;
688 *fLog << err
689 << "You can try to calibrate using the MCalibrationCalc::SkipBlindPixelFit()" << endl;
690 fCalibrations->SetBlindPixelMethodValid(kFALSE);
691 }
692 }
693 else
694 *fLog << inf << GetDescriptor() << ": Skipping Blind Pixel Calibration! " << endl;
695
696
697 if (TESTBIT(fFlags,kUsePinDiodeFit) && fCalibrations->IsPINDiodeMethodValid())
698 {
699
700 if (!fCalibrations->CalcNumPhotInsidePlexiglass())
701 {
702 *fLog << err
703 << "Could not calculate the number of photons from the blind pixel " << endl;
704 *fLog << err
705 << "You can try to calibrate using the MCalibrationCalc::SkipPINDiodeFit()" << endl;
706 fCalibrations->SetPINDiodeMethodValid(kFALSE);
707 }
708 }
709 else
710 *fLog << inf << GetDescriptor() << ": Skipping PIN Diode Calibration! " << endl;
711
712 fCalibrations->SetReadyToSave();
713
714 if (GetNumExecutions()==0)
715 return kTRUE;
716
717 *fLog << inf << endl;
718 *fLog << dec << setfill(' ') << fCosmics << " Events presumably cosmics" << endl;
719
720 return kTRUE;
721}
722
Note: See TracBrowser for help on using the repository browser.