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

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