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

Last change on this file since 4595 was 4590, checked in by gaug, 21 years ago
*** empty log message ***
File size: 17.5 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!
24! Copyright: MAGIC Software Development, 2000-2004
25!
26!
27\* ======================================================================== */
28
29///////////////////////////////////////////////////////////////////////////////////
30//
31// MCalibrateData
32//
33// This task takes the integrated charge from MExtractedSignal and applies
34// the calibration constants from MCalibrationCam to convert the summed FADC
35// slices into photons. The number of photons obtained is stored in MCerPhotEvt.
36// Optionally, the calibration of pedestals from an MPedestalCam container into
37// an MPedPhotCam container can be chosen with the member functions
38// SetPedestalType(). Default is 'kRun', i.e. calibration of pedestals from a
39// dedicated pedestal run.
40// In case, the chosen pedestal type is kRun or kEvt, in ReInit() the MPedPhotCam
41// container is filled using the information from MPedestalCam, MExtractedSignalCam,
42// MCalibrationChargeCam and MCalibrationQECam
43//
44// Selection of different calibration methods is allowed through the
45// SetCalibrationMode() member function (default: kFfactor)
46//
47// The calibration modes which exclude non-valid pixels are the following:
48//
49// kFfactor: calibrates using the F-Factor method
50// kBlindpixel: calibrates using the BlindPixel method
51// kBlindpixel: calibrates using the BlindPixel method
52// kFlatCharge: perform a charge flat-flatfielding. Outer pixels are area-corrected.
53// kDummy: calibrates with fixed conversion factors of 1 and errors of 0.
54//
55// The calibration modes which include all pixels regardless of their validity is:
56//
57// kNone: calibrates with fixed conversion factors of 1 and errors of 0.
58//
59// Use the kDummy and kNone methods ONLY FOR DEBUGGING!
60//
61// Input Containers:
62// MPedestalCam
63// MExtractedSignalCam
64// MCalibrationChargeCam
65// MCalibrationQECam
66//
67// Output Containers:
68// MPedPhotCam
69// MCerPhotEvt
70//
71// See also: MJCalibration, MJPedestal, MJExtractSignal, MJExtractCalibTest
72//
73//////////////////////////////////////////////////////////////////////////////
74#include "MCalibrateData.h"
75
76#include "MLog.h"
77#include "MLogManip.h"
78
79#include "MParList.h"
80#include "MH.h"
81
82#include "MGeomCam.h"
83
84#include "MPedestalCam.h"
85#include "MPedestalPix.h"
86
87#include "MCalibrationChargeCam.h"
88#include "MCalibrationChargePix.h"
89
90#include "MCalibrationQECam.h"
91#include "MCalibrationQEPix.h"
92
93#include "MExtractedSignalCam.h"
94#include "MExtractedSignalPix.h"
95
96#include "MPedPhotCam.h"
97#include "MPedPhotPix.h"
98
99#include "MBadPixelsCam.h"
100#include "MBadPixelsPix.h"
101
102#include "MCerPhotEvt.h"
103
104ClassImp(MCalibrateData);
105
106using namespace std;
107// --------------------------------------------------------------------------
108//
109// Default constructor.
110//
111MCalibrateData::MCalibrateData(CalibrationMode_t calmode,const char *name, const char *title)
112 : fGeomCam(NULL), fPedestal(NULL), fBadPixels(NULL), fCalibrations(NULL), fSignals(NULL),
113 fPedPhot(NULL), fCerPhotEvt(NULL), fCalibrationMode(calmode), fFlags(kRun), fNamePedADCRunContainer("MPedestalCam"), fNamePedADCEventContainer("MPedestalCamFromData"), fNamePedPhotRunContainer("MPedPhotCam"), fNamePedPhotEventContainer("MPedPhotCamFromData")
114{
115 fName = name ? name : "MCalibrateData";
116 fTitle = title ? title : "Task to calculate the number of photons in one event";
117}
118
119// --------------------------------------------------------------------------
120//
121// The PreProcess searches for the following input containers:
122// - MGeomCam
123// - MPedestalCam
124// - MCalibrationChargeCam
125// - MCalibrationQECam
126// - MExtractedSignalCam
127//
128// The following output containers are also searched and created if
129// they were not found:
130//
131// - MPedPhotCam
132// - MCerPhotEvt
133//
134Int_t MCalibrateData::PreProcess(MParList *pList)
135{
136
137 // input containers
138
139
140 if (TestFlag(kRun))
141 {
142 fPedestal = (MPedestalCam*)pList->FindObject(fNamePedADCRunContainer, AddSerialNumber("MPedestalCam"));
143 if (!fPedestal)
144 {
145 *fLog << err << "container 'MPedestalCam'"
146 << " not found ... aborting" << endl;
147 return kFALSE;
148 }
149 }
150
151
152 if (TestFlag(kEvent))
153 {
154 fPedestalFromData = (MPedestalCam*)pList->FindObject(fNamePedADCEventContainer, AddSerialNumber("MPedestalCam"));
155 if (!fPedestalFromData)
156 {
157 *fLog << err << "container 'MPedestalCamFromData'"
158 << " not found ... aborting" << endl;
159 return kFALSE;
160 }
161 }
162
163
164 fSignals = (MExtractedSignalCam*)pList->FindObject(AddSerialNumber("MExtractedSignalCam"));
165 if (!fSignals)
166 {
167 *fLog << err << AddSerialNumber("MExtractedSignalCam") << " not found ... aborting" << endl;
168 return kFALSE;
169 }
170
171 fBadPixels = (MBadPixelsCam*)pList->FindObject(AddSerialNumber("MBadPixelsCam"));
172
173 if (!fBadPixels)
174 {
175 *fLog << err << AddSerialNumber("MBadPixelsCam") << " not found ... aborting" << endl;
176 return kFALSE;
177 }
178
179 if (fCalibrationMode>kNone)
180 {
181 fCalibrations = (MCalibrationChargeCam*)pList->FindObject(AddSerialNumber("MCalibrationChargeCam"));
182 if (!fCalibrations)
183 {
184 *fLog << err << AddSerialNumber("MCalibrationChargeCam") << " not found ... aborting." << endl;
185 return kFALSE;
186 }
187 fQEs = (MCalibrationQECam*)pList->FindObject(AddSerialNumber("MCalibrationQECam"));
188 if (!fQEs)
189 {
190 *fLog << err << AddSerialNumber("MCalibrationQECam") << " not found ... aborting." << endl;
191 return kFALSE;
192 }
193 }
194
195 // output containers
196
197 if (TestFlag(kRun))
198 {
199 fPedPhot = (MPedPhotCam*)pList->FindCreateObj(AddSerialNumber("MPedPhotCam"), fNamePedPhotRunContainer);
200 if (!fPedPhot)
201 {
202 *fLog << err << "container 'MPedPhotCam'"
203 << " not found ... aborting" << endl;
204 return kFALSE;
205 }
206 }
207
208
209 if (TestFlag(kEvent))
210 {
211 fPedPhotFromData = (MPedPhotCam*)pList->FindCreateObj(AddSerialNumber("MPedPhotCam"), fNamePedPhotEventContainer);
212 if (!fPedPhotFromData)
213 {
214 *fLog << err << "container 'MPedPhotCamFromData'"
215 << " not found ... aborting" << endl;
216 return kFALSE;
217 }
218 }
219
220 fCerPhotEvt = (MCerPhotEvt*)pList->FindCreateObj(AddSerialNumber("MCerPhotEvt"));
221 if (!fCerPhotEvt)
222 return kFALSE;
223
224 return kTRUE;
225}
226
227// --------------------------------------------------------------------------
228//
229// Check for validity of the selected calibration method, switch to a
230// different one in case of need
231//
232// fill the MPedPhotCam container using the information from MPedestalCam,
233// MExtractedSignalCam and MCalibrationCam
234//
235//
236Bool_t MCalibrateData::ReInit(MParList *pList)
237{
238
239 if(fCalibrationMode == kBlindPixel && !fQEs->IsBlindPixelMethodValid())
240 {
241 *fLog << warn << GetDescriptor() << "Warning: Blind pixel calibration method not valid, switching to F-factor method" << endl;
242 fCalibrationMode = kFfactor;
243 }
244
245 if(fCalibrationMode == kPinDiode && !fQEs->IsPINDiodeMethodValid())
246 {
247 *fLog << warn << GetDescriptor() << "Warning: PIN diode calibration method not valid, switching to F-factor method" << endl;
248 fCalibrationMode = kFfactor;
249 }
250
251
252 if (TestFlag(kRun))
253 {
254 if (!CalibratePedestal(fPedestal,fPedPhot)) return kFALSE;
255 }
256
257 return kTRUE;
258}
259
260
261
262// --------------------------------------------------------------------------
263//
264// Calibrate the Pedestal values
265//
266//
267Bool_t MCalibrateData::CalibratePedestal(MPedestalCam *pedADCCam, MPedPhotCam *pedPhotCam)
268{
269 //---------------------------------------------
270 // fill MPedPhot container using the informations from
271 // MPedestalCam, MExtractedSignalCam and MCalibrationCam
272
273 fNumUsedHiGainFADCSlices = fSignals->GetNumUsedHiGainFADCSlices();
274
275 // is pixid equal to pixidx ?
276 if ( (Int_t)(pedADCCam->GetSize()) != fSignals->GetSize())
277 {
278 *fLog << err << "MCalibrateData::ReInit(); sizes of MPedestalCam and MCalibrationCam are different"
279 << endl;
280 }
281
282
283
284 for (Int_t pixid=0; pixid<pedADCCam->GetSize(); pixid++)
285 {
286 const MPedestalPix &ped = (*pedADCCam)[pixid];
287
288 // pedestals/(used FADC slices) in [ADC] counts
289 Float_t pedes = ped.GetPedestal() * fNumUsedHiGainFADCSlices;
290 Float_t pedrms = ped.GetPedestalRms() * sqrt(fNumUsedHiGainFADCSlices);
291
292 //----------------------------------
293 // get phe/ADC conversion factor
294
295 Float_t hiloconv;
296 Float_t hiloconverr;
297 Float_t calibConv;
298 Float_t calibConvVar;
299 Float_t calibFFactor;
300
301 if ( !GetConversionFactor(pixid, hiloconv, hiloconverr,
302 calibConv, calibConvVar, calibFFactor ))
303 continue;
304
305 //----------------------------------
306
307 // pedestals/(used FADC slices) in [number of photons]
308 Float_t pedphot = pedes * calibConv;
309 Float_t pedphotrms = pedrms * calibConv;
310
311 (*pedPhotCam)[pixid].Set(pedphot, pedphotrms);
312
313 }
314
315 //---------------------------------------------
316
317 pedPhotCam->SetReadyToSave();
318
319 return kTRUE;
320}
321
322
323
324
325
326
327
328// --------------------------------------------------------------------------
329//
330// Get conversion factor and its error from MCalibrationCam
331//
332//
333Bool_t MCalibrateData::GetConversionFactor(UInt_t pixidx,
334 Float_t &hiloconv, Float_t &hiloconverr,
335 Float_t &calibConv, Float_t &calibConvVar, Float_t &calibFFactor)
336{
337
338 //
339 // For the moment, we use only a dummy zenith for the calibration:
340 //
341 const Float_t zenith = -1.;
342
343 hiloconv = 1.;
344 hiloconverr = 0.;
345 calibConv = 1.;
346 calibConvVar = 0.;
347 calibFFactor = 0.;
348 Float_t calibQE = 1.;
349 Float_t calibQEVar = 0.;
350
351 Float_t avMean = 1.;
352 Float_t avMeanRelVar = 0.;
353
354
355
356 if (fCalibrationMode == kFlatCharge)
357 {
358 MCalibrationChargePix &avpix = (MCalibrationChargePix&)fCalibrations->GetAverageArea(0);
359 avMean = avpix.GetMean();
360 avMeanRelVar = avpix.GetMeanRelVar();
361 }
362
363
364 if(fCalibrationMode!=kNone)
365 {
366
367 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCalibrations)[pixidx];
368
369 hiloconv = pix.GetConversionHiLo ();
370 hiloconverr= pix.GetConversionHiLoErr();
371
372 if (fBadPixels)
373 {
374 MBadPixelsPix &bad = (*fBadPixels)[pixidx];
375 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
376 return kFALSE;
377 }
378
379 calibConv = pix.GetMeanConvFADC2Phe();
380 calibConvVar = pix.GetMeanConvFADC2PheVar();
381 calibFFactor = pix.GetMeanFFactorFADC2Phot();
382
383 MCalibrationQEPix &qe = (MCalibrationQEPix&) (*fQEs)[pixidx];
384
385 switch(fCalibrationMode)
386 {
387 case kFlatCharge:
388 calibConv = avMean / pix.GetMean() / fGeomCam->GetPixRatio(pixidx) ;
389 calibConvVar = (avMeanRelVar + pix.GetMeanRelVar()) * calibConv * calibConv;
390 if (pix.IsFFactorMethodValid())
391 {
392 const Float_t convmin1 = qe.GetQECascadesFFactor(zenith)/pix.GetMeanConvFADC2Phe();
393 if (convmin1 > 0)
394 calibFFactor *= TMath::Sqrt(convmin1);
395 else
396 calibFFactor = -1.;
397 }
398 break;
399 case kBlindPixel:
400 if (qe.IsBlindPixelMethodValid())
401 {
402 calibQE = qe.GetQECascadesBlindPixel ( zenith );
403 calibQEVar = qe.GetQECascadesBlindPixelVar( zenith );
404 }
405 else
406 return kFALSE;
407 break;
408 case kPinDiode:
409 if (qe.IsPINDiodeMethodValid())
410 {
411 calibQE = qe.GetQECascadesPINDiode ( zenith );
412 calibQEVar = qe.GetQECascadesPINDiodeVar( zenith );
413 }
414 else
415 return kFALSE;
416 break;
417 case kFfactor:
418 if (pix.IsFFactorMethodValid())
419 {
420 calibQE = qe.GetQECascadesFFactor ( zenith );
421 calibQEVar = qe.GetQECascadesFFactorVar( zenith );
422 }
423 else
424 return kFALSE;
425 break;
426 case kCombined:
427 if (qe.IsCombinedMethodValid())
428 {
429 calibQE = qe.GetQECascadesCombined ( zenith );
430 calibQEVar = qe.GetQECascadesCombinedVar( zenith );
431 }
432 else
433 return kFALSE;
434 break;
435 case kDummy:
436 hiloconv = 1.;
437 hiloconverr = 0.;
438 calibQE = 1.;
439 calibQEVar = 0.;
440 break;
441
442 } /* switch calibration mode */
443 } /* if(fCalibrationMode!=kNone) */
444 else
445 {
446 hiloconv = 1.;
447 hiloconverr = 0.;
448 calibConv = 1./fGeomCam->GetPixRatio(pixidx);
449 calibConvVar = 0.;
450 calibFFactor = 0.;
451 calibQE = 1.;
452 calibQEVar = 0.;
453 }
454
455 calibConv /= calibQE;
456 calibConvVar /= calibQE;
457
458
459 return kTRUE;
460}
461
462
463// --------------------------------------------------------------------------
464//
465// Apply the calibration factors to the extracted signal according to the
466// selected calibration method
467//
468Int_t MCalibrateData::Process()
469{
470 /*
471 if (fCalibrations->GetNumPixels() != (UInt_t)fSignals->GetSize())
472 {
473 // FIXME: MExtractedSignal must be of variable size -
474 // like MCerPhotEvt - because we must be able
475 // to reduce size by zero supression
476 // For the moment this check could be done in ReInit...
477 *fLog << err << "MExtractedSignal and MCalibrationCam have different sizes... abort." << endl;
478 return kFALSE;
479 }
480 */
481
482 UInt_t npix = fSignals->GetSize();
483
484 Float_t hiloconv;
485 Float_t hiloconverr;
486 Float_t calibrationConversionFactor;
487 Float_t calibrationConversionFactorErr;
488 Float_t calibFFactor;
489
490 for (UInt_t pixidx=0; pixidx<npix; pixidx++)
491 {
492 if ( !GetConversionFactor(pixidx, hiloconv, hiloconverr,
493 calibrationConversionFactor, calibrationConversionFactorErr, calibFFactor) )
494 continue;
495
496 MExtractedSignalPix &sig = (*fSignals)[pixidx];
497
498 Float_t signal;
499 Float_t signalErr = 0.;
500 Float_t nphot,nphotErr;
501
502 if (sig.IsLoGainUsed())
503 {
504 signal = sig.GetExtractedSignalLoGain()*hiloconv;
505 signalErr = signal*hiloconverr;
506 }
507 else
508 {
509 if (sig.GetExtractedSignalHiGain() > 9999.)
510 {
511 signal = 0.;
512 signalErr = 0.;
513 }
514 else
515 signal = sig.GetExtractedSignalHiGain();
516 }
517
518 nphot = signal*calibrationConversionFactor;
519 nphotErr = calibFFactor*TMath::Sqrt(TMath::Abs(nphot));
520 // nphotErr = signal*calibrationConversionFactorErr
521 // *signal*calibrationConversionFactorErr
522 // +signalErr*calibrationConversionFactor
523 // *signalErr*calibrationConversionFactor;
524 // nphotErr = TMath::Sqrt(nphotErr);
525
526 MCerPhotPix *cpix = fCerPhotEvt->AddPixel(pixidx, nphot, nphotErr);
527
528 if (sig.GetNumHiGainSaturated() > 0)
529 cpix->SetPixelHGSaturated();
530
531 if (sig.GetNumLoGainSaturated() > 0)
532 cpix->SetPixelSaturated();
533 }
534
535 fCerPhotEvt->FixSize();
536 fCerPhotEvt->SetReadyToSave();
537
538
539
540 if(TestFlag(kEvent))
541 {
542 if (!CalibratePedestal(fPedestalFromData,fPedPhotFromData)) return kFALSE;
543 }
544
545 /*
546 //---------------------------------------------
547 // fill MPedPhot(FromData) container using the informations from
548 // MPedestalCam, MExtractedSignalCam and MCalibrationCam
549
550 fNumUsedHiGainFADCSlices = fSignals->GetNumUsedHiGainFADCSlices();
551
552 // is pixid equal to pixidx ?
553 if ( (Int_t)(fPedestal->GetSize()) != fSignals->GetSize())
554 {
555 *fLog << err << "MCalibrateData::ReInit(); sizes of MPedestalCam and MCalibrationCam are different"
556 << endl;
557 } */
558
559 /*
560 *fLog << all << "MCalibrateData::ReInit(); fill MPedPhotCam container"
561 << endl;
562 *fLog << all << " fNumUsedHiGainADCSlices = "
563 << fNumUsedHiGainFADCSlices << endl;
564 *fLog << all << " pixid, calibrationConversionFactor, ped, pedRMS, pedphot, pedphotRMS :"
565 << endl;
566 */ /*
567 for (Int_t pixid=0; pixid<fPedestalFromData->GetSize(); pixid++)
568 {
569 const MPedestalPix &ped = (*fPedestalFromData)[pixid];
570
571 // pedestals/(used FADC slices) in [ADC] counts
572 Float_t pedes = ped.GetPedestal() * fNumUsedHiGainFADCSlices;
573 Float_t pedrms = ped.GetPedestalRms() * sqrt(fNumUsedHiGainFADCSlices);
574
575 //----------------------------------
576 // get phe/ADC conversion factor
577
578 Float_t hiloconv;
579 Float_t hiloconverr;
580 Float_t calibConv;
581 Float_t calibFFactor;
582
583 if ( !GetConversionFactor(pixid, hiloconv, hiloconverr,
584 calibConv, calibFFactor ))
585 continue;
586
587 //----------------------------------
588
589 // pedestals/(used FADC slices) in [number of photons]
590 Float_t pedphot = pedes * calibConv;
591 Float_t pedphotrms = pedrms * calibConv;
592
593 (*fPedPhotFromData)[pixid].Set(pedphot, pedphotrms);
594
595 // *fLog << all << pixid << ", " << calibConv << ", "
596 // << ped.GetPedestal() << ", " << ped.GetPedestalRms() << ", "
597 // << pedphot << ", " << pedphotrms << endl;
598 }
599
600 //---------------------------------------------
601
602 fPedPhotFromData->SetReadyToSave();
603
604 }
605 */
606
607 return kTRUE;
608}
Note: See TracBrowser for help on using the repository browser.