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

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