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

Last change on this file since 14181 was 14181, checked in by Jens Buss, 12 years ago
fitt to sum spectrum added
  • Property svn:executable set to *
File size: 22.8 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
739 double fwhm2 = 0;
740 for (int i=1; i<maxbin; i++)
741 if (hSum1.GetBinContent(maxbin-i)+hSum1.GetBinContent(maxbin+i)<ampl*2)
742 {
743 fwhm2 = (2*i+1)*binwidth;
744 break;
745 }
746
747 if (fwhm2==0)
748 {
749 cout << "Could not determine start value for sigma." << endl;
750 }
751
752
753 Double_t par[5] =
754 {
755 ampl, hSum1.GetBinCenter(maxbin), fwhm2*1.4, 0.1, -10
756 };
757
758 func2.SetParameters(par);
759
760 const double min = par[1]-fwhm2*1.4;
761 const double max = par[1]*5;
762
763 hSum1.Fit(&func2, "NOQS", "", min, max);
764 func2.DrawCopy("SAME");
765 //-------------------END OF: fit sum spectrum-----------------------------
766
767
768 TCanvas &c12 = d->AddTab("GainHist");
769 gPad->SetLogy();
770 hSum2.DrawCopy();
771 const Double_t fMaxPos = hSum2.GetBinCenter(maxbin);
772
773 //------------------fit gausses to peaks -----------------------
774
775 UInt_t nfunc = 5;
776 TF1 **gaus = new TF1*[nfunc];
777
778 for (UInt_t nr = 0; nr<nfunc; nr++)
779 {
780 char fname[20];
781 sprintf(fname,"gaus%d",nr);
782 gaus[nr] = new TF1(fname,"gaus", fMaxPos*(nr+1)-0.6, fMaxPos*(nr+1)+0.6);
783 hSum2.Fit(gaus[nr], "N0QS", "", fMaxPos*(nr+1)-0.6, fMaxPos*(nr+1)+0.6);
784 gaus[nr]->DrawCopy("SAME");
785// delete gaus[nr];
786 }
787
788
789 //------------------fit gain corrected sum spectrum-----------------------
790
791 const Int_t maxbin2 = hSum2.GetMaximumBin();
792 const Double_t binwidth2 = hSum2.GetBinWidth(maxbin);
793 const Double_t ampl2 = hSum2.GetBinContent(maxbin);
794 cout << "AMPL: " << ampl2 << endl;
795
796 fwhm2 = 0;
797 for (int i=1; i<maxbin; i++)
798 if (hSum2.GetBinContent(maxbin2-i)+hSum2.GetBinContent(maxbin2+i)<ampl2*2)
799 {
800 fwhm2 = (2*i+1)*binwidth2;
801 break;
802 }
803
804 if (fwhm2==0)
805 {
806 cout << "Could not determine start value for sigma." << endl;
807 }
808
809
810 Double_t par2[5] =
811 {
812 ampl2, hSum2.GetBinCenter(maxbin2), fwhm2*1.4, 0.1, 0
813 };
814
815 cout << "AMPL: " << ampl2 << endl;
816 cout << "maxpos: " << hSum2.GetBinCenter(maxbin2) << endl;
817 cout << "fwhm: " << fwhm2*1.4 << endl;
818
819 func3.SetParameters(par2);
820
821 const double min2 = par2[1]-fwhm2*1.4;
822 const double max2 = par2[1]*5;
823
824 hSum2.Fit(&func3, "N0QS", "", min2, max2);
825 func3.DrawCopy("SAME");
826 //-------------------END OF: fit sum spectrum-----------------------------
827
828 TCanvas &c2 = d->AddTab("Result");
829 c2.Divide(3,2);
830 c2.cd(1);
831 gPad->SetLogy();
832 hAmplitude.DrawCopy();
833
834 c2.cd(2);
835 gPad->SetLogy();
836 hGain.DrawCopy();
837
838 c2.cd(3);
839 gPad->SetLogy();
840 hGainRMS.DrawCopy();
841
842 c2.cd(4);
843 gPad->SetLogy();
844 hCrosstlk.DrawCopy();
845
846 c2.cd(5);
847 gPad->SetLogy();
848 hOffset.DrawCopy();
849
850 c2.cd(6);
851 gPad->SetLogy();
852 hFitProb.DrawCopy();
853
854 /*
855 TCanvas &c3 = d->AddTab("Separation");
856 gPad->SetLogy();
857 hSep.DrawCopy();
858 */
859 return 0;
860}
861
862Double_t fcn(Double_t *xx, Double_t *par)
863{
864 Double_t x = xx[0];
865
866 Double_t ampl = par[0];
867 Double_t gain = par[1];
868 Double_t sigma = par[2];
869 Double_t cross = par[3];
870 Double_t shift = par[4];
871
872 Double_t xTalk = 1;
873
874 Double_t y = 0;
875 for (int order = 1; order < 14; order++)
876 {
877 Double_t arg = (x - order*gain - shift)/sigma;
878
879 Double_t gauss = TMath::Exp(-0.5 * arg * arg/order);
880
881 y += xTalk*gauss;
882
883 xTalk *= cross;
884 }
885
886 return ampl*y;
887}
Note: See TracBrowser for help on using the repository browser.