Context Navigation

source: trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.cc@ 3114

Visit:

Last change on this file since 3114 was 3112, checked in by tbretz, 21 years ago
* empty log message *
File size: 26.9 KB

Line
1	/* ======================================================================== *\
2	!
3	! *
4	! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5	! * Software. It is distributed to you in the hope that it can be a useful
6	! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7	! * It is distributed WITHOUT ANY WARRANTY.
8	! *
9	! * Permission to use, copy, modify and distribute this software and its
10	! * documentation for any purpose is hereby granted without fee,
11	! * provided that the above copyright notice appear in all copies and
12	! * that both that copyright notice and this permission notice appear
13	! * in supporting documentation. It is provided "as is" without express
14	! * or implied warranty.
15	! *
16	!
17	!
18	! Author(s): Markus Gaug 11/2003 <mailto:markus@ifae.es>
19	!
20	! Copyright: MAGIC Software Development, 2000-2002
21	!
22	!
23	\* ======================================================================== */
24
25	//////////////////////////////////////////////////////////////////////////////
26	//
27	// MHCalibrationBlindPixel
28	//
29	// Performs all the Single Photo-Electron Fit to extract
30	// the mean number of photons and to derive the light flux
31	//
32	// The fit result is accepted under condition that:
33	// 1) the Probability is greater than gkProbLimit (default 0.001 == 99.7%)
34	// 2) at least 100 events are in the single Photo-electron peak
35	//
36	// Used numbers are the following:
37	//
38	// Electronic conversion factor:
39	// Assume, we have N_e electrons at the anode,
40	// thus a charge of N_e*e (e = electron charge) Coulomb.
41	//
42	// This charge is AC coupled and runs into a R_pre = 50 Ohm resistency.
43	// The corresponding current is amplified by a gain factor G_pre = 400
44	// (the precision of this value still has to be checked !!!) and again AC coupled to
45	// the output.
46	// The corresponding signal goes through the whole transmission and
47	// amplification chain and is digitized in the FADCs.
48	// The conversion Signal Area to FADC counts (Conv_trans) has been measured
49	// by David and Oscar to be approx. 3.9 pVs^-1
50	//
51	// Thus: Conversion FADC counts to Number of Electrons at Anode:
52	// FADC counts = (1/Conv_tran) * G_pre * R_pre * e * N_e = 8 * 10^-4 N_e.
53	//
54	// Also: FADC counts = 810^-4 GAIN * N_phe
55	//
56	// In the blind pixel, there is an additional pre-amplifier with an amplification of
57	// about 10. Therefore, we have for the blind pixel:
58	//
59	// FADC counts (Blind Pixel) = 810^-3 GAIN * N_phe
60	//
61	//////////////////////////////////////////////////////////////////////////////
62	#include "MHCalibrationBlindPixel.h"
63	#include "MHCalibrationConfig.h"
64
65
66	#include <TStyle.h>
67	#include <TCanvas.h>
68	#include <TPaveText.h>
69
70	#include <TGraph.h>
71	#include <TF1.h>
72	#include <TH1.h>
73	#include <TRandom.h>
74
75	#include "MFFT.h"
76	#include "MLog.h"
77	#include "MLogManip.h"
78
79	ClassImp(MHCalibrationBlindPixel);
80
81	using namespace std;
82
83	const Int_t MHCalibrationBlindPixel::fgBlindPixelChargeNbins = 1000;
84	const Int_t MHCalibrationBlindPixel::fgBlindPixelTimeNbins = 22;
85	const Axis_t MHCalibrationBlindPixel::fgBlindPixelTimeFirst = -9.00;
86	const Axis_t MHCalibrationBlindPixel::fgBlindPixelTimeLast = 12.00;
87	const Double_t MHCalibrationBlindPixel::fgBlindPixelElectronicAmp = 0.008;
88	const Double_t MHCalibrationBlindPixel::fgBlindPixelElectronicAmpError = 0.002;
89
90	const Int_t MHCalibrationBlindPixel::fPSDNbins = 50;
91	const Int_t MHCalibrationBlindPixel::fPulserFrequency = 500;
92	// --------------------------------------------------------------------------
93	//
94	// Default Constructor.
95	//
96	MHCalibrationBlindPixel::MHCalibrationBlindPixel(const char name, const char title)
97	: fHBlindPixelPSD(NULL),
98	fSinglePheFit(NULL),
99	fTimeGausFit(NULL),
100	fSinglePhePedFit(NULL),
101	fPSDHiGain(NULL),
102	fPSDLoGain(NULL),
103	fHPSD(NULL),
104	fPSDExpFit(NULL),
105	fChargeXaxis(NULL),
106	fPSDXaxis(NULL),
107	fCurrentSize(1024),
108	fFitLegend(NULL)
109	{
110
111	fName = name ? name : "MHCalibrationBlindPixel";
112	fTitle = title ? title : "Fill the accumulated charges and times all Blind Pixel events and perform fits";
113
114	// Create a large number of bins, later we will rebin
115	fBlindPixelChargefirst = -400.;
116	fBlindPixelChargelast = 1600.;
117
118	fHBlindPixelCharge = new TH1F("HBlindPixelCharge","Distribution of Summed FADC Slices",
119	fgBlindPixelChargeNbins,fBlindPixelChargefirst,fBlindPixelChargelast);
120	fHBlindPixelCharge->SetXTitle("Sum FADC Slices");
121	fHBlindPixelCharge->SetYTitle("Nr. of events");
122	fHBlindPixelCharge->Sumw2();
123	fHBlindPixelCharge->SetDirectory(NULL);
124
125	fHBlindPixelTime = new TH1F("HBlindPixelTime","Distribution of Mean Arrival Times",
126	fgBlindPixelTimeNbins,fgBlindPixelTimeFirst,fgBlindPixelTimeLast);
127	fHBlindPixelTime->SetXTitle("Mean Arrival Times [FADC slice nr]");
128	fHBlindPixelTime->SetYTitle("Nr. of events");
129	fHBlindPixelTime->Sumw2();
130	fHBlindPixelTime->SetDirectory(NULL);
131
132	fHiGains = new TArrayF(fCurrentSize);
133	fLoGains = new TArrayF(fCurrentSize);
134
135	Clear();
136	}
137
138	MHCalibrationBlindPixel::~MHCalibrationBlindPixel()
139	{
140
141	if (fFitLegend)
142	delete fFitLegend;
143
144	delete fHBlindPixelCharge;
145	delete fHBlindPixelTime;
146
147	delete fHiGains;
148	delete fLoGains;
149
150	if (fHBlindPixelPSD)
151	delete fHBlindPixelPSD;
152	if (fSinglePheFit)
153	delete fSinglePheFit;
154	if (fTimeGausFit)
155	delete fTimeGausFit;
156	if(fSinglePhePedFit)
157	delete fSinglePhePedFit;
158	if (fPSDExpFit)
159	delete fPSDExpFit;
160	if (fHPSD)
161	delete fHPSD;
162	if (fChargeXaxis)
163	delete fChargeXaxis;
164	if (fPSDXaxis)
165	delete fPSDXaxis;
166
167	}
168
169	void MHCalibrationBlindPixel::Clear(Option_t *o)
170	{
171
172	fTotalEntries = 0;
173	fCurrentSize = 1024;
174
175	fBlindPixelChargefirst = -400.;
176	fBlindPixelChargelast = 1600.;
177
178	fLambda = 0.;
179	fMu0 = 0.;
180	fMu1 = 0.;
181	fSigma0 = 0.;
182	fSigma1 = 0.;
183
184	fLambdaErr = 0.;
185	fMu0Err = 0.;
186	fMu1Err = 0.;
187	fSigma0Err = 0.;
188	fSigma1Err = 0.;
189
190	fChisquare = -1.;
191	fProb = -1.;
192	fNdf = -1;
193
194	fMeanTime = -1.;
195	fMeanTimeErr = -1.;
196	fSigmaTime = -1.;
197	fSigmaTimeErr = -1.;
198
199	fLambdaCheck = -1.;
200	fLambdaCheckErr = -1.;
201
202	fMeanPedestal = 0.;
203	fMeanPedestalErr = 0.;
204	fSigmaPedestal = 0.;
205	fSigmaPedestalErr = 0.;
206
207	fFitFunc = kEPoisson4;
208
209	if (fFitLegend)
210	delete fFitLegend;
211	if (fHBlindPixelPSD)
212	delete fHBlindPixelPSD;
213	if (fSinglePheFit)
214	delete fSinglePheFit;
215	if (fTimeGausFit)
216	delete fTimeGausFit;
217	if(fSinglePhePedFit)
218	delete fSinglePhePedFit;
219	if (fPSDExpFit)
220	delete fPSDExpFit;
221	if (fHPSD)
222	delete fHPSD;
223	if (fChargeXaxis)
224	delete fChargeXaxis;
225	if (fPSDXaxis)
226	delete fPSDXaxis;
227	if (fPSDHiGain)
228	delete fPSDHiGain;
229	if (fPSDLoGain)
230	delete fPSDLoGain;
231
232	CLRBIT(fFlags,kFitOK);
233	CLRBIT(fFlags,kOscillating);
234
235	return;
236	}
237
238	void MHCalibrationBlindPixel::Reset()
239	{
240
241	Clear();
242
243	fHBlindPixelCharge->Reset();
244	fHBlindPixelTime->Reset();
245
246	fHiGains->Set(1024);
247	fLoGains->Set(1024);
248
249	fHiGains->Reset(0.);
250	fLoGains->Reset(0.);
251
252
253	}
254
255	Bool_t MHCalibrationBlindPixel::CheckOscillations()
256	{
257
258	if (fPSDExpFit)
259	return IsOscillating();
260
261	// This cuts only the non-used zero's, but MFFT will later cut the rest
262	CutArrayBorder(fHiGains);
263	CutArrayBorder(fLoGains);
264
265	MFFT fourier;
266
267	fPSDLoGain = fourier.PowerSpectrumDensity(fLoGains);
268	fPSDHiGain = fourier.PowerSpectrumDensity(fHiGains);
269
270	fHPSD = ProjectArray(*fPSDHiGain, fPSDNbins,
271	"PSDProjectionBlindPixel",
272	"Power Spectrum Density Projection HiGain Blind Pixel");
273
274	//
275	// First guesses for the fit (should be as close to reality as possible,
276	//
277	const Double_t xmax = fHPSD->GetXaxis()->GetXmax();
278
279	fPSDExpFit = new TF1("PSDExpFit","exp([0]-[1]*x)",0.,xmax);
280
281	const Double_t slope_guess = (TMath::Log(fHPSD->GetEntries())+1.)/xmax;
282	const Double_t offset_guess = slope_guess*xmax;
283
284	fPSDExpFit->SetParameters(offset_guess, slope_guess);
285	fPSDExpFit->SetParNames("Offset","Slope");
286	fPSDExpFit->SetParLimits(0,offset_guess/2.,2.*offset_guess);
287	fPSDExpFit->SetParLimits(1,slope_guess/1.5,1.5*slope_guess);
288	fPSDExpFit->SetRange(0.,xmax);
289
290	fHPSD->Fit(fPSDExpFit,"RQL0");
291
292	fPSDProb = fPSDExpFit->GetProb();
293
294	if (fPSDProb < gkProbLimit)
295	{
296	SETBIT(fFlags,kOscillating);
297	return kTRUE;
298	}
299
300	CLRBIT(fFlags,kOscillating);
301
302	return kFALSE;
303	}
304
305	void MHCalibrationBlindPixel::CreatePSDXaxis(Int_t n)
306	{
307
308	if (fPSDXaxis)
309	return;
310
311	fPSDXaxis = new TArrayF(n);
312
313	for (Int_t i=0;i<n;i++)
314	for (Int_t i=0;i<n;i++)
315	fPSDXaxis->AddAt((Float_t)(fPulserFrequency*i)/2./n,i);
316	}
317
318	void MHCalibrationBlindPixel::CreateChargeXaxis(Int_t n)
319	{
320
321	if (!fChargeXaxis)
322	{
323	fChargeXaxis = new TArrayF(n);
324	for (Int_t i=0;i<n;i++)
325	fChargeXaxis->AddAt((Float_t)i/fPulserFrequency,i);
326	return;
327	}
328
329	if (fChargeXaxis->GetSize() == n)
330	return;
331
332	const Int_t diff = fChargeXaxis->GetSize()-n;
333	fChargeXaxis->Set(n);
334	if (diff < 0)
335	for (Int_t i=n;i<n+diff;i++)
336	fChargeXaxis->AddAt((Float_t)i/fPulserFrequency,i);
337	}
338
339	void MHCalibrationBlindPixel::CutArrayBorder(TArrayF *array)
340	{
341
342	Int_t i;
343
344	for (i=array->GetSize()-1;i>=0;i--)
345	if (array->At(i) != 0)
346	{
347	array->Set(i+1);
348	break;
349	}
350	}
351
352	const Bool_t MHCalibrationBlindPixel::IsFitOK() const
353	{
354	return TESTBIT(fFlags,kFitOK);
355	}
356
357	const Bool_t MHCalibrationBlindPixel::IsOscillating()
358	{
359
360	if (fPSDExpFit)
361	return TESTBIT(fFlags,kOscillating);
362
363	return CheckOscillations();
364
365	}
366
367	Bool_t MHCalibrationBlindPixel::FillGraphs(Float_t qhi,Float_t qlo)
368	{
369
370	if (fTotalEntries >= fCurrentSize)
371	{
372	fCurrentSize *= 2;
373
374	fHiGains->Set(fCurrentSize);
375	fLoGains->Set(fCurrentSize);
376	}
377
378	fHiGains->AddAt(qhi,fTotalEntries);
379	fLoGains->AddAt(qlo,fTotalEntries);
380
381	fTotalEntries++;
382
383	return kTRUE;
384
385	}
386
387
388
389	Bool_t MHCalibrationBlindPixel::FillBlindPixelCharge(Float_t q)
390	{
391	return fHBlindPixelCharge->Fill(q) > -1;
392	}
393
394	Bool_t MHCalibrationBlindPixel::FillBlindPixelTime(Float_t t)
395	{
396	return fHBlindPixelTime->Fill(t) > -1;
397	}
398
399
400	// -------------------------------------------------------------------------
401	//
402	// Draw a legend with the fit results
403	//
404	void MHCalibrationBlindPixel::DrawLegend()
405	{
406
407	fFitLegend = new TPaveText(0.05,0.05,0.95,0.95);
408
409	if (IsFitOK())
410	fFitLegend->SetFillColor(80);
411	else
412	fFitLegend->SetFillColor(2);
413
414	fFitLegend->SetLabel("Results of the single PhE Fit (to k=6):");
415	fFitLegend->SetTextSize(0.05);
416
417	const TString line1 =
418	Form("Mean: #lambda = %2.2f #pm %2.2f",GetLambda(),GetLambdaErr());
419	TText *t1 = fFitLegend->AddText(line1.Data());
420	t1->SetBit(kCanDelete);
421
422	const TString line6 =
423	Form("Mean #lambda (check) = %2.2f #pm %2.2f",GetLambdaCheck(),GetLambdaCheckErr());
424	TText *t2 = fFitLegend->AddText(line6.Data());
425	t2->SetBit(kCanDelete);
426
427	const TString line2 =
428	Form("Pedestal: #mu_{0} = %2.2f #pm %2.2f",GetMu0(),GetMu0Err());
429	TText *t3 = fFitLegend->AddText(line2.Data());
430	t3->SetBit(kCanDelete);
431
432	const TString line3 =
433	Form("Width Pedestal: #sigma_{0} = %2.2f #pm %2.2f",GetSigma0(),GetSigma0Err());
434	TText *t4 = fFitLegend->AddText(line3.Data());
435	t4->SetBit(kCanDelete);
436
437	const TString line4 =
438	Form("1^{st} Phe-peak: #mu_{1} = %2.2f #pm %2.2f",GetMu1(),GetMu1Err());
439	TText *t5 = fFitLegend->AddText(line4.Data());
440	t5->SetBit(kCanDelete);
441
442	const TString line5 =
443	Form("Width 1^{st} Phe-peak: #sigma_{1} = %2.2f #pm %2.2f",GetSigma1(),GetSigma1Err());
444	TText *t6 = fFitLegend->AddText(line5.Data());
445	t6->SetBit(kCanDelete);
446
447	const TString line7 =
448	Form("#chi^{2} / N_{dof}: %4.2f / %3i",GetChiSquare(),GetNdf());
449	TText *t7 = fFitLegend->AddText(line7.Data());
450	t7->SetBit(kCanDelete);
451
452	const TString line8 =
453	Form("Probability: %4.2f ",GetProb());
454	TText *t8 = fFitLegend->AddText(line8.Data());
455	t8->SetBit(kCanDelete);
456
457	if (IsFitOK())
458	{
459	TText *t = fFitLegend->AddText(0.,0.,"Result of the Fit: OK");
460	t->SetBit(kCanDelete);
461	}
462	else
463	{
464	TText *t = fFitLegend->AddText("Result of the Fit: NOT OK");
465	t->SetBit(kCanDelete);
466	}
467
468	fFitLegend->SetBit(kCanDelete);
469	fFitLegend->Draw();
470
471	return;
472	}
473
474	TObject MHCalibrationBlindPixel::DrawClone(Option_t option) const
475	{
476
477	gROOT->SetSelectedPad(NULL);
478
479	MHCalibrationBlindPixel newobj = (MHCalibrationBlindPixel)Clone();
480
481	if (!newobj)
482	return 0;
483	newobj->SetBit(kCanDelete);
484
485	if (strlen(option))
486	newobj->Draw(option);
487	else
488	newobj->Draw(GetDrawOption());
489
490	return newobj;
491	}
492
493
494	// -------------------------------------------------------------------------
495	//
496	// Draw the histogram
497	//
498	void MHCalibrationBlindPixel::Draw(Option_t *opt)
499	{
500
501	gStyle->SetOptFit(1);
502	gStyle->SetOptStat(111111);
503
504	TCanvas *c = MakeDefCanvas(this,550,700);
505
506	c->Divide(2,5);
507
508	gROOT->SetSelectedPad(NULL);
509
510	c->cd(1);
511	gPad->SetLogx(0);
512	gPad->SetLogy(1);
513	gPad->SetTicks();
514
515	fHBlindPixelCharge->Draw(opt);
516
517	if (IsFitOK())
518	fSinglePheFit->SetLineColor(kGreen);
519	else
520	fSinglePheFit->SetLineColor(kRed);
521
522	fSinglePheFit->Draw("same");
523	c->Modified();
524	c->Update();
525
526	fSinglePhePedFit->SetLineColor(kBlue);
527	fSinglePhePedFit->Draw("same");
528
529	c->cd(2);
530	DrawLegend();
531	c->Modified();
532	c->Update();
533
534	c->cd(3);
535	gPad->SetLogy(1);
536	gPad->SetBorderMode(0);
537	fHBlindPixelTime->Draw(opt);
538
539	CutArrayBorder(fHiGains);
540	CreateChargeXaxis(fHiGains->GetSize());
541
542	c->cd(5);
543	gPad->SetTicks();
544	TGraph *gr1 = new TGraph(fChargeXaxis->GetSize(),
545	fChargeXaxis->GetArray(),
546	fHiGains->GetArray());
547	gr1->SetTitle("Evolution of HiGain charges with time");
548	gr1->GetXaxis()->SetTitle("Time [s]");
549	gr1->GetYaxis()->SetTitle("Sum FADC slices");
550	gr1->SetBit(kCanDelete);
551	gr1->Draw("AL");
552
553	CutArrayBorder(fLoGains);
554	CreateChargeXaxis(fHiGains->GetSize());
555
556	c->cd(6);
557	gPad->SetTicks();
558	TGraph *gr2 = new TGraph(fChargeXaxis->GetSize(),
559	fChargeXaxis->GetArray(),
560	fLoGains->GetArray());
561	gr2->SetTitle("Evolution of LoGain charges with time");
562	gr2->GetXaxis()->SetTitle("Time [s]");
563	gr2->GetYaxis()->SetTitle("Sum FADC slices");
564	gr2->SetBit(kCanDelete);
565	gr2->Draw("AL");
566
567	c->Modified();
568	c->Update();
569
570	c->cd(7);
571
572	TArrayF *array;
573
574	if (!fPSDHiGain)
575	return;
576	array = fPSDHiGain;
577
578	if (!fPSDXaxis)
579	CreatePSDXaxis(array->GetSize());
580
581	TGraph *gr3 = new TGraph(fPSDXaxis->GetSize(),fPSDXaxis->GetArray(),array->GetArray());
582	gr3->SetTitle("Power Spectrum Density");
583	gr3->GetXaxis()->SetTitle("Freq. [Hz]");
584	gr3->GetYaxis()->SetTitle("P(f)");
585	gr3->SetBit(kCanDelete);
586	gr3->Draw("AL");
587
588	c->Modified();
589	c->Update();
590
591	c->cd(8);
592
593	gStyle->SetOptStat(111111);
594	gPad->SetTicks();
595
596	if (fHPSD->Integral() > 0)
597	gPad->SetLogy();
598
599	fHPSD->Draw(opt);
600	fHPSD->GetXaxis()->SetTitle("P(f)");
601	fHPSD->GetYaxis()->SetTitle("Counts");
602
603	if (fPSDExpFit)
604	{
605	fPSDExpFit->SetLineColor(IsOscillating() ? kRed : kGreen);
606	fPSDExpFit->Draw("same");
607	}
608
609	c->cd(9);
610
611	array = fPSDLoGain;
612
613	if (!fPSDXaxis)
614	CreatePSDXaxis(array->GetSize());
615
616	TGraph *gr4 = new TGraph(fPSDXaxis->GetSize(),fPSDXaxis->GetArray(),array->GetArray());
617	gr4->SetTitle("Power Spectrum Density");
618	gr4->GetXaxis()->SetTitle("Freq. [Hz]");
619	gr4->GetYaxis()->SetTitle("P(f)");
620	gr4->SetBit(kCanDelete);
621	gr4->Draw("AL");
622
623	c->Modified();
624	c->Update();
625	}
626
627
628
629	Bool_t MHCalibrationBlindPixel::SimulateSinglePhe(Double_t lambda, Double_t mu0, Double_t mu1, Double_t sigma0, Double_t sigma1)
630	{
631	gRandom->SetSeed();
632
633	if (fHBlindPixelCharge->GetIntegral() != 0)
634	{
635	*fLog << err << "Histogram " << fHBlindPixelCharge->GetTitle() << " is already filled. " << endl;
636	*fLog << err << "Create new class MHCalibrationBlindPixel for simulation! " << endl;
637	return kFALSE;
638	}
639
640	if (!InitFit(fBlindPixelChargefirst,fBlindPixelChargelast))
641	return kFALSE;
642
643	for (Int_t i=0;i<10000; i++)
644	fHBlindPixelCharge->Fill(fSinglePheFit->GetRandom());
645
646	return kTRUE;
647	}
648
649	Bool_t MHCalibrationBlindPixel::InitFit(Axis_t min, Axis_t max)
650	{
651
652	//
653	// First guesses for the fit (should be as close to reality as possible,
654	// otherwise the fit goes gaga because of high number of dimensions ...
655	//
656	const Stat_t entries = fHBlindPixelCharge->Integral("width");
657	const Double_t lambda_guess = 0.5;
658	const Double_t maximum_bin = fHBlindPixelCharge->GetBinCenter(fHBlindPixelCharge->GetMaximumBin());
659	const Double_t norm = entries/gkSq2Pi;
660
661	//
662	// Initialize the fit function
663	//
664	switch (fFitFunc)
665	{
666	case kEPoisson4:
667	fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto4,min,max,6);
668	break;
669	case kEPoisson5:
670	fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto5,min,max,6);
671	break;
672	case kEPoisson6:
673	fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto6,min,max,6);
674	break;
675	case kEPolya:
676	fSinglePheFit = new TF1("SinglePheFit",&fPolya,min,max,8);
677	break;
678	case kEMichele:
679	break;
680
681	default:
682	*fLog << warn << "WARNING: Could not find Fit Function for Blind Pixel " << endl;
683	return kFALSE;
684	break;
685	}
686
687	const Double_t mu_0_guess = maximum_bin;
688	const Double_t si_0_guess = 40.;
689	const Double_t mu_1_guess = mu_0_guess + 100.;
690	const Double_t si_1_guess = si_0_guess + si_0_guess;
691	// Michele
692	const Double_t lambda_1cat_guess = 0.5;
693	const Double_t lambda_1dyn_guess = 0.5;
694	const Double_t mu_1cat_guess = mu_0_guess + 50.;
695	const Double_t mu_1dyn_guess = mu_0_guess + 20.;
696	const Double_t si_1cat_guess = si_0_guess + si_0_guess;
697	const Double_t si_1dyn_guess = si_0_guess;
698	// Polya
699	const Double_t excessPoisson_guess = 0.5;
700	const Double_t delta1_guess = 8.;
701	const Double_t delta2_guess = 5.;
702	const Double_t electronicAmp_guess = fgBlindPixelElectronicAmp;
703	const Double_t electronicAmp_limit = fgBlindPixelElectronicAmpError;
704
705	//
706	// Initialize boundaries and start parameters
707	//
708	switch (fFitFunc)
709	{
710
711	case kEPoisson4:
712	if ((fMeanPedestal) && (fSigmaPedestal))
713	fSinglePheFit->SetParameters(lambda_guess,fMeanPedestal,mu_1_guess,fSigmaPedestal,si_1_guess,norm);
714	else
715	fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm);
716
717	fSinglePheFit->SetParNames("#lambda","#mu_{0}","#mu_{1}","#sigma_{0}","#sigma_{1}","Area");
718
719	fSinglePheFit->SetParLimits(0,0.,1.);
720
721	if ((fMeanPedestal) && (fSigmaPedestal))
722	fSinglePheFit->SetParLimits(1,
723	fMeanPedestal-1.*fMeanPedestalErr,
724	fMeanPedestal+1.*fMeanPedestalErr);
725	else
726	fSinglePheFit->SetParLimits(1,-3.,0.);
727
728	fSinglePheFit->SetParLimits(2,(max-min)/2.,max);
729
730	if ((fMeanPedestal) && (fSigmaPedestal))
731	fSinglePheFit->SetParLimits(3,
732	fSigmaPedestal-3.*fSigmaPedestalErr,
733	fSigmaPedestal+3.*fSigmaPedestalErr);
734	else
735	fSinglePheFit->SetParLimits(3,1.0,(max-min)/2.0);
736
737	fSinglePheFit->SetParLimits(4,1.0,(max-min));
738	fSinglePheFit->SetParLimits(5,norm-0.5,norm+0.5);
739	break;
740	case kEPoisson5:
741	case kEPoisson6:
742	fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm);
743	fSinglePheFit->SetParNames("#lambda","#mu_{0}","#mu_{1}","#sigma_{0}","#sigma_{1}","Area");
744	fSinglePheFit->SetParLimits(0,0.,1.);
745	fSinglePheFit->SetParLimits(1,min,(max-min)/1.5);
746	fSinglePheFit->SetParLimits(2,(max-min)/2.,(max-0.05*(max-min)));
747	fSinglePheFit->SetParLimits(3,1.0,(max-min)/2.0);
748	fSinglePheFit->SetParLimits(4,1.0,(max-min)/2.5);
749	fSinglePheFit->SetParLimits(5,norm-0.1,norm+0.1);
750	break;
751
752	case kEPolya:
753	if ((fMeanPedestal) && (fSigmaPedestal))
754	fSinglePheFit->SetParameters(lambda_guess, excessPoisson_guess,
755	delta1_guess,delta2_guess,
756	electronicAmp_guess,
757	fSigmaPedestal,
758	norm,
759	fMeanPedestal);
760	else
761	fSinglePheFit->SetParameters(lambda_guess, excessPoisson_guess,
762	delta1_guess,delta2_guess,
763	electronicAmp_guess,
764	si_0_guess,
765	norm, mu_0_guess);
766	fSinglePheFit->SetParNames("#lambda","b_{tot}",
767	"#delta_{1}","#delta_{2}",
768	"amp_{e}","#sigma_{0}",
769	"Area", "#mu_{0}");
770	fSinglePheFit->SetParLimits(0,0.,1.);
771	fSinglePheFit->SetParLimits(1,0.,1.);
772	fSinglePheFit->SetParLimits(2,6.,12.);
773	fSinglePheFit->SetParLimits(3,3.,8.);
774	fSinglePheFit->SetParLimits(4,electronicAmp_guess-electronicAmp_limit,
775	electronicAmp_guess+electronicAmp_limit);
776	if ((fMeanPedestal) && (fSigmaPedestal))
777	fSinglePheFit->SetParLimits(5,
778	fSigmaPedestal-3.*fSigmaPedestalErr,
779	fSigmaPedestal+3.*fSigmaPedestalErr);
780	else
781	fSinglePheFit->SetParLimits(5,min,(max-min)/1.5);
782
783	fSinglePheFit->SetParLimits(6,norm-0.1,norm+0.1);
784	if ((fMeanPedestal) && (fSigmaPedestal))
785	fSinglePheFit->SetParLimits(7,
786	fMeanPedestal-3.*fMeanPedestalErr,
787	fMeanPedestal+3.*fMeanPedestalErr);
788	else
789	fSinglePheFit->SetParLimits(7,-35.,15.);
790	break;
791
792	case kEMichele:
793
794
795	break;
796
797	default:
798	*fLog << warn << "WARNING: Could not find Fit Function for Blind Pixel " << endl;
799	return kFALSE;
800	break;
801	}
802
803	return kTRUE;
804	}
805
806	void MHCalibrationBlindPixel::ExitFit(TF1 *f)
807	{
808
809
810	//
811	// Finalize
812	//
813	switch (fFitFunc)
814	{
815
816	case kEPoisson4:
817	case kEPoisson5:
818	case kEPoisson6:
819	case kEPoisson7:
820	fLambda = fSinglePheFit->GetParameter(0);
821	fMu0 = fSinglePheFit->GetParameter(1);
822	fMu1 = fSinglePheFit->GetParameter(2);
823	fSigma0 = fSinglePheFit->GetParameter(3);
824	fSigma1 = fSinglePheFit->GetParameter(4);
825
826	fLambdaErr = fSinglePheFit->GetParError(0);
827	fMu0Err = fSinglePheFit->GetParError(1);
828	fMu1Err = fSinglePheFit->GetParError(2);
829	fSigma0Err = fSinglePheFit->GetParError(3);
830	fSigma1Err = fSinglePheFit->GetParError(4);
831	break;
832	case kEPolya:
833	fLambda = fSinglePheFit->GetParameter(0);
834	fMu0 = fSinglePheFit->GetParameter(7);
835	fMu1 = 0.;
836	fSigma0 = fSinglePheFit->GetParameter(5);
837	fSigma1 = 0.;
838
839	fLambdaErr = fSinglePheFit->GetParError(0);
840	fMu0Err = fSinglePheFit->GetParError(7);
841	fMu1Err = 0.;
842	fSigma0Err = fSinglePheFit->GetParError(5);
843	fSigma1Err = 0.;
844	default:
845	break;
846	}
847
848	}
849
850
851	Bool_t MHCalibrationBlindPixel::FitSinglePhe(Axis_t rmin, Axis_t rmax, Option_t *opt)
852	{
853
854	//
855	// Get the fitting ranges
856	//
857	rmin = (rmin != 0.) ? rmin : fBlindPixelChargefirst;
858	rmax = (rmax != 0.) ? rmax : fBlindPixelChargelast;
859
860	if (!InitFit(rmin,rmax))
861	return kFALSE;
862
863	fHBlindPixelCharge->Fit(fSinglePheFit,opt);
864
865	ExitFit(fSinglePheFit);
866
867	fProb = fSinglePheFit->GetProb();
868	fChisquare = fSinglePheFit->GetChisquare();
869	fNdf = fSinglePheFit->GetNDF();
870
871	// Perform the cross-check fitting only the pedestal:
872	fSinglePhePedFit = new TF1("GausPed","gaus",rmin,0.);
873	fHBlindPixelCharge->Fit(fSinglePhePedFit,opt);
874
875	const Stat_t entries = fHBlindPixelCharge->Integral("width");
876
877	Double_t pedarea = fSinglePhePedFit->GetParameter(0)gkSq2PifSinglePhePedFit->GetParameter(2);
878	fLambdaCheck = TMath::Log(entries/pedarea);
879	fLambdaCheckErr = fSinglePhePedFit->GetParError(0)/fSinglePhePedFit->GetParameter(0)
880	+ fSinglePhePedFit->GetParError(2)/fSinglePhePedFit->GetParameter(2);
881
882	*fLog << inf << "Results of the Blind Pixel Fit: " << endl;
883	*fLog << inf << "Chisquare: " << fChisquare << endl;
884	*fLog << inf << "DoF: " << fNdf << endl;
885	*fLog << inf << "Probability: " << fProb << endl;
886
887	//
888	// The fit result is accepted under condition that:
889	// 1) the Probability is greater than gkProbLimit (default 0.001 == 99.7%)
890	// 2) at least 50 events are in the single Photo-electron peak
891	//
892	if (fProb < gkProbLimit)
893	{
894	*fLog << err << "ERROR: Fit Probability " << fProb
895	<< " is smaller than the allowed value: " << gkProbLimit << endl;
896	CLRBIT(fFlags,kFitOK);
897	return kFALSE;
898	}
899
900	Float_t contSinglePhe = TMath::Exp(-1.0fLambda)fLambda*entries;
901
902	if (contSinglePhe < 50.)
903	{
904	*fLog << err << "ERROR: Statistics is too low: Only " << contSinglePhe
905	<< " in the Single Photo-Electron peak " << endl;
906	CLRBIT(fFlags,kFitOK);
907	return kFALSE;
908	}
909	else
910	*fLog << inf << contSinglePhe << " in Single Photo-Electron peak " << endl;
911
912	SETBIT(fFlags,kFitOK);
913
914	return kTRUE;
915	}
916
917
918	void MHCalibrationBlindPixel::CutAllEdges()
919	{
920
921	Int_t nbins = 30;
922
923	CutEdges(fHBlindPixelCharge,nbins);
924
925	fBlindPixelChargefirst = fHBlindPixelCharge->GetBinLowEdge(fHBlindPixelCharge->GetXaxis()->GetFirst());
926	fBlindPixelChargelast = fHBlindPixelCharge->GetBinLowEdge(fHBlindPixelCharge->GetXaxis()->GetLast())+fHBlindPixelCharge->GetBinWidth(0);
927
928	}
929
930	Bool_t MHCalibrationBlindPixel::FitTime(Axis_t rmin, Axis_t rmax, Option_t *opt)
931	{
932
933	rmin = (rmin != 0.) ? rmin : 4.;
934	rmax = (rmax != 0.) ? rmax : 9.;
935
936	const Stat_t entries = fHBlindPixelTime->Integral();
937	const Double_t mu_guess = fHBlindPixelTime->GetBinCenter(fHBlindPixelTime->GetMaximumBin());
938	const Double_t sigma_guess = (rmax - rmin)/2.;
939	const Double_t area_guess = entries/gkSq2Pi;
940
941	fTimeGausFit = new TF1("GausTime","gaus",rmin,rmax);
942	fTimeGausFit->SetParameters(area_guess,mu_guess,sigma_guess);
943	fTimeGausFit->SetParNames("Area","#mu","#sigma");
944	fTimeGausFit->SetParLimits(0,0.,entries);
945	fTimeGausFit->SetParLimits(1,rmin,rmax);
946	fTimeGausFit->SetParLimits(2,0.,rmax-rmin);
947
948	fHBlindPixelTime->Fit(fTimeGausFit,opt);
949	rmin = fTimeGausFit->GetParameter(1) - 2.*fTimeGausFit->GetParameter(2);
950	rmax = fTimeGausFit->GetParameter(1) + 2.*fTimeGausFit->GetParameter(2);
951	fTimeGausFit->SetRange(rmin,rmax);
952
953	fHBlindPixelTime->Fit(fTimeGausFit,opt);
954
955	fMeanTime = fTimeGausFit->GetParameter(2);
956	fSigmaTime = fTimeGausFit->GetParameter(3);
957	fMeanTimeErr = fTimeGausFit->GetParError(2);
958	fSigmaTimeErr = fTimeGausFit->GetParError(3);
959
960	*fLog << inf << "Results of the Times Fit: " << endl;
961	*fLog << inf << "Chisquare: " << fTimeGausFit->GetChisquare() << endl;
962	*fLog << inf << "Ndf: " << fTimeGausFit->GetNDF() << endl;
963
964	return kTRUE;
965
966	}

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

Download in other formats: