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

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