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

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