Context Navigation

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

Visit:

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

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

Download in other formats: