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

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