source: branches/Mars_MC/mcalib/MMcCalibrationCalc.cc@ 17758

Last change on this file since 17758 was 8024, checked in by tbretz, 18 years ago
*** empty log message ***
File size: 11.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-2004
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MMcCalibrationCalc
28//
29// Input Containers:
30// MMcConfigRunHeader
31// MRawRunHeader
32// MMcFadcHeader
33// MHillas
34// MNewImagePar
35// MMcEvt
36//
37// Output Containers:
38// (containers mus exist already, they are filled with new values).
39// MCalibrationChargeCam
40// MCalibrationQECam
41//
42/////////////////////////////////////////////////////////////////////////////
43#include "MMcCalibrationCalc.h"
44
45#include <TH1.h>
46
47#include "MLog.h"
48#include "MLogManip.h"
49
50#include "MParList.h"
51
52#include "MCalibrationChargePix.h"
53#include "MCalibrationChargeCam.h"
54
55#include "MCalibrationQEPix.h"
56#include "MCalibrationQECam.h"
57
58#include "MGeomCam.h"
59#include "MRawRunHeader.h"
60#include "MMcConfigRunHeader.h"
61
62#include "MHillas.h"
63#include "MImagePar.h"
64#include "MNewImagePar.h"
65
66#include "MMcEvt.hxx"
67#include "MMcFadcHeader.hxx"
68
69ClassImp(MMcCalibrationCalc);
70
71using namespace std;
72
73// --------------------------------------------------------------------------
74//
75// Constructor. Default value for fMinSize is 1000 ADC counts. This must be
76// set in general by the user (SetMinSize), since it depends among other things
77// on the signal extractor used.
78//
79MMcCalibrationCalc::MMcCalibrationCalc(const char *name, const char *title): fMinSize(1000)
80{
81 fName = name ? name : "MMcCalibrationCalc";
82 fTitle = title ? title : "Calculate and write conversion factors into MCalibrationChargeCam and MCalibrationQECam containers";
83}
84
85// --------------------------------------------------------------------------
86//
87// Check for the run type. Return kTRUE if it is a MC run or if there
88// is no MC run header (old camera files) kFALSE in case of a different
89// run type
90//
91Bool_t MMcCalibrationCalc::CheckRunType(MParList *pList) const
92{
93 const MRawRunHeader *run = (MRawRunHeader*)pList->FindObject(AddSerialNumber("MRawRunHeader"));
94 if (!run)
95 {
96 *fLog << warn << "Warning - cannot check file type, " << AddSerialNumber("MRawRunHeader")
97 << " not found." << endl;
98 return kTRUE;
99 }
100
101 return run->IsMonteCarloRun();
102}
103
104// --------------------------------------------------------------------------
105//
106// Look for all necessary containers and create histograms
107//
108Int_t MMcCalibrationCalc::PreProcess(MParList *pList)
109{
110 fADC2PhotEl = 0;
111 fPhot2PhotEl = 0;
112
113 fCalCam = (MCalibrationChargeCam*) pList->FindObject(AddSerialNumber("MCalibrationChargeCam"));
114 if (!fCalCam)
115 {
116 *fLog << err << AddSerialNumber("MCalibrationChargeCam") << "not found... aborting." << endl;
117 return kFALSE;
118 }
119
120 fQECam = (MCalibrationQECam*) pList->FindObject(AddSerialNumber("MCalibrationQECam"));
121 if (!fQECam)
122 {
123 *fLog << err << AddSerialNumber("MCalibrationQECam") << "not found... aborting." << endl;
124 return kFALSE;
125 }
126
127 fHillas = (MHillas*) pList->FindObject(AddSerialNumber("MHillas"));
128 if ( !fHillas)
129 {
130 *fLog << err << AddSerialNumber("MHillas") << "not found... aborting." << endl;
131 return kFALSE;
132 }
133
134 fNew = (MNewImagePar*)pList->FindObject(AddSerialNumber("MNewImagePar"));
135 if (!fNew)
136 {
137 *fLog << err << AddSerialNumber("MNewImagePar") << "not found... aborting." << endl;
138 return kFALSE;
139 }
140
141 fPar = (MImagePar*)pList->FindObject(AddSerialNumber("MImagePar"));
142 if (!fPar)
143 {
144 *fLog << err << AddSerialNumber("MImagePar") << "not found... aborting." << endl;
145 return kFALSE;
146 }
147
148 fMcEvt = (MMcEvt*) pList->FindObject(AddSerialNumber("MMcEvt"));
149 if (!fMcEvt)
150 {
151 *fLog << err << AddSerialNumber("MMcEvt") << "not found... aborting." << endl;
152 return kFALSE;
153 }
154
155 //
156 // Create histograms:
157 //
158
159 fHistADC2PhotEl = new TH1F(AddSerialNumber("ADC2PhotEl"), "log10(fPhotElfromShower/fSize)", 1500, -3., 3.);
160 fHistADC2PhotEl->SetXTitle("log_{10}(fPhotElfromShower / fSize) [photel/ADC count]");
161
162
163 fHistPhot2PhotEl = new TH1F(AddSerialNumber("Phot2PhotEl"), "Photon conversion efficiency", 1000, 0., 1.);
164 fHistPhot2PhotEl->SetXTitle("Overall photon conversion efficiency [photoelectron/photon]");
165
166
167 return kTRUE;
168}
169
170// --------------------------------------------------------------------------
171//
172// Check for the runtype.
173// Search for MGeomCam and MMcFadcHeader.
174//
175Bool_t MMcCalibrationCalc::ReInit(MParList *pList)
176{
177 //
178 // If it is no MC file display error and exit
179 //
180 if (!CheckRunType(pList))
181 {
182 *fLog << err << "MMcCalibrationCalc can only be used with MC files... aborting." << endl;
183 return kFALSE;
184 }
185
186 //
187 // Now check the existence of all necessary containers.
188 //
189 fGeom = (MGeomCam*) pList->FindObject(AddSerialNumber("MGeomCam"));
190 if (!fGeom)
191 {
192 *fLog << err << AddSerialNumber("MGeomCam") << " not found... aborting." << endl;
193 return kFALSE;
194 }
195
196 fHeaderFadc = (MMcFadcHeader*)pList->FindObject(AddSerialNumber("MMcFadcHeader"));
197 if (!fHeaderFadc)
198 {
199 *fLog << err << AddSerialNumber("MMcFadcHeader") << " not found... aborting." << endl;
200 return kFALSE;
201 }
202
203 for (UInt_t ipix = 0; ipix < fGeom->GetNumPixels(); ipix++)
204 {
205 if (fHeaderFadc->GetPedestalRmsHigh(ipix) > 0 ||
206 fHeaderFadc->GetPedestalRmsLow(ipix) > 0 )
207 {
208 *fLog << err << "Trying to calibrate the data using a Camera file produced with added noise." << endl;
209 *fLog << "Please use a noiseless file for calibration... aborting." << endl << endl;
210 return kFALSE;
211 }
212 }
213
214 // Now check the light collection for inner and outer pixels to
215 // calculate the ratio between the two. FIXME! Light collection
216 // depends on the incidence angle of the light w.r.t. the camera
217 // plane. For the moment we take the ratio for light impinging
218 // perpendicular to the camera plane.
219 //
220 MMcConfigRunHeader* mcconfig = (MMcConfigRunHeader*) pList->FindObject(AddSerialNumber("MMcConfigRunHeader"));
221 if (!mcconfig)
222 {
223 *fLog << err << AddSerialNumber("MMcConfigRunHeader") << " not found... aborting." << endl;
224 return kFALSE;
225 }
226 const TArrayF &innerlightcoll = mcconfig->GetLightCollectionFactor();
227 const TArrayF &outerlightcoll = mcconfig->GetLightCollectionFactorOuter();
228
229 // In principle outer pixels seem to have a different average light
230 // collection efficiency than outer ones. We set here the factor between
231 // the two.
232
233 fOuterPixelsLightCollection = outerlightcoll[90] / innerlightcoll[90];
234 // (at angle = 90 deg)
235
236
237 return kTRUE;
238}
239
240
241// --------------------------------------------------------------------------
242//
243// Obtain average ratio of photons in camera to image Size.
244//
245Int_t MMcCalibrationCalc::Process()
246{
247 //
248 // Exclude events with some high-gain saturated pixel
249 //
250 if (fPar->GetNumSatPixelsHG()>0)
251 return kTRUE;
252
253 const Float_t size = fHillas->GetSize();
254 const Float_t innersize = fNew->GetInnerSize();
255
256 // Size will at this point be in ADC counts (still uncalibrated)
257 //
258 // Exclude events with low Size (larger fluctuations)
259 //
260
261 if (size < fMinSize)
262 return kTRUE;
263
264 //
265 // PATCH: Convert number of photoelectrons in camera to the approximate number
266 // of photoelectrons that would have been registered if all pixels had the same
267 // light collection efficiency as inner ones (called here "corrected_photel").
268 //
269
270 const Float_t corrected_photel = (Float_t) fMcEvt->GetPhotElfromShower() /
271 (fOuterPixelsLightCollection + innersize / size * (1. - fOuterPixelsLightCollection));
272
273 fHistADC2PhotEl->Fill(TMath::Log10(corrected_photel/size));
274 fHistPhot2PhotEl->Fill( corrected_photel / (Float_t) fMcEvt->GetPassPhotCone() );
275
276 return kTRUE;
277}
278
279// --------------------------------------------------------------------------
280//
281// Fill the MCalibrationCam object
282//
283Int_t MMcCalibrationCalc::PostProcess()
284{
285 const Stat_t n = fHistADC2PhotEl->GetEntries();
286 if (n<1)
287 {
288 *fLog << err << "No events read... aborting." << endl;
289 return kFALSE;
290 }
291
292 fPhot2PhotEl = fHistPhot2PhotEl->GetMean();
293 // Average quantum efficiency. For now we will set this value for all pixels, although
294 // even in MC there may be differences from pixel to pixel, if the .dat file containing
295 // QE vs lambda, input of the camera simulation, has different QE curves for different
296 // pixels. FIXME?
297
298 MCalibrationQEPix &avqepix = (MCalibrationQEPix&)(fQECam->GetAverageArea(0));
299 avqepix.SetAverageQE(fPhot2PhotEl); // Needed by MCalibrateData!
300
301 //
302 // Find peak of log10(photel/Size) histogram:
303 //
304 const Int_t reach = 2;
305 Stat_t summax = 0;
306 Int_t mode = 0;
307 for (Int_t ibin = 1+reach; ibin <= fHistADC2PhotEl->GetNbinsX()-reach; ibin++)
308 {
309 const Stat_t sum = fHistADC2PhotEl->Integral(ibin-reach, ibin+reach);
310
311 if (sum <= summax)
312 continue;
313
314 summax = sum;
315 mode = ibin;
316 }
317
318 fADC2PhotEl = TMath::Power(10, fHistADC2PhotEl->GetBinCenter(mode));
319
320 const Int_t num = fCalCam->GetSize();
321
322 for (int i=0; i<num; i++)
323 {
324 MCalibrationChargePix &calpix = (MCalibrationChargePix&)(*fCalCam)[i];
325 MCalibrationQEPix &qepix = (MCalibrationQEPix&) (*fQECam) [i];
326
327 Float_t qe = fPhot2PhotEl;
328 if (fGeom->GetPixRatio(i) < 1.)
329 qe *= fOuterPixelsLightCollection;
330 qepix.SetAverageQE(qe);
331
332 qepix.SetAvNormFFactor(1.);
333 // This factor should convert the default average QE for different colors to
334 // average QE for a spectrum like that of Cherenkov light (see the documentration
335 // of MCalibrationQEPix).
336 // Here we obtain average QE using already a Cherenkov spectrum so AvNormFFactor
337 // must be 1.
338
339
340 Float_t factor = fADC2PhotEl;
341
342 //
343 // We take into account the (possibly) different gain of outer pixels:
344 // FIXME: we are now assuming that all inner pixels have the same gain, and all
345 // outer pixels have the same gain (different from inner ones though). This can
346 // only be like this in camera 0.7, but might change in future versions of camera.
347 //
348
349 if (fGeom->GetPixRatio(i) < 1.)
350 factor *= fHeaderFadc->GetAmplitud()/fHeaderFadc->GetAmplitudOuter();
351
352 calpix.SetMeanConvFADC2Phe(factor);
353 calpix.SetMeanConvFADC2PheVar(0.);
354
355 calpix.SetMeanFFactorFADC2Phot(0.);
356 }
357
358
359 return kTRUE;
360}
Note: See TracBrowser for help on using the repository browser.