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

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