Context Navigation

source: trunk/Mars/hawc/fit_singles.C@ 19887

Visit:

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

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

Download in other formats: