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

Last change on this file since 7189 was 7189, checked in by tbretz, 19 years ago
*** empty log message ***
File size: 31.1 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-2004
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 "MCalibrationIntensityChargeCam.h"
115#include "MCalibrationChargeCam.h"
116#include "MCalibrationChargePix.h"
117
118#include "MCalibrationIntensityQECam.h"
119#include "MCalibrationQECam.h"
120#include "MCalibrationQEPix.h"
121
122#include "MCalibrationIntensityConstCam.h"
123#include "MCalibConstCam.h"
124#include "MCalibConstPix.h"
125
126#include "MExtractedSignalCam.h"
127#include "MExtractedSignalPix.h"
128
129#include "MPedPhotCam.h"
130#include "MPedPhotPix.h"
131
132#include "MBadPixelsCam.h"
133#include "MBadPixelsPix.h"
134
135#include "MSignalCam.h"
136
137ClassImp(MCalibrateData);
138
139using namespace std;
140
141const Float_t MCalibrateData::gkCalibConvMinLimit = 0.01;
142const Float_t MCalibrateData::gkCalibConvMaxLimit = 5.;
143
144const MCalibrateData::CalibrationMode_t MCalibrateData::gkDefault = kFfactor;
145
146// --------------------------------------------------------------------------
147//
148// Default constructor.
149//
150// Sets all pointers to NULL
151//
152// Initializes:
153// - fCalibrationMode to kDefault
154// - fPedestalFlag to kNo
155//
156MCalibrateData::MCalibrateData(CalibrationMode_t calmode,const char *name, const char *title)
157 : fGeomCam(NULL), fBadPixels(NULL), fCalibrations(NULL), fIntensCalib(NULL),
158 fQEs(NULL), fIntensQE(NULL), fSignals(NULL), fCerPhotEvt(NULL), fCalibConstCam(NULL),
159 fIntensConst(NULL), fPedestalFlag(kNo), fSignalType(kPhot), fRenormFactor(1.),
160 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 fIntensConst = (MCalibrationIntensityConstCam*)pList->FindObject(AddSerialNumber("MCalibrationIntensityConstCam"));
233 if (fIntensConst)
234 *fLog << inf << "Found MCalibrationIntensityConstCam ... " << endl;
235 else
236 {
237 fCalibConstCam = (MCalibConstCam*)pList->FindCreateObj(AddSerialNumber("MCalibConstCam"));
238 if (!fCalibConstCam)
239 return kFALSE;
240 }
241 }
242
243 fCalibrations = 0;
244 fQEs = 0;
245 if (fCalibrationMode>kNone)
246 {
247 fIntensCalib = (MCalibrationIntensityChargeCam*)pList->FindObject(AddSerialNumber("MCalibrationIntensityChargeCam"));
248 if (fIntensCalib)
249 *fLog << inf << "Found MCalibrationIntensityChargeCam ... " << endl;
250
251 fCalibrations = (MCalibrationChargeCam*)pList->FindObject(AddSerialNumber("MCalibrationChargeCam"));
252 if (!fCalibrations)
253 {
254 *fLog << err << AddSerialNumber("MCalibrationChargeCam") << " not found ... aborting." << endl;
255 return kFALSE;
256 }
257
258 fIntensQE = (MCalibrationIntensityQECam*)pList->FindObject(AddSerialNumber("MCalibrationIntensityQECam"));
259 if (fIntensQE)
260 *fLog << inf << "Found MCalibrationIntensityQECam ... " << endl;
261
262 fQEs = (MCalibrationQECam*)pList->FindObject(AddSerialNumber("MCalibrationQECam"));
263 if (!fQEs)
264 {
265 *fLog << err << AddSerialNumber("MCalibrationQECam") << " not found ... aborting." << endl;
266 return kFALSE;
267 }
268
269 }
270
271 if (fNamesPedestal.GetSize()>0 && fPedestalFlag==kNo)
272 {
273 *fLog << warn << "Pedestal list contains entries, but mode is set to kNo... setting to kEvent." << endl;
274 fPedestalFlag = kEvent;
275 }
276
277 if (fPedestalFlag!=kNo)
278 {
279 if (fNamesPedestal.GetSize()==0)
280 {
281 *fLog << inf << "No container names specified... using default: MPedestalCam and MPedPhotCam." << endl;
282 AddPedestal();
283 }
284
285 fPedestalCams.Clear();
286 fPedPhotCams.Clear();
287
288 TIter Next(&fNamesPedestal);
289 TObject *o=0;
290 while ((o=Next()))
291 {
292 TObject *pedcam = pList->FindObject(AddSerialNumber(o->GetName()), "MPedestalCam");
293 if (!pedcam)
294 {
295 *fLog << err << AddSerialNumber(o->GetName()) << " [MPedestalCam] not found ... aborting" << endl;
296 return kFALSE;
297 }
298 TObject *pedphot = pList->FindCreateObj("MPedPhotCam", AddSerialNumber(o->GetTitle()));
299 if (!pedphot)
300 return kFALSE;
301
302 fPedestalCams.Add(pedcam);
303 fPedPhotCams.Add(pedphot);
304 }
305 }
306
307 switch (fSignalType)
308 {
309 case kPhe:
310 //
311 // Average QE for C-photons, for pixels in the inner pixel region ("area 0"),
312 // used later to convert from C-photons into "equivalent phes":
313 //
314 if (fQEs->GetAverageAreas() > 0)
315 {
316 MCalibrationQEPix& avqepix = (MCalibrationQEPix&)(fQEs->GetAverageArea(0));
317 fRenormFactor = avqepix.GetAverageQE();
318 }
319 else // Average areas not yet initialized, use default value
320 fRenormFactor = MCalibrationQEPix::gkDefaultAverageQE;
321
322 fRenormFactor = MCalibrationQEPix::gkDefaultAverageQE;
323 break;
324 case kPhot:
325 fRenormFactor = 1.;
326 break;
327 }
328
329 fCalibConsts.Reset();
330 fCalibFFactors.Reset();
331 fHiLoConv.Reset();
332 fHiLoConvErr.Reset();
333
334 return kTRUE;
335}
336
337// --------------------------------------------------------------------------
338//
339// The ReInit searches for the following input containers:
340//
341// - MGeomCam
342//
343// Check for validity of the selected calibration method, switch to a
344// different one in case of need
345//
346// Fill the MPedPhotCam container using the information from MPedestalCam,
347// MExtractedSignalCam and MCalibrationCam
348//
349Bool_t MCalibrateData::ReInit(MParList *pList)
350{
351 MRawRunHeader *header = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
352 if (!header)
353 {
354 *fLog << err << "MRawRunHeader not found... abort." << endl;
355 return kFALSE;
356 }
357
358 fGeomCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
359 if (!fGeomCam)
360 {
361 *fLog << err << "No MGeomCam found... aborting." << endl;
362 return kFALSE;
363 }
364
365 // Sizes might have changed
366 if (fPedestalFlag!=kNo)
367 {
368 TIter Next(&fPedestalCams);
369 MPedestalCam *cam=0;
370 while ((cam=(MPedestalCam*)Next()))
371 if ((Int_t)cam->GetSize() != fSignals->GetSize())
372 {
373 *fLog << err << "Size mismatch of " << cam->GetDescriptor() << " and MCalibrationCam... abort." << endl;
374 return kFALSE;
375 }
376 }
377
378 const MCalibrationQECam *qecam = fIntensQE ? (MCalibrationQECam*)fIntensQE->GetCam() : fQEs;
379 if(fCalibrationMode == kBlindPixel && !qecam->IsBlindPixelMethodValid())
380 {
381 *fLog << warn << "Blind pixel calibration method not valid, switching to F-factor method" << endl;
382 fCalibrationMode = kFfactor;
383 }
384
385 if(fCalibrationMode == kPinDiode && !qecam->IsPINDiodeMethodValid())
386 {
387 *fLog << warn << "PIN diode calibration method not valid, switching to F-factor method" << endl;
388 fCalibrationMode = kFfactor;
389 }
390
391 if(fCalibrationMode == kCombined && !qecam->IsCombinedMethodValid())
392 {
393 *fLog << warn << "Combined calibration method not valid, switching to F-factor method" << endl;
394 fCalibrationMode = kFfactor;
395 }
396
397 //
398 // output information or warnings:
399 //
400 switch(fCalibrationMode)
401 {
402 case kBlindPixel:
403 break;
404 case kFfactor:
405 break;
406 case kPinDiode:
407 *fLog << err << "PIN Diode Calibration mode not yet available!" << endl;
408 return kFALSE;
409 break;
410 case kCombined:
411 *fLog << err << "Combined Calibration mode not yet available!" << endl;
412 return kFALSE;
413 break;
414 case kFlatCharge:
415 *fLog << warn << "WARNING - Flat-fielding charges - only for muon calibration!" << endl;
416 break;
417 case kDummy:
418 *fLog << warn << "WARNING - Dummy calibration, no calibration applied!" << endl;
419 break;
420 case kNone:
421 *fLog << warn << "WARNING - No calibration applied!" << endl;
422 break;
423 default:
424 *fLog << warn << "WARNING - Calibration mode value (" << fCalibrationMode << ") not known" << endl;
425 return kFALSE;
426 }
427
428 //
429 // output information or warnings:
430 //
431 switch(fSignalType)
432 {
433 case kPhe:
434 *fLog << inf << "Calibrating in units of equivalent (outer/inner=4) photo-electrons." << endl;
435 break;
436 case kPhot:
437 *fLog << inf << "Calibrating in units of photons." << endl;
438 break;
439 }
440
441 if (header->IsMonteCarloRun())
442 {
443 *fLog << "Additional scale factor: 1 (MonteCarloRun)" << endl;
444 fScaleFactor = 1;
445 }
446 else
447 *fLog << "Additional scale factor: " << fScaleFactor << endl;
448
449 const Int_t npixels = fGeomCam->GetNumPixels();
450
451 if (fCalibrationMode > kNone)
452 {
453
454 const MCalibrationCam *chargecam = fIntensCalib ? fIntensCalib->GetCam() : fCalibrations;
455 if (chargecam->GetSize() != npixels)
456 {
457 *fLog << "Size mismatch between MGeomCam and MCalibrationChargeCam... abort!" << endl;
458 return kFALSE;
459 }
460
461 if (fBadPixels->GetSize() != npixels)
462 {
463 *fLog << "Size mismatch between MGeomCam and MBadPixelsCam... abort!" << endl;
464 return kFALSE;
465 }
466 }
467
468 fCalibConsts .Set(npixels);
469 fCalibFFactors.Set(npixels);
470 fHiLoConv .Set(npixels);
471 fHiLoConvErr .Set(npixels);
472
473 if (!UpdateConversionFactors())
474 return kFALSE;
475
476 if (TestPedestalFlag(kRun))
477 Calibrate(kFALSE, kTRUE);
478
479 return kTRUE;
480}
481
482// --------------------------------------------------------------------------
483//
484// Update the conversion factors and F-Factors from MCalibrationCams into
485// the arrays. Later, the here pre-calcualted conversion factors get simply
486// copied from memory.
487//
488// This function can be called from outside in case that the MCalibrationCams
489// have been updated...
490//
491Bool_t MCalibrateData::UpdateConversionFactors( const MCalibrationChargeCam *updatecam)
492{
493
494 *fLog << inf << GetDescriptor() << ": Updating Conversion Factors... " << endl;
495
496 fCalibConsts.Reset();
497 fCalibFFactors.Reset();
498 fHiLoConv.Reset();
499 fHiLoConvErr.Reset();
500
501 MCalibConstCam *constcam = fIntensConst ? fIntensConst->GetCam() : fCalibConstCam;
502
503 MCalibrationChargeCam *chargecam = NULL;
504 MCalibrationQECam *qecam = NULL;
505 if (updatecam)
506 {
507 chargecam = fCalibrations;
508 qecam = fQEs;
509 }
510 else
511 {
512 chargecam = fIntensCalib ? (MCalibrationChargeCam*)fIntensCalib->GetCam() : fCalibrations;
513 qecam = fIntensQE ? (MCalibrationQECam*) fIntensQE->GetCam() : fQEs;
514 }
515
516 //
517 // For the moment, we use only a dummy zenith for the calibration:
518 //
519 const Float_t zenith = -1.;
520
521 UInt_t skip = 0;
522
523 for (UInt_t pixidx=0; pixidx<fGeomCam->GetNumPixels(); pixidx++)
524 {
525
526 Float_t hiloconv = 1.;
527 Float_t hiloconverr = 0.;
528 Float_t calibConv = 1.;
529 Float_t calibConvVar = 0.;
530 Float_t calibFFactor = 0.;
531
532 Float_t calibQE = 1.;
533 Float_t calibQEVar = 0.;
534
535 Float_t calibUpdate = 1.;
536
537 MCalibConstPix &cpix = (*constcam)[pixidx];
538
539 if(fCalibrationMode!=kNone)
540 {
541 if ((*fBadPixels)[pixidx].IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
542 {
543 skip++;
544 calibConv = -1.;
545 continue;
546 }
547
548 const MCalibrationChargePix &pix = (MCalibrationChargePix&)(*chargecam)[pixidx];
549 const MCalibrationChargePix &avpix = (MCalibrationChargePix&)chargecam->GetAverageArea(0);
550
551 hiloconv = pix.GetConversionHiLo ();
552 hiloconverr= pix.GetConversionHiLoSigma();
553
554 calibConv = pix.GetMeanConvFADC2Phe();
555 calibConvVar = pix.GetMeanConvFADC2PheVar();
556 calibFFactor = pix.GetMeanFFactorFADC2Phot();
557
558 const MCalibrationQEPix &qe = (MCalibrationQEPix&)(*qecam)[pixidx];
559
560 if (updatecam)
561 {
562 const MCalibrationChargePix &upix = (MCalibrationChargePix&)(*updatecam)[pixidx];
563
564 //
565 // Correct for the possible change in amplification of the individual pixels chain
566 //
567 const Float_t pixmean = upix.GetConvertedMean();
568 calibUpdate = pixmean==0 ? 1 : pix.GetConvertedMean()/pixmean;
569
570 //
571 // Correct for global shifts in light emission
572 //
573 const MCalibrationChargePix &ugpix = (MCalibrationChargePix&)updatecam->GetAverageArea(0);
574
575 const Float_t globmean = avpix.GetConvertedMean();
576 calibUpdate = globmean==0 ? 1 : ugpix.GetConvertedMean()/globmean;
577
578 MBadPixelsPix &ubad = (MBadPixelsPix&)updatecam->GetAverageBadArea(0);
579 if (ubad.IsUncalibrated(MBadPixelsPix::kChargeIsPedestal))
580 {
581 *fLog << warn << GetDescriptor() << ": Mean charge in inner pixels is smaller than 3 ped. RMS." << endl;
582 *fLog << "Maybe calibration pulses have been lost!" << endl;
583 calibUpdate = 1.;
584
585 }
586 }
587
588 switch(fCalibrationMode)
589 {
590 case kFlatCharge:
591 {
592 calibConv = avpix.GetConvertedMean()
593 / (pix.GetConvertedMean() * fGeomCam->GetPixRatio(pixidx));
594 calibConvVar = (avpix.GetMeanRelVar() + pix.GetMeanRelVar()) * calibConv * calibConv;
595 if (pix.IsFFactorMethodValid())
596 {
597 const Float_t convmin1 = qe.GetQECascadesFFactor(zenith)/pix.GetMeanConvFADC2Phe();
598 if (convmin1 > 0)
599 calibFFactor *= TMath::Sqrt(convmin1);
600 else
601 calibFFactor = -1.;
602 }
603 break;
604 }
605 case kBlindPixel:
606 if (!qe.IsBlindPixelMethodValid())
607 {
608 skip++;
609 continue;
610 }
611 calibQE = qe.GetQECascadesBlindPixel ( zenith );
612 calibQEVar = qe.GetQECascadesBlindPixelVar( zenith );
613 break;
614
615 case kPinDiode:
616 if (!qe.IsPINDiodeMethodValid())
617 {
618 skip++;
619 continue;
620 }
621 calibQE = qe.GetQECascadesPINDiode ( zenith );
622 calibQEVar = qe.GetQECascadesPINDiodeVar( zenith );
623 break;
624
625 case kFfactor:
626 if (!pix.IsFFactorMethodValid())
627 {
628 skip++;
629 continue;
630 }
631 calibQE = qe.GetQECascadesFFactor ( zenith );
632 calibQEVar = qe.GetQECascadesFFactorVar( zenith );
633 break;
634
635 case kCombined:
636 if (!qe.IsCombinedMethodValid())
637 {
638 skip++;
639 continue;
640 }
641 calibQE = qe.GetQECascadesCombined ( zenith );
642 calibQEVar = qe.GetQECascadesCombinedVar( zenith );
643 break;
644
645 case kDummy:
646 hiloconv = 1.;
647 hiloconverr = 0.;
648 calibUpdate = 1.;
649 break;
650
651 default:
652 break;
653 } /* switch calibration mode */
654 } /* if(fCalibrationMode!=kNone) */
655 else
656 {
657 calibConv = 1./fGeomCam->GetPixRatio(pixidx);
658 }
659
660 calibConv /= calibQE;
661
662 if (calibConv > 0.00001 && calibQE > 0.00001)
663 {
664 calibConvVar = calibConvVar/(calibConv*calibConv) + calibQEVar/(calibQE*calibQE);
665 calibConvVar *= (calibConv*calibConv);
666 }
667
668 calibConv *= fRenormFactor*fScaleFactor * calibUpdate;
669 calibFFactor *= TMath::Sqrt(fRenormFactor*fScaleFactor);
670
671 fHiLoConv [pixidx] = hiloconv;
672 fHiLoConvErr [pixidx] = hiloconverr;
673 fCalibConsts [pixidx] = calibConv;
674 fCalibFFactors[pixidx] = calibFFactor;
675
676 if (calibConv < fCalibConvMinLimit || calibConv > fCalibConvMaxLimit)
677 {
678 (*fBadPixels)[pixidx].SetUnsuitable(MBadPixelsPix::kUnsuitableRun);
679 calibConv = -1.;
680 calibFFactor = -1.;
681 *fLog << warn << "Conversion factor of Pixel " << pixidx << " out of range " << endl;
682 }
683 cpix.SetCalibConst(calibConv);
684 cpix.SetCalibFFactor(calibFFactor);
685
686 } /* for (Int_t pixidx=0; pixidx<fGeomCam->GetNumPixels(); pixidx++) */
687
688 if (skip>fGeomCam->GetNumPixels()*0.9)
689 {
690 *fLog << warn << GetDescriptor()
691 << ": WARNING - GetConversionFactor has skipped more than 90% of the pixels... abort." << endl;
692 return kFALSE;
693 }
694
695 return kTRUE;
696}
697
698
699// --------------------------------------------------------------------------
700//
701// Apply the conversion factors and F-Factors from the arrays to the data.
702//
703// The flags 'data' and 'pedestal' decide whether the signal and/or the pedetals
704// shall be calibrated, respectively.
705//
706Int_t MCalibrateData::Calibrate(Bool_t data, Bool_t pedestal) const
707{
708 if (!data && !pedestal)
709 return kTRUE;
710
711 const UInt_t npix = fSignals->GetSize();
712 const Float_t slices = fSignals->GetNumUsedHiGainFADCSlices();
713 const Float_t sqrtslices = TMath::Sqrt(slices);
714
715 Int_t numsatlo=0;
716 Int_t numsathi=0;
717
718 for (UInt_t pixidx=0; pixidx<npix; pixidx++)
719 {
720
721 if (data)
722 {
723 const MExtractedSignalPix &sig = (*fSignals)[pixidx];
724
725 Float_t signal = 0.;
726 Float_t signalErr = 0.;
727
728 if (sig.IsLoGainUsed())
729 {
730 if (fHiLoConv[pixidx] < 0.5)
731 {
732 signal = sig.GetExtractedSignalHiGain()*1.5;
733 signalErr = sig.GetExtractedSignalHiGain()*0.5;
734 }
735 else
736 {
737 const Float_t siglo = sig.GetExtractedSignalLoGain();
738
739 if (siglo > 0.1) // low-gain signal has been extracted successfully
740 {
741 signal = siglo*fHiLoConv [pixidx];
742 signalErr = siglo*fHiLoConvErr[pixidx];
743 }
744 else // low-gain signal has not been extracted successfully, get a rough estimate from the high-gain
745 {
746 signal = sig.GetExtractedSignalHiGain()*1.5;
747 signalErr = sig.GetExtractedSignalHiGain()*0.5;
748 }
749 }
750 }
751 else
752 {
753 if (sig.GetExtractedSignalHiGain() <= 9999.)
754 signal = sig.GetExtractedSignalHiGain();
755 }
756
757 const Float_t nphot = signal * fCalibConsts [pixidx];
758 const Float_t nphotErr = TMath::Sqrt(TMath::Abs(nphot)) * fCalibFFactors[pixidx];
759
760 fCerPhotEvt->AddPixel(pixidx, nphot, nphotErr);
761
762 if (sig.GetNumHiGainSaturated() > 0)
763 numsathi++;
764
765 if (sig.GetNumLoGainSaturated() > 0)
766 numsatlo++;
767 } /* if (data) */
768
769
770 if (pedestal)
771 {
772 TIter NextPed(&fPedestalCams);
773 TIter NextPhot(&fPedPhotCams);
774
775 MPedestalCam *pedestal = 0;
776 MPedPhotCam *pedphot = 0;
777
778 const Float_t pedmeancalib = slices *fCalibConsts[pixidx];
779 const Float_t pedrmscalib = sqrtslices*fCalibConsts[pixidx];
780
781 while ((pedestal=(MPedestalCam*)NextPed()) &&
782 (pedphot =(MPedPhotCam*)NextPhot()))
783 {
784 // pedestals/(used FADC slices) in [number of photons]
785 const Float_t mean = (*pedestal)[pixidx].GetPedestal() *pedmeancalib;
786 const Float_t rms = (*pedestal)[pixidx].GetPedestalRms()*pedrmscalib;
787
788 (*pedphot)[pixidx].Set(mean, rms);
789 pedphot->SetReadyToSave();
790 }
791 } /* if (pedestal) */
792 }
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(ofstream &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 return rc;
934}
935
936void MCalibrateData::Print(Option_t *o) const
937{
938
939 *fLog << all << GetDescriptor() << ":" << endl;
940
941 for (UInt_t pixidx=0; pixidx<fGeomCam->GetNumPixels(); pixidx++)
942 {
943 *fLog << all
944 << "Pixel: " << Form("%3i",pixidx)
945 << " CalibConst: " << Form("%4.2f",fCalibConsts[pixidx])
946 << " F-Factor: " << Form("%4.2f",fCalibFFactors[pixidx])
947 << " HiLoConv: " << Form("%4.2f",fHiLoConv[pixidx])
948 << endl;
949 }
950}
951
Note: See TracBrowser for help on using the repository browser.