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

Last change on this file since 5518 was 5196, checked in by gaug, 21 years ago
*** empty log message ***
File size: 21.4 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 MCerPhotEvt.
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 SetNamePedADCContainer() and/or SetNamePedPhotContainer() you
69// can control the name of the 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.SetNamePedADCContainer("MPedestalFromLoGain");
79// cal2.SetNamePedPhotContainer("MPedPhotFromLoGain")
80// cal2.SetPedestalFlag(MCalibrateData::kEvent);
81//
82//
83// Input Containers:
84// [MPedestalCam]
85// [MExtractedSignalCam]
86// [MCalibrationChargeCam]
87// [MCalibrationQECam]
88// MBadPixelsCam
89//
90// Output Containers:
91// [MPedPhotCam]
92// [MCerPhotEvt]
93//
94// See also: MJCalibration, MJPedestal, MJExtractSignal, MJExtractCalibTest
95//
96//////////////////////////////////////////////////////////////////////////////
97#include "MCalibrateData.h"
98
99#include <fstream>
100
101#include <TEnv.h>
102
103#include "MLog.h"
104#include "MLogManip.h"
105
106#include "MParList.h"
107#include "MH.h"
108
109#include "MGeomCam.h"
110
111#include "MPedestalCam.h"
112#include "MPedestalPix.h"
113
114#include "MCalibrationChargeCam.h"
115#include "MCalibrationChargePix.h"
116
117#include "MCalibrationQECam.h"
118#include "MCalibrationQEPix.h"
119
120#include "MExtractedSignalCam.h"
121#include "MExtractedSignalPix.h"
122
123#include "MPedPhotCam.h"
124#include "MPedPhotPix.h"
125
126#include "MBadPixelsCam.h"
127#include "MBadPixelsPix.h"
128
129#include "MCerPhotEvt.h"
130
131ClassImp(MCalibrateData);
132
133using namespace std;
134
135const TString MCalibrateData::fgNamePedADCContainer = "MPedestalCam";
136const TString MCalibrateData::fgNamePedPhotContainer = "MPedPhotCam";
137// --------------------------------------------------------------------------
138//
139// Default constructor.
140//
141// Sets all pointers to NULL
142//
143// Initializes:
144// - fCalibrationMode to kDefault
145// - fPedestalFlag to kRun
146// - fNamePedADCRunContainer to "MPedestalCam"
147// - fNamePedPhotRunContainer to "MPedPhotCam"
148//
149MCalibrateData::MCalibrateData(CalibrationMode_t calmode,const char *name, const char *title)
150 : fGeomCam(NULL), fPedestal(NULL),
151 fBadPixels(NULL), fCalibrations(NULL), fQEs(NULL), fSignals(NULL),
152 fPedPhot(NULL), fCerPhotEvt(NULL), fPedestalFlag(kNo)
153{
154
155 fName = name ? name : "MCalibrateData";
156 fTitle = title ? title : "Task to calculate the number of photons in one event";
157
158 SetCalibrationMode(calmode);
159
160 SetNamePedADCContainer();
161 SetNamePedPhotContainer();
162}
163
164// --------------------------------------------------------------------------
165//
166// The PreProcess searches for the following input containers:
167//
168// - MGeomCam
169// - MPedestalCam
170// - MCalibrationChargeCam
171// - MCalibrationQECam
172// - MExtractedSignalCam
173// - MBadPixelsCam
174//
175// The following output containers are also searched and created if
176// they were not found:
177//
178// - MPedPhotCam
179// - MCerPhotEvt
180//
181Int_t MCalibrateData::PreProcess(MParList *pList)
182{
183 // input containers
184
185 fBadPixels = (MBadPixelsCam*)pList->FindObject(AddSerialNumber("MBadPixelsCam"));
186 if (!fBadPixels)
187 {
188 *fLog << err << AddSerialNumber("MBadPixelsCam") << " not found ... aborting" << endl;
189 return kFALSE;
190 }
191
192 fSignals = 0;
193 fCerPhotEvt = 0;
194 if (fCalibrationMode>kSkip)
195 {
196 fSignals = (MExtractedSignalCam*)pList->FindObject(AddSerialNumber("MExtractedSignalCam"));
197 if (!fSignals)
198 {
199 *fLog << err << AddSerialNumber("MExtractedSignalCam") << " not found ... aborting" << endl;
200 return kFALSE;
201 }
202
203 fCerPhotEvt = (MCerPhotEvt*)pList->FindCreateObj(AddSerialNumber("MCerPhotEvt"));
204 if (!fCerPhotEvt)
205 return kFALSE;
206 }
207
208 fCalibrations = 0;
209 fQEs = 0;
210 if (fCalibrationMode>kNone)
211 {
212 fCalibrations = (MCalibrationChargeCam*)pList->FindObject(AddSerialNumber("MCalibrationChargeCam"));
213 if (!fCalibrations)
214 {
215 *fLog << err << AddSerialNumber("MCalibrationChargeCam") << " not found ... aborting." << endl;
216 return kFALSE;
217 }
218
219 fQEs = (MCalibrationQECam*)pList->FindObject(AddSerialNumber("MCalibrationQECam"));
220 if (!fQEs)
221 {
222 *fLog << err << AddSerialNumber("MCalibrationQECam") << " not found ... aborting." << endl;
223 return kFALSE;
224 }
225 }
226
227 fPedestal = 0;
228 fPedPhot = 0;
229 if (fPedestalFlag)
230 {
231 fPedestal = (MPedestalCam*)pList->FindObject(AddSerialNumber(fNamePedADCContainer), "MPedestalCam");
232 if (!fPedestal)
233 {
234 *fLog << err << AddSerialNumber(fNamePedADCContainer) << " [MPedestalCam] not found ... aborting" << endl;
235 return kFALSE;
236 }
237
238 fPedPhot = (MPedPhotCam*)pList->FindCreateObj("MPedPhotCam", AddSerialNumber(fNamePedPhotContainer));
239 if (!fPedPhot)
240 return kFALSE;
241 }
242
243 return kTRUE;
244}
245
246// --------------------------------------------------------------------------
247//
248// The ReInit searches for the following input containers:
249//
250// - MGeomCam
251//
252// Check for validity of the selected calibration method, switch to a
253// different one in case of need
254//
255// Fill the MPedPhotCam container using the information from MPedestalCam,
256// MExtractedSignalCam and MCalibrationCam
257//
258Bool_t MCalibrateData::ReInit(MParList *pList)
259{
260 fGeomCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
261 if (!fGeomCam)
262 {
263 *fLog << err << "No MGeomCam found... aborting." << endl;
264 return kFALSE;
265 }
266
267 // Sizes might have changed
268 if (fPedestalFlag && (Int_t)fPedestal->GetSize() != fSignals->GetSize())
269 {
270 *fLog << err << "Size mismatch of MPedestalCam and MCalibrationCam... abort." << endl;
271 return kFALSE;
272 }
273
274 if(fCalibrationMode == kBlindPixel && !fQEs->IsBlindPixelMethodValid())
275 {
276 *fLog << warn << "Blind pixel calibration method not valid, switching to F-factor method" << endl;
277 fCalibrationMode = kFfactor;
278 }
279
280 if(fCalibrationMode == kPinDiode && !fQEs->IsPINDiodeMethodValid())
281 {
282 *fLog << warn << "PIN diode calibration method not valid, switching to F-factor method" << endl;
283 fCalibrationMode = kFfactor;
284 }
285
286 if(fCalibrationMode == kCombined && !fQEs->IsCombinedMethodValid())
287 {
288 *fLog << warn << "Combined calibration method not valid, switching to F-factor method" << endl;
289 fCalibrationMode = kFfactor;
290 }
291
292 //
293 // output information or warnings:
294 //
295 switch(fCalibrationMode)
296 {
297 case kBlindPixel:
298 break;
299 case kFfactor:
300 break;
301 case kPinDiode:
302 *fLog << err << "PIN Diode Calibration mode not yet available!" << endl;
303 return kFALSE;
304 break;
305 case kCombined:
306 *fLog << err << "Combined Calibration mode not yet available!" << endl;
307 return kFALSE;
308 break;
309 case kFlatCharge:
310 *fLog << warn << "WARNING - Flat-fielding charges - only for muon calibration!" << endl;
311 break;
312 case kDummy:
313 *fLog << warn << "WARNING - Dummy calibration, no calibration applied!" << endl;
314 break;
315 case kNone:
316 *fLog << warn << "WARNING - No calibration applied!" << endl;
317 break;
318 default:
319 *fLog << warn << "WARNING - Calibration mode value (" << fCalibrationMode << ") not known" << endl;
320 return kFALSE;
321 }
322
323 if (TestPedestalFlag(kRun))
324 Calibrate(kFALSE, kTRUE);
325
326 return kTRUE;
327}
328
329// --------------------------------------------------------------------------
330//
331// Get conversion factor and its error from MCalibrationCam
332//
333Bool_t MCalibrateData::GetConversionFactor(UInt_t pixidx, Float_t &hiloconv, Float_t &hiloconverr,
334 Float_t &calibConv, Float_t &calibConvVar,
335 Float_t &calibFFactor) const
336{
337 //
338 // For the moment, we use only a dummy zenith for the calibration:
339 //
340 const Float_t zenith = -1.;
341
342 hiloconv = 1.;
343 hiloconverr = 0.;
344 calibConv = 1.;
345 calibConvVar = 0.;
346 calibFFactor = 0.;
347
348 Float_t calibQE = 1.;
349 Float_t calibQEVar = 0.;
350
351 if(fCalibrationMode!=kNone)
352 {
353 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCalibrations)[pixidx];
354
355 hiloconv = pix.GetConversionHiLo ();
356 hiloconverr= pix.GetConversionHiLoErr();
357
358 if ((*fBadPixels)[pixidx].IsUnsuitable())
359 return kFALSE;
360
361 calibConv = pix.GetMeanConvFADC2Phe();
362 calibConvVar = pix.GetMeanConvFADC2PheVar();
363 calibFFactor = pix.GetMeanFFactorFADC2Phot();
364
365 MCalibrationQEPix &qe = (MCalibrationQEPix&) (*fQEs)[pixidx];
366
367 switch(fCalibrationMode)
368 {
369 case kFlatCharge:
370 {
371 MCalibrationChargePix &avpix = (MCalibrationChargePix&)fCalibrations->GetAverageArea(0);
372 calibConv = avpix.GetMean() / (pix.GetMean() * fGeomCam->GetPixRatio(pixidx));
373 calibConvVar = (avpix.GetMeanRelVar() + pix.GetMeanRelVar()) * calibConv * calibConv;
374 if (pix.IsFFactorMethodValid())
375 {
376 const Float_t convmin1 = qe.GetQECascadesFFactor(zenith)/pix.GetMeanConvFADC2Phe();
377 if (convmin1 > 0)
378 calibFFactor *= TMath::Sqrt(convmin1);
379 else
380 calibFFactor = -1.;
381 }
382 }
383 break;
384
385 case kBlindPixel:
386 if (!qe.IsBlindPixelMethodValid())
387 return kFALSE;
388 calibQE = qe.GetQECascadesBlindPixel ( zenith );
389 calibQEVar = qe.GetQECascadesBlindPixelVar( zenith );
390 break;
391
392 case kPinDiode:
393 if (!qe.IsPINDiodeMethodValid())
394 return kFALSE;
395 calibQE = qe.GetQECascadesPINDiode ( zenith );
396 calibQEVar = qe.GetQECascadesPINDiodeVar( zenith );
397 break;
398
399 case kFfactor:
400 if (!pix.IsFFactorMethodValid())
401 return kFALSE;
402 calibQE = qe.GetQECascadesFFactor ( zenith );
403 calibQEVar = qe.GetQECascadesFFactorVar( zenith );
404 break;
405
406 case kCombined:
407 if (!qe.IsCombinedMethodValid())
408 return kFALSE;
409 calibQE = qe.GetQECascadesCombined ( zenith );
410 calibQEVar = qe.GetQECascadesCombinedVar( zenith );
411 break;
412
413 case kDummy:
414 hiloconv = 1.;
415 hiloconverr = 0.;
416 break;
417 } /* switch calibration mode */
418 } /* if(fCalibrationMode!=kNone) */
419 else
420 {
421 calibConv = 1./fGeomCam->GetPixRatio(pixidx);
422 }
423
424 calibConv /= calibQE;
425
426 if (calibConv != 0. && calibQE != 0.)
427 {
428 // Now doing:
429 // calibConvVar = calibConvVar/(calibConv*calibConv) + calibQEVar/(calibQE*calibQE);
430 // calibConvVar *= (calibConv*calibConv);
431 calibConvVar += calibQEVar*(calibConv*calibConv)/(calibQE*calibQE);
432 }
433
434 return kTRUE;
435}
436
437Int_t MCalibrateData::Calibrate(Bool_t data, Bool_t pedestal) const
438{
439 if (!data && !pedestal)
440 return kTRUE;
441
442 const UInt_t npix = fSignals->GetSize();
443 const Float_t slices = fSignals->GetNumUsedHiGainFADCSlices();
444 const Float_t sqrtslices = TMath::Sqrt(slices);
445
446 Float_t hiloconv;
447 Float_t hiloconverr;
448 Float_t calibConv;
449 Float_t calibConvErr;
450 Float_t calibFFactor;
451
452 UInt_t skip = 0;
453 for (UInt_t pixidx=0; pixidx<npix; pixidx++)
454 {
455 if (!GetConversionFactor(pixidx, hiloconv, hiloconverr,
456 calibConv, calibConvErr, calibFFactor))
457 {
458 skip++;
459 continue;
460 }
461
462 if (data)
463 {
464 const MExtractedSignalPix &sig = (*fSignals)[pixidx];
465
466 Float_t signal = 0;
467 Float_t signalErr = 0.;
468
469 if (sig.IsLoGainUsed())
470 {
471 signal = sig.GetExtractedSignalLoGain()*hiloconv;
472 signalErr = signal*hiloconverr;
473 }
474 else
475 {
476 if (sig.GetExtractedSignalHiGain() <= 9999.)
477 signal = sig.GetExtractedSignalHiGain();
478 }
479
480 const Float_t nphot = signal*calibConv;
481 const Float_t nphotErr = calibFFactor*TMath::Sqrt(TMath::Abs(nphot));
482
483 //
484 // The following part is the commented first version of the error calculation
485 // Contact Markus Gaug for questions (or wait for the next documentation update...)
486 //
487 /*
488 nphotErr = signal > 0 ? signalErr*signalErr / (signal * signal) : 0.
489 + calibConv > 0 ? calibConvVar / (calibConv * calibConv ) : 0.;
490 nphotErr = TMath::Sqrt(nphotErr) * nphot;
491 */
492
493 MCerPhotPix *cpix = fCerPhotEvt->AddPixel(pixidx, nphot, nphotErr);
494
495 if (sig.GetNumHiGainSaturated() > 0)
496 cpix->SetPixelHGSaturated();
497
498 if (sig.GetNumLoGainSaturated() > 0)
499 cpix->SetPixelSaturated();
500 }
501 if (pedestal)
502 {
503 const MPedestalPix &ped = (*fPedestal)[pixidx];
504
505 // pedestals/(used FADC slices) in [ADC] counts
506 const Float_t pedes = ped.GetPedestal() * slices;
507 const Float_t pedrms = ped.GetPedestalRms() * sqrtslices;
508
509 //
510 // pedestals/(used FADC slices) in [number of photons]
511 //
512 const Float_t pedphot = pedes * calibConv;
513 const Float_t pedphotrms = pedrms * calibConv;
514
515 (*fPedPhot)[pixidx].Set(pedphot, pedphotrms);
516 }
517 }
518
519 if (skip>npix*0.9)
520 {
521 *fLog << warn << "WARNING - GetConversionFactor has skipped more than 90% of the pixels... skip." << endl;
522 return kCONTINUE;
523 }
524
525 if (pedestal)
526 fPedPhot->SetReadyToSave();
527
528 if (data)
529 {
530 fCerPhotEvt->FixSize();
531 fCerPhotEvt->SetReadyToSave();
532 }
533 return kTRUE;
534}
535
536// --------------------------------------------------------------------------
537//
538// Apply the calibration factors to the extracted signal according to the
539// selected calibration method
540//
541Int_t MCalibrateData::Process()
542{
543 /*
544 if (fCalibrations->GetNumPixels() != (UInt_t)fSignals->GetSize())
545 {
546 // FIXME: MExtractedSignal must be of variable size -
547 // like MCerPhotEvt - because we must be able
548 // to reduce size by zero supression
549 // For the moment this check could be done in ReInit...
550 *fLog << err << "MExtractedSignal and MCalibrationCam have different sizes... abort." << endl;
551 return kFALSE;
552 }
553 */
554
555 return Calibrate(fCalibrationMode!=kSkip, TestPedestalFlag(kEvent));
556}
557
558// --------------------------------------------------------------------------
559//
560// Implementation of SavePrimitive. Used to write the call to a constructor
561// to a macro. In the original root implementation it is used to write
562// gui elements to a macro-file.
563//
564void MCalibrateData::StreamPrimitive(ofstream &out) const
565{
566 out << " " << ClassName() << " " << GetUniqueName() << "(\"";
567 out << "\"" << fName << "\", \"" << fTitle << "\");" << endl;
568
569 if (TestPedestalFlag(kEvent))
570 out << " " << GetUniqueName() << ".EnablePedestalType(MCalibrateData::kEvent)" << endl;
571 if (TestPedestalFlag(kRun))
572 out << " " << GetUniqueName() << ".EnablePedestalType(MCalibrateData::kRun)" << endl;
573
574 if (fCalibrationMode != kDefault)
575 {
576 out << " " << GetUniqueName() << ".SetCalibrationMode(MCalibrateData::";
577 switch (fCalibrationMode)
578 {
579 case kSkip: out << "kSkip"; break;
580 case kNone: out << "kNone"; break;
581 case kFlatCharge: out << "kFlatCharge"; break;
582 case kBlindPixel: out << "kBlindPixel"; break;
583 case kFfactor: out << "kFfactor"; break;
584 case kPinDiode: out << "kPinDiode"; break;
585 case kCombined: out << "kCombined"; break;
586 case kDummy: out << "kDummy"; break;
587 default: out << (int)fCalibrationMode; break;
588 }
589 out << ")" << endl;
590 }
591
592 if (fNamePedADCContainer != fgNamePedADCContainer)
593 {
594 out << " " << GetUniqueName() << ".SetNamePedADCContainer(";
595 out << fNamePedADCContainer.Data() << ");" << endl;
596 }
597
598 if (fNamePedPhotContainer != fgNamePedPhotContainer)
599 {
600 out << " " << GetUniqueName() << ".SetNamePedPhotContainer(";
601 out << fNamePedPhotContainer.Data() << ");" << endl;
602 }
603}
604
605// --------------------------------------------------------------------------
606//
607// Read the setup from a TEnv, eg:
608// MJPedestal.MCalibrateDate.PedestalFlag: no,run,event
609// MJPedestal.MCalibrateDate.CalibrationMode: skip,none,flatcharge,blindpixel,ffactor,pindiode,combined,dummy,default
610//
611Int_t MCalibrateData::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
612{
613 Bool_t rc = kFALSE;
614 if (IsEnvDefined(env, prefix, "PedestalFlag", print))
615 {
616 rc = kTRUE;
617 TString s = GetEnvValue(env, prefix, "PedestalFlag", "");
618 s.ToLower();
619 if (s.BeginsWith("no"))
620 SetPedestalFlag(kNo);
621 if (s.BeginsWith("run"))
622 SetPedestalFlag(kRun);
623 if (s.BeginsWith("event"))
624 SetPedestalFlag(kEvent);
625 }
626
627 if (IsEnvDefined(env, prefix, "CalibrationMode", print))
628 {
629 rc = kTRUE;
630 TString s = GetEnvValue(env, prefix, "CalibrationMode", "");
631 s.ToLower();
632 if (s.BeginsWith("skip"))
633 SetCalibrationMode(kSkip);
634 if (s.BeginsWith("none"))
635 SetCalibrationMode(kNone);
636 if (s.BeginsWith("flatcharge"))
637 SetCalibrationMode(kFlatCharge);
638 if (s.BeginsWith("blindpixel"))
639 SetCalibrationMode(kBlindPixel);
640 if (s.BeginsWith("ffactor"))
641 SetCalibrationMode(kFfactor);
642 if (s.BeginsWith("pindiode"))
643 SetCalibrationMode(kPinDiode);
644 if (s.BeginsWith("combined"))
645 SetCalibrationMode(kCombined);
646 if (s.BeginsWith("dummy"))
647 SetCalibrationMode(kDummy);
648 if (s.BeginsWith("default"))
649 SetCalibrationMode();
650 }
651
652 return rc;
653}
Note: See TracBrowser for help on using the repository browser.