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

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