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

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