source: trunk/Mars/hawc/extract_pulse.C@ 19948

Last change on this file since 19948 was 19876, checked in by tbretz, 5 years ago
New macro for pulse extraction (taken from calibration paper)
File size: 40.2 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 "MStatusArray.h"
36#include "MStatusDisplay.h"
37#include "MTaskInteractive.h"
38#include "MPedestalSubtractedEvt.h"
39#include "MHCamera.h"
40#include "MGeomCamFAMOUS.h"
41#include "MRawRunHeader.h"
42#include "MPedestalCam.h"
43#include "MPedestalPix.h"
44#include "MParameters.h"
45#include "PixelMap.h"
46
47using namespace std;
48using namespace TMath;
49
50MHCamera *hgain = 0;
51MHCamera *hbase = 0;
52
53PixelMap pmap;
54
55// Structure to store the extracted value and the extracted time
56struct Single
57{
58 float fSignal;
59 float fTime;
60};
61
62//Storage class to kKeep a list of the extracted single
63class MSingles : public MParContainer, public MCamEvent
64{
65 Int_t fIntegrationWindow;
66
67 vector<vector<Single>> fData;
68
69public:
70 MSingles(const char *name=NULL, const char *title=NULL) : fIntegrationWindow(30)
71 {
72 fName = name ? name : "MSingles";
73 fName = title ? title : "Storeage for vector of singles";
74 }
75
76 void InitSize(const UInt_t i)
77 {
78 fData.resize(i);
79 }
80
81 vector<Single> &operator[](UInt_t i) { return fData[i]; }
82 vector<Single> &GetVector(UInt_t i) { return fData[i]; }
83
84 UInt_t GetNumPixels() const { return fData.size(); }
85
86 void SetIntegrationWindow(Int_t w) { fIntegrationWindow = w; }
87 Int_t GetIntegrationWindow() const { return fIntegrationWindow; }
88
89 Bool_t GetPixelContent(Double_t &, Int_t , const MGeomCam &, Int_t) const
90 {
91 return kTRUE;
92 }
93 void DrawPixelContent(Int_t) const { }
94
95 ClassDef(MSingles, 1)
96};
97
98// Histogram class to extract the baseline
99class MHBaseline : public MH
100{
101 TH2F fBaseline;
102
103 MPedestalCam *fPedestal;
104
105 // The baseline is only extracted where also the signal is extracted
106 // FIXME: Make sure this is consistent with MExtractSingles
107 UShort_t fSkipStart;
108 UShort_t fSkipEnd;
109
110public:
111 MHBaseline() : fPedestal(0), fSkipStart(20), fSkipEnd(10)
112 {
113 fName = "MHBaseline";
114
115 // Setup the histogram
116 fBaseline.SetName("Baseline");
117 fBaseline.SetTitle("Median spectrum");
118 fBaseline.SetXTitle("Pixel [idx]");
119 fBaseline.SetYTitle("Median baseline [mV]");
120 fBaseline.SetDirectory(NULL);
121 }
122
123 Bool_t ReInit(MParList *plist)
124 {
125 fPedestal = (MPedestalCam*)plist->FindCreateObj("MPedestalCam");
126 if (!fPedestal)
127 return kFALSE;
128
129 const MRawRunHeader *header = (MRawRunHeader*)plist->FindObject("MRawRunHeader");
130 if (!header)
131 {
132 *fLog << err << "MRawRunHeader not found... abort." << endl;
133 return kFALSE;
134 }
135
136 // Bin width should be around 1 dac count which is about 0.5mV
137 MBinning binsx, binsy;
138 binsx.SetEdges(64, -0.5, 64-0.5);
139 binsy.SetEdges(100, -20.5, 29.5);
140
141 // Setup binnning
142 MH::SetBinning(fBaseline, binsx, binsy);
143
144 return kTRUE;
145 }
146
147 // Fill the samples into the histogram
148 Int_t Fill(const MParContainer *par, const Stat_t)
149 {
150 const MPedestalSubtractedEvt *evt = dynamic_cast<const MPedestalSubtractedEvt*>(par);
151
152 const Int_t n = evt->GetNumSamples()-fSkipStart-fSkipEnd;
153
154 // Loop over all pixels
155 for (int pix=0; pix<64; pix++)
156 {
157 // Get samples for each pixel
158 const Float_t *ptr = evt->GetSamples(pix);
159
160 // Average two consecutive samples
161 for (int i=0; i<n; i+=2)
162 {
163 const Double_t val = 0.5*ptr[i+fSkipStart]+0.5*ptr[i+1+fSkipStart];
164 fBaseline.Fill(pix, val);
165 }
166 }
167
168 return kTRUE;
169 }
170
171 // Extract the baseline value from the distrbutions
172 Bool_t Finalize()
173 {
174 if (!fPedestal)
175 return kTRUE;
176
177 fPedestal->InitSize(fBaseline.GetNbinsX());
178 fPedestal->SetNumEvents(GetNumExecutions());
179
180 Int_t cnt = 0;
181 Double_t avg = 0;
182 Double_t rms = 0;
183
184 // Loop over all 'pixels'
185 for (int x=0; x<64; x++)
186 {
187 // Get the corresponding slice from the histogram
188 TH1D *hist = fBaseline.ProjectionY("proj", x+1, x+1);
189
190 // Get the maximum bin
191 const Int_t bin = hist->GetMaximumBin();
192
193 // Calculate a parabola through this and the surrounding points
194 // on logarithmic values (that's a gaussian)
195
196 //
197 // Quadratic interpolation
198 //
199 // calculate the parameters of a parabula such that
200 // y(i) = a + b*x(i) + c*x(i)^2
201 // y'(i) = b + 2*c*x(i)
202 //
203 //
204
205 // -------------------------------------------------------------------------
206 // a = y2;
207 // b = (y3-y1)/2;
208 // c = (y3+y1)/2 - y2;
209
210 const Double_t v1 = hist->GetBinContent(bin-1);
211 const Double_t v2 = hist->GetBinContent(bin);
212 const Double_t v3 = hist->GetBinContent(bin+1);
213 if (v1<=0 || v2<=0 || v3<=0)
214 continue;
215
216 const Double_t y1 = TMath::Log(v1);
217 const Double_t y2 = TMath::Log(v2);
218 const Double_t y3 = TMath::Log(v3);
219
220 //Double_t a = y2;
221 const Double_t b = (y3-y1)/2;
222 const Double_t c = (y1+y3)/2 - y2;
223 if (c>=0)
224 continue;
225
226 const Double_t w = -1./(2*c);
227 const Double_t dx = b*w;
228
229 if (dx<-1 || dx>1)
230 continue;
231
232 // y = exp( - (x-k)^2 / s^2 / 2 )
233 //
234 // -2*s^2 * log(y) = x^2 - 2*k*x + k^2
235 //
236 // a = (k/s0)^2/2
237 // b = k/s^2
238 // c = -1/(2s^2) --> s = sqrt(-1/2c)
239
240 const Double_t binx = hist->GetBinCenter(bin);
241 const Double_t binw = hist->GetBinWidth(bin);
242
243 //const Double_t p = hist->GetBinCenter(hist->GetMaximumBin());
244 const Double_t p = binx + dx*binw;
245
246 avg += p;
247 rms += p*p;
248
249 cnt++;
250
251 // Store baseline and sigma
252 MPedestalPix &pix = (*fPedestal)[x];
253
254 pix.SetPedestal(p);
255 pix.SetPedestalRms(TMath::Sqrt(w)*binw);
256
257 delete hist;
258 }
259
260 avg /= cnt;
261 rms /= cnt;
262
263 cout << "Baseline(" << cnt << "): " << avg << " +- " << sqrt(rms-avg*avg) << endl;
264
265 return kTRUE;
266 }
267
268 // Draw histogram
269 void Draw(Option_t *)
270 {
271 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
272
273 AppendPad("");
274
275 pad->SetBorderMode(0);
276 pad->SetFrameBorderMode(0);
277 fBaseline.Draw("colz");
278 }
279
280
281 ClassDef(MHBaseline, 1);
282};
283
284// Histogram class for the signal and time distribution as
285// well as the pulse shape
286class MHSingles : public MH
287{
288 TH2F fSignalAll;
289 TH2F fSignalOut;
290 TH2F fSignal1;
291 TH2F fSignal2;
292 TH2F fSignal3;
293 TH2F fSignal4;
294 TH2F fSignal5;
295 TH2F fSignal6;
296 TH2F fSignal7;
297 TH2F fSignal8;
298
299 TProfile fProfAll;
300 TProfile fProf1;
301 TProfile fProf2;
302 TProfile fProf3;
303 TProfile fProf4;
304 TProfile fProf5;
305 TProfile fProf6;
306 TProfile fProf7;
307 TProfile fProf8;
308
309 UInt_t fNumEntries;
310
311 MSingles *fSingles; //!
312 MPedestalSubtractedEvt *fEvent; //!
313 MBadPixelsCam *fBadPix; //!
314 MPedestalCam *fPedestal; //!
315
316public:
317 MHSingles() : fNumEntries(0), fSingles(0), fEvent(0), fPedestal(0)
318 {
319 fName = "MHSingles";
320
321 // Setup histograms
322 fSignalAll.SetName("SignalAll");
323 fSignalAll.SetTitle("");
324 fSignalAll.SetXTitle("Time [ns]");
325 fSignalAll.SetYTitle("Amplitude [a.u.]");
326 fSignalAll.GetXaxis()->CenterTitle();
327 fSignalAll.GetXaxis()->SetTitleSize(0.045);
328 fSignalAll.GetYaxis()->SetTitleSize(0.045);
329 fSignalAll.SetDirectory(NULL);
330 fSignalAll.SetStats(kFALSE);
331
332 fSignalOut.SetName("SignalOut");
333 fSignalOut.SetTitle("");
334 fSignalOut.SetXTitle("Time [ns]");
335 fSignalOut.SetYTitle("Amplitude [a.u.]");
336 fSignalOut.GetXaxis()->CenterTitle();
337 fSignalOut.GetXaxis()->SetTitleSize(0.045);
338 fSignalOut.GetYaxis()->SetTitleSize(0.045);
339 fSignalOut.SetDirectory(NULL);
340 fSignalOut.SetStats(kFALSE);
341
342 fSignal1.SetName("Signal1");
343 fSignal1.SetTitle("");
344 fSignal1.SetXTitle("Time [ns]");
345 fSignal1.SetYTitle("Amplitude [a.u.]");
346 fSignal1.GetXaxis()->CenterTitle();
347 fSignal1.GetXaxis()->SetTitleSize(0.045);
348 fSignal1.GetYaxis()->SetTitleSize(0.045);
349 fSignal1.SetDirectory(NULL);
350 fSignal1.SetStats(kFALSE);
351
352 fSignal2.SetName("Signal2");
353 fSignal2.SetTitle("");
354 fSignal2.SetXTitle("Time [ns]");
355 fSignal2.SetYTitle("Amplitude [a.u.]");
356 fSignal2.GetXaxis()->CenterTitle();
357 fSignal2.GetXaxis()->SetTitleSize(0.045);
358 fSignal2.GetYaxis()->SetTitleSize(0.045);
359 fSignal2.SetDirectory(NULL);
360 fSignal2.SetStats(kFALSE);
361
362 fSignal3.SetName("Signal3");
363 fSignal3.SetTitle("");
364 fSignal3.SetXTitle("Time [ns]");
365 fSignal3.SetYTitle("Amplitude [a.u.]");
366 fSignal3.GetXaxis()->CenterTitle();
367 fSignal3.GetXaxis()->SetTitleSize(0.045);
368 fSignal3.GetYaxis()->SetTitleSize(0.045);
369 fSignal3.SetDirectory(NULL);
370 fSignal3.SetStats(kFALSE);
371
372 fSignal4.SetName("Signal4");
373 fSignal4.SetTitle("");
374 fSignal4.SetXTitle("Time [ns]");
375 fSignal4.SetYTitle("Amplitude [a.u.]");
376 fSignal4.GetXaxis()->CenterTitle();
377 fSignal4.GetXaxis()->SetTitleSize(0.045);
378 fSignal4.GetYaxis()->SetTitleSize(0.045);
379 fSignal4.SetDirectory(NULL);
380 fSignal4.SetStats(kFALSE);
381
382 fSignal5.SetName("Signal5");
383 fSignal5.SetTitle("");
384 fSignal5.SetXTitle("Time [ns]");
385 fSignal5.SetYTitle("Amplitude [a.u.]");
386 fSignal5.GetXaxis()->CenterTitle();
387 fSignal5.GetXaxis()->SetTitleSize(0.045);
388 fSignal5.GetYaxis()->SetTitleSize(0.045);
389 fSignal5.SetDirectory(NULL);
390 fSignal5.SetStats(kFALSE);
391
392 fSignal6.SetName("Signal6");
393 fSignal6.SetTitle("");
394 fSignal6.SetXTitle("Time [ns]");
395 fSignal6.SetYTitle("Amplitude [a.u.]");
396 fSignal6.GetXaxis()->CenterTitle();
397 fSignal6.GetXaxis()->SetTitleSize(0.045);
398 fSignal6.GetYaxis()->SetTitleSize(0.045);
399 fSignal6.SetDirectory(NULL);
400 fSignal6.SetStats(kFALSE);
401
402 fSignal7.SetName("Signal7");
403 fSignal7.SetTitle("");
404 fSignal7.SetXTitle("Time [ns]");
405 fSignal7.SetYTitle("Amplitude [a.u.]");
406 fSignal7.GetXaxis()->CenterTitle();
407 fSignal7.GetXaxis()->SetTitleSize(0.045);
408 fSignal7.GetYaxis()->SetTitleSize(0.045);
409 fSignal7.SetDirectory(NULL);
410 fSignal7.SetStats(kFALSE);
411
412 fSignal8.SetName("Signal8");
413 fSignal8.SetTitle("");
414 fSignal8.SetXTitle("Time [ns]");
415 fSignal8.SetYTitle("Amplitude [a.u.]");
416 fSignal8.GetXaxis()->CenterTitle();
417 fSignal8.GetXaxis()->SetTitleSize(0.045);
418 fSignal8.GetYaxis()->SetTitleSize(0.045);
419 fSignal8.SetDirectory(NULL);
420 fSignal8.SetStats(kFALSE);
421
422 fProfAll.SetName("ProfAll");
423 fProf1.SetName("Prof1");
424 fProf2.SetName("Prof2");
425 fProf3.SetName("Prof3");
426 fProf4.SetName("Prof4");
427 fProf5.SetName("Prof5");
428 fProf6.SetName("Prof6");
429 fProf7.SetName("Prof7");
430 fProf8.SetName("Prof8");
431
432 fProfAll.SetErrorOption("s");
433 fProf1.SetErrorOption("s");
434 fProf2.SetErrorOption("s");
435 fProf3.SetErrorOption("s");
436 fProf4.SetErrorOption("s");
437 fProf5.SetErrorOption("s");
438 fProf6.SetErrorOption("s");
439 fProf7.SetErrorOption("s");
440 fProf8.SetErrorOption("s");
441
442 fProfAll.GetYaxis()->SetTitleOffset(0.8);
443 fProf1.GetYaxis()->SetTitleOffset(0.8);
444 fProf2.GetYaxis()->SetTitleOffset(0.8);
445 fProf3.GetYaxis()->SetTitleOffset(0.8);
446 fProf4.GetYaxis()->SetTitleOffset(0.8);
447 fProf5.GetYaxis()->SetTitleOffset(0.8);
448 fProf6.GetYaxis()->SetTitleOffset(0.8);
449 fProf7.GetYaxis()->SetTitleOffset(0.8);
450 fProf8.GetYaxis()->SetTitleOffset(0.8);
451 }
452
453 Bool_t SetupFill(const MParList *plist)
454 {
455 fSingles = (MSingles*)plist->FindObject("MSingles");
456 if (!fSingles)
457 {
458 *fLog << err << "MSingles not found... abort." << endl;
459 return kFALSE;
460 }
461
462 fEvent = (MPedestalSubtractedEvt*)plist->FindObject("MPedestalSubtractedEvt");
463 if (!fEvent)
464 {
465 *fLog << err << "MPedestalSubtractedEvt not found... abort." << endl;
466 return kFALSE;
467 }
468
469 fBadPix = (MBadPixelsCam*)plist->FindObject("MBadPixelsCam");
470 if (!fBadPix)
471 {
472 *fLog << err << "MBadPixelsCam not found... abort." << endl;
473 return kFALSE;
474 }
475
476 fPedestal = (MPedestalCam*)plist->FindObject("MPedestalCam");
477 if (!fPedestal)
478 {
479 *fLog << err << "MPedestalCam not found... abort." << endl;
480 return kFALSE;
481 }
482
483 fNumEntries = 0;
484
485 return kTRUE;
486 }
487
488 Bool_t ReInit(MParList *plist)
489 {
490 const MRawRunHeader *header = (MRawRunHeader*)plist->FindObject("MRawRunHeader");
491 if (!header)
492 {
493 *fLog << err << "MRawRunHeader not found... abort." << endl;
494 return kFALSE;
495 }
496
497 // Setup binning
498 //const Int_t w = fSingles->GetIntegrationWindow();
499
500 MBinning binsx, binsy;
501 binsx.SetEdges(1024, -512, 512);
502 binsy.SetEdges( 800, -20, 20);
503
504 MH::SetBinning(fSignalAll, binsx, binsy);
505 MH::SetBinning(fSignalOut, binsx, binsy);
506 MH::SetBinning(fSignal1, binsx, binsy);
507 MH::SetBinning(fSignal2, binsx, binsy);
508 MH::SetBinning(fSignal3, binsx, binsy);
509 MH::SetBinning(fSignal4, binsx, binsy);
510 MH::SetBinning(fSignal5, binsx, binsy);
511 MH::SetBinning(fSignal6, binsx, binsy);
512 MH::SetBinning(fSignal7, binsx, binsy);
513 MH::SetBinning(fSignal8, binsx, binsy);
514 MH::SetBinning(fProfAll, binsx);
515 MH::SetBinning(fProf1, binsx);
516 MH::SetBinning(fProf2, binsx);
517 MH::SetBinning(fProf3, binsx);
518 MH::SetBinning(fProf4, binsx);
519 MH::SetBinning(fProf5, binsx);
520 MH::SetBinning(fProf6, binsx);
521 MH::SetBinning(fProf7, binsx);
522 MH::SetBinning(fProf8, binsx);
523
524 fSignal1.GetXaxis()->SetRangeUser(-25, 300);
525 fSignal1.GetYaxis()->SetRangeUser(-1.5, 2.5);
526
527 fSignal2.GetXaxis()->SetRangeUser(-25, 300);
528 fSignal2.GetYaxis()->SetRangeUser(-1.5, 3.5);
529
530 fSignal3.GetXaxis()->SetRangeUser(-25, 300);
531 fSignal3.GetYaxis()->SetRangeUser(-1.5, 4.5);
532
533 fSignal4.GetXaxis()->SetRangeUser(-25, 300);
534 fSignal4.GetYaxis()->SetRangeUser(-1.5, 5.5);
535
536 fSignal5.GetXaxis()->SetRangeUser(-25, 300);
537 fSignal5.GetYaxis()->SetRangeUser(-1.5, 6.5);
538
539 fSignal6.GetXaxis()->SetRangeUser(-25, 300);
540 fSignal6.GetYaxis()->SetRangeUser(-1.5, 7.5);
541
542 fSignal7.GetXaxis()->SetRangeUser(-25, 300);
543 fSignal7.GetYaxis()->SetRangeUser(-1.5, 8.5);
544
545 fSignal8.GetXaxis()->SetRangeUser(-25, 300);
546 fSignal8.GetYaxis()->SetRangeUser(-1.5, 9.5);
547
548 return kTRUE;
549 }
550
551 // Fill singles into histograms
552 Int_t Fill(const MParContainer *, const Stat_t)
553 {
554 // Get pointer to samples to fill samples
555 const Float_t *ptr = fEvent->GetSamples(0);
556 const Short_t roi = fEvent->GetNumSamples();
557
558 if (!ptr)
559 return kTRUE;
560
561 const Int_t nx = fSignalAll.GetNbinsX();
562 const Int_t ny = fSignalAll.GetNbinsY();
563
564 const double x0 = fSignalAll.GetXaxis()->GetBinCenter(1);
565 const double y0 = fSignalAll.GetYaxis()->GetBinCenter(1);
566
567 const double wx = fSignalAll.GetXaxis()->GetBinCenter(nx)-x0;
568 const double wy = fSignalAll.GetYaxis()->GetBinCenter(ny)-y0;
569
570 // Loop over all pixels
571 for (unsigned int i=0; i<64; i++, ptr+=roi)
572 {
573 int idx = i;//pmap.hw(i).index;
574
575 if ((*fBadPix)[i].IsUnsuitable())
576 continue;
577
578 if (!hgain->IsUsed(idx))
579 continue;
580
581 if (idx!=0)
582 continue;
583
584 const double gain = hgain->GetBinContent(idx+1)/30;
585 const double base = hbase->GetBinContent(idx+1);
586
587 const double ped = (*fPedestal)[i].GetPedestal();
588
589 // loop over all singles
590 const vector<Single> &result = fSingles->GetVector(i);
591 if (result.size()!=1)
592 continue;
593
594 // Fill pulse shape
595 const double tm = result[0].fTime;
596 const double sig = result[0].fSignal/30 - base;
597
598 bool ok = true;//tm>3 && (ptr[int(tm)-3]-ped-base)/gain > -5;
599
600 double mean = 0;
601
602 for (int s=20; s<roi-20; s++)
603 {
604 double x = (s - tm)/2; // [ns]
605 double y = (ptr[s]-ped-base)/gain;
606
607 if (x>-100 && x<-10)
608 mean += y/90;
609
610 if (y<-3)
611 ok = false;
612 if (x>5 && x<40 && y<=x/40*-0.5)
613 ok = false;
614 if (x>-150 && x<-5 && y>1.5)
615 ok = false;
616 if (y>2.5)
617 ok = false;
618 }
619
620 // Nothing at or bleow zero from 5 to 40machen
621
622 bool ok1 = sig>gain*0.5 && sig<gain*1.5;
623 bool ok2 = sig>gain*1.5 && sig<gain*2.5;
624 bool ok3 = sig>gain*2.5 && sig<gain*3.5;
625 bool ok4 = sig>gain*3.5 && sig<gain*4.5;
626 bool ok5 = sig>gain*4.5 && sig<gain*5.5;
627 bool ok6 = sig>gain*5.5 && sig<gain*6.5;
628 bool ok7 = sig>gain*6.5 && sig<gain*7.5;
629 bool ok8 = sig>gain*7.5 && sig<gain*8.5;
630
631
632 for (int s=20; s<roi-20; s++)
633 {
634 double x = (s - tm)/2; // [ns]
635 double y = (ptr[s]-ped-base)/gain - mean;
636
637 int ix = TMath::Nint((x-x0)/wx*nx);
638 int iy = TMath::Nint((y-y0)/wy*ny);
639
640 if (ix<0 || iy<0 || ix>=nx || iy>=ny)
641 continue;
642
643 int bin = (iy+1)*(nx+2)+(ix+1);
644
645 if (ok)
646 {
647 fSignalAll.AddBinContent(bin);
648 fProfAll.Fill(x, y);
649
650 if (ok1)
651 {
652 fSignal1.AddBinContent(bin);
653 fProf1.Fill(x, y);
654 }
655 if (ok2)
656 {
657 fSignal2.AddBinContent(bin);
658 fProf2.Fill(x, y);
659 }
660 if (ok3)
661 {
662 fSignal3.AddBinContent(bin);
663 fProf3.Fill(x, y);
664 }
665 if (ok4)
666 {
667 fSignal4.AddBinContent(bin);
668 fProf4.Fill(x, y);
669 }
670 if (ok5)
671 {
672 fSignal5.AddBinContent(bin);
673 fProf5.Fill(x, y);
674 }
675 if (ok6)
676 {
677 fSignal6.AddBinContent(bin);
678 fProf6.Fill(x, y);
679 }
680 if (ok7)
681 {
682 fSignal7.AddBinContent(bin);
683 fProf7.Fill(x, y);
684 }
685 if (ok8)
686 {
687 fSignal8.AddBinContent(bin);
688 fProf8.Fill(x, y);
689 }
690 }
691 else
692 fSignalOut.AddBinContent(bin);
693 }
694
695 if (ok)
696 {
697 fSignalAll.SetEntries(fSignalAll.GetEntries()+1);
698
699 if (ok1)
700 fSignal1.SetEntries(fSignal1.GetEntries()+1);
701 if (ok2)
702 fSignal2.SetEntries(fSignal2.GetEntries()+1);
703 if (ok3)
704 fSignal3.SetEntries(fSignal3.GetEntries()+1);
705 if (ok4)
706 fSignal4.SetEntries(fSignal4.GetEntries()+1);
707 if (ok5)
708 fSignal5.SetEntries(fSignal5.GetEntries()+1);
709 if (ok6)
710 fSignal6.SetEntries(fSignal6.GetEntries()+1);
711 if (ok7)
712 fSignal7.SetEntries(fSignal7.GetEntries()+1);
713 if (ok8)
714 fSignal8.SetEntries(fSignal8.GetEntries()+1);
715 }
716 else
717 fSignalOut.SetEntries(fSignalOut.GetEntries()+1);
718 }
719
720 fNumEntries++;
721
722 return kTRUE;
723 }
724
725 // Getter for histograms
726 TH2 *GetSignalAll() { return &fSignalAll; }
727 TH2 *GetSignalOut() { return &fSignalOut; }
728 TH2 *GetSignal1() { return &fSignal1; }
729 TH2 *GetSignal2() { return &fSignal2; }
730 TH2 *GetSignal3() { return &fSignal3; }
731 TH2 *GetSignal4() { return &fSignal4; }
732 TH2 *GetSignal5() { return &fSignal5; }
733 TH2 *GetSignal6() { return &fSignal6; }
734 TH2 *GetSignal7() { return &fSignal7; }
735 TH2 *GetSignal8() { return &fSignal8; }
736 TH1 *GetProf1() { return &fProf1; }
737 TH1 *GetProf2() { return &fProf2; }
738 TH1 *GetProf3() { return &fProf3; }
739 TH1 *GetProf4() { return &fProf4; }
740 TH1 *GetProf5() { return &fProf5; }
741 TH1 *GetProf6() { return &fProf6; }
742 TH1 *GetProf7() { return &fProf7; }
743 TH1 *GetProf8() { return &fProf8; }
744
745 UInt_t GetNumEntries() const { return fNumEntries; }
746
747 void Draw(Option_t *)
748 {
749 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
750
751 AppendPad("");
752
753 pad->Divide(2,2);
754
755 pad->cd(1);
756 gPad->SetBorderMode(0);
757 gPad->SetFrameBorderMode(0);
758 fSignalAll.Draw("colz");
759 fProfAll.Draw("same");
760
761 pad->cd(2);
762 gPad->SetBorderMode(0);
763 gPad->SetFrameBorderMode(0);
764 fSignalOut.Draw("colz");
765
766 pad->cd(3);
767 gPad->SetBorderMode(0);
768 gPad->SetFrameBorderMode(0);
769 fSignal1.Draw("colz");
770 fProf1.Draw("same");
771
772 pad->cd(4);
773 gPad->SetBorderMode(0);
774 gPad->SetFrameBorderMode(0);
775 fSignal8.Draw("colz");
776 //fProf8.Draw("same");
777 }
778
779 void DrawCopy(Option_t *)
780 {
781 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
782
783 AppendPad("");
784
785 pad->Divide(2,2);
786
787 pad->cd(1);
788 gPad->SetBorderMode(0);
789 gPad->SetFrameBorderMode(0);
790 fSignalAll.DrawCopy("colz");
791
792 pad->cd(2);
793 gPad->SetBorderMode(0);
794 gPad->SetFrameBorderMode(0);
795 fSignalOut.DrawCopy("colz");
796
797 pad->cd(3);
798 gPad->SetBorderMode(0);
799 gPad->SetFrameBorderMode(0);
800 fSignal1.DrawCopy("colz");
801
802 pad->cd(4);
803 gPad->SetBorderMode(0);
804 gPad->SetFrameBorderMode(0);
805 fSignal8.DrawCopy("colz");
806 }
807
808 ClassDef(MHSingles, 2)
809};
810
811// Task to extract the singles
812class MExtractSingles : public MTask
813{
814 MSingles *fSingles;
815 MPedestalCam *fPedestal;
816 MPedestalSubtractedEvt *fEvent;
817
818 // On time for each pixel in samples
819 MArrayI fExtractionRange;
820
821 // Number of samples for sliding average
822 size_t fNumAverage;
823
824 // The range in which the singles have to fit in
825 // FIXME: Make sure this is consistent with MHBaseline
826 UInt_t fStartSkip;
827 UInt_t fEndSkip;
828
829 UInt_t fIntegrationSize;
830 UInt_t fMaxSearchWindow;
831
832 Int_t fMaxDist;
833 Double_t fThreshold;
834
835public:
836 MExtractSingles() : fSingles(0), fPedestal(0), fEvent(0),
837 fExtractionRange(64),
838 fNumAverage(10), fStartSkip(20), fEndSkip(10),
839 fIntegrationSize(3*10), fMaxSearchWindow(30)
840 {
841 }
842
843 void SetMaxDist(Int_t m) { fMaxDist = m; }
844 void SetThreshold(Int_t th) { fThreshold = th; }
845
846 UInt_t GetIntegrationSize() const { return fIntegrationSize; }
847
848 const MArrayI &GetExtractionRange() const { return fExtractionRange; }
849
850
851 Int_t PreProcess(MParList *plist)
852 {
853 fSingles = (MSingles*)plist->FindCreateObj("MSingles");
854 if (!fSingles)
855 return kFALSE;
856
857 fEvent = (MPedestalSubtractedEvt*)plist->FindObject("MPedestalSubtractedEvt");
858 if (!fEvent)
859 {
860 *fLog << err << "MPedestalSubtractedEvt not found... abort." << endl;
861 return kFALSE;
862 }
863
864 fPedestal = (MPedestalCam*)plist->FindObject("MPedestalCam");
865 if (!fPedestal)
866 {
867 *fLog << err << "MPedestalCam not found... abort." << endl;
868 return kFALSE;
869 }
870
871 return kTRUE;
872 }
873
874 Int_t Process()
875 {
876 const UInt_t roi = fEvent->GetNumSamples();
877
878 const size_t navg = fNumAverage;
879 const float threshold = fThreshold;
880 const UInt_t start_skip = fStartSkip;
881 const UInt_t end_skip = fEndSkip;
882 const UInt_t integration_size = fIntegrationSize;
883 const UInt_t max_search_window = fMaxSearchWindow;
884 const UInt_t max_dist = fMaxDist;
885
886 vector<float> val(roi-navg);//result of first sliding average
887 for (int pix=0; pix<64; pix++)
888 {
889 // Total number of samples as 'on time'
890 fExtractionRange[pix] += roi-start_skip-navg-end_skip-integration_size;
891
892 // Clear singles for this pixel
893 vector<Single> &result = fSingles->GetVector(pix);
894 result.clear();
895
896 // Get pointer to samples
897 const Float_t *ptr = fEvent->GetSamples(pix);
898
899 // Get Baseline
900 const Double_t ped = (*fPedestal)[pix].GetPedestal();
901
902 // Subtract baseline and do a sliding average over
903 // the samples to remove noise (mainly the 200MHz noise)
904 for (size_t i=0; i<roi-navg; i++)
905 {
906 val[i] = 0;
907 for (size_t j=i; j<i+navg; j++)
908 val[i] += ptr[j];
909 val[i] /= navg;
910 val[i] -= ped;
911 }
912
913 // peak finding via threshold
914 UInt_t imax = roi-navg-end_skip-integration_size;
915 for (UInt_t i=start_skip; i<imax; i++)
916 {
917 //search for threshold crossings
918 if (val[i+0]>threshold ||
919 val[i+4]>threshold ||
920
921 val[i+5]<threshold ||
922 val[i+9]<threshold)
923 continue;
924
925 //search for maximum after threshold crossing
926 UInt_t k_max = i+5;
927 for (UInt_t k=i+6; k<i+max_search_window; k++)
928 {
929 if (val[k] > val[k_max])
930 k_max = k;
931 }
932
933 // Check if the findings make sense
934 if (k_max == i+5 || k_max == i + max_search_window-1)
935 continue;
936
937 //search for half maximum before maximum
938 UInt_t k_half_max = 0;
939 for (UInt_t k=k_max; k>k_max-25; k--)
940 {
941 if (val[k-1] < val[k_max]/2 &&
942 val[k] >= val[k_max]/2 )
943 {
944 k_half_max = k;
945 break;
946 }
947 }
948 // Check if the finding makes sense
949 if (k_half_max+integration_size > roi-navg-end_skip)
950 continue;
951 if (k_half_max == 0)
952 continue;
953 if (k_max - k_half_max > max_dist)
954 continue;
955
956 // FIXME: Evaluate arrival time more precisely!!!
957 // FIXME: Get a better integration window
958
959 // Compile "single"
960 Single single;
961 single.fSignal = 0;
962 single.fTime = k_half_max;
963
964 // Crossing of the threshold is at 5
965 for (UInt_t j=k_half_max; j<k_half_max+integration_size; j++)
966 single.fSignal += ptr[j];
967
968 single.fSignal -= integration_size*ped;
969
970 // Add single to list
971 result.push_back(single);
972
973 // skip falling edge
974 i += 5+integration_size;
975
976 // Remove skipped samples from 'on time'
977 fExtractionRange[pix] -= i>imax ? 5+integration_size-(i-imax) : 5+integration_size;
978 }
979 }
980 return kTRUE;
981 }
982};
983
984int extract_pulse(
985 Int_t max_dist = 14,
986 Double_t threshold = 5,
987 const char *runfile = "/media/tbretz/d84f26d3-fc9b-4b30-8cfd-85246846dd1a/hawcseye01/raw/2019/10/27/20191027_",
988 int firstRunID = 6,
989 int lastRunID = 6,
990 const char *drsfile = "/media/tbretz/d84f26d3-fc9b-4b30-8cfd-85246846dd1a/hawcseye01/raw/2019/10/27/20191027_003.drs.fits",
991 const char *outpath = "."
992 )
993{
994 // ======================================================
995 if (!pmap.Read("../hawc/FAMOUSmap171218.txt"))
996 {
997 gLog << err << "../hawc/FAMOUSmap171218.txt not found." << endl;
998 return 1;
999 }
1000
1001 // ======================================================
1002
1003
1004 MDirIter iter;
1005
1006 // built output filename and path
1007 TString outfile = "20191027_006_006-pulse.root";
1008 TString infile = "20190227_006_006-fit.root";
1009
1010 //if (!gSystem->AccessPathName(outfile))
1011 // return 0;
1012
1013 // add input files to directory iterator
1014 for (Int_t fileNr = firstRunID; fileNr <= lastRunID; fileNr++)
1015 {
1016 TString filename = runfile;
1017 filename += Form("%03d.fits", fileNr);
1018 iter.AddDirectory(gSystem->DirName(runfile), gSystem->BaseName(filename+".fz"));
1019 iter.AddDirectory(gSystem->DirName(runfile), gSystem->BaseName(filename+".gz"));
1020 iter.AddDirectory(gSystem->DirName(runfile), gSystem->BaseName(filename));
1021 }
1022
1023 // ======================================================
1024
1025 TFile file(infile);
1026 MStatusArray arr;
1027 if (arr.Read()<=0)
1028 {
1029 gLog << "ERROR - Cannot open gain file '" << infile << "'" << endl;
1030 return 1;
1031 }
1032
1033 hgain = (MHCamera*)arr.FindObjectInCanvas("Gain_copy", "MHCamera", "Cams1");
1034 hbase = (MHCamera*)arr.FindObjectInCanvas("Baseline_copy", "MHCamera", "Cams1");
1035 if (!hgain || !hbase)
1036 {
1037 gLog << "ERROR - Cannot find Gain/Baseline" << endl;
1038 return 1;
1039 }
1040
1041 // ======================================================
1042
1043 // true: Display correctly mapped pixels in the camera displays
1044 // but the value-vs-index plot is in software/spiral indices
1045 // false: Display pixels in hardware/linear indices,
1046 // but the order is the camera display is distorted
1047 bool usemap = true;
1048
1049 // map file to use (get that from La Palma!)
1050 const char *map = usemap ? "../hawc/FAMOUSmap171218.txt" : NULL;
1051
1052 // ------------------------------------------------------
1053
1054 Long_t max0 = 0;
1055 Long_t max1 = 0;
1056
1057 // ======================================================
1058
1059 if (map && gSystem->AccessPathName(map, kFileExists))
1060 {
1061 gLog << "ERROR - Cannot access mapping file '" << map << "'" << endl;
1062 return 1;
1063 }
1064 if (gSystem->AccessPathName(drsfile, kFileExists))
1065 {
1066 gLog << "ERROR - Cannot access drsfile file '" << drsfile << "'" << endl;
1067 return 2;
1068 }
1069
1070 // --------------------------------------------------------------------------------
1071
1072 gLog.Separator("Pulses");
1073 gLog << all << "Extract pulses with '" << drsfile << "'" << endl;
1074 gLog << endl;
1075
1076 // ------------------------------------------------------
1077
1078 MDrsCalibration drscalib300;
1079 if (!drscalib300.ReadFits(drsfile))
1080 return 3;
1081
1082 // ------------------------------------------------------
1083
1084 iter.Sort();
1085 iter.Print();
1086
1087 TString title;
1088 title = iter.Next();
1089 title += "; ";
1090 title += max1;
1091 iter.Reset();
1092
1093 // ======================================================
1094
1095 MStatusDisplay *d = new MStatusDisplay;
1096
1097 MBadPixelsCam badpixels;
1098 badpixels.InitSize(64);
1099 badpixels[7].SetUnsuitable(MBadPixelsPix::kUnsuitableRun);
1100 badpixels[17].SetUnsuitable(MBadPixelsPix::kUnsuitableRun);
1101
1102 MH::SetPalette();
1103
1104 MContinue cont("(MRawEvtHeader.GetTriggerID&0xff00)!=0x100", "SelectCal");
1105
1106 MGeomApply apply;
1107
1108 MDrsCalibApply drsapply;
1109 drsapply.SetMaxNumPrevEvents(0); // Switched off step-correction -> crashes due to the requirement of N*9 ch
1110
1111 MRawFitsRead read;
1112 read.LoadMap(map);
1113 read.AddFiles(iter);
1114
1115 // ======================================================
1116
1117 gLog << endl;
1118 gLog.Separator("Calculating baseline");
1119
1120 MTaskList tlist0;
1121
1122 MRawRunHeader header;
1123 MPedestalCam pedcam;
1124
1125 MParList plist;
1126 plist.AddToList(&tlist0);
1127 plist.AddToList(&drscalib300);
1128 plist.AddToList(&header);
1129 plist.AddToList(&pedcam);
1130
1131 // ------------------ Setup the tasks ---------------
1132
1133 MFillH fill0("MHBaseline", "MPedestalSubtractedEvt");
1134
1135 drsapply.SetRemoveSpikes(4);
1136
1137 // ------------------ Setup eventloop and run analysis ---------------
1138
1139 tlist0.AddToList(&read);
1140 tlist0.AddToList(&apply);
1141 tlist0.AddToList(&drsapply);
1142 tlist0.AddToList(&fill0);
1143
1144 // ------------------ Setup and run eventloop ---------------
1145
1146 MEvtLoop loop0(title);
1147 loop0.SetDisplay(d);
1148 loop0.SetParList(&plist);
1149
1150 if (!loop0.Eventloop(max0))
1151 return 4;
1152
1153 if (!loop0.GetDisplay())
1154 return 5;
1155
1156 // ----------------------------------------------------------------
1157
1158 MGeomCamFAMOUS fact;
1159 MHCamera hped(fact);
1160 hped.SetCamContent(pedcam);
1161 hped.SetCamError(pedcam, 4);
1162 hped.SetAllUsed();
1163 MHCamEvent display;
1164 display.SetHist(hped);
1165
1166 d->AddTab("Baseline");
1167 display.Clone()->Draw();
1168
1169 // ======================================================
1170
1171 gLog << endl;
1172 gLog.Separator("Extracting singles");
1173
1174 MTaskList tlist1;
1175
1176 MSingles singles;
1177
1178 plist.Replace(&tlist1);
1179 plist.AddToList(&badpixels);
1180 plist.AddToList(&singles);
1181
1182 // ------------------ Setup the tasks ---------------
1183
1184 MExtractSingles extract;
1185 extract.SetMaxDist(max_dist);
1186 extract.SetThreshold(threshold);
1187
1188 MFillH fill("MHSingles");
1189
1190 // ------------------ Setup eventloop and run analysis ---------------
1191
1192 tlist1.AddToList(&read);
1193 tlist1.AddToList(&apply);
1194 tlist1.AddToList(&drsapply);
1195 tlist1.AddToList(&extract);
1196 tlist1.AddToList(&fill);
1197
1198 // ------------------ Setup and run eventloop ---------------
1199
1200 MEvtLoop loop1(title);
1201 loop1.SetDisplay(d);
1202 loop1.SetParList(&plist);
1203
1204 if (!loop1.Eventloop(max1))
1205 return 6;
1206
1207 if (!loop1.GetDisplay())
1208 return 7;
1209
1210 if (fill.GetNumExecutions()==0)
1211 return 8;
1212
1213 // =============================================================
1214
1215 MHSingles* h = (MHSingles*) plist.FindObject("MHSingles");
1216 if (!h)
1217 return 9;
1218
1219 gStyle->SetOptStat(10);
1220
1221 TH1 *h1 = h->GetSignal1();
1222 TH1 *p1 = h->GetProf1();
1223
1224 TH1 *h2 = h->GetSignal2();
1225 TH1 *p2 = h->GetProf2();
1226
1227 TH1 *h3 = h->GetSignal3();
1228 TH1 *p3 = h->GetProf3();
1229
1230 TH1 *h4 = h->GetSignal4();
1231 TH1 *p4 = h->GetProf4();
1232
1233 TH1 *h5 = h->GetSignal5();
1234 TH1 *p5 = h->GetProf5();
1235
1236 TH1 *h6 = h->GetSignal6();
1237 TH1 *p6 = h->GetProf6();
1238
1239 TH1 *h7 = h->GetSignal7();
1240 TH1 *p7 = h->GetProf7();
1241
1242 TH1 *h8 = h->GetSignal8();
1243 TH1 *p8 = h->GetProf8();
1244
1245 TF1 pulse("pulse", "[0]*(1-exp(-x*x/[1]))*exp(-x/[2])", 0, 300);
1246 pulse.SetLineColor(kBlack);
1247 pulse.SetParameters(0, 1./0.094, 1./0.05);
1248
1249 // =====
1250
1251 h1->SetMinimum(h1->GetEntries()/1000);
1252 h2->SetMinimum(h2->GetEntries()/1000);
1253 h3->SetMinimum(h3->GetEntries()/1000);
1254 h4->SetMinimum(h4->GetEntries()/1000);
1255 h5->SetMinimum(h5->GetEntries()/1000);
1256 h6->SetMinimum(h6->GetEntries()/1000);
1257 h7->SetMinimum(h7->GetEntries()/1000);
1258 h8->SetMinimum(h8->GetEntries()/1000);
1259
1260 // =====
1261
1262 d->AddTab("1");
1263 gPad->SetGrid();
1264 gPad->SetLeftMargin(0.07);
1265 gPad->SetTopMargin(0.01);
1266 gPad->SetRightMargin(0.15);
1267 h1->DrawCopy("colz");
1268 p1->DrawCopy("same");
1269
1270 pulse.SetParameter(0, 1*1.3);
1271 p1->Fit(&pulse, "N0", "", 1.5, 25);
1272 TF1 *f1 = pulse.DrawCopy("same");
1273
1274 // =====
1275
1276 d->AddTab("2");
1277 gPad->SetGrid();
1278 gPad->SetLeftMargin(0.07);
1279 gPad->SetTopMargin(0.01);
1280 gPad->SetRightMargin(0.15);
1281 h2->DrawCopy("colz");
1282 p2->DrawCopy("same");
1283
1284 pulse.SetParameter(0, 2*1.3);
1285 p2->Fit(&pulse, "N0", "", 1.5, 25);
1286 TF1 *f2 = pulse.DrawCopy("same");
1287
1288 // =====
1289
1290 d->AddTab("3");
1291 gPad->SetGrid();
1292 gPad->SetLeftMargin(0.07);
1293 gPad->SetTopMargin(0.01);
1294 gPad->SetRightMargin(0.15);
1295 h3->DrawCopy("colz");
1296 p3->DrawCopy("same");
1297
1298 pulse.SetParameter(0, 3*1.3);
1299 p3->Fit(&pulse, "N0", "", 1.5, 25);
1300 TF1 *f3 = pulse.DrawCopy("same");
1301
1302 // =====
1303
1304 d->AddTab("4");
1305 gPad->SetGrid();
1306 gPad->SetLeftMargin(0.07);
1307 gPad->SetTopMargin(0.01);
1308 gPad->SetRightMargin(0.15);
1309 h4->DrawCopy("colz");
1310 p4->DrawCopy("same");
1311
1312 pulse.SetParameter(0, 4*1.3);
1313 p4->Fit(&pulse, "N0", "", 1.5, 25);
1314 TF1 *f4 = pulse.DrawCopy("same");
1315
1316 // =====
1317
1318 d->AddTab("5");
1319 gPad->SetGrid();
1320 gPad->SetLeftMargin(0.07);
1321 gPad->SetTopMargin(0.01);
1322 gPad->SetRightMargin(0.15);
1323 h5->DrawCopy("colz");
1324 p5->DrawCopy("same");
1325
1326 pulse.SetParameter(0, 5*1.3);
1327 p5->Fit(&pulse, "N0", "", 1.5, 25);
1328 TF1 *f5 = pulse.DrawCopy("same");
1329
1330 // =====
1331
1332 d->AddTab("6");
1333 gPad->SetGrid();
1334 gPad->SetLeftMargin(0.07);
1335 gPad->SetTopMargin(0.01);
1336 gPad->SetRightMargin(0.15);
1337 h6->DrawCopy("colz");
1338 p6->DrawCopy("same");
1339
1340 pulse.SetParameter(0, 6*1.3);
1341 p6->Fit(&pulse, "N0", "", 1.5, 25);
1342 TF1 *f6 = pulse.DrawCopy("same");
1343
1344 // =====
1345
1346 d->AddTab("7");
1347 gPad->SetGrid();
1348 gPad->SetLeftMargin(0.07);
1349 gPad->SetTopMargin(0.01);
1350 gPad->SetRightMargin(0.15);
1351 h7->DrawCopy("colz");
1352 p7->DrawCopy("same");
1353
1354 pulse.SetParameter(0, 7*1.3);
1355 p7->Fit(&pulse, "N0", "", 1.5, 25);
1356 TF1 *f7 = pulse.DrawCopy("same");
1357
1358 // =====
1359
1360 d->AddTab("8");
1361 gPad->SetGrid();
1362 gPad->SetLeftMargin(0.07);
1363 gPad->SetTopMargin(0.01);
1364 gPad->SetRightMargin(0.15);
1365 h8->DrawCopy("colz");
1366 p8->DrawCopy("same");
1367
1368 pulse.SetParameter(0, 6*1.3);
1369 p8->Fit(&pulse, "N0", "", 1.5, 25);
1370 TF1 *f8 = pulse.DrawCopy("same");
1371
1372 // =====
1373
1374 d->AddTab("Pulse");
1375 gPad->SetGrid();
1376 gPad->SetLeftMargin(0.07);
1377 gPad->SetTopMargin(0.01);
1378 gPad->SetRightMargin(0.15);
1379 p8->SetStats(kFALSE);
1380 p8->GetXaxis()->SetRangeUser(-25, 300);
1381 p8->GetXaxis()->CenterTitle();
1382 p8->GetXaxis()->SetTitleSize(0.045);
1383 p8->GetYaxis()->SetTitleSize(0.045);
1384 p8->GetYaxis()->SetTitleOffset(0.8);
1385 p8->SetYTitle("Amplitude");
1386 p8->SetXTitle("Time [ns]");
1387 p8->DrawCopy();
1388 p7->DrawCopy("same");
1389 p6->DrawCopy("same");
1390 p5->DrawCopy("same");
1391 p4->DrawCopy("same");
1392 p3->DrawCopy("same");
1393 p2->DrawCopy("same");
1394 p1->DrawCopy("same");
1395
1396 f1->Draw("same");
1397 f2->Draw("same");
1398 f3->Draw("same");
1399 f4->Draw("same");
1400 f5->Draw("same");
1401 f6->Draw("same");
1402 f7->Draw("same");
1403 f8->Draw("same");
1404
1405 d->SaveAs("paper_pulse.root");
1406
1407
1408 /*
1409 d->SaveAs(outfile);
1410
1411 TFile f(outfile, "UPDATE");
1412
1413 MParameterI par("NumEvents");
1414 par.SetVal(fill.GetNumExecutions());
1415 par.Write();
1416
1417 MParameterI win("IntegrationWindow");
1418 win.SetVal(extract.GetIntegrationSize());
1419 win.Write();
1420
1421 extract.GetExtractionRange().Write("ExtractionRange");
1422
1423 if (firstRunID==lastRunID)
1424 header.Write();
1425 */
1426 return 0;
1427}
Note: See TracBrowser for help on using the repository browser.