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

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