source: branches/Mars_MC/mcalib/MCalibrateData.cc@ 17038

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