Context Navigation

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

Visit:

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

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

Download in other formats: