Context Navigation

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

Visit:

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

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

Download in other formats: