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

Last change on this file since 5478 was 5030, checked in by tbretz, 20 years ago
*** empty log message ***
File size: 12.7 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 << "WARNING - Blind pixel calibration method not valid, switching to F-factor method" << endl;
166 fCalibrationMode = kFfactor;
167 }
168
169 if(fCalibrationMode == kPinDiode && !fQEs->IsPINDiodeMethodValid())
170 {
171 *fLog << warn << "WARNING - PIN diode calibration method not valid, switching to F-factor method" << endl;
172 fCalibrationMode = kFfactor;
173 }
174
175 if(fCalibrationMode == kCombined && !fQEs->IsCombinedMethodValid())
176 {
177 *fLog << warn << "WARNING - Combined calibration method not valid, switching to F-factor method" << endl;
178 fCalibrationMode = kFfactor;
179 }
180
181
182 if (fCalibrationMode == kFlatCharge)
183 {
184 fGeomCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
185 if (!fGeomCam)
186 {
187 *fLog << err << AddSerialNumber("MGeomCam") << " not found ... aborting." << endl;
188 return kFALSE;
189 }
190 }
191
192
193 switch(fCalibrationMode)
194 {
195 case kBlindPixel:
196 break;
197 case kFfactor:
198 break;
199 case kPinDiode:
200 *fLog << err << "PIN Diode Calibration mode not yet available " << endl;
201 return kFALSE;
202 break;
203 case kCombined:
204 *fLog << err << "Combined Calibration mode not yet available " << endl;
205 return kFALSE;
206 break;
207 case kFlatCharge:
208 *fLog << warn << "WARNING - Flat-fielding charges - only for Keiichi!!" << endl;
209 break;
210 case kDummy:
211 *fLog << warn << "WARNING - Dummy calibration, no calibration applied!!" << endl;
212 break;
213 case kNone:
214 *fLog << warn << "WARNING - No calibration applied!!" << endl;
215 break;
216 default:
217 *fLog << warn << "WARNING - Calibration mode value (" <<fCalibrationMode<<") not known" << endl;
218 return kFALSE;
219 }
220
221 return kTRUE;
222}
223// --------------------------------------------------------------------------
224//
225// Apply the calibration factors to the extracted signal according to the
226// selected calibration method
227//
228Int_t MCalibrate::Process()
229{
230
231 //
232 // For the moment, we use only a dummy zenith for the calibration:
233 //
234 const Float_t zenith = 0;
235
236 /*
237 if (fCalibrations->GetNumPixels() != (UInt_t)fSignals->GetSize())
238 {
239 // FIXME: MExtractedSignal must be of variable size -
240 // like MCerPhotEvt - because we must be able
241 // to reduce size by zero supression
242 // For the moment this check could be done in ReInit...
243 *fLog << err << "MExtractedSignal and MCalibrationCam have different sizes... abort." << endl;
244 return kFALSE;
245 }
246 */
247
248 UInt_t npix = fSignals->GetSize();
249
250 Float_t hiloconv = 0.;
251 Float_t hiloconverr = 0.;
252 Float_t calibConv = 0.;
253 Float_t calibConvVar = 0.;
254 Float_t calibFFactor = 0.;
255 Float_t calibQE = 1.;
256 Float_t calibQEVar = 0.;
257 Float_t avMean = 1.;
258 Float_t avMeanRelVar = 0.;
259
260
261 if (fCalibrationMode == kFlatCharge)
262 {
263 MCalibrationChargePix &avpix = (MCalibrationChargePix&)fCalibrations->GetAverageArea(0);
264 avMean = avpix.GetMean();
265 avMeanRelVar = avpix.GetMeanRelVar();
266 }
267
268
269 for (UInt_t pixidx=0; pixidx<npix; pixidx++)
270 {
271
272 if(fCalibrationMode!=kNone)
273 {
274
275 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCalibrations)[pixidx];
276
277 hiloconv = pix.GetConversionHiLo ();
278 hiloconverr= pix.GetConversionHiLoErr();
279
280 if (fBadPixels)
281 {
282 MBadPixelsPix &bad = (*fBadPixels)[pixidx];
283 if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
284 continue;
285 }
286
287 calibConv = pix.GetMeanConvFADC2Phe();
288 calibConvVar = pix.GetMeanConvFADC2PheVar();
289 calibFFactor = pix.GetMeanFFactorFADC2Phot();
290
291 MCalibrationQEPix &qe = (MCalibrationQEPix&) (*fQEs)[pixidx];
292
293 switch(fCalibrationMode)
294 {
295 case kFlatCharge:
296 calibConv = avMean / pix.GetMean() / fGeomCam->GetPixRatio(pixidx) ;
297 calibConvVar = (avMeanRelVar + pix.GetMeanRelVar()) * calibConv * calibConv;
298 if (pix.IsFFactorMethodValid())
299 {
300 const Float_t convmin1 = qe.GetQECascadesFFactor(zenith)/pix.GetMeanConvFADC2Phe();
301 if (convmin1 > 0)
302 calibFFactor *= TMath::Sqrt(convmin1);
303 else
304 calibFFactor = -1.;
305 }
306 break;
307 case kBlindPixel:
308 if (qe.IsBlindPixelMethodValid())
309 {
310 calibQE = qe.GetQECascadesBlindPixel ( zenith );
311 calibQEVar = qe.GetQECascadesBlindPixelVar( zenith );
312 }
313 else
314 continue;
315 break;
316 case kPinDiode:
317 if (qe.IsPINDiodeMethodValid())
318 {
319 calibQE = qe.GetQECascadesPINDiode ( zenith );
320 calibQEVar = qe.GetQECascadesPINDiodeVar( zenith );
321 }
322 else
323 continue;
324 break;
325 case kFfactor:
326 if (pix.IsFFactorMethodValid())
327 {
328 calibQE = qe.GetQECascadesFFactor ( zenith );
329 calibQEVar = qe.GetQECascadesFFactorVar( zenith );
330 }
331 else
332 continue;
333 break;
334 case kCombined:
335 if (qe.IsCombinedMethodValid())
336 {
337 calibQE = qe.GetQECascadesCombined ( zenith );
338 calibQEVar = qe.GetQECascadesCombinedVar( zenith );
339 }
340 else
341 continue;
342 break;
343 case kDummy:
344 hiloconv = 1.;
345 hiloconverr = 0.;
346 calibQE = 1.;
347 calibQEVar = 0.;
348 break;
349
350 } /* switch calibration mode */
351 } /* if(fCalibrationMode!=kNone) */
352 else
353 {
354 hiloconv = 1.;
355 hiloconverr = 0.;
356 calibConv = 1./fGeomCam->GetPixRatio(pixidx);
357 calibConvVar = 0.;
358 calibFFactor = 0.;
359 calibQE = 1.;
360 calibQEVar = 0.;
361 }
362
363 MExtractedSignalPix &sig = (*fSignals)[pixidx];
364
365 Float_t signal;
366 Float_t signalErr = 0.;
367 Float_t nphot,nphotErr;
368
369 if (sig.IsLoGainUsed())
370 {
371 signal = sig.GetExtractedSignalLoGain()*hiloconv;
372 signalErr = signal*hiloconverr;
373 }
374 else
375 {
376 if (sig.GetExtractedSignalHiGain() > 9999.)
377 {
378 signal = 0.;
379 signalErr = 0.;
380 }
381 else
382 signal = sig.GetExtractedSignalHiGain();
383 }
384
385 nphot = signal*calibConv/calibQE;
386 nphotErr = calibFFactor*TMath::Sqrt(TMath::Abs(nphot));
387
388 //
389 // The following part is the commented first version of the error calculation
390 // Contact Markus Gaug for questions (or wait for the next documentation update...)
391 //
392 /*
393 nphotErr = signal > 0 ? signalErr*signalErr / (signal * signal) : 0.
394 + calibConv > 0 ? calibConvVar / (calibConv * calibConv ) : 0.
395 + calibQE > 0 ? calibQEVar / (calibQE * calibQE ) : 0.;
396 nphotErr = TMath::Sqrt(nphotErr) * nphot;
397 */
398
399 MCerPhotPix *cpix = fCerPhotEvt->AddPixel(pixidx, nphot, nphotErr);
400
401 if (sig.GetNumLoGainSaturated() > 0)
402 cpix->SetPixelSaturated();
403
404 if (sig.GetNumHiGainSaturated() > 0)
405 cpix->SetPixelHGSaturated();
406
407
408 } /* for (UInt_t pixidx=0; pixidx<npix; pixidx++) */
409
410 fCerPhotEvt->FixSize();
411 fCerPhotEvt->SetReadyToSave();
412
413 return kTRUE;
414}
Note: See TracBrowser for help on using the repository browser.