source: trunk/Mars/fact/analysis/callisto.C@ 17873

Last change on this file since 17873 was 17873, checked in by tbretz, 12 years ago
Turned off usage of light-pulser calculated delays as default.
File size: 30.7 KB
Line 
1#include "MLogManip.h"
2
3int callisto(const char *seqfile="seq/2012/01/23/20120123_023.seq", const char *outpath = "output")
4{
5 // ======================================================
6
7 // true: Display correctly mapped pixels in the camera displays
8 // but the value-vs-index plot is in software/spiral indices
9 // false: Display pixels in hardware/linear indices,
10 // but the order is the camera display is distorted
11 bool usemap = true;
12
13 // map file to use (get that from La Palma!)
14 const char *map = usemap ? "/home/fact/FACT++/FACTmap111030.txt" : NULL;
15
16 Bool_t maximum = kTRUE;
17
18 const char *lp_template = maximum ?
19 "template-lp-extractor-maximum.root" :
20 "template-lp-extractor-leading-edge.root";
21
22 const char *pulse_template = "template-pulse.root";
23
24 // ------------------------------------------------------
25
26 bool use_delays=false;
27
28 // ------------------------------------------------------
29
30 // Extraction range in slices. It will always(!) contain the full range
31 // of integration
32 const int first_slice = 20; // 10ns
33 const int last_slice = 250; // 125ns
34
35 // Note that rise and fall time mean different things whether you use IntegralFixed or IntegraRel:
36 //
37 // IntegralFixed:
38 // * fRiseTime: Number of slices left from arrival time
39 // * fFallTime: Number of slices right from arrival time
40 // IntegralRel:
41 // * fRiseTime: Number of slices left from maximum time
42 // * fFallTime: Number of slices right from maximum time
43 //
44 const int rise_time_cal = maximum ? 40 : 10; // was 13; 5ns
45 const int fall_time_cal = maximum ? 120 : 160; // was 23; 80ns
46
47 const int rise_time_dat = maximum ? 10 : 2; // was 13; was 10; 1ns
48 const int fall_time_dat = maximum ? 40 : 48; // was 23; was 40; 24ns
49
50 // Extraction type: Extract integral and half leading edge
51
52 const int type = maximum ? (MExtralgoSpline::kIntegralRel) : (MExtralgoSpline::kIntegralFixed);
53 //const int type = MExtralgoSpline::kIntegralFixed;
54
55
56 const double heighttm = 0.5; // IntegralAbs { 1.5pe * 9.6mV/pe } / IntegralRel { 0.5 }
57
58 Long_t max = 0; // All
59 Long_t max0 = max; // Time marker
60 Long_t max1 = max; // Light pulser
61 //Long_t max2 = 3000; // Calibration ratio
62 Long_t max3 = max; // Pedestal Rndm
63 Long_t max4 = max; // Pedestal Ext
64 Long_t max5 = max; // Data
65
66 // ======================================================
67
68 if (map && gSystem->AccessPathName(map, kFileExists))
69 {
70 gLog << err << "ERROR - Cannot access mapping file '" << map << "'" << endl;
71 return 1;
72 }
73
74 // The sequence file which defines the files for the analysis
75 MSequence seq(seqfile);
76 if (!seq.IsValid())
77 {
78 gLog << err << "ERROR - Sequence '" << seqfile << "' invalid!" << endl;
79 return 2;
80 }
81
82 // --------------------------------------------------------------------------------
83
84 gLog.Separator("Callisto");
85 gLog << all << "Calibrate data of sequence '" << seq.GetFileName() << "'" << endl;
86 gLog << endl;
87
88 // ------------------------------------------------------
89
90 ostringstream drsname;
91 drsname << gSystem->DirName(seqfile) << "/";
92 drsname << seq.GetNight().GetNightAsInt() << "_";
93 drsname << Form("%03d", seq.GetDrsSequence()) << ".drs.seq";
94
95 MSequence drs(drsname.str().c_str());
96 if (!drs.IsValid())
97 {
98 gLog << err << "ERROR - DRS sequence invalid!" << endl;
99 return 3;
100 }
101
102 gLog << all << "DRS sequence file: " << drsname.str() << '\n' << endl;
103
104 TString drsfile = seq.GetFileName(0, MSequence::kRawDrs);
105 if (drsfile.IsNull())
106 {
107 gLog << err << "No DRS file available in sequence." << endl;
108 return 4;
109 }
110
111 TString timfile = drs.GetFileName(0, MSequence::kFitsDat);
112 TString drs1024 = drs.GetFileName(0, MSequence::kFitsDrs);
113 TString pedfile = seq.GetFileName(0, MSequence::kFitsPed);
114 TString calfile = seq.GetFileName(0, MSequence::kFitsCal);
115
116 gLog << all;
117 gLog << "DRS calib 300: " << drsfile << '\n';
118 gLog << "DRS calib 1024: " << drs1024 << "\n\n";
119
120 MDrsCalibration drscalib300;
121 if (!drscalib300.ReadFits(drsfile.Data()))
122 return 5;
123
124 MDrsCalibration drscalib1024;
125 if (!drscalib1024.ReadFits(drs1024.Data()))
126 return 6;
127
128 gLog << all;
129 gLog << "Time calibration : " << timfile << '\n';
130 gLog << "Pedestal file: " << pedfile << '\n';
131 gLog << "Light Pulser file: " << calfile << '\n' << endl;
132
133 // ------------------------------------------------------
134
135 MDirIter iter;
136 if (seq.GetRuns(iter, MSequence::kFitsDat)<=0)
137 {
138 gLog << err << "ERROR - Sequence valid but without files." << endl;
139 return 7;
140 }
141 iter.Print();
142
143 // ======================================================
144
145 MStatusArray arrt, arrp;
146
147 TFile ft(lp_template);
148 if (arrt.Read()<=0)
149 {
150 gLog << err << "ERROR - Reading LP template from " << lp_template << endl;
151 return 100;
152 }
153
154 MHCamera *lpref = (MHCamera*)arrt.FindObjectInCanvas("ExtCalSig;avg", "MHCamera", "Cam");
155 if (!lpref)
156 {
157 gLog << err << "ERROR - LP Template not found in " << lp_template << endl;
158 return 101;
159 }
160 lpref->SetDirectory(0);
161
162 MHCamera *gain = (MHCamera*)arrt.FindObjectInCanvas("gain", "MHCamera", "Gain");
163 if (!gain)
164 {
165 gLog << err << "ERROR - Gain not found in " << lp_template << endl;
166 return 101;
167 }
168 gain->SetDirectory(0);
169
170 TFile fp(pulse_template);
171 if (arrp.Read()<=0)
172 {
173 gLog << err << "ERROR - Reading Pulse template from " << pulse_template << endl;
174 return 102;
175 }
176
177 TH1F *hpulse = (TH1F*)arrp.FindObjectInCanvas("hPixelEdgeMean0_0", "TH1F", "cgpPixelPulses0");
178 if (!hpulse)
179 {
180 gLog << err << "ERROR - Pulse Template not found in " << pulse_template << endl;
181 return 103;
182 }
183 hpulse->SetDirectory(0);
184
185 // ======================================================
186
187 MStatusDisplay *d = new MStatusDisplay;
188
189 MBadPixelsCam badpixels;
190 badpixels.InitSize(1440);
191 badpixels[ 424].SetUnsuitable(MBadPixelsPix::kUnsuitable);
192 badpixels[ 583].SetUnsuitable(MBadPixelsPix::kUnsuitable);
193 badpixels[ 830].SetUnsuitable(MBadPixelsPix::kUnsuitable);
194 badpixels[ 923].SetUnsuitable(MBadPixelsPix::kUnsuitable);
195 badpixels[1208].SetUnsuitable(MBadPixelsPix::kUnsuitable);
196 badpixels[1399].SetUnsuitable(MBadPixelsPix::kUnsuitable);
197
198 // Twin pixel
199 // 113
200 // 115
201 // 354
202 // 423
203 // 1195
204 // 1393
205
206 MDrsCalibrationTime timecam;
207
208 // Plot the trigger pattern rates vs. run-number
209 MH3 hrate("MRawRunHeader.GetFileID", "MRawEvtHeader.GetTriggerID&0xff00");
210 hrate.SetWeight("1./TMath::Max(MRawRunHeader.GetRunLength,1)");
211 hrate.SetName("Rate");
212 hrate.SetTitle("Event rate [Hz];File Id;Trigger Type;");
213 hrate.InitLabels(MH3::kLabelsXY);
214 hrate.DefineLabelY( 0, "Data"); // What if TriggerID==0 already???
215 hrate.DefineLabelY(0x100, "Cal");
216 hrate.DefineLabelY(0x400, "Ped");
217 // hrate.DefaultLabelY("ERROR");
218
219 Bool_t isinteg =
220 (type&MExtralgoSpline::kIntegral) ||
221 (type&MExtralgoSpline::kFixedWidth) ||
222 (type&MExtralgoSpline::kDynWidth)
223 ? kTRUE : kFALSE;
224
225 gStyle->SetOptFit(kTRUE);
226
227 // ======================================================
228
229 gLog << endl;
230 gLog.Separator("Processing DRS timing calibration run");
231
232 MTaskList tlist0;
233
234 MParList plist0;
235 plist0.AddToList(&tlist0);
236 plist0.AddToList(&drscalib1024);
237 plist0.AddToList(&badpixels);
238 plist0.AddToList(&timecam);
239
240 MEvtLoop loop0("DetermineTimeCal");
241 loop0.SetDisplay(d);
242 loop0.SetParList(&plist0);
243
244 // ------------------ Setup the tasks ---------------
245
246 MRawFitsRead read0(timfile);
247
248 MContinue cont0("MRawEvtHeader.GetTriggerID!=33792", "SelectTim");
249
250 MGeomApply apply0;
251
252 MDrsCalibApply drsapply0;
253 drsapply0.SetRemoveSpikes(4);
254
255 MFillH fill0("MHDrsCalibrationTime");
256 fill0.SetNameTab("DeltaT");
257
258 tlist0.AddToList(&read0);
259 tlist0.AddToList(&apply0);
260 tlist0.AddToList(&drsapply0);
261 tlist0.AddToList(&cont0);
262 tlist0.AddToList(&fill0);
263
264 if (!loop0.Eventloop(max0))
265 return 8;
266
267 if (!loop0.GetDisplay())
268 return 9;
269
270 /*
271 MHDrsCalibrationT *t = (MHDrsCalibrationT*)plist4.FindObject("MHDrsCalibrationT");
272 t->SetDisplay(d);
273 t->PlotAll();
274 */
275
276 // ======================================================
277
278 gLog << endl;
279 gLog.Separator("Processing external light pulser run");
280
281 MTaskList tlist1;
282
283 MParList plist1;
284 plist1.AddToList(&tlist1);
285 plist1.AddToList(&drscalib300);
286 plist1.AddToList(&badpixels);
287 plist1.AddToList(&timecam);
288
289 MEvtLoop loop1("DetermineCalConst");
290 loop1.SetDisplay(d);
291 loop1.SetParList(&plist1);
292
293 // ------------------ Setup the tasks ---------------
294
295 MRawFitsRead read1;
296 read1.LoadMap(map);
297 read1.AddFile(calfile);
298
299 MContinue cont1("(MRawEvtHeader.GetTriggerID&0xff00)!=0x100", "SelectCal");
300
301 MGeomApply apply1;
302
303 MDrsCalibApply drsapply1;
304 drsapply1.SetRemoveSpikes(4);
305
306 /*
307 MPedestalCam fPedestalCamOut1a;
308 MPedestalCam fPedestalCamOut1b;
309 MPedestalCam fPedestalCamOut1c;
310
311 MPedCalcPedRun pedcalc1a;
312 MPedCalcPedRun pedcalc1b;
313 MPedCalcPedRun pedcalc1c;
314 pedcalc1a.SetPedestalsOut(&fPedestalCamOut1a);
315 pedcalc1b.SetPedestalsOut(&fPedestalCamOut1b);
316 pedcalc1c.SetPedestalsOut(&fPedestalCamOut1c);
317
318 MExtractTimeAndChargeSpline extractor1ab;
319 extractor1a.SetRange(first_slice, last_slice);
320 extractor1a.SetRiseTimeHiGain(rise_time);
321 extractor1a.SetFallTimeHiGain(fall_time);
322 extractor1a.SetChargeType(type);
323 extractor1a.SetSaturationLimit(600000);
324 extractor1a.SetNoiseCalculation(kFALSE);
325
326 pedcalc1a.SetRandomCalculation(kTRUE);
327 pedcalc1b.SetRandomCalculation(kFALSE);
328 pedcalc1a.SetExtractor(&extractor1a);
329 pedcalc1b.SetExtractor($extractor1a);
330 pedcalc1c.SetRangeFromExtractor(&extractor1a);
331 */
332
333 // ---
334
335 MExtractTimeAndChargeSpline extractor1b("ExtractPulse");
336 extractor1b.SetRange(first_slice, last_slice);
337 extractor1b.SetRiseTimeHiGain(rise_time_cal);
338 extractor1b.SetFallTimeHiGain(fall_time_cal);
339 extractor1b.SetHeightTm(heighttm);
340 extractor1b.SetChargeType(type);
341 extractor1b.SetSaturationLimit(600000);
342 extractor1b.SetNoiseCalculation(kFALSE);
343
344 MExtractTimeAndChargeSpline extractor1c("ExtractAmplitude");
345 extractor1c.SetRange(first_slice, last_slice);
346 extractor1c.SetChargeType(MExtralgoSpline::kAmplitude);
347 extractor1c.SetSaturationLimit(600000);
348 extractor1c.SetNoiseCalculation(kFALSE);
349 extractor1c.SetNameSignalCam("Amplitude");
350 extractor1c.SetNameTimeCam("AmplitudePos");
351
352 // ---
353
354 MHCamEvent evt1a(5, "CalRatio", "Ratio per slice between integrated signal and amplitude;; r [1/n]");
355 evt1a.SetNameSub("Amplitude", kTRUE);
356 MFillH fill1a(&evt1a, "MExtractedSignalCam", "FillRatio");
357 fill1a.SetDrawOption("gaus");
358
359 MParameterD ratio1a;
360 ratio1a.SetVal(1./(fall_time_cal+rise_time_cal));
361 fill1a.SetWeight(&ratio1a);
362
363 // ---
364
365 MHCamEvent evt1f(0, "ExtCalSig", "Extracted calibration signal;;S [mV·sl]");
366 MHCamEvent evt1g(4, "ExtCalTm", "Extracted arrival times;;T [sl]");
367 MHCamEvent evt1h(6, "CalCalTm", "Calibrated arrival times;;T [sl]");
368
369 MHSectorVsTime hist1rmsb("ExtSigVsTm");
370 MHSectorVsTime hist1tmb("CalTmVsTm");
371 hist1rmsb.SetTitle("Extracted calibration vs event number;;S [mV·sl]");
372 hist1rmsb.SetType(0);
373 hist1tmb.SetTitle("Extracted arrival time vs event number;;T [sl]");
374 //hist1tmb.SetType(4);
375 hist1tmb.SetType(6);
376
377 MFillH fill1f(&evt1f, "MExtractedSignalCam", "FillExtSig");
378 MFillH fill1g(&evt1g, "MArrivalTimeCam", "FillExtTm");
379 MFillH fill1h(&evt1h, "MSignalCam", "FillCalTm");
380 MFillH fill1r(&hist1rmsb, "MExtractedSignalCam", "FillExtSigVsTm");
381 //MFillH fill1j(&hist1tmb, "MArrivalTimeCam", "FillExtTmVsTm");
382 MFillH fill1j(&hist1tmb, "MSignalCam", "FillCalTmVsTm");
383
384 fill1f.SetDrawOption("gaus");
385 fill1h.SetDrawOption("gaus");
386
387 // ---
388
389 MCalibrateDrsTimes calctm1a("CalibrateCalEvents");
390 calctm1a.SetNameUncalibrated("UncalibratedTimes");
391
392 MBadPixelsTreat treat1;
393 treat1.SetProcessPedestalRun(kFALSE);
394 treat1.SetProcessPedestalEvt(kFALSE);
395
396 // ---
397
398 MHCamEvent evt1c(6, "ExtCalTmShift", "Relative extracted arrival time of calibration pulse (w.r.t. event-median);;\\Delta T [ns]");
399 MHCamEvent evt1d(6, "CalCalTmShift", "Relative calibrated arrival time of calibration pulse (w.r.t. event-median);;\\Delta T [ns]");
400
401 evt1c.SetMedianShift();
402 evt1d.SetMedianShift();
403
404 MFillH fill1c(&evt1c, "UncalibratedTimes", "FillExtCalTm");
405 MFillH fill1d(&evt1d, "MSignalCam", "FillCalCalTm");
406 fill1d.SetDrawOption("gaus");
407
408 // ------------------ Setup eventloop and run analysis ---------------
409
410 tlist1.AddToList(&read1);
411 tlist1.AddToList(&apply1);
412 tlist1.AddToList(&drsapply1);
413 tlist1.AddToList(&cont1);
414 tlist1.AddToList(&extractor1b);
415 if (isinteg)
416 {
417 tlist1.AddToList(&extractor1c);
418 tlist1.AddToList(&fill1a);
419 }
420 tlist1.AddToList(&calctm1a);
421 tlist1.AddToList(&treat1);
422 tlist1.AddToList(&fill1f);
423 tlist1.AddToList(&fill1g);
424 tlist1.AddToList(&fill1h);
425 tlist1.AddToList(&fill1r);
426 tlist1.AddToList(&fill1j);
427 tlist1.AddToList(&fill1c);
428 tlist1.AddToList(&fill1d);
429
430 if (!loop1.Eventloop(max1))
431 return 10;
432
433 if (!loop1.GetDisplay())
434 return 11;
435
436 if (use_delays)
437 timecam.SetDelays(*evt1h.GetHist());
438
439 // ========================= Result ==================================
440
441 Double_t avgS = evt1f.GetHist()->GetMean();
442 Double_t medS = evt1f.GetHist()->GetMedian();
443 Double_t rmsS = evt1f.GetHist()->GetRMS();
444 Double_t maxS = evt1f.GetHist()->GetMaximum();
445
446 MArrayF der1(hpulse->GetNbinsX());
447 MArrayF der2(hpulse->GetNbinsX());
448
449 MExtralgoSpline spline(hpulse->GetArray()+1, hpulse->GetNbinsX(),
450 der1.GetArray(), der2.GetArray());
451 spline.SetRiseFallTime(rise_time_dat, fall_time_dat);
452 spline.SetExtractionType(type);
453 spline.SetHeightTm(heighttm);
454
455 spline.Extract(hpulse->GetMaximumBin()-1);
456
457 // The pulser signal is most probably around 400mV/9.5mV
458 // IntegraFixed 2/48 corresponds to roughly 215mV*50slices
459 Double_t scale = 1./spline.GetSignal();
460
461 MArrayD calib(1440);
462 for (int i=0; i<1440; i++)
463 {
464 Double_t g = gain->GetBinContent(i+1)>0.5 ? gain->GetBinContent(i+1) : 1;
465 if (evt1f.GetHist()->GetBinContent(i+1)>0 && !badpixels[i].IsUnsuitable())
466 calib[i] = lpref->GetBinContent(i+1) / evt1f.GetHist()->GetBinContent(i+1) / g;
467 }
468
469 gROOT->SetSelectedPad(0);
470 d->AddTab("PulseTemp");
471 gPad->SetGrid();
472 hpulse->SetNameTitle("Pulse", "Single p.e. pulse template");
473 hpulse->SetDirectory(0);
474 hpulse->SetLineColor(kBlack);
475 hpulse->DrawCopy();
476
477 TAxis *ax = hpulse->GetXaxis();
478
479 Double_t w = hpulse->GetBinWidth(1);
480 Double_t T = w*(spline.GetTime()+0.5) +ax->GetXmin();
481 Double_t H = w*(hpulse->GetMaximumBin()+0.5)+ax->GetXmin();
482
483 TLine line;
484 line.SetLineColor(kRed);
485 line.DrawLine(T-rise_time_dat*w, spline.GetHeight(),
486 T+fall_time_dat*w, spline.GetHeight());
487 line.DrawLine(T, spline.GetHeight()/4, T, 3*spline.GetHeight()/4);
488 line.DrawLine(T-rise_time_dat*w, 0,
489 T-rise_time_dat*w, spline.GetHeight());
490 line.DrawLine(T+fall_time_dat*w, 0,
491 T+fall_time_dat*w, spline.GetHeight());
492
493 TGraph gg;
494 for (int ix=1; ix<=hpulse->GetNbinsX(); ix++)
495 for (int i=0; i<10; i++)
496 {
497 Double_t x = hpulse->GetBinLowEdge(ix)+i*hpulse->GetBinWidth(ix)/10.;
498 gg.SetPoint(gg.GetN(), x+w/2, spline.EvalAt(ix-1+i/10.));
499 }
500
501 gg.SetLineColor(kBlue);
502 gg.SetMarkerColor(kBlue);
503 gg.SetMarkerStyle(kFullDotMedium);
504 gg.DrawClone("L");
505
506 gROOT->SetSelectedPad(0);
507 d->AddTab("CalConst");
508 MGeomCamFACT fact;
509 MHCamera hcalco(fact);
510 hcalco.SetName("CalConst");
511 hcalco.SetTitle(Form("Relative calibration constant [%.0f/pe]", 1./scale));
512 hcalco.SetCamContent(calib);
513 hcalco.SetAllUsed();
514 //hcalco.Scale(scale);
515 hcalco.DrawCopy();
516
517 // ======================================================
518
519 gLog << endl;
520 gLog.Separator("Extracting random pedestal");
521
522 MTaskList tlist3;
523
524 MParList plist3;
525 plist3.AddToList(&tlist3);
526 plist3.AddToList(&drscalib300);
527 plist3.AddToList(&badpixels);
528 plist3.AddToList(&timecam);
529
530 MEvtLoop loop3("DetermineRndmPed");
531 loop3.SetDisplay(d);
532 loop3.SetParList(&plist3);
533
534 // ------------------ Setup the tasks ---------------
535
536 MRawFitsRead read3;
537 read3.LoadMap(map);
538 read3.AddFile(pedfile);
539
540 MFillH fill3a(&hrate);
541
542 MContinue cont3("(MRawEvtHeader.GetTriggerID&0xff00)!=0x400", "SelectPed");
543
544 MGeomApply apply3;
545
546 MDrsCalibApply drsapply3;
547 drsapply3.SetRemoveSpikes(4);
548
549 //---
550
551 MExtractTimeAndChargeSpline extractor3;
552 extractor3.SetRange(first_slice, last_slice);
553 extractor3.SetRiseTimeHiGain(rise_time_dat);
554 extractor3.SetFallTimeHiGain(fall_time_dat);
555 extractor3.SetHeightTm(heighttm);
556 extractor3.SetChargeType(type);
557 extractor3.SetSaturationLimit(600000);
558 extractor3.SetNoiseCalculation(kTRUE);
559
560// MHCamEvent evt2a(0, "PedRdm", "Extracted Pedestal Signal;;S");
561
562// MFillH fill2a(&evt2a, "MExtractedSignalCam", "FillPedRndm");
563
564 // Use this for data, but not for calibration events
565// evt2a.SetErrorSpread(kFALSE);
566
567 /*
568 MCalibrateData conv3;
569 conv3.SetCalibrationMode(MCalibrateData::kNone);
570 conv3.SetPedestalFlag(MCalibrateData::kNo);
571 conv3.SetCalibConvMinLimit(0);
572 conv3.SetCalibConvMaxLimit(10000);
573 conv3.SetScaleFactor(scale);
574 */
575
576 MCalibrateFact conv3;
577 conv3.SetScale(scale);
578 conv3.SetCalibConst(calib);
579
580 MBadPixelsTreat treat3;
581 treat3.SetProcessPedestalRun(kFALSE);
582 treat3.SetProcessPedestalEvt(kFALSE);
583 treat3.SetProcessTimes(kFALSE);
584
585 MHCamEvent evt3b(0, "PedRdm","Interpolated random pedestal;;Signal [~phe]");
586 //evt2b.SetErrorSpread(kFALSE);
587
588 MFillH fill3b(&evt3b, "MSignalCam", "FillPedRdm");
589 fill3b.SetDrawOption("gaus");
590
591 // ------------------ Setup eventloop and run analysis ---------------
592
593 tlist3.AddToList(&read3);
594 tlist3.AddToList(&apply3);
595 tlist3.AddToList(&drsapply3);
596 tlist3.AddToList(&cont3);
597 tlist3.AddToList(&extractor3);
598// tlist3.AddToList(&fill3a);
599 tlist3.AddToList(&conv3);
600 tlist3.AddToList(&treat3);
601 tlist3.AddToList(&fill3b);
602
603 if (!loop3.Eventloop(max3))
604 return 14;
605
606 if (!loop3.GetDisplay())
607 return 15;
608
609 // ======================================================
610
611 gLog << endl;
612 gLog.Separator("Extracting pedestal");
613
614 MTaskList tlist4;
615
616 MParList plist4;
617 plist4.AddToList(&tlist4);
618 plist4.AddToList(&drscalib300);
619 plist4.AddToList(&badpixels);
620 plist4.AddToList(&timecam);
621
622 MEvtLoop loop4("DetermineExtractedPed");
623 loop4.SetDisplay(d);
624 loop4.SetParList(&plist4);
625
626 // ------------------ Setup the tasks ---------------
627
628 MRawFitsRead read4;
629 read4.LoadMap(map);
630 read4.AddFile(pedfile);
631
632 MContinue cont4("(MRawEvtHeader.GetTriggerID&0xff00)!=0x400", "SelectPed");
633
634 MGeomApply apply4;
635
636 MDrsCalibApply drsapply4;
637 drsapply4.SetRemoveSpikes(4);
638
639 MExtractTimeAndChargeSpline extractor4;
640 extractor4.SetRange(first_slice, last_slice);
641 extractor4.SetRiseTimeHiGain(rise_time_dat);
642 extractor4.SetFallTimeHiGain(fall_time_dat);
643 extractor4.SetHeightTm(heighttm);
644 extractor4.SetChargeType(type);
645 extractor4.SetSaturationLimit(600000);
646 extractor4.SetNoiseCalculation(kFALSE);
647
648 // MHCamEvent evt3a(0, "PedExt", "Extracted Pedestal Signal;;S");
649
650 // MFillH fill3a(&evt3a, "MExtractedSignalCam", "FillPedExt");
651
652 // Use this for data, but not for calibration events
653// evt3a.SetErrorSpread(kFALSE);
654/*
655 MCalibrateData conv4;
656 conv4.SetCalibrationMode(MCalibrateData::kNone);
657 conv4.SetPedestalFlag(MCalibrateData::kNo);
658 conv4.SetCalibConvMinLimit(0);
659 conv4.SetCalibConvMaxLimit(10000);
660 conv4.SetScaleFactor(scale);
661*/
662 MCalibrateFact conv4;
663 conv4.SetScale(scale);
664 conv4.SetCalibConst(calib);
665
666 MBadPixelsTreat treat4;
667 treat4.SetProcessPedestalRun(kFALSE);
668 treat4.SetProcessPedestalEvt(kFALSE);
669
670 MHCamEvent evt4b(0, "PedExt","Interpolated extracted pedestal;;Signal [~phe]");
671 //evt4b.SetErrorSpread(kFALSE);
672
673 MFillH fill4b(&evt4b, "MSignalCam", "FillPedExt");
674 fill4b.SetDrawOption("gaus");
675
676 // ------------------ Setup eventloop and run analysis ---------------
677
678 tlist4.AddToList(&read4);
679 tlist4.AddToList(&apply4);
680 tlist4.AddToList(&drsapply4);
681 tlist4.AddToList(&cont4);
682 tlist4.AddToList(&extractor4);
683// tlist4.AddToList(&fill4a);
684 tlist4.AddToList(&conv4);
685 tlist4.AddToList(&treat4);
686 tlist4.AddToList(&fill4b);
687
688 if (!loop4.Eventloop(max4))
689 return 15;
690
691 if (!loop4.GetDisplay())
692 return 16;
693
694 // ===================================================================
695
696 gLog << endl;
697 gLog.Separator("Extracting and calibration data");
698
699 MTaskList tlist5;
700
701 MParList plist5;
702 plist5.AddToList(&tlist5);
703 plist5.AddToList(&drscalib300);
704 plist5.AddToList(&badpixels);
705 plist5.AddToList(&timecam);
706
707 MEvtLoop loop5("CalibratingData");
708 loop5.SetDisplay(d);
709 loop5.SetParList(&plist5);
710
711 // ------------------ Setup the tasks ---------------
712
713 MRawFitsRead read5;
714 read5.LoadMap(map);
715 read5.AddFiles(iter);
716
717 MFillH fill5a(&hrate);
718
719 MGeomApply apply5;
720
721 MDrsCalibApply drsapply5;
722 drsapply5.SetRemoveSpikes(4);
723
724 MFDataPhrase filterdat("(MRawEvtHeader.GetTriggerID&0xff00)==0", "SelectDat");
725 MFDataPhrase filtercal("(MRawEvtHeader.GetTriggerID&0xff00)==0x100", "SelectCal");
726 MFDataPhrase filterped("(MRawEvtHeader.GetTriggerID&0xff00)==0x400", "SelectPed");
727 MFDataPhrase filterncl("(MRawEvtHeader.GetTriggerID&0xff00)!=0x100", "SelectNonCal");
728
729 //MContinue cont4("MRawEvtHeader.GetTriggerID!=4", "SelectData");
730
731 // ---
732
733 MExtractTimeAndChargeSpline extractor5dat;
734 extractor5dat.SetRange(first_slice, last_slice);
735 extractor5dat.SetRiseTimeHiGain(rise_time_dat);
736 extractor5dat.SetFallTimeHiGain(fall_time_dat);
737 extractor5dat.SetHeightTm(heighttm);
738 extractor5dat.SetChargeType(type);
739 extractor5dat.SetSaturationLimit(600000);
740 extractor5dat.SetNoiseCalculation(kFALSE);
741
742 MExtractTimeAndChargeSpline extractor5cal;
743 extractor5cal.SetRange(first_slice, last_slice);
744 extractor5cal.SetRiseTimeHiGain(rise_time_cal);
745 extractor5cal.SetFallTimeHiGain(fall_time_cal);
746 extractor5cal.SetHeightTm(heighttm);
747 extractor5cal.SetChargeType(type);
748 extractor5cal.SetSaturationLimit(600000);
749 extractor5cal.SetNoiseCalculation(kFALSE);
750
751 MExtractTimeAndChargeSpline extractor5tm("ExtractTM");
752 extractor5tm.SetRange(last_slice, 294);
753 extractor5tm.SetRiseTimeHiGain(1);
754 extractor5tm.SetFallTimeHiGain(1);
755 extractor5tm.SetHeightTm(0.5);
756 extractor5tm.SetChargeType(MExtralgoSpline::kAmplitudeRel);
757 extractor5tm.SetSaturationLimit(600000);
758 extractor5tm.SetNoiseCalculation(kFALSE);
759 extractor5tm.SetNameSignalCam("TimeMarkerAmplitude");
760 extractor5tm.SetNameTimeCam("TimeMarkerTime");
761
762 extractor5dat.SetFilter(&filterncl);
763 extractor5cal.SetFilter(&filtercal);
764 //extractor4tm.SetFilter(&filtercal);
765
766 // ---
767/*
768 MCalibrateData conv5;
769 conv5.SetCalibrationMode(MCalibrateData::kNone);
770 conv5.SetPedestalFlag(MCalibrateData::kNo);
771 conv5.SetCalibConvMinLimit(0);
772 conv5.SetCalibConvMaxLimit(10000);
773 conv5.SetScaleFactor(scale);
774*/
775 MCalibrateFact conv5;
776 conv5.SetScale(scale);
777 conv5.SetCalibConst(calib);
778
779 MCalibrateDrsTimes calctm5;
780 calctm5.SetNameUncalibrated("UncalibratedTimes");
781
782 MCalibrateDrsTimes calctm5tm("CalibrateTimeMarker");
783 calctm5tm.SetNameArrivalTime("TimeMarkerTime");
784 calctm5tm.SetNameUncalibrated("UncalTimeMarker");
785 calctm5tm.SetNameCalibrated("TimeMarker");
786 calctm5tm.SetTimeMarker();
787 //calctm4tm.SetFilter(&filtercal);
788
789 MHCamEvent evt5m(6, "ExtTm", "Extracted arrival times of calibration pulse;;\\Delta T [ns]");
790 MHCamEvent evt5n(6, "CalTm", "Calibrated arrival times of calibration pulse;;\\Delta T [ns]");
791 MHCamEvent evt5q(6, "ExtTmShift", "Relative extracted arrival times of calibration pulse (w.r.t. event-median);;\\Delta T [ns]");
792 MHCamEvent evt5r(6, "CalTmShift", "Relative calibrated arrival times of calibration pulse (w.r.t. event-median);;\\Delta T [ns]");
793 MHCamEvent evt5s(6, "ExtTM", "Extracted absolute time marker position;;T [sl]");
794 MHCamEvent evt5t(6, "CalTM", "Calibrated absolute time marker position;;T [ns]");
795 MHCamEvent evt5u(6, "ExtTMshift", "Relative extracted time marker position (w.r.t. event-median);;\\Delta T [ns]");
796 MHCamEvent evt5v(6, "CalTMshift", "Relative calibrated time marker position (w.r.t. event-median);;\\Delta T [ns]");
797 MHCamEvent evt5w(6, "ExtDiff", "Difference between extracted arrival time of time marker and calibration pulse;;\\Delta T [ns]");
798 MHCamEvent evt5x(6, "CalDiff", "Difference between calibrated arrival time of time marker and calibration pulse;;\\Delta T [ns]");
799
800 evt5w.SetNameSub("UncalibratedTimes");
801 evt5x.SetNameSub("MSignalCam");
802
803 evt5q.SetMedianShift();
804 evt5r.SetMedianShift();
805 evt5u.SetMedianShift();
806 evt5v.SetMedianShift();
807 //evt4w.SetMedianShift();
808 //evt4x.SetMedianShift();
809
810 MFillH fill5m(&evt5m, "UncalibratedTimes", "FillExtTm");
811 MFillH fill5n(&evt5n, "MSignalCam", "FillCalTm");
812 MFillH fill5q(&evt5q, "UncalibratedTimes", "FillExtTmShift");
813 MFillH fill5r(&evt5r, "MSignalCam" , "FillCalTmShift");
814 MFillH fill5s(&evt5s, "UncalTimeMarker", "FillExtTM");
815 MFillH fill5t(&evt5t, "TimeMarker", "FillCalTM");
816 MFillH fill5u(&evt5u, "UncalTimeMarker", "FillExtTMshift");
817 MFillH fill5v(&evt5v, "TimeMarker", "FillCalTMshift");
818 MFillH fill5w(&evt5w, "UncalTimeMarker", "FillExtDiff");
819 MFillH fill5x(&evt5x, "TimeMarker", "FillCalDiff");
820
821 fill5m.SetDrawOption("gaus");
822 fill5n.SetDrawOption("gaus");
823 fill5q.SetDrawOption("gaus");
824 fill5r.SetDrawOption("gaus");
825 //fill5s.SetDrawOption("gaus");
826 //fill5t.SetDrawOption("gaus");
827 //fill5u.SetDrawOption("gaus");
828 //fill5v.SetDrawOption("gaus");
829 //fill5w.SetDrawOption("gaus");
830 //fill5x.SetDrawOption("gaus");
831
832
833 MBadPixelsTreat treat5;
834 treat5.SetProcessPedestalRun(kFALSE);
835 treat5.SetProcessPedestalEvt(kFALSE);
836
837 MHSectorVsTime hist5cal("CalVsTm");
838 MHSectorVsTime hist5ped("PedVsTm");
839 hist5cal.SetTitle("Median calibrated calibration signal vs event number;;Signal [~phe]");
840 hist5ped.SetTitle("Median calibrated pedestal signal vs event number;;Signal [~phe]");
841 hist5cal.SetType(0);
842 hist5ped.SetType(0);
843 hist5cal.SetMinimum(0);
844 hist5ped.SetMinimum(0);
845 hist5cal.SetUseMedian();
846 hist5ped.SetUseMedian();
847 hist5cal.SetNameTime("MTime");
848 hist5ped.SetNameTime("MTime");
849
850 MFillH fill5cal(&hist5cal, "MSignalCam", "FillCalVsTm");
851 MFillH fill5ped(&hist5ped, "MSignalCam", "FillPedVsTm");
852 fill5cal.SetFilter(&filtercal);
853 fill5ped.SetFilter(&filterped);
854
855 MHCamEvent evt5b(0, "ExtSig", "Extracted signal;;S [mV·sl]");
856 MHCamEvent evt5c(0, "CalSig", "Calibrated and interpolated signal;;S [~phe]");
857 MHCamEvent evt5d(4, "ExtSigTm", "Extracted time;;T [sl]");
858 MHCamEvent evt5e(6, "CalSigTm", "Calibrated and interpolated time;;T [ns]");
859
860 MFillH fill5b(&evt5b, "MExtractedSignalCam", "FillExtSig");
861 MFillH fill5c(&evt5c, "MSignalCam", "FillCalSig");
862 MFillH fill5d(&evt5d, "MArrivalTimeCam", "FillExtTm");
863 MFillH fill5e(&evt5e, "MSignalCam", "FillCalTm");
864
865 fill5c.SetDrawOption("gaus");
866 fill5d.SetDrawOption("gaus");
867 fill5e.SetDrawOption("gaus");
868
869 /*
870 fill4b.SetFilter(&filterdat);
871 fill4c.SetFilter(&filterdat);
872 fill4d.SetFilter(&filterdat);
873 fill4e.SetFilter(&filterdat);
874 */
875
876 //MFSoftwareTrigger swtrig;
877 //MContinue contsw(&swtrig, "FilterSwTrigger", "Software trigger");
878 //contsw.SetInverted();
879
880 const TString fname(Form("s/([0-9]+_[0-9]+)[.]fits([.][fg]z)?$/%s\\/$1_C.root/",
881 MJob::Esc(outpath).Data()));
882
883 // The second rule is for the case reading raw-files!
884 MWriteRootFile write5(2, fname, "RECREATE", "Calibrated Data");
885 write5.AddContainer("MRawRunHeader", "RunHeaders");
886 write5.AddContainer("MGeomCam", "RunHeaders");
887 write5.AddContainer("MSignalCam", "Events");
888 write5.AddContainer("MTime", "Events");
889 write5.AddContainer("MRawEvtHeader", "Events");
890 //write.AddContainer("MTriggerPattern", "Events");
891
892 // ------------------ Setup histograms and fill tasks ----------------
893
894 MContinue test;
895 test.SetFilter(&filterncl);
896
897 MTaskList tlist5tm;
898 tlist5tm.AddToList(&extractor5tm);
899 tlist5tm.AddToList(&calctm5tm);
900 tlist5tm.AddToList(&fill5m);
901 tlist5tm.AddToList(&fill5n);
902 tlist5tm.AddToList(&fill5q);
903 tlist5tm.AddToList(&fill5r);
904 //tlist5tm.AddToList(&fill5s);
905 //tlist5tm.AddToList(&fill5t);
906 tlist5tm.AddToList(&fill5u);
907 tlist5tm.AddToList(&fill5v);
908 tlist5tm.AddToList(&fill5w);
909 tlist5tm.AddToList(&fill5x);
910 tlist5tm.SetFilter(&filtercal);
911
912 MTaskList tlist5dat;
913 tlist5dat.AddToList(&fill5b);
914 tlist5dat.AddToList(&fill5c);
915 tlist5dat.AddToList(&fill5d);
916 tlist5dat.AddToList(&fill5e);
917 tlist5dat.SetFilter(&filterdat);
918
919 tlist5.AddToList(&read5);
920 tlist5.AddToList(&apply5);
921 tlist5.AddToList(&drsapply5);
922 tlist5.AddToList(&filterncl);
923 //tlist5.AddToList(&test);
924 tlist5.AddToList(&filterdat);
925 tlist5.AddToList(&filtercal);
926 tlist5.AddToList(&filterped);
927 tlist5.AddToList(&fill5a);
928 tlist5.AddToList(&extractor5dat);
929 tlist5.AddToList(&extractor5cal);
930 tlist5.AddToList(&calctm5);
931 tlist5.AddToList(&tlist5tm);
932 tlist5.AddToList(&conv5);
933 tlist5.AddToList(&treat5);
934 tlist5.AddToList(&fill5ped);
935 tlist5.AddToList(&fill5cal);
936 tlist5.AddToList(&tlist5dat);
937 tlist5.AddToList(&write5);
938
939 if (!loop5.Eventloop(max4))
940 return 18;
941
942 if (!loop5.GetDisplay())
943 return 19;
944
945 TString title = "-- Calibrated signal #";
946 title += seq.GetSequence();
947 title += " (";
948 title += drsfile;
949 title += ") --";
950 d->SetTitle(title, kFALSE);
951
952 TString path;
953 path += Form("%s/20%6d_%03d-calibration.root", outpath,
954 seq.GetSequence()/1000, seq.GetSequence()%1000);
955
956 d->SaveAs(path);
957
958 return 0;
959}
Note: See TracBrowser for help on using the repository browser.