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

Last change on this file since 4882 was 4605, 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#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 fQEs = (MCalibrationQECam*)pList->FindObject(AddSerialNumber("MCalibrationQECam"));
138 if (!fQEs)
139 {
140 *fLog << err << AddSerialNumber("MCalibrationQECam") << " not found ... aborting." << endl;
141 return kFALSE;
142 }
143 }
144
145 fCerPhotEvt = (MCerPhotEvt*)pList->FindCreateObj(AddSerialNumber("MCerPhotEvt"));
146 if (!fCerPhotEvt)
147 return kFALSE;
148
149 fGeomCam = (MGeomCam*)pList->FindCreateObj(AddSerialNumber("MGeomCam"));
150 if (!fGeomCam)
151 return kFALSE;
152
153 return kTRUE;
154}
155
156// --------------------------------------------------------------------------
157//
158// Check for validity of the selected calibration method, switch to a
159// different one in case of need
160//
161Bool_t MCalibrate::ReInit(MParList *pList)
162{
163
164 if(fCalibrationMode == kBlindPixel && !fQEs->IsBlindPixelMethodValid())
165 {
166 *fLog << warn << GetDescriptor()
167 << "Warning: Blind pixel calibration method not valid, switching to F-factor method" << endl;
168 fCalibrationMode = kFfactor;
169 }
170
171 if(fCalibrationMode == kPinDiode && !fQEs->IsPINDiodeMethodValid())
172 {
173 *fLog << warn << GetDescriptor()
174 << "Warning: PIN diode calibration method not valid, switching to F-factor method" << endl;
175 fCalibrationMode = kFfactor;
176 }
177
178 if(fCalibrationMode == kCombined && !fQEs->IsCombinedMethodValid())
179 {
180 *fLog << warn << GetDescriptor()
181 << "Warning: Combined calibration method not valid, switching to F-factor method" << endl;
182 fCalibrationMode = kFfactor;
183 }
184
185
186 if (fCalibrationMode == kFlatCharge)
187 {
188 fGeomCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
189 if (!fGeomCam)
190 {
191 *fLog << err << AddSerialNumber("MGeomCam") << " not found ... aborting." << endl;
192 return kFALSE;
193 }
194 }
195
196
197 switch(fCalibrationMode)
198 {
199 case kBlindPixel:
200 break;
201 case kFfactor:
202 break;
203 case kPinDiode:
204 *fLog << err << GetDescriptor()
205 << ": PIN Diode Calibration mode not yet available " << endl;
206 return kFALSE;
207 break;
208 case kCombined:
209 *fLog << err << GetDescriptor()
210 << ": Combined Calibration mode not yet available " << endl;
211 return kFALSE;
212 break;
213 case kFlatCharge:
214 *fLog << warn << GetDescriptor()
215 << ": WARNING: Flat-fielding charges - only for Keiichi!!" << endl;
216 break;
217 case kDummy:
218 *fLog << warn << GetDescriptor()
219 << ": WARNING: Dummy calibration, no calibration applied!!" << endl;
220 break;
221 case kNone:
222 *fLog << warn << GetDescriptor()
223 << ": WARNING: No calibration applied!!" << endl;
224 break;
225 default:
226 *fLog << warn << GetDescriptor()
227 << ": WARNING: Calibration mode value ("
228 <<fCalibrationMode<<") not known" << endl;
229 return kFALSE;
230 }
231
232 return kTRUE;
233}
234// --------------------------------------------------------------------------
235//
236// Apply the calibration factors to the extracted signal according to the
237// selected calibration method
238//
239Int_t MCalibrate::Process()
240{
241
242 //
243 // For the moment, we use only a dummy zenith for the calibration:
244 //
245 const Float_t zenith = 0;
246
247 /*
248 if (fCalibrations->GetNumPixels() != (UInt_t)fSignals->GetSize())
249 {
250 // FIXME: MExtractedSignal must be of variable size -
251 // like MCerPhotEvt - because we must be able
252 // to reduce size by zero supression
253 // For the moment this check could be done in ReInit...
254 *fLog << err << "MExtractedSignal and MCalibrationCam have different sizes... abort." << endl;
255 return kFALSE;
256 }
257 */
258
259 UInt_t npix = fSignals->GetSize();
260
261 Float_t hiloconv = 0.;
262 Float_t hiloconverr = 0.;
263 Float_t calibConv = 0.;
264 Float_t calibConvVar = 0.;
265 Float_t calibFFactor = 0.;
266 Float_t calibQE = 1.;
267 Float_t calibQEVar = 0.;
268 Float_t avMean = 1.;
269 Float_t avMeanRelVar = 0.;
270
271
272 if (fCalibrationMode == kFlatCharge)
273 {
274 MCalibrationChargePix &avpix = (MCalibrationChargePix&)fCalibrations->GetAverageArea(0);
275 avMean = avpix.GetMean();
276 avMeanRelVar = avpix.GetMeanRelVar();
277 }
278
279
280 for (UInt_t pixidx=0; pixidx<npix; pixidx++)
281 {
282
283 if(fCalibrationMode!=kNone)
284 {
285
286 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCalibrations)[pixidx];
287
288 hiloconv = pix.GetConversionHiLo ();
289 hiloconverr= pix.GetConversionHiLoErr();
290
291 if (fBadPixels)
292 {
293 MBadPixelsPix &bad = (*fBadPixels)[pixidx];
294 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
295 continue;
296 }
297
298 calibConv = pix.GetMeanConvFADC2Phe();
299 calibConvVar = pix.GetMeanConvFADC2PheVar();
300 calibFFactor = pix.GetMeanFFactorFADC2Phot();
301
302 MCalibrationQEPix &qe = (MCalibrationQEPix&) (*fQEs)[pixidx];
303
304 switch(fCalibrationMode)
305 {
306 case kFlatCharge:
307 calibConv = avMean / pix.GetMean() / fGeomCam->GetPixRatio(pixidx) ;
308 calibConvVar = (avMeanRelVar + pix.GetMeanRelVar()) * calibConv * calibConv;
309 if (pix.IsFFactorMethodValid())
310 {
311 const Float_t convmin1 = qe.GetQECascadesFFactor(zenith)/pix.GetMeanConvFADC2Phe();
312 if (convmin1 > 0)
313 calibFFactor *= TMath::Sqrt(convmin1);
314 else
315 calibFFactor = -1.;
316 }
317 break;
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 } /* if(fCalibrationMode!=kNone) */
363 else
364 {
365 hiloconv = 1.;
366 hiloconverr = 0.;
367 calibConv = 1./fGeomCam->GetPixRatio(pixidx);
368 calibConvVar = 0.;
369 calibFFactor = 0.;
370 calibQE = 1.;
371 calibQEVar = 0.;
372 }
373
374 MExtractedSignalPix &sig = (*fSignals)[pixidx];
375
376 Float_t signal;
377 Float_t signalErr = 0.;
378 Float_t nphot,nphotErr;
379
380 if (sig.IsLoGainUsed())
381 {
382 signal = sig.GetExtractedSignalLoGain()*hiloconv;
383 signalErr = signal*hiloconverr;
384 }
385 else
386 {
387 if (sig.GetExtractedSignalHiGain() > 9999.)
388 {
389 signal = 0.;
390 signalErr = 0.;
391 }
392 else
393 signal = sig.GetExtractedSignalHiGain();
394 }
395
396 nphot = signal*calibConv/calibQE;
397 nphotErr = calibFFactor*TMath::Sqrt(TMath::Abs(nphot));
398
399 //
400 // The following part is the commented first version of the error calculation
401 // Contact Markus Gaug for questions (or wait for the next documentation update...)
402 //
403 /*
404 nphotErr = signal > 0 ? signalErr*signalErr / (signal * signal) : 0.
405 + calibConv > 0 ? calibConvVar / (calibConv * calibConv ) : 0.
406 + calibQE > 0 ? calibQEVar / (calibQE * calibQE ) : 0.;
407 nphotErr = TMath::Sqrt(nphotErr) * nphot;
408 */
409
410 MCerPhotPix *cpix = fCerPhotEvt->AddPixel(pixidx, nphot, nphotErr);
411
412 if (sig.GetNumLoGainSaturated() > 0)
413 cpix->SetPixelSaturated();
414
415 if (sig.GetNumHiGainSaturated() > 0)
416 cpix->SetPixelHGSaturated();
417
418
419 } /* for (UInt_t pixidx=0; pixidx<npix; pixidx++) */
420
421 fCerPhotEvt->FixSize();
422 fCerPhotEvt->SetReadyToSave();
423
424 return kTRUE;
425}
Note: See TracBrowser for help on using the repository browser.