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

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