Context Navigation

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

Visit:

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

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

Download in other formats: