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

Last change on this file since 17033 was 17033, checked in by tbretz, 12 years ago
First version of files.
File size: 23.1 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 "MStatusArray.h"
12#include "MStatusDisplay.h"
13#include "MHCamera.h"
14#include "MGeomCamFACT.h"
15#include "MParameters.h"
16#include "MArrayI.h"
17
18using namespace std;
19
20// --------------------------------------------------------------------------
21
22// Fit function for a single pe spectrum
23Double_t fcn_g(Double_t *xx, Double_t *par)
24{
25 const Double_t ampl = par[0];
26 const Double_t gain = par[1];
27 const Double_t sigma = par[2]*gain;
28 const Double_t cross = par[3];
29 const Double_t shift = par[4];
30 const Double_t noise = par[5];
31 const Double_t expo = par[6];
32
33 Double_t y = 0;
34 for (int N=1; N<14; N++)
35 {
36 const Double_t muN = N*gain + shift;
37 const Double_t sigN = TMath::Sqrt(N*sigma*sigma + noise*noise);
38
39 const Double_t p = TMath::Power(cross, N-1) * TMath::Power(N, -expo);
40
41 y += TMath::Gaus(xx[0], muN, sigN) * p / sigN;
42 }
43
44 const Double_t sig1 = TMath::Sqrt(sigma*sigma + noise*noise);
45 return ampl*sig1*y;
46}
47
48// Calculate the crosstalk from the function parameters
49Double_t xtalk(TF1 &f)
50{
51 Double_t cross = f.GetParameter(3);
52 Double_t expo = f.GetParameter(6);
53
54 Double_t y = 0;
55 for (int N=2; N<14; N++)
56 y += TMath::Power(cross, N-1) * TMath::Power(N, -expo);
57
58 return y / (y + 1);
59}
60
61// calculate the integral in units per millisecond
62double integral(TF1 &func, TH1 &hist)
63{
64 const Double_t sigma = func.GetParameter(2)*func.GetParameter(1);
65 const Double_t cross = func.GetParameter(3);
66 const Double_t noise = func.GetParameter(5);
67 const Double_t expo = func.GetParameter(6);
68
69 Double_t sum = 0;
70 for (int N=1; N<14; N++)
71 sum += TMath::Power(cross, N-1) * TMath::Power(N, -expo);
72
73 const Double_t scale = hist.GetBinWidth(1);
74 const Double_t sig1 = TMath::Sqrt(sigma*sigma + noise*noise);
75 const Double_t integ = func.GetParameter(0)*sum*sig1*sqrt(TMath::TwoPi())/scale;
76
77 return integ/1e-9/1e6;
78}
79
80
81// --------------------------------------------------------------------------
82
83// Print function parameters
84void PrintFunc(TF1 &f, float integration_window=30)
85{
86 //cout << "--------------------" << endl;
87 cout << "Ampl: " << setw(8) << f.GetParameter(0) << " +/- " << f.GetParError(0) << endl;
88 cout << "Gain: " << setw(8) << f.GetParameter(1) << " +/- " << f.GetParError(1) << endl;
89 cout << "Rel.sigma: " << setw(8) << f.GetParameter(2) << " +/- " << f.GetParError(2) << endl;
90 cout << "Baseline: " << setw(8) << f.GetParameter(4)/integration_window << " +/- " << f.GetParError(4)/integration_window << endl;
91 cout << "Crosstalk: " << setw(8) << f.GetParameter(3) << " +/- " << f.GetParError(3) << endl;
92 cout << "Pcrosstalk: " << setw(8) << xtalk(f) << endl;
93 cout << "Noise: " << setw(8) << f.GetParameter(5)/sqrt(integration_window) << " +/- " << f.GetParError(5)/sqrt(integration_window) << endl;
94 cout << "Expo: " << setw(8) << f.GetParameter(6) << " +/- " << f.GetParError(6) << endl;
95 //cout << "--------------------" << endl;
96
97}
98
99// --------------------------------------------------------------------------
100
101// The parameters for the function are the filenam from the gain extraction
102// and the output filename
103int fit_spectra(const char *filename = "20130222_018_018.root",
104 const char *outfile = "20120802_247_247-fit.root" )
105{
106 // It is more correct to fit the integral, but this is super
107 // slow, especially for 1440 pixel. To get that, one would have
108 // to analytically integrate and fit this function.
109 // (Currently the integral is switched off for the 1440 individual
110 // spectra in all cases)
111 bool fast = false; // Switch off using integral
112
113 // Map for which pixel shall be plotted and which not
114 TArrayC usePixel(1440);
115 memset(usePixel.GetArray(), 1, 1440);
116
117 // List of Pixel that should be ignored in camera view
118 usePixel[424] = 0;
119 usePixel[923] = 0;
120 usePixel[1208] = 0;
121 usePixel[583] = 0;
122 usePixel[830] = 0;
123 usePixel[1399] = 0;
124 usePixel[113] = 0;
125 usePixel[115] = 0;
126 usePixel[354] = 0;
127 usePixel[423] = 0;
128 usePixel[1195] = 0;
129 usePixel[1393] = 0;
130
131 cout << setprecision(3);
132
133 // ======================================================
134 // Read data and histograms from file
135
136 TFile file(filename);
137 if (file.IsZombie())
138 {
139 cout << "Opening file '" << filename << "' failed." << endl;
140 return 1;
141 }
142
143 MStatusArray arr;
144 if (arr.Read()<=0)
145 {
146 cout << "Reading of MStatusArray from '" << filename << "' failed." << endl;
147 return 2;
148 }
149
150 TH2 *hsignal = (TH2*)arr.FindObjectInCanvas("Signal", "TH2F", "MHSingles");
151 if (!hsignal)
152 {
153 cout << "Histogram Signal not found in '" << filename << "'." << endl;
154 return 3;
155 }
156 TH2 *htime = (TH2*)arr.FindObjectInCanvas("Time", "TH2F", "MHSingles");
157 if (!htime)
158 {
159 cout << "Histogram Time not found in '" << filename << "'." << endl;
160 return 3;
161 }
162 TProfile2D *hpulse = (TProfile2D*)arr.FindObjectInCanvas("Pulse", "TProfile2D", "MHSingles");
163 if (!hpulse)
164 {
165 cout << "Histogram Pulse not found in '" << filename << "'." << endl;
166 return 3;
167 }
168 TH2F *hbase = (TH2F*)arr.FindObjectInCanvas("Baseline", "TH2F", "MHBaseline");
169 if (!hbase)
170 {
171 cout << "Histogram Baseline not found in '" << filename << "'." << endl;
172 return 3;
173 }
174
175 MParameterI par("NumEvents");
176 if (par.Read()<=0)
177 {
178 cout << "NumEvents not found in '" << filename << "'." << endl;
179 return 4;
180 }
181
182 MParameterI win("IntegrationWindow");
183 if (win.Read()<=0)
184 {
185 cout << "IntegrationWindow not found in '" << filename << "'." << endl;
186 return 5;
187 }
188
189 Int_t integration_window = win.GetVal();
190
191 MArrayI ext;
192 if (ext.Read("ExtractionRange")<=0)
193 {
194 cout << "ExtractionRange not found in '" << filename << "'." << endl;
195 return 6;
196
197 }
198
199 // ======================================================
200
201 MStatusDisplay *d = new MStatusDisplay;
202
203 // Camera geometry for displays
204 MGeomCamFACT fact;
205
206 // ------------------ Spectrum Fit ---------------
207 // Instantiate the display histograms
208 MHCamera cRate(fact);
209 MHCamera cGain(fact);
210 MHCamera cRelSigma(fact);
211 MHCamera cCrosstalk(fact);
212 MHCamera cBaseline(fact);
213 MHCamera cNoise(fact);
214 MHCamera cChi2(fact);
215 MHCamera cNormGain(fact);
216 MHCamera cFitProb(fact);
217
218 // Set name and title for the histograms
219 cRate.SetNameTitle ("Rate", "Dark count rate");
220 cGain.SetNameTitle ("Gain", "Gain distribution");
221 cRelSigma.SetNameTitle ("RelSigma", "Rel. Sigma");
222 cCrosstalk.SetNameTitle("Crosstalk", "Crosstalk probability");
223 cBaseline.SetNameTitle ("Baseline", "Baseline per sample");
224 cNoise.SetNameTitle ("Noise", "Noise per sample");
225 cChi2.SetNameTitle ("Chi2", "\\Chi^2");
226 cNormGain.SetNameTitle ("NormGain", "Normalized gain");
227 cFitProb.SetNameTitle ("FitProb", "Root's fit probability");
228
229 // Instantiate 1D histograms for the distributions
230 TH1F hRate ("Rate", "Dark count rate", 200, 0, 10);
231 TH1F hGain ("Gain", "Gain distribution", 100, 0, 400);
232 TH1F hRelSigma ("RelSigma", "Rel. Sigma", 160, 0, 0.40);
233 TH1F hCrosstalk("Crosstalk", "Crosstalk probability", 120, 0, 0.30);
234 TH1F hBaseline ("Baseline", "Baseline per sample", 75, -7.5, 7.5);
235 TH1F hNoise ("Noise", "Noise per sample", 60, 0, 30);
236 TH1F hChi2 ("Chi2", "\\Chi^2", 200, 0, 4);
237 TH1F hNormGain ("NormGain", "Normalized gain", 51, 0.5, 1.5);
238 TH1F hFitProb ("FitProb", "FitProb distribution", 100, 0, 1);
239
240 // Histigram for the sum of all spectrums
241 TH1F hSum("Sum1", "Signal spectrum of all pixels",
242 hsignal->GetNbinsY(), hsignal->GetYaxis()->GetXmin(), hsignal->GetYaxis()->GetXmax());
243 hSum.SetXTitle("pulse integral [mV sample]");
244 hSum.SetYTitle("Counts");
245 hSum.SetStats(false);
246 hSum.Sumw2();
247
248 // Histogram for the sum of all pixels excluding the ones with faulty fits
249 TH1F hSumClear("SumC", "Signal spectrum of all pixels",
250 hsignal->GetNbinsY(), hsignal->GetYaxis()->GetXmin(), hsignal->GetYaxis()->GetXmax());
251 hSumClear.SetXTitle("pulse integral [mV sample]");
252 hSumClear.SetYTitle("Counts");
253 hSumClear.SetStats(false);
254 hSumClear.SetLineColor(kBlue);
255 hSumClear.Sumw2();
256
257 // Arrival time spectrum of the extracted pulses
258 TH1F hTime("Time", "Arrival time spectrum", htime->GetNbinsY(), htime->GetYaxis()->GetXmin(), htime->GetYaxis()->GetXmax());
259 hTime.SetXTitle("pulse arrival time [sample]");
260 hTime.SetYTitle("Counts");
261 hTime.SetStats(false);
262
263 // average pulse shape of the extracted pulses
264 TH1F hPulse("Puls", "Average pulse", hpulse->GetNbinsY(), hpulse->GetYaxis()->GetXmin(), hpulse->GetYaxis()->GetXmax());
265 hPulse.SetXTitle("pulse arrival time [sample]");
266 hPulse.SetYTitle("Counts");
267 hPulse.SetStats(false);
268
269 // Spektrum for the normalized individual spectra
270 TH1F hSumScale("Sum2", "Signal spectrum of all pixels", 120, 0.05, 12.05);
271 hSumScale.SetXTitle("pulse integral [pe]");
272 hSumScale.SetYTitle("Rate");
273 hSumScale.SetStats(false);
274 hSumScale.Sumw2();
275
276 // define fit function for Amplitudespectrum
277 TF1 func("spektrum", fcn_g, 0, hSum.GetXaxis()->GetXmax(), 7);
278 func.SetNpx(2000);
279 func.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset", "Noise", "Expo");
280 func.SetLineColor(kRed);
281
282 //--------------------- fill sum spectrum --------------------------------
283
284 d->SetStatusLine("Calculating sum spectrum", 0);
285
286 // Begin of Loop over Pixels
287 for (Int_t pixel = 0; pixel < hsignal->GetNbinsX(); pixel++)
288 {
289 //jump to next pixel if the current is marked as faulty
290 if (usePixel[pixel]==0)
291 continue;
292
293 TH1D *hist = hsignal->ProjectionY("proj", pixel+1, pixel+1);
294 hSum.Add(hist);
295 delete hist;
296 }
297
298 //----------------- get starting values -------------------------------
299
300 hSum.GetXaxis()->SetRangeUser(150, hSum.GetXaxis()->GetXmax());
301
302 const Int_t maxbin = hSum.GetMaximumBin();
303 const Double_t maxpos = hSum.GetBinCenter(maxbin);
304 const Double_t binwidth = hSum.GetBinWidth(maxbin);
305 const Double_t ampl = hSum.GetBinContent(maxbin);
306
307 double fwhmSum = 0;
308
309 //Calculate full width half Maximum
310 for (int i=1; i<maxbin; i++)
311 if (hSum.GetBinContent(maxbin-i)+hSum.GetBinContent(maxbin+i)<ampl)
312 {
313 fwhmSum = 2*(i-0.5)*binwidth;
314 break;
315 }
316
317 if (fwhmSum==0)
318 {
319 cout << "Could not determine start value for sigma." << endl;
320 }
321
322 Double_t sigma_est = fwhmSum/2.3548; // FWHM = 2*sqrt(2*ln(2))*sigma
323
324 Double_t fitmin = maxpos-3*sigma_est; // was 3
325 Double_t fitmax = hSum.GetXaxis()->GetXmax();
326
327 // ------------------- fit --------------------------------
328
329 //Fit and draw spectrum
330 func.SetParLimits(0, 0, 2*ampl); // Amplitude
331 func.SetParLimits(1, 0, 2*maxpos); // Gain
332 func.SetParLimits(2, 0, 1); // Sigma
333 func.SetParLimits(3, 0, 1); // Crosstalk
334 func.SetParLimits(5, 0, 150); // Noise
335
336 func.SetParameter(0, ampl); // Amplitude
337 func.SetParameter(1, maxpos); // Gain
338 func.SetParameter(2, 0.1); // Sigma
339 func.SetParameter(3, 0.15); // Crosstalk
340 func.SetParameter(4, 0*integration_window); // Baseline
341 func.SetParameter(5, 1.*sqrt(integration_window)); // Noise
342 func.SetParameter(6, 0.4); // Expo
343
344 func.SetRange(fitmin, fitmax);
345 hSum.Fit(&func, fast?"N0QSR":"IN0QSR");
346
347 Double_t res_par[7];
348 func.GetParameters(res_par);
349
350 // ------------------ display result -------------------------------
351
352 cout << "--------------------" << endl;
353 cout << "AmplEst: " << ampl << endl;
354 cout << "GainEst: " << maxpos << endl;
355 cout << "SigmaEst: " << sigma_est << endl;
356 PrintFunc(func, integration_window);
357 cout << "--------------------" << endl;
358
359 gROOT->SetSelectedPad(0);
360 TCanvas &c11 = d->AddTab("SumHist");
361 c11.cd();
362 gPad->SetLogy();
363 gPad->SetGridx();
364 gPad->SetGridy();
365 hSum.GetXaxis()->SetRange();
366 hSum.DrawCopy("hist");
367 func.DrawCopy("same");
368
369 // ===================================================================
370 // Gain Calculation for each Pixel
371 // ===================================================================
372
373 // counter for number of processed pixel
374 int count_ok = 0;
375
376 // Begin of Loop over Pixels
377 for (Int_t pixel=0; pixel<hsignal->GetNbinsX(); pixel++)
378 {
379 // User information
380 d->SetStatusLine(Form("Fitting pixel #%d", pixel), 0);
381 d->SetProgressBarPosition((pixel+1.)/hsignal->GetNbinsX(), 1);
382
383 // Skip pixels known to be faulty
384 if (usePixel[pixel]==0)
385 continue;
386
387 //Projectipon of a certain Pixel out of Ampl.Spectrum
388 TH1D *hist = hsignal->ProjectionY("proj", pixel+1, pixel+1);
389 hist->SetDirectory(0);
390
391 if (hist->GetEntries()<100)
392 {
393 cout << pixel << " ...histogram empty." << endl;
394 usePixel[pixel] = 0;
395 delete hist;
396 continue;
397 }
398
399 //Rebin Projection
400 hist->Rebin(2);
401
402 // Fit range
403 hist->GetXaxis()->SetRangeUser(150, hist->GetXaxis()->GetXmax());
404
405 // Determine start values
406 const Int_t maxBin = hist->GetMaximumBin();
407 const Double_t maxPos = hist->GetBinCenter(maxBin);
408
409 const Double_t gain = res_par[1];
410 const Double_t GainRMS = res_par[2];
411
412 const double fit_min = maxPos-GainRMS*gain*2.5;
413 const double fit_max = fitmax;//maxPos+gain*(maxOrder-0.5);
414
415 TArrayD cpy_par(7, res_par);
416
417 cpy_par[0] = hist->GetBinContent(maxBin);
418 cpy_par[1] = maxPos-res_par[4]; // correct position for avg baseline
419
420 func.SetParameters(cpy_par.GetArray());
421 func.SetParLimits(0, 0, 2*cpy_par[0]);
422 func.SetParLimits(1, 0, 2*cpy_par[1]);
423
424 // For individual spectra, the average fit yields 1 anyway
425 func.FixParameter(6, 1); // Expo
426
427 // ----------- Fit Pixels spectrum ---------------
428
429 const TFitResultPtr rc = hist->Fit(&func, /*fast?*/"LLN0QS"/*:"LLIN0QS"*/, "", fit_min, fit_max);
430
431 // ----------- Calculate quality parameter ---------------
432
433 Int_t b1 = hist->GetXaxis()->FindFixBin(fit_min);
434 Int_t b2 = hist->GetXaxis()->FindFixBin(fit_max);
435
436 Double_t chi2 = 0;
437 Int_t cnt = 0;
438 for (int i=b1; i<=b2; i++)
439 {
440 if (hist->GetBinContent(i)<1.5 || func.Eval(hist->GetBinCenter(i))<1.5)
441 continue;
442
443 const Double_t y = func.Integral(hist->GetBinLowEdge(i), hist->GetBinLowEdge(i+1));
444 const Double_t v = hist->GetBinContent(i)*hist->GetBinWidth(i);
445
446 const Double_t chi = (v-y)/v;
447
448 chi2 += chi*chi;
449 cnt ++;
450 }
451
452 chi2 = cnt==0 ? 0 : sqrt(chi2/cnt);
453
454 // ----------------- Fit result --------------------
455
456 const double fit_prob = rc->Prob();
457
458 const float fRate = integral(func, *hist)/(ext[pixel]*0.5);
459 const float fGain = func.GetParameter(1);
460 const float fGainRMS = func.GetParameter(2);
461 const float fCrosstlk = xtalk(func);
462 const float fOffset = func.GetParameter(4);
463 const float fNoise = func.GetParameter(5)/sqrt(integration_window);
464
465 // Fill histograms with result values
466 cRate.SetBinContent( pixel+1, fRate);
467 cGain.SetBinContent( pixel+1, fGain);
468 cRelSigma.SetBinContent( pixel+1, fGainRMS);
469 cCrosstalk.SetBinContent( pixel+1, fCrosstlk);
470 cBaseline.SetBinContent( pixel+1, fOffset/integration_window);
471 cNoise.SetBinContent( pixel+1, fNoise);
472 cChi2.SetBinContent( pixel+1, chi2);
473 cNormGain.SetBinContent( pixel+1, fGain/gain);
474 cFitProb.SetBinContent( pixel+1, fit_prob);
475
476 // ======================================================
477
478 // Try to determine faulty fits
479
480 bool ok = int(rc)==0;
481
482 // mark pixels suspicious with failed fit
483 if (!ok)
484 cout << pixel << " ...fit failed!" << endl;
485
486 // mark pixels suspicious with negative GainRMS
487 if (fabs(fGain/gain-1)>0.3)
488 {
489 cout << pixel << " ...gain deviates more than 30% from sum-gain." << endl;
490 ok = 0;
491 }
492
493 if (fabs(fOffset/integration_window)>3)
494 {
495 cout << pixel << " ...baseline deviates." << endl;
496 ok = 0;
497 }
498
499 // cancel out pixel where the fit was not succsessfull
500 usePixel[pixel] = ok;
501
502 // Plot pixel 0 and all faulty fits
503 if (pixel==0 || !ok)
504 {
505 TCanvas &c = d->AddTab(Form("Pix%d", pixel));
506 c.cd();
507 gPad->SetLogy();
508 gPad->SetGridx();
509 gPad->SetGridy();
510
511 hist->SetStats(false);
512 hist->SetXTitle("Extracted signal");
513 hist->SetYTitle("Counts");
514 hist->SetName(Form("Pix%d", pixel));
515 hist->GetXaxis()->SetRange();
516 hist->DrawCopy("hist")->SetDirectory(0);
517 func.DrawCopy("SAME")->SetRange(fit_min, fit_max);
518
519 cout << "--------------------" << endl;
520 cout << "Pixel: " << pixel << endl;
521 cout << "fit prob: " << fit_prob << endl;
522 cout << "AmplEst: " << cpy_par[0] << endl;
523 cout << "GainEst: " << cpy_par[1] << endl;
524 PrintFunc(func, integration_window);
525 cout << "--------------------" << endl;
526 }
527
528 if (!ok)
529 {
530 delete hist;
531 continue;
532 }
533
534 // Fill Parameters into histograms
535 hRate.Fill( fRate);
536 hGain.Fill( fGain);
537 hRelSigma.Fill( fGainRMS);
538 hCrosstalk.Fill( fCrosstlk);
539 hBaseline.Fill( fOffset/integration_window);
540 hNoise.Fill( fNoise);
541 hChi2.Fill( chi2);
542 hNormGain.Fill( fGain/gain);
543 hFitProb.Fill( fit_prob);
544
545 // Fill sum spectrum
546 for (int b=1; b<=hist->GetNbinsX(); b++)
547 hSumScale.Fill((hist->GetBinCenter(b)-fOffset)/fGain, hist->GetBinContent(b));
548 delete hist;
549
550 // Because of the rebinning...
551 hist = hsignal->ProjectionY("proj", pixel+1, pixel+1);
552 hSumClear.Add(hist);
553 delete hist;
554
555 hist = htime->ProjectionY("proj", pixel+1, pixel+1);
556 hTime.Add(hist);
557 delete hist;
558
559 hist = hpulse->ProjectionY("proj", pixel+1, pixel+1);
560 hPulse.Add(hist);
561 delete hist;
562
563 count_ok++;
564 }
565
566 //------------------fill histograms-----------------------
567 // Display only pixels used and with valid fits
568
569 cRate.SetUsed(usePixel);
570 cGain.SetUsed(usePixel);
571 cRelSigma.SetUsed(usePixel);
572 cCrosstalk.SetUsed(usePixel);
573 cBaseline.SetUsed(usePixel);
574 cNoise.SetUsed(usePixel);
575 cChi2.SetUsed(usePixel);
576 cNormGain.SetUsed(usePixel);
577 cFitProb.SetUsed(usePixel);
578
579 // --------------------------------------------------------
580 // Display data
581
582 TCanvas *canv = &d->AddTab("Cams");
583 canv->Divide(4,2);
584
585 canv->cd(1);
586 cRate.DrawCopy();
587
588 canv->cd(2);
589 cGain.DrawCopy();
590
591 canv->cd(3);
592 cRelSigma.DrawCopy();
593
594 canv->cd(4);
595 cCrosstalk.DrawCopy();
596
597 canv->cd(5);
598 cBaseline.DrawCopy();
599
600 canv->cd(6);
601 cNoise.DrawCopy();
602
603 canv->cd(7);
604 cFitProb.DrawCopy();
605
606 canv->cd(8);
607 cChi2.DrawCopy();
608
609 // --------------------------------------------------------
610
611 gStyle->SetOptFit(1);
612
613 canv = &d->AddTab("Hists");
614 canv->Divide(4,2);
615
616 TH1 *hh = 0;
617
618 canv->cd(1);
619 hh = hRate.DrawCopy();
620
621 canv->cd(2);
622 hh = hGain.DrawCopy();
623 hh->Fit("gaus");
624
625 canv->cd(3);
626 hh = hRelSigma.DrawCopy();
627 hh->Fit("gaus");
628
629 canv->cd(4);
630 hh = hCrosstalk.DrawCopy();
631 hh->Fit("gaus");
632
633 canv->cd(5);
634 hh = hBaseline.DrawCopy();
635 hh->Fit("gaus");
636
637 canv->cd(6);
638 hh = hNoise.DrawCopy();
639 hh->Fit("gaus");
640
641 canv->cd(7);
642 gPad->SetLogy();
643 hFitProb.DrawCopy();
644
645 canv->cd(8);
646 hChi2.DrawCopy();
647
648 // --------------------------------------------------------
649
650 canv = &d->AddTab("NormGain");
651 canv->Divide(2,1);
652
653 canv->cd(1);
654 cNormGain.SetMinimum(0.8);
655 cNormGain.SetMaximum(1.2);
656 cNormGain.DrawCopy();
657
658 canv->cd(2);
659 gPad->SetLogy();
660 hh = hNormGain.DrawCopy();
661 hh->Fit("gaus");
662
663 //------------------ Draw gain corrected sum specetrum -------------------
664 gROOT->SetSelectedPad(0);
665 c11.cd();
666 hSumClear.DrawCopy("hist same");
667
668 //-------------------- fit gain corrected sum spectrum -------------------
669
670 gROOT->SetSelectedPad(0);
671 TCanvas &c11b = d->AddTab("CleanHist");
672 c11b.cd();
673 gPad->SetLogy();
674 gPad->SetGridx();
675 gPad->SetGridy();
676
677 const Int_t maxbin1 = hSumClear.GetMaximumBin();
678 const Double_t ampl1 = hSumClear.GetBinContent(maxbin1);
679
680 func.SetParameters(res_par);
681 func.SetParLimits(0, 0, 2*ampl1);
682 func.SetParameter(0, ampl1);
683 func.ReleaseParameter(6);
684
685 hSumClear.Fit(&func, fast?"LN0QSR":"LIN0QSR");
686
687 hSumClear.DrawCopy();
688 func.DrawCopy("same");
689
690 cout << "--------------------" << endl;
691 PrintFunc(func, integration_window);
692 cout << "--------------------" << endl;
693
694 //-------------------- fit gain corrected sum spectrum -------------------
695
696 gROOT->SetSelectedPad(0);
697 TCanvas &c12 = d->AddTab("GainHist");
698 c12.cd();
699 gPad->SetLogy();
700 gPad->SetGridx();
701 gPad->SetGridy();
702
703 const Int_t maxbin2 = hSumScale.GetMaximumBin();
704 const Double_t ampl2 = hSumScale.GetBinContent(maxbin2);
705
706 //Set fit parameters
707 Double_t par2[7] =
708 {
709 ampl2, 1, 0.1, res_par[3], 0, res_par[5]/res_par[1], res_par[6]
710 };
711
712 func.SetParameters(par2);
713 func.SetParLimits(0, 0, 2*ampl2);
714 func.FixParameter(1, 1);
715 func.FixParameter(4, 0);
716
717 func.SetRange(0.62, 9);
718 hSumScale.Fit(&func, fast?"LN0QSR":"LIN0QSR");
719
720 hSumScale.DrawCopy();
721 func.DrawCopy("same");
722
723 cout << "--------------------" << endl;
724 PrintFunc(func, integration_window);
725 cout << "--------------------" << endl;
726
727 //--------fit gausses to peaks of gain corrected sum specetrum -----------
728
729 d->AddTab("ArrTime");
730 gPad->SetGrid();
731 hTime.DrawCopy();
732
733 // -----------------------------------------------------------------
734
735 d->AddTab("Pulse");
736 gPad->SetGrid();
737 hPulse.DrawCopy();
738
739 // ================================================================
740
741 cout << "saving results to rootfile" << endl;
742 d->SaveAs(outfile);
743 cout << "..success!" << endl;
744
745 TFile f(outfile, "UPDATE");
746 par.Write();
747 ext.Write("ExtractionRange");
748
749 return 0;
750}
751
Note: See TracBrowser for help on using the repository browser.