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

Last change on this file since 3939 was 3907, checked in by gaug, 21 years ago
*** empty log message ***
File size: 12.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): 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 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
197
198 switch(fCalibrationMode)
199 {
200 case kBlindPixel:
201 break;
202 case kFfactor:
203 break;
204 case kPinDiode:
205 *fLog << err << GetDescriptor()
206 << ": PIN Diode Calibration mode not yet available " << endl;
207 return kFALSE;
208 break;
209 case kCombined:
210 *fLog << err << GetDescriptor()
211 << ": Combined Calibration mode not yet available " << endl;
212 return kFALSE;
213 break;
214 case kFlatCharge:
215 *fLog << warn << GetDescriptor()
216 << ": WARNING: Flat-fielding charges - only for Keiichi!!" << endl;
217 break;
218 case kDummy:
219 *fLog << warn << GetDescriptor()
220 << ": WARNING: Dummy calibration, no calibration applied!!" << endl;
221 break;
222 case kNone:
223 *fLog << warn << GetDescriptor()
224 << ": WARNING: No calibration applied!!" << endl;
225 break;
226 default:
227 *fLog << warn << GetDescriptor()
228 << ": WARNING: Calibration mode value ("
229 <<fCalibrationMode<<") not known" << endl;
230 return kFALSE;
231 }
232
233 return kTRUE;
234}
235// --------------------------------------------------------------------------
236//
237// Apply the calibration factors to the extracted signal according to the
238// selected calibration method
239//
240Int_t MCalibrate::Process()
241{
242
243 //
244 // For the moment, we use only a dummy zenith for the calibration:
245 //
246 const Float_t zenith = 0;
247
248 /*
249 if (fCalibrations->GetNumPixels() != (UInt_t)fSignals->GetSize())
250 {
251 // FIXME: MExtractedSignal must be of variable size -
252 // like MCerPhotEvt - because we must be able
253 // to reduce size by zero supression
254 // For the moment this check could be done in ReInit...
255 *fLog << err << "MExtractedSignal and MCalibrationCam have different sizes... abort." << endl;
256 return kFALSE;
257 }
258 */
259
260 UInt_t npix = fSignals->GetSize();
261
262 Float_t hiloconv = 0.;
263 Float_t hiloconverr = 0.;
264 Float_t calibConv = 0.;
265 Float_t calibConvVar = 0.;
266 Float_t calibFFactor = 0.;
267 Float_t calibQE = 1.;
268 Float_t calibQEVar = 0.;
269 Float_t avMean = 1.;
270 Float_t avMeanRelVar = 0.;
271
272
273 if (fCalibrationMode == kFlatCharge)
274 {
275 MCalibrationChargePix &avpix = (MCalibrationChargePix&)fCalibrations->GetAverageArea(0);
276 avMean = avpix.GetMean();
277 avMeanRelVar = avpix.GetMeanRelVar();
278 }
279
280
281 for (UInt_t pixidx=0; pixidx<npix; pixidx++)
282 {
283
284 if(fCalibrationMode!=kNone)
285 {
286
287 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCalibrations)[pixidx];
288
289 hiloconv = pix.GetConversionHiLo ();
290 hiloconverr= pix.GetConversionHiLoErr();
291
292 if (fBadPixels)
293 {
294 MBadPixelsPix &bad = (*fBadPixels)[pixidx];
295 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
296 continue;
297 }
298
299 calibConv = pix.GetMeanConvFADC2Phe();
300 calibConvVar = pix.GetMeanConvFADC2PheVar();
301 calibFFactor = pix.GetMeanFFactorFADC2Phot();
302
303 if (fCalibrationMode== kFlatCharge)
304 {
305 calibConv = avMean / pix.GetMean() / fGeomCam->GetPixRatio(pixidx) ;
306 calibConvVar = (avMeanRelVar + pix.GetMeanRelVar()) * calibConv * calibConv;
307 calibFFactor = pix.GetRSigmaPerCharge();
308 }
309 else
310 {
311
312 MCalibrationQEPix &qe = (MCalibrationQEPix&) (*fQEs)[pixidx];
313
314 switch(fCalibrationMode)
315 {
316 case kBlindPixel:
317 if (qe.IsBlindPixelMethodValid())
318 {
319 calibQE = qe.GetQECascadesBlindPixel ( zenith );
320 calibQEVar = qe.GetQECascadesBlindPixelVar( zenith );
321 }
322 else
323 continue;
324 break;
325 case kPinDiode:
326 if (qe.IsPINDiodeMethodValid())
327 {
328 calibQE = qe.GetQECascadesPINDiode ( zenith );
329 calibQEVar = qe.GetQECascadesPINDiodeVar( zenith );
330 }
331 else
332 continue;
333 break;
334 case kFfactor:
335 if (pix.IsFFactorMethodValid())
336 {
337 calibQE = qe.GetQECascadesFFactor ( zenith );
338 calibQEVar = qe.GetQECascadesFFactorVar( zenith );
339 }
340 else
341 continue;
342 break;
343 case kCombined:
344 if (qe.IsCombinedMethodValid())
345 {
346 calibQE = qe.GetQECascadesCombined ( zenith );
347 calibQEVar = qe.GetQECascadesCombinedVar( zenith );
348 }
349 else
350 continue;
351 break;
352 case kDummy:
353 hiloconv = 1.;
354 hiloconverr = 0.;
355 calibQE = 1.;
356 calibQEVar = 0.;
357 break;
358
359 } /* switch calibration mode */
360 } /* else fCalibrationMode == kFlatCharge */
361 } /* if(fCalibrationMode!=kNone) */
362 else
363 {
364 hiloconv = 1.;
365 hiloconverr = 0.;
366 calibConv = 1.;
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 if (fCalibrationMode == kFlatCharge)
399 nphotErr = calibFFactor * signal;
400
401
402 //
403 // The following part is the outcommented first version of the error calculation
404 // Contact Markus Gaug for questions (or wait for the next documentation update...)
405 //
406 /*
407 nphotErr = signal > 0 ? signalErr*signalErr / (signal * signal) : 0.
408 + calibConv > 0 ? calibConvVar / (calibConv * calibConv ) : 0.
409 + calibQE > 0 ? calibQEVar / (calibQE * calibQE ) : 0.;
410 nphotErr = TMath::Sqrt(nphotErr) * nphot;
411 */
412
413 MCerPhotPix *cpix = fCerPhotEvt->AddPixel(pixidx, nphot, nphotErr);
414
415 if (sig.GetNumLoGainSaturated() > 0)
416 cpix->SetPixelSaturated();
417
418 if (sig.GetNumHiGainSaturated() > 0)
419 cpix->SetPixelHGSaturated();
420
421
422 } /* for (UInt_t pixidx=0; pixidx<npix; pixidx++) */
423
424 fCerPhotEvt->FixSize();
425 fCerPhotEvt->SetReadyToSave();
426
427 return kTRUE;
428}
Note: See TracBrowser for help on using the repository browser.