source: tags/Mars-V0.9/mcalib/MCalibrationIntensityChargeCam.cc

Last change on this file was 6412, checked in by gaug, 20 years ago
*** empty log message ***
File size: 35.6 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("abstime"))
597 {
598 var [i] = pix.GetAbsTimeMean();
599 varerr[i] = pix.GetAbsTimeRms();
600 }
601 if (option.Contains("blackout"))
602 {
603 var [i] = pix.GetNumBlackout();
604 varerr[i] = 0.;
605 }
606 if (option.Contains("pickup"))
607 {
608 var [i] = pix.GetNumPickup();
609 varerr[i] = 0.;
610 }
611 if (option.Contains("outlier"))
612 {
613 var [i] = pix.GetNumPickup() + pix.GetNumBlackout();
614 varerr[i] = 0.;
615 }
616 if (option.Contains("conversionhilo"))
617 {
618 var [i] = pix.GetConversionHiLo();
619 varerr[i] = pix.GetConversionHiLoErr();
620 }
621 if (option.Contains("convertedmean"))
622 {
623 var [i] = pix.GetConvertedMean();
624 varerr[i] = pix.GetConvertedMeanErr();
625 }
626 if (option.Contains("convertedsigma"))
627 {
628 var [i] = pix.GetConvertedSigma();
629 varerr[i] = pix.GetConvertedSigmaErr();
630 }
631 if (option.Contains("convertedrsigma"))
632 {
633 var [i] = pix.GetConvertedRSigma();
634 varerr[i] = pix.GetConvertedRSigmaErr();
635 }
636 if (option.Contains("meanconvfadc2phe"))
637 {
638 var [i] = pix.GetMeanConvFADC2Phe();
639 varerr[i] = pix.GetMeanConvFADC2PheErr();
640 }
641 if (option.Contains("meanffactorfadc2phot"))
642 {
643 var [i] = pix.GetMeanFFactorFADC2Phot();
644 varerr[i] = pix.GetMeanFFactorFADC2PhotErr();
645 }
646 if (option.Contains("ped"))
647 {
648 var [i] = pix.GetPed();
649 varerr[i] = pix.GetPedErr();
650 }
651 if (option.Contains("pedrms"))
652 {
653 var [i] = pix.GetPedRms();
654 varerr[i] = pix.GetPedRmsErr();
655 }
656 if (option.Contains("pheffactormethod"))
657 {
658 var [i] = pix.GetPheFFactorMethod();
659 varerr[i] = pix.GetPheFFactorMethodErr();
660 }
661 if (option.Contains("rsigmapercharge"))
662 {
663 var [i] = pix.GetRSigmaPerCharge();
664 varerr[i] = pix.GetRSigmaPerChargeErr();
665 }
666 }
667
668
669 TGraphErrors *gr = new TGraphErrors(size,
670 nr.GetArray(),var.GetArray(),
671 nrerr.GetArray(),varerr.GetArray());
672 gr->SetTitle(Form("%s%3i","Pixel ",pixid));
673 gr->GetXaxis()->SetTitle("Camera Nr.");
674 // gr->GetYaxis()->SetTitle("<Q> [FADC cnts]");
675 return gr;
676}
677
678// --------------------------------------------------------------------------------
679//
680// Returns a TGraphErrors with a pre-defined variable with name (handed over in 'opt')
681// per area index 'aidx' vs. the calibration camera number
682//
683TGraphErrors *MCalibrationIntensityChargeCam::GetVarPerAreaVsTime( const Int_t aidx, const MGeomCam &geom, const Option_t *varname)
684{
685
686 const Int_t size = GetSize();
687
688 if (size == 0)
689 return NULL;
690
691 TString option(varname);
692 option.ToLower();
693
694 TArrayF vararea(size);
695 TArrayF varareaerr(size);
696 TArrayF nr(size);
697 TArrayF nrerr(size);
698
699 TH1D *h = 0;
700
701 for (Int_t i=0;i<GetSize();i++)
702 {
703 //
704 // Get the calibration cam from the intensity cam
705 //
706 MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);
707
708 //
709 // Get the calibration pix from the calibration cam
710 //
711 Double_t variab = 0.;
712 Double_t variab2 = 0.;
713 Double_t variance = 0.;
714 Int_t num = 0;
715 Float_t pvar = 0.;
716
717 MHCamera camcharge(geom,"CamCharge","Variable;;channels");
718 //
719 // Get the area calibration pix from the calibration cam
720 //
721 for (Int_t j=0; j<cam->GetSize(); j++)
722 {
723 const MCalibrationChargePix &pix = (MCalibrationChargePix&)(*cam)[j];
724 //
725 // Don't use bad pixels
726 //
727 if (!pix.IsFFactorMethodValid())
728 continue;
729 //
730 //
731 if (aidx != geom[j].GetAidx())
732 continue;
733
734 pvar = 0.;
735
736 if (option.Contains("rsigma"))
737 pvar = pix.GetRSigma();
738 if (option.Contains("abstime"))
739 pvar = pix.GetAbsTimeMean();
740 if (option.Contains("conversionhilo"))
741 pvar = pix.GetConversionHiLo();
742 if (option.Contains("convertedmean"))
743 pvar = pix.GetConvertedMean();
744 if (option.Contains("convertedsigma"))
745 pvar = pix.GetConvertedSigma();
746 if (option.Contains("convertedrsigma"))
747 pvar = pix.GetConvertedRSigma();
748 if (option.Contains("meanconvfadc2phe"))
749 pvar = pix.GetMeanConvFADC2Phe();
750 if (option.Contains("meanffactorfadc2phot"))
751 pvar = pix.GetMeanFFactorFADC2Phot();
752 if (option.Contains("ped"))
753 pvar = pix.GetPed();
754 if (option.Contains("pedrms"))
755 pvar = pix.GetPedRms();
756 if (option.Contains("pheffactormethod"))
757 pvar = pix.GetPheFFactorMethod();
758 if (option.Contains("rsigmapercharge"))
759 pvar = pix.GetRSigmaPerCharge();
760
761 variab += pvar;
762 variab2 += pvar*pvar;
763 num++;
764
765 camcharge.Fill(j,pvar);
766 camcharge.SetUsed(j);
767 }
768
769 if (num > 1)
770 {
771 variab /= num;
772 variance = (variab2 - variab*variab*num) / (num-1);
773
774 vararea[i] = variab;
775 if (variance > 0.)
776 varareaerr[i] = TMath::Sqrt(variance);
777 else
778 varareaerr[i] = 999999999.;
779
780 //
781 // Make also a Gauss-fit to the distributions. The RMS can be determined by
782 // outlier, thus we look at the sigma and the RMS and take the smaller one, afterwards.
783 //
784 h = camcharge.ProjectionS(TArrayI(),TArrayI(1,&aidx),"_py",750);
785 h->SetDirectory(NULL);
786 h->Fit("gaus","QL");
787 TF1 *fit = h->GetFunction("gaus");
788
789 Float_t ci2 = fit->GetChisquare();
790 Float_t sigma = fit->GetParameter(2);
791
792 if (ci2 > 500. || sigma > varareaerr[i])
793 {
794 h->Fit("gaus","QLM");
795 fit = h->GetFunction("gaus");
796
797 ci2 = fit->GetChisquare();
798 sigma = fit->GetParameter(2);
799 }
800
801 const Float_t mean = fit->GetParameter(1);
802 const Float_t ndf = fit->GetNDF();
803
804 *fLog << inf << "Camera Nr: " << i << endl;
805 *fLog << inf << option.Data() << " area idx: " << aidx << " Results: " << endl;
806 *fLog << inf << "Mean: " << Form("%4.3f",mean)
807 << "+-" << Form("%4.3f",fit->GetParError(1))
808 << " Sigma: " << Form("%4.3f",sigma) << "+-" << Form("%4.3f",fit->GetParError(2))
809 << " Chisquare: " << Form("%4.3f",fit->GetChisquare()) << " NDF : " << ndf << endl;
810 delete h;
811 gROOT->GetListOfFunctions()->Remove(fit);
812
813 if (sigma < varareaerr[i] && ndf > 2)
814 {
815 vararea [i] = mean;
816 varareaerr[i] = sigma;
817 }
818 }
819 else
820 {
821 vararea[i] = -1.;
822 varareaerr[i] = 0.;
823 }
824
825 nr[i] = i;
826 nrerr[i] = 0.;
827 }
828
829 TGraphErrors *gr = new TGraphErrors(size,
830 nr.GetArray(),vararea.GetArray(),
831 nrerr.GetArray(),varareaerr.GetArray());
832 gr->SetTitle(Form("%s Area %3i Average",option.Data(),aidx));
833 gr->GetXaxis()->SetTitle("Camera Nr.");
834 // gr->GetYaxis()->SetTitle("<Q> [1]");
835 return gr;
836}
837
838
839// -------------------------------------------------------------------
840//
841// Returns a TGraphErrors with the mean effective number of photon
842// vs. the calibration camera number. With the string 'method', different
843// calibration methods can be called.
844//
845TGraphErrors *MCalibrationIntensityChargeCam::GetPhotVsTime( const Option_t *method )
846{
847
848 const Int_t size = GetSize();
849
850 if (size == 0)
851 return NULL;
852
853 TString option(method);
854
855 TArrayF photarr(size);
856 TArrayF photarrerr(size);
857 TArrayF nr(size);
858 TArrayF nrerr(size);
859
860 for (Int_t i=0;i<GetSize();i++)
861 {
862 //
863 // Get the calibration cam from the intensity cam
864 //
865 MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);
866
867 //
868 // Get the calibration pix from the calibration cam
869 //
870 Float_t phot = 0.;
871 Float_t photerr = 0.;
872
873 if (option.Contains("BlindPixel"))
874 {
875 phot = cam->GetNumPhotonsBlindPixelMethod();
876 photerr = cam->GetNumPhotonsBlindPixelMethodErr();
877 }
878 if (option.Contains("FFactor"))
879 {
880 phot = cam->GetNumPhotonsFFactorMethod();
881 photerr = cam->GetNumPhotonsFFactorMethodErr();
882 }
883 if (option.Contains("PINDiode"))
884 {
885 phot = cam->GetNumPhotonsPINDiodeMethod();
886 photerr = cam->GetNumPhotonsPINDiodeMethodErr();
887 }
888
889 photarr[i] = phot;
890 photarrerr[i] = photerr;
891
892 nr[i] = i;
893 nrerr[i] = 0.;
894 }
895
896 TGraphErrors *gr = new TGraphErrors(size,
897 nr.GetArray(),photarr.GetArray(),
898 nrerr.GetArray(),photarrerr.GetArray());
899 gr->SetTitle("Photons Average");
900 gr->GetXaxis()->SetTitle("Camera Nr.");
901 gr->GetYaxis()->SetTitle("<N_phot> [1]");
902 return gr;
903}
904
905// -------------------------------------------------------------------
906//
907// Returns a TGraphErrors with the mean effective number of photo-electrons per
908// area index 'aidx' vs. the calibration camera number
909//
910TGraphErrors *MCalibrationIntensityChargeCam::GetPhePerAreaVsTime( const Int_t aidx, const MGeomCam &geom)
911{
912
913 const Int_t size = GetSize();
914
915 if (size == 0)
916 return NULL;
917
918 TArrayF phearea(size);
919 TArrayF pheareaerr(size);
920 TArrayF time(size);
921 TArrayF timeerr(size);
922
923 for (Int_t i=0;i<GetSize();i++)
924 {
925 //
926 // Get the calibration cam from the intensity cam
927 //
928 MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);
929
930 //
931 // Get the calibration pix from the calibration cam
932 //
933 const MCalibrationChargePix &apix = (MCalibrationChargePix&)cam->GetAverageArea(aidx);
934 const Float_t phe = apix.GetPheFFactorMethod();
935 const Float_t pheerr = apix.GetPheFFactorMethodErr();
936
937 phearea[i] = phe;
938 pheareaerr[i] = pheerr;
939
940 time[i] = i;
941 timeerr[i] = 0.;
942 }
943
944 TGraphErrors *gr = new TGraphErrors(size,
945 time.GetArray(),phearea.GetArray(),
946 timeerr.GetArray(),pheareaerr.GetArray());
947 gr->SetTitle(Form("Phes Area %d Average",aidx));
948 gr->GetXaxis()->SetTitle("Camera Nr.");
949 gr->GetYaxis()->SetTitle("<N_phes> [1]");
950 return gr;
951}
952
953// -------------------------------------------------------------------
954//
955// Returns a TGraphErrors with the event-by-event averaged charge per
956// area index 'aidx' vs. the calibration camera number
957//
958TGraphErrors *MCalibrationIntensityChargeCam::GetChargePerAreaVsTime( const Int_t aidx, const MGeomCam &geom)
959{
960
961 const Int_t size = GetSize();
962
963 if (size == 0)
964 return NULL;
965
966 TArrayF chargearea(size);
967 TArrayF chargeareaerr(size);
968 TArrayF nr(size);
969 TArrayF nrerr(size);
970
971 for (Int_t i=0;i<GetSize();i++)
972 {
973 //
974 // Get the calibration cam from the intensity cam
975 //
976 MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);
977
978 //
979 // Get the calibration pix from the calibration cam
980 //
981 const MCalibrationChargePix &apix = (MCalibrationChargePix&)cam->GetAverageArea(aidx);
982 const Float_t charge = apix.GetConvertedMean();
983 const Float_t chargeerr = apix.GetConvertedSigma();
984
985 chargearea[i] = charge;
986 chargeareaerr[i] = chargeerr;
987
988 nr[i] = i;
989 nrerr[i] = 0.;
990 }
991
992 TGraphErrors *gr = new TGraphErrors(size,
993 nr.GetArray(),chargearea.GetArray(),
994 nrerr.GetArray(),chargeareaerr.GetArray());
995 gr->SetTitle(Form("Averaged Charges Area Idx %d",aidx));
996 gr->GetXaxis()->SetTitle("Camera Nr.");
997 gr->GetYaxis()->SetTitle("<Q> [FADC cnts]");
998 return gr;
999}
1000
1001TH1F *MCalibrationIntensityChargeCam::GetVarFluctuations( const Int_t aidx, const MGeomCam &geom, const Option_t *varname )
1002{
1003
1004 const Int_t size = GetSize();
1005
1006 if (size == 0)
1007 return NULL;
1008
1009 TString option(varname);
1010
1011 TH1F *hist = new TH1F("hist",Form("%s - Rel. Fluctuations %s Pixel",option.Data(),aidx ? "Outer" : "Inner"),
1012 200,0.,100.);
1013 hist->SetXTitle("Relative Fluctuation [%]");
1014 hist->SetYTitle("Nr. channels [1]");
1015 hist->SetFillColor(kRed+aidx);
1016
1017 MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam();
1018
1019 //
1020 // Loop over pixels
1021 //
1022 for (Int_t npix=0;npix<cam->GetSize();npix++)
1023 {
1024 if (geom[npix].GetAidx() != aidx)
1025 continue;
1026
1027 Double_t variab = 0.;
1028 Double_t variab2 = 0.;
1029 Double_t variance = 0.;
1030 Int_t num = 0;
1031 Float_t pvar = 0.;
1032 Float_t relrms = 99.9;
1033 //
1034 // Loop over the Cams for each pixel
1035 //
1036 for (Int_t i=0; i<GetSize(); i++)
1037 {
1038 MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);
1039 //
1040 // Get the calibration pix from the calibration cam
1041 //
1042 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*cam)[npix];
1043 //
1044 // Don't use bad pixels
1045 //
1046 if (!pix.IsFFactorMethodValid())
1047 continue;
1048
1049 if (option.Contains("RSigma"))
1050 pvar = pix.GetRSigma();
1051 if (option.Contains("AbsTime"))
1052 pvar = pix.GetAbsTimeMean();
1053 if (option.Contains("ConversionHiLo"))
1054 pvar = pix.GetConversionHiLo();
1055 if (option.Contains("ConvertedMean"))
1056 pvar = pix.GetConvertedMean();
1057 if (option.Contains("ConvertedSigma"))
1058 pvar = pix.GetConvertedSigma();
1059 if (option.Contains("ConvertedRSigma"))
1060 pvar = pix.GetConvertedRSigma();
1061 if (option.Contains("MeanConvFADC2Phe"))
1062 pvar = pix.GetMeanConvFADC2Phe();
1063 if (option.Contains("MeanFFactorFADC2Phot"))
1064 pvar = pix.GetMeanFFactorFADC2Phot();
1065 if (option.Contains("Ped"))
1066 pvar = pix.GetPed();
1067 if (option.Contains("PedRms"))
1068 pvar = pix.GetPedRms();
1069 if (option.Contains("PheFFactorMethod"))
1070 pvar = pix.GetPheFFactorMethod();
1071 if (option.Contains("RSigmaPerCharge"))
1072 pvar = pix.GetRSigmaPerCharge();
1073
1074 variab += pvar;
1075 variab2 += pvar*pvar;
1076 num++;
1077 }
1078
1079 if (num > 1)
1080 {
1081 variab /= num;
1082 variance = (variab2 - variab*variab*num) / (num-1);
1083
1084 if (variance > 0.)
1085 relrms = TMath::Sqrt(variance)/variab * 100.;
1086 }
1087 hist->Fill(relrms);
1088 }
1089 return hist;
1090}
1091
1092void MCalibrationIntensityChargeCam::DrawRazmikPlot( const UInt_t pixid )
1093{
1094 TGraphErrors *gr = GetRazmikPlot(pixid );
1095 gr->SetBit(kCanDelete);
1096 gr->Draw("A*");
1097
1098}
1099void MCalibrationIntensityChargeCam::DrawPheVsCharge( const UInt_t pixid, const MCalibrationCam::PulserColor_t col)
1100{
1101 TGraphErrors *gr = GetPheVsCharge(pixid,col);
1102 gr->SetBit(kCanDelete);
1103 gr->Draw("A*");
1104}
1105void MCalibrationIntensityChargeCam::DrawPhePerCharge( const UInt_t pixid, const MCalibrationCam::PulserColor_t col)
1106{
1107 TGraphErrors *gr = GetPhePerCharge(pixid,MGeomCamMagic(),col);
1108 gr->SetBit(kCanDelete);
1109 gr->Draw("A*");
1110}
1111void MCalibrationIntensityChargeCam::DrawPhePerChargePerArea( const Int_t aidx, const MCalibrationCam::PulserColor_t col)
1112{
1113 TGraphErrors *gr = GetPhePerChargePerArea(aidx,MGeomCamMagic(),col);
1114 gr->SetBit(kCanDelete);
1115 gr->Draw("A*");
1116}
1117void MCalibrationIntensityChargeCam::DrawPheVsChargePerArea( const Int_t aidx, const MCalibrationCam::PulserColor_t col)
1118{
1119 TGraphErrors *gr = GetPheVsChargePerArea(aidx,col);
1120 gr->SetBit(kCanDelete);
1121 gr->Draw("A*");
1122}
1123void MCalibrationIntensityChargeCam::DrawRazmikPlotResults( const Int_t aidx)
1124{
1125 TH2F *h = GetRazmikPlotResults(aidx,MGeomCamMagic());
1126 h->SetBit(kCanDelete);
1127 h->Draw();
1128}
1129
1130void MCalibrationIntensityChargeCam::DrawChargePerAreaVsTime( const Int_t aidx)
1131{
1132 TGraphErrors *gr = GetChargePerAreaVsTime(aidx,MGeomCamMagic());
1133 gr->SetBit(kCanDelete);
1134 gr->Draw("A*");
1135}
1136void MCalibrationIntensityChargeCam::DrawPhePerAreaVsTime( const Int_t aidx)
1137{
1138 TGraphErrors *gr = GetPhePerAreaVsTime(aidx,MGeomCamMagic());
1139 gr->SetBit(kCanDelete);
1140 gr->Draw("A*");
1141}
1142void MCalibrationIntensityChargeCam::DrawPhotVsTime( const Option_t *method)
1143{
1144 TGraphErrors *gr = GetPhotVsTime(method);
1145 gr->SetBit(kCanDelete);
1146 gr->Draw("A*");
1147}
1148
1149void MCalibrationIntensityChargeCam::DrawVarPerAreaVsTime( const Int_t aidx, const Option_t *varname )
1150{
1151 TGraphErrors *gr = GetVarPerAreaVsTime(aidx,MGeomCamMagic(),varname );
1152 gr->SetBit(kCanDelete);
1153 gr->Draw("A*");
1154}
1155void MCalibrationIntensityChargeCam::DrawVarVsTime( const Int_t pixid , const Option_t *varname )
1156{
1157 TGraphErrors *gr = GetVarVsTime(pixid,varname );
1158 gr->SetBit(kCanDelete);
1159 gr->Draw("A*");
1160}
1161void MCalibrationIntensityChargeCam::DrawVarFluctuations( const Int_t aidx, const Option_t *varname)
1162{
1163 TH1F *h = GetVarFluctuations( aidx,MGeomCamMagic(),varname);
1164 h->SetBit(kCanDelete);
1165 h->Draw();
1166}
Note: See TracBrowser for help on using the repository browser.