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

Last change on this file since 8462 was 8106, checked in by tbretz, 19 years ago
*** empty log message ***
File size: 37.1 KB
Line 
1/* ======================================================================== *\
2! $Name: not supported by cvs2svn $:$Id: MCalibrationIntensityChargeCam.cc,v 1.23 2006-10-17 17:15:59 tbretz Exp $
3! --------------------------------------------------------------------------
4!
5! *
6! * This file is part of MARS, the MAGIC Analysis and Reconstruction
7! * Software. It is distributed to you in the hope that it can be a useful
8! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
9! * It is distributed WITHOUT ANY WARRANTY.
10! *
11! * Permission to use, copy, modify and distribute this software and its
12! * documentation for any purpose is hereby granted without fee,
13! * provided that the above copyright notice appear in all copies and
14! * that both that copyright notice and this permission notice appear
15! * in supporting documentation. It is provided "as is" without express
16! * or implied warranty.
17! *
18!
19!
20! Author(s): Markus Gaug 11/2003 <mailto:markus@ifae.es>
21!
22! Copyright: MAGIC Software Development, 2000-2004
23!
24!
25\* ======================================================================== */
26/////////////////////////////////////////////////////////////////////////////
27//
28// MCalibrationIntensityChargeCam
29//
30// Storage container for intensity charge calibration results.
31//
32// Individual MCalibrationChargeCam's can be retrieved with:
33// - GetCam() yielding the current cam.
34// - GetCam("name") yielding the current camera with name "name".
35// - GetCam(i) yielding the i-th camera.
36//
37// See also: MCalibrationIntensityCam, MCalibrationChargeCam,
38// MCalibrationChargePix, MCalibrationChargeCalc, MCalibrationQECam
39// MCalibrationBlindCam, MCalibrationChargePINDiode
40// MHCalibrationChargePix, MHCalibrationChargeCam
41//
42/////////////////////////////////////////////////////////////////////////////
43#include "MCalibrationIntensityChargeCam.h"
44
45#include <TF1.h>
46#include <TH2.h>
47#include <TGraphErrors.h>
48#include <TOrdCollection.h>
49#include <TH1.h>
50
51#include "MLog.h"
52#include "MLogManip.h"
53
54#include "MHCamera.h"
55
56#include "MGeomCamMagic.h"
57#include "MGeomCam.h"
58#include "MGeomPix.h"
59
60#include "MCalibrationChargeCam.h"
61#include "MCalibrationChargePix.h"
62
63ClassImp(MCalibrationIntensityChargeCam);
64
65using namespace std;
66
67// --------------------------------------------------------------------------
68//
69// Default constructor.
70//
71MCalibrationIntensityChargeCam::MCalibrationIntensityChargeCam(const char *name, const char *title)
72{
73
74 fName = name ? name : "MCalibrationIntensityChargeCam";
75 fTitle = title ? title : "Results of the Intensity Calibration";
76
77 InitSize(1);
78}
79
80// -------------------------------------------------------------------
81//
82// Add MCalibrationChargeCam's in the ranges from - to.
83//
84void MCalibrationIntensityChargeCam::Add(const UInt_t from, const UInt_t to)
85{
86 for (UInt_t i=from; i<to; i++)
87 fCams->AddAt(new MCalibrationChargeCam,i);
88}
89
90// -------------------------------------------------------------------
91//
92// Returns a TGraphErrors with the number of photo-electrons vs.
93// the extracted signal of pixel "pixid".
94//
95TGraphErrors *MCalibrationIntensityChargeCam::GetPheVsCharge( const UInt_t pixid, const MCalibrationCam::PulserColor_t col)
96{
97
98 Int_t size = CountNumEntries(col);
99
100 if (size == 0)
101 return NULL;
102
103 TArrayF phe(size);
104 TArrayF pheerr(size);
105 TArrayF sig(size);
106 TArrayF sigerr(size);
107
108 Int_t cnt = 0;
109
110 for (Int_t i=0;i<GetSize();i++)
111 {
112 //
113 // Get the calibration cam from the intensity cam
114 //
115 MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);
116
117 if (col != MCalibrationCam::kNONE)
118 if (cam->GetPulserColor() != col)
119 continue;
120 //
121 // Get the calibration pix from the calibration cam
122 //
123 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*cam)[pixid];
124 //
125 // Don't use bad pixels
126 //
127 if (!pix.IsFFactorMethodValid())
128 continue;
129 //
130 phe[cnt] = pix.GetPheFFactorMethod();
131 pheerr[cnt] = pix.GetPheFFactorMethodErr();
132 //
133 // For the calculation of Q, we have to use the
134 // converted value!
135 //
136 sig [cnt] = pix.GetConvertedMean();
137 sigerr[cnt] = pix.GetConvertedMeanErr();
138 cnt++;
139 }
140
141 TGraphErrors *gr = new TGraphErrors(size,
142 sig.GetArray(),phe.GetArray(),
143 sigerr.GetArray(),pheerr.GetArray());
144 gr->SetTitle(Form("%s%3i","Pixel ",pixid));
145 gr->GetXaxis()->SetTitle("Q [FADC counts]");
146 gr->GetYaxis()->SetTitle("photo-electrons [1]");
147 return gr;
148}
149
150// -------------------------------------------------------------------
151//
152// Returns a TGraphErrors with the mean effective number of photo-electrons divided by
153// the mean charge of that pixel vs. the mean number of photo-electrons.
154//
155TGraphErrors *MCalibrationIntensityChargeCam::GetPhePerCharge( const UInt_t pixid, const MGeomCam &geom, const MCalibrationCam::PulserColor_t col)
156{
157
158 Int_t size = CountNumValidEntries(pixid,col);
159
160 if (size == 0)
161 return NULL;
162
163 TArrayF phepersig(size);
164 TArrayF phepersigerr(size);
165 TArrayF sig(size);
166 TArrayF sigerr(size);
167
168 Int_t cnt = 0;
169
170 for (Int_t i=0;i<GetSize();i++)
171 {
172 //
173 // Get the calibration cam from the intensity cam
174 //
175 MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);
176
177 if (col != MCalibrationCam::kNONE)
178 if (cam->GetPulserColor() != col)
179 continue;
180 //
181 // Get the calibration pix from the calibration cam
182 //
183 const MCalibrationChargePix &pix = (MCalibrationChargePix&)(*cam)[pixid];
184 //
185 // Don't use bad pixels
186 //
187 if (!pix.IsFFactorMethodValid())
188 continue;
189 //
190 // For the calculation of Q, we have to use the
191 // converted value!
192 //
193 const Int_t aidx = geom[pixid].GetAidx();
194 const MCalibrationChargePix &apix = (MCalibrationChargePix&)cam->GetAverageArea(aidx);
195
196 const Float_t q = pix.GetConvertedMean();
197 const Float_t qerr = pix.GetConvertedMeanErr();
198 //
199 const Float_t phe = apix.GetPheFFactorMethod();
200 const Float_t pheerr = apix.GetPheFFactorMethodErr();
201
202 sig[cnt] = phe;
203 sigerr[cnt] = pheerr;
204
205
206 phepersig[cnt] = q > 0.00001 ? phe/q : -1.;
207
208 Float_t var = 0.;
209
210 if (q > 0.00001 && phe > 0.00001)
211 {
212 var = pheerr * pheerr / phe / phe + qerr*qerr/q/q;
213 if (var > 0.00001)
214 var = TMath::Sqrt(var)*phepersig[cnt];
215 }
216 phepersigerr[cnt] = var;
217 cnt++;
218 }
219
220 TGraphErrors *gr = new TGraphErrors(size,
221 sig.GetArray(),phepersig.GetArray(),
222 sigerr.GetArray(),phepersigerr.GetArray());
223 gr->SetTitle(Form("%s%3i","Pixel ",pixid));
224 gr->GetXaxis()->SetTitle("<photo-electrons> [1]");
225 gr->GetYaxis()->SetTitle("<phes> / <Q> [FADC cts^{-1}]");
226 return gr;
227}
228
229// -------------------------------------------------------------------
230//
231// Returns a TGraphErrors with the mean effective number of photo-electrons divided by
232// the mean charge of that pixel vs. the mean number of photo-electrons.
233//
234TGraphErrors *MCalibrationIntensityChargeCam::GetPhePerChargePerArea( const Int_t aidx, const MGeomCam &geom, const MCalibrationCam::PulserColor_t col)
235{
236
237 Int_t size = CountNumEntries(col);
238
239 if (size == 0)
240 return NULL;
241
242 TArrayF phepersig(size);
243 TArrayF phepersigerr(size);
244 TArrayF sig(size);
245 TArrayF sigerr(size);
246
247 Int_t cnt = 0;
248
249 for (Int_t i=0;i<GetSize();i++)
250 {
251 //
252 // Get the calibration cam from the intensity cam
253 //
254 MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);
255
256 if (col != MCalibrationCam::kNONE)
257 if (cam->GetPulserColor() != col)
258 continue;
259 //
260 // Get the calibration pix from the calibration cam
261 //
262 const MCalibrationChargePix &apix = (MCalibrationChargePix&)cam->GetAverageArea(aidx);
263 const Float_t phe = apix.GetPheFFactorMethod();
264 const Float_t pherelvar = apix.GetPheFFactorMethodRelVar();
265 const Float_t pheerr = apix.GetPheFFactorMethodErr();
266
267 sig[cnt] = phe;
268 sigerr[cnt] = pheerr;
269
270 Double_t sig1 = 0.;
271 Double_t sig2 = 0.;
272 Int_t num = 0;
273
274 for (Int_t j=0; j<cam->GetSize(); j++)
275 {
276 const MCalibrationChargePix &pix = (MCalibrationChargePix&)(*cam)[j];
277 //
278 // Don't use bad pixels
279 //
280 if (!pix.IsFFactorMethodValid())
281 continue;
282 //
283 //
284 if (aidx != geom[j].GetAidx())
285 continue;
286
287 sig1 += pix.GetConvertedMean();
288 sig2 += pix.GetConvertedMean() * pix.GetConvertedMean();
289 num++;
290 }
291
292 if (num > 1)
293 {
294 sig1 /= num;
295
296 Double_t var = (sig2 - sig1*sig1*num) / (num-1);
297 var /= sig1*sig1;
298 var += pherelvar;
299
300 phepersig[cnt] = phe/sig1;
301 if (var > 0.)
302 phepersigerr[cnt] = TMath::Sqrt(var) * phepersig[cnt];
303 else
304 phepersigerr[cnt] = 0.;
305 }
306 else
307 {
308 phepersig[cnt] = -1.;
309 phepersigerr[cnt] = 0.;
310 }
311 cnt++;
312 }
313
314 TGraphErrors *gr = new TGraphErrors(size,
315 sig.GetArray(),phepersig.GetArray(),
316 sigerr.GetArray(),phepersigerr.GetArray());
317 gr->SetTitle(Form("Conv. Factors Area %d Average",aidx));
318 gr->GetXaxis()->SetTitle("<photo-electrons> [1]");
319 gr->GetYaxis()->SetTitle("<phes> / <Q> [FADC cts^{-1}]");
320 return gr;
321}
322
323// -------------------------------------------------------------------
324//
325// Returns a TGraphErrors with the number of photo-electrons vs.
326// the extracted signal over all pixels with area index "aidx".
327//
328// The points represent the means of the pixels values, while the error bars
329// the sigma of the pixels values.
330//
331TGraphErrors *MCalibrationIntensityChargeCam::GetPheVsChargePerArea( const Int_t aidx, const MCalibrationCam::PulserColor_t col)
332{
333
334 Int_t size = CountNumEntries(col);
335
336 TArrayF phe(size);
337 TArrayF pheerr(size);
338 TArrayF sig(size);
339 TArrayF sigerr(size);
340
341 Int_t cnt = 0;
342
343 for (Int_t i=0;i<GetSize();i++)
344 {
345 //
346 // Get the calibration cam from the intensity cam
347 //
348 MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);
349
350 if (col != MCalibrationCam::kNONE)
351 if (cam->GetPulserColor() != col)
352 continue;
353
354 //
355 // Get the area calibration pix from the calibration cam
356 //
357 MCalibrationChargePix &pix = (MCalibrationChargePix&)(cam->GetAverageArea(aidx));
358
359 phe[cnt] = pix.GetPheFFactorMethod();
360 pheerr[cnt] = pix.GetPheFFactorMethodErr();
361 //
362 // For the calculation of Q, we have to use the
363 // converted value!
364 //
365 sig [cnt] = pix.GetConvertedMean();
366 sigerr[cnt] = pix.GetConvertedMeanErr();
367
368 cnt++;
369 }
370
371 TGraphErrors *gr = new TGraphErrors(size,
372 sig.GetArray(),phe.GetArray(),
373 sigerr.GetArray(),pheerr.GetArray());
374 gr->SetTitle(Form("%s%3i","Area Index ",aidx));
375 gr->GetXaxis()->SetTitle("Q [FADC counts]");
376 gr->GetYaxis()->SetTitle("photo-electrons [1]");
377 return gr;
378}
379
380// -------------------------------------------------------------------
381//
382// Returns a TGraphErrors with the 'Razmik plot' of pixel "pixid".
383// The Razmik plot shows the value of 'R' vs. 1/Q where:
384//
385// sigma^2 F^2
386// R = ------- = ------
387// <Q>^2 <m_pe>
388//
389// and 1/Q is the inverse (mean) extracted signal
390//
391TGraphErrors *MCalibrationIntensityChargeCam::GetRazmikPlot( const UInt_t pixid )
392{
393
394 const Int_t size = GetSize();
395
396 TArrayF r(size);
397 TArrayF rerr(size);
398 TArrayF oneoverq(size);
399 TArrayF oneoverqerr(size);
400
401 for (Int_t i=0;i<size;i++)
402 {
403 //
404 // Get the calibration cam from the intensity cam
405 //
406 MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);
407 //
408 // Get the calibration pix from the calibration cam
409 //
410 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*cam)[pixid];
411 //
412 // Don't use bad pixels
413 //
414 if (!pix.IsFFactorMethodValid())
415 continue;
416 //
417 // For the calculation of R, use the un-converted values, like
418 // in the calibration, since:
419 // C^2*sigma^2 sigma^2
420 // R(lowgain) = ----------- = ------ = R
421 // C^2*<Q>^2 <Q>^2
422 //
423 const Float_t mean = pix.GetMean();
424 const Float_t meanerr = pix.GetMeanErr();
425 const Float_t rsigma = pix.GetRSigma();
426 const Float_t rsigmaerr = pix.GetRSigmaErr();
427 r[i] = rsigma*rsigma/mean/mean;
428 const Float_t rrelvar = 4.*rsigmaerr*rsigmaerr/rsigma/rsigma + 4.*meanerr*meanerr/mean/mean;
429 rerr[i] = rrelvar * r[i] * r[i];
430 rerr[i] = rerr[i] <= 0 ? 0. : TMath::Sqrt(rerr[i]);
431 //
432 // For the calculation of 1/Q, we have to use the
433 // converted value!
434 //
435 const Float_t q = pix.GetConvertedMean();
436 const Float_t qe = pix.GetConvertedMeanErr();
437 oneoverq [i] = 1./q;
438 oneoverqerr[i] = qe / (q * q);
439 }
440
441 TGraphErrors *gr = new TGraphErrors(size,
442 oneoverq.GetArray(),r.GetArray(),
443 oneoverqerr.GetArray(),rerr.GetArray());
444 gr->SetTitle(Form("%s%3i","Pixel ",pixid));
445 gr->GetXaxis()->SetTitle("1/Q [FADC counts^{-1}]");
446 gr->GetYaxis()->SetTitle("\\sigma_{red}^{2}/Q^{2} [1]");
447 return gr;
448}
449
450// -------------------------------------------------------------------
451//
452// Returns a 2-dimensional histogram with the fit results of the
453// 'Razmik plot' for each pixel of area index "aidx" (see GetRazmikPlot())
454//
455// The results of the polynomial fit of grade 1 are:
456//
457// x-axis: Offset (Parameter 0 of the polynomial)
458// y-axis: Slope (Parameter 1 of the polynomial)
459//
460// The offset is a measure of how well-known the supposed additional contributions
461// to the value "reduced sigma" are. Because a photo-multiplier is a linear instrument,
462// the excess fluctuations are linear w.r.t. the signal amplitude and can be expressed by
463// the proportionality constant F (the "F-Factor").
464// Adding noise from outside (e.g. night sky background) modifies the recorded noise, but
465// not the mean extracted signal, due to the AC-coupling. Thus, noise contributions from outside
466// (e.g. calculating the pedestal RMS)have to be subtracted from the recorded signal fluctuations
467// in order to retrieve the linearity relation:
468//
469// sigma(signal)^2 / mean(signal)^2 = sigma^2 / <Q>^2 = F^2 / <n_phe> (1)
470//
471// Any systematic offset in the sigma(signal) will produce an offset in the "Razmik plot"),
472// characterized by the Offset of the polynomial fit. Thus, in an ideal case, all pixels have their
473// "offset" centered very closely around zero.
474//
475// The "slope" is the proportionality constant F^2, multiplied with the conversion factor
476// phe's to mean signal (because the "Razmik plot" plots the left side of eq. (1) w.r.t.
477// 1/<Q> instead of 1/<n_phe>. However, the mean number of photo-electrons <n_phe> can be
478// expressed by <Q> with the relation:
479//
480// <n_phe> = c_phe * <Q> (2)
481//
482// Thus:
483//
484// 1/<n_phe> = 1/c_phe * 1/<Q> (3)
485//
486// and:
487//
488// Slope = F^2 / c_phe
489//
490// In the ideal case of having equal photo-multipliers and a perfectly flat-fielded camera,
491// the "slope" -values should thus all be closely centered around F^2/c_phe.
492//
493TH2F *MCalibrationIntensityChargeCam::GetRazmikPlotResults( const Int_t aidx, const MGeomCam &geom)
494{
495
496 TH2F *hist = new TH2F("hist","R vs. Inverse Charges - Fit results",45,-0.02,0.02,45,0.,30.);
497 hist->SetXTitle("Offset [FADC counts^{-1}]");
498 hist->SetYTitle("F^{2} / <n_phe>/<Q> [FADC count / phe]");
499 hist->SetFillColor(kRed+aidx);
500
501 MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam();
502
503 for (Int_t npix=0;npix<cam->GetSize();npix++)
504 {
505
506 if (geom[npix].GetAidx() == aidx)
507 {
508 TGraph *gr = GetRazmikPlot(npix);
509 gr->Fit("pol1","Q");
510 hist->Fill(gr->GetFunction("pol1")->GetParameter(0),gr->GetFunction("pol1")->GetParameter(1));
511 }
512 }
513 return hist;
514}
515
516
517// --------------------------------------------------------------------
518//
519// Returns the number of camera entries matching the required colour
520// and the requirement that pixel "pixid" has been correctly calibrated
521//
522Int_t MCalibrationIntensityChargeCam::CountNumValidEntries(const UInt_t pixid, const MCalibrationCam::PulserColor_t col) const
523{
524
525 Int_t nvalid = 0;
526
527 for (Int_t i=0;i<GetSize();i++)
528 {
529 const MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);
530 const MCalibrationChargePix &pix = (MCalibrationChargePix&)(*cam)[pixid];
531
532 if (col == MCalibrationCam::kNONE)
533 {
534 if (pix.IsFFactorMethodValid())
535 nvalid++;
536 }
537 else
538 {
539 if (cam->GetPulserColor() == col)
540 {
541 if (pix.IsFFactorMethodValid())
542 nvalid++;
543 }
544 }
545 }
546
547 return nvalid;
548}
549
550
551// -------------------------------------------------------------------
552//
553// Returns a TGraphErrors with the development of the number of
554// photo-electrons vs. camera number for pixel 'pixid'
555//
556TGraphErrors *MCalibrationIntensityChargeCam::GetVarVsTime( const Int_t pixid , const Option_t *varname )
557{
558
559 const Int_t size = GetSize();
560
561 if (size == 0)
562 return NULL;
563
564 TString option(varname);
565 option.ToLower();
566
567 TArrayF nr(size);
568 TArrayF nrerr(size);
569 TArrayF var (size);
570 TArrayF varerr(size);
571
572 for (Int_t i=0;i<size;i++)
573 {
574 //
575 // Get the calibration cam from the intensity cam
576 //
577 MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);
578 //
579 // Get the calibration pix from the calibration cam
580 //
581 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*cam)[pixid];
582 //
583 nr[i] = i;
584 nrerr[i] = 0.;
585 var[i] = -1.;
586 varerr[i] = -1.;
587 //
588 // Don't use bad pixels
589 //
590 if (!pix.IsFFactorMethodValid())
591 continue;
592 //
593 if (option.Contains("rsigma"))
594 {
595 var [i] = pix.GetRSigma();
596 varerr[i] = pix.GetRSigmaErr();
597 }
598 if (option.Contains("abstimemean"))
599 {
600 var [i] = pix.GetAbsTimeMean();
601 varerr[i] = pix.GetAbsTimeRms();
602 }
603 if (option.Contains("abstimerms"))
604 {
605 var [i] = pix.GetAbsTimeRms();
606 varerr[i] = pix.GetAbsTimeRms()/2.;
607 }
608 if (option.Contains("blackout"))
609 {
610 var [i] = pix.GetNumBlackout();
611 varerr[i] = 0.;
612 }
613 if (option.Contains("pickup"))
614 {
615 var [i] = pix.GetNumPickup();
616 varerr[i] = 0.;
617 }
618 if (option.Contains("outlier"))
619 {
620 var [i] = pix.GetNumPickup() + pix.GetNumBlackout();
621 varerr[i] = 0.;
622 }
623 if (option.Contains("conversionhilo"))
624 {
625 var [i] = pix.GetConversionHiLo();
626 varerr[i] = pix.GetConversionHiLoErr();
627 }
628 if (option.Contains("convertedmean"))
629 {
630 var [i] = pix.GetConvertedMean();
631 varerr[i] = pix.GetConvertedMeanErr();
632 }
633 if (option.Contains("convertedsigma"))
634 {
635 var [i] = pix.GetConvertedSigma();
636 varerr[i] = pix.GetConvertedSigmaErr();
637 }
638 if (option.Contains("convertedrsigma"))
639 {
640 var [i] = pix.GetConvertedRSigma();
641 varerr[i] = pix.GetConvertedRSigmaErr();
642 }
643 if (option.Contains("meanconvfadc2phe"))
644 {
645 var [i] = pix.GetMeanConvFADC2Phe();
646 varerr[i] = pix.GetMeanConvFADC2PheErr();
647 }
648 if (option.Contains("meanffactorfadc2phot"))
649 {
650 var [i] = pix.GetMeanFFactorFADC2Phot();
651 varerr[i] = pix.GetMeanFFactorFADC2PhotErr();
652 }
653 if (option.Contains("ped"))
654 {
655 var [i] = pix.GetPed();
656 varerr[i] = pix.GetPedErr();
657 }
658 if (option.Contains("pedrms"))
659 {
660 var [i] = pix.GetPedRms();
661 varerr[i] = pix.GetPedRmsErr();
662 }
663 if (option.Contains("pheffactormethod"))
664 {
665 var [i] = pix.GetPheFFactorMethod();
666 varerr[i] = pix.GetPheFFactorMethodErr();
667 }
668 if (option.Contains("rsigmapercharge"))
669 {
670 var [i] = pix.GetRSigmaPerCharge();
671 varerr[i] = pix.GetRSigmaPerChargeErr();
672 }
673 if (option.Contains("conversionfactor"))
674 {
675 const MCalibrationChargePix &apix = (MCalibrationChargePix&)cam->GetAverageArea(0);
676 const Float_t mean = pix.GetConvertedMean();
677 const Float_t phe = apix.GetPheFFactorMethod();
678
679 var[i] = phe/mean;
680 varerr[i] = TMath::Sqrt(apix.GetPheFFactorMethodErr()*apix.GetPheFFactorMethodErr()/mean/mean
681 + phe*phe/mean/mean/mean/mean*pix.GetConvertedMeanErr()*pix.GetConvertedMeanErr());
682 }
683 }
684
685
686 TGraphErrors *gr = new TGraphErrors(size,
687 nr.GetArray(),var.GetArray(),
688 nrerr.GetArray(),varerr.GetArray());
689 gr->SetTitle(Form("%s%3i","Pixel ",pixid));
690 gr->GetXaxis()->SetTitle("Camera Nr.");
691 // gr->GetYaxis()->SetTitle("<Q> [FADC cnts]");
692 return gr;
693}
694
695// --------------------------------------------------------------------------------
696//
697// Returns a TGraphErrors with a pre-defined variable with name (handed over in 'opt')
698// per area index 'aidx' vs. the calibration camera number
699//
700TGraphErrors *MCalibrationIntensityChargeCam::GetVarPerAreaVsTime( const Int_t aidx, const MGeomCam &geom, const Option_t *varname)
701{
702
703 const Int_t size = GetSize();
704
705 if (size == 0)
706 return NULL;
707
708 TString option(varname);
709 option.ToLower();
710
711 TArrayF vararea(size);
712 TArrayF varareaerr(size);
713 TArrayF nr(size);
714 TArrayF nrerr(size);
715
716 TH1D *h = 0;
717
718 for (Int_t i=0;i<GetSize();i++)
719 {
720 //
721 // Get the calibration cam from the intensity cam
722 //
723 MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);
724
725 //
726 // Get the calibration pix from the calibration cam
727 //
728 Double_t variab = 0.;
729 Double_t variab2 = 0.;
730 Double_t variance = 0.;
731 Int_t num = 0;
732 Float_t pvar = 0.;
733
734 MHCamera camcharge(geom,"CamCharge","Variable;;channels");
735 //
736 // Get the area calibration pix from the calibration cam
737 //
738 for (Int_t j=0; j<cam->GetSize(); j++)
739 {
740 const MCalibrationChargePix &pix = (MCalibrationChargePix&)(*cam)[j];
741 //
742 // Don't use bad pixels
743 //
744 if (!pix.IsFFactorMethodValid())
745 continue;
746 //
747 //
748 if (aidx != geom[j].GetAidx())
749 continue;
750
751 pvar = 0.;
752
753 if (option.Contains("rsigma"))
754 pvar = pix.GetRSigma();
755 if (option.Contains("abstimemean"))
756 pvar = pix.GetAbsTimeMean();
757 if (option.Contains("abstimerms"))
758 pvar = pix.GetAbsTimeRms();
759 if (option.Contains("conversionhilo"))
760 pvar = pix.GetConversionHiLo();
761 if (option.Contains("convertedmean"))
762 pvar = pix.GetConvertedMean();
763 if (option.Contains("convertedsigma"))
764 pvar = pix.GetConvertedSigma();
765 if (option.Contains("convertedrsigma"))
766 pvar = pix.GetConvertedRSigma();
767 if (option.Contains("meanconvfadc2phe"))
768 pvar = pix.GetMeanConvFADC2Phe();
769 if (option.Contains("meanffactorfadc2phot"))
770 pvar = pix.GetMeanFFactorFADC2Phot();
771 if (option.Contains("ped"))
772 pvar = pix.GetPed();
773 if (option.Contains("pedrms"))
774 pvar = pix.GetPedRms();
775 if (option.Contains("pheffactormethod"))
776 pvar = pix.GetPheFFactorMethod();
777 if (option.Contains("rsigmapercharge"))
778 pvar = pix.GetRSigmaPerCharge();
779 if (option.Contains("conversionfactor"))
780 {
781 const MCalibrationChargePix &apix = (MCalibrationChargePix&)cam->GetAverageArea(aidx);
782 pvar = apix.GetPheFFactorMethod()/pix.GetConvertedMean();
783 }
784
785
786 variab += pvar;
787 variab2 += pvar*pvar;
788 num++;
789
790 camcharge.Fill(j,pvar);
791 camcharge.SetUsed(j);
792 }
793
794 if (num > 1)
795 {
796 variab /= num;
797 variance = (variab2 - variab*variab*num) / (num-1);
798
799 vararea[i] = variab;
800 varareaerr[i] = variance>0 ? TMath::Sqrt(variance/num) : 999999999.;
801
802 //
803 // Make also a Gauss-fit to the distributions. The RMS can be determined by
804 // outlier, thus we look at the sigma and the RMS and take the smaller one, afterwards.
805 //
806 h = camcharge.ProjectionS(TArrayI(),TArrayI(1,&aidx),"_py",750);
807 h->SetDirectory(NULL);
808 h->Fit("gaus","QL");
809 TF1 *fit = h->GetFunction("gaus");
810
811 Float_t ci2 = fit->GetChisquare();
812 Float_t sigma = fit->GetParameter(2);
813
814 if (ci2 > 500. || sigma > varareaerr[i])
815 {
816 h->Fit("gaus","QLM");
817 fit = h->GetFunction("gaus");
818
819 ci2 = fit->GetChisquare();
820 sigma = fit->GetParameter(2);
821 }
822
823 const Float_t mean = fit->GetParameter(1);
824 const Float_t ndf = fit->GetNDF();
825
826 *fLog << inf << "Camera Nr: " << i << endl;
827 *fLog << inf << option.Data() << " area idx: " << aidx << " Results: " << endl;
828 *fLog << inf << "Mean: " << Form("%4.3f",mean)
829 << "+-" << Form("%4.3f",fit->GetParError(1))
830 << " Sigma: " << Form("%4.3f",sigma) << "+-" << Form("%4.3f",fit->GetParError(2))
831 << " Chisquare: " << Form("%4.3f",ci2) << " NDF : " << ndf << endl;
832 delete h;
833 gROOT->GetListOfFunctions()->Remove(fit);
834
835 if (sigma<varareaerr[i] && ndf>2 && ci2<500.)
836 {
837 vararea [i] = mean;
838 varareaerr[i] = sigma/TMath::Sqrt((Float_t)num);
839 }
840 }
841 else
842 {
843 vararea[i] = -1.;
844 varareaerr[i] = 0.;
845 }
846
847 nr[i] = i;
848 nrerr[i] = 0.;
849 }
850
851 TGraphErrors *gr = new TGraphErrors(size,
852 nr.GetArray(),vararea.GetArray(),
853 nrerr.GetArray(),varareaerr.GetArray());
854 gr->SetTitle(Form("%s Area %3i Average",option.Data(),aidx));
855 gr->GetXaxis()->SetTitle("Camera Nr.");
856 // gr->GetYaxis()->SetTitle("<Q> [1]");
857 return gr;
858}
859
860
861// -------------------------------------------------------------------
862//
863// Returns a TGraphErrors with the mean effective number of photon
864// vs. the calibration camera number. With the string 'method', different
865// calibration methods can be called.
866//
867TGraphErrors *MCalibrationIntensityChargeCam::GetPhotVsTime( const Option_t *method )
868{
869
870 const Int_t size = GetSize();
871
872 if (size == 0)
873 return NULL;
874
875 TString option(method);
876
877 TArrayF photarr(size);
878 TArrayF photarrerr(size);
879 TArrayF nr(size);
880 TArrayF nrerr(size);
881
882 for (Int_t i=0;i<GetSize();i++)
883 {
884 //
885 // Get the calibration cam from the intensity cam
886 //
887 MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);
888
889 //
890 // Get the calibration pix from the calibration cam
891 //
892 Float_t phot = 0.;
893 Float_t photerr = 0.;
894
895 if (option.Contains("BlindPixel"))
896 {
897 phot = cam->GetNumPhotonsBlindPixelMethod();
898 photerr = cam->GetNumPhotonsBlindPixelMethodErr();
899 }
900 if (option.Contains("FFactor"))
901 {
902 phot = cam->GetNumPhotonsFFactorMethod();
903 photerr = cam->GetNumPhotonsFFactorMethodErr();
904 }
905 if (option.Contains("PINDiode"))
906 {
907 phot = cam->GetNumPhotonsPINDiodeMethod();
908 photerr = cam->GetNumPhotonsPINDiodeMethodErr();
909 }
910
911 photarr[i] = phot;
912 photarrerr[i] = photerr;
913
914 nr[i] = i;
915 nrerr[i] = 0.;
916 }
917
918 TGraphErrors *gr = new TGraphErrors(size,
919 nr.GetArray(),photarr.GetArray(),
920 nrerr.GetArray(),photarrerr.GetArray());
921 gr->SetTitle("Photons Average");
922 gr->GetXaxis()->SetTitle("Camera Nr.");
923 gr->GetYaxis()->SetTitle("<N_{phot}> [1]");
924 return gr;
925}
926
927// -------------------------------------------------------------------
928//
929// Returns a TGraphErrors with the mean effective number of photo-electrons per
930// area index 'aidx' vs. the calibration camera number
931//
932TGraphErrors *MCalibrationIntensityChargeCam::GetPhePerAreaVsTime( const Int_t aidx, const MGeomCam &geom)
933{
934
935 const Int_t size = GetSize();
936
937 if (size == 0)
938 return NULL;
939
940 TArrayF phearea(size);
941 TArrayF pheareaerr(size);
942 TArrayF time(size);
943 TArrayF timeerr(size);
944
945 for (Int_t i=0;i<GetSize();i++)
946 {
947 //
948 // Get the calibration cam from the intensity cam
949 //
950 MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);
951
952 //
953 // Get the calibration pix from the calibration cam
954 //
955 const MCalibrationChargePix &apix = (MCalibrationChargePix&)cam->GetAverageArea(aidx);
956 const Float_t phe = apix.GetPheFFactorMethod();
957 const Float_t pheerr = apix.GetPheFFactorMethodErr();
958
959 phearea[i] = phe;
960 pheareaerr[i] = pheerr;
961
962 time[i] = i;
963 timeerr[i] = 0.;
964 }
965
966 TGraphErrors *gr = new TGraphErrors(size,
967 time.GetArray(),phearea.GetArray(),
968 timeerr.GetArray(),pheareaerr.GetArray());
969 gr->SetTitle(Form("Phes Area %d Average",aidx));
970 gr->GetXaxis()->SetTitle("Camera Nr.");
971 gr->GetYaxis()->SetTitle("<N_{phe}> [1]");
972 return gr;
973}
974
975// -------------------------------------------------------------------
976//
977// Returns a TGraphErrors with the event-by-event averaged charge per
978// area index 'aidx' vs. the calibration camera number
979//
980TGraphErrors *MCalibrationIntensityChargeCam::GetChargePerAreaVsTime( const Int_t aidx, const MGeomCam &geom)
981{
982
983 const Int_t size = GetSize();
984
985 if (size == 0)
986 return NULL;
987
988 TArrayF chargearea(size);
989 TArrayF chargeareaerr(size);
990 TArrayF nr(size);
991 TArrayF nrerr(size);
992
993 for (Int_t i=0;i<GetSize();i++)
994 {
995 //
996 // Get the calibration cam from the intensity cam
997 //
998 MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);
999
1000 //
1001 // Get the calibration pix from the calibration cam
1002 //
1003 const MCalibrationChargePix &apix = (MCalibrationChargePix&)cam->GetAverageArea(aidx);
1004 const Float_t charge = apix.GetConvertedMean();
1005 const Float_t chargeerr = apix.GetConvertedSigma();
1006
1007 chargearea[i] = charge;
1008 chargeareaerr[i] = chargeerr;
1009
1010 nr[i] = i;
1011 nrerr[i] = 0.;
1012 }
1013
1014 TGraphErrors *gr = new TGraphErrors(size,
1015 nr.GetArray(),chargearea.GetArray(),
1016 nrerr.GetArray(),chargeareaerr.GetArray());
1017 gr->SetTitle(Form("Averaged Charges Area Idx %d",aidx));
1018 gr->GetXaxis()->SetTitle("Camera Nr.");
1019 gr->GetYaxis()->SetTitle("<Q> [FADC cnts]");
1020 return gr;
1021}
1022
1023TH1F *MCalibrationIntensityChargeCam::GetVarFluctuations( const Int_t aidx, const MGeomCam &geom, const Option_t *varname )
1024{
1025
1026 const Int_t size = GetSize();
1027
1028 if (size == 0)
1029 return NULL;
1030
1031 TString option(varname);
1032 option.ToLower();
1033
1034 TH1F *hist = new TH1F("hist",Form("%s - Rel. Fluctuations %s Pixel",option.Data(),aidx ? "Outer" : "Inner"),
1035 200,0.,100.);
1036 hist->SetXTitle("Relative Fluctuation [%]");
1037 hist->SetYTitle("Nr. channels [1]");
1038 hist->SetFillColor(kRed+aidx);
1039
1040 MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam();
1041
1042 //
1043 // Loop over pixels
1044 //
1045 for (Int_t npix=0;npix<cam->GetSize();npix++)
1046 {
1047 if (geom[npix].GetAidx() != aidx)
1048 continue;
1049
1050 Double_t variab = 0.;
1051 Double_t variab2 = 0.;
1052 Double_t variance = 0.;
1053 Int_t num = 0;
1054 Float_t pvar = 0.;
1055 Float_t relrms = 99.9;
1056 //
1057 // Loop over the Cams for each pixel
1058 //
1059 for (Int_t i=0; i<GetSize(); i++)
1060 {
1061 MCalibrationChargeCam *ccam = (MCalibrationChargeCam*)GetCam(i);
1062 //
1063 // Get the calibration pix from the calibration cam
1064 //
1065 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*ccam)[npix];
1066 //
1067 // Don't use bad pixels
1068 //
1069 if (!pix.IsFFactorMethodValid())
1070 continue;
1071
1072 if (option.Contains("rsigma"))
1073 pvar = pix.GetRSigma();
1074 if (option.Contains("abstimemean"))
1075 pvar = pix.GetAbsTimeMean();
1076 if (option.Contains("abstimerms"))
1077 pvar = pix.GetAbsTimeRms();
1078 if (option.Contains("conversionhilo"))
1079 pvar = pix.GetConversionHiLo();
1080 if (option.Contains("convertedmean"))
1081 pvar = pix.GetConvertedMean();
1082 if (option.Contains("convertedsigma"))
1083 pvar = pix.GetConvertedSigma();
1084 if (option.Contains("convertedrsigma"))
1085 pvar = pix.GetConvertedRSigma();
1086 if (option.Contains("meanconvfadc2phe"))
1087 pvar = pix.GetMeanConvFADC2Phe();
1088 if (option.Contains("meanffactorfadc2phot"))
1089 pvar = pix.GetMeanFFactorFADC2Phot();
1090 if (option.Contains("ped"))
1091 pvar = pix.GetPed();
1092 if (option.Contains("pedrms"))
1093 pvar = pix.GetPedRms();
1094 if (option.Contains("pheffactormethod"))
1095 pvar = pix.GetPheFFactorMethod();
1096 if (option.Contains("rsigmapercharge"))
1097 pvar = pix.GetRSigmaPerCharge();
1098 if (option.Contains("conversionfactor"))
1099 {
1100 const MCalibrationChargePix &apix = (MCalibrationChargePix&)ccam->GetAverageArea(0);
1101 pvar = apix.GetPheFFactorMethod()/pix.GetConvertedMean();
1102 }
1103
1104
1105 variab += pvar;
1106 variab2 += pvar*pvar;
1107 num++;
1108 }
1109
1110 if (num > 1)
1111 {
1112 variab /= num;
1113 variance = (variab2 - variab*variab*num) / (num-1);
1114
1115 if (variance > 0.)
1116 relrms = TMath::Sqrt(variance)/variab * 100.;
1117 }
1118 hist->Fill(relrms);
1119 }
1120 return hist;
1121}
1122
1123void MCalibrationIntensityChargeCam::DrawRazmikPlot( const UInt_t pixid )
1124{
1125 TGraphErrors *gr = GetRazmikPlot(pixid );
1126 gr->SetBit(kCanDelete);
1127 gr->Draw("A*");
1128
1129}
1130void MCalibrationIntensityChargeCam::DrawPheVsCharge( const UInt_t pixid, const MCalibrationCam::PulserColor_t col)
1131{
1132 TGraphErrors *gr = GetPheVsCharge(pixid,col);
1133 gr->SetBit(kCanDelete);
1134 gr->Draw("A*");
1135}
1136void MCalibrationIntensityChargeCam::DrawPhePerCharge( const UInt_t pixid, const MCalibrationCam::PulserColor_t col)
1137{
1138 TGraphErrors *gr = GetPhePerCharge(pixid,MGeomCamMagic(),col);
1139 gr->SetBit(kCanDelete);
1140 gr->Draw("A*");
1141}
1142void MCalibrationIntensityChargeCam::DrawPhePerChargePerArea( const Int_t aidx, const MCalibrationCam::PulserColor_t col)
1143{
1144 TGraphErrors *gr = GetPhePerChargePerArea(aidx,MGeomCamMagic(),col);
1145 gr->SetBit(kCanDelete);
1146 gr->Draw("A*");
1147}
1148void MCalibrationIntensityChargeCam::DrawPheVsChargePerArea( const Int_t aidx, const MCalibrationCam::PulserColor_t col)
1149{
1150 TGraphErrors *gr = GetPheVsChargePerArea(aidx,col);
1151 gr->SetBit(kCanDelete);
1152 gr->Draw("A*");
1153}
1154void MCalibrationIntensityChargeCam::DrawRazmikPlotResults( const Int_t aidx)
1155{
1156 TH2F *h = GetRazmikPlotResults(aidx,MGeomCamMagic());
1157 h->SetBit(kCanDelete);
1158 h->Draw();
1159}
1160
1161void MCalibrationIntensityChargeCam::DrawChargePerAreaVsTime( const Int_t aidx)
1162{
1163 TGraphErrors *gr = GetChargePerAreaVsTime(aidx,MGeomCamMagic());
1164 gr->SetBit(kCanDelete);
1165 gr->Draw("A*");
1166}
1167void MCalibrationIntensityChargeCam::DrawPhePerAreaVsTime( const Int_t aidx)
1168{
1169 TGraphErrors *gr = GetPhePerAreaVsTime(aidx,MGeomCamMagic());
1170 gr->SetBit(kCanDelete);
1171 gr->Draw("A*");
1172}
1173void MCalibrationIntensityChargeCam::DrawPhotVsTime( const Option_t *method)
1174{
1175 TGraphErrors *gr = GetPhotVsTime(method);
1176 gr->SetBit(kCanDelete);
1177 gr->Draw("A*");
1178}
1179
1180void MCalibrationIntensityChargeCam::DrawVarPerAreaVsTime( const Int_t aidx, const Option_t *varname )
1181{
1182 TGraphErrors *gr = GetVarPerAreaVsTime(aidx,MGeomCamMagic(),varname );
1183 gr->SetBit(kCanDelete);
1184 gr->Draw("A*");
1185}
1186void MCalibrationIntensityChargeCam::DrawVarVsTime( const Int_t pixid , const Option_t *varname )
1187{
1188 TGraphErrors *gr = GetVarVsTime(pixid,varname );
1189 gr->SetBit(kCanDelete);
1190 gr->Draw("A*");
1191}
1192void MCalibrationIntensityChargeCam::DrawVarFluctuations( const Int_t aidx, const Option_t *varname)
1193{
1194 TH1F *h = GetVarFluctuations( aidx,MGeomCamMagic(),varname);
1195 h->SetBit(kCanDelete);
1196 h->Draw();
1197}
Note: See TracBrowser for help on using the repository browser.