source: trunk/Mars/fact/analysis/gain/extract_singles.C@ 17033

Last change on this file since 17033 was 17033, checked in by tbretz, 11 years ago
First version of files.
File size: 24.1 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 <TGraph.h>
9#include <TFitResultPtr.h>
10#include <TFitResult.h>
11#include <TFile.h>
12
13#include <cstdio>
14#include <stdio.h>
15#include <stdint.h>
16
17#include "MH.h"
18#include "MArrayI.h"
19#include "MLog.h"
20#include "MLogManip.h"
21#include "MDirIter.h"
22#include "MFillH.h"
23#include "MEvtLoop.h"
24#include "MCamEvent.h"
25#include "MHCamEvent.h"
26#include "MGeomApply.h"
27#include "MTaskList.h"
28#include "MParList.h"
29#include "MContinue.h"
30#include "MBinning.h"
31#include "MDrsCalibApply.h"
32#include "MDrsCalibration.h"
33#include "MRawFitsRead.h"
34#include "MBadPixelsCam.h"
35#include "MStatusDisplay.h"
36#include "MTaskInteractive.h"
37#include "MPedestalSubtractedEvt.h"
38#include "MHCamera.h"
39#include "MGeomCamFACT.h"
40#include "MRawRunHeader.h"
41#include "MPedestalCam.h"
42#include "MPedestalPix.h"
43#include "MParameters.h"
44
45using namespace std;
46using namespace TMath;
47
48// Structure to store the extracted value and the extracted time
49struct Single
50{
51 float fSignal;
52 float fTime;
53};
54
55//Storage class to kKeep a list of the extracted single
56class MSingles : public MParContainer, public MCamEvent
57{
58 Int_t fIntegrationWindow;
59
60 vector<vector<Single>> fData;
61
62public:
63 MSingles(const char *name=NULL, const char *title=NULL) : fIntegrationWindow(30)
64 {
65 fName = name ? name : "MSingles";
66 fName = title ? title : "Storeage for vector of singles";
67 }
68
69 void InitSize(const UInt_t i)
70 {
71 fData.resize(i);
72 }
73
74 vector<Single> &operator[](UInt_t i) { return fData[i]; }
75 vector<Single> &GetVector(UInt_t i) { return fData[i]; }
76
77 UInt_t GetNumPixels() const { return fData.size(); }
78
79 void SetIntegrationWindow(Int_t w) { fIntegrationWindow = w; }
80 Int_t GetIntegrationWindow() const { return fIntegrationWindow; }
81
82 Bool_t GetPixelContent(Double_t &, Int_t , const MGeomCam &, Int_t) const
83 {
84 return kTRUE;
85 }
86 void DrawPixelContent(Int_t) const { }
87
88 ClassDef(MSingles, 1)
89};
90
91// Histogram class to extract the baseline
92class MHBaseline : public MH
93{
94 TH2F fBaseline;
95
96 MPedestalCam *fPedestal;
97
98 // The baseline is only extracted where also the signal is extracted
99 // FIXME: Make sure this is consistent with MExtractSingles
100 UShort_t fSkipStart;
101 UShort_t fSkipEnd;
102
103public:
104 MHBaseline() : fPedestal(0), fSkipStart(20), fSkipEnd(10)
105 {
106 fName = "MHBaseline";
107
108 // Setup the histogram
109 fBaseline.SetName("Baseline");
110 fBaseline.SetTitle("Median spectrum");
111 fBaseline.SetXTitle("Pixel [idx]");
112 fBaseline.SetYTitle("Median baseline [mV]");
113 fBaseline.SetDirectory(NULL);
114 }
115
116 Bool_t ReInit(MParList *plist)
117 {
118 fPedestal = (MPedestalCam*)plist->FindCreateObj("MPedestalCam");
119 if (!fPedestal)
120 return kFALSE;
121
122 const MRawRunHeader *header = (MRawRunHeader*)plist->FindObject("MRawRunHeader");
123 if (!header)
124 {
125 *fLog << err << "MRawRunHeader not found... abort." << endl;
126 return kFALSE;
127 }
128
129 // Bin width should be around 1 dac count which is about 0.5mV
130 MBinning binsx, binsy;
131 binsx.SetEdges(header->GetNumNormalPixels(), -0.5, header->GetNumNormalPixels()-0.5);
132 binsy.SetEdges(100, -20.5, 29.5);
133
134 // Setup binnning
135 MH::SetBinning(fBaseline, binsx, binsy);
136
137 return kTRUE;
138 }
139
140 // Fill the samples into the histogram
141 Int_t Fill(const MParContainer *par, const Stat_t)
142 {
143 const MPedestalSubtractedEvt *evt = dynamic_cast<const MPedestalSubtractedEvt*>(par);
144
145 const Int_t n = evt->GetNumSamples()-fSkipStart-fSkipEnd;
146
147 // Loop over all pixels
148 for (int pix=0; pix<evt->GetNumPixels(); pix++)
149 {
150 // Get samples for each pixel
151 const Float_t *ptr = evt->GetSamples(pix);
152
153 // Average two consecutive samples
154 for (int i=0; i<n; i+=2)
155 {
156 const Double_t val = 0.5*ptr[i+fSkipStart]+0.5*ptr[i+1+fSkipStart];
157 fBaseline.Fill(pix, val);
158 }
159 }
160
161 return kTRUE;
162 }
163
164 // Extract the baseline value from the distrbutions
165 Bool_t Finalize()
166 {
167 fPedestal->InitSize(fBaseline.GetNbinsX());
168 fPedestal->SetNumEvents(GetNumExecutions());
169
170 Int_t cnt = 0;
171 Double_t avg = 0;
172 Double_t rms = 0;
173
174 // Loop over all 'pixels'
175 for (int x=0; x<fBaseline.GetNbinsX(); x++)
176 {
177 // Get the corresponding slice from the histogram
178 TH1D *hist = fBaseline.ProjectionY("proj", x+1, x+1);
179
180 // Get the maximum bin
181 const Int_t bin = hist->GetMaximumBin();
182
183 // Calculate a parabola through this and the surrounding points
184 // on logarithmic values (that's a gaussian)
185
186 //
187 // Quadratic interpolation
188 //
189 // calculate the parameters of a parabula such that
190 // y(i) = a + b*x(i) + c*x(i)^2
191 // y'(i) = b + 2*c*x(i)
192 //
193 //
194
195 // -------------------------------------------------------------------------
196 // a = y2;
197 // b = (y3-y1)/2;
198 // c = (y3+y1)/2 - y2;
199
200 const Double_t v1 = hist->GetBinContent(bin-1);
201 const Double_t v2 = hist->GetBinContent(bin);
202 const Double_t v3 = hist->GetBinContent(bin+1);
203 if (v1<=0 || v2<=0 || v3<=0)
204 continue;
205
206 const Double_t y1 = TMath::Log(v1);
207 const Double_t y2 = TMath::Log(v2);
208 const Double_t y3 = TMath::Log(v3);
209
210 //Double_t a = y2;
211 const Double_t b = (y3-y1)/2;
212 const Double_t c = (y1+y3)/2 - y2;
213 if (c>=0)
214 continue;
215
216 const Double_t w = -1./(2*c);
217 const Double_t dx = b*w;
218
219 if (dx<-1 || dx>1)
220 continue;
221
222 // y = exp( - (x-k)^2 / s^2 / 2 )
223 //
224 // -2*s^2 * log(y) = x^2 - 2*k*x + k^2
225 //
226 // a = (k/s0)^2/2
227 // b = k/s^2
228 // c = -1/(2s^2) --> s = sqrt(-1/2c)
229
230 const Double_t binx = hist->GetBinCenter(bin);
231 const Double_t binw = hist->GetBinWidth(bin);
232
233 //const Double_t p = hist->GetBinCenter(hist->GetMaximumBin());
234 const Double_t p = binx + dx*binw;
235
236 avg += p;
237 rms += p*p;
238
239 cnt++;
240
241 // Store baseline and sigma
242 MPedestalPix &pix = (*fPedestal)[x];
243
244 pix.SetPedestal(p);
245 pix.SetPedestalRms(TMath::Sqrt(w)*binw);
246
247 delete hist;
248 }
249
250 avg /= cnt;
251 rms /= cnt;
252
253 cout << "Baseline(" << cnt << "): " << avg << " +- " << sqrt(rms-avg*avg) << endl;
254
255 return kTRUE;
256 }
257
258 // Draw histogram
259 void Draw(Option_t *)
260 {
261 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
262
263 AppendPad("");
264
265 pad->SetBorderMode(0);
266 pad->SetFrameBorderMode(0);
267 fBaseline.Draw("colz");
268 }
269
270
271 ClassDef(MHBaseline, 1);
272};
273
274// Histogram class for the signal and time distribution as
275// well as the pulse shape
276class MHSingles : public MH
277{
278 TH2F fSignal;
279 TH2F fTime;
280 TProfile2D fPulse;
281
282 UInt_t fNumEntries;
283
284 MSingles *fSingles; //!
285 MPedestalSubtractedEvt *fEvent; //!
286 MBadPixelsCam *fBadPix; //!
287
288public:
289 MHSingles() : fNumEntries(0), fSingles(0), fEvent(0)
290 {
291 fName = "MHSingles";
292
293 // Setup histograms
294 fSignal.SetName("Signal");
295 fSignal.SetTitle("pulse integral spectrum");
296 fSignal.SetXTitle("pixel [SoftID]");
297 fSignal.SetYTitle("time [a.u.]");
298 fSignal.SetDirectory(NULL);
299
300 fTime.SetName("Time");
301 fTime.SetTitle("arival time spectrum");
302 fTime.SetXTitle("pixel [SoftID]");
303 fTime.SetYTitle("time [a.u.]");
304 fTime.SetDirectory(NULL);
305
306 fPulse.SetName("Pulse");
307 fPulse.SetTitle("average pulse shape spectrum");
308 fPulse.SetXTitle("pixel [SoftID]");
309 fPulse.SetYTitle("time [a.u.]");
310 fPulse.SetDirectory(NULL);
311 }
312
313 Bool_t SetupFill(const MParList *plist)
314 {
315 fSingles = (MSingles*)plist->FindObject("MSingles");
316 if (!fSingles)
317 {
318 *fLog << err << "MSingles not found... abort." << endl;
319 return kFALSE;
320 }
321
322 fEvent = (MPedestalSubtractedEvt*)plist->FindObject("MPedestalSubtractedEvt");
323 if (!fEvent)
324 {
325 *fLog << err << "MPedestalSubtractedEvt not found... abort." << endl;
326 return kFALSE;
327 }
328
329 fBadPix = (MBadPixelsCam*)plist->FindObject("MBadPixelsCam");
330 if (!fBadPix)
331 {
332 *fLog << err << "MBadPixelsCam not found... abort." << endl;
333 return kFALSE;
334 }
335
336 fNumEntries = 0;
337
338 return kTRUE;
339 }
340
341 Bool_t ReInit(MParList *plist)
342 {
343 const MRawRunHeader *header = (MRawRunHeader*)plist->FindObject("MRawRunHeader");
344 if (!header)
345 {
346 *fLog << err << "MRawRunHeader not found... abort." << endl;
347 return kFALSE;
348 }
349
350 // Setup binning
351 const Int_t w = fSingles->GetIntegrationWindow();
352
353 MBinning binsx, binsy;
354 binsx.SetEdges(fSingles->GetNumPixels(), -0.5, fSingles->GetNumPixels()-0.5);
355 binsy.SetEdges(22*w, -10*w, 100*w);
356
357 MH::SetBinning(fSignal, binsx, binsy);
358
359 const UShort_t roi = header->GetNumSamples();
360
361 // Correct for DRS timing!!
362 MBinning binst(roi, -0.5, roi-0.5);
363 MH::SetBinning(fTime, binsx, binst);
364
365 MBinning binspy(2*roi, -roi-0.5, roi-0.5);
366 MH::SetBinning(fPulse, binsx, binspy);
367
368 return kTRUE;
369 }
370
371 // Fill singles into histograms
372 Int_t Fill(const MParContainer *, const Stat_t)
373 {
374 // Get pointer to samples to fill samples
375 const Float_t *ptr = fEvent->GetSamples(0);
376 const Short_t roi = fEvent->GetNumSamples();
377
378 // Loop over all pixels
379 for (unsigned int i=0; i<fSingles->GetNumPixels(); i++)
380 {
381 if ((*fBadPix)[i].IsUnsuitable())
382 continue;
383
384 // loop over all singles
385 const vector<Single> &result = fSingles->GetVector(i);
386
387 for (unsigned int j=0; j<result.size(); j++)
388 {
389 // Fill signal and time
390 fSignal.Fill(i, result[j].fSignal);
391 fTime.Fill (i, result[j].fTime);
392
393 if (!ptr)
394 continue;
395
396 // Fill pulse shape
397 const Short_t tm = floor(result[j].fTime);
398
399 for (int s=0; s<roi; s++)
400 fPulse.Fill(i, s-tm, ptr[s]);
401 }
402
403 ptr += roi;
404 }
405
406 fNumEntries++;
407
408 return kTRUE;
409 }
410
411 // Getter for histograms
412 const TH2 *GetSignal() const { return &fSignal; }
413 const TH2 *GetTime() const { return &fTime; }
414 const TH2 *GetPulse() const { return &fPulse; }
415
416 UInt_t GetNumEntries() const { return fNumEntries; }
417
418 void Draw(Option_t *)
419 {
420 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
421
422 AppendPad("");
423
424 pad->Divide(2,2);
425
426 pad->cd(1);
427 gPad->SetBorderMode(0);
428 gPad->SetFrameBorderMode(0);
429 fSignal.Draw("colz");
430
431 pad->cd(2);
432 gPad->SetBorderMode(0);
433 gPad->SetFrameBorderMode(0);
434 fTime.Draw("colz");
435
436 pad->cd(3);
437 gPad->SetBorderMode(0);
438 gPad->SetFrameBorderMode(0);
439 fPulse.Draw("colz");
440 }
441
442 void DrawCopy(Option_t *)
443 {
444 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
445
446 AppendPad("");
447
448 pad->Divide(2,2);
449
450 pad->cd(1);
451 gPad->SetBorderMode(0);
452 gPad->SetFrameBorderMode(0);
453 fSignal.DrawCopy("colz");
454
455 pad->cd(2);
456 gPad->SetBorderMode(0);
457 gPad->SetFrameBorderMode(0);
458 fTime.DrawCopy("colz");
459
460 pad->cd(3);
461 gPad->SetBorderMode(0);
462 gPad->SetFrameBorderMode(0);
463 fPulse.DrawCopy("colz");
464 }
465
466 ClassDef(MHSingles, 1)
467};
468
469// Task to extract the singles
470class MExtractSingles : public MTask
471{
472 MSingles *fSingles;
473 MPedestalCam *fPedestal;
474 MPedestalSubtractedEvt *fEvent;
475
476 // On time for each pixel in samples
477 MArrayI fExtractionRange;
478
479 // Number of samples for sliding average
480 size_t fNumAverage;
481
482 // The range in which the singles have to fit in
483 // FIXME: Make sure this is consistent with MHBaseline
484 UInt_t fStartSkip;
485 UInt_t fEndSkip;
486
487 UInt_t fIntegrationSize;
488 UInt_t fMaxSearchWindow;
489
490 Int_t fMaxDist;
491 Double_t fThreshold;
492
493public:
494 MExtractSingles() : fSingles(0), fPedestal(0), fEvent(0),
495 fExtractionRange(1440),
496 fNumAverage(10), fStartSkip(20), fEndSkip(10),
497 fIntegrationSize(3*10), fMaxSearchWindow(30)
498 {
499 }
500
501 void SetMaxDist(Int_t m) { fMaxDist = m; }
502 void SetThreshold(Int_t th) { fThreshold = th; }
503
504 UInt_t GetIntegrationSize() const { return fIntegrationSize; }
505
506 const MArrayI &GetExtractionRange() const { return fExtractionRange; }
507
508
509 Int_t PreProcess(MParList *plist)
510 {
511 fSingles = (MSingles*)plist->FindCreateObj("MSingles");
512 if (!fSingles)
513 return kFALSE;
514
515 fEvent = (MPedestalSubtractedEvt*)plist->FindObject("MPedestalSubtractedEvt");
516 if (!fEvent)
517 {
518 *fLog << err << "MPedestalSubtractedEvt not found... abort." << endl;
519 return kFALSE;
520 }
521
522 fPedestal = (MPedestalCam*)plist->FindObject("MPedestalCam");
523 if (!fPedestal)
524 {
525 *fLog << err << "MPedestalCam not found... abort." << endl;
526 return kFALSE;
527 }
528
529 return kTRUE;
530 }
531
532 Int_t Process()
533 {
534 const UInt_t roi = fEvent->GetNumSamples();
535
536 const size_t navg = fNumAverage;
537 const float threshold = fThreshold;
538 const UInt_t start_skip = fStartSkip;
539 const UInt_t end_skip = fEndSkip;
540 const UInt_t integration_size = fIntegrationSize;
541 const UInt_t max_search_window = fMaxSearchWindow;
542 const UInt_t max_dist = fMaxDist;
543
544 vector<float> val(roi-navg);//result of first sliding average
545 for (int pix=0; pix<fEvent->GetNumPixels(); pix++)
546 {
547 // Total number of samples as 'on time'
548 fExtractionRange[pix] += roi-start_skip-navg-end_skip-integration_size;
549
550 // Clear singles for this pixel
551 vector<Single> &result = fSingles->GetVector(pix);
552 result.clear();
553
554 // Get pointer to samples
555 const Float_t *ptr = fEvent->GetSamples(pix);
556
557 // Get Baseline
558 const Double_t ped = (*fPedestal)[pix].GetPedestal();
559
560 // Subtract baseline and do a sliding average over
561 // the samples to remove noise (mainly the 200MHz noise)
562 for (size_t i=0; i<roi-navg; i++)
563 {
564 val[i] = 0;
565 for (size_t j=i; j<i+navg; j++)
566 val[i] += ptr[j];
567 val[i] /= navg;
568 val[i] -= ped;
569 }
570
571 // peak finding via threshold
572 UInt_t imax = roi-navg-end_skip-integration_size;
573 for (UInt_t i=start_skip; i<imax; i++)
574 {
575 //search for threshold crossings
576 if (val[i+0]>threshold ||
577 val[i+4]>threshold ||
578
579 val[i+5]<threshold ||
580 val[i+9]<threshold)
581 continue;
582
583 //search for maximum after threshold crossing
584 UInt_t k_max = i+5;
585 for (UInt_t k=i+6; k<i+max_search_window; k++)
586 {
587 if (val[k] > val[k_max])
588 k_max = k;
589 }
590
591 // Check if the findings make sense
592 if (k_max == i+5 || k_max == i + max_search_window-1)
593 continue;
594
595 //search for half maximum before maximum
596 UInt_t k_half_max = 0;
597 for (UInt_t k=k_max; k>k_max-25; k--)
598 {
599 if (val[k-1] < val[k_max]/2 &&
600 val[k] >= val[k_max]/2 )
601 {
602 k_half_max = k;
603 break;
604 }
605 }
606 // Check if the finding makes sense
607 if (k_half_max+integration_size > roi-navg-end_skip)
608 continue;
609 if (k_half_max == 0)
610 continue;
611 if (k_max - k_half_max > max_dist)
612 continue;
613
614 // FIXME: Evaluate arrival time more precisely!!!
615 // FIXME: Get a better integration window
616
617 // Compile "single"
618 Single single;
619 single.fSignal = 0;
620 single.fTime = k_half_max;
621
622 // Crossing of the threshold is at 5
623 for (UInt_t j=k_half_max; j<k_half_max+integration_size; j++)
624 single.fSignal += ptr[j];
625
626 single.fSignal -= integration_size*ped;
627
628 // Add single to list
629 result.push_back(single);
630
631 // skip falling edge
632 i += 5+integration_size;
633
634 // Remove skipped samples from 'on time'
635 fExtractionRange[pix] -= i>imax ? 5+integration_size-(i-imax) : 5+integration_size;
636 }
637 }
638 return kTRUE;
639 }
640};
641
642int extract_singles(
643 Int_t max_dist = 14,
644 Double_t threshold = 5,
645 const char *runfile = "/fact/raw/2012/07/23/20120723_",
646 int firstRunID = 6,
647 int lastRunID = 6,
648 const char *drsfile = "/fact/raw/2012/07/23/20120723_004.drs.fits.gz",
649 const char *outpath = "."
650 )
651{
652 // ======================================================
653
654 MDirIter iter;
655
656 // built outpu filename and path
657 TString outfile = Form("%s/", outpath);
658 outfile += gSystem->BaseName(runfile);
659 outfile += Form("%03d_%03d.root", firstRunID, lastRunID);
660
661
662 if (!gSystem->AccessPathName(outfile))
663 return 0;
664
665 // add input files to directory iterator
666 for (Int_t fileNr = firstRunID; fileNr <= lastRunID; fileNr++)
667 {
668 TString filename = runfile;
669 filename += Form("%03d.fits", fileNr);
670 iter.AddDirectory(gSystem->DirName(runfile), gSystem->BaseName(filename+".fz"));
671 iter.AddDirectory(gSystem->DirName(runfile), gSystem->BaseName(filename+".gz"));
672 }
673
674 // ======================================================
675
676 // true: Display correctly mapped pixels in the camera displays
677 // but the value-vs-index plot is in software/spiral indices
678 // false: Display pixels in hardware/linear indices,
679 // but the order is the camera display is distorted
680 bool usemap = true;
681
682 // map file to use (get that from La Palma!)
683 const char *map = usemap ? "/scratch/fact/FACTmap111030.txt" : NULL;
684
685 // ------------------------------------------------------
686
687 Long_t max0 = 0;
688 Long_t max1 = 0;
689
690 // ======================================================
691
692 if (map && gSystem->AccessPathName(map, kFileExists))
693 {
694 gLog << "ERROR - Cannot access mapping file '" << map << "'" << endl;
695 return 1;
696 }
697 if (gSystem->AccessPathName(drsfile, kFileExists))
698 {
699 gLog << "ERROR - Cannot access drsfile file '" << drsfile << "'" << endl;
700 return 2;
701 }
702
703 // --------------------------------------------------------------------------------
704
705 gLog.Separator("Singles");
706 gLog << "Extract singles with '" << drsfile << "'" << endl;
707 gLog << endl;
708
709 // ------------------------------------------------------
710
711 MDrsCalibration drscalib300;
712 if (!drscalib300.ReadFits(drsfile))
713 return 3;
714
715 // ------------------------------------------------------
716
717 iter.Sort();
718 iter.Print();
719
720 TString title;
721 title = iter.Next();
722 title += "; ";
723 title += max1;
724 iter.Reset();
725
726 // ======================================================
727
728 MStatusDisplay *d = new MStatusDisplay;
729
730 MBadPixelsCam badpixels;
731 badpixels.InitSize(1440);
732 badpixels[ 424].SetUnsuitable(MBadPixelsPix::kUnsuitable);
733 badpixels[ 583].SetUnsuitable(MBadPixelsPix::kUnsuitable);
734 badpixels[ 830].SetUnsuitable(MBadPixelsPix::kUnsuitable);
735 badpixels[ 923].SetUnsuitable(MBadPixelsPix::kUnsuitable);
736 badpixels[1208].SetUnsuitable(MBadPixelsPix::kUnsuitable);
737 badpixels[1399].SetUnsuitable(MBadPixelsPix::kUnsuitable);
738
739 // Twin pixel
740 // 113
741 // 115
742 // 354
743 // 423
744 // 1195
745 // 1393
746
747 MH::SetPalette();
748
749 MContinue cont("(MRawEvtHeader.GetTriggerID&0xff00)!=0x100", "SelectCal");
750
751 MGeomApply apply;
752
753 MDrsCalibApply drsapply;
754
755 MRawFitsRead read;
756 read.LoadMap(map);
757 read.AddFiles(iter);
758
759 // ======================================================
760
761 gLog << endl;
762 gLog.Separator("Calculating baseline");
763
764 MTaskList tlist0;
765
766 MRawRunHeader header;
767 MPedestalCam pedcam;
768
769 MParList plist;
770 plist.AddToList(&tlist0);
771 plist.AddToList(&drscalib300);
772 plist.AddToList(&header);
773 plist.AddToList(&pedcam);
774
775 // ------------------ Setup the tasks ---------------
776
777 MFillH fill0("MHBaseline", "MPedestalSubtractedEvt");
778
779 drsapply.SetRemoveSpikes(3);
780
781 // ------------------ Setup eventloop and run analysis ---------------
782
783 tlist0.AddToList(&read);
784 tlist0.AddToList(&apply);
785 tlist0.AddToList(&drsapply);
786 tlist0.AddToList(&fill0);
787
788 // ------------------ Setup and run eventloop ---------------
789
790 MEvtLoop loop0(title);
791 loop0.SetDisplay(d);
792 loop0.SetParList(&plist);
793
794 if (!loop0.Eventloop(max0))
795 return 4;
796
797 if (!loop0.GetDisplay())
798 return 5;
799
800 // ----------------------------------------------------------------
801
802 MGeomCamFACT fact;
803 MHCamera hped(fact);
804 hped.SetCamContent(pedcam);
805 hped.SetCamError(pedcam, 4);
806 hped.SetAllUsed();
807 MHCamEvent display;
808 display.SetHist(hped);
809
810 d->AddTab("Baseline");
811 display.Clone()->Draw();
812
813 // ======================================================
814
815 gLog << endl;
816 gLog.Separator("Extracting singles");
817
818 MTaskList tlist1;
819
820 MSingles singles;
821
822 plist.Replace(&tlist1);
823 plist.AddToList(&badpixels);
824 plist.AddToList(&singles);
825
826 // ------------------ Setup the tasks ---------------
827
828 MExtractSingles extract;
829 extract.SetMaxDist(max_dist);
830 extract.SetThreshold(threshold);
831
832 MFillH fill("MHSingles");
833
834 // ------------------ Setup eventloop and run analysis ---------------
835
836 tlist1.AddToList(&read);
837 tlist1.AddToList(&apply);
838 tlist1.AddToList(&drsapply);
839 tlist1.AddToList(&extract);
840 tlist1.AddToList(&fill);
841
842 // ------------------ Setup and run eventloop ---------------
843
844 MEvtLoop loop1(title);
845 loop1.SetDisplay(d);
846 loop1.SetParList(&plist);
847
848 if (!loop1.Eventloop(max1))
849 return 6;
850
851 if (!loop1.GetDisplay())
852 return 7;
853
854 // =============================================================
855
856 gStyle->SetOptFit(kTRUE);
857
858 MHSingles* h = (MHSingles*) plist.FindObject("MHSingles");
859 if (!h)
860 return 8;
861
862 const TH2 *htime = h->GetTime();
863 const TH2 *hpulse = h->GetPulse();
864
865 d->AddTab("Time");
866 TH1 *ht = htime->ProjectionY();
867 ht->SetYTitle("Counts [a.u.]");
868 ht->Scale(1./1440);
869 ht->Draw();
870
871 d->AddTab("Pulse");
872 TH1 *hp = hpulse->ProjectionY();
873 hp->SetYTitle("Counts [a.u.]");
874 hp->Scale(1./1440);
875 hp->Draw();
876
877 d->SaveAs(outfile);
878
879 TFile f(outfile, "UPDATE");
880
881 MParameterI par("NumEvents");
882 par.SetVal(fill.GetNumExecutions());
883 par.Write();
884
885 MParameterI win("IntegrationWindow");
886 win.SetVal(extract.GetIntegrationSize());
887 win.Write();
888
889 extract.GetExtractionRange().Write("ExtractionRange");
890
891 return 0;
892}
Note: See TracBrowser for help on using the repository browser.