source: fact/tools/marsmacros/recalcSinglepe.C@ 16858

Last change on this file since 16858 was 14277, checked in by Jens Buss, 12 years ago
Initial Commit: Recalculate gainfits for an allready calculated Amplitude/Integral Spectrum
  • Property svn:executable set to *
File size: 28.6 KB
Line 
1#include <TStyle.h>
2#include <TCanvas.h>
3#include <TSystem.h>
4#include <TF1.h>
5#include <TFile.h>
6#include <TH1.h>
7#include <TH2.h>
8#include <TProfile.h>
9#include <TProfile2D.h>
10#include <TMath.h>
11#include <TFitResultPtr.h>
12#include <TFitResult.h>
13#include <TGProgressBar.h> // TGHProgressBar
14#include <TVirtualFitter.h>
15#include <Getline.h>
16
17#include <cstdio>
18#include <stdio.h>
19#include <stdint.h>
20
21#include "MH.h"
22#include "MLog.h"
23#include "MLogManip.h"
24#include "MDirIter.h"
25#include "MFillH.h"
26#include "MEvtLoop.h"
27#include "MCamEvent.h"
28#include "MGeomApply.h"
29#include "MTaskList.h"
30#include "MParList.h"
31#include "MContinue.h"
32#include "MBinning.h"
33#include "MHDrsCalib.h"
34#include "MStatusArray.h"
35#include "MDrsCalibApply.h"
36#include "MRawFitsRead.h"
37#include "MBadPixelsCam.h"
38#include "MStatusDisplay.h"
39#include "MTaskInteractive.h"
40#include "MPedestalSubtractedEvt.h"
41#include "MHCamera.h"
42#include "MGeomCamFACT.h"
43#include "MRawRunHeader.h"
44
45using namespace std;
46using namespace TMath;
47
48MPedestalSubtractedEvt *fEvent = 0;
49MLog *fLog = &gLog;
50
51struct Single
52{
53 float fSignal;
54 float fTime;
55};
56
57class MSingles : public MParContainer, public MCamEvent
58{
59 Int_t fIntegrationWindow;
60
61 vector<vector<Single>> fData;
62
63public:
64 MSingles(const char *name=NULL, const char *title=NULL) : fIntegrationWindow(0)
65 {
66 fName = name ? name : "MSingles";
67 }
68
69 void InitSize(const UInt_t i)
70 {
71 fData.resize(i);
72 }
73
74 vector<Single> &operator[](UInt_t i) { return fData[i]; }
75 vector<Single> &GetVector(UInt_t i) { return fData[i]; }
76
77 UInt_t GetNumPixels() const { return fData.size(); }
78
79 void SetIntegrationWindow(Int_t w) { fIntegrationWindow = w; }
80 Int_t GetIntegrationWindow() const { return fIntegrationWindow; }
81
82 Bool_t GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const
83 {
84 return kTRUE;
85 }
86 void DrawPixelContent(Int_t num) const { }
87
88 ClassDef(MSingles, 1)
89};
90/*
91class MH2F : public TH2F
92{
93public:
94 Int_t Fill(Double_t x,Double_t y)
95 {
96 Int_t binx, biny, bin;
97 fEntries++;
98 binx = TMath::Nint(x+1);
99 biny = fYaxis.FindBin(y);
100 bin = biny*(fXaxis.GetNbins()+2) + binx;
101 AddBinContent(bin);
102 if (fSumw2.fN) ++fSumw2.fArray[bin];
103 if (biny == 0 || biny > fYaxis.GetNbins()) {
104 if (!fgStatOverflows) return -1;
105 }
106 ++fTsumw;
107 ++fTsumw2;
108 fTsumwy += y;
109 fTsumwy2 += y*y;
110 return bin;
111 }
112 ClassDef(MH2F, 1);
113};
114
115class MProfile2D : public TProfile2D
116{
117public:
118 Int_t Fill(Double_t x, Double_t y, Double_t z)
119 {
120 Int_t bin,binx,biny;
121 fEntries++;
122 binx =TMath::Nint(x+1);
123 biny =fYaxis.FindBin(y);
124 bin = GetBin(binx, biny);
125 fArray[bin] += z;
126 fSumw2.fArray[bin] += z*z;
127 fBinEntries.fArray[bin] += 1;
128 if (fBinSumw2.fN) fBinSumw2.fArray[bin] += 1;
129 if (biny == 0 || biny > fYaxis.GetNbins()) {
130 if (!fgStatOverflows) return -1;
131 }
132 ++fTsumw;
133 ++fTsumw2;
134 fTsumwy += y;
135 fTsumwy2 += y*y;
136 fTsumwz += z;
137 fTsumwz2 += z*z;
138 return bin;
139 }
140 ClassDef(MProfile2D, 1);
141};
142*/
143
144
145class MHSingles : public MH
146{
147 TH2F fSignal;
148 TH2F fTime;
149 TProfile2D fPulse;
150
151 UInt_t fNumEntries;
152
153 MSingles *fSingles; //!
154 MPedestalSubtractedEvt *fEvent; //!
155
156public:
157 MHSingles() : fNumEntries(0), fSingles(0), fEvent(0)
158 {
159 fName = "MHSingles";
160
161 fSignal.SetName("Signal");
162 fSignal.SetTitle("pulse integral spectrum");
163 fSignal.SetXTitle("pixel [SoftID]");
164 fSignal.SetYTitle("time [a.u.]");
165 fSignal.SetDirectory(NULL);
166
167 fTime.SetName("Time");
168 fTime.SetTitle("arival time spectrum");
169 fTime.SetXTitle("pixel [SoftID]");
170 fTime.SetYTitle("time [a.u.]");
171 fTime.SetDirectory(NULL);
172
173 fPulse.SetName("Pulse");
174 fPulse.SetTitle("average pulse shape spectrum");
175 fPulse.SetXTitle("pixel [SoftID]");
176 fPulse.SetYTitle("time [a.u.]");
177 fPulse.SetDirectory(NULL);
178 }
179
180 Bool_t SetupFill(const MParList *plist)
181 {
182 fSingles = (MSingles*)plist->FindObject("MSingles");
183 if (!fSingles)
184 {
185 *fLog << /* err << */ "MSingles not found... abort." << endl;
186 return kFALSE;
187 }
188
189 fEvent = (MPedestalSubtractedEvt*)plist->FindObject("MPedestalSubtractedEvt");
190 if (!fEvent)
191 {
192 *fLog << /* err << */ "MPedestalSubtractedEvt not found... abort." << endl;
193 return kFALSE;
194 }
195
196 fNumEntries = 0;
197
198 return kTRUE;
199 }
200 Bool_t ReInit(MParList *plist)
201 {
202 const MRawRunHeader *header = (MRawRunHeader*)plist->FindObject("MRawRunHeader");
203 if (!header)
204 {
205 *fLog << /* err << */ "MRawRunHeader not found... abort." << endl;
206 return kFALSE;
207 }
208
209 const Int_t w = fSingles->GetIntegrationWindow();
210
211 MBinning binsx, binsy;
212 binsx.SetEdges(fSingles->GetNumPixels(), -0.5, fSingles->GetNumPixels()-0.5);
213 binsy.SetEdges(2*150, -2*7.5*w, 2*(60-7.5)*w);
214
215 MH::SetBinning(fSignal, binsx, binsy);
216
217 const UShort_t roi = header->GetNumSamples();
218
219 // Correct for DRS timing!!
220 MBinning binst(roi, -0.5, roi-0.5);
221 MH::SetBinning(fTime, binsx, binst);
222
223 MBinning binspy(2*roi, -roi-0.5, roi-0.5);
224 MH::SetBinning(fPulse, binsx, binspy);
225
226 return kTRUE;
227 }
228
229 Int_t Fill(const MParContainer *par, const Stat_t weight=1)
230 {
231 const Float_t *ptr = fEvent->GetSamples(0);
232 const Short_t roi = fEvent->GetNumSamples();
233
234 for (unsigned int i=0; i<fSingles->GetNumPixels(); i++)
235 {
236 const vector<Single> &result = fSingles->GetVector(i);
237
238 for (unsigned int j=0; j<result.size(); j++)
239 {
240 fSignal.Fill(i, result[j].fSignal);
241 fTime.Fill (i, result[j].fTime);
242
243 if (!ptr)
244 continue;
245
246 const Short_t tm = floor(result[j].fTime);
247
248 for (int s=0; s<roi; s++)
249 fPulse.Fill(i, s-tm, ptr[s]);
250 }
251
252 ptr += roi;
253 }
254
255 fNumEntries++;
256
257 return kTRUE;
258 }
259
260 TH2 *GetSignal() { return &fSignal; }
261 TH2 *GetTime() { return &fTime; }
262 TH2 *GetPulse() { return &fPulse; }
263
264 UInt_t GetNumEntries() const { return fNumEntries; }
265
266 void Draw(Option_t *opt)
267 {
268 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
269
270 AppendPad("");
271
272 pad->Divide(2,2);
273
274 pad->cd(1);
275 gPad->SetBorderMode(0);
276 gPad->SetFrameBorderMode(0);
277 fSignal.Draw("colz");
278
279 pad->cd(2);
280 gPad->SetBorderMode(0);
281 gPad->SetFrameBorderMode(0);
282 fTime.Draw("colz");
283
284 pad->cd(3);
285 gPad->SetBorderMode(0);
286 gPad->SetFrameBorderMode(0);
287 fPulse.Draw("colz");
288 }
289
290 void DrawCopy(Option_t *opt)
291 {
292 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
293
294 AppendPad("");
295
296 pad->Divide(2,2);
297
298 pad->cd(1);
299 gPad->SetBorderMode(0);
300 gPad->SetFrameBorderMode(0);
301 fSignal.DrawCopy("colz");
302
303 pad->cd(2);
304 gPad->SetBorderMode(0);
305 gPad->SetFrameBorderMode(0);
306 fTime.DrawCopy("colz");
307
308 pad->cd(3);
309 gPad->SetBorderMode(0);
310 gPad->SetFrameBorderMode(0);
311 fPulse.DrawCopy("colz");
312 }
313
314 ClassDef(MHSingles, 1)
315};
316
317MStatusDisplay *d = new MStatusDisplay;
318
319MSingles *fSingles;
320
321Int_t PreProcess(MParList *plist)
322{
323// fSingles = (MSingles*)plist->FindCreateObj("MSingles");
324// if (!fSingles)
325// return kFALSE;
326
327// fEvent = (MPedestalSubtractedEvt*)plist->FindObject("MPedestalSubtractedEvt");
328// if (!fEvent)
329// {
330// *fLog << /* err << */ "MPedestalSubtractedEvt not found... abort." << endl;
331// return kFALSE;
332// }
333
334// d->AddTab("Fluct");
335// fluct2->Draw();
336// fluct1->Draw("same");
337
338// fSingles->SetIntegrationWindow(20);
339
340// //if (weights.GetN()==0)
341// // return kFALSE;
342
343 return kTRUE;
344}
345
346Int_t Process()
347{
348 return kTRUE;
349}
350
351Int_t PostProcess()
352{
353 return kTRUE;
354}
355
356Double_t fcn(Double_t *x, Double_t *par);
357
358int recalcSinglepe(const char *file, const char *outfile="", TString option="")
359{
360 bool save = 0;
361
362 // ======================================================
363 //save options
364
365 if ( option.Contains("S") ){
366 save = 1;
367 cout << "...will save" << endl;
368 }
369
370 //Maximum pe order to which the sepctrum will be fitted
371 Int_t maxOrder = 5;
372
373 // ======================================================
374
375 cout << file << endl;
376 TFile fi(file, "READ");
377 if (!fi.IsOpen())
378 {
379 cout << "ERROR - Could not find file " << file << endl;
380 return 2;
381 }
382
383 MStatusArray arr;
384 if (arr.Read()<=0)
385 {
386 cout << "ERROR - could not read MStatusArray." << endl;
387 return 2;
388 }
389
390 // ======================================================
391
392 // Load histogramm for fluctuations of all Pixel
393 TH1 *fluct = (TH1*)arr.FindObjectInCanvas("", "TH1F", "Fluct");
394 if (!fluct)
395 {
396 cout << "WARNING - Reading of fluct failed." << endl;
397 return 2;
398 }
399
400 d->AddTab("Fluct");
401 fluct->DrawCopy();
402
403
404 MHSingles* singles = (MHSingles*) arr.FindObjectInCanvas("MHSingles", "MHSingles", "MHSingles");
405 if (!singles)
406 {
407 cout << "WARNING - Reading of singles failed." << endl;
408 return 2;
409 }
410
411 d->AddTab("MHSingles");
412 singles->DrawCopy("colz");
413
414 TH1 *time = (TH1*)arr.FindObjectInCanvas("Time_py", "TH1D", "Time");
415 if (!fluct)
416 {
417 cout << "WARNING - Reading of time failed." << endl;
418 return 2;
419 }
420
421 d->AddTab("Time");
422 time->DrawCopy();
423
424 TH1 *pulse = (TH1*)arr.FindObjectInCanvas("Pulse_py", "TH1D", "Pulse");
425 if (!fluct)
426 {
427 cout << "WARNING - Reading of pulse failed." << endl;
428 return 2;
429 }
430
431 d->AddTab("Pulse");
432 pulse->DrawCopy();
433
434 // ======================================================
435
436 // Load histogramm for Sum Spectrum of all Pixel
437 TH1 *hSum = (TH1*)arr.FindObjectInCanvas("Sum1", "TH1F", "SumHist");
438 if (!hSum)
439 {
440 cout << "WARNING - Reading of hsum failed." << endl;
441 return 2;
442 }
443
444 // Load fit function for Sum Spectrum of all Pixel
445 TF1 *func1 = (TF1*)arr.FindObjectInCanvas("spektrum", "TF1", "SumHist");
446 if (!func1)
447 {
448 cout << "WARNING - Reading of func failed." << endl;
449 return 2;
450 }
451
452 //Draw histogramm for Sum Spectrum of all Pixel
453// d->AddTab("SumHist");
454// gPad->SetLogy();
455// gPad->SetGrid();
456// hSum->DrawCopy("hist");
457// func1->DrawCopy("SAME");
458
459 // Load TH2 histogramm for Pixel Spectra
460 TH2F *hsignal = (TH2F*)arr.FindObjectInCanvas("Signal", "TH2F", "MHSingles");
461 if (!hsignal)
462 {
463 cout << "WARNING - Reading of hSpectra failed." << endl;
464 return 2;
465 }
466
467 // Draw TH2 histogramm for Pixel Spectra
468 d->AddTab("Signal");
469 hsignal->SetYTitle("Counts [a.u.]");
470 hsignal->DrawCopy("colz");
471
472
473 // ======================================================
474
475 // true: Display correctly mapped pixels in the camera displays
476 // but the value-vs-index plot is in software/spiral indices
477 // false: Display pixels in hardware/linear indices,
478 // but the order is the camera display is distorted
479 bool usemap = true;
480
481 // map file to use (get that from La Palma!)
482 const char *map = usemap ? "/scratch_nfs/FACTmap111030.txt" : NULL;
483
484
485
486 // ======================================================
487
488 if (map && gSystem->AccessPathName(map, kFileExists))
489 {
490 gLog << "ERROR - Cannot access mapping file '" << map << "'" << endl;
491 return 1;
492 }
493
494 // ======================================================
495
496 //MStatusDisplay *d = new MStatusDisplay;
497
498 MBadPixelsCam badpixels;
499 badpixels.InitSize(1440);
500 badpixels[ 424].SetUnsuitable(MBadPixelsPix::kUnsuitable);
501 badpixels[ 583].SetUnsuitable(MBadPixelsPix::kUnsuitable);
502 badpixels[ 830].SetUnsuitable(MBadPixelsPix::kUnsuitable);
503 badpixels[ 923].SetUnsuitable(MBadPixelsPix::kUnsuitable);
504 badpixels[1208].SetUnsuitable(MBadPixelsPix::kUnsuitable);
505 badpixels[1399].SetUnsuitable(MBadPixelsPix::kUnsuitable);
506
507 // Twin pixel
508 // 113
509 // 115
510 // 354
511 // 423
512 // 1195
513 // 1393
514
515 // ======================================================
516
517 gLog << endl;
518 gLog.Separator("recalculation of gain from pixel spectra");
519
520 MTaskList tlist1;
521
522 MParList plist1;
523 plist1.AddToList(&tlist1);
524 plist1.AddToList(&badpixels);
525
526 MEvtLoop loop1(file);
527 loop1.SetDisplay(d);
528 loop1.SetParList(&plist1);
529
530 // ------------------ Setup the tasks ---------------
531
532 MRawFitsRead read1;
533 read1.LoadMap(map);
534
535 MGeomApply apply1;
536
537 MTaskInteractive mytask;
538 mytask.SetPreProcess(PreProcess);
539 mytask.SetProcess(Process);
540 mytask.SetPostProcess(PostProcess);
541
542 // ------------------ Setup eventloop and run analysis ---------------
543
544 tlist1.AddToList(&mytask);
545
546// if (!loop1.Eventloop(max1))
547// return 10;
548
549 if (!loop1.GetDisplay())
550 return 11;
551
552 gStyle->SetOptFit(kTRUE);
553
554 // ======================================================
555 // ----------------------- Spectrum Fit -----------------
556 // AFTER this the Code for the spektrum fit follows
557 // access the spektrumhistogram by the pointer *hist
558
559 MGeomCamFACT fact;
560 MHCamera hcam(fact);
561 MHCamera hcam2(fact);
562 MHCamera hNormCam(fact);
563
564 // ======================================================
565 // create histograms to plot spectra and parameter distributions
566
567 TH1F hAmplitude ("Rate", "Average number of dark counts per event",
568 200, 0, 20);
569 TH1F hGain ("Gain", "Gain distribution",
570 50, 150, 350);
571 TH1F hGain2 ("NormGain", "Normalized gain distribution",
572 51, 0.5, 1.5);
573 TH1F hGainRMS ("RelSigma", "Distribution of rel. Sigma",
574 100, 0, 30);
575 TH1F hCrosstlk ("Crosstalk", "Distribution of crosstalk probability",
576 100, 0, 25);
577 TH1F hOffset ("Offset", "Baseline offset distribution",
578 50, -7.5, 2.5);
579 TH1F hFitProb ("FitProb", "FitProb distribution",
580 50, 0, 100);
581 TH1F hNoise ("Noise", "Noise distribution",
582 80, -10, 30);
583 TH1F hXtalkPot ("XtalkPot", "Crosstalk Potential Distribution",
584 200, 0, 2);
585
586 hAmplitude.SetXTitle("Amplitude [a.u.]");
587 hGain.SetXTitle("gain [a.u.]");
588 hGain2.SetXTitle("gain [a.u.]");
589 hGainRMS.SetXTitle("gainRMS [a.u.]");
590 //hGainRMS.SetStats(false);
591 hCrosstlk.SetXTitle("crosstalk [a.u.]");
592 hOffset.SetXTitle("baseline offset [a.u.]");
593 hFitProb.SetXTitle("FitProb p-value [a.u.]");
594 //hFitProb.SetStats(false);
595
596 TH1F hSumScale("Sum2", "Signal spectrum of all pixels", 100, 0.05, 10.05);
597 hSumScale.SetXTitle("pulse integral [pe]");
598 hSumScale.SetYTitle("Rate");
599 hSumScale.SetStats(false);
600 hSumScale.Sumw2();
601
602 // define fit function for Amplitudespectrum
603 TF1 func("spektrum", fcn, 0, hSum->GetXaxis()->GetXmax(), 7);
604 func.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset", "Noise", "XTalkPot");
605 func.SetLineColor(kBlue);
606
607 // ======================================================
608
609 // Map for which pixel shall be plotted and which not
610 TArrayC usePixel(1440);
611 memset(usePixel.GetArray(), 1, 1440);
612
613 // List of Pixel that should be ignored in camera view
614 usePixel[424] = 0;
615 usePixel[923] = 0;
616 usePixel[1208] = 0;
617 usePixel[583] = 0;
618 usePixel[830] = 0;
619 usePixel[1399] = 0;
620 usePixel[113] = 0;
621 usePixel[115] = 0;
622 usePixel[354] = 0;
623 usePixel[423] = 0;
624 usePixel[1195] = 0;
625 usePixel[1393] = 0;
626
627 // ======================================================
628 //-------------------_fit sum spectrum-------------------
629
630 gLog << endl;
631 gLog.Separator("Fitting Sum Spectrum");
632
633 //Set Status Display Output and Prograssbar position
634 d->SetStatusLine( "...fitting Sum Spectrum", 0);
635 d->SetProgressBarPosition( 0.0 , 1);
636
637 const Int_t maxbin = hSum->GetMaximumBin();
638 const Double_t maxpos = hSum->GetBinCenter(maxbin);
639 const Double_t binwidth = hSum->GetBinWidth(maxbin);
640 const Double_t ampl = hSum->GetBinContent(maxbin);
641
642 // ======================================================
643
644 // full width half maximum of sumspectra
645 double fwhmSum = 0;
646
647 //Calculate full width half Maximum
648 for (int i=1; i<maxbin; i++)
649 if (hSum->GetBinContent(maxbin-i)+hSum->GetBinContent(maxbin+i)<ampl*2)
650 {
651 fwhmSum = (2*i+1)*binwidth;
652 break;
653 }
654
655 if (fwhmSum==0)
656 {
657 cout << "Could not determine start value for sigma." << endl;
658 }
659
660 // ======================================================
661
662 //Set fit parameters
663 Double_t sum_par[7] =
664 {
665 ampl, maxpos, fwhmSum*1.4, 0.1, -1, 0, 1
666 };
667
668 Double_t sum_par_err[7];
669
670 Double_t fitmin = maxpos-fwhmSum*3;
671 Double_t fitmax = hSum->GetXaxis()->GetXmax();
672
673 //set status display prograssbar position
674 d->SetProgressBarPosition(0.05 , 1);
675
676 //Fit and draw spectrum
677 func.SetParameters(sum_par);
678 func.SetRange(fitmin, fitmax);
679 hSum->Fit(&func, "INOQSR", "", fitmin, fitmax);
680 func.GetParameters(sum_par);
681
682 //Save Parameter Errors for output
683 for (int i = 0; i < 7; i++){
684 sum_par_err[i] = func.GetParError(i);
685 }
686
687 //Set Parameters for following fits
688 const float gain = sum_par[1];
689 const float GainRMS = sum_par[2];
690 const float Crosstlk = sum_par[3];
691 const float Offset = sum_par[4];
692 const float Noise = sum_par[5];
693 const float XtalkPot = sum_par[6];
694
695 cout << "--------------------" << endl;
696 cout << "MaxPos: " << maxpos << " \t +/- " << sum_par_err[0] << endl;
697 cout << "FWHM: " << fwhmSum << endl;
698 cout << "GAIN: " << gain << " \t +/- " << sum_par_err[1] << endl;
699 cout << "RMS: " << GainRMS << " \t +/- " << sum_par_err[2] << endl;
700 cout << "xTalk: " << Crosstlk << " \t +/- " << sum_par_err[3] << endl;
701 cout << "XtPot: " << XtalkPot << " \t +/- " << sum_par_err[6] << endl;
702 cout << "Noise: " << Noise << " \t +/- " << sum_par_err[5] << endl;
703 cout << "Offset: " << Offset << " \t +/- " << sum_par_err[4] << endl;
704 cout << "--------------------" << endl;
705
706 gROOT->SetSelectedPad(0);
707
708 //set status display prograssbar position
709 d->SetProgressBarPosition(0.1 , 1);
710
711 // ======================================================
712
713 TCanvas &c11 = d->AddTab("SumHist");
714 c11.cd();
715 gPad->SetLogy();
716 gPad->SetGridx();
717 gPad->SetGridy();
718 hSum->DrawCopy("hist");
719 func.DrawCopy("SAME");
720
721 // ======================================================
722 //----------Gain Calculation for each Pixel--------------
723
724 gLog << endl;
725 gLog.Separator("pixelwise gain calculation");
726
727 // counter for number of processed pixel
728 int count_ok = 0;
729
730 // marker for pixel with abnormal behavior
731 bool suspicous = 0;
732
733 // number of pixels to process
734 UInt_t nPixels = hsignal->GetNbinsX();
735
736 //set status display status
737 d->SetStatusLine( "...fitting pixels", 0);
738
739 // ======================================================
740
741 // Loop over Pixels
742 for (UInt_t pixel = 0; pixel < nPixels; pixel++)
743 {
744 //set status display prograssbar position
745 d->SetProgressBarPosition(0.1 + 0.8*( (double)(pixel+1)/(double)nPixels ), 1);
746 suspicous = 0;
747
748 // jump to next pixel if the current is marked as broken
749 if (usePixel[pixel]==0)
750 continue;
751
752 //Projectipon of a certain Pixel out of Ampl.Spectrum
753 TH1D *hist = hsignal->ProjectionY("proj", pixel+1, pixel+1);
754 hist->SetDirectory(0);
755
756 //Rebin Projection
757 hist->Rebin(2);
758
759 // Callculate values of the bin with maximum entries
760 const Int_t maxbin = hist->GetMaximumBin();
761 const Double_t maxPos = hist->GetBinCenter(maxbin);
762
763 // Calculate range for the fit
764 const double fit_min = maxPos-GainRMS*2;
765 const double fit_max = maxPos+gain*(maxOrder-0.5);
766
767 // ======================================================
768
769 // Set fit parameter start values
770 sum_par[0] = hist->GetBinContent(maxbin);
771// sum_par[2] = first_gaus_RMS;
772// sum_par[3] = second_gaus_ampl/first_gaus_ampl;
773// sum_par[4] = -1*(second_gaus_ampl-2*first_gaus_ampl);
774 func.SetParameters(sum_par);
775// func.FixParameter(0, first_gaus_ampl);
776// func.SetParLimits(0, 0.8*sum_par[0], 1.1*sum_par[0]);
777// func.SetParLimits(1, maxPos - 2*GainRMS, maxPos + 2*GainRMS);
778// func.SetParLimits(2, -first_gaus_RMS, +first_gaus_RMS);
779// func.FixParameter(5, sum_par[5]);
780// func.SetParLimits(5, 0.7, 1.1);
781// func.FixParameter(6, sum_par[6]);
782// func.SetParLimits(6, 0.1, 1.9);
783
784 // ======================================================
785 // -------------- Fit Pixels spectrum -------------------
786
787 // Save fit result to a variable
788 const TFitResultPtr rc = hist->Fit(&func, "LLN0QS", "", fit_min, fit_max);
789
790 // marker to mark if fit was successfull or not
791 const bool ok = int(rc)==0;
792
793 // Save Parametervalues for statistics and bug fixing
794 const double fit_prob = rc->Prob();
795// const float fMaxAmpl = func.GetParameter(0);
796 const float fGain = func.GetParameter(1);
797 const float fAmplitude = func.Integral(Offset, maxOrder*fGain+Offset); //func.GetParameter(0);
798 const float f2GainRMS = func.GetParameter(2);
799 const float fCrosstlk = func.GetParameter(3);
800 const float fOffset = func.GetParameter(4)/30;
801 const float fNoise = func.GetParameter(5);
802 const float fXtlkPot = func.GetParameter(6);
803
804 // Fill Parameter value hgistograms
805 hAmplitude.Fill(fAmplitude);
806 hGain.Fill(fGain);
807 hGain2.Fill(fGain/gain);
808 hCrosstlk.Fill(fCrosstlk*100);
809 hOffset.Fill(fOffset);
810 hFitProb.Fill(100*fit_prob);
811 hGainRMS.Fill(100*f2GainRMS/fGain);
812 hNoise.Fill(fNoise);
813 hXtalkPot.Fill(fXtlkPot);
814
815 // Plot calculated Gain values to Camera Display
816 hcam.SetBinContent(pixel+1, fGain);
817 hcam2.SetBinContent(pixel+1, 100*f2GainRMS/fGain);
818 hNormCam.SetBinContent(pixel+1, fGain/gain);
819
820 // ======================================================
821
822 // cancel out pixel where the fit was not succsessfull
823 usePixel[pixel] = ok;
824
825 // mark pixels suspicious with negative noise
826 if ( fNoise < 0)
827 suspicous = 1;
828
829 // mark pixels suspicious with negative GainRMS
830 if ( f2GainRMS < 0)
831 suspicous = 1;
832
833 // mark pixels suspicious with with large deviation from mean gain
834 if ( fGain < gain - 1.4*GainRMS || fGain > gain + 1.4*GainRMS)
835 suspicous = 1;
836
837 //plott especialy marked pixel
838 if (count_ok<3 || suspicous == 1 || !ok)
839 {
840 hist->SetName("Pix%d");
841 TCanvas &c = d->AddTab(Form("Pix%d", pixel));
842 c.cd();
843 gPad->SetLogy();
844 gPad->SetGridx();
845 gPad->SetGridy();
846
847 rc->Print("V");
848
849// hist->SetStats(false);
850 hist->SetYTitle("Counts [a.u.]");
851 hist->SetName(Form("Pix%d", pixel));
852 hist->DrawCopy()->SetDirectory(0);
853 func.DrawCopy("SAME")->SetRange(fit_min, fit_max);
854 }
855
856 if (!ok)
857 {
858 cout << "Pixel " << setw(4) << pixel << ": failed!" << endl;
859 delete hist;
860 continue;
861 }
862
863 // ======================================================
864
865 //Fill sum spectrum
866 for (int b=1; b<=hist->GetNbinsX(); b++)
867 hSumScale.Fill(
868 (hist->GetBinCenter(b)-fOffset)/fGain,
869 hist->GetBinContent(b)
870 );
871
872 count_ok++;
873
874 delete hist;
875 }
876
877 // define a general gaus function
878 TF1 fgaus("fgaus", "gaus");
879 hGain.Fit(&fgaus, "N0QS");
880
881 // ======================================================
882 //------------------Draw histograms----------------------
883
884 gLog << endl;
885 gLog.Separator("Drawing Histograms");
886
887 d->SetStatusLine( "...drawing Histograms", 0);
888 d->SetProgressBarPosition(0.9 , 1);
889
890 // cancel out pixel in camera view
891 hcam.SetUsed(usePixel);
892 hcam2.SetUsed(usePixel);
893
894 // ======================================================
895
896 TCanvas &canv = d->AddTab("Gain");
897 canv.cd();
898 canv.Divide(2,2);
899
900 canv.cd(1);
901 hcam.SetNameTitle("gain", "Gain distribution over whole camera");
902 hcam.DrawCopy();
903
904 canv.cd(2);
905 hcam2.SetNameTitle("gainRMS", "GainRMS distribution over whole camera");
906 hcam2.DrawCopy();
907
908 TCanvas &cam_canv = d->AddTab("Camera_Gain");
909
910 // ======================================================
911 //---------- Draw gain corrected sum specetrum ----------
912
913 gROOT->SetSelectedPad(0);
914 TCanvas &c12 = d->AddTab("GainHist");
915 c12.cd();
916 gPad->SetLogy();
917 gPad->SetGridx();
918 gPad->SetGridy();
919
920 hSumScale.DrawCopy("hist");
921
922 //---------fit gain corrected sum spectrum ---------------
923
924 // Callculate values of the bin with maximum entries
925 const Int_t maxbin2 = hSumScale.GetMaximumBin();
926 const Double_t ampl2 = hSumScale.GetBinContent(maxbin2);
927
928 //Set fit start parameters
929 Double_t par2[7] =
930 {
931 ampl2, 1, GainRMS/fgaus.GetParameter(1), Crosstlk, 0, 0, 1
932 };
933 func.SetParameters(par2);
934
935 // fix gain to 1 as it was normalized allready
936 func.FixParameter(1, 1);
937
938 // fix offset parameter
939 // func.FixParameter(4, 0);
940
941 // set progressbar possition
942 d->SetStatusLine( "...fitting corr. spectrum", 0);
943
944 // Fit corrected Spectrum with Likelyhood Method
945 const TFitResultPtr rc = hSumScale.Fit(&func, "LLEN0QS", "", 0.75, 10);
946
947 func.DrawCopy("SAME")->SetRange(0.5, 10);
948
949 // Print fit results
950 rc->Print("V");
951
952 // set progressbar possition
953 d->SetProgressBarPosition(0.95 , 1);
954
955 // ======================================================
956
957 TCanvas &c2 = d->AddTab("Result");
958 c2.Divide(3,2);
959 c2.cd(1);
960 gPad->SetLogy();
961 hXtalkPot.DrawCopy();
962
963 c2.cd(2);
964 gPad->SetLogy();
965 hGain.DrawCopy();
966
967 //plot normailzed camera view
968 cam_canv.cd();
969 hNormCam.SetStats(false);
970
971 hNormCam.SetNameTitle("GainCam", "Normalized gain distribution");
972 hNormCam.SetUsed(usePixel);
973 hNormCam.DrawCopy();
974
975 c2.cd(3);
976 gPad->SetLogy();
977 hGainRMS.DrawCopy();
978
979 c2.cd(4);
980 gPad->SetLogy();
981 hCrosstlk.DrawCopy();
982
983 c2.cd(5);
984 gPad->SetLogy();
985 hOffset.DrawCopy();
986
987 c2.cd(6);
988 gPad->SetLogy();
989 hNoise.DrawCopy();
990
991 // ======================================================
992
993 TCanvas &c3 =d->AddTab("Prop");
994 c3.Divide(2,1);
995 c3.cd(1);
996 gPad->SetGrid();
997 gPad->SetLogy();
998 hFitProb.DrawCopy();
999
1000 c3.cd(2);
1001 gPad->SetGrid();
1002 gPad->SetLogy();
1003 hAmplitude.DrawCopy();
1004
1005 // ======================================================
1006
1007 d->AddTab("Norm");
1008 gPad->SetGrid();
1009 gPad->SetLogy();
1010 TH1 *hh = hGain2.DrawCopy();
1011 hh->Fit("gaus");
1012
1013 // ======================================================
1014
1015 //save results to file
1016 if (save){
1017 d->SetStatusLine( "...saving results to rootfile", 0);
1018 d->SaveAs(outfile);
1019 cout << "..success!" << endl;
1020 }
1021 d->SetProgressBarPosition(1 , 1);
1022 d->SetStatusLine( "...done", 0);
1023 return 0;
1024}
1025
1026Double_t fcn(Double_t *xx, Double_t *par)
1027{
1028 Double_t x = xx[0];
1029
1030 Double_t ampl = par[0];
1031 Double_t gain = par[1];
1032 Double_t sigma = par[2];
1033 Double_t cross = par[3];
1034 Double_t shift = par[4];
1035 Double_t noise = par[5];
1036 Double_t crossPot = par[6];
1037
1038 Double_t xTalk = 1;
1039
1040 Double_t y = 0;
1041 for (int order = 1; order < 14; order++)
1042 {
1043 Double_t arg = (x - order*gain - shift)/(order==1?sigma+noise:sigma);
1044
1045 Double_t gauss = TMath::Exp(-0.5 * arg * arg/order);
1046
1047 y += xTalk*gauss;
1048 xTalk = Power(cross, Power(order, crossPot) );
1049 }
1050
1051 return ampl*y;
1052}
1053
1054// end of PlotMedianEachSliceOfPulse
1055//----------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.