source: branches/Mars_MC/manalysis/MMcCalibrationUpdate.cc@ 17054

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