source: trunk/Mars/hawc/extract_singles.C@ 19856

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