Context Navigation

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

Visit:

Last change on this file since 4382 was 4380, checked in by gaug, 22 years ago
* empty log message *
File size: 34.4 KB

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

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

Download in other formats: