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

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