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

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