Context Navigation

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

Visit:

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: