Context Navigation

source: trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeBlindPix.cc@ 4204

Visit:

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

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

Download in other formats: