source: trunk/MagicSoft/Mars/mcalib/MCalibrate.cc@ 5020

Last change on this file since 5020 was 4934, checked in by gaug, 20 years ago
*** empty log message ***
File size: 13.1 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): Markus Gaug 04/2004 <mailto:markus@ifae.es>
21!
22! Copyright: MAGIC Software Development, 2000-2004
23!
24!
25\* ======================================================================== */
26
27//////////////////////////////////////////////////////////////////////////////
28//
29// MCalibrate
30//
31// This task takes the integrated charge from MExtractedSignal and apply
32// the calibration constants from MCalibraitionCam to convert the summed FADC
33// slices to photons. The number of photons obtained is stored in MCerPhotEvt.
34//
35// Selection of different calibration methods is possible through the
36// SetCalibrationMode member function
37//
38// The calibration modes which exclude non-valid pixels are the following:
39//
40// kFfactor: calibrates using the F-Factor method
41// kBlindpixel: calibrates using the BlindPixel method
42// kBlindpixel: calibrates using the BlindPixel method
43// kFlatCharge: perform a charge flat-flatfielding. Outer pixels are area-corrected.
44// kDummy: calibrates with fixed conversion factors of 1 and errors of 0.
45//
46// The calibration modes which include all pixels regardless of their validity is:
47//
48// kNone: calibrates with fixed conversion factors of 1 and errors of 0.
49//
50// Use the kDummy and kNone methods ONLY FOR DEBUGGING!
51//
52// Input Containers:
53// MExtractedSingal
54// MCalibrationChargeCam
55//
56// Output Containers:
57// MCerPhotEvt
58//
59//////////////////////////////////////////////////////////////////////////////
60#include "MCalibrate.h"
61
62#include "MLog.h"
63#include "MLogManip.h"
64
65#include "MParList.h"
66
67#include "MGeomCam.h"
68
69#include "MCalibrationChargeCam.h"
70#include "MCalibrationChargePix.h"
71
72#include "MCalibrationQECam.h"
73#include "MCalibrationQEPix.h"
74
75#include "MExtractedSignalCam.h"
76#include "MExtractedSignalPix.h"
77
78#include "MBadPixelsCam.h"
79#include "MBadPixelsPix.h"
80
81#include "MCerPhotEvt.h"
82
83ClassImp(MCalibrate);
84
85using namespace std;
86// --------------------------------------------------------------------------
87//
88// Default constructor.
89//
90MCalibrate::MCalibrate(CalibrationMode_t calmode,const char *name, const char *title)
91 : fGeomCam(NULL), fCalibrations(NULL), fQEs(NULL), fBadPixels(NULL), fSignals(NULL),
92 fCerPhotEvt(NULL), fCalibrationMode(calmode)
93{
94 fName = name ? name : "MCalibrate";
95 fTitle = title ? title : "Task to calculate the number of photons in one event";
96}
97
98// --------------------------------------------------------------------------
99//
100// The PreProcess searches for the following input containers:
101// - MGeomCam
102// - MCalibrationChargeCam
103// - MExtractedSignalCam
104// - MBadPixelsCam
105//
106// The following output containers are also searched and created if
107// they were not found:
108//
109// - MCerPhotEvt
110//
111Int_t MCalibrate::PreProcess(MParList *pList)
112{
113
114 fSignals = (MExtractedSignalCam*)pList->FindObject(AddSerialNumber("MExtractedSignalCam"));
115
116 if (!fSignals)
117 {
118 *fLog << err << AddSerialNumber("MExtractedSignalCam") << " not found ... aborting" << endl;
119 return kFALSE;
120 }
121
122 fBadPixels = (MBadPixelsCam*)pList->FindObject(AddSerialNumber("MBadPixelsCam"));
123 if (!fBadPixels)
124 *fLog << warn << AddSerialNumber("MBadPixelsCam") << " not found ... no action" << endl;
125
126 if(fCalibrationMode>kNone)
127 {
128
129 fCalibrations = (MCalibrationChargeCam*)pList->FindObject(AddSerialNumber("MCalibrationChargeCam"));
130 if (!fCalibrations)
131 {
132 *fLog << err << AddSerialNumber("MCalibrationChargeCam") << " not found ... aborting." << endl;
133 return kFALSE;
134 }
135
136 fQEs = (MCalibrationQECam*)pList->FindObject(AddSerialNumber("MCalibrationQECam"));
137 if (!fQEs)
138 {
139 *fLog << err << AddSerialNumber("MCalibrationQECam") << " not found ... aborting." << endl;
140 return kFALSE;
141 }
142 }
143
144 fCerPhotEvt = (MCerPhotEvt*)pList->FindCreateObj(AddSerialNumber("MCerPhotEvt"));
145 if (!fCerPhotEvt)
146 return kFALSE;
147
148 fGeomCam = (MGeomCam*)pList->FindCreateObj(AddSerialNumber("MGeomCam"));
149 if (!fGeomCam)
150 return kFALSE;
151
152 return kTRUE;
153}
154
155// --------------------------------------------------------------------------
156//
157// Check for validity of the selected calibration method, switch to a
158// different one in case of need
159//
160Bool_t MCalibrate::ReInit(MParList *pList)
161{
162
163 if(fCalibrationMode == kBlindPixel && !fQEs->IsBlindPixelMethodValid())
164 {
165 *fLog << warn << GetDescriptor()
166 << "Warning: Blind pixel calibration method not valid, switching to F-factor method" << endl;
167 fCalibrationMode = kFfactor;
168 }
169
170 if(fCalibrationMode == kPinDiode && !fQEs->IsPINDiodeMethodValid())
171 {
172 *fLog << warn << GetDescriptor()
173 << "Warning: PIN diode calibration method not valid, switching to F-factor method" << endl;
174 fCalibrationMode = kFfactor;
175 }
176
177 if(fCalibrationMode == kCombined && !fQEs->IsCombinedMethodValid())
178 {
179 *fLog << warn << GetDescriptor()
180 << "Warning: Combined calibration method not valid, switching to F-factor method" << endl;
181 fCalibrationMode = kFfactor;
182 }
183
184
185 if (fCalibrationMode == kFlatCharge)
186 {
187 fGeomCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
188 if (!fGeomCam)
189 {
190 *fLog << err << AddSerialNumber("MGeomCam") << " not found ... aborting." << endl;
191 return kFALSE;
192 }
193 }
194
195
196 switch(fCalibrationMode)
197 {
198 case kBlindPixel:
199 break;
200 case kFfactor:
201 break;
202 case kPinDiode:
203 *fLog << err << GetDescriptor()
204 << ": PIN Diode Calibration mode not yet available " << endl;
205 return kFALSE;
206 break;
207 case kCombined:
208 *fLog << err << GetDescriptor()
209 << ": Combined Calibration mode not yet available " << endl;
210 return kFALSE;
211 break;
212 case kFlatCharge:
213 *fLog << warn << GetDescriptor()
214 << ": WARNING: Flat-fielding charges - only for Keiichi!!" << endl;
215 break;
216 case kDummy:
217 *fLog << warn << GetDescriptor()
218 << ": WARNING: Dummy calibration, no calibration applied!!" << endl;
219 break;
220 case kNone:
221 *fLog << warn << GetDescriptor()
222 << ": WARNING: No calibration applied!!" << endl;
223 break;
224 default:
225 *fLog << warn << GetDescriptor()
226 << ": WARNING: Calibration mode value ("
227 <<fCalibrationMode<<") not known" << endl;
228 return kFALSE;
229 }
230
231 return kTRUE;
232}
233// --------------------------------------------------------------------------
234//
235// Apply the calibration factors to the extracted signal according to the
236// selected calibration method
237//
238Int_t MCalibrate::Process()
239{
240
241 //
242 // For the moment, we use only a dummy zenith for the calibration:
243 //
244 const Float_t zenith = 0;
245
246 /*
247 if (fCalibrations->GetNumPixels() != (UInt_t)fSignals->GetSize())
248 {
249 // FIXME: MExtractedSignal must be of variable size -
250 // like MCerPhotEvt - because we must be able
251 // to reduce size by zero supression
252 // For the moment this check could be done in ReInit...
253 *fLog << err << "MExtractedSignal and MCalibrationCam have different sizes... abort." << endl;
254 return kFALSE;
255 }
256 */
257
258 UInt_t npix = fSignals->GetSize();
259
260 Float_t hiloconv = 0.;
261 Float_t hiloconverr = 0.;
262 Float_t calibConv = 0.;
263 Float_t calibConvVar = 0.;
264 Float_t calibFFactor = 0.;
265 Float_t calibQE = 1.;
266 Float_t calibQEVar = 0.;
267 Float_t avMean = 1.;
268 Float_t avMeanRelVar = 0.;
269
270
271 if (fCalibrationMode == kFlatCharge)
272 {
273 MCalibrationChargePix &avpix = (MCalibrationChargePix&)fCalibrations->GetAverageArea(0);
274 avMean = avpix.GetMean();
275 avMeanRelVar = avpix.GetMeanRelVar();
276 }
277
278
279 for (UInt_t pixidx=0; pixidx<npix; pixidx++)
280 {
281
282 if(fCalibrationMode!=kNone)
283 {
284
285 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCalibrations)[pixidx];
286
287 hiloconv = pix.GetConversionHiLo ();
288 hiloconverr= pix.GetConversionHiLoErr();
289
290 if (fBadPixels)
291 {
292 MBadPixelsPix &bad = (*fBadPixels)[pixidx];
293 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
294 continue;
295 }
296
297 calibConv = pix.GetMeanConvFADC2Phe();
298 calibConvVar = pix.GetMeanConvFADC2PheVar();
299 calibFFactor = pix.GetMeanFFactorFADC2Phot();
300
301 MCalibrationQEPix &qe = (MCalibrationQEPix&) (*fQEs)[pixidx];
302
303 switch(fCalibrationMode)
304 {
305 case kFlatCharge:
306 calibConv = avMean / pix.GetMean() / fGeomCam->GetPixRatio(pixidx) ;
307 calibConvVar = (avMeanRelVar + pix.GetMeanRelVar()) * calibConv * calibConv;
308 if (pix.IsFFactorMethodValid())
309 {
310 const Float_t convmin1 = qe.GetQECascadesFFactor(zenith)/pix.GetMeanConvFADC2Phe();
311 if (convmin1 > 0)
312 calibFFactor *= TMath::Sqrt(convmin1);
313 else
314 calibFFactor = -1.;
315 }
316 break;
317 case kBlindPixel:
318 if (qe.IsBlindPixelMethodValid())
319 {
320 calibQE = qe.GetQECascadesBlindPixel ( zenith );
321 calibQEVar = qe.GetQECascadesBlindPixelVar( zenith );
322 }
323 else
324 continue;
325 break;
326 case kPinDiode:
327 if (qe.IsPINDiodeMethodValid())
328 {
329 calibQE = qe.GetQECascadesPINDiode ( zenith );
330 calibQEVar = qe.GetQECascadesPINDiodeVar( zenith );
331 }
332 else
333 continue;
334 break;
335 case kFfactor:
336 if (pix.IsFFactorMethodValid())
337 {
338 calibQE = qe.GetQECascadesFFactor ( zenith );
339 calibQEVar = qe.GetQECascadesFFactorVar( zenith );
340 }
341 else
342 continue;
343 break;
344 case kCombined:
345 if (qe.IsCombinedMethodValid())
346 {
347 calibQE = qe.GetQECascadesCombined ( zenith );
348 calibQEVar = qe.GetQECascadesCombinedVar( zenith );
349 }
350 else
351 continue;
352 break;
353 case kDummy:
354 hiloconv = 1.;
355 hiloconverr = 0.;
356 calibQE = 1.;
357 calibQEVar = 0.;
358 break;
359
360 } /* switch calibration mode */
361 } /* if(fCalibrationMode!=kNone) */
362 else
363 {
364 hiloconv = 1.;
365 hiloconverr = 0.;
366 calibConv = 1./fGeomCam->GetPixRatio(pixidx);
367 calibConvVar = 0.;
368 calibFFactor = 0.;
369 calibQE = 1.;
370 calibQEVar = 0.;
371 }
372
373 MExtractedSignalPix &sig = (*fSignals)[pixidx];
374
375 Float_t signal;
376 Float_t signalErr = 0.;
377 Float_t nphot,nphotErr;
378
379 if (sig.IsLoGainUsed())
380 {
381 signal = sig.GetExtractedSignalLoGain()*hiloconv;
382 signalErr = signal*hiloconverr;
383 }
384 else
385 {
386 if (sig.GetExtractedSignalHiGain() > 9999.)
387 {
388 signal = 0.;
389 signalErr = 0.;
390 }
391 else
392 signal = sig.GetExtractedSignalHiGain();
393 }
394
395 nphot = signal*calibConv/calibQE;
396 nphotErr = calibFFactor*TMath::Sqrt(TMath::Abs(nphot));
397
398 //
399 // The following part is the commented first version of the error calculation
400 // Contact Markus Gaug for questions (or wait for the next documentation update...)
401 //
402 /*
403 nphotErr = signal > 0 ? signalErr*signalErr / (signal * signal) : 0.
404 + calibConv > 0 ? calibConvVar / (calibConv * calibConv ) : 0.
405 + calibQE > 0 ? calibQEVar / (calibQE * calibQE ) : 0.;
406 nphotErr = TMath::Sqrt(nphotErr) * nphot;
407 */
408
409 MCerPhotPix *cpix = fCerPhotEvt->AddPixel(pixidx, nphot, nphotErr);
410
411 if (sig.GetNumLoGainSaturated() > 0)
412 cpix->SetPixelSaturated();
413
414 if (sig.GetNumHiGainSaturated() > 0)
415 cpix->SetPixelHGSaturated();
416
417
418 } /* for (UInt_t pixidx=0; pixidx<npix; pixidx++) */
419
420 fCerPhotEvt->FixSize();
421 fCerPhotEvt->SetReadyToSave();
422
423 return kTRUE;
424}
Note: See TracBrowser for help on using the repository browser.