source: trunk/MagicSoft/Mars/mcalib/MCalibrateData.cc@ 8462

Last change on this file since 8462 was 8452, checked in by tbretz, 19 years ago
*** empty log message ***
File size: 31.0 KB
Line 
1/* ======================================================================== *\
2! $Name: not supported by cvs2svn $:$Id: MCalibrateData.cc,v 1.68 2007-04-27 10:04:46 tbretz Exp $
3! --------------------------------------------------------------------------
4!
5! *
6! * This file is part of MARS, the MAGIC Analysis and Reconstruction
7! * Software. It is distributed to you in the hope that it can be a useful
8! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
9! * It is distributed WITHOUT ANY WARRANTY.
10! *
11! * Permission to use, copy, modify and distribute this software and its
12! * documentation for any purpose is hereby granted without fee,
13! * provided that the above copyright notice appear in all copies and
14! * that both that copyright notice and this permission notice appear
15! * in supporting documentation. It is provided "as is" without express
16! * or implied warranty.
17! *
18!
19!
20! Author(s): Javier Lopez 12/2003 <mailto:jlopez@ifae.es>
21! Author(s): Javier Rico 01/2004 <mailto:jrico@ifae.es>
22! Author(s): Wolfgang Wittek 02/2004 <mailto:wittek@mppmu.mpg.de>
23! Author(s): Markus Gaug 04/2004 <mailto:markus@ifae.es>
24! Author(s): Hendrik Bartko 08/2004 <mailto:hbartko@mppmu.mpg.de>
25! Author(s): Thomas Bretz 08/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
26!
27! Copyright: MAGIC Software Development, 2000-2006
28!
29!
30\* ======================================================================== */
31
32///////////////////////////////////////////////////////////////////////////////////
33//
34// MCalibrateData
35//
36// This task takes the integrated charge from MExtractedSignal and applies
37// the calibration constants from MCalibrationCam to convert the summed FADC
38// slices into photons. The number of photons obtained is stored in MSignalCam.
39// Optionally, the calibration of pedestals from an MPedestalCam container into
40// an MPedPhotCam container can be chosen with the member functions
41// SetPedestalType(). Default is 'kRun', i.e. calibration of pedestals from a
42// dedicated pedestal run.
43// In case, the chosen pedestal type is kRun or kEvent, in ReInit() the MPedPhotCam
44// container is filled using the information from MPedestalCam, MExtractedSignalCam,
45// MCalibrationChargeCam and MCalibrationQECam
46//
47// Selection of different calibration methods is allowed through the
48// SetCalibrationMode() member function (default: kFfactor)
49//
50// The calibration modes which exclude non-valid pixels are the following:
51//
52// kFfactor: calibrates using the F-Factor method
53// kBlindpixel: calibrates using the BlindPixel method
54// kBlindpixel: calibrates using the BlindPixel method
55// kFlatCharge: perform a charge flat-flatfielding. Outer pixels are area-corrected.
56// kDummy: calibrates with fixed conversion factors of 1 and errors of 0.
57//
58// The calibration modes which include all pixels regardless of their validity is:
59//
60// kNone: calibrates with fixed conversion factors of 1 and errors of 0.
61//
62// Use the kDummy and kNone methods ONLY FOR DEBUGGING!
63//
64//
65// This class can calibrate data and/or pedestals. To switch off calibration of data
66// set the Calibration Mode to kSkip. To switch on pedestal calibration call either
67// SetPedestalFlag(MCalibrateData::kRun) (calibration is done once in ReInit)
68// SetPedestalFlag(MCalibrateData::kEvent) (calibration is done for each event)
69//
70// By calling AddPedestal() you can control the name of the
71// MPedestalCam and/or MPedPhotCam container which is used.
72//
73// Assume you want to calibrate "MPedestalCam" once and "MPedestalFromLoGain" [MPedestalCam]
74// event-by-event, so:
75// MCalibrateData cal1;
76// cal1.SetCalibrationMode(MCalibrateData::kSkip);
77// cal1.SetPedestalFlag(MCalibrateData::kRun);
78// MCalibrateData cal2;
79// cal2.SetCalibrationMode(MCalibrateData::kSkip);
80// cal2.AddPedestal("FromLoGain");
81// cal2.SetPedestalFlag(MCalibrateData::kEvent);
82//
83//
84// Input Containers:
85// [MPedestalCam]
86// [MExtractedSignalCam]
87// [MCalibrationChargeCam]
88// [MCalibrationQECam]
89// MBadPixelsCam
90//
91// Output Containers:
92// [MPedPhotCam]
93// [MSignalCam]
94//
95// See also: MJCalibration, MJPedestal, MJExtractSignal, MJExtractCalibTest
96//
97//////////////////////////////////////////////////////////////////////////////
98#include "MCalibrateData.h"
99
100#include <fstream>
101
102#include <TEnv.h>
103
104#include "MLog.h"
105#include "MLogManip.h"
106
107#include "MParList.h"
108#include "MH.h"
109
110#include "MGeomCam.h"
111#include "MRawRunHeader.h"
112
113#include "MPedestalCam.h"
114#include "MPedestalPix.h"
115
116#include "MCalibrationChargeCam.h"
117#include "MCalibrationChargePix.h"
118
119#include "MCalibrationQECam.h"
120#include "MCalibrationQEPix.h"
121
122#include "MCalibConstCam.h"
123#include "MCalibConstPix.h"
124
125#include "MExtractedSignalCam.h"
126#include "MExtractedSignalPix.h"
127
128#include "MPedPhotCam.h"
129#include "MPedPhotPix.h"
130
131#include "MBadPixelsCam.h"
132#include "MBadPixelsPix.h"
133
134#include "MSignalCam.h"
135
136ClassImp(MCalibrateData);
137
138using namespace std;
139
140const Float_t MCalibrateData::gkCalibConvMinLimit = 0.01;
141const Float_t MCalibrateData::gkCalibConvMaxLimit = 5.;
142
143const MCalibrateData::CalibrationMode_t MCalibrateData::gkDefault = kFfactor;
144
145// --------------------------------------------------------------------------
146//
147// Default constructor.
148//
149// Sets all pointers to NULL
150//
151// Initializes:
152// - fCalibrationMode to kDefault
153// - fPedestalFlag to kNo
154//
155MCalibrateData::MCalibrateData(CalibrationMode_t calmode,const char *name, const char *title)
156 : fGeomCam(NULL), fBadPixels(NULL), fCalibrations(NULL),
157 fQEs(NULL), fSignals(NULL), fCerPhotEvt(NULL), fCalibConstCam(NULL),
158 fPedestalFlag(kNo), fSignalType(kPhot), fRenormFactor(1.), fScaleFactor(1.)
159{
160
161 fName = name ? name : "MCalibrateData";
162 fTitle = title ? title : "Task to calculate the number of photons in one event";
163
164 SetCalibrationMode(calmode);
165
166 SetCalibConvMinLimit();
167 SetCalibConvMaxLimit();
168
169 fNamesPedestal.SetOwner();
170}
171
172void MCalibrateData::AddPedestal(const char *name)
173{
174 TString ped(name);
175 TString pho(name);
176 ped.Prepend("MPedestal");
177 pho.Prepend("MPedPhot");
178
179 fNamesPedestal.Add(new TNamed(ped, pho));
180}
181
182void MCalibrateData::AddPedestal(const char *pedestal, const char *pedphot)
183{
184 fNamesPedestal.Add(new TNamed(pedestal, pedphot));
185}
186
187// --------------------------------------------------------------------------
188//
189// The PreProcess searches for the following input containers:
190//
191// - MGeomCam
192// - MPedestalCam
193// - MCalibrationChargeCam
194// - MCalibrationQECam
195// - MExtractedSignalCam
196// - MBadPixelsCam
197//
198// The following output containers are also searched and created if
199// they were not found:
200//
201// - MPedPhotCam
202// - MSignalCam
203//
204Int_t MCalibrateData::PreProcess(MParList *pList)
205{
206 // input containers
207
208 fBadPixels = (MBadPixelsCam*)pList->FindObject(AddSerialNumber("MBadPixelsCam"));
209 if (!fBadPixels)
210 {
211 *fLog << err << AddSerialNumber("MBadPixelsCam") << " not found ... aborting" << endl;
212 return kFALSE;
213 }
214
215 fSignals = 0;
216 fCerPhotEvt = 0;
217 if (fCalibrationMode>kSkip)
218 {
219 fSignals = (MExtractedSignalCam*)pList->FindObject(AddSerialNumber("MExtractedSignalCam"));
220 if (!fSignals)
221 {
222 *fLog << err << AddSerialNumber("MExtractedSignalCam") << " not found ... aborting" << endl;
223 return kFALSE;
224 }
225
226 fCerPhotEvt = (MSignalCam*)pList->FindCreateObj(AddSerialNumber("MSignalCam"));
227 if (!fCerPhotEvt)
228 return kFALSE;
229
230 fCalibConstCam = (MCalibConstCam*)pList->FindCreateObj(AddSerialNumber("MCalibConstCam"));
231 if (!fCalibConstCam)
232 return kFALSE;
233 }
234
235 fCalibrations = 0;
236 fQEs = 0;
237 if (fCalibrationMode>kNone)
238 {
239 fCalibrations = (MCalibrationChargeCam*)pList->FindObject(AddSerialNumber("MCalibrationChargeCam"));
240 if (!fCalibrations)
241 {
242 *fLog << err << AddSerialNumber("MCalibrationChargeCam") << " not found ... aborting." << endl;
243 return kFALSE;
244 }
245
246 fQEs = (MCalibrationQECam*)pList->FindObject(AddSerialNumber("MCalibrationQECam"));
247 if (!fQEs)
248 {
249 *fLog << err << AddSerialNumber("MCalibrationQECam") << " not found ... aborting." << endl;
250 return kFALSE;
251 }
252
253 }
254
255 if (fNamesPedestal.GetSize()>0 && fPedestalFlag==kNo)
256 {
257 *fLog << warn << "Pedestal list contains entries, but mode is set to kNo... setting to kEvent." << endl;
258 fPedestalFlag = kEvent;
259 }
260
261 if (fPedestalFlag!=kNo)
262 {
263 if (fNamesPedestal.GetSize()==0)
264 {
265 *fLog << inf << "No container names specified... using default: MPedestalCam and MPedPhotCam." << endl;
266 AddPedestal();
267 }
268
269 fPedestalCams.Clear();
270 fPedPhotCams.Clear();
271
272 TIter Next(&fNamesPedestal);
273 TObject *o=0;
274 while ((o=Next()))
275 {
276 TObject *pedcam = pList->FindObject(AddSerialNumber(o->GetName()), "MPedestalCam");
277 if (!pedcam)
278 {
279 *fLog << err << AddSerialNumber(o->GetName()) << " [MPedestalCam] not found ... aborting" << endl;
280 return kFALSE;
281 }
282 TObject *pedphot = pList->FindCreateObj("MPedPhotCam", AddSerialNumber(o->GetTitle()));
283 if (!pedphot)
284 return kFALSE;
285
286 fPedestalCams.Add(pedcam);
287 fPedPhotCams.Add(pedphot);
288 }
289 }
290
291 switch (fSignalType)
292 {
293 case kPhe:
294 //
295 // Average QE for C-photons, for pixels in the inner pixel region ("area 0"),
296 // used later to convert from C-photons into "equivalent phes":
297 //
298
299 // Average areas not yet initialized? use default value
300 fRenormFactor = MCalibrationQEPix::gkDefaultAverageQE;
301 if (fQEs->GetAverageAreas() > 0)
302 {
303 MCalibrationQEPix& avqepix = (MCalibrationQEPix&)(fQEs->GetAverageArea(0));
304 fRenormFactor = avqepix.GetAverageQE();
305 }
306 break;
307
308 case kPhot:
309 fRenormFactor = 1.;
310 break;
311 }
312
313 fCalibConsts.Reset();
314 fCalibFFactors.Reset();
315 fHiLoConv.Reset();
316 fHiLoConvErr.Reset();
317
318 return kTRUE;
319}
320
321// --------------------------------------------------------------------------
322//
323// The ReInit searches for the following input containers:
324//
325// - MGeomCam
326//
327// Check for validity of the selected calibration method, switch to a
328// different one in case of need
329//
330// Fill the MPedPhotCam container using the information from MPedestalCam,
331// MExtractedSignalCam and MCalibrationCam
332//
333Bool_t MCalibrateData::ReInit(MParList *pList)
334{
335 MRawRunHeader *header = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
336 if (!header)
337 {
338 *fLog << err << "MRawRunHeader not found... abort." << endl;
339 return kFALSE;
340 }
341
342 fGeomCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
343 if (!fGeomCam)
344 {
345 *fLog << err << "No MGeomCam found... aborting." << endl;
346 return kFALSE;
347 }
348
349 // Sizes might have changed
350 if (fPedestalFlag!=kNo)
351 {
352 TIter Next(&fPedestalCams);
353 MPedestalCam *cam=0;
354 while ((cam=(MPedestalCam*)Next()))
355 if ((Int_t)cam->GetSize() != fSignals->GetSize())
356 {
357 *fLog << err << "Size mismatch of " << cam->GetDescriptor() << " and MCalibrationCam... abort." << endl;
358 return kFALSE;
359 }
360 }
361
362 if(fCalibrationMode == kBlindPixel && !fQEs->IsBlindPixelMethodValid())
363 {
364 *fLog << warn << "Blind pixel calibration method not valid, switching to F-factor method" << endl;
365 fCalibrationMode = kFfactor;
366 }
367
368 if(fCalibrationMode == kPinDiode && !fQEs->IsPINDiodeMethodValid())
369 {
370 *fLog << warn << "PIN diode calibration method not valid, switching to F-factor method" << endl;
371 fCalibrationMode = kFfactor;
372 }
373
374 if(fCalibrationMode == kCombined && !fQEs->IsCombinedMethodValid())
375 {
376 *fLog << warn << "Combined calibration method not valid, switching to F-factor method" << endl;
377 fCalibrationMode = kFfactor;
378 }
379
380 //
381 // output information or warnings:
382 //
383 switch(fCalibrationMode)
384 {
385 case kBlindPixel:
386 break;
387 case kFfactor:
388 break;
389 case kPinDiode:
390 *fLog << err << "PIN Diode Calibration mode not yet available!" << endl;
391 return kFALSE;
392 break;
393 case kCombined:
394 *fLog << err << "Combined Calibration mode not yet available!" << endl;
395 return kFALSE;
396 break;
397 case kFlatCharge:
398 *fLog << warn << "WARNING - Flat-fielding charges - only for muon calibration!" << endl;
399 break;
400 case kDummy:
401 *fLog << warn << "WARNING - Dummy calibration, no calibration applied!" << endl;
402 break;
403 case kNone:
404 *fLog << warn << "WARNING - No calibration applied!" << endl;
405 break;
406 default:
407 *fLog << warn << "WARNING - Calibration mode value (" << fCalibrationMode << ") not known" << endl;
408 return kFALSE;
409 }
410
411 //
412 // output information or warnings:
413 //
414 switch (fSignalType)
415 {
416 case kPhe:
417 *fLog << inf << "Calibrating in units of equivalent (outer/inner=4) photo-electrons." << endl;
418 break;
419 case kPhot:
420 *fLog << inf << "Calibrating in units of photons." << endl;
421 break;
422 }
423
424 if (header->IsMonteCarloRun())
425 {
426 *fLog << inf << "Additional scale factor forced to: 1 (MonteCarloRun)" << endl;
427 fScaleFactor = 1;
428 }
429 else
430 {
431 if (!fFileNameScale.IsNull())
432 {
433 if (gSystem->AccessPathName(fFileNameScale, kFileExists))
434 {
435 *fLog << err << "ERROR - Configuration file '" << fFileNameScale << "' doesn't exist... abort." << endl;
436 return kFALSE;
437 }
438
439 const Int_t p = header->GetRunStart().GetMagicPeriod();
440 const TString per = Form("%2d", p);
441
442 TEnv rc(fFileNameScale);
443 const Double_t scale = rc.GetValue(per, -1.);
444
445 if (scale<=0)
446 {
447 *fLog << err << "ERROR - No valid entry for scale factor found for period " << p << " in " << fFileNameScale << "... abort." << endl;
448 return kFALSE;
449 }
450
451 *fLog << inf << "New additional scale factor for period " << p << ": " << scale << endl;
452 fScaleFactor = scale;
453 }
454 else
455 *fLog << inf << "Additional scale factor set to: " << fScaleFactor << endl;
456 }
457
458 const Int_t npixels = fGeomCam->GetNumPixels();
459
460 if (fCalibrationMode > kNone)
461 {
462 if (fCalibrations->GetSize() != npixels)
463 {
464 *fLog << "Size mismatch between MGeomCam and MCalibrationChargeCam... abort!" << endl;
465 return kFALSE;
466 }
467
468 if (fBadPixels->GetSize() != npixels)
469 {
470 *fLog << "Size mismatch between MGeomCam and MBadPixelsCam... abort!" << endl;
471 return kFALSE;
472 }
473 }
474
475 fCalibConsts .Set(npixels);
476 fCalibFFactors.Set(npixels);
477 fHiLoConv .Set(npixels);
478 fHiLoConvErr .Set(npixels);
479
480 if (!UpdateConversionFactors())
481 return kFALSE;
482
483 if (TestPedestalFlag(kRun))
484 Calibrate(kFALSE, kTRUE);
485
486 return kTRUE;
487}
488
489// --------------------------------------------------------------------------
490//
491// Update the conversion factors and F-Factors from MCalibrationCams into
492// the arrays. Later, the here pre-calcualted conversion factors get simply
493// copied from memory.
494//
495// This function can be called from outside in case that the MCalibrationCams
496// have been updated...
497//
498Bool_t MCalibrateData::UpdateConversionFactors( const MCalibrationChargeCam *updatecam)
499{
500
501 *fLog << inf << GetDescriptor() << ": Updating Conversion Factors... " << endl;
502
503 fCalibConsts.Reset();
504 fCalibFFactors.Reset();
505 fHiLoConv.Reset();
506 fHiLoConvErr.Reset();
507
508 //
509 // For the moment, we use only a dummy zenith for the calibration:
510 //
511 UInt_t skip = 0;
512
513 for (UInt_t pixidx=0; pixidx<fGeomCam->GetNumPixels(); pixidx++)
514 {
515
516 Float_t hiloconv = 1.;
517 Float_t hiloconverr = 0.;
518 Float_t calibConv = 1.;
519 Float_t calibConvVar = 0.;
520 Float_t calibFFactor = 0.;
521
522 Float_t calibQE = 1.;
523 Float_t calibQEVar = 0.;
524
525 Float_t calibUpdate = 1.;
526
527 if(fCalibrationMode!=kNone)
528 {
529 if ((*fBadPixels)[pixidx].IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
530 {
531 skip++;
532 continue; // calibConv will remain 0
533 }
534
535 const MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCalibrations)[pixidx];
536 const MCalibrationChargePix &avpix = (MCalibrationChargePix&)fCalibrations->GetAverageArea(0);
537
538 hiloconv = pix.GetConversionHiLo ();
539 hiloconverr= pix.GetConversionHiLoSigma();
540
541 calibConv = pix.GetMeanConvFADC2Phe();
542 calibConvVar = pix.GetMeanConvFADC2PheVar();
543 calibFFactor = pix.GetMeanFFactorFADC2Phot();
544
545 const MCalibrationQEPix &qe = (MCalibrationQEPix&)(*fQEs)[pixidx];
546
547 if (updatecam)
548 {
549 const MCalibrationChargePix &upix = (MCalibrationChargePix&)(*updatecam)[pixidx];
550
551 //
552 // Correct for the possible change in amplification of the individual pixels chain
553 //
554 const Float_t pixmean = upix.GetConvertedMean();
555 calibUpdate = pixmean==0 ? 1 : pix.GetConvertedMean()/pixmean;
556
557 //
558 // Correct for global shifts in light emission
559 //
560 const MCalibrationChargePix &ugpix = (MCalibrationChargePix&)updatecam->GetAverageArea(0);
561
562 const Float_t globmean = avpix.GetConvertedMean();
563 calibUpdate = globmean==0 ? 1 : ugpix.GetConvertedMean()/globmean;
564
565 MBadPixelsPix &ubad = (MBadPixelsPix&)updatecam->GetAverageBadArea(0);
566 if (ubad.IsUncalibrated(MBadPixelsPix::kChargeIsPedestal))
567 {
568 *fLog << warn << GetDescriptor() << ": Mean charge in inner pixels is smaller than 3 ped. RMS." << endl;
569 *fLog << "Maybe calibration pulses have been lost!" << endl;
570 calibUpdate = 1.;
571
572 }
573 }
574
575 switch(fCalibrationMode)
576 {
577 case kFlatCharge:
578 {
579 calibConv = avpix.GetConvertedMean()
580 / (pix.GetConvertedMean() * fGeomCam->GetPixRatio(pixidx));
581 calibConvVar = (avpix.GetMeanRelVar() + pix.GetMeanRelVar()) * calibConv * calibConv;
582 if (pix.IsFFactorMethodValid())
583 {
584 const Float_t convmin1 = qe.GetQECascadesFFactor()/pix.GetMeanConvFADC2Phe();
585 if (convmin1 > 0)
586 calibFFactor *= TMath::Sqrt(convmin1);
587 else
588 calibFFactor = -1.;
589 }
590 break;
591 }
592 case kBlindPixel:
593 if (!qe.IsBlindPixelMethodValid())
594 {
595 skip++;
596 continue;
597 }
598 calibQE = qe.GetQECascadesBlindPixel();
599 calibQEVar = qe.GetQECascadesBlindPixelVar();
600 break;
601
602 case kPinDiode:
603 if (!qe.IsPINDiodeMethodValid())
604 {
605 skip++;
606 continue;
607 }
608 calibQE = qe.GetQECascadesPINDiode();
609 calibQEVar = qe.GetQECascadesPINDiodeVar();
610 break;
611
612 case kFfactor:
613 if (!pix.IsFFactorMethodValid())
614 {
615 skip++;
616 continue;
617 }
618 calibQE = qe.GetQECascadesFFactor();
619 calibQEVar = qe.GetQECascadesFFactorVar();
620 break;
621
622 case kCombined:
623 if (!qe.IsCombinedMethodValid())
624 {
625 skip++;
626 continue;
627 }
628 calibQE = qe.GetQECascadesCombined();
629 calibQEVar = qe.GetQECascadesCombinedVar();
630 break;
631
632 case kDummy:
633 hiloconv = 1.;
634 hiloconverr = 0.;
635 calibUpdate = 1.;
636 break;
637
638 default:
639 break;
640 } /* switch calibration mode */
641 } /* if(fCalibrationMode!=kNone) */
642 else
643 {
644 calibConv = 1./fGeomCam->GetPixRatio(pixidx);
645 }
646
647 calibConv /= calibQE;
648
649 if (calibConv > 0.00001 && calibQE > 0.00001)
650 {
651 calibConvVar = calibConvVar/(calibConv*calibConv) + calibQEVar/(calibQE*calibQE);
652 calibConvVar *= (calibConv*calibConv);
653 }
654
655 calibConv *= fRenormFactor*fScaleFactor * calibUpdate;
656 calibFFactor *= TMath::Sqrt(fRenormFactor*fScaleFactor);
657
658 fHiLoConv [pixidx] = hiloconv;
659 fHiLoConvErr [pixidx] = hiloconverr;
660 fCalibConsts [pixidx] = calibConv;
661 fCalibFFactors[pixidx] = calibFFactor;
662
663 if (calibConv < fCalibConvMinLimit || calibConv > fCalibConvMaxLimit)
664 {
665 (*fBadPixels)[pixidx].SetUnsuitable(MBadPixelsPix::kUnsuitableRun);
666 calibConv = -1.;
667 calibFFactor = -1.;
668 *fLog << warn << GetDescriptor() << ": WARNING - ";
669 *fLog << "Conversion factor " << calibConv << " of Pixel " << pixidx << " out of range ]";
670 *fLog << fCalibConvMinLimit << "," << fCalibConvMaxLimit << "[... set to 0. " << endl;
671 }
672
673 MCalibConstPix &cpix = (*fCalibConstCam)[pixidx];
674
675 cpix.SetCalibConst(calibConv);
676 cpix.SetCalibFFactor(calibFFactor);
677
678 } /* for (Int_t pixidx=0; pixidx<fGeomCam->GetNumPixels(); pixidx++) */
679
680 fCalibConstCam->SetReadyToSave();
681
682 if (skip>fGeomCam->GetNumPixels()*0.9)
683 {
684 *fLog << err << GetDescriptor() << ": ERROR - ";
685 *fLog << "GetConversionFactor has skipped more than 90% of the pixels... abort." << endl;
686 return kFALSE;
687 }
688
689 return kTRUE;
690}
691
692
693// --------------------------------------------------------------------------
694//
695// Apply the conversion factors and F-Factors from the arrays to the data.
696//
697// The flags 'data' and 'pedestal' decide whether the signal and/or the pedetals
698// shall be calibrated, respectively.
699//
700Int_t MCalibrateData::Calibrate(Bool_t data, Bool_t pedestal) const
701{
702 if (!data && !pedestal)
703 return kTRUE;
704
705 const UInt_t npix = fSignals->GetSize();
706
707 // The number of used slices are just a mean value
708 // the real number might change from event to event.
709 // (up to 50%!)
710 const Float_t slices = fSignals->GetNumUsedHiGainFADCSlices();
711 const Float_t sqrtslices = TMath::Sqrt(slices);
712
713 Int_t numsatlo=0;
714 Int_t numsathi=0;
715
716 for (UInt_t pixidx=0; pixidx<npix; pixidx++)
717 {
718 MBadPixelsPix &bpix = (*fBadPixels)[pixidx];
719
720 if (data)
721 {
722 const MExtractedSignalPix &sig = (*fSignals)[pixidx];
723
724 Float_t signal = 0.;
725 Float_t signalErr = 0.;
726 Bool_t ok = kFALSE;
727
728 // If hi-gain is not saturated and has successfully been
729 // extracted use the hi-gain arrival time
730 if (!sig.IsHiGainSaturated() && sig.IsHiGainValid())
731 {
732 signal = sig.GetExtractedSignalHiGain();
733 ok = kTRUE;
734 }
735 else
736 {
737 // Use lo-gain if it can be used
738 if (sig.IsLoGainValid() && fHiLoConv[pixidx]>0.5)
739 {
740 signal = sig.GetExtractedSignalLoGain()*fHiLoConv[pixidx];
741 signalErr = sig.GetExtractedSignalLoGain()*fHiLoConvErr[pixidx];
742 ok = kTRUE;
743 }
744 }
745
746 // In this cases the signal extraction has failed.
747 if (!ok)
748 bpix.SetUnsuitable(MBadPixelsPix::kUnsuitableEvt);
749
750 const Float_t nphot = signal * fCalibConsts [pixidx];
751 const Float_t nphotErr = TMath::Sqrt(TMath::Abs(nphot)) * fCalibFFactors[pixidx];
752
753 fCerPhotEvt->AddPixel(pixidx, nphot, nphotErr);
754
755 if (!bpix.IsUnsuitable())
756 {
757 if (sig.IsHiGainSaturated())
758 numsathi++;
759
760 if (sig.IsLoGainSaturated())
761 numsatlo++;
762 }
763 } /* if (data) */
764
765 if (pedestal)
766 {
767 TIter NextPed(&fPedestalCams);
768 TIter NextPhot(&fPedPhotCams);
769
770 MPedestalCam *pedestal = 0;
771 MPedPhotCam *pedphot = 0;
772
773 const Float_t pedmeancalib = slices *fCalibConsts[pixidx];
774 const Float_t pedrmscalib = sqrtslices*fCalibConsts[pixidx];
775
776 while ((pedestal=(MPedestalCam*)NextPed()) &&
777 (pedphot =(MPedPhotCam*)NextPhot()))
778 {
779 // pedestals/(used FADC slices) in [number of photons]
780 const Float_t mean = (*pedestal)[pixidx].GetPedestal() *pedmeancalib;
781 const Float_t rms = (*pedestal)[pixidx].GetPedestalRms()*pedrmscalib;
782
783 (*pedphot)[pixidx].Set(mean, rms);
784 pedphot->SetReadyToSave();
785 //break;
786 }
787 } /* if (pedestal) */
788 }
789
790 // Now we should take the bias (MPedPhotExtractor/Mean) and
791 // pedestal rms (MPedPhotExtractorRndm/Rms) and store it
792 // into MSignalPix
793
794 if (data)
795 {
796 fCerPhotEvt->SetNumPixelsSaturated(numsathi, numsatlo);
797 fCerPhotEvt->SetReadyToSave();
798 }
799 return kTRUE;
800}
801
802// --------------------------------------------------------------------------
803//
804// Apply the calibration factors to the extracted signal according to the
805// selected calibration method
806//
807Int_t MCalibrateData::Process()
808{
809 return Calibrate(fCalibrationMode!=kSkip, TestPedestalFlag(kEvent));
810}
811
812// --------------------------------------------------------------------------
813//
814// Implementation of SavePrimitive. Used to write the call to a constructor
815// to a macro. In the original root implementation it is used to write
816// gui elements to a macro-file.
817//
818void MCalibrateData::StreamPrimitive(ostream &out) const
819{
820 out << " " << ClassName() << " " << GetUniqueName() << "(\"";
821 out << "\"" << fName << "\", \"" << fTitle << "\");" << endl;
822
823 if (TestPedestalFlag(kEvent))
824 out << " " << GetUniqueName() << ".EnablePedestalType(MCalibrateData::kEvent)" << endl;
825 if (TestPedestalFlag(kRun))
826 out << " " << GetUniqueName() << ".EnablePedestalType(MCalibrateData::kRun)" << endl;
827
828 if (fCalibrationMode != gkDefault)
829 {
830 out << " " << GetUniqueName() << ".SetCalibrationMode(MCalibrateData::";
831 switch (fCalibrationMode)
832 {
833 case kSkip: out << "kSkip"; break;
834 case kNone: out << "kNone"; break;
835 case kFlatCharge: out << "kFlatCharge"; break;
836 case kBlindPixel: out << "kBlindPixel"; break;
837 case kFfactor: out << "kFfactor"; break;
838 case kPinDiode: out << "kPinDiode"; break;
839 case kCombined: out << "kCombined"; break;
840 case kDummy: out << "kDummy"; break;
841 default: out << (int)fCalibrationMode; break;
842 }
843 out << ")" << endl;
844 }
845
846 TIter Next(&fNamesPedestal);
847 TObject *o=0;
848 while ((o=Next()))
849 {
850 out << " " << GetUniqueName() << ".AddPedestal(\"";
851 out << o->GetName() << "\", \"" << o->GetTitle() << "\");" << endl;
852 }
853}
854
855// --------------------------------------------------------------------------
856//
857// Read the setup from a TEnv, eg:
858// MJPedestal.MCalibrateDate.PedestalFlag: no,run,event
859// MJPedestal.MCalibrateDate.CalibrationMode: skip,none,flatcharge,blindpixel,ffactor,pindiode,combined,dummy,default
860//
861Int_t MCalibrateData::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
862{
863 Bool_t rc = kFALSE;
864 if (IsEnvDefined(env, prefix, "PedestalFlag", print))
865 {
866 rc = kTRUE;
867 TString s = GetEnvValue(env, prefix, "PedestalFlag", "");
868 s.ToLower();
869 if (s.BeginsWith("no"))
870 SetPedestalFlag(kNo);
871 if (s.BeginsWith("run"))
872 SetPedestalFlag(kRun);
873 if (s.BeginsWith("event"))
874 SetPedestalFlag(kEvent);
875 }
876
877 if (IsEnvDefined(env, prefix, "CalibrationMode", print))
878 {
879 rc = kTRUE;
880 TString s = GetEnvValue(env, prefix, "CalibrationMode", "");
881 s.ToLower();
882 if (s.BeginsWith("skip"))
883 SetCalibrationMode(kSkip);
884 if (s.BeginsWith("none"))
885 SetCalibrationMode(kNone);
886 if (s.BeginsWith("flatcharge"))
887 SetCalibrationMode(kFlatCharge);
888 if (s.BeginsWith("blindpixel"))
889 SetCalibrationMode(kBlindPixel);
890 if (s.BeginsWith("ffactor"))
891 SetCalibrationMode(kFfactor);
892 if (s.BeginsWith("pindiode"))
893 SetCalibrationMode(kPinDiode);
894 if (s.BeginsWith("combined"))
895 SetCalibrationMode(kCombined);
896 if (s.BeginsWith("dummy"))
897 SetCalibrationMode(kDummy);
898 if (s.BeginsWith("default"))
899 SetCalibrationMode();
900 }
901
902 if (IsEnvDefined(env, prefix, "SignalType", print))
903 {
904 rc = kTRUE;
905 TString s = GetEnvValue(env, prefix, "SignalType", "");
906 s.ToLower();
907 if (s.BeginsWith("phot"))
908 SetSignalType(kPhot);
909 if (s.BeginsWith("phe"))
910 SetSignalType(kPhe);
911 if (s.BeginsWith("default"))
912 SetSignalType();
913 }
914
915 if (IsEnvDefined(env, prefix, "CalibConvMinLimit", print))
916 {
917 fCalibConvMinLimit = GetEnvValue(env, prefix, "CalibConvMinLimit", fCalibConvMinLimit);
918 rc = kTRUE;
919 }
920
921 if (IsEnvDefined(env, prefix, "CalibConvMaxLimit", print))
922 {
923 fCalibConvMaxLimit = GetEnvValue(env, prefix, "CalibConvMaxLimit", fCalibConvMaxLimit);
924 rc = kTRUE;
925 }
926
927 if (IsEnvDefined(env, prefix, "ScaleFactor", print))
928 {
929 fScaleFactor = GetEnvValue(env, prefix, "ScaleFactor", fScaleFactor);
930 rc = kTRUE;
931 }
932
933 if (IsEnvDefined(env, prefix, "FileNameScale", print))
934 {
935 fFileNameScale = GetEnvValue(env, prefix, "FileNameScale", fFileNameScale);
936 rc = kTRUE;
937 }
938
939 return rc;
940}
941
942void MCalibrateData::Print(Option_t *o) const
943{
944
945 *fLog << all << GetDescriptor() << ":" << endl;
946
947 for (UInt_t pixidx=0; pixidx<fGeomCam->GetNumPixels(); pixidx++)
948 {
949 *fLog << all
950 << "Pixel: " << Form("%3i",pixidx)
951 << " CalibConst: " << Form("%4.2f",fCalibConsts[pixidx])
952 << " F-Factor: " << Form("%4.2f",fCalibFFactors[pixidx])
953 << " HiLoConv: " << Form("%4.2f",fHiLoConv[pixidx])
954 << endl;
955 }
956}
957
Note: See TracBrowser for help on using the repository browser.