source: tags/Mars-V0.9/manalysis/MMcCalibrationUpdate.cc

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