Context Navigation

fit_spectra.C@ 20115

Visit:

Last change on this file since 20115 was 19874, checked in by tbretz, 5 years ago
Just found out that this was the wrong version of the macro. The correct one (which also contains the correct function from the paper) is this one.
File size: 29.7 KB

Line
1	#include <iomanip>
2	#include <iostream>
3
4	#include <TF1.h>
5	#include <TH2.h>
6	#include <TProfile2D.h>
7	#include <TFile.h>
8	#include <TFitResult.h>
9	#include <TStyle.h>
10
11	#include "MLog.h"
12	#include "MLogManip.h"
13	#include "MStatusArray.h"
14	#include "MStatusDisplay.h"
15	#include "MHCamera.h"
16	#include "MGeomCamFACT.h"
17	#include "MParameters.h"
18	#include "MArrayI.h"
19	#include "MRawRunHeader.h"
20	#include "PixelMap.h"
21
22	using namespace std;
23
24	// --------------------------------------------------------------------------
25	// After fit_spectra_bt2b.C
26
27	// Fit function for a single pe spectrum
28	Double_t fcn_g(Double_t xx, Double_t par)
29	{
30	const Double_t ampl = par[0];
31	const Double_t gain = par[1];
32	const Double_t sigma = par[2]*gain;
33	const Double_t cross = par[3];
34	const Double_t shift = par[4];
35	const Double_t noise = par[5]<0 ? sigma : par[5];
36	const Double_t expo = par[6];
37
38	const Double_t P = cross*TMath::Exp(-cross);
39
40	Double_t y = 0;
41	for (int N=1; N<14; N++)
42	{
43	const Double_t muN = N*gain + shift;
44	const Double_t sigN = TMath::Sqrt(Nsigmasigma + noise*noise);
45
46	const Double_t p =
47	TMath::Power(N*P, N-1)/TMath::Power(TMath::Factorial(N-1), expo);
48
49	y += TMath::Gaus(xx[0], muN, sigN) * p / sigN;
50	}
51
52	const Double_t sig1 = TMath::Sqrt(sigmasigma + noisenoise);
53	return amplsig1y;
54	}
55
56	// Calculate the crosstalk from the function parameters
57	Double_t xtalk(TF1 &f)
58	{
59	Double_t cross = f.GetParameter(3);
60	Double_t expo = f.GetParameter(6);
61
62	const Double_t P = cross*TMath::Exp(-cross);
63
64	Double_t y = 0;
65	for (int N=2; N<14; N++)
66	y +=
67	TMath::Power(N*P, N-1)/TMath::Power(TMath::Factorial(N-1), expo);
68
69	return y / (y + 1);
70	}
71
72	// calculate the integral in units per millisecond
73	double integral(TF1 &func, TH1 &hist)
74	{
75	const Double_t sigma = func.GetParameter(2)*func.GetParameter(1);
76	const Double_t cross = func.GetParameter(3);
77	const Double_t noise = func.GetParameter(5)<0 ? sigma : func.GetParameter(5);
78	const Double_t expo = func.GetParameter(6);
79
80	const Double_t P = cross*TMath::Exp(-cross);
81
82	Double_t sum = 0;
83	for (int N=1; N<14; N++)
84	sum +=
85	TMath::Power(N*P, N-1)/TMath::Power(TMath::Factorial(N-1), expo);
86
87	const Double_t scale = hist.GetBinWidth(1);
88	const Double_t sig1 = TMath::Sqrt(sigmasigma + noisenoise);
89	const Double_t integ = func.GetParameter(0)sumsig1*sqrt(TMath::TwoPi())/scale;
90
91	return integ/1e-9/1e6;
92	}
93
94
95	// --------------------------------------------------------------------------
96
97	// Print function parameters
98	void PrintFunc(TF1 &f, float integration_window=30)
99	{
100	//cout << "--------------------" << endl;
101	cout << "Ampl: " << setw(8) << f.GetParameter(0) << " +/- " << f.GetParError(0) << endl;
102	cout << "Gain: " << setw(8) << f.GetParameter(1) << " +/- " << f.GetParError(1) << endl;
103	cout << "Rel.sigma: " << setw(8) << f.GetParameter(2) << " +/- " << f.GetParError(2) << endl;
104	cout << "Baseline: " << setw(8) << f.GetParameter(4)/integration_window << " +/- " << f.GetParError(4)/integration_window << endl;
105	cout << "Crosstalk: " << setw(8) << f.GetParameter(3) << " +/- " << f.GetParError(3) << endl;
106	cout << "Pcrosstalk: " << setw(8) << xtalk(f) << endl;
107	if (f.GetParameter(5)>=0)
108	cout << "Noise: " << setw(8) << f.GetParameter(5)/sqrt(integration_window) << " +/- " << f.GetParError(5)/sqrt(integration_window) << endl;
109	cout << "Expo: " << setw(8) << f.GetParameter(6) << " +/- " << f.GetParError(6) << endl;
110	//cout << "--------------------" << endl;
111
112	}
113
114	// --------------------------------------------------------------------------
115
116	// The parameters for the function are the filenam from the gain extraction
117	// and the output filename
118	int fit_spectra(const char *filename = "20130222_018_018.root",
119	const char *outfile = "20130222_018_018-fit.root",
120	bool fixednoise=true)
121	{
122	// Read the mapping between bias channels and hardware channels
123	PixelMap pmap;
124	if (!pmap.Read("FACTmap111030.txt"))
125	{
126	cout << "FACTmap111030.txt not found." << endl;
127	return 1;
128	}
129
130	// It is more correct to fit the integral, but this is super
131	// slow, especially for 1440 pixel. To get that, one would have
132	// to analytically integrate and fit this function.
133	// (Currently the integral is switched off for the 1440 individual
134	// spectra in all cases)
135	bool fast = false; // Switch off using integral
136
137	// Values which should be read from the file but not available atm
138	Int_t integration_window = 30;
139
140	// Map for which pixel shall be plotted and which not
141	TArrayC usePixel(1440);
142	memset(usePixel.GetArray(), 1, 1440);
143
144	// List of Pixel that should be ignored in camera view
145	usePixel[424] = 0;
146	usePixel[923] = 0;
147	usePixel[1208] = 0;
148	usePixel[583] = 0;
149	usePixel[830] = 0;
150	usePixel[1399] = 0;
151	usePixel[113] = 0;
152	usePixel[115] = 0;
153	usePixel[354] = 0;
154	usePixel[423] = 0;
155	usePixel[1195] = 0;
156	usePixel[1393] = 0;
157
158	cout << setprecision(3);
159
160	// ======================================================
161	// Read data and histograms from file
162
163	TFile file(filename);
164	if (file.IsZombie())
165	{
166	gLog << err << "Opening file '" << filename << "' failed." << endl;
167	return 1;
168	}
169
170	MStatusArray arr;
171	if (arr.Read()<=0)
172	{
173	gLog << err << "Reading of MStatusArray from '" << filename << "' failed." << endl;
174	return 2;
175	}
176
177	TH2 hsignal = (TH2)arr.FindObjectInCanvas("Signal", "TH2F", "MHSingles");
178	if (!hsignal)
179	{
180	gLog << err << "Histogram Signal not found in '" << filename << "'." << endl;
181	return 3;
182	}
183	TH2 htime = (TH2)arr.FindObjectInCanvas("Time", "TH2F", "MHSingles");
184	if (!htime)
185	{
186	gLog << err << "Histogram Time not found in '" << filename << "'." << endl;
187	return 4;
188	}
189	TProfile2D hpulse = (TProfile2D)arr.FindObjectInCanvas("Pulse", "TProfile2D", "MHSingles");
190	if (!hpulse)
191	{
192	gLog << err << "Histogram Pulse not found in '" << filename << "'." << endl;
193	return 5;
194	}
195	TH2F hbase = (TH2F)arr.FindObjectInCanvas("Baseline", "TH2F", "MHBaseline");
196	if (!hbase)
197	{
198	gLog << err << "Histogram Baseline not found in '" << filename << "'." << endl;
199	return 6;
200	}
201
202	MRawRunHeader header;
203	if (header.Read()<=0)
204	{
205	gLog << err << "MRawRunheader not found in '" << filename << "'." << endl;
206	return 7;
207	}
208
209	MParameterI par("NumEvents");
210	if (par.Read()<=0)
211	{
212	gLog << err << "NumEvents not found in '" << filename << "'." << endl;
213	return 8;
214	}
215
216	MArrayI ext;
217	if (ext.Read("ExtractionRange")<=0)
218	{
219	gLog << err << "ExtractionRange not found in '" << filename << "'." << endl;
220	return 9;
221
222	}
223
224	// ======================================================
225
226	MStatusDisplay *d = new MStatusDisplay;
227
228	// Camera geometry for displays
229	MGeomCamFACT fact;
230
231	// ------------------ Spectrum Fit ---------------
232	// Instantiate the display histograms
233	MHCamera cRate(fact);
234	MHCamera cGain(fact);
235	MHCamera cRelSigma(fact);
236	MHCamera cCrosstalk(fact);
237	MHCamera cBaseline(fact);
238	MHCamera cNoise(fact);
239	MHCamera cChi2(fact);
240	MHCamera cNormGain(fact);
241	MHCamera cFitProb(fact);
242	MHCamera cCrosstalkP(fact);
243	MHCamera cCoeffR(fact);
244
245	// Set name and title for the histograms
246	cRate.SetNameTitle ("Rate", "Dark count rate");
247	cGain.SetNameTitle ("Gain", "Gain distribution");
248	cRelSigma.SetNameTitle ("RelSigma", "Rel. Sigma");
249	cCrosstalk.SetNameTitle ("Crosstalk", "Crosstalk probability");
250	cBaseline.SetNameTitle ("Baseline", "Baseline per sample");
251	cNoise.SetNameTitle ("Noise", "Noise per sample");
252	cChi2.SetNameTitle ("Chi2", "\\Chi^2");
253	cNormGain.SetNameTitle ("NormGain", "Normalized gain");
254	cFitProb.SetNameTitle ("FitProb", "Root's fit probability");
255	cCrosstalkP.SetNameTitle("Pxtalk", "Crosstalk coeff. P");
256	cCoeffR.SetNameTitle ("CoeffR", "Coefficient R");
257
258	// Instantiate 1D histograms for the distributions
259	// including TM channels
260	TH1F hRate1 ("Rate1", "Dark count rate", 150, 0, 15);
261	TH1F hGain1 ("Gain1", "Gain distribution", 100, 0, 400);
262	TH1F hRelSigma1 ("RelSigma1", "Rel. Sigma", 160, 0, 0.40);
263	TH1F hCrosstalk1 ("Crosstalk1", "Crosstalk probability", 90, 0, 0.30);
264	TH1F hBaseline1 ("Baseline1", "Baseline per sample", 75, -7.5, 7.5);
265	TH1F hNoise1 ("Noise1", "Noise per sample", 60, 0, 30);
266	TH1F hChiSq1 ("ChiSq1", "\\Chi^2", 200, 0, 4);
267	TH1F hNormGain1 ("NormGain1", "Normalized gain", 51, 0.5, 1.5);
268	TH1F hFitProb1 ("FitProb1", "FitProb distribution", 100, 0, 1);
269	TH1F hCrosstalkP1("Pxtalk1", "Crosstalk coeff.", 90, 0, 0.3);
270	TH1F hCoeffR1 ("CoeffR1", "Coefficient R", 90, -1, 2);
271
272	// excluding TM channels
273	TH1F hRate2 ("Rate2", "Dark count rate", 150, 0, 15);
274	TH1F hGain2 ("Gain2", "Gain distribution", 100, 0, 400);
275	TH1F hRelSigma2 ("RelSigma2", "Rel. Sigma", 160, 0, 0.40);
276	TH1F hCrosstalk2 ("Crosstalk2", "Crosstalk probability", 90, 0, 0.30);
277	TH1F hBaseline2 ("Baseline2", "Baseline per sample", 75, -7.5, 7.5);
278	TH1F hNoise2 ("Noise2", "Noise per sample", 60, 0, 30);
279	TH1F hChiSq2 ("ChiSq2", "\\Chi^2", 200, 0, 4);
280	TH1F hNormGain2 ("NormGain2", "Normalized gain", 51, 0.5, 1.5);
281	TH1F hFitProb2 ("FitProb2", "FitProb distribution", 100, 0, 1);
282	TH1F hCrosstalkP2("Pxtalk2", "Crosstalk coeff.", 90, 0, 0.3);
283	TH1F hCoeffR2 ("CoeffR2", "Coefficient R", 90, -1, 2);
284
285	// Histigram for the sum of all spectrums
286	TH1F hSum("Sum1", "Signal spectrum of all pixels",
287	hsignal->GetNbinsY(), hsignal->GetYaxis()->GetXmin(), hsignal->GetYaxis()->GetXmax());
288	hSum.SetXTitle("pulse integral [mV sample]");
289	hSum.SetYTitle("Counts");
290	hSum.SetStats(false);
291	hSum.Sumw2();
292
293	// Histogram for the sum of all pixels excluding the ones with faulty fits
294	TH1F hSumClear1("SumC1", "Signal spectrum of all pixels (incl TM)",
295	hsignal->GetNbinsY(), hsignal->GetYaxis()->GetXmin(), hsignal->GetYaxis()->GetXmax());
296	hSumClear1.SetXTitle("pulse integral [mV sample]");
297	hSumClear1.SetYTitle("Counts");
298	hSumClear1.SetStats(false);
299	hSumClear1.SetLineColor(kBlue);
300	hSumClear1.Sumw2();
301
302	TH1F hSumClear2("SumC2", "Signal spectrum of all pixels (excp TM)",
303	hsignal->GetNbinsY(), hsignal->GetYaxis()->GetXmin(), hsignal->GetYaxis()->GetXmax());
304	hSumClear2.SetXTitle("pulse integral [mV sample]");
305	hSumClear2.SetYTitle("Counts");
306	hSumClear2.SetStats(false);
307	hSumClear2.SetLineColor(kBlue);
308	hSumClear2.Sumw2();
309
310	// Arrival time spectrum of the extracted pulses
311	TH1F hTime("Time", "Arrival time spectrum", htime->GetNbinsY(), htime->GetYaxis()->GetXmin(), htime->GetYaxis()->GetXmax());
312	hTime.SetXTitle("pulse arrival time [sample]");
313	hTime.SetYTitle("Counts");
314	hTime.SetStats(false);
315
316	// average pulse shape of the extracted pulses
317	TH1F hPulse("Puls", "Average pulse", hpulse->GetNbinsY(), hpulse->GetYaxis()->GetXmin(), hpulse->GetYaxis()->GetXmax());
318	hPulse.SetXTitle("pulse arrival time [sample]");
319	hPulse.SetYTitle("Counts");
320	hPulse.SetStats(false);
321
322	// Spektrum for the normalized individual spectra
323	TH1F hSumScale1("SumScale1", "Signal spectrum of all pixels (incl TM)", 120, 0.05, 12.05);
324	hSumScale1.SetXTitle("pulse integral [pe]");
325	hSumScale1.SetYTitle("Rate");
326	hSumScale1.SetStats(false);
327	hSumScale1.Sumw2();
328
329	TH1F hSumScale2("SumScale2", "Signal spectrum of all pixels (excl TM)", 120, 0.05, 12.05);
330	hSumScale2.SetXTitle("pulse integral [pe]");
331	hSumScale2.SetYTitle("Rate");
332	hSumScale2.SetStats(false);
333	hSumScale2.Sumw2();
334
335	// define fit function for Amplitudespectrum
336	TF1 func("spektrum", fcn_g, 0, hSum.GetXaxis()->GetXmax(), 7);
337	func.SetNpx(2000);
338	func.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset", "Noise", "Expo");
339	func.SetLineColor(kRed);
340
341	//--------------------- fill sum spectrum --------------------------------
342
343	d->SetStatusLine1("Calculating sum spectrum");
344
345	// Begin of Loop over Pixels
346	for (Int_t pixel = 0; pixel < hsignal->GetNbinsX(); pixel++)
347	{
348	//jump to next pixel if the current is marked as faulty
349	if (usePixel[pixel]==0)
350	continue;
351
352	TH1D *hist = hsignal->ProjectionY("proj", pixel+1, pixel+1);
353	hSum.Add(hist);
354	delete hist;
355	}
356
357	//----------------- get starting values -------------------------------
358
359	hSum.GetXaxis()->SetRangeUser(150, hSum.GetXaxis()->GetXmax());
360
361	const Int_t maxbin = hSum.GetMaximumBin();
362	const Double_t maxpos = hSum.GetBinCenter(maxbin);
363	const Double_t binwidth = hSum.GetBinWidth(maxbin);
364	const Double_t ampl = hSum.GetBinContent(maxbin);
365
366	double fwhmSum = 0;
367
368	//Calculate full width half Maximum
369	for (int i=1; i<maxbin; i++)
370	if (hSum.GetBinContent(maxbin-i)+hSum.GetBinContent(maxbin+i)<ampl)
371	{
372	fwhmSum = 2(i-0.5)binwidth;
373	break;
374	}
375
376	if (fwhmSum==0)
377	{
378	gLog << warn << "Could not determine start value for sigma." << endl;
379	}
380
381	Double_t sigma_est = fwhmSum/2.3548; // FWHM = 2sqrt(2ln(2))*sigma
382
383	Double_t fitmin = maxpos-3*sigma_est; // was 3
384	Double_t fitmax = hSum.GetXaxis()->GetXmax();
385
386	// ------------------- fit --------------------------------
387
388	//Fit and draw spectrum
389	func.SetParLimits(0, 0, 2*ampl); // Amplitude
390	func.SetParLimits(1, 0, 2*maxpos); // Gain
391	func.SetParLimits(2, 0, 1); // Sigma
392	func.SetParLimits(3, 0, 1); // Crosstalk
393	if (!fixednoise)
394	func.SetParLimits(5, 0, 150); // Noise
395	func.SetParLimits(6, 0, 2); // Expo
396
397	func.SetParameter(0, ampl); // Amplitude
398	func.SetParameter(1, maxpos); // Gain
399	func.SetParameter(2, 0.1); // Sigma
400	func.SetParameter(3, 0.16); // Crosstalk
401	func.SetParameter(4, 0*integration_window); // Baseline
402	if (fixednoise)
403	func.FixParameter(5, -1); // Noise
404	else
405	func.SetParameter(5, 0.1*maxpos); // Noise
406
407	func.SetParameter(6, 0.95); // Expo
408
409	func.SetRange(fitmin, fitmax);
410	hSum.Fit(&func, fast?"N0QSR":"IN0QSR");
411
412	Double_t res_par[7];
413	func.GetParameters(res_par);
414
415	//func.FixParameter(6, func.GetParameter(6)); // Expo
416
417	// ------------------ display result -------------------------------
418
419	cout << "--------------------" << endl;
420	cout << "AmplEst: " << ampl << endl;
421	cout << "GainEst: " << maxpos << endl;
422	cout << "SigmaEst: " << sigma_est << endl;
423	PrintFunc(func, integration_window);
424	cout << "--------------------" << endl;
425
426	gROOT->SetSelectedPad(0);
427	TCanvas &c11 = d->AddTab("SumHist");
428	c11.cd();
429	gPad->SetLogy();
430	gPad->SetGridx();
431	gPad->SetGridy();
432	hSum.GetXaxis()->SetRange();
433	hSum.DrawCopy("hist");
434	func.DrawCopy("same");
435
436	// ===================================================================
437	// Gain Calculation for each Pixel
438	// ===================================================================
439
440	// counter for number of processed pixel
441	int count_ok = 0;
442
443	// Begin of Loop over Pixels
444	for (Int_t pixel=0; pixel<hsignal->GetNbinsX(); pixel++)
445	{
446	// User information
447	d->SetStatusLine1(Form("Fitting pixel #%d", pixel));
448	d->SetProgressBarPosition((pixel+1.)/hsignal->GetNbinsX(), 1);
449
450	// Skip pixels known to be faulty
451	if (usePixel[pixel]==0)
452	continue;
453
454	//Projectipon of a certain Pixel out of Ampl.Spectrum
455	TH1D *hist = hsignal->ProjectionY("proj", pixel+1, pixel+1);
456	hist->SetDirectory(0);
457
458	if (hist->GetEntries()<100)
459	{
460	gLog << warn << pixel << " ...histogram empty." << endl;
461	usePixel[pixel] = 0;
462	delete hist;
463	continue;
464	}
465
466	//Rebin Projection
467	hist->Rebin(2);
468
469	// Fit range
470	hist->GetXaxis()->SetRangeUser(150, hist->GetXaxis()->GetXmax());
471
472	// Determine start values
473	const Int_t maxBin = hist->GetMaximumBin();
474	const Double_t maxPos = hist->GetBinCenter(maxBin);
475
476	const Double_t gain = res_par[1];
477	const Double_t GainRMS = res_par[2];
478
479	const double fit_min = maxPos-GainRMSgain2.5;
480	const double fit_max = fitmax;//maxPos+gain*(maxOrder-0.5);
481
482	TArrayD cpy_par(7, res_par);
483
484	cpy_par[0] = hist->GetBinContent(maxBin);
485	cpy_par[1] = maxPos-res_par[4]; // correct position for avg baseline
486
487	func.SetParameters(cpy_par.GetArray());
488	func.SetParLimits(0, 0, 2*cpy_par[0]);
489	func.SetParLimits(1, 0, 2*cpy_par[1]);
490
491	// For individual spectra, the average fit yields 1 anyway
492	//func.SetParameter(6, 0); // Expo
493
494	// ----------- Fit Pixels spectrum ---------------
495
496	const TFitResultPtr rc = hist->Fit(&func, /fast?/"LLN0QS"/:"LLIN0QS"/, "", fit_min, fit_max);
497
498	// ----------- Calculate quality parameter ---------------
499
500	Int_t b1 = hist->GetXaxis()->FindFixBin(fit_min);
501	Int_t b2 = hist->GetXaxis()->FindFixBin(fit_max);
502
503	Double_t chi2 = 0;
504	Int_t cnt = 0;
505	for (int i=b1; i<=b2; i++)
506	{
507	if (hist->GetBinContent(i)<1.5 \|\| func.Eval(hist->GetBinCenter(i))<1.5)
508	continue;
509
510	const Double_t y = func.Integral(hist->GetBinLowEdge(i), hist->GetBinLowEdge(i+1));
511	const Double_t v = hist->GetBinContent(i)*hist->GetBinWidth(i);
512
513	const Double_t chi = (v-y)/v;
514
515	chi2 += chi*chi;
516	cnt ++;
517	}
518
519	chi2 = cnt==0 ? 0 : sqrt(chi2/cnt);
520
521	// ----------------- Fit result --------------------
522
523	const double fit_prob = rc->Prob();
524
525	const float fRate = integral(func, hist)/(ext[pixel]0.5);
526	const float fGain = func.GetParameter(1);
527	const float fGainRMS = func.GetParameter(2);
528	const float fCrosstalkP= func.GetParameter(3);
529	const float fCrosstlk = xtalk(func);
530	const float fOffset = func.GetParameter(4);
531	const float fNoise = func.GetParameter(5)<0 ? fGainRMS*fGain/sqrt(integration_window) : func.GetParameter(5)/sqrt(integration_window);
532	const float fCoeffR = func.GetParameter(6);
533
534	// Fill histograms with result values
535	cRate.SetBinContent( pixel+1, fRate);
536	cGain.SetBinContent( pixel+1, fGain);
537	cRelSigma.SetBinContent( pixel+1, fGainRMS);
538	cCrosstalk.SetBinContent( pixel+1, fCrosstlk);
539	cBaseline.SetBinContent( pixel+1, fOffset/integration_window);
540	cNoise.SetBinContent( pixel+1, fNoise);
541	cChi2.SetBinContent( pixel+1, chi2);
542	cNormGain.SetBinContent( pixel+1, fGain/gain);
543	cFitProb.SetBinContent( pixel+1, fit_prob);
544	cCrosstalkP.SetBinContent(pixel+1, fCrosstalkP);
545	cCoeffR.SetBinContent( pixel+1, fCoeffR);
546
547	// ======================================================
548
549	// Try to determine faulty fits
550
551	bool ok = int(rc)==0;
552
553	// mark pixels suspicious with failed fit
554	if (!ok)
555	gLog << warn << pixel << " ...fit failed!" << endl;
556
557	// mark pixels suspicious with negative GainRMS
558	if (fabs(fGain/gain-1)>0.3)
559	{
560	gLog << warn << pixel << " ...gain deviates more than 30% from sum-gain." << endl;
561	ok = 0;
562	}
563
564	if (fabs(fOffset/integration_window)>3)
565	{
566	gLog << warn << pixel << " ...baseline deviates." << endl;
567	ok = 0;
568	}
569
570	// cancel out pixel where the fit was not succsessfull
571	usePixel[pixel] = ok;
572
573	// Plot pixel 0 and 5 (TM) and all faulty fits
574	if (pixel==0 \|\| pixel==5 \|\| !ok)
575	{
576	TCanvas &c = d->AddTab(Form("Pix%d", pixel));
577	c.cd();
578	gPad->SetLogy();
579	gPad->SetGridx();
580	gPad->SetGridy();
581
582	hist->SetStats(false);
583	hist->SetXTitle("Extracted signal");
584	hist->SetYTitle("Counts");
585	hist->SetName(Form("Pix%d", pixel));
586	hist->GetXaxis()->SetRange();
587	hist->DrawCopy("hist")->SetDirectory(0);
588	func.DrawCopy("SAME")->SetRange(fit_min, fit_max);
589
590	cout << "--------------------" << endl;
591	cout << "Pixel: " << pixel << endl;
592	cout << "fit prob: " << fit_prob << endl;
593	cout << "AmplEst: " << cpy_par[0] << endl;
594	cout << "GainEst: " << cpy_par[1] << endl;
595	PrintFunc(func, integration_window);
596	cout << "--------------------" << endl;
597	}
598
599	if (!ok)
600	{
601	delete hist;
602	continue;
603	}
604
605	// Fill Parameters into histograms
606	hRate1.Fill( fRate);
607	hGain1.Fill( fGain);
608	hRelSigma1.Fill( fGainRMS);
609	hCrosstalk1.Fill( fCrosstlk);
610	hBaseline1.Fill( fOffset/integration_window);
611	hNoise1.Fill( fNoise);
612	hChiSq1.Fill( chi2);
613	hNormGain1.Fill( fGain/gain);
614	hFitProb1.Fill( fit_prob);
615	hCrosstalkP1.Fill(fCrosstalkP);
616	hCoeffR1.Fill( fCoeffR);
617
618	if (pmap.index(pixel).pixel()!=8)
619	{
620	hRate2.Fill( fRate);
621	hGain2.Fill( fGain);
622	hRelSigma2.Fill( fGainRMS);
623	hCrosstalk2.Fill( fCrosstlk);
624	hBaseline2.Fill( fOffset/integration_window);
625	hNoise2.Fill( fNoise);
626	hChiSq2.Fill( chi2);
627	hNormGain2.Fill( fGain/gain);
628	hFitProb2.Fill( fit_prob);
629	hCrosstalkP2.Fill(fCrosstalkP);
630	hCoeffR2.Fill( fCoeffR);
631	}
632
633	// Fill sum spectrum
634	for (int b=1; b<=hist->GetNbinsX(); b++)
635	hSumScale1.Fill((hist->GetBinCenter(b)-fOffset)/fGain, hist->GetBinContent(b));
636
637	if (pmap.index(pixel).pixel()!=8)
638	for (int b=1; b<=hist->GetNbinsX(); b++)
639	hSumScale2.Fill((hist->GetBinCenter(b)-fOffset)/fGain, hist->GetBinContent(b));
640
641	delete hist;
642
643	// Because of the rebinning...
644	hist = hsignal->ProjectionY("proj", pixel+1, pixel+1);
645	hSumClear1.Add(hist);
646	if (pmap.index(pixel).pixel()!=8)
647	hSumClear2.Add(hist);
648	delete hist;
649
650	hist = htime->ProjectionY("proj", pixel+1, pixel+1);
651	hTime.Add(hist);
652	delete hist;
653
654	hist = hpulse->ProjectionY("proj", pixel+1, pixel+1);
655	hPulse.Add(hist);
656	delete hist;
657
658	count_ok++;
659	}
660
661	//------------------fill histograms-----------------------
662	// Display only pixels used and with valid fits
663
664	cRate.SetUsed(usePixel);
665	cGain.SetUsed(usePixel);
666	cRelSigma.SetUsed(usePixel);
667	cCrosstalk.SetUsed(usePixel);
668	cBaseline.SetUsed(usePixel);
669	cNoise.SetUsed(usePixel);
670	cChi2.SetUsed(usePixel);
671	cNormGain.SetUsed(usePixel);
672	cFitProb.SetUsed(usePixel);
673	cCrosstalkP.SetUsed(usePixel);
674	cCoeffR.SetUsed(usePixel);
675
676	// --------------------------------------------------------
677	// Display data
678
679	TCanvas *canv = &d->AddTab("Cams1");
680	canv->Divide(3,2);
681
682	canv->cd(1);
683	cRate.DrawCopy();
684
685	canv->cd(2);
686	cGain.DrawCopy();
687
688	canv->cd(3);
689	cBaseline.DrawCopy();
690
691	canv->cd(4);
692	cRelSigma.DrawCopy();
693
694	canv->cd(5);
695	cCrosstalk.DrawCopy();
696
697	canv->cd(6);
698	cNoise.DrawCopy();
699
700
701	canv = &d->AddTab("Cams2");
702	canv->Divide(3,2);
703
704	canv->cd(1);
705	cFitProb.DrawCopy();
706
707	canv->cd(2);
708	cChi2.DrawCopy();
709
710	canv->cd(4);
711	cCoeffR.DrawCopy();
712
713	canv->cd(5);
714	cCrosstalkP.DrawCopy();
715
716	// --------------------------------------------------------
717
718	gStyle->SetOptFit(1);
719
720	canv = &d->AddTab("Hists1");
721	canv->Divide(3,2);
722
723	TH1 *hh = 0;
724
725	canv->cd(1);
726	hh = hRate1.DrawCopy();
727	hh = hRate2.DrawCopy("same");
728
729	canv->cd(2);
730	hh = hGain1.DrawCopy();
731	hh = hGain2.DrawCopy("same");
732	hh->Fit("gaus");
733
734	canv->cd(3);
735	hh = hBaseline1.DrawCopy();
736	hh = hBaseline2.DrawCopy("same");
737	hh->Fit("gaus");
738
739	canv->cd(4);
740	hh = hRelSigma1.DrawCopy();
741	hh = hRelSigma2.DrawCopy("same");
742	hh->Fit("gaus");
743
744	canv->cd(5);
745	hh = hCrosstalk1.DrawCopy();
746	hh = hCrosstalk2.DrawCopy("same");
747	hh->Fit("gaus");
748
749	canv->cd(6);
750	hh = hNoise1.DrawCopy();
751	hh = hNoise2.DrawCopy("same");
752	hh->Fit("gaus");
753
754	// --------------------------------------------------------
755
756	canv = &d->AddTab("Hists2");
757	canv->Divide(3,2);
758
759	canv->cd(1);
760	gPad->SetLogy();
761	hh = hFitProb1.DrawCopy();
762	hh = hFitProb2.DrawCopy("same");
763	hh->Fit("gaus");
764
765	canv->cd(2);
766	hChiSq1.DrawCopy();
767	hChiSq2.DrawCopy("same");
768
769	canv->cd(4);
770	hh = hCoeffR1.DrawCopy();
771	hh = hCoeffR2.DrawCopy("same");
772	hh->Fit("gaus");
773
774	canv->cd(5);
775	hh = hCrosstalkP1.DrawCopy();
776	hh = hCrosstalkP2.DrawCopy("same");
777	hh->Fit("gaus");
778
779	// --------------------------------------------------------
780
781	canv = &d->AddTab("NormGain");
782	canv->Divide(2,1);
783
784	canv->cd(1);
785	cNormGain.SetMinimum(0.8);
786	cNormGain.SetMaximum(1.2);
787	cNormGain.DrawCopy();
788
789	canv->cd(2);
790	gPad->SetLogy();
791	hh = hNormGain1.DrawCopy();
792	hh = hNormGain2.DrawCopy("same");
793	hh->Fit("gaus");
794
795	//------------------ Draw gain corrected sum specetrum -------------------
796	gROOT->SetSelectedPad(0);
797	c11.cd();
798	hSumClear1.DrawCopy("hist same");
799
800	//-------------------- fit gain corrected sum spectrum -------------------
801
802	gROOT->SetSelectedPad(0);
803	TCanvas &c11b = d->AddTab("CleanHist1");
804	c11b.cd();
805	gPad->SetLogy();
806	gPad->SetGridx();
807	gPad->SetGridy();
808
809	const Int_t maxbin1 = hSumClear1.GetMaximumBin();
810	const Double_t ampl1 = hSumClear1.GetBinContent(maxbin1);
811
812	func.SetParameters(res_par);
813	func.SetParLimits(0, 0, 2*ampl1);
814	func.SetParameter(0, ampl1);
815	func.ReleaseParameter(6);
816
817	hSumClear1.Fit(&func, fast?"LN0QSR":"LIN0QSR");
818
819	hSumClear1.DrawCopy();
820	func.DrawCopy("same");
821
822	cout << "--------------------" << endl;
823	PrintFunc(func, integration_window);
824	cout << "--------------------" << endl;
825
826	//-------------------- fit gain corrected sum spectrum -------------------
827
828	gROOT->SetSelectedPad(0);
829	TCanvas &c11c = d->AddTab("CleanHist2");
830	c11c.cd();
831	gPad->SetLogy();
832	gPad->SetGridx();
833	gPad->SetGridy();
834
835	const Int_t maxbin1b = hSumClear2.GetMaximumBin();
836	const Double_t ampl1b = hSumClear2.GetBinContent(maxbin1b);
837
838	func.SetParameters(res_par);
839	func.SetParLimits(0, 0, 2*ampl1b);
840	func.SetParameter(0, ampl1b);
841	func.ReleaseParameter(6);
842
843	hSumClear2.Fit(&func, fast?"LN0QSR":"LIN0QSR");
844
845	hSumClear2.DrawCopy();
846	func.DrawCopy("same");
847
848	cout << "--------------------" << endl;
849	PrintFunc(func, integration_window);
850	cout << "--------------------" << endl;
851
852	//-------------------- fit gain corrected sum spectrum -------------------
853
854	gROOT->SetSelectedPad(0);
855	TCanvas &c12 = d->AddTab("GainHist1");
856	c12.cd();
857	gPad->SetLogy();
858	gPad->SetGridx();
859	gPad->SetGridy();
860
861	const Int_t maxbin2 = hSumScale1.GetMaximumBin();
862	const Double_t ampl2 = hSumScale1.GetBinContent(maxbin2);
863
864	//Set fit parameters
865	Double_t par2[7] =
866	{
867	ampl2, 1, 0.1, res_par[3], 0, res_par[5]<0 ? -1 : res_par[5]/res_par[1], res_par[6]
868	};
869
870	func.SetParameters(par2);
871	func.SetParLimits(0, 0, 2*ampl2);
872	func.FixParameter(1, 1);
873	func.FixParameter(4, 0);
874
875	func.SetRange(0.62, 9);
876	hSumScale1.Fit(&func, fast?"LN0QSR":"LIN0QSR");
877
878	hSumScale1.DrawCopy();
879	func.DrawCopy("same");
880
881	cout << "--------------------" << endl;
882	PrintFunc(func, integration_window);
883	cout << "--------------------" << endl;
884
885	//-------------------- fit gain corrected sum spectrum -------------------
886
887	gROOT->SetSelectedPad(0);
888	TCanvas &c12b = d->AddTab("GainHist2");
889	c12b.cd();
890	gPad->SetLogy();
891	gPad->SetGridx();
892	gPad->SetGridy();
893
894	const Int_t maxbin2b = hSumScale2.GetMaximumBin();
895	const Double_t ampl2b = hSumScale2.GetBinContent(maxbin2b);
896
897	//Set fit parameters
898	Double_t par2b[7] =
899	{
900	ampl2b, 1, 0.1, res_par[3], 0, res_par[5]<0 ? -1 : res_par[5]/res_par[1], res_par[6]
901	};
902
903	func.SetParameters(par2b);
904	func.SetParLimits(0, 0, 2*ampl2b);
905	func.FixParameter(1, 1);
906	func.FixParameter(4, 0);
907
908	func.SetRange(0.62, 9);
909	hSumScale2.Fit(&func, fast?"LN0QSR":"LIN0QSR");
910
911	hSumScale2.DrawCopy();
912	func.DrawCopy("same");
913
914	cout << "--------------------" << endl;
915	PrintFunc(func, integration_window);
916	cout << "--------------------" << endl;
917
918	//--------fit gausses to peaks of gain corrected sum specetrum -----------
919
920	d->AddTab("ArrTime");
921	gPad->SetGrid();
922	hTime.DrawCopy();
923
924	// -----------------------------------------------------------------
925
926	d->AddTab("Pulse");
927	gPad->SetGrid();
928	hPulse.DrawCopy();
929
930	// ================================================================
931
932	cout << "Saving results to '" << outfile << "'" << endl;
933	d->SaveAs(outfile);
934	cout << "..success!" << endl;
935
936	TFile f(outfile, "UPDATE");
937	par.Write();
938	ext.Write("ExtractionRange");
939	header.Write();
940
941	return 0;
942	}

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

Context Navigation

source: trunk/Mars/fact/analysis/gain/fit_spectra.C@ 20115

Download in other formats: