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

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