Context Navigation

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

Visit:

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

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

Download in other formats: