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

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