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

Last change on this file since 16858 was 14246, checked in by Jens Buss, 12 years ago
bug fix in scaled spectrum fit, took of likelihood fit
  • Property svn:executable set to *
File size: 32.7 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#include <Getline.h>
11
12#include <cstdio>
13#include <stdio.h>
14#include <stdint.h>
15
16#include "MH.h"
17#include "MLog.h"
18#include "MLogManip.h"
19#include "MDirIter.h"
20#include "MFillH.h"
21#include "MEvtLoop.h"
22#include "MCamEvent.h"
23#include "MGeomApply.h"
24#include "MTaskList.h"
25#include "MParList.h"
26#include "MContinue.h"
27#include "MBinning.h"
28#include "MHDrsCalib.h"
29#include "MDrsCalibApply.h"
30#include "MRawFitsRead.h"
31#include "MBadPixelsCam.h"
32#include "MStatusDisplay.h"
33#include "MTaskInteractive.h"
34#include "MPedestalSubtractedEvt.h"
35#include "MHCamera.h"
36#include "MGeomCamFACT.h"
37#include "MRawRunHeader.h"
38
39using namespace std;
40
41MPedestalSubtractedEvt *fEvent = 0;
42MLog *fLog = &gLog;
43
44struct Single
45{
46 float fSignal;
47 float fTime;
48};
49
50class MSingles : public MParContainer, public MCamEvent
51{
52 Int_t fIntegrationWindow;
53
54 vector<vector<Single>> fData;
55
56 public:
57 MSingles(const char *name=NULL, const char *title=NULL) : fIntegrationWindow(0)
58 {
59 fName = name ? name : "MSingles";
60 }
61
62 void InitSize(const UInt_t i)
63 {
64 fData.resize(i);
65 }
66
67 vector<Single> &operator[](UInt_t i)
68 {
69 return fData[i];
70 }
71 vector<Single> &GetVector(UInt_t i)
72 {
73 return fData[i];
74 }
75
76 UInt_t GetNumPixels() const
77 {
78 return fData.size();
79 }
80
81 void SetIntegrationWindow(Int_t w)
82 {
83 fIntegrationWindow = w;
84 }
85 Int_t GetIntegrationWindow() const
86 {
87 return fIntegrationWindow;
88 }
89
90 Bool_t GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const
91 {
92 return kTRUE;
93 }
94 void DrawPixelContent(Int_t num) const { }
95
96 ClassDef(MSingles, 1)
97};
98
99class MH2F : public TH2F
100{
101 public:
102 Int_t Fill(Double_t x,Double_t y)
103 {
104 Int_t binx, biny, bin;
105 fEntries++;
106 binx = TMath::Nint(x+1);
107 biny = fYaxis.FindBin(y);
108 bin = biny*(fXaxis.GetNbins()+2) + binx;
109 AddBinContent(bin);
110
111 if (fSumw2.fN) ++fSumw2.fArray[bin];
112
113 if (biny == 0 || biny > fYaxis.GetNbins())
114 {
115 if (!fgStatOverflows) return -1;
116 }
117
118 ++fTsumw;
119 ++fTsumw2;
120 fTsumwy += y;
121 fTsumwy2 += y*y;
122 return bin;
123 }
124 ClassDef(MH2F, 1);
125};
126
127class MProfile2D : public TProfile2D
128{
129 public:
130 Int_t Fill(Double_t x, Double_t y, Double_t z)
131 {
132 Int_t bin,binx,biny;
133 fEntries++;
134 binx =TMath::Nint(x+1);
135 biny =fYaxis.FindBin(y);
136 bin = GetBin(binx, biny);
137 fArray[bin] += z;
138 fSumw2.fArray[bin] += z*z;
139 fBinEntries.fArray[bin] += 1;
140
141 if (fBinSumw2.fN) fBinSumw2.fArray[bin] += 1;
142
143 if (biny == 0 || biny > fYaxis.GetNbins())
144 {
145 if (!fgStatOverflows) return -1;
146 }
147
148 ++fTsumw;
149 ++fTsumw2;
150 fTsumwy += y;
151 fTsumwy2 += y*y;
152 fTsumwz += z;
153 fTsumwz2 += z*z;
154 return bin;
155 }
156 ClassDef(MProfile2D, 1);
157};
158
159
160
161class MHSingles : public MH
162{
163 TH2F fSignal; //changed from MH2F to TH2F
164 TH2F fTime; //changed from MH2F to TH2F
165 TProfile2D fPulse; //changed from MProfile2D to TProfile2D
166
167
168 UInt_t fNumEntries;
169
170 MSingles *fSingles; //!
171 MPedestalSubtractedEvt *fEvent; //!
172
173 public:
174 MHSingles() : fNumEntries(0), fSingles(0), fEvent(0)
175 {
176 fName = "MHSingles";
177
178 fSignal.SetName("Signal");
179 fSignal.SetTitle("pulse integral spectrum");
180 fSignal.SetXTitle("pixel [SoftID]");
181 fSignal.SetYTitle("time [a.u.]");
182 fSignal.SetDirectory(NULL);
183
184 fTime.SetName("Time");
185 fTime.SetTitle("arival time spectrum");
186 fTime.SetXTitle("pixel [SoftID]");
187 fTime.SetYTitle("time [a.u.]");
188 fTime.SetDirectory(NULL);
189
190 fPulse.SetName("Pulse");
191 fPulse.SetTitle("average pulse shape spectrum");
192 fPulse.SetXTitle("pixel [SoftID]");
193 fPulse.SetYTitle("time [a.u.]");
194 fPulse.SetDirectory(NULL);
195 }
196
197 Bool_t SetupFill(const MParList *plist)
198 {
199 fSingles = (MSingles*)plist->FindObject("MSingles");
200
201 if (!fSingles)
202 {
203 *fLog << /* err << */ "MSingles not found... abort." << endl;
204 return kFALSE;
205 }
206
207 fEvent = (MPedestalSubtractedEvt*)plist->FindObject("MPedestalSubtractedEvt");
208
209 if (!fEvent)
210 {
211 *fLog << /* err << */ "MPedestalSubtractedEvt not found... abort." << endl;
212 return kFALSE;
213 }
214
215 fNumEntries = 0;
216
217 return kTRUE;
218 }
219 Bool_t ReInit(MParList *plist)
220 {
221 const MRawRunHeader *header = (MRawRunHeader*)plist->FindObject("MRawRunHeader");
222
223 if (!header)
224 {
225 *fLog << /* err << */ "MRawRunHeader not found... abort." << endl;
226 return kFALSE;
227 }
228
229 const Int_t w = fSingles->GetIntegrationWindow();
230
231 MBinning binsx, binsy;
232
233 binsx.SetEdges(fSingles->GetNumPixels(), -0.5, fSingles->GetNumPixels()-0.5);
234
235 binsy.SetEdges(2*150, -2*7.5*w, 2*(60-7.5)*w);
236
237 MH::SetBinning(fSignal, binsx, binsy);
238
239 const UShort_t roi = header->GetNumSamples();
240
241 // Correct for DRS timing!!
242 MBinning binst(roi, -0.5, roi-0.5);
243
244 MH::SetBinning(fTime, binsx, binst);
245
246 MBinning binspy(2*roi, -roi-0.5, roi-0.5);
247
248 MH::SetBinning(fPulse, binsx, binspy);
249
250 return kTRUE;
251 }
252
253 Int_t Fill(const MParContainer *par, const Stat_t weight=1)
254 {
255 const Float_t *ptr = fEvent->GetSamples(0);
256 const Short_t roi = fEvent->GetNumSamples();
257
258 for (unsigned int i=0; i<fSingles->GetNumPixels(); i++)
259 {
260 const vector<Single> &result = fSingles->GetVector(i);
261
262 for (unsigned int j=0; j<result.size(); j++)
263 {
264 fSignal.Fill(i, result[j].fSignal);
265 fTime.Fill (i, result[j].fTime);
266
267 if (!ptr)
268 continue;
269
270 const Short_t tm = floor(result[j].fTime);
271
272 for (int s=0; s<roi; s++)
273 fPulse.Fill(i, s-tm, ptr[s]);
274 }
275
276 ptr += roi;
277 }
278
279 fNumEntries++;
280
281 return kTRUE;
282 }
283
284 TH2 *GetSignal()
285 {
286 return &fSignal;
287 }
288 TH2 *GetTime()
289 {
290 return &fTime;
291 }
292 TH2 *GetPulse()
293 {
294 return &fPulse;
295 }
296
297 UInt_t GetNumEntries() const
298 {
299 return fNumEntries;
300 }
301
302 void Draw(Option_t *opt)
303 {
304 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
305
306 AppendPad("");
307
308 pad->Divide(2,2);
309
310 pad->cd(1);
311 gPad->SetBorderMode(0);
312 gPad->SetFrameBorderMode(0);
313 fSignal.Draw("colz");
314
315 pad->cd(2);
316 gPad->SetBorderMode(0);
317 gPad->SetFrameBorderMode(0);
318 fTime.Draw("colz");
319
320 pad->cd(3);
321 gPad->SetBorderMode(0);
322 gPad->SetFrameBorderMode(0);
323 fPulse.Draw("colz");
324 }
325
326 ClassDef(MHSingles, 1)
327};
328
329MStatusDisplay *d = new MStatusDisplay;
330
331MSingles *fSingles;
332
333TH1F *fluct1 = new TH1F("", "", 200, -200, 200);
334TH1F *fluct2 = new TH1F("", "", 100, -50, 50);
335
336TGraph weights("weights.txt");
337
338Int_t PreProcess(MParList *plist)
339{
340 fSingles = (MSingles*)plist->FindCreateObj("MSingles");
341
342 if (!fSingles)
343 return kFALSE;
344
345 fEvent = (MPedestalSubtractedEvt*)plist->FindObject("MPedestalSubtractedEvt");
346
347 if (!fEvent)
348 {
349 *fLog << /* err << */ "MPedestalSubtractedEvt not found... abort." << endl;
350 return kFALSE;
351 }
352
353 d->AddTab("Fluct");
354 fluct2->Draw();
355 fluct1->Draw("same");
356
357 fSingles->SetIntegrationWindow(20);
358
359 //if (weights.GetN()==0)
360 // return kFALSE;
361
362 return kTRUE;
363}
364
365Int_t Process()
366{
367 const UInt_t roi = fEvent->GetNumSamples();
368
369 const size_t navg = 10;
370 const float threshold = 5;
371 const UInt_t start_skip = 20;
372 const UInt_t end_skip = 10;
373 const UInt_t integration_size = 10*3;
374 const UInt_t max_search_window = 30;
375
376 vector<float> val(roi-navg);//result of first sliding average
377
378 for (int pix=0; pix<fEvent->GetNumPixels(); pix++)
379 {
380 vector<Single> &result = fSingles->GetVector(pix);
381 result.clear();
382
383 const Float_t *ptr = fEvent->GetSamples(pix);
384
385 //Sliding window
386 for (size_t i=0; i<roi-navg; i++)
387 {
388 val[i] = 0;
389
390 for (size_t j=i; j<i+navg; j++)
391 val[i] += ptr[j];
392
393 val[i] /= navg;
394 }
395
396 //peak finding via threshold
397 for (UInt_t i=start_skip; i<roi-navg-end_skip-30; i++)
398 {
399 //search for threshold crossings
400 if (val[i+0]>threshold ||
401 val[i+4]>threshold ||
402
403 val[i+5]<threshold ||
404 val[i+9]<threshold)
405 continue;
406
407 //search for maximum after threshold crossing
408 UInt_t k_max = i+5;
409
410 for (UInt_t k=i+5; k<i+max_search_window; k++)
411 {
412 if (val[k] > val[k_max])
413 k_max = k;
414 }
415
416 if (k_max == i+5 || k_max == i + max_search_window-1) continue;
417
418 //search for half maximum before maximum
419 UInt_t k_half_max = 0;
420
421 for (UInt_t k=k_max; k>k_max-25; k--)
422 {
423 if (val[k-1] < val[k_max]/2 &&
424 val[k] >= val[k_max]/2 )
425 {
426 k_half_max = k;
427 break;
428 }
429 }
430
431 if (k_half_max > roi-navg-end_skip-max_search_window - integration_size)
432 continue;
433
434 if (k_half_max == 0)
435 continue;
436
437 if (k_max - k_half_max > 14)
438 continue;
439
440 fluct2->Fill(k_max - k_half_max);
441
442 // Evaluate arrival time more precisely!!!
443 // Get a better integration window
444
445 Single single;
446 single.fSignal = 0;
447 single.fTime = k_half_max;
448
449
450 // Crossing of the threshold is at 5
451 for (UInt_t j=k_half_max; j<k_half_max+integration_size; j++)
452 single.fSignal += ptr[j];
453
454 result.push_back(single);
455
456 i += 5+integration_size;
457 }
458 }
459
460 return kTRUE;
461}
462
463Int_t PostProcess()
464{
465 return kTRUE;
466}
467
468Double_t fcn(Double_t *x, Double_t *par);
469Double_t fcn_cross(Double_t *x, Double_t *par);
470
471void FitGaussiansToSpectrum( UInt_t nfunc, TH1 &hSource, TH1 &hDest,
472 Int_t maxbin, double fwhm, Double_t gain );
473
474Double_t MedianOfH1 (
475 TH1* inputHisto
476);
477
478int singlepe(
479 const char *runfile, const char *drsfile, const char *outfile
480)
481{
482
483 // ======================================================
484
485// const char *drsfile = "/fact/raw/2012/06/25/20120625_112.drs.fits.gz";
486
487 MDirIter iter;
488 iter.AddDirectory(gSystem->DirName(runfile), gSystem->BaseName(runfile));
489
490 // ======================================================
491
492 // true: Display correctly mapped pixels in the camera displays
493 // but the value-vs-index plot is in software/spiral indices
494 // false: Display pixels in hardware/linear indices,
495 // but the order is the camera display is distorted
496 bool usemap = true;
497
498 // map file to use (get that from La Palma!)
499 const char *map = usemap ? "/scratch_nfs/FACTmap111030.txt" : NULL;
500
501 // ------------------------------------------------------
502
503 Long_t max1 = 0;
504
505 // ======================================================
506
507 if (map && gSystem->AccessPathName(map, kFileExists))
508 {
509 gLog << "ERROR - Cannot access mapping file '" << map << "'" << endl;
510 return 1;
511 }
512
513 if (gSystem->AccessPathName(drsfile, kFileExists))
514 {
515 gLog << "ERROR - Cannot access drsfile file '" << drsfile << "'" << endl;
516 return 2;
517 }
518
519 // --------------------------------------------------------------------------------
520
521 gLog.Separator("Singles");
522 gLog << "Extract singles with '" << drsfile << "'" << endl;
523 gLog << endl;
524
525 // ------------------------------------------------------
526
527 MDrsCalibration drscalib300;
528
529 if (!drscalib300.ReadFits(drsfile))
530 return 4;
531
532 // ------------------------------------------------------
533
534 iter.Sort();
535 iter.Print();
536
537 // ======================================================
538
539 //MStatusDisplay *d = new MStatusDisplay;
540
541 MBadPixelsCam badpixels;
542 badpixels.InitSize(1440);
543 badpixels[ 424].SetUnsuitable(MBadPixelsPix::kUnsuitable);
544 badpixels[ 583].SetUnsuitable(MBadPixelsPix::kUnsuitable);
545 badpixels[ 830].SetUnsuitable(MBadPixelsPix::kUnsuitable);
546 badpixels[ 923].SetUnsuitable(MBadPixelsPix::kUnsuitable);
547 badpixels[1208].SetUnsuitable(MBadPixelsPix::kUnsuitable);
548 badpixels[1399].SetUnsuitable(MBadPixelsPix::kUnsuitable);
549
550 // Twin pixel
551 // 113
552 // 115
553 // 354
554 // 423
555 // 1195
556 // 1393
557
558 //MDrsCalibrationTime timecam;
559
560 MH::SetPalette();
561
562 // ======================================================
563
564 gLog << endl;
565 gLog.Separator("Processing external light pulser run");
566
567 MTaskList tlist1;
568
569 MSingles singles;
570 MRawRunHeader header;
571
572 MParList plist1;
573 plist1.AddToList(&tlist1);
574 plist1.AddToList(&drscalib300);
575 plist1.AddToList(&badpixels);
576 plist1.AddToList(&singles);
577 plist1.AddToList(&header);
578
579 TString Title;
580 Title = iter.Next();
581 iter.Reset();
582 Title += "; ";
583 Title += max1;
584
585 MEvtLoop loop1(Title);
586 loop1.SetDisplay(d);
587 loop1.SetParList(&plist1);
588
589 // ------------------ Setup the tasks ---------------
590
591 MRawFitsRead read1;
592 read1.LoadMap(map);
593 read1.AddFiles(iter);
594
595 MContinue cont1("(MRawEvtHeader.GetTriggerID&0xff00)!=0x100", "SelectCal");
596
597 MGeomApply apply1;
598
599 MDrsCalibApply drsapply1;
600
601 MTaskInteractive mytask;
602 mytask.SetPreProcess(PreProcess);
603 mytask.SetProcess(Process);
604 mytask.SetPostProcess(PostProcess);
605
606 MFillH fill("MHSingles");
607
608 // ------------------ Setup eventloop and run analysis ---------------
609
610 tlist1.AddToList(&read1);
611// tlist1.AddToList(&cont1);
612 tlist1.AddToList(&apply1);
613 tlist1.AddToList(&drsapply1);
614 tlist1.AddToList(&mytask);
615 tlist1.AddToList(&fill);
616
617 if (!loop1.Eventloop(max1))
618 return 10;
619
620 if (!loop1.GetDisplay())
621 return 11;
622
623 gStyle->SetOptFit(kTRUE);
624
625 MHSingles* h = (MHSingles*) plist1.FindObject("MHSingles");
626
627 if (!h)
628 return 0;
629
630 TH2 *hsignal = h->GetSignal();
631 TH2 *htime = h->GetTime();
632 TH2 *hpulse = h->GetPulse();
633
634 d->AddTab("Time");
635 TH1 *ht = htime->ProjectionY();
636 ht->SetYTitle("Counts [a.u.]");
637 ht->Scale(1./1440);
638 ht->Draw();
639
640 d->AddTab("Pulse");
641 TH1 *hp = hpulse->ProjectionY();
642 hp->SetYTitle("Counts [a.u.]");
643 hp->Scale(1./1440);
644 hp->Draw();
645
646
647 // ------------------ Spectrum Fit ---------------
648 // AFTER this the Code for the spektrum fit follows
649 // access the spektrumhistogram by the pointer *hist
650 UInt_t maxOrder = 8;
651
652
653 MGeomCamFACT fact;
654 MHCamera hcam(fact);
655 MHCamera hcam2(fact);
656
657 //built an array of Amplitude spectra
658 TH1F hAmplitude ("Rate", "Average number of dark counts per event",
659 200, 0, 20);
660 hAmplitude.SetXTitle("Amplitude [a.u.]");
661 hAmplitude.SetYTitle("Rate");
662
663 TH1F hGain ("Gain", "Gain distribution",
664 100, 0, 400);
665 hGain.SetXTitle("gain [a.u.]");
666 hGain.SetYTitle("Rate");
667
668 TH1F hGainRMS ("RelSigma", "Distribution of rel. Sigma",
669 100, 0, 30);
670 hGainRMS.SetXTitle("gainRMS [a.u.]");
671 hGainRMS.SetYTitle("Rate");
672 hGainRMS.SetStats(false);
673
674 TH1F hCrosstlk ("Crosstalk", "Distribution of crosstalk probability",
675 100, 0, 25);
676 hCrosstlk.SetXTitle("crosstalk [a.u.]");
677 hCrosstlk.SetYTitle("Rate");
678
679 TH1F hOffset ("Offset", "Baseline offset distribution",
680 50, -7.5, 2.5);
681 hOffset.SetXTitle("baseline offset [a.u.]");
682 hOffset.SetYTitle("Rate");
683
684 TH1F hFitProb ("FitProb", "FitProb distribution",
685 50, 0, 100);
686 hFitProb.SetXTitle("FitProb p-value [a.u.]");
687 hFitProb.SetYTitle("Rate");
688 hFitProb.SetStats(false);
689
690
691 TH1F hSum ("Sum1", "Sum of all pixel's' integral spectra",
692 100, 0, 2000);
693 hSum.SetXTitle("pulse integral [a.u.]");
694 hSum.SetYTitle("Rate");
695 hSum.SetStats(false);
696 hSum.Sumw2();
697
698
699 TH1F hSumScale ("Sum2",
700 "Sum of all pixel's' integral spectra (scaled with gain)",
701 100, 0.05, 10.05);
702 hSumScale.SetXTitle("pulse integral [pe]");
703 hSumScale.SetYTitle("Rate");
704 hSumScale.SetStats(false);
705 hSumScale.Sumw2();
706
707 // define fit function for Amplitudespectrum
708 TF1 func("spektrum", fcn, 0, 1000, 5);
709 func.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset");
710 func.SetLineColor(kBlue);
711
712 // define fit function for Amplitudespectrum
713 TF1 funcSum("sum_spektrum", fcn, 0, 2000, 5);
714 funcSum.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset", "PowOfXtalk");
715 funcSum.SetLineColor(kBlue);
716
717 // define fit function for Amplitudespectrum
718 TF1 funcScaled("gain_sum_spektrum", fcn, 0, 10, 5);
719 funcScaled.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset");
720 funcScaled.SetLineColor(kBlue);
721
722 /*
723 TF1 funcScaledBL("gain_sum_spektrum", fcn, 0, 10, 5);
724 funcScaledBL.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset");
725 funcScaledBL.SetLineColor(kRed);
726 */
727
728
729 // Map for which pixel shall be plotted and which not
730 TArrayC usePixel(1440);
731 memset(usePixel.GetArray(), 1, 1440);
732
733 // List of Pixel that should be ignored in camera view
734 usePixel[424] = 0;
735 usePixel[923] = 0;
736 usePixel[1208] = 0;
737 usePixel[583] = 0;
738 usePixel[830] = 0;
739 usePixel[1399] = 0;
740 usePixel[113] = 0;
741 usePixel[115] = 0;
742 usePixel[354] = 0;
743 usePixel[423] = 0;
744 usePixel[1195] = 0;
745 usePixel[1393] = 0;
746
747// usePixel[990] = 0;
748
749 // Map for which pixel which were suspicous
750 bool suspicous[1440];
751
752 for (int i=0; i<1440; i++) suspicous[i]=false;
753
754 // List of Pixel that should be ignored in camera view
755 // suspicous[802] = true;
756
757 //------------------------------------------------------------------------
758 // Fill and calculate sum spectrum
759
760 // Begin of Loop over Pixels
761 for (UInt_t pixel = 0; pixel < singles.GetNumPixels(); pixel++)
762 {
763 //jump to next pixel if the current is marked as faulty
764 if (usePixel[pixel]==0)
765 continue;
766
767 TH1D *hist = hsignal->ProjectionY("proj", pixel+1, pixel+1);
768 hist->SetDirectory(0);
769
770 //loop over slices
771 for (int b=1; b<=hist->GetNbinsX(); b++)
772 {
773 //Fill sum spectrum
774 hSum.Fill(hist->GetBinCenter(b), hist->GetBinContent(b));
775 }
776 }
777
778 for (UInt_t bin = 1; bin < (UInt_t)hSum.GetNbinsX()+1; bin++)
779 {
780 hSum.SetBinError(bin, hSum.GetBinContent(bin)*0.1);
781 }
782
783 gROOT->SetSelectedPad(0);
784 TCanvas &c11 = d->AddTab("SumHist");
785 c11.cd();
786 gPad->SetLogy();
787 gPad->SetGridx();
788 gPad->SetGridy();
789 hSum.DrawCopy("hist");
790 //------------------------fit sum spectrum-------------------------------
791 const Int_t maxbin = hSum.GetMaximumBin();
792 const Double_t binwidth = hSum.GetBinWidth(maxbin);
793 const Double_t ampl = hSum.GetBinContent(maxbin);
794
795 double fwhmSum = 0;
796
797 //Calculate full width half Maximum
798 for (int i=1; i<maxbin; i++)
799 if (hSum.GetBinContent(maxbin-i)+hSum.GetBinContent(maxbin+i)<ampl*2)
800 {
801 fwhmSum = (2*i+1)*binwidth;
802 break;
803 }
804
805 if (fwhmSum==0)
806 {
807 cout << "Could not determine start value for sigma." << endl;
808 }
809
810 //Set fit parameters
811 Double_t sum_par[6] =
812 {
813 ampl, hSum.GetBinCenter(maxbin), fwhmSum*1.4, 0.1, -10, 1
814 };
815 funcSum.SetParameters(sum_par);
816
817 //calculate fit range
818 const double min = sum_par[1]-fwhmSum*2;
819 const double max = sum_par[1]*(maxOrder);
820 funcSum.SetRange(min, max);
821 funcSum.FixParameter(5,0.1);
822
823 //Fit and draw spectrum
824 hSum.Fit(&funcSum, "NOQS", "", min, max);
825 funcSum.DrawCopy("SAME");
826 funcSum.GetParameters(sum_par);
827
828 //Set Parameters for following fits
829// const float Amplitude = sum_par[0];
830 const float gain = sum_par[1];
831// const float GainRMS = sum_par[2];
832 const float Crosstlk = sum_par[3];
833 const float Offset = sum_par[4];
834// const float PowCrosstlk = sum_par[5];
835
836
837 //--------fit gausses to peaks of sum specetrum -----------
838
839 // create histo for height of Maximum amplitudes vs. pe
840 double min_bin = hSum.GetBinCenter(maxbin);
841 double max_bin = hSum.GetBinCenter((maxOrder-2)*(maxbin));
842 double bin_width= (max_bin - min_bin)/10;
843
844 TH1D hMaxHeightsSum("MaxHeightSum", "peak heights",
845 10, min_bin-bin_width, max_bin*2-bin_width/2 );
846
847 FitGaussiansToSpectrum
848 (
849 maxOrder-2,
850 hSum,
851 hMaxHeightsSum,
852 maxbin,
853 fwhmSum,
854 gain
855 );
856
857 //EOF: fit sum spectrum-----------------------------
858
859 hMaxHeightsSum.SetLineColor(kRed);
860// hMaxHeightsSum.DrawCopy("SAME");
861
862 //Fit Peak heights
863 TF1 fExpo1( "fExpo1","expo", min, max);
864 fExpo1.SetLineColor(kRed);
865 hMaxHeightsSum.Fit(&fExpo1, "N0QS" );
866// hMaxHeightsSum.DrawCopy("SAME");
867// fExpo1.DrawCopy("SAME");
868
869 // EOF: Fill and calculate sum spectrum
870 //------------------------------------------------------------------------
871
872
873
874 //------------------------------------------------------------------------
875 // Gain Calculation for each Pixel
876
877 int count_ok = 0;
878
879 // Begin of Loop over Pixels
880 for (UInt_t pixel = 0; pixel < singles.GetNumPixels(); pixel++)
881 {
882
883 if (usePixel[pixel]==0)
884 continue;
885
886 //Projectipon of a certain Pixel out of Ampl.Spectrum
887 TH1D *hist = hsignal->ProjectionY("proj", pixel+1, pixel+1);
888 hist->SetDirectory(0);
889
890 for (UInt_t bin = 1; bin < (UInt_t)hSum.GetNbinsX()+1; bin++)
891 {
892 hist->SetBinError(bin, hSum.GetBinContent(bin)*0.1);
893 }
894
895 const Int_t maxbin = hist->GetMaximumBin();
896
897 const Double_t binwidth = hist->GetBinWidth(maxbin);
898
899 const Double_t ampl = hist->GetBinContent(maxbin);
900
901 const Double_t maxPos = hist->GetBinCenter(maxbin);
902
903 double fwhm = 0;
904
905 for (int i=1; i<maxbin; i++)
906 if (hist->GetBinContent(maxbin-i)+hist->GetBinContent(maxbin+i)<ampl*2)
907 {
908 fwhm = (2*i+1)*binwidth;
909 break;
910 }
911
912 if (fwhm==0)
913 {
914 cout << "Could not determine start value for sigma." << endl;
915 continue;
916 }
917
918 // Amplitude, Gain, Sigma, Crosstalk Prob, Shift
919 // par[1] + par[4] is the position of the first maximum
920 Double_t par[5] =
921 {
922 ampl,
923 hist->GetBinCenter(maxbin),
924// gain,
925 fwhm*2,
926// GainRMS,
927 Crosstlk,
928 Offset
929 };
930
931
932 func.SetParameters(par);
933 const double fit_min = par[1]-fwhm*1.4;
934 const double fit_max = par[1]*maxOrder;
935 func.SetRange(fit_min-fwhm*2, fit_max);
936 func.SetParLimits(1, maxPos - fwhm, maxPos + fwhm);
937
938 if (suspicous[pixel] == true)
939 {
940 cout << "Maxbin\t" << maxbin << endl;
941 cout << "min\t" << fit_min << endl;
942 cout << "max\t" << fit_max << endl;
943 cout << "Amplitude, 1.Max, Sigma, fwhm:\t"
944 << par[0] << "\t"
945 << par[1] << "\t"
946 << par[2] << "\t"
947 << fwhm << "\t" << endl;
948 }
949
950 //Rebin Projection
951 hist->Rebin(2);
952
953 //Fit Pixels spectrum
954 const TFitResultPtr rc = hist->Fit(&func, "N0QS", "", fit_min, fit_max);
955
956 const bool ok = int(rc)==0;
957
958 const double fit_prob = rc->Prob();
959 const float fGain = func.GetParameter(1);
960 const float fAmplitude = func.Integral(0, 7*fGain)/h->GetNumEntries(); //func.GetParameter(0);
961 const float f2GainRMS = func.GetParameter(2);
962 const float fCrosstlk = func.GetParameter(3);
963 const float fOffset = func.GetParameter(4)/singles.GetIntegrationWindow();
964// const float fCrossOrd = func.GetParameter(5);
965
966 hAmplitude.Fill(fAmplitude);
967 hGain.Fill(fGain);
968 //hGainRMS.Fill(f2GainRMS);
969 hCrosstlk.Fill(fCrosstlk*100);
970 hOffset.Fill(fOffset);
971 hFitProb.Fill(100*fit_prob);
972 hGainRMS.Fill(100*f2GainRMS/fGain);
973
974 hcam.SetBinContent(pixel+1, fGain);
975 hcam2.SetBinContent(pixel+1, 100*f2GainRMS/fGain);
976
977 usePixel[pixel] = ok;
978
979 if (!ok)
980 {
981 cout << "Pixel " << setw(4) << pixel << ": failed!" << endl;
982 continue;
983 }
984
985 //Fill sum spectrum
986 for (int b=1; b<=hist->GetNbinsX(); b++)
987 {
988// hSum.Fill(hist->GetBinCenter(b), hist->GetBinContent(b));
989 hSumScale.Fill((hist->GetBinCenter(b)-fOffset)/fGain, (hist->GetBinContent(b)));
990 }
991
992 //plott special pixel
993 if (
994 count_ok==0 ||
995 suspicous[pixel]
996 )
997 {
998 TCanvas &c = d->AddTab(Form("Pix%d", pixel));
999 c.cd();
1000 gPad->SetLogy();
1001 gPad->SetGridx();
1002 gPad->SetGridy();
1003
1004 hist->SetStats(false);
1005 hist->SetYTitle("Counts [a.u.]");
1006 hist->SetName(Form("Pix%d", pixel));
1007 hist->DrawCopy("hist")->SetDirectory(0);
1008// hist->Draw();
1009 func.DrawCopy("SAME");
1010 }
1011
1012 count_ok++;
1013
1014 delete hist;
1015// if (suspicous[pixel] == true) usePixel[pixel]=0;
1016 }
1017
1018 // End of Loop over Pixel
1019
1020 //------------------fill histograms-----------------------
1021
1022 hcam.SetUsed(usePixel);
1023 hcam2.SetUsed(usePixel);
1024
1025 TCanvas &canv = d->AddTab("Gain");
1026 canv.cd();
1027 canv.Divide(2,2);
1028
1029 canv.cd(1);
1030 hcam.SetNameTitle( "gain","Gain distribution over whole camera");
1031 hcam.DrawCopy();
1032
1033 canv.cd(2);
1034 hcam2.SetNameTitle( "gainRMS","GainRMS distribution over whole camera");
1035 hcam2.DrawCopy();
1036
1037 TCanvas &cam_canv = d->AddTab("Camera_Gain");
1038
1039 //------------------ Draw gain corrected sum specetrum -------------------
1040 gROOT->SetSelectedPad(0);
1041 TCanvas &c12 = d->AddTab("GainHist");
1042 c12.cd();
1043 gPad->SetLogy();
1044 gPad->SetGridx();
1045 gPad->SetGridy();
1046
1047 for (UInt_t bin = 1; bin < (UInt_t)hSum.GetNbinsX()+1; bin++)
1048 {
1049 hSumScale.SetBinError(bin, hSumScale.GetBinContent(bin)*0.1);
1050 }
1051
1052 hSumScale.DrawCopy("hist");
1053
1054
1055 //-------------------- fit gain corrected sum spectrum -------------------
1056 const Int_t maxbin2 = hSumScale.GetMaximumBin();
1057 const Double_t binwidth2 = hSumScale.GetBinWidth(maxbin2);
1058 const Double_t ampl2 = hSumScale.GetBinContent(maxbin2);
1059
1060 //Calculate full width half Maximum
1061 double fwhmScaled = 0;
1062
1063 for (int i=1; i<maxbin; i++)
1064 if (hSumScale.GetBinContent(maxbin2-i)+hSumScale.GetBinContent(maxbin2+i)<ampl2*2)
1065 {
1066 fwhmScaled = (2*i+1)*binwidth2;
1067 break;
1068 }
1069
1070 if (fwhmScaled==0)
1071 {
1072 cout << "Could not determine start value for sigma." << endl;
1073 }
1074
1075 //Set fit parameters
1076 Double_t par2[6] =
1077 {
1078 ampl2, 0.9, fwhmScaled*1.4, Crosstlk, 0, 1.1
1079 };
1080
1081 funcScaled.SetParameters(par2);
1082// funcScaledBL.SetParameters(par2);
1083// funcScaled.FixParameter(1,0.9);
1084 funcScaled.FixParameter(4,0);
1085// funcScaledBL.FixParameter(4,0);
1086 funcScaled.FixParameter(5,Crosstlk);
1087// funcScaledBL.FixParameter(5,Crosstlk);
1088
1089 const double minScaled = par2[1]-fwhmScaled;
1090 const double maxScaled = par2[1]*maxOrder;
1091 cout << "par2[1]" << par2[1] << endl;
1092 cout << "maxScaled\t" << maxScaled
1093 << " \t\t minScaled\t" << minScaled << endl;
1094
1095 funcScaled.SetRange(minScaled, maxScaled);
1096// funcScaledBL.SetRange(minScaled, maxScaled);
1097 hSumScale.Fit(&funcScaled, "N0QS", "", minScaled, maxScaled);
1098// hSumScale.Fit(&funcScaledBL, "WLN0QS", "", minScaled, maxScaled);
1099 funcScaled.DrawCopy("SAME");
1100// funcScaledBL.DrawCopy("SAME");
1101 //--------fit gausses to peaks of gain corrected sum specetrum -----------
1102
1103 // create histo for height of Maximum amplitudes vs. pe
1104 TH1D hMaxHeights("MaxHeights", "peak heights",
1105 10, hSumScale.GetBinCenter(maxbin-1)-0.5, 10+hSumScale.GetBinCenter(maxbin-1)-0.5);
1106
1107 FitGaussiansToSpectrum
1108 (
1109 maxOrder-2,
1110 hSumScale,
1111 hMaxHeights,
1112 maxbin2,
1113 fwhmScaled,
1114 1
1115 );
1116
1117
1118
1119// TCanvas &c16 = d->AddTab("PeakHeights");
1120// gPad->SetLogy();
1121 hMaxHeights.SetLineColor(kRed);
1122// hMaxHeights.DrawCopy("SAME");
1123 TF1 fExpo( "fExpo","expo", 0, maxOrder-2);
1124 fExpo.SetLineColor(kRed);
1125 hMaxHeights.Fit(&fExpo, "N0QS" );
1126// fExpo.DrawCopy("SAME");
1127// c12.cd();
1128// fExpo.DrawCopy("SAME");
1129
1130 TCanvas &c2 = d->AddTab("Result");
1131 c2.Divide(3,2);
1132 c2.cd(1);
1133 gPad->SetLogy();
1134 hAmplitude.DrawCopy();
1135
1136 c2.cd(2);
1137 gPad->SetLogy();
1138// hGain.DrawCopy();
1139
1140
1141 TF1 GainGaus( "GainGaus", "gaus", 0, hGain.GetXaxis()->GetXmax());
1142 hGain.Rebin(2);
1143 hGain.Fit(&GainGaus, "N0QS");
1144
1145// float gain_mean = GainGaus.GetParameter(1);
1146 float gain_median = MedianOfH1(&hGain);
1147 TH1F hNormGain ("NormGain", "median normalized Gain distribution",
1148 hGain.GetNbinsX(),
1149 (hGain.GetXaxis()->GetXmin())/gain_median,
1150 (hGain.GetXaxis()->GetXmax())/gain_median
1151 );
1152 hNormGain.SetXTitle("gain [fraction of median gain value]");
1153 hNormGain.SetYTitle("Counts");
1154
1155
1156 hNormGain.Add(&hGain);
1157// hNormGain.SetStats(false);
1158 hNormGain.GetXaxis()->SetRangeUser(
1159 1-(hNormGain.GetXaxis()->GetXmax()-1),
1160 hNormGain.GetXaxis()->GetXmax()
1161 );
1162 hNormGain.DrawCopy();
1163 hNormGain.Fit(&GainGaus, "N0QS");
1164 GainGaus.DrawCopy("SAME");
1165
1166 //plot normailzed camera view
1167 cam_canv.cd();
1168 MHCamera hNormCam(fact);
1169 hNormCam.SetStats(false);
1170
1171 for (UInt_t pixel =1; pixel <= 1440; pixel++)
1172 {
1173 hNormCam.SetBinContent(pixel, hcam.GetBinContent(pixel)/gain_median);
1174 }
1175
1176 hNormCam.SetNameTitle( "gain","normalized [Median] Gain distribution over whole camera");
1177 hNormCam.SetZTitle("Horst");
1178 hNormCam.SetUsed(usePixel);
1179 hNormCam.DrawCopy();
1180// hGain.Scale(1/GainGaus.GetParameter(1));
1181// GainGaus.DrawCopy("SAME");
1182
1183 c2.cd(3);
1184 gPad->SetLogy();
1185 hGainRMS.DrawCopy();
1186
1187 c2.cd(4);
1188 gPad->SetLogy();
1189 hCrosstlk.DrawCopy();
1190
1191 c2.cd(5);
1192 gPad->SetLogy();
1193 hOffset.DrawCopy();
1194
1195 c2.cd(6);
1196 gPad->SetLogy();
1197// hFitProb.DrawCopy();
1198 hGain.DrawCopy();
1199
1200/*
1201 TCanvas &c25 = d->AddTab("Spectra");
1202 c25.Divide(2,1);
1203 c25.cd(1);
1204 gPad->SetLogy();
1205 gPad->SetGrid();
1206 hSum.DrawCopy("hist");
1207 funcSum.DrawCopy("SAME");
1208
1209 c25.cd(2);
1210 gPad->SetLogy();
1211 gPad->SetGrid();
1212 hSumScale.DrawCopy("hist");
1213 funcScaled.DrawCopy("SAME");
1214 funcScaledBL.DrawCopy("SAME");
1215*/
1216
1217 /*
1218 TCanvas &c3 = d->AddTab("Separation");
1219 gPad->SetLogy();
1220 hSep.DrawCopy();
1221 */
1222
1223 d->SaveAs(outfile);
1224 return 0;
1225}
1226
1227Double_t fcn(Double_t *xx, Double_t *par)
1228{
1229 Double_t x = xx[0];
1230
1231 Double_t ampl = par[0];
1232 Double_t gain = par[1];
1233 Double_t sigma = par[2];
1234 Double_t cross = par[3];
1235 Double_t shift = par[4];
1236
1237 Double_t xTalk = 1;
1238
1239 Double_t y = 0;
1240
1241 for (int order = 1; order < 14; order++)
1242 {
1243 Double_t arg = (x - order*gain - shift)/sigma;
1244
1245 Double_t gauss = TMath::Exp(-0.5 * arg * arg/order);
1246
1247 y += xTalk*gauss;
1248
1249 xTalk *= cross;
1250
1251 }
1252
1253 return ampl*y;
1254}
1255
1256Double_t fcn_cross(Double_t *xx, Double_t *par)
1257{
1258 Double_t x = xx[0];
1259
1260 Double_t ampl = par[0];
1261 Double_t gain = par[1];
1262 Double_t sigma = par[2];
1263 Double_t cross = par[3];
1264 Double_t shift = par[4];
1265// Double_t powOfX = par[5];
1266
1267 Double_t xTalk = 1;
1268
1269 Double_t y = 0;
1270
1271 for (int order = 1; order < 14; order++)
1272 {
1273 Double_t arg = (x - order*gain - shift)/sigma;
1274
1275 Double_t gauss = TMath::Exp(-0.5 * arg * arg/order);
1276
1277 y += xTalk*gauss;
1278
1279// xTalk = pow(cross, order*powOfX);
1280// for (int j = 1; j < powOfX; j++)
1281// {
1282// xTalk *= cross;
1283// }
1284// cross *= TMath::Exp(-powOfX*cross);
1285 xTalk *= cross;
1286 }
1287
1288 return ampl*y;
1289}
1290
1291
1292void
1293FitGaussiansToSpectrum(
1294 UInt_t maxOrder,
1295 TH1 &hSource,
1296 TH1 &hDest,
1297 Int_t maxbin,
1298 double fwhm,
1299 Double_t gain
1300)
1301{
1302
1303 Double_t peakposition = hSource.GetBinCenter(maxbin);
1304 char fname[20];
1305
1306 // fit gauss functions to spectrum
1307 for (UInt_t nr = 1; nr<maxOrder; nr++)
1308 {
1309 sprintf(fname,"gaus%d",nr);
1310
1311 //create gauss functions
1312 TF1 gaus( fname,"gaus", peakposition-fwhm, peakposition+fwhm);
1313 gaus.SetParLimits(1, peakposition-fwhm, peakposition+fwhm);
1314 gaus.SetParLimits(2, -fwhm, fwhm );
1315 //fit and draw gauss functions
1316 hSource.Fit( &gaus, "N0QS", "", peakposition-fwhm, peakposition+fwhm);
1317// gaus.DrawCopy("SAME");
1318
1319 peakposition = gaus.GetParameter(1);
1320
1321 //fill histo for height of Maximum amplitudes vs. pe
1322 hDest.SetBinContent(nr, gaus.GetParameter(0));
1323
1324 peakposition += gain;
1325 }
1326
1327 return;
1328}
1329
1330Double_t MedianOfH1 (
1331 TH1* inputHisto
1332)
1333{
1334 //compute the median for 1-d histogram h1
1335 Int_t nbins = inputHisto->GetXaxis()->GetNbins();
1336 Double_t *x = new Double_t[nbins];
1337 Double_t *y = new Double_t[nbins];
1338
1339 for (Int_t i=0; i<nbins; i++)
1340 {
1341 x[i] = inputHisto->GetXaxis()->GetBinCenter(i+1);
1342 y[i] = inputHisto->GetBinContent(i+1);
1343 }
1344
1345 Double_t median = TMath::Median(nbins,x,y);
1346 delete [] x;
1347 delete [] y;
1348 return median;
1349}
1350// end of PlotMedianEachSliceOfPulse
1351//----------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.