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

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