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

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