source: trunk/MagicSoft/Mars/mcalib/MCalibrationIntensityChargeCam.cc@ 5791

Last change on this file since 5791 was 5690, checked in by tbretz, 21 years ago
*** empty log message ***
File size: 17.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): Markus Gaug 11/2003 <mailto:markus@ifae.es>
19!
20! Copyright: MAGIC Software Development, 2000-2004
21!
22!
23\* ======================================================================== */
24/////////////////////////////////////////////////////////////////////////////
25//
26// MCalibrationIntensityChargeCam
27//
28// Storage container for intensity charge calibration results.
29//
30// Individual MCalibrationChargeCam's can be retrieved with:
31// - GetCam() yielding the current cam.
32// - GetCam("name") yielding the current camera with name "name".
33// - GetCam(i) yielding the i-th camera.
34//
35// See also: MCalibrationIntensityCam, MCalibrationChargeCam,
36// MCalibrationChargePix, MCalibrationChargeCalc, MCalibrationQECam
37// MCalibrationBlindCam, MCalibrationChargePINDiode
38// MHCalibrationChargePix, MHCalibrationChargeCam
39//
40/////////////////////////////////////////////////////////////////////////////
41#include "MCalibrationIntensityChargeCam.h"
42
43#include <TF1.h>
44#include <TH2.h>
45#include <TGraphErrors.h>
46#include <TOrdCollection.h>
47
48#include "MLog.h"
49
50#include "MGeomCam.h"
51#include "MGeomPix.h"
52
53#include "MCalibrationChargeCam.h"
54#include "MCalibrationChargePix.h"
55
56ClassImp(MCalibrationIntensityChargeCam);
57
58using namespace std;
59
60// --------------------------------------------------------------------------
61//
62// Default constructor.
63//
64MCalibrationIntensityChargeCam::MCalibrationIntensityChargeCam(const char *name, const char *title)
65{
66
67 fName = name ? name : "MCalibrationIntensityChargeCam";
68 fTitle = title ? title : "Results of the Intensity Calibration";
69
70 InitSize(1);
71}
72
73// -------------------------------------------------------------------
74//
75// Add MCalibrationChargeCam's in the ranges from - to.
76//
77void MCalibrationIntensityChargeCam::Add(const UInt_t from, const UInt_t to)
78{
79 for (UInt_t i=from; i<to; i++)
80 fCams->AddAt(new MCalibrationChargeCam,i);
81}
82
83// -------------------------------------------------------------------
84//
85// Returns a TGraphErrors with the number of photo-electrons vs.
86// the extracted signal of pixel "pixid".
87//
88TGraphErrors *MCalibrationIntensityChargeCam::GetPheVsCharge( const UInt_t pixid, const MCalibrationCam::PulserColor_t col)
89{
90
91 Int_t size = CountNumEntries(col);
92
93 if (size == 0)
94 return NULL;
95
96 TArrayF phe(size);
97 TArrayF pheerr(size);
98 TArrayF sig(size);
99 TArrayF sigerr(size);
100
101 Int_t cnt = 0;
102
103 for (Int_t i=0;i<GetSize();i++)
104 {
105 //
106 // Get the calibration cam from the intensity cam
107 //
108 MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);
109
110 if (col != MCalibrationCam::kNONE)
111 if (cam->GetPulserColor() != col)
112 continue;
113 //
114 // Get the calibration pix from the calibration cam
115 //
116 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*cam)[pixid];
117 //
118 // Don't use bad pixels
119 //
120 if (!pix.IsFFactorMethodValid())
121 continue;
122 //
123 phe[cnt] = pix.GetPheFFactorMethod();
124 pheerr[cnt] = pix.GetPheFFactorMethodErr();
125 //
126 // For the calculation of Q, we have to use the
127 // converted value!
128 //
129 sig [cnt] = pix.GetConvertedMean();
130 sigerr[cnt] = pix.GetConvertedMeanErr();
131 cnt++;
132 }
133
134 TGraphErrors *gr = new TGraphErrors(size,
135 sig.GetArray(),phe.GetArray(),
136 sigerr.GetArray(),pheerr.GetArray());
137 gr->SetTitle(Form("%s%3i","Pixel ",pixid));
138 gr->GetXaxis()->SetTitle("Q [FADC counts]");
139 gr->GetYaxis()->SetTitle("photo-electrons [1]");
140 return gr;
141}
142
143// -------------------------------------------------------------------
144//
145// Returns a TGraphErrors with the mean effective number of photo-electrons divided by
146// the mean charge of that pixel vs. the mean number of photo-electrons.
147//
148TGraphErrors *MCalibrationIntensityChargeCam::GetPhePerCharge( const UInt_t pixid, const MGeomCam &geom, const MCalibrationCam::PulserColor_t col)
149{
150
151 Int_t size = CountNumValidEntries(pixid,col);
152
153 if (size == 0)
154 return NULL;
155
156 TArrayF phepersig(size);
157 TArrayF phepersigerr(size);
158 TArrayF sig(size);
159 TArrayF sigerr(size);
160
161 Int_t cnt = 0;
162
163 for (Int_t i=0;i<GetSize();i++)
164 {
165 //
166 // Get the calibration cam from the intensity cam
167 //
168 MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);
169
170 if (col != MCalibrationCam::kNONE)
171 if (cam->GetPulserColor() != col)
172 continue;
173 //
174 // Get the calibration pix from the calibration cam
175 //
176 const MCalibrationChargePix &pix = (MCalibrationChargePix&)(*cam)[pixid];
177 //
178 // Don't use bad pixels
179 //
180 if (!pix.IsFFactorMethodValid())
181 continue;
182 //
183 // For the calculation of Q, we have to use the
184 // converted value!
185 //
186 const Int_t aidx = geom[pixid].GetAidx();
187 const MCalibrationChargePix &apix = (MCalibrationChargePix&)cam->GetAverageArea(aidx);
188
189 const Float_t q = pix.GetConvertedMean();
190 const Float_t qerr = pix.GetConvertedMeanErr();
191 //
192 const Float_t phe = apix.GetPheFFactorMethod();
193 const Float_t pheerr = apix.GetPheFFactorMethodErr();
194
195 sig[cnt] = phe;
196 sigerr[cnt] = pheerr;
197
198
199 phepersig[cnt] = q > 0.00001 ? phe/q : -1.;
200
201 Float_t var = 0.;
202
203 if (q > 0.00001 && phe > 0.00001)
204 {
205 var = pheerr * pheerr / phe / phe + qerr*qerr/q/q;
206 if (var > 0.00001)
207 var = TMath::Sqrt(var)*phepersig[cnt];
208 }
209 phepersigerr[cnt] = var;
210 cnt++;
211 }
212
213 TGraphErrors *gr = new TGraphErrors(size,
214 sig.GetArray(),phepersig.GetArray(),
215 sigerr.GetArray(),phepersigerr.GetArray());
216 gr->SetTitle(Form("%s%3i","Pixel ",pixid));
217 gr->GetXaxis()->SetTitle("<photo-electrons> [1]");
218 gr->GetYaxis()->SetTitle("<phes> / <Q> [FADC cts^{-1}]");
219 return gr;
220}
221
222// -------------------------------------------------------------------
223//
224// Returns a TGraphErrors with the mean effective number of photo-electrons divided by
225// the mean charge of that pixel vs. the mean number of photo-electrons.
226//
227TGraphErrors *MCalibrationIntensityChargeCam::GetPhePerChargePerArea( const Int_t aidx, const MGeomCam &geom, const MCalibrationCam::PulserColor_t col)
228{
229
230 Int_t size = CountNumEntries(col);
231
232 if (size == 0)
233 return NULL;
234
235 TArrayF phepersig(size);
236 TArrayF phepersigerr(size);
237 TArrayF sig(size);
238 TArrayF sigerr(size);
239
240 Int_t cnt = 0;
241
242 for (Int_t i=0;i<GetSize();i++)
243 {
244 //
245 // Get the calibration cam from the intensity cam
246 //
247 MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);
248
249 if (col != MCalibrationCam::kNONE)
250 if (cam->GetPulserColor() != col)
251 continue;
252 //
253 // Get the calibration pix from the calibration cam
254 //
255 const MCalibrationChargePix &apix = (MCalibrationChargePix&)cam->GetAverageArea(aidx);
256 const Float_t phe = apix.GetPheFFactorMethod();
257 const Float_t pherelvar = apix.GetPheFFactorMethodRelVar();
258 const Float_t pheerr = apix.GetPheFFactorMethodErr();
259
260 sig[cnt] = phe;
261 sigerr[cnt] = pheerr;
262
263 Double_t sig = 0.;
264 Double_t sig2 = 0.;
265 Int_t num = 0;
266
267 for (Int_t i=0; i<cam->GetSize(); i++)
268 {
269 const MCalibrationChargePix &pix = (MCalibrationChargePix&)(*cam)[i];
270 //
271 // Don't use bad pixels
272 //
273 if (!pix.IsFFactorMethodValid())
274 continue;
275 //
276 //
277 if (aidx != geom[i].GetAidx())
278 continue;
279
280 sig += pix.GetConvertedMean();
281 sig2 += pix.GetConvertedMean() * pix.GetConvertedMean();
282 num++;
283 }
284
285 if (num > 1)
286 {
287 sig /= num;
288
289 Double_t var = (sig2 - sig*sig*num) / (num-1);
290 var /= sig*sig;
291 var += pherelvar;
292
293 phepersig[cnt] = phe/sig;
294 if (var > 0.)
295 phepersigerr[cnt] = TMath::Sqrt(var) * phepersig[cnt];
296 else
297 phepersigerr[cnt] = 0.;
298 }
299 else
300 {
301 phepersig[cnt] = -1.;
302 phepersigerr[cnt] = 0.;
303 }
304 cnt++;
305 }
306
307 TGraphErrors *gr = new TGraphErrors(size,
308 sig.GetArray(),phepersig.GetArray(),
309 sigerr.GetArray(),phepersigerr.GetArray());
310 gr->SetTitle(Form("%s%3i","Conv. Factors Area %d Average",aidx));
311 gr->GetXaxis()->SetTitle("<photo-electrons> [1]");
312 gr->GetYaxis()->SetTitle("<phes> / <Q> [FADC cts^{-1}]");
313 return gr;
314}
315
316// -------------------------------------------------------------------
317//
318// Returns a TGraphErrors with the number of photo-electrons vs.
319// the extracted signal over all pixels with area index "aidx".
320//
321// The points represent the means of the pixels values, while the error bars
322// the sigma of the pixels values.
323//
324TGraphErrors *MCalibrationIntensityChargeCam::GetPheVsChargePerArea( const Int_t aidx, const MCalibrationCam::PulserColor_t col)
325{
326
327 Int_t size = CountNumEntries(col);
328
329 TArrayF phe(size);
330 TArrayF pheerr(size);
331 TArrayF sig(size);
332 TArrayF sigerr(size);
333
334 Int_t cnt = 0;
335
336 for (Int_t i=0;i<GetSize();i++)
337 {
338 //
339 // Get the calibration cam from the intensity cam
340 //
341 MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);
342
343 if (col != MCalibrationCam::kNONE)
344 if (cam->GetPulserColor() != col)
345 continue;
346
347 //
348 // Get the area calibration pix from the calibration cam
349 //
350 MCalibrationChargePix &pix = (MCalibrationChargePix&)(cam->GetAverageArea(aidx));
351
352 phe[cnt] = pix.GetPheFFactorMethod();
353 pheerr[cnt] = pix.GetPheFFactorMethodErr();
354 //
355 // For the calculation of Q, we have to use the
356 // converted value!
357 //
358 sig [cnt] = pix.GetConvertedMean();
359 sigerr[cnt] = pix.GetConvertedMeanErr();
360
361 cnt++;
362 }
363
364 TGraphErrors *gr = new TGraphErrors(size,
365 sig.GetArray(),phe.GetArray(),
366 sigerr.GetArray(),pheerr.GetArray());
367 gr->SetTitle(Form("%s%3i","Area Index ",aidx));
368 gr->GetXaxis()->SetTitle("Q [FADC counts]");
369 gr->GetYaxis()->SetTitle("photo-electrons [1]");
370 return gr;
371}
372
373// -------------------------------------------------------------------
374//
375// Returns a TGraphErrors with the 'Razmik plot' of pixel "pixid".
376// The Razmik plot shows the value of 'R' vs. 1/Q where:
377//
378// sigma^2 F^2
379// R = ------- = ------
380// <Q>^2 <m_pe>
381//
382// and 1/Q is the inverse (mean) extracted signal
383//
384TGraphErrors *MCalibrationIntensityChargeCam::GetRazmikPlot( const UInt_t pixid )
385{
386
387 const Int_t size = GetSize();
388
389 TArrayF r(size);
390 TArrayF rerr(size);
391 TArrayF oneoverq(size);
392 TArrayF oneoverqerr(size);
393
394 for (Int_t i=0;i<size;i++)
395 {
396 //
397 // Get the calibration cam from the intensity cam
398 //
399 MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);
400 //
401 // Get the calibration pix from the calibration cam
402 //
403 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*cam)[pixid];
404 //
405 // Don't use bad pixels
406 //
407 if (!pix.IsFFactorMethodValid())
408 continue;
409 //
410 // For the calculation of R, use the un-converted values, like
411 // in the calibration, since:
412 // C^2*sigma^2 sigma^2
413 // R(lowgain) = ----------- = ------ = R
414 // C^2*<Q>^2 <Q>^2
415 //
416 const Float_t mean = pix.GetMean();
417 const Float_t meanerr = pix.GetMeanErr();
418 const Float_t rsigma = pix.GetRSigma();
419 const Float_t rsigmaerr = pix.GetRSigmaErr();
420 r[i] = rsigma*rsigma/mean/mean;
421 const Float_t rrelvar = 4.*rsigmaerr*rsigmaerr/rsigma/rsigma + 4.*meanerr*meanerr/mean/mean;
422 rerr[i] = rrelvar * r[i] * r[i];
423 rerr[i] = rerr[i] <= 0 ? 0. : TMath::Sqrt(rerr[i]);
424 //
425 // For the calculation of 1/Q, we have to use the
426 // converted value!
427 //
428 const Float_t q = pix.GetConvertedMean();
429 const Float_t qe = pix.GetConvertedMeanErr();
430 oneoverq [i] = 1./q;
431 oneoverqerr[i] = qe / (q * q);
432 }
433
434 TGraphErrors *gr = new TGraphErrors(size,
435 oneoverq.GetArray(),r.GetArray(),
436 oneoverqerr.GetArray(),rerr.GetArray());
437 gr->SetTitle(Form("%s%3i","Pixel ",pixid));
438 gr->GetXaxis()->SetTitle("1/Q [FADC counts^{-1}]");
439 gr->GetYaxis()->SetTitle("\\sigma_{red}^{2}/Q^{2} [1]");
440 return gr;
441}
442
443// -------------------------------------------------------------------
444//
445// Returns a 2-dimensional histogram with the fit results of the
446// 'Razmik plot' for each pixel of area index "aidx" (see GetRazmikPlot())
447//
448// The results of the polynomial fit of grade 1 are:
449//
450// x-axis: Offset (Parameter 0 of the polynomial)
451// y-axis: Slope (Parameter 1 of the polynomial)
452//
453// The offset is a measure of how well-known the supposed additional contributions
454// to the value "reduced sigma" are. Because a photo-multiplier is a linear instrument,
455// the excess fluctuations are linear w.r.t. the signal amplitude and can be expressed by
456// the proportionality constant F (the "F-Factor").
457// Adding noise from outside (e.g. night sky background) modifies the recorded noise, but
458// not the mean extracted signal, due to the AC-coupling. Thus, noise contributions from outside
459// (e.g. calculating the pedestal RMS)have to be subtracted from the recorded signal fluctuations
460// in order to retrieve the linearity relation:
461//
462// sigma(signal)^2 / mean(signal)^2 = sigma^2 / <Q>^2 = F^2 / <n_phe> (1)
463//
464// Any systematic offset in the sigma(signal) will produce an offset in the "Razmik plot"),
465// characterized by the Offset of the polynomial fit. Thus, in an ideal case, all pixels have their
466// "offset" centered very closely around zero.
467//
468// The "slope" is the proportionality constant F^2, multiplied with the conversion factor
469// phe's to mean signal (because the "Razmik plot" plots the left side of eq. (1) w.r.t.
470// 1/<Q> instead of 1/<n_phe>. However, the mean number of photo-electrons <n_phe> can be
471// expressed by <Q> with the relation:
472//
473// <n_phe> = c_phe * <Q> (2)
474//
475// Thus:
476//
477// 1/<n_phe> = 1/c_phe * 1/<Q> (3)
478//
479// and:
480//
481// Slope = F^2 / c_phe
482//
483// In the ideal case of having equal photo-multipliers and a perfectly flat-fielded camera,
484// the "slope" -values should thus all be closely centered around F^2/c_phe.
485//
486TH2F *MCalibrationIntensityChargeCam::GetRazmikPlotResults( const Int_t aidx, const MGeomCam &geom)
487{
488
489 TH2F *hist = new TH2F("hist","R vs. Inverse Charges - Fit results",45,-0.02,0.02,45,0.,30.);
490 hist->SetXTitle("Offset [FADC counts^{-1}]");
491 hist->SetYTitle("F^{2} / <n_phe>/<Q> [FADC count / phe]");
492 hist->SetFillColor(kRed+aidx);
493
494 MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam();
495
496 for (Int_t npix=0;npix<cam->GetSize();npix++)
497 {
498
499 if (geom[npix].GetAidx() == aidx)
500 {
501 TGraph *gr = GetRazmikPlot(npix);
502 gr->Fit("pol1","Q");
503 hist->Fill(gr->GetFunction("pol1")->GetParameter(0),gr->GetFunction("pol1")->GetParameter(1));
504 }
505 }
506 return hist;
507}
508
509
510// --------------------------------------------------------------------
511//
512// Returns the number of camera entries matching the required colour
513// and the requirement that pixel "pixid" has been correctly calibrated
514//
515Int_t MCalibrationIntensityChargeCam::CountNumValidEntries(const UInt_t pixid, const MCalibrationCam::PulserColor_t col) const
516{
517
518 Int_t nvalid = 0;
519
520 for (Int_t i=0;i<GetSize();i++)
521 {
522 const MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);
523 const MCalibrationChargePix &pix = (MCalibrationChargePix&)(*cam)[pixid];
524
525 if (col == MCalibrationCam::kNONE)
526 {
527 if (pix.IsFFactorMethodValid())
528 nvalid++;
529 }
530 else
531 {
532 if (cam->GetPulserColor() == col)
533 {
534 if (pix.IsFFactorMethodValid())
535 nvalid++;
536 }
537 }
538 }
539
540 return nvalid;
541}
Note: See TracBrowser for help on using the repository browser.