Context Navigation

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

Visit:

Last change on this file since 5092 was 5084, checked in by gaug, 22 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	// - 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	#if 0
234	// --------------------------------------------------------------------------
235	//
236	// ATTENTION: This nasty Clone function is necessary since the ROOT cloning
237	// lead to crashes on SOME machines (unfortunately not mine...).
238	// This function is a workaround in order to achieve the correct
239	// DrawClone() behaviour.
240	//
241	TObject MHCalibrationChargeBlindPix::Clone(const char name) const
242	{
243
244	MHCalibrationChargeBlindPix &pix =
245	*new MHCalibrationChargeBlindPix(name ? name : fName.Data(),fTitle.Data());
246
247	//
248	// Copy MHGausEvents data members
249	//
250	pix.fBinsAfterStripping = fBinsAfterStripping;
251	pix.fCurrentSize = fCurrentSize;
252	pix.fFlags = fFlags;
253	pix.fPowerProbabilityBins = fPowerProbabilityBins;
254
255	if (fHPowerProbability)
256	pix.fHPowerProbability=(TH1I*)fHPowerProbability->Clone();
257
258	if (fPowerSpectrum)
259	pix.fPowerSpectrum = new TArrayF(*fPowerSpectrum);
260
261	pix.fEvents = fEvents;
262
263	if (fFGausFit)
264	pix.fFGausFit=(TF1*)fFGausFit->Clone();
265	if (fFExpFit)
266	pix.fFExpFit=(TF1*)fFExpFit->Clone();
267
268	pix.fFirst = fFirst;
269
270	if (fGraphEvents)
271	pix.fGraphEvents=(TGraph*)fGraphEvents->Clone();
272	if (fGraphPowerSpectrum)
273	pix.fGraphPowerSpectrum=(TGraph*)fGraphPowerSpectrum->Clone();
274
275	pix.fHGausHist = fHGausHist;
276
277	pix.fLast = fLast;
278	pix.fMean = fMean;
279	pix.fMeanErr = fMeanErr;
280	pix.fNbins = fNbins;
281	pix.fNDFLimit = fNDFLimit;
282	pix.fSigma = fSigma;
283	pix.fSigmaErr = fSigmaErr;
284	pix.fProb = fProb;
285	pix.fProbLimit = fProbLimit;
286
287	//
288	// Copy MHCalibrationPix data members
289	//
290	pix.fEventFrequency = fEventFrequency;
291	pix.fBlackoutLimit = fBlackoutLimit;
292	pix.fSaturated = fSaturated;
293	pix.fPickupLimit = fPickupLimit;
294	pix.fPixId = fPixId;
295
296	//
297	// Copy MHCalibrationChargeBlindPix data members
298	//
299	pix.fSinglePheCut = fSinglePheCut;
300	pix.fNumSinglePheLimit = fNumSinglePheLimit;
301
302	fASinglePheFADCSlices.Copy(pix.fASinglePheFADCSlices);
303	fAPedestalFADCSlices.Copy( pix.fAPedestalFADCSlices);
304
305	pix.fSinglePheFit = (TF1*)fSinglePheFit->Clone();
306
307	pix.fNumSinglePhes = fNumSinglePhes;
308	pix.fNumPedestals = fNumPedestals;
309	pix.fLambda = fLambda;
310	pix.fLambdaCheck = fLambdaCheck;
311	pix.fMu0 = fMu0;
312	pix.fMu1 = fMu1;
313	pix.fSigma0 = fSigma0;
314	pix.fSigma1 = fSigma1;
315	pix.fLambdaErr = fLambdaErr;
316	pix.fLambdaCheckErr = fLambdaCheckErr;
317	pix.fMu0Err = fMu0Err;
318	pix.fMu1Err = fMu1Err;
319	pix.fSigma0Err = fSigma0Err;
320	pix.fSigma1Err = fSigma1Err;
321	pix.fChisquare = fChisquare;
322	pix.fNDF = fNDF;
323	pix.fProb = fProb;
324
325	return &pix;
326	}
327	#endif
328
329	/*
330	// --------------------------------------------------------------------------
331	//
332	// Our own clone function is necessary since root 3.01/06 or Mars 0.4
333	// I don't know the reason.
334	//
335	// Creates new MHCalibrationCam
336	//
337	TObject MHCalibrationChargeBlindPix::Clone(const char ) const
338	{
339
340	MHCalibrationChargeBlindPix *pix = new MHCalibrationChargeBlindPix();
341	this->Copy(*pix);
342
343	this->fHGausHist.Copy(pix->fHGausHist);
344	this->fSinglePheFit->Copy(*(pix->fSinglePheFit));
345	this->fHSinglePheFADCSlices->Copy(*(pix->fHSinglePheFADCSlices));
346	this->fHPedestalFADCSlices->Copy(*(pix->fHPedestalFADCSlices));
347
348
349	return pix;
350	}
351	*/
352
353	// --------------------------------------------------------------------------
354	//
355	// Set bit kSinglePheFitOK from outside
356	//
357	void MHCalibrationChargeBlindPix::SetSinglePheFitOK (const Bool_t b )
358	{
359	b ? SETBIT(fFlags,kSinglePheFitOK) : CLRBIT(fFlags,kSinglePheFitOK);
360	}
361
362	// --------------------------------------------------------------------------
363	//
364	// Set bit kPedestalFitOK from outside
365	//
366	void MHCalibrationChargeBlindPix::SetPedestalFitOK(const Bool_t b)
367	{
368	b ? SETBIT(fFlags,kPedestalFitOK) : CLRBIT(fFlags,kPedestalFitOK);
369	}
370
371	// --------------------------------------------------------------------------
372	//
373	// Ask for status of bit kSinglePheFitOK
374	//
375	const Bool_t MHCalibrationChargeBlindPix::IsSinglePheFitOK() const
376	{
377	return TESTBIT(fFlags,kSinglePheFitOK);
378	}
379
380	// --------------------------------------------------------------------------
381	//
382	// Ask for status of bit kPedestalFitOK
383	//
384	const Bool_t MHCalibrationChargeBlindPix::IsPedestalFitOK() const
385	{
386	return TESTBIT(fFlags,kPedestalFitOK);
387	}
388
389	// --------------------------------------------------------------------------
390	//
391	// Normalizes fASinglePheFADCSlices and fAPedestalFADCSlices
392	//
393	void MHCalibrationChargeBlindPix::FinalizeSinglePheSpectrum()
394	{
395
396	if (fNumSinglePhes > 1)
397	for (Int_t i=0;i<fASinglePheFADCSlices.GetNrows();i++)
398	fASinglePheFADCSlices[i] = fASinglePheFADCSlices[i]/fNumSinglePhes;
399	if (fNumPedestals > 1)
400	for (Int_t i=0;i<fAPedestalFADCSlices.GetNrows();i++)
401	fAPedestalFADCSlices[i] = fAPedestalFADCSlices[i]/fNumPedestals;
402	}
403
404	// --------------------------------------------------------------------------
405	//
406	// Checks again for the size and fills fASinglePheFADCSlices with the FADC slice entries
407	//
408	void MHCalibrationChargeBlindPix::FillSinglePheFADCSlices(const MRawEvtPixelIter &iter)
409	{
410
411	const Int_t n = iter.GetNumHiGainSamples() + iter.GetNumLoGainSamples();
412
413	if (fASinglePheFADCSlices.GetNrows() < n)
414	fASinglePheFADCSlices.ResizeTo(n);
415
416	Int_t i=0;
417
418	Byte_t *start = iter.GetHiGainSamples();
419	Byte_t *end = start + iter.GetNumHiGainSamples();
420
421	for (Byte_t *ptr = start; ptr < end; ptr++, i++)
422	fASinglePheFADCSlices(i) = fASinglePheFADCSlices(i) + (Float_t)*ptr;
423
424	start = iter.GetLoGainSamples();
425	end = start + iter.GetNumLoGainSamples();
426
427	for (Byte_t *ptr = start; ptr < end; ptr++, i++)
428	fASinglePheFADCSlices(i) = fASinglePheFADCSlices(i) + (Float_t)*ptr;
429
430	fNumSinglePhes++;
431	}
432
433	// --------------------------------------------------------------------------
434	//
435	// Checks again for the size and fills fAPedestalFADCSlices with the FADC slice entries
436	//
437	void MHCalibrationChargeBlindPix::FillPedestalFADCSlices(const MRawEvtPixelIter &iter)
438	{
439
440	const Int_t n = iter.GetNumHiGainSamples() + iter.GetNumLoGainSamples();
441
442	if (fAPedestalFADCSlices.GetNrows() < n)
443	fAPedestalFADCSlices.ResizeTo(n);
444
445	Int_t i = 0;
446	Byte_t *start = iter.GetHiGainSamples();
447	Byte_t *end = start + iter.GetNumHiGainSamples();
448
449	for (Byte_t *ptr = start; ptr < end; ptr++, i++)
450	fAPedestalFADCSlices(i) = fAPedestalFADCSlices(i)+ (Float_t)*ptr;
451
452	start = iter.GetLoGainSamples();
453	end = start + iter.GetNumLoGainSamples();
454
455	for (Byte_t *ptr = start; ptr < end; ptr++, i++)
456	fAPedestalFADCSlices(i) = fAPedestalFADCSlices(i)+ (Float_t)*ptr;
457
458	fNumPedestals++;
459	}
460
461
462	// --------------------------------------------------------------------------
463	//
464	// Task to simulate single phe spectrum with the given parameters
465	//
466	Bool_t MHCalibrationChargeBlindPix::SimulateSinglePhe(Double_t lambda, Double_t mu0, Double_t mu1, Double_t sigma0, Double_t sigma1)
467	{
468
469	gRandom->SetSeed();
470
471	if (fHGausHist.GetIntegral() != 0)
472	{
473	*fLog << err << "Histogram " << fHGausHist.GetTitle() << " is already filled. " << endl;
474	*fLog << err << "Create new class MHCalibrationBlindPixel for simulation! " << endl;
475	return kFALSE;
476	}
477
478	if (!InitFit())
479	return kFALSE;
480
481	for (Int_t i=0;i<10000; i++)
482	fHGausHist.Fill(fSinglePheFit->GetRandom());
483
484	return kTRUE;
485	}
486
487	// --------------------------------------------------------------------------
488	//
489	// - Get the ranges from the stripped histogram
490	// - choose reasonable start values for the fit
491	// - initialize the fit function depending on fFitFunc
492	// - initialize parameter names and limits depending on fFitFunc
493	//
494	Bool_t MHCalibrationChargeBlindPix::InitFit()
495	{
496
497	//
498	// Get the fitting ranges
499	//
500	Axis_t rmin = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetFirst());
501	Axis_t rmax = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetLast());
502
503	if (rmin < 0.)
504	rmin = 0.;
505
506	//
507	// First guesses for the fit (should be as close to reality as possible,
508	// otherwise the fit goes gaga because of high number of dimensions ...
509	//
510	const Stat_t entries = fHGausHist.Integral("width");
511	const Double_t lambda_guess = 0.5;
512	//const Double_t maximum_bin = fHGausHist.GetBinCenter(fHGausHist.GetMaximumBin());
513	const Double_t norm = entries/TMath::Sqrt(TMath::TwoPi());
514
515	//
516	// Initialize the fit function
517	//
518	switch (fFitFunc)
519	{
520	case kEPoisson4:
521	fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto4,rmin,rmax,6);
522	rmin += 6.5;
523	break;
524	case kEPoisson5:
525	fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto5,rmin,rmax,6);
526	rmin = 0.;
527	break;
528	case kEPoisson6:
529	fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto6,rmin,rmax,6);
530	break;
531	case kEPolya:
532	fSinglePheFit = new TF1("SinglePheFit",&fPolya,rmin,rmax,8);
533	break;
534	case kEMichele:
535	fSinglePheFit = new TF1("SinglePheFit",&fFitFuncMichele,rmin,rmax,9);
536	break;
537	default:
538	*fLog << warn << "WARNING: Could not find Fit Function for Blind Pixel " << endl;
539	return kFALSE;
540	break;
541	}
542
543	if (!fSinglePheFit)
544	{
545	*fLog << warn << dbginf << "WARNING: Could not create fit function for Single Phe fit" << endl;
546	return kFALSE;
547	}
548
549	//
550	// For the fits, we have to take special care since ROOT
551	// has stored the function pointer in a global list which
552	// lead to removing the object twice. We have to take out
553	// the following functions of the global list of functions
554	// as well:
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: