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

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