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

Last change on this file since 7374 was 7374, checked in by tbretz, 19 years ago
*** empty log message ***
File size: 32.3 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 << inf << "Additional scale factor forced to: 1 (MonteCarloRun)" << endl;
444 fScaleFactor = 1;
445 }
446 else
447 {
448 if (!fFileNameScale.IsNull())
449 {
450 if (gSystem->AccessPathName(fFileNameScale, kFileExists))
451 {
452 *fLog << err << "ERROR - Configuration file '" << fFileNameScale << "' doesn't exist... abort." << endl;
453 return kFALSE;
454 }
455
456 const Int_t p = header->GetRunStart().GetMagicPeriod();
457 const TString per = Form("%2d", p);
458
459 TEnv rc(fFileNameScale);
460 const Double_t scale = rc.GetValue(per, -1.);
461
462 if (scale<=0)
463 {
464 *fLog << err << "ERROR - No valid entry for scale factor found for period " << p << " in " << fFileNameScale << "... abort." << endl;
465 return kFALSE;
466 }
467
468 *fLog << inf << "New additional scale factor for period " << p << ": " << scale << endl;
469 fScaleFactor = scale;
470 }
471 else
472 *fLog << inf << "Additional scale factor set to: " << fScaleFactor << endl;
473 }
474
475 const Int_t npixels = fGeomCam->GetNumPixels();
476
477 if (fCalibrationMode > kNone)
478 {
479
480 const MCalibrationCam *chargecam = fIntensCalib ? fIntensCalib->GetCam() : fCalibrations;
481 if (chargecam->GetSize() != npixels)
482 {
483 *fLog << "Size mismatch between MGeomCam and MCalibrationChargeCam... abort!" << endl;
484 return kFALSE;
485 }
486
487 if (fBadPixels->GetSize() != npixels)
488 {
489 *fLog << "Size mismatch between MGeomCam and MBadPixelsCam... abort!" << endl;
490 return kFALSE;
491 }
492 }
493
494 fCalibConsts .Set(npixels);
495 fCalibFFactors.Set(npixels);
496 fHiLoConv .Set(npixels);
497 fHiLoConvErr .Set(npixels);
498
499 if (!UpdateConversionFactors())
500 return kFALSE;
501
502 if (TestPedestalFlag(kRun))
503 Calibrate(kFALSE, kTRUE);
504
505 return kTRUE;
506}
507
508// --------------------------------------------------------------------------
509//
510// Update the conversion factors and F-Factors from MCalibrationCams into
511// the arrays. Later, the here pre-calcualted conversion factors get simply
512// copied from memory.
513//
514// This function can be called from outside in case that the MCalibrationCams
515// have been updated...
516//
517Bool_t MCalibrateData::UpdateConversionFactors( const MCalibrationChargeCam *updatecam)
518{
519
520 *fLog << inf << GetDescriptor() << ": Updating Conversion Factors... " << endl;
521
522 fCalibConsts.Reset();
523 fCalibFFactors.Reset();
524 fHiLoConv.Reset();
525 fHiLoConvErr.Reset();
526
527 MCalibConstCam *constcam = fIntensConst ? fIntensConst->GetCam() : fCalibConstCam;
528
529 MCalibrationChargeCam *chargecam = NULL;
530 MCalibrationQECam *qecam = NULL;
531 if (updatecam)
532 {
533 chargecam = fCalibrations;
534 qecam = fQEs;
535 }
536 else
537 {
538 chargecam = fIntensCalib ? (MCalibrationChargeCam*)fIntensCalib->GetCam() : fCalibrations;
539 qecam = fIntensQE ? (MCalibrationQECam*) fIntensQE->GetCam() : fQEs;
540 }
541
542 //
543 // For the moment, we use only a dummy zenith for the calibration:
544 //
545 const Float_t zenith = -1.;
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(zenith)/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 ( zenith );
637 calibQEVar = qe.GetQECascadesBlindPixelVar( zenith );
638 break;
639
640 case kPinDiode:
641 if (!qe.IsPINDiodeMethodValid())
642 {
643 skip++;
644 continue;
645 }
646 calibQE = qe.GetQECascadesPINDiode ( zenith );
647 calibQEVar = qe.GetQECascadesPINDiodeVar( zenith );
648 break;
649
650 case kFfactor:
651 if (!pix.IsFFactorMethodValid())
652 {
653 skip++;
654 continue;
655 }
656 calibQE = qe.GetQECascadesFFactor ( zenith );
657 calibQEVar = qe.GetQECascadesFFactorVar( zenith );
658 break;
659
660 case kCombined:
661 if (!qe.IsCombinedMethodValid())
662 {
663 skip++;
664 continue;
665 }
666 calibQE = qe.GetQECascadesCombined ( zenith );
667 calibQEVar = qe.GetQECascadesCombinedVar( zenith );
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 const Float_t slices = fSignals->GetNumUsedHiGainFADCSlices();
739 const Float_t sqrtslices = TMath::Sqrt(slices);
740
741 Int_t numsatlo=0;
742 Int_t numsathi=0;
743
744 for (UInt_t pixidx=0; pixidx<npix; pixidx++)
745 {
746
747 if (data)
748 {
749 const MExtractedSignalPix &sig = (*fSignals)[pixidx];
750
751 Float_t signal = 0.;
752 Float_t signalErr = 0.;
753
754 if (sig.IsLoGainUsed())
755 {
756 if (fHiLoConv[pixidx] < 0.5)
757 {
758 signal = sig.GetExtractedSignalHiGain()*1.5;
759 signalErr = sig.GetExtractedSignalHiGain()*0.5;
760 }
761 else
762 {
763 const Float_t siglo = sig.GetExtractedSignalLoGain();
764
765 if (siglo > 0.1) // low-gain signal has been extracted successfully
766 {
767 signal = siglo*fHiLoConv [pixidx];
768 signalErr = siglo*fHiLoConvErr[pixidx];
769 }
770 else // low-gain signal has not been extracted successfully, get a rough estimate from the high-gain
771 {
772 signal = sig.GetExtractedSignalHiGain()*1.5;
773 signalErr = sig.GetExtractedSignalHiGain()*0.5;
774 }
775 }
776 }
777 else
778 {
779 if (sig.GetExtractedSignalHiGain() <= 9999.)
780 signal = sig.GetExtractedSignalHiGain();
781 }
782
783 const Float_t nphot = signal * fCalibConsts [pixidx];
784 const Float_t nphotErr = TMath::Sqrt(TMath::Abs(nphot)) * fCalibFFactors[pixidx];
785
786 fCerPhotEvt->AddPixel(pixidx, nphot, nphotErr);
787
788 if (sig.GetNumHiGainSaturated() > 0)
789 numsathi++;
790
791 if (sig.GetNumLoGainSaturated() > 0)
792 numsatlo++;
793 } /* if (data) */
794
795
796 if (pedestal)
797 {
798 TIter NextPed(&fPedestalCams);
799 TIter NextPhot(&fPedPhotCams);
800
801 MPedestalCam *pedestal = 0;
802 MPedPhotCam *pedphot = 0;
803
804 const Float_t pedmeancalib = slices *fCalibConsts[pixidx];
805 const Float_t pedrmscalib = sqrtslices*fCalibConsts[pixidx];
806
807 while ((pedestal=(MPedestalCam*)NextPed()) &&
808 (pedphot =(MPedPhotCam*)NextPhot()))
809 {
810 // pedestals/(used FADC slices) in [number of photons]
811 const Float_t mean = (*pedestal)[pixidx].GetPedestal() *pedmeancalib;
812 const Float_t rms = (*pedestal)[pixidx].GetPedestalRms()*pedrmscalib;
813
814 (*pedphot)[pixidx].Set(mean, rms);
815 pedphot->SetReadyToSave();
816 }
817 } /* if (pedestal) */
818 }
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(ofstream &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.