Context Navigation

MCalibrationChargeCam.cc@ 19365

Visit:

Last change on this file since 19365 was 10166, checked in by tbretz, 14 years ago
Removed the old obsolete cvs header line.
File size: 37.1 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-2008
21	!
22	!
23	\* ======================================================================== */
24	/////////////////////////////////////////////////////////////////////////////
25	//
26	// MCalibrationChargeCam
27	//
28	// Storage container for charge calibration results from the signal distribution
29	// fits (see MHCalibrationChargeCam and MHCalibrationChargePix), the calculation
30	// of reduced sigmas and number of photo-electrons (this class) and conversion
31	// factors sum FADC slices to photo-electrons (see MCalibrationChargeCalc)
32	//
33	// Individual pixels have to be cast when retrieved e.g.:
34	// MCalibrationChargePix &avpix = (MCalibrationChargePix&)(*fChargeCam)[i]
35	//
36	// Averaged values over one whole area index (e.g. inner or outer pixels for
37	// the MAGIC camera), can be retrieved via:
38	// MCalibrationChargePix &avpix = (MCalibrationChargePix&)fChargeCam->GetAverageArea(i)
39	//
40	// Averaged values over one whole camera sector can be retrieved via:
41	// MCalibrationChargePix &avpix = (MCalibrationChargePix&)fChargeCam->GetAverageSector(i)
42	//
43	// Note the averageing has been done on an event-by-event basis. Resulting
44	// Sigma's of the Gauss fit have been multiplied with the square root of the number
45	// of involved pixels in order to make a direct comparison possible with the mean of
46	// sigmas.
47	//
48	// Final numbers of uncalibrated or unreliable pixels can be retrieved via the commands:
49	// GetNumUncalibrated(aidx) and GetNumUnreliable(aidx) where aidx is the area index (0 for
50	// inner and 1 for outer pixels in the MAGIC camera).
51	//
52	// The following "calibration" constants are used for the calibration of each pixel
53	// (see MCalibrate):
54	//
55	// - MCalibrationQEPix::GetMeanConvFADC2Phe(): The mean conversion factor from the
56	// summed FADC slices to the number of photo-electrons (in first order independent
57	// on colour and intensity)
58	// - MCalibrationQEPix::GetMeanFFactorFADC2Phot(): The mean F-Factor of the total
59	// readout chain dividing the signal to noise of the incoming number of photons
60	// (= sqrt(number photons)) by the signal to noise of the outgoing summed FADC slices
61	// signal (= MCalibrationPix::GetMean() / MCalibrationChargePix::GetRSigma() )
62	//
63	// The following calibration constants can be retrieved directly from this class:
64	//
65	// - GetConversionFactorFFactor ( Int_t idx, Float_t &mean, Float_t &err, Float_t &sigma );
66	//
67	// where:
68	// - idx is the pixel software ID
69	// - "mean" is the mean conversion constant, to be multiplied with the retrieved signal
70	// in order to get a calibrated number of PHOTONS.
71	// - "err" is the pure statistical uncertainty about the MEAN
72	// - "sigma", if mulitplied with the square root of signal, gives the approximate sigma of the
73	// retrieved mean number of incident Cherenkov photons.
74	//
75	// Note, Conversion is ALWAYS (included the F-Factor method!) from summed FADC slices to PHOTONS.
76	//
77	// See also: MCalibrationChargePix, MCalibrationChargeCalc, MCalibrationQECam
78	// MHCalibrationChargePix, MHCalibrationChargeCam
79	// MCalibrationBlindPix, MCalibrationChargePINDiode
80	//
81	/////////////////////////////////////////////////////////////////////////////
82	#include "MCalibrationChargeCam.h"
83	#include "MCalibrationChargePix.h"
84
85	#include <TMath.h>
86
87	#include <TH1.h>
88	#include <TF1.h>
89	#include <TOrdCollection.h>
90
91	#include "MLog.h"
92	#include "MLogManip.h"
93
94	#include "MGeomCam.h"
95	#include "MGeom.h"
96
97	#include "MBadPixelsCam.h"
98	#include "MBadPixelsPix.h"
99
100	#include "MCalibrationQECam.h"
101	#include "MCalibrationQEPix.h"
102
103	#include "MHCamera.h"
104
105	ClassImp(MCalibrationChargeCam);
106
107	using namespace std;
108	// --------------------------------------------------------------------------
109	//
110	// Default constructor.
111	//
112	// Calls:
113	// - Clear()
114	//
115	MCalibrationChargeCam::MCalibrationChargeCam(const char name, const char title)
116	{
117	fName = name ? name : "MCalibrationChargeCam";
118	fTitle = title ? title : "Storage container for the Calibration Information in the camera";
119
120	Clear();
121	}
122
123	/*
124	// --------------------------------------------------------------------------
125	//
126	// Creates new MCalibrationCam only for the averaged areas:
127	// the rest has to be retrieved directly, e.g. via:
128	//
129	TObject MCalibrationChargeCam::Clone(const char ) const
130	{
131
132	//
133	// FIXME, this might be done faster and more elegant, by direct copy.
134	//
135	MCalibrationChargeCam *cam = new MCalibrationChargeCam(fName,fTitle);
136
137	cam->fNumUnsuitable = fNumUnsuitable;
138	cam->fNumUnreliable = fNumUnreliable;
139	cam->fNumHiGainFADCSlices = fNumHiGainFADCSlices;
140	cam->fNumLoGainFADCSlices = fNumLoGainFADCSlices;
141	cam->fPulserColor = fPulserColor;
142
143	cam->fFlags = fFlags;
144
145	cam->fNumPhotonsBlindPixelMethod = fNumPhotonsBlindPixelMethod;
146	cam->fNumPhotonsFFactorMethod = fNumPhotonsFFactorMethod;
147	cam->fNumPhotonsPINDiodeMethod = fNumPhotonsPINDiodeMethod;
148	cam->fNumPhotonsBlindPixelMethodErr = fNumPhotonsBlindPixelMethodErr;
149	cam->fNumPhotonsFFactorMethodErr = fNumPhotonsFFactorMethodErr;
150	cam->fNumPhotonsPINDiodeMethodErr = fNumPhotonsPINDiodeMethodErr;
151
152	for (Int_t i=0; i<GetSize(); i++)
153	cam->fPixels->AddAt((*this)[i].Clone(),i);
154
155	for (Int_t i=0; i<GetAverageAreas(); i++)
156	{
157	cam->fAverageAreas->AddAt(GetAverageArea(i).Clone(),i);
158	cam->fAverageBadAreas->AddAt(GetAverageBadArea(i).Clone(),i);
159	}
160	for (Int_t i=0; i<GetAverageSectors(); i++)
161	{
162	cam->fAverageSectors->AddAt(GetAverageSector(i).Clone(),i);
163	cam->fAverageBadSectors->AddAt(GetAverageBadSector(i).Clone(),i);
164	}
165
166	return cam;
167	}
168	*/
169
170	// -------------------------------------------------------------------
171	//
172	// Add MCalibrationChargePix's in the ranges from - to to fPixels
173	//
174	void MCalibrationChargeCam::Add(const UInt_t a, const UInt_t b)
175	{
176	for (UInt_t i=a; i<b; i++)
177	{
178	fPixels->AddAt(new MCalibrationChargePix,i);
179	(*this)[i].SetPixId(i);
180	}
181	}
182
183	// -------------------------------------------------------------------
184	//
185	// Add MCalibrationChargePix's in the ranges from - to to fAverageAreas
186	//
187	void MCalibrationChargeCam::AddArea(const UInt_t a, const UInt_t b)
188	{
189	for (UInt_t i=a; i<b; i++)
190	{
191	fAverageAreas->AddAt(new MCalibrationChargePix,i);
192	GetAverageArea(i).SetPixId(i);
193	}
194	}
195
196	// -------------------------------------------------------------------
197	//
198	// Add MCalibrationChargePix's in the ranges from - to to fAverageSectors
199	//
200	void MCalibrationChargeCam::AddSector(const UInt_t a, const UInt_t b)
201	{
202	for (UInt_t i=a; i<b; i++)
203	{
204	fAverageSectors->AddAt(new MCalibrationChargePix,i);
205	GetAverageSector(i).SetPixId(i);
206	}
207	}
208
209
210	// --------------------------------------
211	//
212	// Sets all variable to 0.
213	// Sets all flags to kFALSE
214	// Calls MCalibrationCam::Clear()
215	//
216	void MCalibrationChargeCam::Clear(Option_t *o)
217	{
218	// DO NOT RESET THE HI-/LO-GAIN CONSTANTS
219	SetFFactorMethodValid ( kFALSE );
220
221	fNumPhotonsBlindPixelMethod = 0.;
222	fNumPhotonsFFactorMethod = 0.;
223	fNumPhotonsPINDiodeMethod = 0.;
224	fNumPhotonsBlindPixelMethodErr = 0.;
225	fNumPhotonsFFactorMethodErr = 0.;
226	fNumPhotonsPINDiodeMethodErr = 0.;
227
228	MCalibrationCam::Clear();
229	}
230
231	// -----------------------------------------------
232	//
233	// Sets the kFFactorMethodValid bit from outside
234	//
235	void MCalibrationChargeCam::SetFFactorMethodValid(const Bool_t b)
236	{
237	b ? SETBIT(fFlags, kFFactorMethodValid) : CLRBIT(fFlags, kFFactorMethodValid);
238	}
239
240
241	// --------------------------------------------------------------------------
242	//
243	// Test bit kFFactorMethodValid
244	//
245	Bool_t MCalibrationChargeCam::IsFFactorMethodValid() const
246	{
247	return TESTBIT(fFlags,kFFactorMethodValid);
248	}
249
250	// --------------------------------------------------------------------------
251	//
252	// Copy High-gain vs. low-gain conversion factors from cam to this.
253	//
254	Bool_t MCalibrationChargeCam::MergeHiLoConversionFactors(const MCalibrationChargeCam &cam) const
255	{
256
257	if (GetSize() != cam.GetSize())
258	{
259	*fLog << warn << "Sizes mismatch! Cannot merge high-gain vs. low-gain convertion factors" << endl;
260	return kFALSE;
261	}
262
263	for (Int_t i=0; i<GetSize(); i++)
264	{
265	((MCalibrationChargePix&)(*this)[i]).SetConversionHiLo (((MCalibrationChargePix&)cam[i]).GetConversionHiLo());
266	((MCalibrationChargePix&)(*this)[i]).SetConversionHiLoErr(((MCalibrationChargePix&)cam[i]).GetConversionHiLoErr());
267	}
268
269	return kTRUE;
270	}
271
272	// --------------------------------------------------------------------------
273	//
274	// Print first the well fitted pixels
275	// and then the ones which are not FitValid
276	//
277	void MCalibrationChargeCam::Print(Option_t *o) const
278	{
279
280	*fLog << all << GetDescriptor() << ":" << endl;
281	int id = 0;
282
283	*fLog << all << "Calibrated pixels:" << endl;
284	*fLog << all << endl;
285
286	TIter Next(fPixels);
287	MCalibrationChargePix *pix;
288	while ((pix=(MCalibrationChargePix*)Next()))
289	{
290
291	if (pix->IsFFactorMethodValid())
292	{
293
294	*fLog << all
295	<< Form("%s%3i","Pixel: ",pix->GetPixId())
296	<< Form("%s%4.2f%s%4.2f"," Ped.Rms: ",pix->GetPedRms(),"+-",pix->GetPedRmsErr())
297	<< Form("%s%4.2f%s%4.2f"," Charge: " ,pix->GetConvertedMean(),"+-",pix->GetConvertedSigma())
298	<< Form("%s%4.2f%s%4.2f"," Red.Sigma: ",pix->GetConvertedRSigma(),"+-",pix->GetConvertedRSigmaErr())
299	<< Form("%s%4.2f%s%4.2f"," Num.Phes: ",pix->GetPheFFactorMethod(),"+-",pix->GetPheFFactorMethodErr())
300	<< Form("%s%4.2f%s%4.2f"," Conv.FADC2Phe: ",pix->GetMeanConvFADC2Phe(),"+-",pix->GetMeanConvFADC2PheErr())
301	<< " Saturated? :" << pix->IsHiGainSaturation()
302	<< Form("%s%4.2f%s%4.2f"," Conv.HiLo: ",pix->GetConversionHiLo(),"+-",pix->GetConversionHiLoErr())
303	<< endl;
304	id++;
305	}
306	}
307
308	*fLog << all << id << " pixels" << endl;
309	id = 0;
310
311
312	*fLog << all << endl;
313	*fLog << all << "Previously Excluded pixels:" << endl;
314	*fLog << all << endl;
315
316	id = 0;
317
318	TIter Next4(fPixels);
319	while ((pix=(MCalibrationChargePix*)Next4()))
320	{
321	if (pix->IsExcluded())
322	{
323	*fLog << all << pix->GetPixId() << " ";
324	id++;
325
326	if (!(id % 25))
327	*fLog << endl;
328	}
329	}
330
331	*fLog << all << endl;
332	*fLog << all << "New Excluded pixels:" << endl;
333	*fLog << all << endl;
334
335	TIter Next5(fPixels);
336	while ((pix=(MCalibrationChargePix*)Next5()))
337	{
338	if (!pix->IsFFactorMethodValid() && !pix->IsExcluded())
339	{
340	*fLog << all << pix->GetPixId() << " ";
341	id++;
342
343	if (!(id % 25))
344	*fLog << endl;
345	}
346	}
347
348	*fLog << endl;
349	*fLog << all << id << " Excluded pixels " << endl;
350	*fLog << endl;
351
352	*fLog << all << endl;
353	*fLog << all << "Averaged Areas:" << endl;
354	*fLog << all << endl;
355
356	TIter Next6(fAverageAreas);
357	while ((pix=(MCalibrationChargePix*)Next6()))
358	{
359	*fLog << all << Form("%s%3i","Area Idx: ",pix->GetPixId())
360	<< " Ped. Rms: " << pix->GetPedRms() << " +- " << pix->GetPedRmsErr()
361	<< " Mean signal: " << pix->GetMean() << " +- " << pix->GetMeanErr()
362	<< " Sigma signal: " << pix->GetSigma() << " +- "<< pix->GetSigmaErr()
363	<< " Reduced Sigma: " << pix->GetRSigma()
364	<< " Nr Phe's: " << pix->GetPheFFactorMethod()
365	<< endl;
366	}
367
368	*fLog << all << endl;
369	*fLog << all << "Averaged Sectors:" << endl;
370	*fLog << all << endl;
371
372	TIter Next7(fAverageSectors);
373	while ((pix=(MCalibrationChargePix*)Next7()))
374	{
375	*fLog << all << Form("%s%3i","Sector: ",pix->GetPixId())
376	<< " Ped. Rms: " << pix->GetPedRms() << " +- " << pix->GetPedRmsErr()
377	<< " Mean signal: " << pix->GetMean() << " +- " << pix->GetMeanErr()
378	<< " Sigma signal: " << pix->GetSigma() << " +- "<< pix->GetSigmaErr()
379	<< " Reduced Sigma: " << pix->GetRSigma()
380	<< " Nr Phe's: " << pix->GetPheFFactorMethod()
381	<< endl;
382	}
383	*fLog << all << endl;
384	}
385
386
387	// --------------------------------------------------------------------------
388	//
389	// The types are as follows:
390	//
391	// Fitted values:
392	// ==============
393	//
394	// 0: Fitted Charge (see MCalibrationPix::GetMean())
395	// 1: Error of fitted Charge (see MCalibrationPix::GetMeanErr())
396	// 2: Sigma of fitted Charge (see MCalibrationPix::GetSigma())
397	// 3: Error of Sigma of fitted Charge (see MCalibrationPix::GetSigmaErr())
398	//
399	// Useful variables derived from the fit results:
400	// =============================================
401	//
402	// 4: Probability Gauss fit Charge distribution (see MCalibrationPix::GetProb())
403	// 5: Reduced Sigma of fitted Charge (see MCalibrationChargePix::GetRSigma())
404	// 6: Error Reduced Sigma of fitted Charge (see MCalibrationChargePix::GetRSigmaErr())
405	// 7: Reduced Sigma per Charge (see MCalibrationChargePix::GetRSigmaPerCharge())
406	// 8: Error of Reduced Sigma per Charge (see MCalibrationChargePix::GetRSigmaPerChargeErr())
407	//
408	// Results of the F-Factor calibration Method:
409	// ===========================================
410	//
411	// 9: Nr. Photo-electrons from F-Factor Method (see MCalibrationChargePix::GetPheFFactorMethod())
412	// 10: Error Nr. Photo-el. from F-Factor Method (see MCalibrationChargePix::GetPheFFactorMethodErr())
413	// 11: Conversion factor from F-Factor Method (see MCalibrationChargePix::GetMeanConvFADC2Phe()
414	// 12: Error conv. factor from F-Factor Method (see MCalibrationChargePix::GetMeanConvFADC2PheErr()
415	// 13: Overall F-Factor from F-Factor Method (see MCalibrationChargePix::GetMeanFFactorFADC2Phot()
416	// 14: Error F-Factor from F-Factor Method (see MCalibrationChargePix::GetMeanFFactorFADC2PhotErr()
417	// 15: Pixels valid calibration F-Factor-Method (see MCalibrationChargePix::IsFFactorMethodValid())
418	//
419	// Results of the Low-Gain vs. High-Gain Conversion:
420	// =================================================
421	//
422	// 16: Mean Signal Hi-Gain / Mean Signal Lo-Gain (see MCalibrationPix::GetHiLoMeansDivided())
423	// 17: Error Signal High-Gain / Signal Low Gain (see MCalibrationPix::GetHiLoMeansDividedErr())
424	// 18: Sigma High-Gain / Sigma Low Gain (see MCalibrationPix::GetHiLoSigmasDivided())
425	// 19: Error Sigma High-Gain / Sigma Low Gain (see MCalibrationPix::GetHiLoSigmasDividedErr())
426	//
427	// Localized defects:
428	// ==================
429	//
430	// 20: Excluded Pixels
431	// 21: Number of pickup events in the Hi Gain (see MCalibrationPix::GetHiGainNumPickup())
432	// 22: Number of pickup events in the Lo Gain (see MCalibrationPix::GetLoGainNumPickup())
433	// 23: Number of blackout events in the Hi Gain (see MCalibrationPix::GetHiGainNumBlackout())
434	// 24: Number of blackout events in the Lo Gain (see MCalibrationPix::GetLoGainNumBlackout())
435	//
436	// Other classifications of pixels:
437	// ================================
438	//
439	// 25: Pixels with saturated High-Gain (see MCalibrationPix::IsHiGainSaturation())
440	//
441	// Calculated absolute arrival times (very low precision!):
442	// ========================================================
443	//
444	// 26: Absolute Arrival time of the signal (see MCalibrationChargePix::GetAbsTimeMean())
445	// 27: RMS Ab. Arrival time of the signal (see MCalibrationChargePix::GetAbsTimeRms())
446	//
447	// Used Pedestals:
448	// ===============
449	//
450	// 28: Pedestal for entire signal extr. range (see MCalibrationChargePix::Ped())
451	// 29: Error Pedestal entire signal extr. range (see MCalibrationChargePix::PedErr())
452	// 30: Ped. RMS entire signal extraction range (see MCalibrationChargePix::PedRms())
453	// 31: Error Ped. RMS entire signal extr. range (see MCalibrationChargePix::PedRmsErr())
454	//
455	// Special variables (for data check):
456	// ====================================
457	//
458	// 32: HiGain RMS divided by Mean for every pixel (with inclusion of the excluded pixels)
459	//
460	Bool_t MCalibrationChargeCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const
461	{
462	if (idx > GetSize())
463	return kFALSE;
464
465	const Float_t area = cam[idx].GetA();
466	if (area == 0)
467	return kFALSE;
468
469	const MCalibrationChargePix &pix = (MCalibrationChargePix&)(*this)[idx];
470	switch (type)
471	{
472	case 0:
473	if (pix.IsExcluded() \|\| pix.GetConvertedMeanErr()<0)
474	return kFALSE;
475	val = pix.GetConvertedMean();
476	return kTRUE;
477	case 1:
478	if (pix.IsExcluded() \|\| pix.GetConvertedMeanErr()<0)
479	return kFALSE;
480	val = pix.GetConvertedMeanErr();
481	return kTRUE;
482	case 2:
483	if (pix.IsExcluded() \|\| pix.GetConvertedSigmaErr()<0)
484	return kFALSE;
485	val = pix.GetConvertedSigma();
486	return kTRUE;
487	case 3:
488	if (pix.IsExcluded() \|\| pix.GetConvertedSigmaErr()<0)
489	return kFALSE;
490	val = pix.GetConvertedSigmaErr();
491	return kTRUE;
492	case 4:
493	if (pix.IsExcluded())
494	return kFALSE;
495	val = pix.GetProb();
496	return val>=0;
497	case 5:
498	if (!pix.IsFFactorMethodValid())
499	return kFALSE;
500	if (pix.GetRSigma() == -1.)
501	return kFALSE;
502	val = pix.GetConvertedRSigma();
503	break;
504	case 6:
505	if (!pix.IsFFactorMethodValid())
506	return kFALSE;
507	if (pix.GetRSigma() == -1.)
508	return kFALSE;
509	val = pix.GetConvertedRSigmaErr();
510	break;
511	case 7:
512	if (!pix.IsFFactorMethodValid())
513	return kFALSE;
514	val = pix.GetRSigmaPerCharge();
515	break;
516	case 8:
517	if (!pix.IsFFactorMethodValid())
518	return kFALSE;
519	val = pix.GetRSigmaPerChargeErr();
520	break;
521	case 9:
522	// if (!pix.IsFFactorMethodValid())
523	// return kFALSE;
524	if (pix.IsExcluded() \|\| pix.GetPheFFactorMethodErr()<0)
525	return kFALSE;
526	val = pix.GetPheFFactorMethod();
527	return kTRUE;
528	case 10:
529	// if (!pix.IsFFactorMethodValid())
530	// return kFALSE;
531	if (pix.IsExcluded() \|\| pix.GetPheFFactorMethodErr()<=0)
532	return kFALSE;
533	val = pix.GetPheFFactorMethodErr();
534	return kTRUE;
535	case 11:
536	// if (!pix.IsFFactorMethodValid())
537	// return kFALSE;
538	if (pix.IsExcluded() \|\| pix.GetMeanConvFADC2PheErr()<0)
539	return kFALSE;
540	val = pix.GetMeanConvFADC2Phe();
541	return kTRUE;
542	case 12:
543	// if (!pix.IsFFactorMethodValid())
544	// return kFALSE;
545	if (pix.IsExcluded() \|\| pix.GetMeanConvFADC2PheErr()<=0)
546	return kFALSE;
547	val = pix.GetMeanConvFADC2PheErr();
548	return kTRUE;
549	case 13:
550	// if (!pix.IsFFactorMethodValid())
551	// return kFALSE;
552	if (pix.IsExcluded())
553	return kFALSE;
554	val = pix.GetMeanFFactorFADC2Phot();
555	break;
556	case 14:
557	// if (!pix.IsFFactorMethodValid())
558	// return kFALSE;
559	if (pix.IsExcluded())
560	return kFALSE;
561	val = pix.GetMeanFFactorFADC2PhotErr();
562	if (val <= 0.)
563	val = 0.00001;
564	break;
565	case 15:
566	if (pix.IsExcluded())
567	return kFALSE;
568	if (!pix.IsFFactorMethodValid())
569	return kFALSE;
570	val = 1;
571	return kTRUE;
572	case 16:
573	if (pix.IsExcluded())
574	return kFALSE;
575	val = pix.GetHiLoMeansDivided();
576	break;
577	case 17:
578	if (pix.IsExcluded())
579	return kFALSE;
580	val = pix.GetHiLoMeansDividedErr();
581	break;
582	case 18:
583	if (pix.IsExcluded())
584	return kFALSE;
585	val = pix.GetHiLoSigmasDivided();
586	break;
587	case 19:
588	if (pix.IsExcluded())
589	return kFALSE;
590	val = pix.GetHiLoSigmasDividedErr();
591	break;
592	case 20:
593	if (!pix.IsExcluded())
594	return kFALSE;
595	val = 1.;
596	return kTRUE;
597	case 21:
598	if (pix.IsExcluded())
599	return kFALSE;
600	val = pix.GetHiGainNumPickup();
601	break;
602	case 22:
603	if (pix.IsExcluded())
604	return kFALSE;
605	val = pix.GetLoGainNumPickup();
606	break;
607	case 23:
608	if (pix.IsExcluded())
609	return kFALSE;
610	val = pix.GetHiGainNumBlackout();
611	break;
612	case 24:
613	if (pix.IsExcluded())
614	return kFALSE;
615	val = pix.GetLoGainNumBlackout();
616	break;
617	case 25:
618	if (pix.IsExcluded())
619	return kFALSE;
620	val = pix.IsHiGainSaturation();
621	break;
622	case 26:
623	if (pix.IsExcluded())
624	return kFALSE;
625	val = pix.GetAbsTimeMean();
626	break;
627	case 27:
628	if (pix.IsExcluded())
629	return kFALSE;
630	val = pix.GetAbsTimeRms();
631	break;
632	case 28:
633	if (pix.IsExcluded())
634	return kFALSE;
635	val = pix.GetPed();
636	break;
637	case 29:
638	if (pix.IsExcluded())
639	return kFALSE;
640	val = pix.GetPedErr();
641	break;
642	case 30:
643	if (pix.IsExcluded())
644	return kFALSE;
645	val = pix.GetPedRms();
646	break;
647	case 31:
648	if (pix.IsExcluded())
649	return kFALSE;
650	val = pix.GetPedErr()/2.;
651	break;
652	case 32:
653	if (pix.IsExcluded() \|\| pix.GetRms()<0)
654	return kFALSE;
655	val = pix.GetMean() == 0. ? 0. : pix.GetRms()/pix.GetMean();
656	return kTRUE;
657	case 33:
658	if (pix.IsExcluded() \|\| pix.GetRms()<0)
659	return kFALSE;
660	if (pix.GetMean() == 0.)
661	val = 0.;
662	else
663	val = pix.GetSigmaErr()/pix.GetMean() + pix.GetRms()*pix.GetMeanErr()/pix.GetMean()/pix.GetMean();
664	return kTRUE;
665	default:
666	return kFALSE;
667	}
668
669	return val!=-1.;
670	}
671
672	Bool_t MCalibrationChargeCam::GetConversionFactorFFactor(Int_t ipx, Float_t &mean, Float_t &ferr, Float_t &ffactor)
673	{
674
675	MCalibrationChargePix &pix = (MCalibrationChargePix&)(*this)[ipx];
676
677	Float_t conv = pix.GetMeanConvFADC2Phe();
678
679	if (conv < 0.)
680	return kFALSE;
681
682	mean = conv;
683	ferr = pix.GetMeanConvFADC2PheErr();
684	ffactor = pix.GetMeanFFactorFADC2Phot();
685
686	return kTRUE;
687	}
688
689
690	// --------------------------------------------------------------------------
691	//
692	// Calculates the average conversion factor FADC counts to photons for pixel sizes.
693	// The geometry container is used to get the necessary
694	// geometry information (area index). The MCalibrationQECam container is necessary for
695	// the quantum efficiency information.
696	// If the bad pixel container is given all pixels which have the flag 'kUnsuitableRun' are ignored
697	// in the calculation of the size average.
698	//
699	// Returns a TArrayF of dimension 2:
700	// arr[0]: averaged conversion factors (default: -1.)
701	// arr[1]: Error (rms) of averaged conversion factors (default: 0.)
702	//
703	TArrayF MCalibrationChargeCam::GetAveragedConvFADC2PhotPerArea ( const MGeomCam &geom, const MCalibrationQECam &qecam,
704	const UInt_t ai, MBadPixelsCam *bad)
705	{
706
707	const Int_t np = GetSize();
708
709	Double_t mean = 0.;
710	Double_t mean2 = 0.;
711	Int_t nr = 0;
712
713	MHCamera convcam(geom,"ConvFactors","Conversion Factors;Conv Factor [phot/FADC cnts];channels");
714
715	for (int i=0; i<np; i++)
716	{
717	if (bad && (*bad)[i].IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
718	continue;
719
720	if (bad && (*bad)[i].IsUncalibrated(MBadPixelsPix::kDeviatingNumPhes))
721	continue;
722
723	const UInt_t aidx = geom[i].GetAidx();
724
725	if (ai != aidx)
726	continue;
727
728	const MCalibrationQEPix &qepix = (MCalibrationQEPix&)qecam[i];
729	if (!qepix.IsAverageQEFFactorAvailable())
730	continue;
731
732	const MCalibrationChargePix &pix = (MCalibrationChargePix&)(*this)[i];
733	const Float_t conv = pix.GetMeanConvFADC2Phe()/qepix.GetQECascadesFFactor();
734
735	mean += conv;
736	mean2 += conv*conv;
737	nr ++;
738
739	convcam.Fill(i,conv);
740	convcam.SetUsed(i);
741	}
742
743	Float_t mn = nr ? mean/nr : -1.;
744	Float_t sg = nr>1 ? TMath::Sqrt((mean2 - mean*mean/nr)/(nr-1)) : 0.;
745
746	const Int_t aidx = (Int_t)ai;
747
748	TH1D *h = convcam.ProjectionS(TArrayI(),TArrayI(1,&aidx),"_py",750);
749	h->SetDirectory(NULL);
750
751	TF1 *fit = NULL;
752
753	if (geom.InheritsFrom("MGeomCamMagic"))
754	{
755
756	fit = new TF1("fit","gaus",0.4,5.);
757
758	// Fix the ranges, as found by Nadia
759	if(aidx == 0)
760	h->Fit(fit, "REQ0", "",0.4,1.5);
761	else
762	h->Fit(fit ,"REQ0", "",1.,5.);
763	}
764	else
765	{
766	h->Fit("gaus","Q0");
767	fit = h->GetFunction("gaus");
768	}
769
770	Float_t ci2 = fit->GetChisquare();
771	Float_t sigma = fit->GetParameter(2);
772
773	if (ci2 > 500. \|\| sigma > sg)
774	{
775	if (geom.InheritsFrom("MGeomCamMagic"))
776	{
777	// Fix the ranges, as found by Nadia
778	if(aidx == 0)
779	h->Fit(fit, "REQ0", "",0.4,1.5);
780	else
781	h->Fit(fit, "REQ0", "",1.,5.);
782	}
783	else
784	{
785	h->Fit("gaus","MREQ0");
786	fit = h->GetFunction("gaus");
787	}
788
789	ci2 = fit->GetChisquare();
790	sigma = fit->GetParameter(2);
791	}
792
793	const Int_t ndf = fit->GetNDF();
794
795	if (ci2 < 500. && sigma < sg && ndf > 2)
796	{
797	mn = fit->GetParameter(1);
798	sg = sigma;
799	}
800
801	*fLog << inf << "Conversion Factors to photons area idx: " << aidx << ":" << endl;
802	*fLog << inf << "Mean: " << Form("%4.3f",mn)
803	<< "+-" << Form("%4.3f",fit->GetParError(1))
804	<< " Sigma: " << Form("%4.3f",sg) << "+-" << Form("%4.3f",fit->GetParError(2))
805	<< " Chisquare: " << Form("%4.3f",fit->GetChisquare()) << " NDF : " << ndf << endl;
806
807	delete fit;
808	delete h;
809	gROOT->GetListOfFunctions()->Remove(fit);
810
811	TArrayF arr(2);
812	arr[0] = mn;
813	arr[1] = sg;
814
815	return arr;
816	}
817
818	// --------------------------------------------------------------------------
819	//
820	// Calculates the average conversion factor FADC counts to equiv. photo-electrons for pixel sizes.
821	// The geometry container is used to get the necessary
822	// geometry information (area index). The MCalibrationQECam container is necessary for
823	// the quantum efficiency information.
824	// If the bad pixel container is given all pixels which have the flag 'kUnsuitableRun' are ignored
825	// in the calculation of the size average.
826	//
827	// Returns a TArrayF of dimension 2:
828	// arr[0]: averaged conversion factors (default: -1.)
829	// arr[1]: Error (rms) of averaged conversion factors (default: 0.)
830	//
831	TArrayF MCalibrationChargeCam::GetAveragedConvFADC2PhePerArea ( const MGeomCam &geom, const MCalibrationQECam &qecam,
832	const UInt_t ai, MBadPixelsCam *bad)
833	{
834
835	const Int_t np = GetSize();
836
837	Double_t mean = 0.;
838	Double_t mean2 = 0.;
839	Int_t nr = 0;
840
841	MHCamera convcam(geom,"ConvFactors","Conversion Factors;Conv Factor [phe/FADC cnts];channels");
842
843	for (int i=0; i<np; i++)
844	{
845	if (bad && (*bad)[i].IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
846	continue;
847
848	if (bad && (*bad)[i].IsUncalibrated(MBadPixelsPix::kDeviatingNumPhes))
849	continue;
850
851	const UInt_t aidx = geom[i].GetAidx();
852
853	if (ai != aidx)
854	continue;
855
856	const MCalibrationQEPix &qepix = (MCalibrationQEPix&)qecam[i];
857	if (!qepix.IsAverageQEFFactorAvailable())
858	continue;
859
860	const MCalibrationChargePix &pix = (MCalibrationChargePix&)(*this)[i];
861	const Float_t conv = pix.GetMeanConvFADC2Phe()/qepix.GetQECascadesFFactor()*MCalibrationQEPix::gkDefaultAverageQE;
862
863	mean += conv;
864	mean2 += conv*conv;
865	nr ++;
866
867	convcam.Fill(i,conv);
868	convcam.SetUsed(i);
869	}
870
871	Float_t mn = nr ? mean/nr : -1.;
872	Float_t sg = nr>1 ? TMath::Sqrt((mean2 - mean*mean/nr)/(nr-1)) : 0.;
873
874	const Int_t aidx = (Int_t)ai;
875
876	TH1D *h = convcam.ProjectionS(TArrayI(),TArrayI(1,&aidx),"_py",750);
877	h->SetDirectory(NULL);
878
879	TF1 *fit = NULL;
880
881	if (geom.InheritsFrom("MGeomCamMagic"))
882	{
883
884	fit = new TF1("fit","gaus",0.01,1.);
885
886	// Fix the ranges, as found by Nadia
887	if(aidx == 0)
888	{h->Fit("fit","REQ0", "",0.07,0.3);}
889	else
890	{h->Fit("fit","REQ0", "",0.15,1.0);}
891	}
892	else
893	{
894	h->Fit("gaus","Q0");
895	fit = h->GetFunction("gaus");
896	}
897
898	Float_t ci2 = fit->GetChisquare();
899	Float_t sigma = fit->GetParameter(2);
900
901	if (ci2 > 500. \|\| sigma > sg)
902	{
903	if (geom.InheritsFrom("MGeomCamMagic"))
904	{
905	// Fix the ranges, as found by Nadia
906	if(aidx == 0)
907	{h->Fit("fit","REQ0", "",0.07,0.3);}
908	else
909	{h->Fit("fit","REQ0", "",0.15,1.0);}
910	}
911	else
912	{
913	h->Fit("gaus","MREQ0");
914	fit = h->GetFunction("gaus");
915	}
916
917	ci2 = fit->GetChisquare();
918	sigma = fit->GetParameter(2);
919	}
920
921	const Int_t ndf = fit->GetNDF();
922
923	if (ci2 < 500. && sigma < sg && ndf > 2)
924	{
925	mn = fit->GetParameter(1);
926	sg = sigma;
927	}
928
929	*fLog << inf << "Conversion Factors to equiv. photo-electrons area idx: " << aidx << ":" << endl;
930	*fLog << inf << "Mean: " << Form("%4.3f",mn)
931	<< "+-" << Form("%4.3f",fit->GetParError(1))
932	<< " Sigma: " << Form("%4.3f",sg) << "+-" << Form("%4.3f",fit->GetParError(2))
933	<< " Chisquare: " << Form("%4.3f",fit->GetChisquare()) << " NDF : " << ndf << endl;
934
935	delete fit;
936	delete h;
937	gROOT->GetListOfFunctions()->Remove(fit);
938
939	TArrayF arr(2);
940	arr[0] = mn;
941	arr[1] = sg;
942
943	return arr;
944	}
945
946	// --------------------------------------------------------------------------
947	//
948	// Calculates the average conversion factor FADC counts to photons for camera sectors.
949	// The geometry container is used to get the necessary
950	// geometry information (area index). The MCalibrationQECam container is necessary for
951	// the quantum efficiency information.
952	// If the bad pixel container is given all pixels which have the flag 'kUnsuitableRun' are ignored
953	// in the calculation of the size average.
954	//
955	// Returns a TArrayF of dimension 2:
956	// arr[0]: averaged conversion factors (default: -1.)
957	// arr[1]: Error (rms) of averaged conversion factors (default: 0.)
958	//
959	TArrayF MCalibrationChargeCam::GetAveragedConvFADC2PhotPerSector( const MGeomCam &geom, const MCalibrationQECam &qecam,
960	const UInt_t sec, MBadPixelsCam *bad)
961	{
962	const Int_t np = GetSize();
963
964	Double_t mean = 0.;
965	Double_t mean2 = 0.;
966	Int_t nr = 0;
967
968	for (int i=0; i<np; i++)
969	{
970	if (bad && (*bad)[i].IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
971	continue;
972
973	const UInt_t sector = geom[i].GetSector();
974
975	if (sector != sec)
976	continue;
977
978	const MCalibrationQEPix &qepix = (MCalibrationQEPix&)qecam[i];
979	if (!qepix.IsAverageQEFFactorAvailable())
980	continue;
981
982	const MCalibrationChargePix &pix = (MCalibrationChargePix&)(*this)[i];
983	const Float_t conv = pix.GetMeanConvFADC2Phe();
984	const Float_t qe = qepix.GetQECascadesFFactor();
985
986	mean += conv/qe;
987	mean2 += conv*conv/qe/qe;
988	nr ++;
989
990	}
991
992	TArrayF arr(2);
993	arr[0] = nr ? mean/nr : -1;
994	arr[1] = nr>1 ? TMath::Sqrt((mean2 - mean*mean/nr)/(nr-1)) : 0;
995	return arr;
996	}
997
998	// --------------------------------------------------------------------------
999	//
1000	// Calculates the average mean arrival times for pixel sizes.
1001	// The geometry container is used to get the necessary
1002	// geometry information (area index).
1003	// If the bad pixel container is given all pixels which have the flag 'kUnsuitableRun' are ignored
1004	// in the calculation of the size average.
1005	//
1006	// Returns a TArrayF of dimension 2:
1007	// arr[0]: averaged mean arrival times (default: -1.)
1008	// arr[1]: Error (rms) of averaged mean arrival times (default: 0.)
1009	//
1010	TArrayF MCalibrationChargeCam::GetAveragedArrivalTimeMeanPerArea(const MGeomCam &geom,
1011	const UInt_t ai, MBadPixelsCam *bad)
1012	{
1013
1014	const Int_t np = GetSize();
1015
1016	Double_t mean = 0.;
1017	Double_t mean2 = 0.;
1018	Int_t nr = 0;
1019
1020	for (int i=0; i<np; i++)
1021	{
1022	if (bad && (*bad)[i].IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
1023	continue;
1024
1025	const UInt_t aidx = geom[i].GetAidx();
1026
1027	if (ai != aidx)
1028	continue;
1029
1030	const MCalibrationChargePix &pix = (MCalibrationChargePix&)(*this)[i];
1031	const Float_t time = pix.GetAbsTimeMean();
1032
1033	mean += time ;
1034	mean2 += time*time;
1035	nr ++;
1036
1037	}
1038
1039	TArrayF arr(2);
1040	arr[0] = nr ? mean/nr : -1;
1041	arr[1] = nr>1 ? TMath::Sqrt((mean2 - mean*mean/nr)/(nr-1)) : 0;
1042	return arr;
1043	}
1044
1045	// --------------------------------------------------------------------------
1046	//
1047	// Calculates the average mean arrival times for camera sectors.
1048	// The geometry container is used to get the necessary
1049	// geometry information (area index).
1050	// If the bad pixel container is given all pixels which have the flag 'kUnsuitableRun' are ignored
1051	// in the calculation of the size average.
1052	//
1053	// Returns a TArrayF of dimension 2:
1054	// arr[0]: averaged mean arrival times (default: -1.)
1055	// arr[1]: Error (rms) of averaged mean arrival times (default: 0.)
1056	//
1057	TArrayF MCalibrationChargeCam::GetAveragedArrivalTimeMeanPerSector(const MGeomCam &geom,
1058	const UInt_t sec, MBadPixelsCam *bad)
1059	{
1060	const Int_t np = GetSize();
1061
1062	Double_t mean = 0.;
1063	Double_t mean2 = 0.;
1064	Int_t nr = 0;
1065
1066	for (int i=0; i<np; i++)
1067	{
1068	if (bad && (*bad)[i].IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
1069	continue;
1070
1071	const UInt_t sector = geom[i].GetSector();
1072
1073	if (sector != sec)
1074	continue;
1075
1076	const MCalibrationChargePix &pix = (MCalibrationChargePix&)(*this)[i];
1077	const Float_t time = pix.GetAbsTimeMean();
1078
1079	mean += time;
1080	mean2 += time*time;
1081	nr ++;
1082
1083	}
1084
1085	TArrayF arr(2);
1086	arr[0] = nr ? mean/nr : -1;
1087	arr[1] = nr>1 ? TMath::Sqrt((mean2 - mean*mean/nr)/(nr-1)) : 0;
1088	return arr;
1089	}
1090
1091	// --------------------------------------------------------------------------
1092	//
1093	// Calculates the average arrival time RMSs for pixel sizes.
1094	// The geometry container is used to get the necessary
1095	// geometry information (area index).
1096	// If the bad pixel container is given all pixels which have the flag 'kUnsuitableRun' are ignored
1097	// in the calculation of the size average.
1098	//
1099	// Returns a TArrayF of dimension 2:
1100	// arr[0]: averaged arrival time RMSs (default: -1.)
1101	// arr[1]: Error (rms) of averaged arrival time RMSs (default: 0.)
1102	//
1103	TArrayF MCalibrationChargeCam::GetAveragedArrivalTimeRmsPerArea ( const MGeomCam &geom,
1104	const UInt_t ai, MBadPixelsCam *bad)
1105	{
1106
1107	const Int_t np = GetSize();
1108
1109	Double_t mean = 0.;
1110	Double_t mean2 = 0.;
1111	Int_t nr = 0;
1112
1113	for (int i=0; i<np; i++)
1114	{
1115	if (bad && (*bad)[i].IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
1116	continue;
1117
1118	const UInt_t aidx = geom[i].GetAidx();
1119
1120	if (ai != aidx)
1121	continue;
1122
1123	const MCalibrationChargePix &pix = (MCalibrationChargePix&)(*this)[i];
1124	const Float_t rms = pix.GetAbsTimeRms();
1125
1126	mean += rms;
1127	mean2 += rms*rms;
1128	nr ++;
1129
1130	}
1131
1132	TArrayF arr(2);
1133	arr[0] = nr ? mean/nr : -1;
1134	arr[1] = nr>1 ? TMath::Sqrt((mean2 - mean*mean/nr)/(nr-1)) : 0;
1135	return arr;
1136	}
1137
1138	// --------------------------------------------------------------------------
1139	//
1140	// Calculates the average arrival time RMSs for camera sectors.
1141	// The geometry container is used to get the necessary
1142	// geometry information (area index).
1143	// If the bad pixel container is given all pixels which have the flag 'kUnsuitableRun' are ignored
1144	// in the calculation of the size average.
1145	//
1146	// Returns a TArrayF of dimension 2:
1147	// arr[0]: averaged arrival time RMSs (default: -1.)
1148	// arr[1]: Error (rms) of averaged arrival time RMSs (default: 0.)
1149	//
1150	TArrayF MCalibrationChargeCam::GetAveragedArrivalTimeRmsPerSector( const MGeomCam &geom, const UInt_t sec, MBadPixelsCam *bad)
1151	{
1152	const Int_t np = GetSize();
1153
1154	Double_t mean = 0.;
1155	Double_t mean2 = 0.;
1156	Int_t nr = 0;
1157
1158	for (int i=0; i<np; i++)
1159	{
1160	if (bad && (*bad)[i].IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
1161	continue;
1162
1163	const UInt_t sector = geom[i].GetSector();
1164
1165	if (sector != sec)
1166	continue;
1167
1168	const MCalibrationChargePix &pix = (MCalibrationChargePix&)(*this)[i];
1169	const Float_t rms = pix.GetAbsTimeRms();
1170
1171	mean += rms;
1172	mean2 += rms*rms;
1173	nr ++;
1174	}
1175
1176	TArrayF arr(2);
1177	arr[0] = nr ? mean/nr : -1;
1178	arr[1] = nr>1 ? TMath::Sqrt((mean2 - mean*mean/nr)/(nr-1)) : 0;
1179
1180	return arr;
1181	}

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

Context Navigation

source: branches/Corsika7405Compatibility/mcalib/MCalibrationChargeCam.cc@ 19365

Download in other formats: