Context Navigation

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

Visit:

Last change on this file since 4955 was 4943, checked in by gaug, 21 years ago
* empty log message *
File size: 37.3 KB

Line
1	/* ======================================================================== *\
2	!
3	! *
4	! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5	! * Software. It is distributed to you in the hope that it can be a useful
6	! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7	! * It is distributed WITHOUT ANY WARRANTY.
8	! *
9	! * Permission to use, copy, modify and distribute this software and its
10	! * documentation for any purpose is hereby granted without fee,
11	! * provided that the above copyright notice appear in all copies and
12	! * that both that copyright notice and this permission notice appear
13	! * in supporting documentation. It is provided "as is" without express
14	! * or implied warranty.
15	! *
16	!
17	!
18	! Author(s): Markus Gaug 02/2004 <mailto:markus@ifae.es>
19	!
20	! Copyright: MAGIC Software Development, 2000-2004
21	!
22	!
23	\* ======================================================================== */
24
25	//////////////////////////////////////////////////////////////////////////////
26	//
27	// MHCalibrationChargeBlindPix
28	//
29	// Histogram class for the charge calibration of the Blind Pixel.
30	// Stores and fits the charges and stores the averaged assumed pedestal and
31	// single-phe FADC slice entries. Charges are taken from MExtractedSignalBlindPix.
32	// Performs the Single Photo-electron fit to extract the Poisson mean and its errors.
33	//
34	// Different fits can be chosen with the function 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	MHCalibrationPix::Clear();
263	return;
264	}
265
266	// --------------------------------------------------------------------------
267	//
268	// Empty function to overload MHCalibrationChargePix::Reset()
269	//
270	void MHCalibrationChargeBlindPix::Reset()
271	{
272	}
273
274
275	// --------------------------------------------------------------------------
276	//
277	// Use the MHGausEvents::Clone function and clone additionally the rest of the
278	// data members.
279	//
280	TObject MHCalibrationChargeBlindPix::Clone(const char name) const
281	{
282
283	MHCalibrationChargeBlindPix &pix = (MHCalibrationChargeBlindPix&)*MHCalibrationPix::Clone(name);
284	//
285	// Copy data members
286	//
287	pix.fSinglePheCut = fSinglePheCut;
288	pix.fNumSinglePheLimit = fNumSinglePheLimit;
289
290	pix.fASinglePheFADCSlices = fASinglePheFADCSlices;
291	pix.fAPedestalFADCSlices = fAPedestalFADCSlices;
292
293	pix.fSinglePheFit = (TF1*)fSinglePheFit->Clone();
294
295	pix.fNumSinglePhes = fNumSinglePhes;
296	pix.fNumPedestals = fNumPedestals;
297	pix.fLambda = fLambda;
298	pix.fLambdaCheck = fLambdaCheck;
299	pix.fMu0 = fMu0;
300	pix.fMu1 = fMu1;
301	pix.fSigma0 = fSigma0;
302	pix.fSigma1 = fSigma1;
303	pix.fLambdaErr = fLambdaErr;
304	pix.fLambdaCheckErr = fLambdaCheckErr;
305	pix.fMu0Err = fMu0Err;
306	pix.fMu1Err = fMu1Err;
307	pix.fSigma0Err = fSigma0Err;
308	pix.fSigma1Err = fSigma1Err;
309	pix.fChisquare = fChisquare;
310	pix.fNDF = fNDF;
311	pix.fProb = fProb;
312	pix.fMeanPedestal = fMeanPedestal;
313	pix.fSigmaPedestal = fSigmaPedestal;
314	pix.fMeanPedestalErr = fMeanPedestalErr;
315	pix.fSigmaPedestalErr = fSigmaPedestalErr;
316	pix.fFlags = fFlags;
317
318	return &pix;
319	}
320
321	/*
322	// --------------------------------------------------------------------------
323	//
324	// Our own clone function is necessary since root 3.01/06 or Mars 0.4
325	// I don't know the reason.
326	//
327	// Creates new MHCalibrationCam
328	//
329	TObject MHCalibrationChargeBlindPix::Clone(const char ) const
330	{
331
332	MHCalibrationChargeBlindPix *pix = new MHCalibrationChargeBlindPix();
333	this->Copy(*pix);
334
335	this->fHGausHist.Copy(pix->fHGausHist);
336	this->fSinglePheFit->Copy(*(pix->fSinglePheFit));
337	this->fHSinglePheFADCSlices->Copy(*(pix->fHSinglePheFADCSlices));
338	this->fHPedestalFADCSlices->Copy(*(pix->fHPedestalFADCSlices));
339
340
341	return pix;
342	}
343	*/
344
345	// --------------------------------------------------------------------------
346	//
347	// Set bit kSinglePheFitOK from outside
348	//
349	void MHCalibrationChargeBlindPix::SetSinglePheFitOK (const Bool_t b )
350	{
351	b ? SETBIT(fFlags,kSinglePheFitOK) : CLRBIT(fFlags,kSinglePheFitOK);
352	}
353
354	// --------------------------------------------------------------------------
355	//
356	// Set bit kPedestalFitOK from outside
357	//
358	void MHCalibrationChargeBlindPix::SetPedestalFitOK(const Bool_t b)
359	{
360	b ? SETBIT(fFlags,kPedestalFitOK) : CLRBIT(fFlags,kPedestalFitOK);
361	}
362
363	// --------------------------------------------------------------------------
364	//
365	// Ask for status of bit kSinglePheFitOK
366	//
367	const Bool_t MHCalibrationChargeBlindPix::IsSinglePheFitOK() const
368	{
369	return TESTBIT(fFlags,kSinglePheFitOK);
370	}
371
372	// --------------------------------------------------------------------------
373	//
374	// Ask for status of bit kPedestalFitOK
375	//
376	const Bool_t MHCalibrationChargeBlindPix::IsPedestalFitOK() const
377	{
378	return TESTBIT(fFlags,kPedestalFitOK);
379	}
380
381	// --------------------------------------------------------------------------
382	//
383	// Gets the pointers to:
384	// - MRawEvtData
385	// - MExtractedSignalBlindPixel
386	//
387	Bool_t MHCalibrationChargeBlindPix::SetupFill(const MParList *pList)
388	{
389
390	fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
391	if (!fRawEvt)
392	{
393	*fLog << err << "MRawEvtData not found... aborting." << endl;
394	return kFALSE;
395	}
396
397	fSignal = (MExtractedSignalBlindPixel*)pList->FindObject("MExtractedSignalBlindPixel");
398	if (!fSignal)
399	{
400	*fLog << err << "MExtractedSignalBlindPixel not found... aborting " << endl;
401	return kFALSE;
402	}
403
404	return kTRUE;
405	}
406
407	// --------------------------------------------------------------------------
408	//
409	// Gets or creates the pointers to:
410	// - MCalibrationChargeBlindPix
411	//
412	// Calls:
413	// - MHGausHist::InitBins()
414	//
415	Bool_t MHCalibrationChargeBlindPix::ReInit(MParList *pList)
416	{
417
418	fBlindPix = (MCalibrationChargeBlindPix*)pList->FindCreateObj("MCalibrationChargeBlindPix");
419	if (!fBlindPix)
420	return kFALSE;
421
422	const Int_t samples = fSignal->GetNumFADCSamples();
423	const Int_t integ = fSignal->IsExtractionType( MExtractBlindPixel::kIntegral );
424
425	//
426	// Modify the histogram size in case, integrals have been used
427	//
428	if ( fLast < samplesintegfgChargeLast )
429	{
430	SetLast ( samples * (fgChargeLast+0.5) - 0.5 );
431	SetSinglePheCut( samples * fgSinglePheCut );
432	}
433
434	MHCalibrationPix::InitBins();
435
436	return kTRUE;
437	}
438
439	// --------------------------------------------------------------------------
440	//
441	// Retrieves from MExtractedSignalBlindPixel:
442	// - number of FADC samples
443	// - extracted signal
444	// - blind Pixel ID
445	//
446	// Resizes (if necessary):
447	// - fASinglePheFADCSlices to sum of HiGain and LoGain samples
448	// - fAPedestalFADCSlices to sum of HiGain and LoGain samples
449	//
450	// Fills the following histograms:
451	// - MHGausEvents::FillHistAndArray(signal)
452	//
453	// Creates MRawEvtPixelIter, jumps to blind pixel ID,
454	// fills the vectors fASinglePheFADCSlices and fAPedestalFADCSlices
455	// with the full FADC slices, depending on the size of the signal w.r.t. fSinglePheCut
456	//
457	Bool_t MHCalibrationChargeBlindPix::Fill(const MParContainer *par, const Stat_t w)
458	{
459
460	const Int_t samples = (Int_t)fRawEvt->GetNumHiGainSamples()
461	+(Int_t)fRawEvt->GetNumLoGainSamples();
462
463	if (!fASinglePheFADCSlices.IsValid())
464	{
465	fASinglePheFADCSlices.ResizeTo(samples);
466	fAPedestalFADCSlices.ResizeTo(samples);
467	}
468
469	if (fASinglePheFADCSlices.GetNrows() != samples)
470	{
471	fASinglePheFADCSlices.ResizeTo(samples);
472	fAPedestalFADCSlices.ResizeTo(samples);
473	}
474
475	Float_t slices = (Float_t)fSignal->GetNumFADCSamples();
476
477	if (slices == 0.)
478	{
479	*fLog << err
480	<< "Number of used signal slices in MExtractedSignalBlindPix "
481	<< "is zero ... abort."
482	<< endl;
483	return kFALSE;
484	}
485
486	//
487	// Signal extraction and histogram filling
488	//
489	const Float_t signal = fSignal->GetExtractedSignal(fPixId);
490
491	if (signal > -0.5)
492	FillHist(signal);
493	else
494	return kTRUE;
495
496	//
497	// In order to study the single-phe posistion, we extract the slices
498	//
499	const Int_t blindpixIdx = fSignal->GetBlindPixelIdx(fPixId);
500
501	MRawEvtPixelIter pixel(fRawEvt);
502	pixel.Jump(blindpixIdx);
503
504	if (signal > fSinglePheCut)
505	FillSinglePheFADCSlices(pixel);
506	else
507	FillPedestalFADCSlices(pixel);
508
509	return kTRUE;
510	}
511
512	// --------------------------------------------------------------------------
513	//
514	// Returns kFALSE, if empty
515	//
516	// - Creates the fourier spectrum and sets bit MHGausEvents::IsFourierSpectrumOK()
517	// - Retrieves the pedestals from MExtractedSignalBlindPixel
518	// - Normalizes fASinglePheFADCSlices and fAPedestalFADCSlices
519	// - Executes FitPedestal()
520	// - Executes FitSinglePhe()
521	// - Retrieves fit results and stores them in MCalibrationChargeBlindPix
522	//
523	Bool_t MHCalibrationChargeBlindPix::Finalize()
524	{
525
526	if (IsEmpty())
527	{
528	*fLog << err << GetDescriptor() << " ID: " << fPixId
529	<< " My histogram has not been filled !! " << endl;
530	return kFALSE;
531	}
532
533	fMeanPedestal = fSignal->GetPed();
534	fMeanPedestalErr = fSignal->GetPedErr();
535	fSigmaPedestal = fSignal->GetPedRms();
536	fSigmaPedestalErr = fSignal->GetPedRmsErr();
537
538	if (fNumSinglePhes > 1)
539	for (Int_t i=0;i<fASinglePheFADCSlices.GetNrows();i++)
540	fASinglePheFADCSlices[i] = fASinglePheFADCSlices[i]/fNumSinglePhes;
541	if (fNumPedestals > 1)
542	for (Int_t i=0;i<fAPedestalFADCSlices.GetNrows();i++)
543	fAPedestalFADCSlices[i] = fAPedestalFADCSlices[i]/fNumPedestals;
544
545	FitPedestal();
546
547	fBlindPix->SetValid(kTRUE);
548
549	if (FitSinglePhe())
550	fBlindPix->SetSinglePheFitOK();
551	else
552	fBlindPix->SetValid(IsPedestalFitOK());
553
554	fBlindPix->SetLambda ( fLambda );
555	fBlindPix->SetLambdaVar ( fLambdaErr*fLambdaErr );
556	fBlindPix->SetMu0 ( fMu0 );
557	fBlindPix->SetMu0Err ( fMu0Err );
558	fBlindPix->SetMu1 ( fMu1 );
559	fBlindPix->SetMu1Err ( fMu1Err );
560	fBlindPix->SetSigma0 ( fSigma0 );
561	fBlindPix->SetSigma0Err ( fSigma0Err );
562	fBlindPix->SetSigma1 ( fSigma1 );
563	fBlindPix->SetSigma1Err ( fSigma1Err );
564	fBlindPix->SetProb ( fProb );
565
566	fBlindPix->SetLambdaCheck ( fLambdaCheck );
567	fBlindPix->SetLambdaCheckErr ( fLambdaCheckErr );
568
569	return kTRUE;
570	}
571
572
573	// --------------------------------------------------------------------------
574	//
575	// Checks again for the size and fills fASinglePheFADCSlices with the FADC slice entries
576	//
577	void MHCalibrationChargeBlindPix::FillSinglePheFADCSlices(const MRawEvtPixelIter &iter)
578	{
579
580	const Int_t n = iter.GetNumHiGainSamples() + iter.GetNumLoGainSamples();
581
582	if (fASinglePheFADCSlices.GetNrows() < n)
583	fASinglePheFADCSlices.ResizeTo(n);
584
585	Int_t i=0;
586
587	Byte_t *start = iter.GetHiGainSamples();
588	Byte_t *end = start + iter.GetNumHiGainSamples();
589
590	for (Byte_t *ptr = start; ptr < end; ptr++, i++)
591	fASinglePheFADCSlices(i) = fASinglePheFADCSlices(i) + (Float_t)*ptr;
592
593	start = iter.GetLoGainSamples();
594	end = start + iter.GetNumLoGainSamples();
595
596	for (Byte_t *ptr = start; ptr < end; ptr++, i++)
597	fASinglePheFADCSlices(i) = fASinglePheFADCSlices(i) + (Float_t)*ptr;
598
599	fNumSinglePhes++;
600	}
601
602	// --------------------------------------------------------------------------
603	//
604	// Checks again for the size and fills fAPedestalFADCSlices with the FADC slice entries
605	//
606	void MHCalibrationChargeBlindPix::FillPedestalFADCSlices(const MRawEvtPixelIter &iter)
607	{
608
609	const Int_t n = iter.GetNumHiGainSamples() + iter.GetNumLoGainSamples();
610
611	if (fAPedestalFADCSlices.GetNrows() < n)
612	fAPedestalFADCSlices.ResizeTo(n);
613
614	Int_t i = 0;
615	Byte_t *start = iter.GetHiGainSamples();
616	Byte_t *end = start + iter.GetNumHiGainSamples();
617
618	for (Byte_t *ptr = start; ptr < end; ptr++, i++)
619	fAPedestalFADCSlices(i) = fAPedestalFADCSlices(i)+ (Float_t)*ptr;
620
621	start = iter.GetLoGainSamples();
622	end = start + iter.GetNumLoGainSamples();
623
624	for (Byte_t *ptr = start; ptr < end; ptr++, i++)
625	fAPedestalFADCSlices(i) = fAPedestalFADCSlices(i)+ (Float_t)*ptr;
626
627	fNumPedestals++;
628	}
629
630
631	// --------------------------------------------------------------------------
632	//
633	// Task to simulate single phe spectrum with the given parameters
634	//
635	Bool_t MHCalibrationChargeBlindPix::SimulateSinglePhe(Double_t lambda, Double_t mu0, Double_t mu1, Double_t sigma0, Double_t sigma1)
636	{
637
638	gRandom->SetSeed();
639
640	if (fHGausHist.GetIntegral() != 0)
641	{
642	*fLog << err << "Histogram " << fHGausHist.GetTitle() << " is already filled. " << endl;
643	*fLog << err << "Create new class MHCalibrationBlindPixel for simulation! " << endl;
644	return kFALSE;
645	}
646
647	if (!InitFit())
648	return kFALSE;
649
650	for (Int_t i=0;i<10000; i++)
651	fHGausHist.Fill(fSinglePheFit->GetRandom());
652
653	return kTRUE;
654	}
655
656	// --------------------------------------------------------------------------
657	//
658	// - Get the ranges from the stripped histogram
659	// - choose reasonable start values for the fit
660	// - initialize the fit function depending on fFitFunc
661	// - initialize parameter names and limits depending on fFitFunc
662	//
663	Bool_t MHCalibrationChargeBlindPix::InitFit()
664	{
665
666	//
667	// Get the fitting ranges
668	//
669	Axis_t rmin = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetFirst());
670	Axis_t rmax = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetLast());
671
672	if (rmin < 0.)
673	rmin = 0.;
674
675	//
676	// First guesses for the fit (should be as close to reality as possible,
677	// otherwise the fit goes gaga because of high number of dimensions ...
678	//
679	const Stat_t entries = fHGausHist.Integral("width");
680	const Double_t lambda_guess = 0.5;
681	//const Double_t maximum_bin = fHGausHist.GetBinCenter(fHGausHist.GetMaximumBin());
682	const Double_t norm = entries/TMath::Sqrt(TMath::TwoPi());
683
684	//
685	// Initialize the fit function
686	//
687	switch (fFitFunc)
688	{
689	case kEPoisson4:
690	fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto4,rmin,rmax,6);
691	rmin += 6.5;
692	break;
693	case kEPoisson5:
694	fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto5,rmin,rmax,6);
695	rmin = 0.;
696	break;
697	case kEPoisson6:
698	fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto6,rmin,rmax,6);
699	break;
700	case kEPolya:
701	fSinglePheFit = new TF1("SinglePheFit",&fPolya,rmin,rmax,8);
702	break;
703	case kEMichele:
704	fSinglePheFit = new TF1("SinglePheFit",&fFitFuncMichele,rmin,rmax,9);
705	break;
706	default:
707	*fLog << warn << "WARNING: Could not find Fit Function for Blind Pixel " << endl;
708	return kFALSE;
709	break;
710	}
711
712	if (!fSinglePheFit)
713	{
714	*fLog << warn << dbginf << "WARNING: Could not create fit function for Single Phe fit" << endl;
715	return kFALSE;
716	}
717
718	const Double_t mu_0_guess = 13.5;
719	const Double_t si_0_guess = 2.5;
720	const Double_t mu_1_guess = 30.;
721	const Double_t si_1_guess = si_0_guess + si_0_guess;
722	// Michele
723	const Double_t lambda_1cat_guess = 1.00;
724	const Double_t lambda_1dyn_guess = lambda_1cat_guess/10.;
725	const Double_t mu_1cat_guess = 50.;
726	const Double_t mu_1dyn_guess = 17.;
727	const Double_t si_1cat_guess = si_0_guess + si_0_guess;
728	const Double_t si_1dyn_guess = si_0_guess + si_0_guess/2.;
729	// Polya
730	const Double_t excessPoisson_guess = 0.5;
731	const Double_t delta1_guess = 8.;
732	const Double_t delta2_guess = 5.;
733	const Double_t electronicAmp_guess = gkElectronicAmp;
734	const Double_t electronicAmp_limit = gkElectronicAmpErr;
735
736	//
737	// Initialize boundaries and start parameters
738	//
739	switch (fFitFunc)
740	{
741
742	case kEPoisson4:
743	fSinglePheFit->SetParNames( "#lambda", "#mu_{0}", "#mu_{1}", "#sigma_{0}", "#sigma_{1}","Area");
744	fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm);
745	fSinglePheFit->SetParLimits(0,0.,2.);
746	fSinglePheFit->SetParLimits(1,10.,17.);
747	fSinglePheFit->SetParLimits(2,17.,50.);
748	fSinglePheFit->SetParLimits(3,1.,5.);
749	fSinglePheFit->SetParLimits(4,5.,30.);
750	fSinglePheFit->SetParLimits(5,norm-(0.5norm),norm+(0.7norm));
751	break;
752	case kEPoisson5:
753	case kEPoisson6:
754	fSinglePheFit->SetParNames("#lambda","#mu_{0}","#mu_{1}","#sigma_{0}","#sigma_{1}","Area");
755	fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,800.,si_0_guess,500.,norm);
756	fSinglePheFit->SetParLimits(0,0.,2.);
757	fSinglePheFit->SetParLimits(1,0.,100.);
758	fSinglePheFit->SetParLimits(2,300.,1500.);
759	fSinglePheFit->SetParLimits(3,30.,250.);
760	fSinglePheFit->SetParLimits(4,100.,1000.);
761	fSinglePheFit->SetParLimits(5,norm/1.5,norm*1.5);
762	break;
763
764	case kEPolya:
765	fSinglePheFit->SetParameters(lambda_guess, excessPoisson_guess,
766	delta1_guess,delta2_guess,
767	electronicAmp_guess,
768	fSigmaPedestal,
769	norm,
770	fMeanPedestal);
771	fSinglePheFit->SetParNames("#lambda","b_{tot}",
772	"#delta_{1}","#delta_{2}",
773	"amp_{e}","#sigma_{0}",
774	"Area", "#mu_{0}");
775	fSinglePheFit->SetParLimits(0,0.,1.);
776	fSinglePheFit->SetParLimits(1,0.,1.);
777	fSinglePheFit->SetParLimits(2,6.,12.);
778	fSinglePheFit->SetParLimits(3,3.,8.);
779	fSinglePheFit->SetParLimits(4,electronicAmp_guess-electronicAmp_limit,
780	electronicAmp_guess+electronicAmp_limit);
781	fSinglePheFit->SetParLimits(5,
782	fSigmaPedestal-3.*fSigmaPedestalErr,
783	fSigmaPedestal+3.*fSigmaPedestalErr);
784	fSinglePheFit->SetParLimits(6,norm-0.1,norm+0.1);
785	fSinglePheFit->SetParLimits(7,
786	fMeanPedestal-3.*fMeanPedestalErr,
787	fMeanPedestal+3.*fMeanPedestalErr);
788	break;
789	case kEMichele:
790	fSinglePheFit->SetParNames("#lambda_{cat}","#lambda_{dyn}",
791	"#mu_{0}","#mu_{1cat}","#mu_{1dyn}",
792	"#sigma_{0}","#sigma_{1cat}","#sigma_{1dyn}",
793	"Area");
794	fSinglePheFit->SetParameters(lambda_1cat_guess, lambda_1dyn_guess,
795	mu_0_guess, mu_1cat_guess,mu_1dyn_guess,
796	si_0_guess, si_1cat_guess,si_1dyn_guess,
797	norm);
798	fSinglePheFit->SetParLimits(0,0.01,2.0);
799	fSinglePheFit->SetParLimits(1,0.,0.5);
800	fSinglePheFit->SetParLimits(2,10.,16.);
801	fSinglePheFit->SetParLimits(3,25.,50.);
802	fSinglePheFit->SetParLimits(4,16.,18.5);
803	fSinglePheFit->SetParLimits(5,1.,5.);
804	fSinglePheFit->SetParLimits(6,10.,50.);
805	fSinglePheFit->SetParLimits(7,5.,10.);
806	fSinglePheFit->SetParLimits(8,norm/2.,norm*2.5);
807	break;
808
809	default:
810	*fLog << warn << "WARNING: Could not find Fit Function for Blind Pixel " << endl;
811	return kFALSE;
812	break;
813	}
814
815	fSinglePheFit->SetRange(rmin,rmax);
816
817	return kTRUE;
818	}
819
820	// --------------------------------------------------------------------------
821	//
822	// - Retrieve the parameters depending on fFitFunc
823	// - Retrieve probability, Chisquare and NDF
824	//
825	void MHCalibrationChargeBlindPix::ExitFit()
826	{
827
828
829	//
830	// Finalize
831	//
832	switch (fFitFunc)
833	{
834
835	case kEPoisson4:
836	case kEPoisson5:
837	case kEPoisson6:
838	case kEPoisson7:
839	fLambda = fSinglePheFit->GetParameter(0);
840	fMu0 = fSinglePheFit->GetParameter(1);
841	fMu1 = fSinglePheFit->GetParameter(2);
842	fSigma0 = fSinglePheFit->GetParameter(3);
843	fSigma1 = fSinglePheFit->GetParameter(4);
844
845	fLambdaErr = fSinglePheFit->GetParError(0);
846	fMu0Err = fSinglePheFit->GetParError(1);
847	fMu1Err = fSinglePheFit->GetParError(2);
848	fSigma0Err = fSinglePheFit->GetParError(3);
849	fSigma1Err = fSinglePheFit->GetParError(4);
850	break;
851	case kEPolya:
852	fLambda = fSinglePheFit->GetParameter(0);
853	fMu0 = fSinglePheFit->GetParameter(7);
854	fMu1 = 0.;
855	fSigma0 = fSinglePheFit->GetParameter(5);
856	fSigma1 = 0.;
857
858	fLambdaErr = fSinglePheFit->GetParError(0);
859	fMu0Err = fSinglePheFit->GetParError(7);
860	fMu1Err = 0.;
861	fSigma0Err = fSinglePheFit->GetParError(5);
862	fSigma1Err = 0.;
863	case kEMichele:
864	fLambda = fSinglePheFit->GetParameter(0);
865	fMu0 = fSinglePheFit->GetParameter(2);
866	fMu1 = fSinglePheFit->GetParameter(3);
867	fSigma0 = fSinglePheFit->GetParameter(5);
868	fSigma1 = fSinglePheFit->GetParameter(6);
869
870	fLambdaErr = fSinglePheFit->GetParError(0);
871	fMu0Err = fSinglePheFit->GetParError(2);
872	fMu1Err = fSinglePheFit->GetParError(3);
873	fSigma0Err = fSinglePheFit->GetParError(5);
874	fSigma1Err = fSinglePheFit->GetParError(6);
875	break;
876	default:
877	break;
878	}
879
880	fProb = fSinglePheFit->GetProb();
881	fChisquare = fSinglePheFit->GetChisquare();
882	fNDF = fSinglePheFit->GetNDF();
883
884	*fLog << all << "Results of the Blind Pixel Fit: " << endl;
885	*fLog << all << "Chisquare: " << fChisquare << endl;
886	*fLog << all << "DoF: " << fNDF << endl;
887	*fLog << all << "Probability: " << fProb << endl;
888
889	}
890
891	// --------------------------------------------------------------------------
892	//
893	// - Executes InitFit()
894	// - Fits the fHGausHist with fSinglePheFit
895	// - Executes ExitFit()
896	//
897	// The fit result is accepted under condition:
898	// 1) The results are not nan's
899	// 2) The NDF is not smaller than fNDFLimit (5)
900	// 3) The Probability is greater than fProbLimit (default 0.001 == 99.9%)
901	// 4) at least fNumSinglePheLimit events are in the single Photo-electron peak
902	//
903	Bool_t MHCalibrationChargeBlindPix::FitSinglePhe(Option_t *opt)
904	{
905
906	if (!InitFit())
907	return kFALSE;
908
909	fHGausHist.Fit(fSinglePheFit,opt);
910
911	ExitFit();
912
913	//
914	// The fit result is accepted under condition:
915	// 1) The results are not nan's
916	// 2) The NDF is not smaller than fNDFLimit (5)
917	// 3) The Probability is greater than fProbLimit (default 0.001 == 99.9%)
918	// 4) at least fNumSinglePheLimit events are in the single Photo-electron peak
919	//
920	if ( TMath::IsNaN(fLambda)
921	\|\| TMath::IsNaN(fLambdaErr)
922	\|\| TMath::IsNaN(fProb)
923	\|\| TMath::IsNaN(fMu0)
924	\|\| TMath::IsNaN(fMu0Err)
925	\|\| TMath::IsNaN(fMu1)
926	\|\| TMath::IsNaN(fMu1Err)
927	\|\| TMath::IsNaN(fSigma0)
928	\|\| TMath::IsNaN(fSigma0Err)
929	\|\| TMath::IsNaN(fSigma1)
930	\|\| TMath::IsNaN(fSigma1Err)
931	\|\| fNDF < fNDFLimit
932	\|\| fProb < fProbLimit )
933	return kFALSE;
934
935	const Stat_t entries = fHGausHist.Integral("width");
936	const Float_t numSinglePhe = TMath::Exp(-1.0fLambda)fLambda*entries;
937
938	if (numSinglePhe < fNumSinglePheLimit)
939	{
940	*fLog << warn << "WARNING - Statistics is too low: Only " << numSinglePhe
941	<< " in the Single Photo-Electron peak " << endl;
942	return kFALSE;
943	}
944	else
945	*fLog << all << numSinglePhe << " in Single Photo-Electron peak " << endl;
946
947	SetSinglePheFitOK();
948	return kTRUE;
949	}
950
951	// --------------------------------------------------------------------------
952	//
953	// - Retrieves limits for the fit
954	// - Fits the fHGausHist with Gauss
955	// - Retrieves the results to fLambdaCheck and fLambdaCheckErr
956	// - Sets a flag IsPedestalFitOK()
957	//
958	void MHCalibrationChargeBlindPix::FitPedestal (Option_t *opt)
959	{
960
961	// Perform the cross-check fitting only the pedestal:
962	const Axis_t rmin = 0.;
963	// const Axis_t rmax = fHGausHist.GetBinCenter(fHGausHist.GetMaximumBin());
964	const Axis_t rmax = fSinglePheCut;
965
966	FitGaus(opt, rmin, rmax);
967
968	const Stat_t entries = fHGausHist.Integral("width");
969	const Double_t pedarea = fFGausFit->Integral(0.,fSinglePheCut);
970
971	fLambdaCheck = TMath::Log(entries/pedarea);
972	// estimate the error by the error of the obtained area from the Gauss-function:
973	fLambdaCheckErr = fFGausFit->GetParError(0)/fFGausFit->GetParameter(0);
974
975	SetPedestalFitOK(IsGausFitOK());
976	return;
977	}
978
979
980	// -------------------------------------------------------------------------
981	//
982	// Draw a legend with the fit results
983	//
984	void MHCalibrationChargeBlindPix::DrawLegend(Option_t *opt)
985	{
986
987	TString option(opt);
988
989	if (!fFitLegend)
990	{
991	fFitLegend = new TPaveText(0.05,0.05,0.95,0.95);
992	fFitLegend->SetLabel(Form("%s%s", "Results of the single PhE Fit (",
993	(fFitFunc == kEPoisson4) ? "Poisson(k=4))" :
994	(fFitFunc == kEPoisson5) ? "Poisson(k=5))" :
995	(fFitFunc == kEPoisson6) ? "Poisson(k=6))" :
996	(fFitFunc == kEPolya ) ? "Polya(k=4))" :
997	(fFitFunc == kEMichele ) ? "Michele)" : " none )" ));
998	fFitLegend->SetTextSize(0.05);
999	}
1000	else
1001	fFitLegend->Clear();
1002
1003	const TString line1 =
1004	Form("Mean: #lambda = %2.2f #pm %2.2f",fLambda,fLambdaErr);
1005	TText *t1 = fFitLegend->AddText(line1.Data());
1006	t1->SetBit(kCanDelete);
1007
1008	const TString line6 =
1009	Form("Mean #lambda_{check} = %2.2f #pm %2.2f",fLambdaCheck,fLambdaCheckErr);
1010	TText *t2 = fFitLegend->AddText(line6.Data());
1011	t2->SetBit(kCanDelete);
1012
1013	if (option.Contains("datacheck"))
1014	{
1015	if (fLambda + 3.fLambdaErr < fLambdaCheck - 3.fLambdaCheckErr
1016	\|\|
1017	fLambda - 3.fLambdaErr > fLambdaCheck + 3.fLambdaCheckErr )
1018	{
1019	TText *t = fFitLegend->AddText("#lambda and #lambda_{check} more than 3#sigma apart!");
1020	t->SetBit(kCanDelete);
1021	}
1022	}
1023	else
1024	{
1025
1026	const TString line2 =
1027	Form("Pedestal: #mu_{0} = %2.2f #pm %2.2f",fMu0,fMu0Err);
1028	TText *t3 = fFitLegend->AddText(line2.Data());
1029	t3->SetBit(kCanDelete);
1030
1031	const TString line3 =
1032	Form("Width Pedestal: #sigma_{0} = %2.2f #pm %2.2f",fSigma0,fSigma0Err);
1033	TText *t4 = fFitLegend->AddText(line3.Data());
1034	t4->SetBit(kCanDelete);
1035
1036	const TString line4 =
1037	Form("1^{st} Phe-peak: #mu_{1} = %2.2f #pm %2.2f",fMu1,fMu1Err);
1038	TText *t5 = fFitLegend->AddText(line4.Data());
1039	t5->SetBit(kCanDelete);
1040
1041	const TString line5 =
1042	Form("Width 1^{st} Phe-peak: #sigma_{1} = %2.2f #pm %2.2f",fSigma1,fSigma1Err);
1043	TText *t6 = fFitLegend->AddText(line5.Data());
1044	t6->SetBit(kCanDelete);
1045	}
1046
1047	const TString line7 =
1048	Form("#chi^{2} / N_{dof}: %4.2f / %3i",fChisquare,fNDF);
1049	TText *t7 = fFitLegend->AddText(line7.Data());
1050	t7->SetBit(kCanDelete);
1051
1052	const TString line8 =
1053	Form("Probability: %6.4f ",fProb);
1054	TText *t8 = fFitLegend->AddText(line8.Data());
1055	t8->SetBit(kCanDelete);
1056
1057	if (IsSinglePheFitOK())
1058	{
1059	TText *t = fFitLegend->AddText("Result of the Fit: OK");
1060	t->SetBit(kCanDelete);
1061	}
1062	else
1063	{
1064	TText *t = fFitLegend->AddText("Result of the Fit: NOT OK");
1065	t->SetBit(kCanDelete);
1066	}
1067
1068	fFitLegend->SetFillColor(IsSinglePheFitOK() ? 80 : 2);
1069	fFitLegend->Draw();
1070
1071	return;
1072	}
1073
1074
1075	// -------------------------------------------------------------------------
1076	//
1077	// Draw the histogram
1078	//
1079	// The following options can be chosen:
1080	//
1081	// "": displays the fHGausHist, the fits, the legend and fASinglePheFADCSlices and fAPedestalFADCSlices
1082	// "all": executes additionally MHGausEvents::Draw(), with option "fourierevents"
1083	// "datacheck" display the fHGausHist, the fits and the legend
1084	//
1085	void MHCalibrationChargeBlindPix::Draw(Option_t *opt)
1086	{
1087
1088	TString option(opt);
1089	option.ToLower();
1090
1091	Int_t win = 1;
1092
1093	TVirtualPad *oldpad = gPad ? gPad : MH::MakeDefCanvas(this,900, 600);
1094	TVirtualPad *pad = NULL;
1095
1096	if (option.Contains("all"))
1097	{
1098	option.ReplaceAll("all","");
1099	oldpad->Divide(2,1);
1100	win = 2;
1101	oldpad->cd(1);
1102	TVirtualPad *newpad = gPad;
1103	pad = newpad;
1104	pad->Divide(2,2);
1105	pad->cd(1);
1106	}
1107	else if (option.Contains("datacheck"))
1108	{
1109	pad = oldpad;
1110	pad->Divide(1,2);
1111	pad->cd(1);
1112	fHGausHist.SetStats(0);
1113	}
1114	else
1115	{
1116	pad = oldpad;
1117	pad->Divide(2,2);
1118	pad->cd(1);
1119	}
1120
1121	if (!IsEmpty() && !IsOnlyOverflow() && !IsOnlyUnderflow())
1122	gPad->SetLogy();
1123
1124	gPad->SetTicks();
1125
1126	fHGausHist.Draw();
1127	if (fFGausFit )
1128	{
1129	fFGausFit->SetLineColor(kBlue);
1130	fFGausFit->Draw("same");
1131	if (!option.Contains("datacheck"))
1132	{
1133	TLine *line = new TLine(fSinglePheCut, 0., fSinglePheCut, 10.);
1134	line->SetBit(kCanDelete);
1135	line->SetLineColor(kBlue);
1136	line->SetLineWidth(3);
1137	line->DrawLine(fSinglePheCut, 0., fSinglePheCut, 2.);
1138	}
1139	}
1140
1141	if (fSinglePheFit)
1142	{
1143	fSinglePheFit->SetFillStyle(0);
1144	fSinglePheFit->SetLineWidth(3);
1145	fSinglePheFit->SetLineColor(IsSinglePheFitOK() ? kGreen : kRed);
1146	fSinglePheFit->Draw("same");
1147	}
1148
1149	pad->cd(2);
1150	DrawLegend(option.Data());
1151
1152	if (option.Contains("datacheck"))
1153	return;
1154
1155	pad->cd(3);
1156
1157	if (fASinglePheFADCSlices.GetNrows()!=1)
1158	{
1159	if (fHSinglePheFADCSlices)
1160	delete fHSinglePheFADCSlices;
1161
1162	fHSinglePheFADCSlices = new TH1F(fASinglePheFADCSlices);
1163	fHSinglePheFADCSlices->SetName("SinglePheFADCSlices");
1164	fHSinglePheFADCSlices->SetTitle(Form("%s%4.1f","Assumed Single Phe FADC Slices, Sum > ",fSinglePheCut));
1165	fHSinglePheFADCSlices->SetXTitle("FADC slice number");
1166	fHSinglePheFADCSlices->SetYTitle("FADC counts");
1167	const Int_t nbins = fHSinglePheFADCSlices->GetNbinsX();
1168	TH2D *nulls = new TH2D("Nulls",fHSinglePheFADCSlices->GetTitle(),nbins,0.,
1169	fHSinglePheFADCSlices->GetXaxis()->GetBinCenter(nbins),
1170	100,0.,50.);
1171	nulls->SetDirectory(NULL);
1172	nulls->SetBit(kCanDelete);
1173	nulls->GetXaxis()->SetTitle(fHSinglePheFADCSlices->GetXaxis()->GetTitle());
1174	nulls->GetYaxis()->SetTitle(fHSinglePheFADCSlices->GetYaxis()->GetTitle());
1175	nulls->GetXaxis()->CenterTitle();
1176	nulls->GetYaxis()->CenterTitle();
1177	nulls->SetStats(0);
1178	nulls->Draw();
1179	fHSinglePheFADCSlices->Draw("same");
1180	}
1181
1182	pad->cd(4);
1183	if (fAPedestalFADCSlices.GetNrows()!=1)
1184	{
1185
1186	if (fHPedestalFADCSlices)
1187	delete fHPedestalFADCSlices;
1188
1189	fHPedestalFADCSlices = new TH1F(fAPedestalFADCSlices);
1190	fHPedestalFADCSlices->SetName("PedestalFADCSlices");
1191	fHPedestalFADCSlices->SetTitle(Form("%s%4.1f","Pedestal FADC Slices, Sum < ",fSinglePheCut));
1192	fHPedestalFADCSlices->SetXTitle("FADC slice number");
1193	fHPedestalFADCSlices->SetYTitle("FADC counts");
1194	const Int_t nbins = fHPedestalFADCSlices->GetNbinsX();
1195	TH2D *nullp = new TH2D("Nullp",fHPedestalFADCSlices->GetTitle(),nbins,0.,
1196	fHPedestalFADCSlices->GetXaxis()->GetBinCenter(nbins),
1197	100,0.,50.);
1198	nullp->SetDirectory(NULL);
1199	nullp->SetBit(kCanDelete);
1200	nullp->GetXaxis()->SetTitle(fHPedestalFADCSlices->GetXaxis()->GetTitle());
1201	nullp->GetYaxis()->SetTitle(fHPedestalFADCSlices->GetYaxis()->GetTitle());
1202	nullp->GetXaxis()->CenterTitle();
1203	nullp->GetYaxis()->CenterTitle();
1204	nullp->SetStats(0);
1205	nullp->Draw();
1206	fHPedestalFADCSlices->Draw("same");
1207	}
1208
1209	if (win < 2)
1210	return;
1211
1212	oldpad->cd(2);
1213	MHCalibrationPix::Draw("fourierevents");
1214	}
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225

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

Download in other formats: