Context Navigation

source: trunk/MagicSoft/Mars/mcalib/MCalibrationCam.cc@ 3065

Visit:

Last change on this file since 3065 was 3065, checked in by gaug, 22 years ago
* empty log message *
File size: 32.4 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	!
19	! Author(s): Markus Gaug 11/2003 <mailto:markus@ifae.es>
20	!
21	! Copyright: MAGIC Software Development, 2000-2001
22	!
23	!
24	\* ======================================================================== */
25
26	/////////////////////////////////////////////////////////////////////////////
27	//
28	// MCalibrationCam
29	//
30	// Hold the whole Calibration results of the camera:
31	//
32	// 1) MCalibrationCam initializes a TClonesArray whose elements are
33	// pointers to MCalibrationPix Containers
34	// 2) It initializes a pointer to an MCalibrationBlindPix container
35	// 3) It initializes a pointer to an MCalibrationPINDiode container
36	//
37	/////////////////////////////////////////////////////////////////////////////
38	#include "MCalibrationCam.h"
39
40	#include <TH2.h>
41	#include <TCanvas.h>
42	#include <TClonesArray.h>
43
44	#include "MLog.h"
45	#include "MLogManip.h"
46
47	#include "MGeomCam.h"
48	#include "MGeomPix.h"
49
50	#include "MCalibrationPix.h"
51	#include "MCalibrationConfig.h"
52	#include "MCalibrationBlindPix.h"
53	#include "MCalibrationPINDiode.h"
54
55	#include "MHCalibrationPixel.h"
56
57	ClassImp(MCalibrationCam);
58
59	using namespace std;
60
61	const Int_t MCalibrationCam::gkBlindPixelId = 559;
62	const Int_t MCalibrationCam::gkPINDiodeId = 9999;
63
64	// --------------------------------------------------------------------------
65	//
66	// Default constructor.
67	//
68	// Creates a TClonesArray of MCalibrationPix containers, initialized to 1 entry
69	// Later, a call to MCalibrationCam::InitSize(Int_t size) has to be performed
70	//
71	// Creates an MCalibrationBlindPix container
72	// Creates an MCalibrationPINDiode container
73	//
74	MCalibrationCam::MCalibrationCam(const char name, const char title)
75	: fOffsets(NULL),
76	fSlopes(NULL),
77	fOffvsSlope(NULL)
78	{
79	fName = name ? name : "MCalibrationCam";
80	fTitle = title ? title : "Storage container for the Calibration Information in the camera";
81
82	fPixels = new TClonesArray("MCalibrationPix",1);
83	fBlindPixel = new MCalibrationBlindPix();
84	fPINDiode = new MCalibrationPINDiode();
85
86	Clear();
87	}
88
89	// --------------------------------------------------------------------------
90	//
91	// Delete the TClonesArray of MCalibrationPix containers
92	// Delete the MCalibrationPINDiode and the MCalibrationBlindPix
93	//
94	// Delete the histograms if they exist
95	//
96	MCalibrationCam::~MCalibrationCam()
97	{
98
99	//
100	// delete fPixels should delete all Objects stored inside
101	//
102	delete fPixels;
103	delete fBlindPixel;
104	delete fPINDiode;
105
106	if (fOffsets)
107	delete fOffsets;
108	if (fSlopes)
109	delete fSlopes;
110	if (fOffvsSlope)
111	delete fOffvsSlope;
112
113	}
114
115	// -------------------------------------------------------------------
116	//
117	// This function simply allocates memory via the ROOT command:
118	// (TObject*) TStorage::ReAlloc(fCont, newSize sizeof(TObject*),
119	// fSize * sizeof(TObject*));
120	// newSize corresponds to size in our case
121	// fSize is the old size (in most cases: 1)
122	//
123	void MCalibrationCam::InitSize(const UInt_t i)
124	{
125
126	//
127	// check if we have already initialized to size
128	//
129	if (CheckBounds(i))
130	return;
131
132	fPixels->ExpandCreate(i);
133
134	}
135
136	// --------------------------------------------------------------------------
137	//
138	// This function returns the current size of the TClonesArray
139	// independently if the MCalibrationPix is filled with values or not.
140	//
141	// It is the size of the array fPixels.
142	//
143	Int_t MCalibrationCam::GetSize() const
144	{
145	return fPixels->GetEntriesFast();
146	}
147
148	// --------------------------------------------------------------------------
149	//
150	// Check if position i is inside the current bounds of the TClonesArray
151	//
152	Bool_t MCalibrationCam::CheckBounds(Int_t i) const
153	{
154	return i < GetSize();
155	}
156
157
158	// --------------------------------------------------------------------------
159	//
160	// Get i-th pixel (pixel number)
161	//
162	MCalibrationPix &MCalibrationCam::operator[](Int_t i)
163	{
164	return static_cast<MCalibrationPix>(fPixels->UncheckedAt(i));
165	}
166
167	// --------------------------------------------------------------------------
168	//
169	// Get i-th pixel (pixel number)
170	//
171	MCalibrationPix &MCalibrationCam::operator[](Int_t i) const
172	{
173	return static_cast<MCalibrationPix>(fPixels->UncheckedAt(i));
174	}
175
176
177	// --------------------------------------
178	//
179	void MCalibrationCam::Clear(Option_t *o)
180	{
181
182	fPixels->ForEach(TObject, Clear)();
183	fBlindPixel->Clear();
184	fPINDiode->Clear();
185
186	fMeanPhotInsidePlexiglass = -1.;
187	fMeanPhotErrInsidePlexiglass = -1.;
188	fMeanPhotOutsidePlexiglass = -1.;
189	fMeanPhotErrOutsidePlexiglass = -1.;
190
191	fNumExcludedPixels = 0;
192
193	CLRBIT(fFlags,kBlindPixelMethodValid);
194	CLRBIT(fFlags,kPINDiodeMethodValid);
195	CLRBIT(fFlags,kNumPhotInsidePlexiglassAvailable);
196	CLRBIT(fFlags,kNumPhotOutsidePlexiglassAvailable);
197
198	return;
199	}
200
201	void MCalibrationCam::SetBlindPixelMethodValid(const Bool_t b)
202	{
203
204	if (b)
205	SETBIT(fFlags, kBlindPixelMethodValid);
206	else
207	CLRBIT(fFlags, kBlindPixelMethodValid);
208
209	}
210
211	void MCalibrationCam::SetPINDiodeMethodValid(const Bool_t b)
212	{
213
214	if (b)
215	SETBIT(fFlags, kPINDiodeMethodValid);
216	else
217	CLRBIT(fFlags, kPINDiodeMethodValid);
218
219
220	}
221
222	Bool_t MCalibrationCam::IsBlindPixelMethodValid() const
223	{
224	return TESTBIT(fFlags,kBlindPixelMethodValid);
225	}
226
227	Bool_t MCalibrationCam::IsPINDiodeMethodValid() const
228	{
229	return TESTBIT(fFlags,kPINDiodeMethodValid);
230	}
231
232
233	Bool_t MCalibrationCam::IsNumPhotInsidePlexiglassAvailable() const
234	{
235	return TESTBIT(fFlags,kNumPhotInsidePlexiglassAvailable);
236	}
237
238	Bool_t MCalibrationCam::IsNumPhotOutsidePlexiglassAvailable() const
239	{
240	return TESTBIT(fFlags,kNumPhotOutsidePlexiglassAvailable);
241	}
242
243
244
245	// --------------------------------------------------------------------------
246	//
247	// Print first the well fitted pixels
248	// and then the ones which are not FitValid
249	//
250	void MCalibrationCam::Print(Option_t *o) const
251	{
252
253	*fLog << all << GetDescriptor() << ":" << endl;
254	int id = 0;
255
256	*fLog << all << "Succesfully calibrated pixels:" << endl;
257	*fLog << all << endl;
258
259	TIter Next(fPixels);
260	MCalibrationPix *pix;
261	while ((pix=(MCalibrationPix*)Next()))
262	{
263
264	if (pix->IsChargeValid() && !pix->IsExcluded())
265	{
266
267	*fLog << all << pix->GetPixId() << " Pedestals: " << pix->GetPed() << " +- "
268	<< pix->GetPedRms() << " Reduced Charge: " << pix->GetCharge() << " +- "
269	<< pix->GetSigmaCharge() << " Reduced Sigma: " << pix->GetRSigmaCharge()
270	<< " Nr Phe's: " << pix->GetPheFFactorMethod() << endl;
271	id++;
272	}
273	}
274
275	*fLog << all << id << " succesful pixels :-))" << endl;
276	id = 0;
277
278	*fLog << all << endl;
279	*fLog << all << "Pixels with errors:" << endl;
280	*fLog << all << endl;
281
282	TIter Next2(fPixels);
283	while ((pix=(MCalibrationPix*)Next2()))
284	{
285
286	if (!pix->IsChargeValid() && !pix->IsExcluded())
287	{
288
289	*fLog << all << pix->GetPixId() << " Pedestals: " << pix->GetPed() << " +- "
290	<< pix->GetPedRms() << " Reduced Charge: " << pix->GetCharge() << " +- "
291	<< pix->GetSigmaCharge() << " Reduced Sigma: " << pix->GetRSigmaCharge() << endl;
292	id++;
293	}
294	}
295	*fLog << all << id << " pixels with errors :-((" << endl;
296
297	*fLog << all << endl;
298	*fLog << all << "Excluded pixels:" << endl;
299	*fLog << all << endl;
300
301	TIter Next3(fPixels);
302	while ((pix=(MCalibrationPix*)Next3()))
303	if (pix->IsExcluded())
304	*fLog << all << pix->GetPixId() << endl;
305
306	*fLog << all << fNumExcludedPixels << " excluded pixels " << endl;
307	}
308
309	// --------------------------------------------------------------------------
310	//
311	// Return true if pixel is inside bounds of the TClonesArray fPixels
312	//
313	Bool_t MCalibrationCam::IsPixelUsed(Int_t idx) const
314	{
315	if (!CheckBounds(idx))
316	return kFALSE;
317
318	return kTRUE;
319	}
320
321	// --------------------------------------------------------------------------
322	//
323	// Return true if pixel has already been fitted once (independent of the result)
324	//
325	Bool_t MCalibrationCam::IsPixelFitted(Int_t idx) const
326	{
327
328	if (!CheckBounds(idx))
329	return kFALSE;
330
331	return (*this)[idx].IsFitted();
332	}
333
334	// --------------------------------------------------------------------------
335	//
336	// Sets the user ranges of all histograms such that
337	// empty bins at the edges are not used. Additionally, it rebins the
338	// histograms such that in total, 50 bins are used.
339	//
340	void MCalibrationCam::CutEdges()
341	{
342
343	fBlindPixel->GetHist()->CutAllEdges();
344	fPINDiode->GetHist()->CutAllEdges();
345
346	TIter Next(fPixels);
347	MCalibrationPix *pix;
348	while ((pix=(MCalibrationPix*)Next()))
349	{
350	pix->GetHist()->CutAllEdges();
351	}
352
353	return;
354	}
355
356
357	// The types are as follows:
358	//
359	// 0: Fitted Charge
360	// 1: Error of fitted Charge
361	// 2: Sigma of fitted Charge
362	// 3: Error of Sigma of fitted Charge
363	// 4: Returned probability of Gauss fit to Charge distribution
364	// 5: Mean arrival time
365	// 6: Sigma of the arrival time
366	// 7: Chi-square of the Gauss fit to the arrival times
367	// 8: Pedestal
368	// 9: Pedestal RMS
369	// 10: Reduced Sigma Square
370	// 11: Number of Photo-electrons after the F-Factor method
371	// 12: Error on the Number of Photo-electrons after the F-Factor method
372	// 13: Mean conversion factor after the F-Factor method
373	// 14: Error on the conversion factor after the F-Factor method
374	// 15: Number of Photons after the Blind Pixel method
375	// 16: Mean conversion factor after the Blind Pixel method
376	//
377	Bool_t MCalibrationCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const
378	{
379
380	if (idx > GetSize())
381	return kFALSE;
382
383	Float_t arearatio = cam.GetPixRatio(idx);
384
385	if (arearatio == 0)
386	return kFALSE;
387
388	switch (type)
389	{
390	case 0:
391	if ((*this)[idx].IsExcluded())
392	return kFALSE;
393	val = (*this)[idx].GetCharge();
394	break;
395	case 1:
396	if ((*this)[idx].IsExcluded())
397	return kFALSE;
398	val = (*this)[idx].GetErrCharge();
399	break;
400	case 2:
401	if ((*this)[idx].IsExcluded())
402	return kFALSE;
403	val = (*this)[idx].GetSigmaCharge();
404	break;
405	case 3:
406	if ((*this)[idx].IsExcluded())
407	return kFALSE;
408	val = (*this)[idx].GetErrSigmaCharge();
409	break;
410	case 4:
411	if ((*this)[idx].IsExcluded())
412	return kFALSE;
413	val = (*this)[idx].GetChargeProb();
414	break;
415	case 5:
416	if ((*this)[idx].IsExcluded())
417	return kFALSE;
418	val = (*this)[idx].GetRSigmaCharge();
419	break;
420	case 6:
421	if ((*this)[idx].IsExcluded())
422	return kFALSE;
423	val = (*this)[idx].GetErrRSigmaCharge();
424	break;
425	case 7:
426	if ((*this)[idx].IsExcluded())
427	return kFALSE;
428	val = (this)[idx].GetRSigmaCharge() / (this)[idx].GetCharge();
429	break;
430	case 8:
431	if ((*this)[idx].IsExcluded())
432	return kFALSE;
433	// relative error RsigmaCharge square
434	val = (this)[idx].GetErrRSigmaCharge() (*this)[idx].GetErrRSigmaCharge()
435	/ ((this)[idx].GetRSigmaCharge() (*this)[idx].GetRSigmaCharge() );
436	// relative error Charge square
437	val += (this)[idx].GetErrCharge() (*this)[idx].GetErrCharge()
438	/ ((this)[idx].GetCharge() (*this)[idx].GetCharge() );
439	// calculate relative error out of squares
440	val = TMath::Sqrt(val) ;
441	// multiply with value to get absolute error
442	val = (this)[idx].GetRSigmaCharge() / (*this)[idx].GetCharge();
443	break;
444	case 9:
445	if ((*this)[idx].IsExcluded())
446	return kFALSE;
447	val = (*this)[idx].GetPheFFactorMethod();
448	break;
449	case 10:
450	if ((*this)[idx].IsExcluded())
451	return kFALSE;
452	val = (*this)[idx].GetPheFFactorMethodError();
453	break;
454	case 11:
455	if ((*this)[idx].IsExcluded())
456	return kFALSE;
457	val = (*this)[idx].GetMeanConversionFFactorMethod();
458	break;
459	case 12:
460	if ((*this)[idx].IsExcluded())
461	return kFALSE;
462	val = (*this)[idx].GetErrorConversionFFactorMethod();
463	break;
464	case 13:
465	if ((*this)[idx].IsExcluded())
466	return kFALSE;
467	val = (*this)[idx].GetTotalFFactorFFactorMethod();
468	break;
469	case 14:
470	if ((*this)[idx].IsExcluded())
471	return kFALSE;
472	val = (*this)[idx].GetTotalFFactorErrorFFactorMethod();
473	break;
474	case 15:
475	if ((*this)[idx].IsExcluded())
476	return kFALSE;
477	if (arearatio == 1.)
478	val = GetMeanPhotInsidePlexiglass();
479	else
480	val = GetMeanPhotInsidePlexiglass()*gkCalibrationOutervsInnerPixelArea;
481	break;
482	case 16:
483	if ((*this)[idx].IsExcluded())
484	return kFALSE;
485	if (arearatio == 1.)
486	val = (double)fMeanPhotInsidePlexiglass;
487	else
488	val = (double)fMeanPhotInsidePlexiglass*gkCalibrationOutervsInnerPixelArea;
489	break;
490	case 17:
491	if ((*this)[idx].IsExcluded())
492	return kFALSE;
493	if (arearatio == 1.)
494	val = (*this)[idx].GetMeanConversionBlindPixelMethod();
495	else
496	val = (this)[idx].GetMeanConversionBlindPixelMethod()gkCalibrationOutervsInnerPixelArea;
497	break;
498	case 18:
499	if ((*this)[idx].IsExcluded())
500	return kFALSE;
501	if (arearatio == 1.)
502	val = (*this)[idx].GetErrorConversionBlindPixelMethod();
503	else
504	{
505	val = (this)[idx].GetErrorConversionBlindPixelMethod()(*this)[idx].GetErrorConversionBlindPixelMethod()
506	* gkCalibrationOutervsInnerPixelArea * gkCalibrationOutervsInnerPixelArea;
507	val += gkCalibrationOutervsInnerPixelAreaError * gkCalibrationOutervsInnerPixelAreaError
508	* (this)[idx].GetMeanConversionBlindPixelMethod() (*this)[idx].GetMeanConversionBlindPixelMethod();
509	val = TMath::Sqrt(val);
510	}
511	break;
512	case 19:
513	if ((*this)[idx].IsExcluded())
514	return kFALSE;
515	val = (*this)[idx].GetTotalFFactorBlindPixelMethod();
516	break;
517	case 20:
518	if ((*this)[idx].IsExcluded())
519	return kFALSE;
520	val = (*this)[idx].GetTotalFFactorErrorBlindPixelMethod();
521	break;
522	case 21:
523	if ((*this)[idx].IsExcluded())
524	return kFALSE;
525	if (arearatio == 1.)
526	val = GetMeanPhotOutsidePlexiglass();
527	else
528	val = GetMeanPhotOutsidePlexiglass()*gkCalibrationOutervsInnerPixelArea;
529	break;
530	case 22:
531	if ((*this)[idx].IsExcluded())
532	return kFALSE;
533	if (arearatio == 1.)
534	val = (double)fMeanPhotOutsidePlexiglass;
535	else
536	val = (double)fMeanPhotOutsidePlexiglass*gkCalibrationOutervsInnerPixelArea;
537	break;
538	case 23:
539	if ((*this)[idx].IsExcluded())
540	return kFALSE;
541	if (arearatio == 1.)
542	val = (*this)[idx].GetMeanConversionPINDiodeMethod();
543	else
544	val = (this)[idx].GetMeanConversionPINDiodeMethod()gkCalibrationOutervsInnerPixelArea;
545	break;
546	case 24:
547	if ((*this)[idx].IsExcluded())
548	return kFALSE;
549	if (arearatio == 1.)
550	val = (*this)[idx].GetErrorConversionPINDiodeMethod();
551	else
552	{
553	val = (this)[idx].GetErrorConversionPINDiodeMethod()(*this)[idx].GetErrorConversionPINDiodeMethod()
554	* gkCalibrationOutervsInnerPixelArea * gkCalibrationOutervsInnerPixelArea;
555	val += gkCalibrationOutervsInnerPixelAreaError * gkCalibrationOutervsInnerPixelAreaError
556	* (this)[idx].GetMeanConversionPINDiodeMethod() (*this)[idx].GetMeanConversionPINDiodeMethod();
557	val = TMath::Sqrt(val);
558	}
559	break;
560	case 25:
561	if ((*this)[idx].IsExcluded())
562	return kFALSE;
563	val = (*this)[idx].GetTotalFFactorBlindPixelMethod();
564	break;
565	case 26:
566	if ((*this)[idx].IsExcluded())
567	return kFALSE;
568	val = (*this)[idx].GetTotalFFactorErrorBlindPixelMethod();
569	break;
570	case 27:
571	if ((*this)[idx].IsExcluded())
572	val = 1.;
573	else
574	return kFALSE;
575	break;
576	case 28:
577	if ((*this)[idx].IsExcluded())
578	return kFALSE;
579	if (!(*this)[idx].IsFitted())
580	val = 1;
581	else
582	return kFALSE;
583	break;
584	case 29:
585	if ((*this)[idx].IsExcluded())
586	return kFALSE;
587	if (!(*this)[idx].IsFitted())
588	return kFALSE;
589	if (!(*this)[idx].IsChargeValid())
590	val = 1;
591	else
592	return kFALSE;
593	break;
594	case 30:
595	if ((*this)[idx].IsExcluded())
596	return kFALSE;
597	if ((*this)[idx].IsOscillating())
598	val = 1;
599	else
600	return kFALSE;
601	break;
602	case 31:
603	if ((*this)[idx].IsExcluded())
604	return kFALSE;
605	if ((*this)[idx].IsHiGainSaturation())
606	val = 1;
607	else
608	return kFALSE;
609	break;
610	case 32:
611	if ((*this)[idx].IsExcluded())
612	return kFALSE;
613	if ((*this)[idx].IsFFactorMethodValid())
614	val = 1;
615	else
616	return kFALSE;
617	break;
618	case 33:
619	if ((*this)[idx].IsExcluded())
620	return kFALSE;
621	if ((*this)[idx].IsBlindPixelMethodValid())
622	val = 1;
623	else
624	return kFALSE;
625	break;
626	case 34:
627	if ((*this)[idx].IsExcluded())
628	return kFALSE;
629	if ((*this)[idx].IsPINDiodeMethodValid())
630	val = 1;
631	else
632	return kFALSE;
633	break;
634	case 35:
635	if ((*this)[idx].IsExcluded())
636	return kFALSE;
637	val = (*this)[idx].GetPed();
638	break;
639	case 36:
640	if ((*this)[idx].IsExcluded())
641	return kFALSE;
642	val = 1.;
643	// val = (*this)[idx].GetPedError();
644	break;
645	case 37:
646	if ((*this)[idx].IsExcluded())
647	return kFALSE;
648	val = (*this)[idx].GetPedRms();
649	break;
650	case 38:
651	if ((*this)[idx].IsExcluded())
652	return kFALSE;
653	val = 1.;
654	// val = (*this)[idx].GetPedRmsError();
655	break;
656	case 39:
657	if ((*this)[idx].IsExcluded())
658	return kFALSE;
659	val = (*this)[idx].GetAbsTimeMean();
660	break;
661	case 40:
662	if ((*this)[idx].IsExcluded())
663	return kFALSE;
664	val = (*this)[idx].GetAbsTimeMeanErr();
665	break;
666	case 41:
667	if ((*this)[idx].IsExcluded())
668	return kFALSE;
669	val = (*this)[idx].GetAbsTimeRms();
670	break;
671	case 42:
672	if ((*this)[idx].IsExcluded())
673	return kFALSE;
674	val = (*this)[idx].GetAbsTimeMeanErr()/TMath::Sqrt(2.);
675	break;
676	default:
677	return kFALSE;
678	}
679	return val!=-1.;
680	}
681
682	// --------------------------------------------------------------------------
683	//
684	// What MHCamera needs in order to draw an individual pixel in the camera
685	//
686	void MCalibrationCam::DrawPixelContent(Int_t idx) const
687	{
688	(*this)[idx].Draw();
689	}
690
691
692	// --------------------------------------------------------------------------
693	//
694	//
695	//
696	Bool_t MCalibrationCam::CalcNumPhotInsidePlexiglass()
697	{
698
699	if (!fBlindPixel->IsFitOK())
700	return kFALSE;
701
702	const Float_t mean = fBlindPixel->GetLambda();
703	const Float_t merr = fBlindPixel->GetErrLambda();
704
705	//
706	// Start calculation of number of photons
707	//
708	// The blind pixel has exactly one cm^2 area (with negligible error),
709	// thus we omi division by 1.
710	//
711	fMeanPhotInsidePlexiglass = mean * gkCalibrationInnerPixelArea;
712
713	// Start calculation of number of photons relative Variance (!!)
714	fMeanPhotErrInsidePlexiglass = merr*merr/mean/mean;
715	fMeanPhotErrInsidePlexiglass += gkCalibrationInnerPixelAreaError*gkCalibrationInnerPixelAreaError
716	/ gkCalibrationInnerPixelArea/gkCalibrationInnerPixelArea;
717
718	switch (fColor)
719	{
720	case kECGreen:
721	fMeanPhotInsidePlexiglass /= gkCalibrationBlindPixelQEGreen;
722	fMeanPhotErrInsidePlexiglass += gkCalibrationBlindPixelQEGreenError*gkCalibrationBlindPixelQEGreenError
723	/ gkCalibrationBlindPixelQEGreen / gkCalibrationBlindPixelQEGreen;
724
725	fMeanPhotInsidePlexiglass *= TMath::Power(10,gkCalibrationBlindPixelAttGreen); // correct for absorption
726	// attenuation has negligible error
727	break;
728	case kECBlue:
729	fMeanPhotInsidePlexiglass /= gkCalibrationBlindPixelQEBlue;
730	fMeanPhotErrInsidePlexiglass += gkCalibrationBlindPixelQEBlueError*gkCalibrationBlindPixelQEBlueError
731	/ gkCalibrationBlindPixelQEBlue / gkCalibrationBlindPixelQEBlue;
732
733	fMeanPhotInsidePlexiglass *= TMath::Power(10,gkCalibrationBlindPixelAttBlue); // correct for absorption
734	// attenuation has negligible error
735	break;
736	case kECUV:
737	fMeanPhotInsidePlexiglass /= gkCalibrationBlindPixelQEUV;
738	fMeanPhotErrInsidePlexiglass += gkCalibrationBlindPixelQEUVError*gkCalibrationBlindPixelQEUVError
739	/ gkCalibrationBlindPixelQEUV / gkCalibrationBlindPixelQEUV;
740
741	fMeanPhotInsidePlexiglass *= TMath::Power(10,gkCalibrationBlindPixelAttUV); // correct for absorption
742	// attenuation has negligible error
743	break;
744	case kECCT1:
745	default:
746	fMeanPhotInsidePlexiglass /= gkCalibrationBlindPixelQECT1;
747	fMeanPhotErrInsidePlexiglass += gkCalibrationBlindPixelQECT1Error*gkCalibrationBlindPixelQECT1Error
748	/ gkCalibrationBlindPixelQECT1 / gkCalibrationBlindPixelQECT1;
749
750	fMeanPhotInsidePlexiglass *= TMath::Power(10,gkCalibrationBlindPixelAttCT1); // correct for absorption
751	// attenuation has negligible error
752	break;
753	}
754
755	*fLog << inf << endl;
756	*fLog << inf << " Mean number of Photons for an Inner Pixel (inside Plexiglass): "
757	<< fMeanPhotInsidePlexiglass << endl;
758
759	if (fMeanPhotInsidePlexiglass > 0.)
760	SETBIT(fFlags,kNumPhotInsidePlexiglassAvailable);
761	else
762	{
763	CLRBIT(fFlags,kNumPhotInsidePlexiglassAvailable);
764	return kFALSE;
765	}
766
767	if (fMeanPhotErrInsidePlexiglass < 0.)
768	{
769	*fLog << warn << " Relative Variance on number of Photons for an Inner Pixel (inside Plexiglass): "
770	<< fMeanPhotErrInsidePlexiglass << endl;
771	CLRBIT(fFlags,kNumPhotInsidePlexiglassAvailable);
772	return kFALSE;
773	}
774
775	// Finish calculation of errors -> convert from relative variance to absolute error
776	fMeanPhotErrInsidePlexiglass = TMath::Sqrt(fMeanPhotErrInsidePlexiglass);
777	fMeanPhotErrInsidePlexiglass *= fMeanPhotInsidePlexiglass;
778
779	*fLog << inf << " Error on number of Photons for an Inner Pixel (inside Plexiglass): "
780	<< fMeanPhotErrInsidePlexiglass << endl;
781	*fLog << inf << endl;
782
783
784
785	TIter Next(fPixels);
786	MCalibrationPix *pix;
787	while ((pix=(MCalibrationPix*)Next()))
788	{
789
790	if(pix->IsChargeValid())
791	{
792
793	const Float_t charge = pix->GetCharge();
794	const Float_t ratio = fGeomCam->GetPixRatio(pix->GetPixId());
795	const Float_t chargeerr = pix->GetErrCharge();
796
797	Float_t nphot;
798	Float_t nphoterr;
799
800	if (ratio == 1.)
801	{
802	nphot = fMeanPhotInsidePlexiglass;
803	nphoterr = fMeanPhotErrInsidePlexiglass;
804	}
805	else
806	{
807	nphot = fMeanPhotInsidePlexiglass*gkCalibrationOutervsInnerPixelArea;
808	nphoterr = fMeanPhotErrInsidePlexiglass*gkCalibrationOutervsInnerPixelArea
809	fMeanPhotErrInsidePlexiglassgkCalibrationOutervsInnerPixelArea;
810	nphoterr += fMeanPhotInsidePlexiglass*gkCalibrationOutervsInnerPixelAreaError
811	fMeanPhotInsidePlexiglassgkCalibrationOutervsInnerPixelAreaError;
812	nphoterr = TMath::Sqrt(nphoterr);
813	}
814
815
816
817	const Float_t conversion = nphot/charge;
818	Float_t conversionerr;
819
820	conversionerr = nphoterr/charge
821	* nphoterr/charge ;
822	conversionerr += chargeerr/charge
823	* chargeerr/charge
824	* conversion*conversion;
825	conversionerr = TMath::Sqrt(conversionerr);
826
827	const Float_t conversionsigma = 0.;
828
829	pix->SetConversionBlindPixelMethod(conversion, conversionerr, conversionsigma);
830
831	if (conversionerr/conversion < 0.1)
832	pix->SetBlindPixelMethodValid();
833	}
834	}
835	return kTRUE;
836	}
837
838
839	Bool_t MCalibrationCam::CalcNumPhotOutsidePlexiglass()
840	{
841
842	if (!fPINDiode->IsChargeFitValid())
843	return kFALSE;
844
845	const Float_t mean = fPINDiode->GetCharge();
846	const Float_t merr = fPINDiode->GetErrCharge();
847
848	// Start calculation of number of photons
849	fMeanPhotOutsidePlexiglass = mean * gkCalibrationInnerPixelvsPINDiodeArea;
850
851	// Start calculation of number of photons relative Variance (!!)
852	fMeanPhotErrOutsidePlexiglass = merr*merr/mean/mean;
853	fMeanPhotErrOutsidePlexiglass += gkCalibrationInnerPixelvsPINDiodeAreaError*gkCalibrationInnerPixelvsPINDiodeAreaError
854	/ gkCalibrationInnerPixelvsPINDiodeArea/gkCalibrationInnerPixelvsPINDiodeArea;
855
856	switch (fColor)
857	{
858	case kECGreen:
859	fMeanPhotOutsidePlexiglass /= gkCalibrationPINDiodeQEGreen;
860	fMeanPhotErrOutsidePlexiglass += gkCalibrationPINDiodeQEGreenError*gkCalibrationPINDiodeQEGreenError
861	/ gkCalibrationPINDiodeQEGreen/gkCalibrationPINDiodeQEGreen;
862	break;
863	case kECBlue:
864	fMeanPhotOutsidePlexiglass /= gkCalibrationPINDiodeQEBlue;
865	fMeanPhotErrOutsidePlexiglass += gkCalibrationPINDiodeQEBlueError*gkCalibrationPINDiodeQEBlueError
866	/ gkCalibrationPINDiodeQEBlue/gkCalibrationPINDiodeQEBlue;
867	break;
868	case kECUV:
869	fMeanPhotOutsidePlexiglass /= gkCalibrationPINDiodeQEUV;
870	fMeanPhotErrOutsidePlexiglass += gkCalibrationPINDiodeQEUVError*gkCalibrationPINDiodeQEUVError
871	/ gkCalibrationPINDiodeQEUV/gkCalibrationPINDiodeQEUV;
872	break;
873	case kECCT1:
874	default:
875	fMeanPhotOutsidePlexiglass /= gkCalibrationPINDiodeQECT1;
876	fMeanPhotErrOutsidePlexiglass += gkCalibrationPINDiodeQECT1Error*gkCalibrationPINDiodeQECT1Error
877	/ gkCalibrationPINDiodeQECT1/gkCalibrationPINDiodeQECT1;
878	break;
879	}
880
881
882	*fLog << inf << endl;
883	*fLog << inf << " Mean number of Photons for an Inner Pixel (outside Plexiglass): "
884	<< fMeanPhotOutsidePlexiglass << endl;
885
886	if (fMeanPhotOutsidePlexiglass > 0.)
887	SETBIT(fFlags,kNumPhotOutsidePlexiglassAvailable);
888	else
889	{
890	CLRBIT(fFlags,kNumPhotOutsidePlexiglassAvailable);
891	return kFALSE;
892	}
893
894	if (fMeanPhotErrOutsidePlexiglass < 0.)
895	{
896	*fLog << warn << "Relative Variance on number of Photons for an Inner Pixel (outside Plexiglass): "
897	<< fMeanPhotErrOutsidePlexiglass << endl;
898	CLRBIT(fFlags,kNumPhotOutsidePlexiglassAvailable);
899	return kFALSE;
900	}
901
902	// Finish calculation of errors -> convert from relative variance to absolute error
903	fMeanPhotErrOutsidePlexiglass = TMath::Sqrt(fMeanPhotErrOutsidePlexiglass);
904	fMeanPhotErrOutsidePlexiglass *= fMeanPhotOutsidePlexiglass;
905
906	*fLog << inf << " Error on number of Photons for an Inner Pixel (outside Plexiglass): "
907	<< fMeanPhotErrOutsidePlexiglass << endl;
908	*fLog << inf << endl;
909
910	TIter Next(fPixels);
911	MCalibrationPix *pix;
912	while ((pix=(MCalibrationPix*)Next()))
913	{
914
915	if (pix->IsChargeValid())
916	{
917
918	const Float_t charge = pix->GetCharge();
919	const Float_t ratio = fGeomCam->GetPixRatio(pix->GetPixId());
920	const Float_t chargeerr = pix->GetErrCharge();
921
922	Float_t nphot;
923	Float_t nphoterr;
924
925	if (ratio == 1.)
926	{
927	nphot = fMeanPhotInsidePlexiglass;
928	nphoterr = fMeanPhotErrInsidePlexiglass;
929	}
930	else
931	{
932	nphot = fMeanPhotInsidePlexiglass*gkCalibrationOutervsInnerPixelArea;
933	nphoterr = fMeanPhotErrInsidePlexiglass*gkCalibrationOutervsInnerPixelArea
934	fMeanPhotErrInsidePlexiglassgkCalibrationOutervsInnerPixelArea;
935	nphoterr += fMeanPhotInsidePlexiglass*gkCalibrationOutervsInnerPixelAreaError
936	fMeanPhotInsidePlexiglassgkCalibrationOutervsInnerPixelAreaError;
937	nphoterr = TMath::Sqrt(nphoterr);
938	}
939
940	const Float_t conversion = nphot/charge;
941	Float_t conversionerr;
942
943	conversionerr = nphoterr/charge
944	* nphoterr/charge ;
945	conversionerr += chargeerr/charge
946	* chargeerr/charge
947	* conversion*conversion;
948	if (conversionerr > 0.)
949	conversionerr = TMath::Sqrt(conversionerr);
950
951	const Float_t conversionsigma = 0.;
952
953	pix->SetConversionBlindPixelMethod(conversion, conversionerr, conversionsigma);
954
955	if (conversionerr/conversion < 0.1)
956	pix->SetBlindPixelMethodValid();
957
958	}
959	}
960	return kTRUE;
961	}
962
963
964
965	Bool_t MCalibrationCam::GetConversionFactorBlindPixel(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma)
966	{
967
968	if (ipx < 0 \|\| !IsPixelFitted(ipx))
969	return kFALSE;
970
971	if (!IsNumPhotInsidePlexiglassAvailable())
972	if (!CalcNumPhotInsidePlexiglass())
973	return kFALSE;
974
975	mean = (*this)[ipx].GetMeanConversionBlindPixelMethod();
976	err = (*this)[ipx].GetErrorConversionBlindPixelMethod();
977	sigma = (*this)[ipx].GetSigmaConversionBlindPixelMethod();
978
979	return kTRUE;
980	}
981
982
983	Bool_t MCalibrationCam::GetConversionFactorFFactor(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma)
984	{
985
986	if (ipx < 0 \|\| !IsPixelFitted(ipx))
987	return kFALSE;
988
989	Float_t conv = (*this)[ipx].GetMeanConversionFFactorMethod();
990
991	if (conv < 0.)
992	return kFALSE;
993
994	mean = conv;
995	err = (*this)[ipx].GetErrorConversionFFactorMethod();
996	sigma = (*this)[ipx].GetSigmaConversionFFactorMethod();
997
998	return kTRUE;
999	}
1000
1001
1002	//-----------------------------------------------------------------------------------
1003	//
1004	// Calculates the conversion factor between the integral of FADCs slices
1005	// (as defined in the signal extractor MExtractSignal.cc)
1006	// and the number of photons reaching the plexiglass for one Inner Pixel
1007	//
1008	// FIXME: The PINDiode is still not working and so is the code
1009	//
1010	Bool_t MCalibrationCam::GetConversionFactorPINDiode(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma)
1011	{
1012
1013	if (ipx < 0 \|\| !IsPixelFitted(ipx))
1014	return kFALSE;
1015
1016	if (!IsNumPhotOutsidePlexiglassAvailable())
1017	if (!CalcNumPhotOutsidePlexiglass())
1018	return kFALSE;
1019
1020	mean = (*this)[ipx].GetMeanConversionPINDiodeMethod();
1021	err = (*this)[ipx].GetErrorConversionPINDiodeMethod();
1022	sigma = (*this)[ipx].GetSigmaConversionPINDiodeMethod();
1023
1024	return kFALSE;
1025
1026	}
1027
1028	//-----------------------------------------------------------------------------------
1029	//
1030	// Calculates the best combination of the three used methods possible
1031	// between the integral of FADCs slices
1032	// (as defined in the signal extractor MExtractSignal.cc)
1033	// and the number of photons reaching one Inner Pixel.
1034	// The procedure is not yet defined.
1035	//
1036	// FIXME: The PINDiode is still not working and so is the code
1037	//
1038	Bool_t MCalibrationCam::GetConversionFactorCombined(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma)
1039	{
1040
1041	if (ipx < 0 \|\| !IsPixelFitted(ipx))
1042	return kFALSE;
1043
1044	return kFALSE;
1045
1046	}
1047
1048
1049	void MCalibrationCam::DrawHiLoFits()
1050	{
1051
1052	if (!fOffsets)
1053	fOffsets = new TH1D("pp","Offsets of the HiGain LoGain Fit",100,-600.,400.);
1054	if (!fSlopes)
1055	fSlopes = new TH1D("mm","Slopes of the HiGain LoGain Fit",100,-2.,2.);
1056	if (!fOffvsSlope)
1057	fOffvsSlope = new TH2D("aa","Slopes vs Offsets of the HiGain LoGain Fit",100,-600.,400.,100,-2.,2.);
1058
1059	TIter Next(fPixels);
1060	MCalibrationPix *pix;
1061	MHCalibrationPixel *hist;
1062	while ((pix=(MCalibrationPix*)Next()))
1063	{
1064	hist = pix->GetHist();
1065	hist->FitHiGainvsLoGain();
1066	fOffsets->Fill(hist->GetOffset(),1.);
1067	fSlopes->Fill(hist->GetSlope(),1.);
1068	fOffvsSlope->Fill(hist->GetOffset(),hist->GetSlope(),1.);
1069	}
1070
1071	TCanvas *c1 = new TCanvas();
1072
1073	c1->Divide(1,3);
1074	c1->cd(1);
1075	fOffsets->Draw();
1076	gPad->Modified();
1077	gPad->Update();
1078
1079	c1->cd(2);
1080	fSlopes->Draw();
1081	gPad->Modified();
1082	gPad->Update();
1083
1084	c1->cd(3);
1085	fOffvsSlope->Draw("col1");
1086	gPad->Modified();
1087	gPad->Update();
1088	}
1089

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

Download in other formats: