Context Navigation

source: branches/Corsika7405Compatibility/fact/analysis/gain/fit_spectra.C@ 19365

Visit:

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

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

Download in other formats: