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

Last change on this file since 17227 was 17074, checked in by tbretz, 12 years ago
Fixed a crash when no events were read; do not write a file when no events were processed; write the header to the output when only one file was processed.
File size: 24.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 <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 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(1440),
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 = "/fact/raw/2012/07/23/20120723_",
649 int firstRunID = 6,
650 int lastRunID = 6,
651 const char *drsfile = "/fact/raw/2012/07/23/20120723_004.drs.fits.gz",
652 const char *outpath = "."
653 )
654{
655 // ======================================================
656
657 MDirIter iter;
658
659 // built outpu 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 return 0;
667
668 // add input files to directory iterator
669 for (Int_t fileNr = firstRunID; fileNr <= lastRunID; fileNr++)
670 {
671 TString filename = runfile;
672 filename += Form("%03d.fits", fileNr);
673 iter.AddDirectory(gSystem->DirName(runfile), gSystem->BaseName(filename+".fz"));
674 iter.AddDirectory(gSystem->DirName(runfile), gSystem->BaseName(filename+".gz"));
675 }
676
677 // ======================================================
678
679 // true: Display correctly mapped pixels in the camera displays
680 // but the value-vs-index plot is in software/spiral indices
681 // false: Display pixels in hardware/linear indices,
682 // but the order is the camera display is distorted
683 bool usemap = true;
684
685 // map file to use (get that from La Palma!)
686 const char *map = usemap ? "/scratch/fact/FACTmap111030.txt" : NULL;
687
688 // ------------------------------------------------------
689
690 Long_t max0 = 0;
691 Long_t max1 = 0;
692
693 // ======================================================
694
695 if (map && gSystem->AccessPathName(map, kFileExists))
696 {
697 gLog << "ERROR - Cannot access mapping file '" << map << "'" << endl;
698 return 1;
699 }
700 if (gSystem->AccessPathName(drsfile, kFileExists))
701 {
702 gLog << "ERROR - Cannot access drsfile file '" << drsfile << "'" << endl;
703 return 2;
704 }
705
706 // --------------------------------------------------------------------------------
707
708 gLog.Separator("Singles");
709 gLog << "Extract singles with '" << drsfile << "'" << endl;
710 gLog << endl;
711
712 // ------------------------------------------------------
713
714 MDrsCalibration drscalib300;
715 if (!drscalib300.ReadFits(drsfile))
716 return 3;
717
718 // ------------------------------------------------------
719
720 iter.Sort();
721 iter.Print();
722
723 TString title;
724 title = iter.Next();
725 title += "; ";
726 title += max1;
727 iter.Reset();
728
729 // ======================================================
730
731 MStatusDisplay *d = new MStatusDisplay;
732
733 MBadPixelsCam badpixels;
734 badpixels.InitSize(1440);
735 badpixels[ 424].SetUnsuitable(MBadPixelsPix::kUnsuitable);
736 badpixels[ 583].SetUnsuitable(MBadPixelsPix::kUnsuitable);
737 badpixels[ 830].SetUnsuitable(MBadPixelsPix::kUnsuitable);
738 badpixels[ 923].SetUnsuitable(MBadPixelsPix::kUnsuitable);
739 badpixels[1208].SetUnsuitable(MBadPixelsPix::kUnsuitable);
740 badpixels[1399].SetUnsuitable(MBadPixelsPix::kUnsuitable);
741
742 // Twin pixel
743 // 113
744 // 115
745 // 354
746 // 423
747 // 1195
748 // 1393
749
750 MH::SetPalette();
751
752 MContinue cont("(MRawEvtHeader.GetTriggerID&0xff00)!=0x100", "SelectCal");
753
754 MGeomApply apply;
755
756 MDrsCalibApply drsapply;
757
758 MRawFitsRead read;
759 read.LoadMap(map);
760 read.AddFiles(iter);
761
762 // ======================================================
763
764 gLog << endl;
765 gLog.Separator("Calculating baseline");
766
767 MTaskList tlist0;
768
769 MRawRunHeader header;
770 MPedestalCam pedcam;
771
772 MParList plist;
773 plist.AddToList(&tlist0);
774 plist.AddToList(&drscalib300);
775 plist.AddToList(&header);
776 plist.AddToList(&pedcam);
777
778 // ------------------ Setup the tasks ---------------
779
780 MFillH fill0("MHBaseline", "MPedestalSubtractedEvt");
781
782 drsapply.SetRemoveSpikes(3);
783
784 // ------------------ Setup eventloop and run analysis ---------------
785
786 tlist0.AddToList(&read);
787 tlist0.AddToList(&apply);
788 tlist0.AddToList(&drsapply);
789 tlist0.AddToList(&fill0);
790
791 // ------------------ Setup and run eventloop ---------------
792
793 MEvtLoop loop0(title);
794 loop0.SetDisplay(d);
795 loop0.SetParList(&plist);
796
797 if (!loop0.Eventloop(max0))
798 return 4;
799
800 if (!loop0.GetDisplay())
801 return 5;
802
803 // ----------------------------------------------------------------
804
805 MGeomCamFACT fact;
806 MHCamera hped(fact);
807 hped.SetCamContent(pedcam);
808 hped.SetCamError(pedcam, 4);
809 hped.SetAllUsed();
810 MHCamEvent display;
811 display.SetHist(hped);
812
813 d->AddTab("Baseline");
814 display.Clone()->Draw();
815
816 // ======================================================
817
818 gLog << endl;
819 gLog.Separator("Extracting singles");
820
821 MTaskList tlist1;
822
823 MSingles singles;
824
825 plist.Replace(&tlist1);
826 plist.AddToList(&badpixels);
827 plist.AddToList(&singles);
828
829 // ------------------ Setup the tasks ---------------
830
831 MExtractSingles extract;
832 extract.SetMaxDist(max_dist);
833 extract.SetThreshold(threshold);
834
835 MFillH fill("MHSingles");
836
837 // ------------------ Setup eventloop and run analysis ---------------
838
839 tlist1.AddToList(&read);
840 tlist1.AddToList(&apply);
841 tlist1.AddToList(&drsapply);
842 tlist1.AddToList(&extract);
843 tlist1.AddToList(&fill);
844
845 // ------------------ Setup and run eventloop ---------------
846
847 MEvtLoop loop1(title);
848 loop1.SetDisplay(d);
849 loop1.SetParList(&plist);
850
851 if (!loop1.Eventloop(max1))
852 return 6;
853
854 if (!loop1.GetDisplay())
855 return 7;
856
857 if (fill.GetNumExecutions()==0)
858 return 8;
859
860 // =============================================================
861
862 gStyle->SetOptFit(kTRUE);
863
864 MHSingles* h = (MHSingles*) plist.FindObject("MHSingles");
865 if (!h)
866 return 9;
867
868 const TH2 *htime = h->GetTime();
869 const TH2 *hpulse = h->GetPulse();
870
871 d->AddTab("Time");
872 TH1 *ht = htime->ProjectionY();
873 ht->SetYTitle("Counts [a.u.]");
874 ht->Scale(1./1440);
875 ht->Draw();
876
877 d->AddTab("Pulse");
878 TH1 *hp = hpulse->ProjectionY();
879 hp->SetYTitle("Counts [a.u.]");
880 hp->Scale(1./1440);
881 hp->Draw();
882
883 d->SaveAs(outfile);
884
885 TFile f(outfile, "UPDATE");
886
887 MParameterI par("NumEvents");
888 par.SetVal(fill.GetNumExecutions());
889 par.Write();
890
891 MParameterI win("IntegrationWindow");
892 win.SetVal(extract.GetIntegrationSize());
893 win.Write();
894
895 extract.GetExtractionRange().Write("ExtractionRange");
896
897 if (firstRunID==lastRunID)
898 header.Write();
899
900 return 0;
901}
Note: See TracBrowser for help on using the repository browser.