source: trunk/MagicSoft/Mars/manalysis/MMcCalibrationUpdate.cc@ 7150

Last change on this file since 7150 was 7005, checked in by tbretz, 20 years ago
*** empty log message ***
File size: 15.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): Abelardo Moralejo, 12/2003 <mailto:moralejo@pd.infn.it>
19!
20! Copyright: MAGIC Software Development, 2000-2003
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MMcCalibrationUpdate
28//
29// This task looks for the ìnformation about FADC pedestals in
30// MMcFadcHeader and translates it to the pedestal mean and rms (in adc counts).
31// If not already existing in the parameter list, an MCalibrationCam object
32// is created, with the conversion factor between ADC counts and photons or
33// photoelectrons (depending on fSignalType) is set to 1 to allow the analysis
34// to proceed.
35//
36// Then it creates and fills also the MPedPhotCam object containing the pedestal
37// mean and rms in units of photons or photoelectrons.
38//
39// Input Containers:
40// MMcFadcHeader
41// MRawRunHeader
42// [MCalibrationChargeCam] (if it existed previously)
43// [MCalibrationQECam] (if it existed previously)
44//
45// Output Containers:
46// MPedPhotCam
47// [MCalibrationChargeCam] (if it did not exist previously)
48// [MCalibrationQECam] (if it did not exist previously)
49//
50/////////////////////////////////////////////////////////////////////////////
51#include "MMcCalibrationUpdate.h"
52
53#include "MParList.h"
54
55#include "MLog.h"
56#include "MLogManip.h"
57
58#include "MCalibrationChargePix.h"
59#include "MCalibrationChargeCam.h"
60
61#include "MCalibrationQEPix.h"
62#include "MCalibrationQECam.h"
63
64#include "MExtractedSignalCam.h"
65#include "MExtractedSignalPix.h"
66#include "MGeomCam.h"
67#include "MPedPhotCam.h"
68#include "MPedPhotPix.h"
69
70#include "MRawRunHeader.h"
71#include "MMcRunHeader.hxx"
72#include "MMcFadcHeader.hxx"
73#include "MMcConfigRunHeader.h"
74#include "MCalibrateData.h"
75
76ClassImp(MMcCalibrationUpdate);
77
78using namespace std;
79
80MMcCalibrationUpdate::MMcCalibrationUpdate(const char *name, const char *title)
81 : fFillCalibrationCam(kTRUE), fOuterPixelsGainScaling(kTRUE), fAmplitude(-1.),
82 fAmplitudeOuter(-1.), fConversionHiLo(-1.), fUserLow2HiGainFactor(-1.),
83 fSignalType(MCalibrateData::kPhe)
84{
85 fName = name ? name : "MMcCalibrationUpdate";
86 fTitle = title ? title : "Write MC pedestals and conversion factors into MCalibration Container";
87}
88
89// --------------------------------------------------------------------------
90//
91// Check for the run type. Return kTRUE if it is a MC run or if there
92// is no MC run header (old camera files) kFALSE in case of a different
93// run type
94//
95Bool_t MMcCalibrationUpdate::CheckRunType(MParList *pList) const
96{
97 const MRawRunHeader *run = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
98 if (!run)
99 {
100 *fLog << warn << dbginf << "Warning - cannot check file type, MRawRunHeader not found." << endl;
101 return kTRUE;
102 }
103
104 return run->IsMonteCarloRun();
105}
106
107// --------------------------------------------------------------------------
108//
109// Make sure, that there is an MCalibrationCam Object in the Parameter List.
110//
111Int_t MMcCalibrationUpdate::PreProcess(MParList *pList)
112{
113 fCalCam = (MCalibrationChargeCam*) pList->FindObject(AddSerialNumber("MCalibrationChargeCam"));
114 fQECam = (MCalibrationQECam*) pList->FindObject(AddSerialNumber("MCalibrationQECam"));
115
116 if (!fCalCam || !fQECam)
117 {
118 fCalCam = (MCalibrationChargeCam*) pList->FindCreateObj(AddSerialNumber("MCalibrationChargeCam"));
119 fQECam = (MCalibrationQECam*) pList->FindCreateObj(AddSerialNumber("MCalibrationQECam"));
120 if (!fCalCam || !fQECam)
121 return kFALSE;
122 }
123 else
124 {
125 fFillCalibrationCam = kFALSE;
126 *fLog << inf << AddSerialNumber("MCalibrationChargeCam") << " and " <<
127 AddSerialNumber("MCalibrationQECam") << " already exist... " << endl;
128 }
129
130 fPedPhotCam = (MPedPhotCam*) pList->FindCreateObj(AddSerialNumber("MPedPhotCam"));
131 if (!fPedPhotCam)
132 return kFALSE;
133
134 fSignalCam = (MExtractedSignalCam*) pList->FindObject(AddSerialNumber("MExtractedSignalCam"));
135 if (!fSignalCam)
136 {
137 *fLog << err << AddSerialNumber("MExtractedSignalCam") << " not found... aborting." << endl;
138 return kFALSE;
139 }
140
141 return kTRUE;
142}
143
144// --------------------------------------------------------------------------
145//
146// Check for the runtype.
147// Search for MGeomCam and MMcFadcHeader.
148// Fill the MCalibrationCam object.
149//
150Bool_t MMcCalibrationUpdate::ReInit(MParList *pList)
151{
152 //
153 // If it is no MC file skip this function...
154 //
155 if (!CheckRunType(pList))
156 {
157 *fLog << inf << "This is no MC file... skipping." << endl;
158 return kTRUE;
159 }
160
161 //
162 // Now check the existence of all necessary containers.
163 //
164 fGeom = (MGeomCam*) pList->FindObject(AddSerialNumber("MGeomCam"));
165 if (!fGeom)
166 {
167 *fLog << err << AddSerialNumber("MGeomCam") << " not found... aborting." << endl;
168 return kFALSE;
169 }
170
171 fHeaderFadc = (MMcFadcHeader*)pList->FindObject(AddSerialNumber("MMcFadcHeader"));
172 if (!fHeaderFadc)
173 {
174 *fLog << err << AddSerialNumber("MMcFadcHeader") << " not found... aborting." << endl;
175 return kFALSE;
176 }
177
178 MMcRunHeader* mcrunh = (MMcRunHeader*) pList->FindObject("MMcRunHeader");
179
180 //
181 // Initialize Fadc simulation parameters:
182 //
183 if (fAmplitude < 0)
184 {
185 fAmplitude = fHeaderFadc->GetAmplitud();
186 if (mcrunh->GetCamVersion() > 60)
187 {
188 fAmplitudeOuter = fHeaderFadc->GetAmplitudOuter();
189
190 fHeaderLow2HiGainFactor = fHeaderFadc->GetLow2HighGain();
191
192 // The high to low gain ratio is stored in MMcFadcHeader.Low2HighGain.
193 // However, this is just the ratio of total pulse integrals. Since the
194 // shape of the low gain pulse is different from that of the high gain,
195 // the factor to be applied to signals extracted from low gain depends
196 // on the type of signal extractor (for instance if we extract the pulse
197 // height, the factor is larger than Low2HighGain, because the low gain
198 // pulse shape is wider and hence lower than the high gain pulse. So the
199 // user can set manually the value of the factor to be applied. If such
200 // value has not been set by the user, then we takes as default Low2HighGain.
201
202 if (fUserLow2HiGainFactor < 0.)
203 fConversionHiLo = fHeaderLow2HiGainFactor;
204 else
205 fConversionHiLo = fUserLow2HiGainFactor;
206
207 }
208 else // old MC files, camera < v0.7
209 {
210 fAmplitudeOuter = fAmplitude;
211 fConversionHiLo = 10; // dummy value
212 }
213
214 }
215 else // Check that the following files have all the same FADC parameters as the first
216 {
217 if ( fabs(fHeaderFadc->GetAmplitud()-fAmplitude) > 1.e-6 )
218 {
219 *fLog << err << "Parameters of MMcFadcHeader are not the same for all files... aborting." << endl;
220 return kFALSE;
221 }
222
223 if (mcrunh->GetCamVersion() > 60) // files from camera 0.7 or newer
224 {
225 if( fabs(fHeaderFadc->GetAmplitudOuter()-fAmplitudeOuter) > 1.e-6 ||
226 fabs(fHeaderLow2HiGainFactor-fHeaderFadc->GetLow2HighGain()) > 1.e-6 )
227 {
228 *fLog << err << "Parameters of MMcFadcHeader are not the same for all files... aborting." << endl;
229 return kFALSE;
230 }
231 }
232 }
233
234 //
235 // If MCalibrationChargeCam and MCalibrationQECam already existed
236 // in the parameter list before MMcCalibrationUpdate::PreProcess was
237 // executed (from a previous calibration loop) we must not fill them,
238 // hence nothing else has to be done in ReInit:
239 //
240 if (!fFillCalibrationCam)
241 return kTRUE;
242
243 // Now check the light collection for inner and outer pixels to
244 // calculate the ratio between the two. FIXME! Light collection
245 // depends on the incidence angle of the light w.r.t. the camera
246 // plane. For the moment we take the ratio for light impinging
247 // perpendicular to the camera plane.
248 //
249 // FIXME! We look for AddSerialNumber("MMcConfigRunHeader") but
250 // for the moment the stereo version of camera does not write one such
251 // header per telescope (it should!)
252 //
253 MMcConfigRunHeader* mcconfig = (MMcConfigRunHeader*) pList->FindObject(AddSerialNumber("MMcConfigRunHeader"));
254 if (!mcconfig)
255 {
256 *fLog << err << AddSerialNumber("MMcConfigRunHeader") <<
257 " not found... aborting." << endl;
258 return kFALSE;
259 }
260 TArrayF innerlightcoll = mcconfig->GetLightCollectionFactor();
261 TArrayF outerlightcoll = mcconfig->GetLightCollectionFactorOuter();
262
263 // In principle outer pixels seem to have a different average light
264 // collection efficiency than outer ones. We set here the factor between
265 // the two.
266
267 fOuterPixelsLightCollection = outerlightcoll[90] / innerlightcoll[90];
268 // (at angle = 90 deg)
269
270 // Set now the default conversion from ADC counts to photoelectrons
271 // (in case no previous calibration existed in the parameter list).
272 //
273 // As default we want to have SIZE in ADC counts, or rather, in "inner pixel
274 // equivalent ADC counts".
275 //
276 // To achieve this:
277 // - In the case fSignalType==kPhot: we set the ADC to photoelectron conversion
278 // equal to the QE, which will later make the ADC to photon conversion factor
279 // (= ADC2PhotEl/QE) to be = 1,
280 //
281 // - In the case fSignalType==kPhe: we set the ADC to photoelectron conversion
282 // equal to 1, since this will be applied directly to the signals...
283
284 if (fSignalType == MCalibrateData::kPhot)
285 fADC2PhElInner = MCalibrationQEPix::gkDefaultAverageQE;
286 else
287 fADC2PhElInner = 1.;
288
289 //
290 // Set the default ADC to "photoelectrons" conversion factor for outer
291 // pixels. One can choose not to apply the known (in MC) gain factor
292 // between inner and outer pixels, (in case fOuterPixelsGainScaling = kFALSE),
293 // which may be useful for display purposes.
294 // If on the contrary we apply the factor, we must take into account the
295 // different gains photoelectrons->ADC counts, given in MC by fAmplitude
296 // and fAmplitudeOuter. This "default" calibration is such that a shower
297 // completely contained in the inner part would have Size in ADC counts,
298 // whereas one partially in the outer part would have Size in "equivalent
299 // inner ADC counts" : the "same" shower (light density distribution) would
300 // have the same Size no matter where in the camera it lies. For this we have
301 // also to set later (see below) the right QE for outer pixels, which may
302 // be different from that of inner pixels.
303 //
304
305 if (fOuterPixelsGainScaling)
306 fADC2PhElOuter = fADC2PhElInner
307 * (fAmplitude / fAmplitudeOuter);
308 else
309 fADC2PhElOuter = fADC2PhElInner;
310
311
312 Int_t num = fCalCam->GetSize();
313
314 fCalCam->SetFFactorMethodValid ( kTRUE );
315 fQECam->SetFFactorMethodValid ( kTRUE );
316 fQECam->SetBlindPixelMethodValid ( kTRUE );
317 fQECam->SetCombinedMethodValid ( kTRUE );
318 fQECam->SetPINDiodeMethodValid ( kTRUE );
319
320 for (Int_t i=0; i<num; i++)
321 {
322 MCalibrationChargePix &calpix = (MCalibrationChargePix&)(*fCalCam)[i];
323
324 calpix.SetFFactorMethodValid();
325
326 calpix.SetConversionHiLo(fConversionHiLo);
327 calpix.SetConversionHiLoErr(0.); // FIXME ?
328
329 //
330 // Write conversion factor ADC to photo-electrons (different for inner
331 // and outer pixels).
332 //
333 Float_t adc2photel = (fGeom->GetPixRatio(i) < fGeom->GetPixRatio(0))?
334 fADC2PhElOuter : fADC2PhElInner;
335
336
337 calpix.SetMeanConvFADC2Phe(adc2photel);
338 calpix.SetMeanConvFADC2PheVar(0.);
339 calpix.SetMeanFFactorFADC2Phot(0.); // Not used for now.
340
341 }
342
343 //
344 // Now set the average QE for each type of pixels. Correct outer pixels
345 // for different light collection efficiency.
346 //
347 num = fQECam->GetSize();
348 for (Int_t i=0; i<num; i++)
349 {
350 MCalibrationQEPix &qepix = (MCalibrationQEPix&)(*fQECam)[i];
351
352 Float_t avqe = MCalibrationQEPix::gkDefaultAverageQE;
353
354 if (fOuterPixelsGainScaling)
355 if (fGeom->GetPixRatio(i) < fGeom->GetPixRatio(0))
356 avqe = MCalibrationQEPix::gkDefaultAverageQE*fOuterPixelsLightCollection;
357
358 qepix.SetAvNormFFactor(1.);
359 // This factor should convert the default average QE to average QE
360 // for a spectrum like that of Cherenkov light. See the documentation
361 // of MCalibrationQEPix. Here it is 1 because we calibrate using
362 // Cherenkov light.
363
364 qepix.SetAverageQE(avqe);
365 }
366
367 return kTRUE;
368}
369
370
371// --------------------------------------------------------------------------
372//
373// Fill the MCerPhotPed object
374//
375// This has to be done on an event by event basis because the (calibrated) pedestal
376// fluctuations depend on whether, for each pixel, we are using the high gain or the
377// low gain branch.
378//
379Int_t MMcCalibrationUpdate::Process()
380{
381 const Int_t num = fCalCam->GetSize();
382
383 for (Int_t i=0; i<num; i++)
384 {
385 MExtractedSignalPix &sigpix = (*fSignalCam)[i];
386
387 //
388 // ped mean and rms per pixel, in ADC counts, according to signal
389 // calculation (hi or low gain and number of integrated slices):
390 //
391 const Float_t pedestmean = sigpix.IsLoGainUsed()?
392 fSignalCam->GetNumUsedLoGainFADCSlices()*fHeaderFadc->GetPedestal(i) :
393 fSignalCam->GetNumUsedHiGainFADCSlices()*fHeaderFadc->GetPedestal(i);
394
395 //
396 // In some cases, depending on the camera simulation parameters, one can
397 // have very little or no noise in the FADC. In the case the rms of
398 // pedestal is zero, the pixel will be cleaned out later in the image
399 // cleaning. To avoid this problem,we set a default value of 0.01 ADC
400 // counts for the RMS per slice:
401 //
402 const Double_t used = (Double_t)(sigpix.IsLoGainUsed() ?
403 fSignalCam->GetNumUsedLoGainFADCSlices() :
404 fSignalCam->GetNumUsedHiGainFADCSlices());
405
406 const Float_t rms0 = sigpix.IsLoGainUsed() ?
407 fHeaderFadc->GetPedestalRmsLow(i) :
408 fHeaderFadc->GetPedestalRmsHigh(i);
409
410 const Float_t pedestrms = TMath::Sqrt(used) * (rms0>0 ? rms0 : 0.01);
411
412 //
413 // Write mean pedestal and pedestal rms per pixel
414 // in number of photons:
415 //
416 MPedPhotPix &pedpix = (*fPedPhotCam)[i];
417
418 MCalibrationChargePix &calpix = (MCalibrationChargePix&)(*fCalCam)[i];
419 MCalibrationQEPix &qepix = (MCalibrationQEPix&)(*fQECam)[i];
420
421 Float_t qe = qepix.GetAverageQE();
422
423
424 Float_t conv = (fSignalType == MCalibrateData::kPhot?
425 calpix.GetMeanConvFADC2Phe() / qe :
426 calpix.GetMeanConvFADC2Phe());
427
428 Float_t hi2lo = calpix.GetConversionHiLo();
429
430 if (sigpix.IsLoGainUsed())
431 pedpix.Set(conv*hi2lo*pedestmean, conv*hi2lo*pedestrms);
432 else
433 pedpix.Set(conv*pedestmean, conv*pedestrms);
434
435 }
436
437 return kTRUE;
438}
439
Note: See TracBrowser for help on using the repository browser.