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

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