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

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