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

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