source: fact/tools/marsmacros/singlepe.C@ 14226

Last change on this file since 14226 was 14182, checked in by Jens Buss, 14 years ago
add array of gaussions to fit to spectrum, calculation pulseheights in spectrum to calculate crosstalk and underlying poisson, plot spektrum height vs pe
  • Property svn:executable set to *
File size: 23.5 KB
Line 
1#include <TStyle.h>
2#include <TCanvas.h>
3#include <TSystem.h>
4#include <TF1.h>
5#include <TProfile.h>
6#include <TProfile2D.h>
7#include <TMath.h>
8#include <TFitResultPtr.h>
9#include <TFitResult.h>
10
11#include "MH.h"
12#include "MLog.h"
13#include "MLogManip.h"
14#include "MDirIter.h"
15#include "MFillH.h"
16#include "MEvtLoop.h"
17#include "MCamEvent.h"
18#include "MGeomApply.h"
19#include "MTaskList.h"
20#include "MParList.h"
21#include "MContinue.h"
22#include "MBinning.h"
23#include "MHDrsCalib.h"
24#include "MDrsCalibApply.h"
25#include "MRawFitsRead.h"
26#include "MBadPixelsCam.h"
27#include "MStatusDisplay.h"
28#include "MTaskInteractive.h"
29#include "MPedestalSubtractedEvt.h"
30#include "MHCamera.h"
31#include "MGeomCamFACT.h"
32#include "MRawRunHeader.h"
33
34using namespace std;
35
36MPedestalSubtractedEvt *fEvent = 0;
37MLog *fLog = &gLog;
38
39struct Single
40{
41 float fSignal;
42 float fTime;
43};
44
45class MSingles : public MParContainer, public MCamEvent
46{
47 Int_t fIntegrationWindow;
48
49 vector<vector<Single>> fData;
50
51public:
52 MSingles(const char *name=NULL, const char *title=NULL) : fIntegrationWindow(0)
53 {
54 fName = name ? name : "MSingles";
55 }
56
57 void InitSize(const UInt_t i)
58 {
59 fData.resize(i);
60 }
61
62 vector<Single> &operator[](UInt_t i) { return fData[i]; }
63 vector<Single> &GetVector(UInt_t i) { return fData[i]; }
64
65 UInt_t GetNumPixels() const { return fData.size(); }
66
67 void SetIntegrationWindow(Int_t w) { fIntegrationWindow = w; }
68 Int_t GetIntegrationWindow() const { return fIntegrationWindow; }
69
70 Bool_t GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const
71 {
72 return kTRUE;
73 }
74 void DrawPixelContent(Int_t num) const { }
75
76 ClassDef(MSingles, 1)
77};
78
79class MH2F : public TH2F
80{
81public:
82 Int_t Fill(Double_t x,Double_t y)
83 {
84 Int_t binx, biny, bin;
85 fEntries++;
86 binx = TMath::Nint(x+1);
87 biny = fYaxis.FindBin(y);
88 bin = biny*(fXaxis.GetNbins()+2) + binx;
89 AddBinContent(bin);
90 if (fSumw2.fN) ++fSumw2.fArray[bin];
91 if (biny == 0 || biny > fYaxis.GetNbins()) {
92 if (!fgStatOverflows) return -1;
93 }
94 ++fTsumw;
95 ++fTsumw2;
96 fTsumwy += y;
97 fTsumwy2 += y*y;
98 return bin;
99 }
100 ClassDef(MH2F, 1);
101};
102
103class MProfile2D : public TProfile2D
104{
105public:
106 Int_t Fill(Double_t x, Double_t y, Double_t z)
107 {
108 Int_t bin,binx,biny;
109 fEntries++;
110 binx =TMath::Nint(x+1);
111 biny =fYaxis.FindBin(y);
112 bin = GetBin(binx, biny);
113 fArray[bin] += z;
114 fSumw2.fArray[bin] += z*z;
115 fBinEntries.fArray[bin] += 1;
116 if (fBinSumw2.fN) fBinSumw2.fArray[bin] += 1;
117 if (biny == 0 || biny > fYaxis.GetNbins()) {
118 if (!fgStatOverflows) return -1;
119 }
120 ++fTsumw;
121 ++fTsumw2;
122 fTsumwy += y;
123 fTsumwy2 += y*y;
124 fTsumwz += z;
125 fTsumwz2 += z*z;
126 return bin;
127 }
128 ClassDef(MProfile2D, 1);
129};
130
131
132
133class MHSingles : public MH
134{
135 MH2F fSignal;
136 MH2F fTime;
137 MProfile2D fPulse;
138
139 UInt_t fNumEntries;
140
141 MSingles *fSingles; //!
142 MPedestalSubtractedEvt *fEvent; //!
143
144public:
145 MHSingles() : fNumEntries(0), fSingles(0), fEvent(0)
146 {
147 fName = "MHSingles";
148
149 fSignal.SetName("Signal");
150 fSignal.SetTitle("Title");
151 fSignal.SetXTitle("X title");
152 fSignal.SetYTitle("Y title");
153 fSignal.SetDirectory(NULL);
154
155 fTime.SetName("Time");
156 fTime.SetTitle("Title");
157 fTime.SetXTitle("X title");
158 fTime.SetYTitle("Y title");
159 fTime.SetDirectory(NULL);
160
161 fPulse.SetName("Pulse");
162 fPulse.SetTitle("Title");
163 fPulse.SetXTitle("X title");
164 fPulse.SetYTitle("Y title");
165 fPulse.SetDirectory(NULL);
166 }
167
168 Bool_t SetupFill(const MParList *plist)
169 {
170 fSingles = (MSingles*)plist->FindObject("MSingles");
171 if (!fSingles)
172 {
173 *fLog << /* err << */ "MSingles not found... abort." << endl;
174 return kFALSE;
175 }
176
177 fEvent = (MPedestalSubtractedEvt*)plist->FindObject("MPedestalSubtractedEvt");
178 if (!fEvent)
179 {
180 *fLog << /* err << */ "MPedestalSubtractedEvt not found... abort." << endl;
181 return kFALSE;
182 }
183
184 fNumEntries = 0;
185
186 return kTRUE;
187 }
188 Bool_t ReInit(MParList *plist)
189 {
190 const MRawRunHeader *header = (MRawRunHeader*)plist->FindObject("MRawRunHeader");
191 if (!header)
192 {
193 *fLog << /* err << */ "MRawRunHeader not found... abort." << endl;
194 return kFALSE;
195 }
196
197 const Int_t w = fSingles->GetIntegrationWindow();
198
199 MBinning binsx, binsy;
200 binsx.SetEdges(fSingles->GetNumPixels(), -0.5, fSingles->GetNumPixels()-0.5);
201 binsy.SetEdges(2*150, -2*7.5*w, 2*(60-7.5)*w);
202
203 MH::SetBinning(fSignal, binsx, binsy);
204
205 const UShort_t roi = header->GetNumSamples();
206
207 // Correct for DRS timing!!
208 MBinning binst(roi, -0.5, roi-0.5);
209 MH::SetBinning(fTime, binsx, binst);
210
211 MBinning binspy(2*roi, -roi-0.5, roi-0.5);
212 MH::SetBinning(fPulse, binsx, binspy);
213
214 return kTRUE;
215 }
216
217 Int_t Fill(const MParContainer *par, const Stat_t weight=1)
218 {
219 const Float_t *ptr = fEvent->GetSamples(0);
220 const Short_t roi = fEvent->GetNumSamples();
221
222 for (unsigned int i=0; i<fSingles->GetNumPixels(); i++)
223 {
224 const vector<Single> &result = fSingles->GetVector(i);
225
226 for (unsigned int j=0; j<result.size(); j++)
227 {
228 fSignal.Fill(i, result[j].fSignal);
229 fTime.Fill (i, result[j].fTime);
230
231 if (!ptr)
232 continue;
233
234 const Short_t tm = floor(result[j].fTime);
235
236 for (int s=0; s<roi; s++)
237 fPulse.Fill(i, s-tm, ptr[s]);
238 }
239
240 ptr += roi;
241 }
242
243 fNumEntries++;
244
245 return kTRUE;
246 }
247
248 TH2 *GetSignal() { return &fSignal; }
249 TH2 *GetTime() { return &fTime; }
250 TH2 *GetPulse() { return &fPulse; }
251
252 UInt_t GetNumEntries() const { return fNumEntries; }
253
254 void Draw(Option_t *opt)
255 {
256 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
257
258 AppendPad("");
259
260 pad->Divide(2,2);
261
262 pad->cd(1);
263 gPad->SetBorderMode(0);
264 gPad->SetFrameBorderMode(0);
265 fSignal.Draw("colz");
266
267 pad->cd(2);
268 gPad->SetBorderMode(0);
269 gPad->SetFrameBorderMode(0);
270 fTime.Draw("colz");
271
272 pad->cd(3);
273 gPad->SetBorderMode(0);
274 gPad->SetFrameBorderMode(0);
275 fPulse.Draw("colz");
276 }
277
278 ClassDef(MHSingles, 1)
279};
280
281MStatusDisplay *d = new MStatusDisplay;
282
283MSingles *fSingles;
284
285TH1F *fluct1 = new TH1F("", "", 200, -200, 200);
286TH1F *fluct2 = new TH1F("", "", 100, -50, 50);
287
288TGraph weights("weights.txt");
289
290Int_t PreProcess(MParList *plist)
291{
292 fSingles = (MSingles*)plist->FindCreateObj("MSingles");
293 if (!fSingles)
294 return kFALSE;
295
296 fEvent = (MPedestalSubtractedEvt*)plist->FindObject("MPedestalSubtractedEvt");
297 if (!fEvent)
298 {
299 *fLog << /* err << */ "MPedestalSubtractedEvt not found... abort." << endl;
300 return kFALSE;
301 }
302
303 d->AddTab("Fluct");
304 fluct2->Draw();
305 fluct1->Draw("same");
306
307 fSingles->SetIntegrationWindow(20);
308
309 //if (weights.GetN()==0)
310 // return kFALSE;
311
312 return kTRUE;
313}
314
315Int_t Process()
316{
317 const UInt_t roi = fEvent->GetNumSamples();
318
319 const size_t navg = 10;
320 const float threshold = 5;
321 const UInt_t start_skip = 20;
322 const UInt_t end_skip = 10;
323 const UInt_t integration_size = 10*3;
324 const UInt_t max_search_window = 30;
325
326 vector<float> val(roi-navg);//result of first sliding average
327 for (int pix=0; pix<fEvent->GetNumPixels(); pix++)
328 {
329 vector<Single> &result = fSingles->GetVector(pix);
330 result.clear();
331
332 const Float_t *ptr = fEvent->GetSamples(pix);
333
334 //Sliding window
335 for (size_t i=0; i<roi-navg; i++)
336 {
337 val[i] = 0;
338 for (size_t j=i; j<i+navg; j++)
339 val[i] += ptr[j];
340 val[i] /= navg;
341 }
342
343 //peak finding via threshold
344 for (UInt_t i=start_skip; i<roi-navg-end_skip-30; i++)
345 {
346 //search for threshold crossings
347 if (val[i+0]>threshold ||
348 val[i+4]>threshold ||
349
350 val[i+5]<threshold ||
351 val[i+9]<threshold)
352 continue;
353
354 //search for maximum after threshold crossing
355 UInt_t k_max = i+5;
356 for (UInt_t k=i+5; k<i+max_search_window; k++)
357 {
358 if (val[k] > val[k_max])
359 k_max = k;
360 }
361
362 if (k_max == i+5 || k_max == i + max_search_window-1) continue;
363
364 //search for half maximum before maximum
365 UInt_t k_half_max = 0;
366 for (UInt_t k=k_max; k>k_max-25; k--)
367 {
368 if (val[k-1] < val[k_max]/2 &&
369 val[k] >= val[k_max]/2 )
370 {
371 k_half_max = k;
372 break;
373 }
374 }
375 if (k_half_max > roi-navg-end_skip-max_search_window - integration_size)
376 continue;
377 if (k_half_max == 0)
378 continue;
379 if (k_max - k_half_max > 14)
380 continue;
381
382 fluct2->Fill(k_max - k_half_max);
383
384 // Evaluate arrival time more precisely!!!
385 // Get a better integration window
386
387 Single single;
388 single.fSignal = 0;
389 single.fTime = k_half_max;
390
391
392 // Crossing of the threshold is at 5
393 for (UInt_t j=k_half_max; j<k_half_max+integration_size; j++)
394 single.fSignal += ptr[j];
395 result.push_back(single);
396
397 i += 5+integration_size;
398 }
399 }
400 return kTRUE;
401}
402
403Int_t PostProcess()
404{
405 return kTRUE;
406}
407
408Double_t fcn(Double_t *x, Double_t *par);
409
410int singlepe()
411{
412 // ======================================================
413
414 const char *drsfile = "/fact/raw/2012/06/01/20120601_003.drs.fits.gz";
415
416 MDirIter iter;
417 iter.AddDirectory("/fact/raw/2012/06/01", "20120601_005.fits.gz");
418
419 // ======================================================
420
421 // true: Display correctly mapped pixels in the camera displays
422 // but the value-vs-index plot is in software/spiral indices
423 // false: Display pixels in hardware/linear indices,
424 // but the order is the camera display is distorted
425 bool usemap = true;
426
427 // map file to use (get that from La Palma!)
428 const char *map = usemap ? "/scratch_nfs/FACTmap111030.txt" : NULL;
429
430 // ------------------------------------------------------
431
432 Long_t max1 = 0;
433
434 // ======================================================
435
436 if (map && gSystem->AccessPathName(map, kFileExists))
437 {
438 gLog << "ERROR - Cannot access mapping file '" << map << "'" << endl;
439 return 1;
440 }
441 if (gSystem->AccessPathName(drsfile, kFileExists))
442 {
443 gLog << "ERROR - Cannot access drsfile file '" << drsfile << "'" << endl;
444 return 2;
445 }
446
447 // --------------------------------------------------------------------------------
448
449 gLog.Separator("Singles");
450 gLog << "Extract singles with '" << drsfile << "'" << endl;
451 gLog << endl;
452
453 // ------------------------------------------------------
454
455 MDrsCalibration drscalib300;
456 if (!drscalib300.ReadFits(drsfile))
457 return 4;
458
459 // ------------------------------------------------------
460
461 iter.Sort();
462 iter.Print();
463
464 // ======================================================
465
466 //MStatusDisplay *d = new MStatusDisplay;
467
468 MBadPixelsCam badpixels;
469 badpixels.InitSize(1440);
470 badpixels[ 424].SetUnsuitable(MBadPixelsPix::kUnsuitable);
471 badpixels[ 583].SetUnsuitable(MBadPixelsPix::kUnsuitable);
472 badpixels[ 830].SetUnsuitable(MBadPixelsPix::kUnsuitable);
473 badpixels[ 923].SetUnsuitable(MBadPixelsPix::kUnsuitable);
474 badpixels[1208].SetUnsuitable(MBadPixelsPix::kUnsuitable);
475 badpixels[1399].SetUnsuitable(MBadPixelsPix::kUnsuitable);
476
477 // Twin pixel
478 // 113
479 // 115
480 // 354
481 // 423
482 // 1195
483 // 1393
484
485 //MDrsCalibrationTime timecam;
486
487 MH::SetPalette();
488
489 // ======================================================
490
491 gLog << endl;
492 gLog.Separator("Processing external light pulser run");
493
494 MTaskList tlist1;
495
496 MSingles singles;
497 MRawRunHeader header;
498
499 MParList plist1;
500 plist1.AddToList(&tlist1);
501 plist1.AddToList(&drscalib300);
502 plist1.AddToList(&badpixels);
503 plist1.AddToList(&singles);
504 plist1.AddToList(&header);
505
506 MEvtLoop loop1("DetermineCalConst");
507 loop1.SetDisplay(d);
508 loop1.SetParList(&plist1);
509
510 // ------------------ Setup the tasks ---------------
511
512 MRawFitsRead read1;
513 read1.LoadMap(map);
514 read1.AddFiles(iter);
515
516 MContinue cont1("(MRawEvtHeader.GetTriggerID&0xff00)!=0x100", "SelectCal");
517
518 MGeomApply apply1;
519
520 MDrsCalibApply drsapply1;
521
522 MTaskInteractive mytask;
523 mytask.SetPreProcess(PreProcess);
524 mytask.SetProcess(Process);
525 mytask.SetPostProcess(PostProcess);
526
527 MFillH fill("MHSingles");
528
529 // ------------------ Setup eventloop and run analysis ---------------
530
531 tlist1.AddToList(&read1);
532// tlist1.AddToList(&cont1);
533 tlist1.AddToList(&apply1);
534 tlist1.AddToList(&drsapply1);
535 tlist1.AddToList(&mytask);
536 tlist1.AddToList(&fill);
537
538 if (!loop1.Eventloop(max1))
539 return 10;
540
541 if (!loop1.GetDisplay())
542 return 11;
543
544 gStyle->SetOptFit(kTRUE);
545
546 MHSingles* h = (MHSingles*) plist1.FindObject("MHSingles");
547 if (!h)
548 return 0;
549
550 TH2 *hsignal = h->GetSignal();
551 TH2 *htime = h->GetTime();
552 TH2 *hpulse = h->GetPulse();
553
554 d->AddTab("Time");
555 TH1 *ht = htime->ProjectionY();
556 ht->Scale(1./1440);
557 ht->Draw();
558
559 d->AddTab("Pulse");
560 TH1 *hp = hpulse->ProjectionY();
561 hp->Scale(1./1440);
562 hp->Draw();
563
564
565 // ------------------ Spectrum Fit ---------------
566 // AFTER this the Code for the spektrum fit follows
567 // access the spektrumhistogram by the pointer *hist
568
569 MGeomCamFACT fact;
570 MHCamera hcam(fact);
571
572 //built an array of Amplitude spectra
573 TH1F hAmplitude("Rate", "Average number of dark counts per event",
574 100, 0, 10);
575 TH1F hGain ("Gain", "Gain distribution",
576 50, 200, 350);
577 TH1F hGainRMS ("RelSigma", "Distribution of rel. Sigma",
578 100, 0, 30);
579 TH1F hCrosstlk ("Crosstalk", "Distribution of crosstalk probability",
580 100, 0, 25);
581 TH1F hOffset ("Offset", "Baseline offset distribution",
582 50, -7.5, 2.5);
583 TH1F hFitProb ("FitProb", "FitProb distribution",
584 50, 0, 100);
585 TH1F hSum1 ("Sum1", "Sum",
586 100, 0, 2000);
587 TH1F hSum2 ("Sum2", "Sum",
588 100, 0, 10);
589
590 // define fit function for Amplitudespectrum
591 TF1 func("spektrum", fcn, 0, 1000, 5);
592 func.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset");
593 func.SetLineColor(kBlue);
594
595 // define fit function for Amplitudespectrum
596 TF1 func2("sum_spektrum", fcn, 0, 2000, 5);
597 func2.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset");
598 func2.SetLineColor(kBlue);
599
600 // define fit function for Amplitudespectrum
601 TF1 func3("gain_sum_spektrum", fcn, 0, 10, 5);
602 func3.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset");
603 func3.SetLineColor(kBlue);
604
605
606
607 // Map for which pixel shall be plotted and which not
608 TArrayC usePixel(1440);
609 memset(usePixel.GetArray(), 1, 1440);
610
611 // List of Pixel that should be ignored in camera view
612 usePixel[424] = 0;
613 usePixel[923] = 0;
614 usePixel[1208] = 0;
615 usePixel[583] = 0;
616 usePixel[830] = 0;
617 usePixel[1399] = 0;
618 usePixel[113] = 0;
619 usePixel[115] = 0;
620 usePixel[354] = 0;
621 usePixel[423] = 0;
622 usePixel[1195] = 0;
623 usePixel[1393] = 0;
624
625 int count_ok = 0;
626
627 // Begin of Loop over Pixels
628 for (UInt_t pixel = 0; pixel < singles.GetNumPixels(); pixel++)
629 {
630 if (usePixel[pixel]==0)
631 continue;
632
633 //Projectipon of a certain Pixel out of Ampl.Spectrum
634 TH1D *hist = hsignal->ProjectionY("proj", pixel+1, pixel+1);
635 hist->SetDirectory(0);
636
637 const Int_t maxbin = hist->GetMaximumBin();
638 const Double_t binwidth = hist->GetBinWidth(maxbin);
639 const Double_t ampl = hist->GetBinContent(maxbin);
640
641 double fwhm = 0;
642 for (int i=1; i<maxbin; i++)
643 if (hist->GetBinContent(maxbin-i)+hist->GetBinContent(maxbin+i)<ampl*2)
644 {
645 fwhm = (2*i+1)*binwidth;
646 break;
647 }
648
649 if (fwhm==0)
650 {
651 cout << "Could not determine start value for sigma." << endl;
652 continue;
653 }
654
655 // Amplitude, Gain, Sigma, Crosstalk Prob, Shift
656 // par[1] + par[4] is the position of the first maximum
657 Double_t par[5] =
658 {
659 ampl, hist->GetBinCenter(maxbin), fwhm*1.4, 0.1, -10
660 };
661
662 func.SetParameters(par);
663
664 const double min = par[1]-fwhm*1.4;
665 const double max = par[1]*5;
666
667 //Fit Pixels spectrum
668 const TFitResultPtr rc = hist->Fit(&func, "N0QS", "", min, max);
669
670 const bool ok = int(rc)==0;
671
672 const double fit_prob = rc->Prob();
673 const float fGain = func.GetParameter(1);
674 const float fAmplitude = func.Integral(0, 7*fGain)/h->GetNumEntries(); //func.GetParameter(0);
675 const float f2GainRMS = func.GetParameter(2);
676 const float fCrosstlk = func.GetParameter(3);
677 const float fOffset = func.GetParameter(4)/singles.GetIntegrationWindow();
678
679 hAmplitude.Fill(fAmplitude);
680 hGain.Fill(fGain);
681 //hGainRMS.Fill(f2GainRMS);
682 hCrosstlk.Fill(fCrosstlk*100);
683 hOffset.Fill(fOffset);
684 hFitProb.Fill(100*fit_prob);
685 hGainRMS.Fill(100*f2GainRMS/fGain);
686
687 hcam.SetBinContent(pixel+1, fGain);
688
689 usePixel[pixel] = ok;
690
691 if (!ok)
692 {
693 cout << "Pixel " << setw(4) << pixel << ": failed!" << endl;
694 continue;
695 }
696
697 for (int b=1; b<=hist->GetNbinsX(); b++)
698 {
699 hSum1.Fill(hist->GetBinCenter(b), hist->GetBinContent(b));
700 hSum2.Fill((hist->GetBinCenter(b)-fOffset)/fGain, hist->GetBinContent(b));
701 }
702
703 if (count_ok==0)
704 {
705 TCanvas &c = d->AddTab(Form("Pix%d", pixel));
706 c.cd();
707 hist->SetName(Form("Pix%d", pixel));
708 hist->DrawCopy()->SetDirectory(0);
709 func.DrawCopy("SAME");
710 }
711
712
713
714 count_ok++;
715
716 delete hist;
717 }
718
719 //------------------fill histograms-----------------------
720
721 hcam.SetUsed(usePixel);
722
723 TCanvas &canv = d->AddTab("Gain");
724 canv.cd();
725 hcam.SetNameTitle( "gain","Gain distribution over whole camera");
726 hcam.DrawCopy();
727
728 gROOT->SetSelectedPad(0);
729 TCanvas &c11 = d->AddTab("SumHist");
730 gPad->SetLogy();
731 hSum1.DrawCopy();
732
733 //------------------------fit sum spectrum-------------------------------
734 const Int_t maxbin = hSum1.GetMaximumBin();
735 const Double_t binwidth = hSum1.GetBinWidth(maxbin);
736 const Double_t ampl = hSum1.GetBinContent(maxbin);
737
738 double fwhm2 = 0;
739
740 //Calculate full width half Maximum
741 for (int i=1; i<maxbin; i++)
742 if (hSum1.GetBinContent(maxbin-i)+hSum1.GetBinContent(maxbin+i)<ampl*2)
743 {
744 fwhm2 = (2*i+1)*binwidth;
745 break;
746 }
747 if (fwhm2==0)
748 {
749 cout << "Could not determine start value for sigma." << endl;
750 }
751
752 //Set fit parameters
753 Double_t par[5] =
754 {
755 ampl, hSum1.GetBinCenter(maxbin), fwhm2*1.4, 0.1, -10
756 };
757 func2.SetParameters(par);
758
759 //calculate fit range
760 const double min = par[1]-fwhm2*1.4;
761 const double max = par[1]*5;
762
763 //Fit and draw spectrum
764 hSum1.Fit(&func2, "NOQS", "", min, max);
765 func2.DrawCopy("SAME");
766
767 //-------------------END OF: fit sum spectrum-----------------------------
768
769 // Draw gain corrected sum specetrum
770 gROOT->SetSelectedPad(0);
771 TCanvas &c12 = d->AddTab("GainHist");
772 gPad->SetLogy();
773 hSum2.DrawCopy();
774 const Double_t fMaxPos = hSum2.GetBinCenter(maxbin);
775
776 //--------fit gausses to peaks of gain corrected sum specetrum -----------
777 UInt_t nfunc = 5;
778
779 // create array for gaus funtions to fit to peaks of spectrum
780 TF1 **gaus = new TF1*[nfunc];
781
782 // create histo for height of Maximum amplitudes vs. pe
783 TH1D hMaxHeights("MaxHeights", "peak heights",
784 20, 0, 20);
785
786 // fit gauss functions to spectrum
787 for (UInt_t nr = 0; nr<nfunc; nr++)
788 {
789 char fname[20];
790 sprintf(fname,"gaus%d",nr);
791
792 //create gauss functions
793 gaus[nr] = new TF1( fname,"gaus",
794 fMaxPos*(nr+1)-0.6, fMaxPos*(nr+1)+0.6);
795 //fit and draw gauss functions
796 hSum2.Fit( gaus[nr], "N0QS", "",
797 fMaxPos*(nr+1)-0.6, fMaxPos*(nr+1)+0.6);
798 gaus[nr]->DrawCopy("SAME");
799
800 //fill histo for height of Maximum amplitudes vs. pe
801 hMaxHeights.SetBinContent(nr, gaus[nr]->GetParameter(0));
802 delete gaus[nr];
803 }
804 delete[] gaus;
805
806 //------------------fit gain corrected sum spectrum-----------------------
807
808 const Int_t maxbin2 = hSum2.GetMaximumBin();
809 const Double_t binwidth2 = hSum2.GetBinWidth(maxbin);
810 const Double_t ampl2 = hSum2.GetBinContent(maxbin);
811
812 //Calculate full width half Maximum
813 fwhm2 = 0;
814 for (int i=1; i<maxbin; i++)
815 if (hSum2.GetBinContent(maxbin2-i)+hSum2.GetBinContent(maxbin2+i)<ampl2*2)
816 {
817 fwhm2 = (2*i+1)*binwidth2;
818 break;
819 }
820 if (fwhm2==0)
821 {
822 cout << "Could not determine start value for sigma." << endl;
823 }
824
825 //Set fit parameters
826 Double_t par2[5] =
827 {
828 ampl2, hSum2.GetBinCenter(maxbin2), fwhm2*1.4, 0.1, 0
829 };
830
831 func3.SetParameters(par2);
832
833 const double min2 = par2[1]-fwhm2*1.4;
834 const double max2 = par2[1]*5;
835
836 hSum2.Fit(&func3, "N0QS", "", min2, max2);
837 func3.DrawCopy("SAME");
838 //-------------------END OF: fit sum spectrum-----------------------------
839
840 TCanvas &c15 = d->AddTab("PeakHeights");
841 gPad->SetLogy();
842 hMaxHeights.DrawCopy();
843
844 TCanvas &c2 = d->AddTab("Result");
845 c2.Divide(3,2);
846 c2.cd(1);
847 gPad->SetLogy();
848 hAmplitude.DrawCopy();
849
850 c2.cd(2);
851 gPad->SetLogy();
852 hGain.DrawCopy();
853
854 c2.cd(3);
855 gPad->SetLogy();
856 hGainRMS.DrawCopy();
857
858 c2.cd(4);
859 gPad->SetLogy();
860 hCrosstlk.DrawCopy();
861
862 c2.cd(5);
863 gPad->SetLogy();
864 hOffset.DrawCopy();
865
866 c2.cd(6);
867 gPad->SetLogy();
868 hFitProb.DrawCopy();
869
870 /*
871 TCanvas &c3 = d->AddTab("Separation");
872 gPad->SetLogy();
873 hSep.DrawCopy();
874 */
875 return 0;
876}
877
878Double_t fcn(Double_t *xx, Double_t *par)
879{
880 Double_t x = xx[0];
881
882 Double_t ampl = par[0];
883 Double_t gain = par[1];
884 Double_t sigma = par[2];
885 Double_t cross = par[3];
886 Double_t shift = par[4];
887
888 Double_t xTalk = 1;
889
890 Double_t y = 0;
891 for (int order = 1; order < 14; order++)
892 {
893 Double_t arg = (x - order*gain - shift)/sigma;
894
895 Double_t gauss = TMath::Exp(-0.5 * arg * arg/order);
896
897 y += xTalk*gauss;
898
899 xTalk *= cross;
900 }
901
902 return ampl*y;
903}
Note: See TracBrowser for help on using the repository browser.