Context Navigation

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

Visit:

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

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

Download in other formats: