Context Navigation

source: trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargeBlindPix.cc@ 5102

Visit:

Last change on this file since 5102 was 5098, checked in by tbretz, 21 years ago
* empty log message *
File size: 29.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	! Author(s): Markus Gaug 02/2004 <mailto:markus@ifae.es>
19	!
20	! Copyright: MAGIC Software Development, 2000-2004
21	!
22	!
23	\* ======================================================================== */
24
25	//////////////////////////////////////////////////////////////////////////////
26	//
27	// MHCalibrationChargeBlindPix
28	//
29	// Histogram class for the charge calibration of the Blind Pixel.
30	// Stores and fits the charges and stores the averaged assumed pedestal and
31	// single-phe FADC slice entries. Charges are taken from MExtractedSignalBlindPix.
32	// Performs the Single Photo-electron fit to extract the Poisson mean and its errors.
33	//
34	// Different fits can be chosen with the function SetFitFunc().
35	//
36	// The fit result is accepted under the condition that:
37	// 1) the Probability is greater than fProbLimit (default 0.001 == 99.7%)
38	// 2) at least fNumSinglePheLimit events are found in the single Photo-electron peak
39	//
40	// The single FADC slice entries are averaged and stored in fASinglePheFADCSlices, if
41	// their sum exceeds fSinglePheCut, otherwise in fAPedestalFADCSlices.
42	//
43	// Used numbers are the following:
44	//
45	// Electronic conversion factor:
46	// Assume, we have N_e electrons at the anode,
47	// thus a charge of N_e*e (e = electron charge) Coulomb.
48	//
49	// This charge is AC coupled and runs into a R_pre = 50 Ohm resistency.
50	// The corresponding current is amplified by a gain factor G_pre = 400
51	// (the precision of this value still has to be checked !!!) and again AC coupled to
52	// the output.
53	// The corresponding signal goes through the whole transmission and
54	// amplification chain and is digitized in the FADCs.
55	// The conversion Signal Area to FADC counts (Conv_trans) has been measured
56	// by David and Oscar to be approx. 3.9 pVs^-1
57	//
58	// Thus: Conversion FADC counts to Number of Electrons at Anode:
59	// FADC counts = (1/Conv_tran) * G_pre * R_pre * e * N_e = 8 * 10^-4 N_e.
60	//
61	// Also: FADC counts = 810^-4 GAIN * N_phe
62	//
63	// In the blind pixel, there is an additional pre-amplifier with an amplification of
64	// about 10. Therefore, we have for the blind pixel:
65	//
66	//
67	// FADC counts (Blind Pixel) = 810^-3 GAIN * N_phe
68	//
69	//////////////////////////////////////////////////////////////////////////////
70	#include "MHCalibrationChargeBlindPix.h"
71
72	#include <TStyle.h>
73	#include <TCanvas.h>
74	#include <TPaveText.h>
75	#include <TPaveStats.h>
76	#include <TGraph.h>
77
78	#include <TVector.h>
79	#include <TF1.h>
80	#include <TH1.h>
81	#include <TH2D.h>
82	#include <TRandom.h>
83
84	#include "MLog.h"
85	#include "MLogManip.h"
86
87	#include "MParList.h"
88
89	#include "MRawEvtPixelIter.h"
90
91	ClassImp(MHCalibrationChargeBlindPix);
92
93	using namespace std;
94
95	const Double_t MHCalibrationChargeBlindPix::gkElectronicAmp = 0.008;
96	const Double_t MHCalibrationChargeBlindPix::gkElectronicAmpErr = 0.002;
97	const Float_t MHCalibrationChargeBlindPix::gkSignalInitializer = -9999.;
98	const Float_t MHCalibrationChargeBlindPix::fgNumSinglePheLimit = 50.;
99	// --------------------------------------------------------------------------
100	//
101	// Default Constructor.
102	//
103	// Sets:
104	// - the default number for fNumSinglePheLimit (fgNumSinglePheLimit)
105	// - the default number of bins after stripping (30)
106	//
107	// Initializes:
108	// - all pointers to NULL
109	// - fASinglePheFADCSlices(0);
110	// - fAPedestalFADCSlices(0);
111	//
112	// Calls:
113	// - Clear()
114	//
115	MHCalibrationChargeBlindPix::MHCalibrationChargeBlindPix(const char name, const char title)
116	: fSinglePheFit(NULL),
117	fFitLegend(NULL),
118	fHSinglePheFADCSlices(NULL), fHPedestalFADCSlices(NULL)
119	{
120
121	fName = name ? name : "MHCalibrationChargeBlindPix";
122	fTitle = title ? title : "Statistics of the FADC sums of Blind Pixel calibration events";
123
124	fASinglePheFADCSlices.ResizeTo(1);
125	fAPedestalFADCSlices.ResizeTo(1);
126
127	SetNumSinglePheLimit();
128	SetProbLimit(0.001);
129	SetBinsAfterStripping(0);
130
131	Clear();
132	}
133
134	// --------------------------------------------------------------------------
135	//
136	// Default Destructor.
137	//
138	// Deletes (if Pointer is not NULL):
139	//
140	// - fFitLegend
141	// - fHSinglePheFADCSlices
142	// - fHPedestalFADCSlices
143	// - fSinglePheFit
144	//
145	MHCalibrationChargeBlindPix::~MHCalibrationChargeBlindPix()
146	{
147
148	if (fSinglePheFit)
149	delete fSinglePheFit;
150
151	if (fFitLegend)
152	delete fFitLegend;
153
154	if (fHSinglePheFADCSlices)
155	delete fHSinglePheFADCSlices;
156
157	if (fHPedestalFADCSlices)
158	delete fHPedestalFADCSlices;
159
160	}
161
162	// --------------------------------------------------------------------------
163	//
164	// Sets:
165	// - all variables to 0., except the fit result variables to gkSignalInitializer
166	// - all flags to kFALSE
167	// - all pointers to NULL
168	// - the default fit function (kEPoisson5)
169	//
170	// Deletes:
171	// - all pointers unequal NULL
172	//
173	// Calls:
174	// - MHCalibrationChargePix::Clear()
175	//
176	void MHCalibrationChargeBlindPix::Clear(Option_t *o)
177	{
178
179	fLambda = gkSignalInitializer;
180	fMu0 = gkSignalInitializer;
181	fMu1 = gkSignalInitializer;
182	fSigma0 = gkSignalInitializer;
183	fSigma1 = gkSignalInitializer;
184	fLambdaErr = gkSignalInitializer;
185	fMu0Err = gkSignalInitializer;
186	fMu1Err = gkSignalInitializer;
187	fSigma0Err = gkSignalInitializer;
188	fSigma1Err = gkSignalInitializer;
189
190	fLambdaCheck = gkSignalInitializer;
191	fLambdaCheckErr = gkSignalInitializer;
192
193	fFitFunc = kEPoisson5;
194
195	fNumSinglePhes = 0;
196	fNumPedestals = 0;
197
198	fChisquare = 0.;
199	fNDF = 0 ;
200	fProb = 0.;
201
202	SetSinglePheFitOK ( kFALSE );
203	SetPedestalFitOK ( kFALSE );
204
205	if (fFitLegend)
206	{
207	delete fFitLegend;
208	fFitLegend = NULL;
209	}
210
211	if (fSinglePheFit)
212	{
213	delete fSinglePheFit;
214	fSinglePheFit = NULL;
215	}
216
217	if (fHSinglePheFADCSlices)
218	{
219	delete fHSinglePheFADCSlices;
220	fHSinglePheFADCSlices = NULL;
221	}
222
223	if (fHPedestalFADCSlices)
224	{
225	delete fHPedestalFADCSlices;
226	fHPedestalFADCSlices = NULL;
227	}
228
229	MHCalibrationPix::Clear();
230	return;
231	}
232
233	/*
234	// --------------------------------------------------------------------------
235	//
236	// Our own clone function is necessary since root 3.01/06 or Mars 0.4
237	// I don't know the reason.
238	//
239	// Creates new MHCalibrationCam
240	//
241	TObject MHCalibrationChargeBlindPix::Clone(const char ) const
242	{
243
244	MHCalibrationChargeBlindPix *pix = new MHCalibrationChargeBlindPix();
245	this->Copy(*pix);
246
247	this->fHGausHist.Copy(pix->fHGausHist);
248	this->fSinglePheFit->Copy(*(pix->fSinglePheFit));
249	this->fHSinglePheFADCSlices->Copy(*(pix->fHSinglePheFADCSlices));
250	this->fHPedestalFADCSlices->Copy(*(pix->fHPedestalFADCSlices));
251
252
253	return pix;
254	}
255	*/
256
257	// --------------------------------------------------------------------------
258	//
259	// Set bit kSinglePheFitOK from outside
260	//
261	void MHCalibrationChargeBlindPix::SetSinglePheFitOK (const Bool_t b )
262	{
263	b ? SETBIT(fFlags,kSinglePheFitOK) : CLRBIT(fFlags,kSinglePheFitOK);
264	}
265
266	// --------------------------------------------------------------------------
267	//
268	// Set bit kPedestalFitOK from outside
269	//
270	void MHCalibrationChargeBlindPix::SetPedestalFitOK(const Bool_t b)
271	{
272	b ? SETBIT(fFlags,kPedestalFitOK) : CLRBIT(fFlags,kPedestalFitOK);
273	}
274
275	// --------------------------------------------------------------------------
276	//
277	// Ask for status of bit kSinglePheFitOK
278	//
279	const Bool_t MHCalibrationChargeBlindPix::IsSinglePheFitOK() const
280	{
281	return TESTBIT(fFlags,kSinglePheFitOK);
282	}
283
284	// --------------------------------------------------------------------------
285	//
286	// Ask for status of bit kPedestalFitOK
287	//
288	const Bool_t MHCalibrationChargeBlindPix::IsPedestalFitOK() const
289	{
290	return TESTBIT(fFlags,kPedestalFitOK);
291	}
292
293	// --------------------------------------------------------------------------
294	//
295	// Normalizes fASinglePheFADCSlices and fAPedestalFADCSlices
296	//
297	void MHCalibrationChargeBlindPix::FinalizeSinglePheSpectrum()
298	{
299
300	if (fNumSinglePhes > 1)
301	for (Int_t i=0;i<fASinglePheFADCSlices.GetNrows();i++)
302	fASinglePheFADCSlices[i] = fASinglePheFADCSlices[i]/fNumSinglePhes;
303	if (fNumPedestals > 1)
304	for (Int_t i=0;i<fAPedestalFADCSlices.GetNrows();i++)
305	fAPedestalFADCSlices[i] = fAPedestalFADCSlices[i]/fNumPedestals;
306	}
307
308	// --------------------------------------------------------------------------
309	//
310	// Checks again for the size and fills fASinglePheFADCSlices with the FADC slice entries
311	//
312	void MHCalibrationChargeBlindPix::FillSinglePheFADCSlices(const MRawEvtPixelIter &iter)
313	{
314
315	const Int_t n = iter.GetNumHiGainSamples() + iter.GetNumLoGainSamples();
316
317	if (fASinglePheFADCSlices.GetNrows() < n)
318	fASinglePheFADCSlices.ResizeTo(n);
319
320	Int_t i=0;
321
322	Byte_t *start = iter.GetHiGainSamples();
323	Byte_t *end = start + iter.GetNumHiGainSamples();
324
325	for (Byte_t *ptr = start; ptr < end; ptr++, i++)
326	fASinglePheFADCSlices(i) = fASinglePheFADCSlices(i) + (Float_t)*ptr;
327
328	start = iter.GetLoGainSamples();
329	end = start + iter.GetNumLoGainSamples();
330
331	for (Byte_t *ptr = start; ptr < end; ptr++, i++)
332	fASinglePheFADCSlices(i) = fASinglePheFADCSlices(i) + (Float_t)*ptr;
333
334	fNumSinglePhes++;
335	}
336
337	// --------------------------------------------------------------------------
338	//
339	// Checks again for the size and fills fAPedestalFADCSlices with the FADC slice entries
340	//
341	void MHCalibrationChargeBlindPix::FillPedestalFADCSlices(const MRawEvtPixelIter &iter)
342	{
343
344	const Int_t n = iter.GetNumHiGainSamples() + iter.GetNumLoGainSamples();
345
346	if (fAPedestalFADCSlices.GetNrows() < n)
347	fAPedestalFADCSlices.ResizeTo(n);
348
349	Int_t i = 0;
350	Byte_t *start = iter.GetHiGainSamples();
351	Byte_t *end = start + iter.GetNumHiGainSamples();
352
353	for (Byte_t *ptr = start; ptr < end; ptr++, i++)
354	fAPedestalFADCSlices(i) = fAPedestalFADCSlices(i)+ (Float_t)*ptr;
355
356	start = iter.GetLoGainSamples();
357	end = start + iter.GetNumLoGainSamples();
358
359	for (Byte_t *ptr = start; ptr < end; ptr++, i++)
360	fAPedestalFADCSlices(i) = fAPedestalFADCSlices(i)+ (Float_t)*ptr;
361
362	fNumPedestals++;
363	}
364
365
366	// --------------------------------------------------------------------------
367	//
368	// Task to simulate single phe spectrum with the given parameters
369	//
370	Bool_t MHCalibrationChargeBlindPix::SimulateSinglePhe(Double_t lambda, Double_t mu0, Double_t mu1, Double_t sigma0, Double_t sigma1)
371	{
372
373	gRandom->SetSeed();
374
375	if (fHGausHist.GetIntegral() != 0)
376	{
377	*fLog << err << "Histogram " << fHGausHist.GetTitle() << " is already filled. " << endl;
378	*fLog << err << "Create new class MHCalibrationBlindPixel for simulation! " << endl;
379	return kFALSE;
380	}
381
382	if (!InitFit())
383	return kFALSE;
384
385	for (Int_t i=0;i<10000; i++)
386	fHGausHist.Fill(fSinglePheFit->GetRandom());
387
388	return kTRUE;
389	}
390
391	// --------------------------------------------------------------------------
392	//
393	// - Get the ranges from the stripped histogram
394	// - choose reasonable start values for the fit
395	// - initialize the fit function depending on fFitFunc
396	// - initialize parameter names and limits depending on fFitFunc
397	//
398	Bool_t MHCalibrationChargeBlindPix::InitFit()
399	{
400
401	//
402	// Get the fitting ranges
403	//
404	Axis_t rmin = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetFirst());
405	Axis_t rmax = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetLast());
406
407	if (rmin < 0.)
408	rmin = 0.;
409
410	//
411	// First guesses for the fit (should be as close to reality as possible,
412	// otherwise the fit goes gaga because of high number of dimensions ...
413	//
414	const Stat_t entries = fHGausHist.Integral("width");
415	const Double_t lambda_guess = 0.5;
416	//const Double_t maximum_bin = fHGausHist.GetBinCenter(fHGausHist.GetMaximumBin());
417	const Double_t norm = entries/TMath::Sqrt(TMath::TwoPi());
418
419	//
420	// Initialize the fit function
421	//
422	switch (fFitFunc)
423	{
424	case kEPoisson4:
425	fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto4,rmin,rmax,6);
426	rmin += 6.5;
427	break;
428	case kEPoisson5:
429	fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto5,rmin,rmax,6);
430	rmin = 0.;
431	break;
432	case kEPoisson6:
433	fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto6,rmin,rmax,6);
434	break;
435	case kEPolya:
436	fSinglePheFit = new TF1("SinglePheFit",&fPolya,rmin,rmax,8);
437	break;
438	case kEMichele:
439	fSinglePheFit = new TF1("SinglePheFit",&fFitFuncMichele,rmin,rmax,9);
440	break;
441	default:
442	*fLog << warn << "WARNING: Could not find Fit Function for Blind Pixel " << endl;
443	return kFALSE;
444	break;
445	}
446
447	if (!fSinglePheFit)
448	{
449	*fLog << warn << dbginf << "WARNING: Could not create fit function for Single Phe fit" << endl;
450	return kFALSE;
451	}
452
453	//
454	// For the fits, we have to take special care since ROOT
455	// has stored the function pointer in a global list which
456	// lead to removing the object twice. We have to take out
457	// the following functions of the global list of functions
458	// as well:
459	//
460	gROOT->GetListOfFunctions()->Remove(fSinglePheFit);
461
462	const Double_t mu_0_guess = 13.5;
463	const Double_t si_0_guess = 2.5;
464	const Double_t mu_1_guess = 30.;
465	const Double_t si_1_guess = si_0_guess + si_0_guess;
466	// Michele
467	const Double_t lambda_1cat_guess = 1.00;
468	const Double_t lambda_1dyn_guess = lambda_1cat_guess/10.;
469	const Double_t mu_1cat_guess = 50.;
470	const Double_t mu_1dyn_guess = 17.;
471	const Double_t si_1cat_guess = si_0_guess + si_0_guess;
472	const Double_t si_1dyn_guess = si_0_guess + si_0_guess/2.;
473	// Polya
474	const Double_t excessPoisson_guess = 0.5;
475	const Double_t delta1_guess = 8.;
476	const Double_t delta2_guess = 5.;
477	const Double_t electronicAmp_guess = gkElectronicAmp;
478	const Double_t electronicAmp_limit = gkElectronicAmpErr;
479
480	//
481	// Initialize boundaries and start parameters
482	//
483	switch (fFitFunc)
484	{
485
486	case kEPoisson4:
487	fSinglePheFit->SetParNames( "#lambda", "#mu_{0}", "#mu_{1}", "#sigma_{0}", "#sigma_{1}","Area");
488	fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm);
489	fSinglePheFit->SetParLimits(0,0.,2.);
490	fSinglePheFit->SetParLimits(1,10.,17.);
491	fSinglePheFit->SetParLimits(2,17.,50.);
492	fSinglePheFit->SetParLimits(3,1.,5.);
493	fSinglePheFit->SetParLimits(4,5.,30.);
494	fSinglePheFit->SetParLimits(5,norm-(0.5norm),norm+(0.7norm));
495	break;
496	case kEPoisson5:
497	case kEPoisson6:
498	fSinglePheFit->SetParNames("#lambda","#mu_{0}","#mu_{1}","#sigma_{0}","#sigma_{1}","Area");
499	fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,800.,si_0_guess,500.,norm);
500	fSinglePheFit->SetParLimits(0,0.,2.);
501	fSinglePheFit->SetParLimits(1,0.,100.);
502	fSinglePheFit->SetParLimits(2,300.,1500.);
503	fSinglePheFit->SetParLimits(3,30.,250.);
504	fSinglePheFit->SetParLimits(4,100.,1000.);
505	fSinglePheFit->SetParLimits(5,norm/1.5,norm*1.5);
506	break;
507
508	case kEPolya:
509	fSinglePheFit->SetParameters(lambda_guess, excessPoisson_guess,
510	delta1_guess,delta2_guess,
511	electronicAmp_guess,
512	10.,
513	norm,
514	0.);
515	fSinglePheFit->SetParNames("#lambda","b_{tot}",
516	"#delta_{1}","#delta_{2}",
517	"amp_{e}","#sigma_{0}",
518	"Area", "#mu_{0}");
519	fSinglePheFit->SetParLimits(0,0.,1.);
520	fSinglePheFit->SetParLimits(1,0.,1.);
521	fSinglePheFit->SetParLimits(2,6.,12.);
522	fSinglePheFit->SetParLimits(3,3.,8.);
523	fSinglePheFit->SetParLimits(4,electronicAmp_guess-electronicAmp_limit,
524	electronicAmp_guess+electronicAmp_limit);
525	fSinglePheFit->SetParLimits(5,0.,40.);
526	fSinglePheFit->SetParLimits(6,norm-0.1,norm+0.1);
527	fSinglePheFit->SetParLimits(7,-10.,10.);
528	break;
529	case kEMichele:
530	fSinglePheFit->SetParNames("#lambda_{cat}","#lambda_{dyn}",
531	"#mu_{0}","#mu_{1cat}","#mu_{1dyn}",
532	"#sigma_{0}","#sigma_{1cat}","#sigma_{1dyn}",
533	"Area");
534	fSinglePheFit->SetParameters(lambda_1cat_guess, lambda_1dyn_guess,
535	mu_0_guess, mu_1cat_guess,mu_1dyn_guess,
536	si_0_guess, si_1cat_guess,si_1dyn_guess,
537	norm);
538	fSinglePheFit->SetParLimits(0,0.01,2.0);
539	fSinglePheFit->SetParLimits(1,0.,0.5);
540	fSinglePheFit->SetParLimits(2,10.,16.);
541	fSinglePheFit->SetParLimits(3,25.,50.);
542	fSinglePheFit->SetParLimits(4,16.,18.5);
543	fSinglePheFit->SetParLimits(5,1.,5.);
544	fSinglePheFit->SetParLimits(6,10.,50.);
545	fSinglePheFit->SetParLimits(7,5.,10.);
546	fSinglePheFit->SetParLimits(8,norm/2.,norm*2.5);
547	break;
548
549	default:
550	*fLog << warn << "WARNING: Could not find Fit Function for Blind Pixel " << endl;
551	return kFALSE;
552	break;
553	}
554
555	fSinglePheFit->SetRange(rmin,rmax);
556
557	return kTRUE;
558	}
559
560	// --------------------------------------------------------------------------
561	//
562	// - Retrieve the parameters depending on fFitFunc
563	// - Retrieve probability, Chisquare and NDF
564	//
565	void MHCalibrationChargeBlindPix::ExitFit()
566	{
567
568
569	//
570	// Finalize
571	//
572	switch (fFitFunc)
573	{
574
575	case kEPoisson4:
576	case kEPoisson5:
577	case kEPoisson6:
578	case kEPoisson7:
579	fLambda = fSinglePheFit->GetParameter(0);
580	fMu0 = fSinglePheFit->GetParameter(1);
581	fMu1 = fSinglePheFit->GetParameter(2);
582	fSigma0 = fSinglePheFit->GetParameter(3);
583	fSigma1 = fSinglePheFit->GetParameter(4);
584
585	fLambdaErr = fSinglePheFit->GetParError(0);
586	fMu0Err = fSinglePheFit->GetParError(1);
587	fMu1Err = fSinglePheFit->GetParError(2);
588	fSigma0Err = fSinglePheFit->GetParError(3);
589	fSigma1Err = fSinglePheFit->GetParError(4);
590	break;
591	case kEPolya:
592	fLambda = fSinglePheFit->GetParameter(0);
593	fMu0 = fSinglePheFit->GetParameter(7);
594	fMu1 = 0.;
595	fSigma0 = fSinglePheFit->GetParameter(5);
596	fSigma1 = 0.;
597
598	fLambdaErr = fSinglePheFit->GetParError(0);
599	fMu0Err = fSinglePheFit->GetParError(7);
600	fMu1Err = 0.;
601	fSigma0Err = fSinglePheFit->GetParError(5);
602	fSigma1Err = 0.;
603	case kEMichele:
604	fLambda = fSinglePheFit->GetParameter(0);
605	fMu0 = fSinglePheFit->GetParameter(2);
606	fMu1 = fSinglePheFit->GetParameter(3);
607	fSigma0 = fSinglePheFit->GetParameter(5);
608	fSigma1 = fSinglePheFit->GetParameter(6);
609
610	fLambdaErr = fSinglePheFit->GetParError(0);
611	fMu0Err = fSinglePheFit->GetParError(2);
612	fMu1Err = fSinglePheFit->GetParError(3);
613	fSigma0Err = fSinglePheFit->GetParError(5);
614	fSigma1Err = fSinglePheFit->GetParError(6);
615	break;
616	default:
617	break;
618	}
619
620	fProb = fSinglePheFit->GetProb();
621	fChisquare = fSinglePheFit->GetChisquare();
622	fNDF = fSinglePheFit->GetNDF();
623
624	*fLog << all << "Results of the Blind Pixel Fit: " << endl;
625	*fLog << all << "Chisquare: " << fChisquare << endl;
626	*fLog << all << "DoF: " << fNDF << endl;
627	*fLog << all << "Probability: " << fProb << endl;
628
629	}
630
631	// --------------------------------------------------------------------------
632	//
633	// - Executes InitFit()
634	// - Fits the fHGausHist with fSinglePheFit
635	// - Executes ExitFit()
636	//
637	// The fit result is accepted under condition:
638	// 1) The results are not nan's
639	// 2) The NDF is not smaller than fNDFLimit (5)
640	// 3) The Probability is greater than fProbLimit (default 0.001 == 99.9%)
641	// 4) at least fNumSinglePheLimit events are in the single Photo-electron peak
642	//
643	Bool_t MHCalibrationChargeBlindPix::FitSinglePhe(Option_t *opt)
644	{
645
646	if (!InitFit())
647	return kFALSE;
648
649	fHGausHist.Fit(fSinglePheFit,opt);
650
651	ExitFit();
652
653	//
654	// The fit result is accepted under condition:
655	// 1) The results are not nan's
656	// 2) The NDF is not smaller than fNDFLimit (5)
657	// 3) The Probability is greater than fProbLimit (default 0.001 == 99.9%)
658	// 4) at least fNumSinglePheLimit events are in the single Photo-electron peak
659	//
660	if ( TMath::IsNaN(fLambda)
661	\|\| TMath::IsNaN(fLambdaErr)
662	\|\| TMath::IsNaN(fProb)
663	\|\| TMath::IsNaN(fMu0)
664	\|\| TMath::IsNaN(fMu0Err)
665	\|\| TMath::IsNaN(fMu1)
666	\|\| TMath::IsNaN(fMu1Err)
667	\|\| TMath::IsNaN(fSigma0)
668	\|\| TMath::IsNaN(fSigma0Err)
669	\|\| TMath::IsNaN(fSigma1)
670	\|\| TMath::IsNaN(fSigma1Err)
671	\|\| fNDF < fNDFLimit
672	\|\| fProb < fProbLimit )
673	return kFALSE;
674
675	const Stat_t entries = fHGausHist.Integral("width");
676	const Float_t numSinglePhe = TMath::Exp(-1.0fLambda)fLambda*entries;
677
678	if (numSinglePhe < fNumSinglePheLimit)
679	{
680	*fLog << warn << "WARNING - Statistics is too low: Only " << numSinglePhe
681	<< " in the Single Photo-Electron peak " << endl;
682	return kFALSE;
683	}
684	else
685	*fLog << all << numSinglePhe << " in Single Photo-Electron peak " << endl;
686
687	SetSinglePheFitOK();
688	return kTRUE;
689	}
690
691	// --------------------------------------------------------------------------
692	//
693	// - Retrieves limits for the fit
694	// - Fits the fHGausHist with Gauss
695	// - Retrieves the results to fLambdaCheck and fLambdaCheckErr
696	// - Sets a flag IsPedestalFitOK()
697	//
698	void MHCalibrationChargeBlindPix::FitPedestal (Option_t *opt)
699	{
700
701	// Perform the cross-check fitting only the pedestal:
702	const Axis_t rmin = 0.;
703	// const Axis_t rmax = fHGausHist.GetBinCenter(fHGausHist.GetMaximumBin());
704	const Axis_t rmax = fSinglePheCut;
705
706	FitGaus(opt, rmin, rmax);
707
708	const Stat_t entries = fHGausHist.Integral("width");
709	const Double_t pedarea = fFGausFit->Integral(0.,fSinglePheCut);
710
711	fLambdaCheck = TMath::Log(entries/pedarea);
712	// estimate the error by the error of the obtained area from the Gauss-function:
713	fLambdaCheckErr = fFGausFit->GetParError(0)/fFGausFit->GetParameter(0);
714
715	SetPedestalFitOK(IsGausFitOK());
716	return;
717	}
718
719
720	// -------------------------------------------------------------------------
721	//
722	// Draw a legend with the fit results
723	//
724	void MHCalibrationChargeBlindPix::DrawLegend(Option_t *opt)
725	{
726
727	TString option(opt);
728
729	if (!fFitLegend)
730	{
731	fFitLegend = new TPaveText(0.05,0.05,0.95,0.95);
732	fFitLegend->SetLabel(Form("%s%s", "Results of the single PhE Fit (",
733	(fFitFunc == kEPoisson4) ? "Poisson(k=4))" :
734	(fFitFunc == kEPoisson5) ? "Poisson(k=5))" :
735	(fFitFunc == kEPoisson6) ? "Poisson(k=6))" :
736	(fFitFunc == kEPolya ) ? "Polya(k=4))" :
737	(fFitFunc == kEMichele ) ? "Michele)"
738	: " none )" ));
739	fFitLegend->SetTextSize(0.05);
740	}
741	else
742	fFitLegend->Clear();
743
744	const TString line1 =
745	Form("Mean: #lambda = %2.2f #pm %2.2f",fLambda,fLambdaErr);
746	TText *t1 = fFitLegend->AddText(line1.Data());
747	t1->SetBit(kCanDelete);
748
749	const TString line6 =
750	Form("Mean #lambda_{check} = %2.2f #pm %2.2f",fLambdaCheck,fLambdaCheckErr);
751	TText *t2 = fFitLegend->AddText(line6.Data());
752	t2->SetBit(kCanDelete);
753
754	if (option.Contains("datacheck"))
755	{
756	if (fLambda + 3.fLambdaErr < fLambdaCheck - 3.fLambdaCheckErr
757	\|\|
758	fLambda - 3.fLambdaErr > fLambdaCheck + 3.fLambdaCheckErr )
759	{
760	TText *t = fFitLegend->AddText("#lambda and #lambda_{check} more than 3#sigma apart!");
761	t->SetBit(kCanDelete);
762	}
763	}
764	else
765	{
766
767	const TString line2 =
768	Form("Pedestal: #mu_{0} = %2.2f #pm %2.2f",fMu0,fMu0Err);
769	TText *t3 = fFitLegend->AddText(line2.Data());
770	t3->SetBit(kCanDelete);
771
772	const TString line3 =
773	Form("Width Pedestal: #sigma_{0} = %2.2f #pm %2.2f",fSigma0,fSigma0Err);
774	TText *t4 = fFitLegend->AddText(line3.Data());
775	t4->SetBit(kCanDelete);
776
777	const TString line4 =
778	Form("1^{st} Phe-peak: #mu_{1} = %2.2f #pm %2.2f",fMu1,fMu1Err);
779	TText *t5 = fFitLegend->AddText(line4.Data());
780	t5->SetBit(kCanDelete);
781
782	const TString line5 =
783	Form("Width 1^{st} Phe-peak: #sigma_{1} = %2.2f #pm %2.2f",fSigma1,fSigma1Err);
784	TText *t6 = fFitLegend->AddText(line5.Data());
785	t6->SetBit(kCanDelete);
786	}
787
788	const TString line7 =
789	Form("#chi^{2} / N_{dof}: %4.2f / %3i",fChisquare,fNDF);
790	TText *t7 = fFitLegend->AddText(line7.Data());
791	t7->SetBit(kCanDelete);
792
793	const TString line8 =
794	Form("Probability: %6.4f ",fProb);
795	TText *t8 = fFitLegend->AddText(line8.Data());
796	t8->SetBit(kCanDelete);
797
798	if (IsSinglePheFitOK())
799	{
800	TText *t = fFitLegend->AddText("Result of the Fit: OK");
801	t->SetBit(kCanDelete);
802	}
803	else
804	{
805	TText *t = fFitLegend->AddText("Result of the Fit: NOT OK");
806	t->SetBit(kCanDelete);
807	}
808
809	fFitLegend->SetFillColor(IsSinglePheFitOK() ? 80 : 2);
810	fFitLegend->Draw();
811
812	return;
813	}
814
815
816	// -------------------------------------------------------------------------
817	//
818	// Draw the histogram
819	//
820	// The following options can be chosen:
821	//
822	// "": displays the fHGausHist, the fits, the legend and fASinglePheFADCSlices and fAPedestalFADCSlices
823	// "all": executes additionally MHGausEvents::Draw(), with option "fourierevents"
824	// "datacheck" display the fHGausHist, the fits and the legend
825	//
826	void MHCalibrationChargeBlindPix::Draw(Option_t *opt)
827	{
828
829	TString option(opt);
830	option.ToLower();
831
832	Int_t win = 1;
833
834	TVirtualPad *oldpad = gPad ? gPad : MH::MakeDefCanvas(this,900, 600);
835	TVirtualPad *pad = NULL;
836
837	if (option.Contains("all"))
838	{
839	option.ReplaceAll("all","");
840	oldpad->Divide(2,1);
841	win = 2;
842	oldpad->cd(1);
843	TVirtualPad *newpad = gPad;
844	pad = newpad;
845	pad->Divide(2,2);
846	pad->cd(1);
847	}
848	else if (option.Contains("datacheck"))
849	{
850	pad = oldpad;
851	pad->Divide(1,2);
852	pad->cd(1);
853	fHGausHist.SetStats(0);
854	}
855	else
856	{
857	pad = oldpad;
858	pad->Divide(2,2);
859	pad->cd(1);
860	}
861
862	if (!IsEmpty() && !IsOnlyOverflow() && !IsOnlyUnderflow())
863	gPad->SetLogy();
864
865	gPad->SetTicks();
866
867	fHGausHist.Draw();
868	if (fFGausFit )
869	{
870	fFGausFit->SetLineColor(kBlue);
871	fFGausFit->Draw("same");
872	if (!option.Contains("datacheck"))
873	{
874	TLine *line = new TLine(fSinglePheCut, 0., fSinglePheCut, 10.);
875	line->SetBit(kCanDelete);
876	line->SetLineColor(kBlue);
877	line->SetLineWidth(3);
878	line->DrawLine(fSinglePheCut, 0., fSinglePheCut, 2.);
879	}
880	}
881
882	if (fSinglePheFit)
883	{
884	fSinglePheFit->SetFillStyle(0);
885	fSinglePheFit->SetLineWidth(3);
886	fSinglePheFit->SetLineColor(IsSinglePheFitOK() ? kGreen : kRed);
887	fSinglePheFit->Draw("same");
888	}
889
890	pad->cd(2);
891	DrawLegend(option.Data());
892
893	if (option.Contains("datacheck"))
894	return;
895
896	pad->cd(3);
897
898	if (fASinglePheFADCSlices.GetNrows()!=1)
899	{
900	if (fHSinglePheFADCSlices)
901	delete fHSinglePheFADCSlices;
902
903	fHSinglePheFADCSlices = new TH1F(fASinglePheFADCSlices);
904	fHSinglePheFADCSlices->SetName("SinglePheFADCSlices");
905	fHSinglePheFADCSlices->SetTitle(Form("%s%4.1f","Assumed Single Phe FADC Slices, Sum > ",fSinglePheCut));
906	fHSinglePheFADCSlices->SetXTitle("FADC slice number");
907	fHSinglePheFADCSlices->SetYTitle("FADC counts");
908	const Int_t nbins = fHSinglePheFADCSlices->GetNbinsX();
909	TH2D *nulls = new TH2D("Nulls",fHSinglePheFADCSlices->GetTitle(),nbins,0.,
910	fHSinglePheFADCSlices->GetXaxis()->GetBinCenter(nbins),
911	100,0.,50.);
912	nulls->SetDirectory(NULL);
913	nulls->SetBit(kCanDelete);
914	nulls->GetXaxis()->SetTitle(fHSinglePheFADCSlices->GetXaxis()->GetTitle());
915	nulls->GetYaxis()->SetTitle(fHSinglePheFADCSlices->GetYaxis()->GetTitle());
916	nulls->GetXaxis()->CenterTitle();
917	nulls->GetYaxis()->CenterTitle();
918	nulls->SetStats(0);
919	nulls->Draw();
920	fHSinglePheFADCSlices->Draw("same");
921	}
922
923	pad->cd(4);
924	if (fAPedestalFADCSlices.GetNrows()!=1)
925	{
926
927	if (fHPedestalFADCSlices)
928	delete fHPedestalFADCSlices;
929
930	fHPedestalFADCSlices = new TH1F(fAPedestalFADCSlices);
931	fHPedestalFADCSlices->SetName("PedestalFADCSlices");
932	fHPedestalFADCSlices->SetTitle(Form("%s%4.1f","Pedestal FADC Slices, Sum < ",fSinglePheCut));
933	fHPedestalFADCSlices->SetXTitle("FADC slice number");
934	fHPedestalFADCSlices->SetYTitle("FADC counts");
935	const Int_t nbins = fHPedestalFADCSlices->GetNbinsX();
936	TH2D *nullp = new TH2D("Nullp",fHPedestalFADCSlices->GetTitle(),nbins,0.,
937	fHPedestalFADCSlices->GetXaxis()->GetBinCenter(nbins),
938	100,0.,50.);
939	nullp->SetDirectory(NULL);
940	nullp->SetBit(kCanDelete);
941	nullp->GetXaxis()->SetTitle(fHPedestalFADCSlices->GetXaxis()->GetTitle());
942	nullp->GetYaxis()->SetTitle(fHPedestalFADCSlices->GetYaxis()->GetTitle());
943	nullp->GetXaxis()->CenterTitle();
944	nullp->GetYaxis()->CenterTitle();
945	nullp->SetStats(0);
946	nullp->Draw();
947	fHPedestalFADCSlices->Draw("same");
948	}
949
950	if (win < 2)
951	return;
952
953	oldpad->cd(2);
954	MHCalibrationPix::Draw("fourierevents");
955	}
956
957
958
959
960
961
962
963
964
965
966

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

Download in other formats: