Context Navigation

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

Visit:

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

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

Download in other formats: