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

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