source: branches/MarsISDCBranchBasedOn17887/fact/analysis/callisto.C@ 19084

Last change on this file since 19084 was 17885, checked in by tbretz, 11 years ago
Added calls which allow to use the sequence number instead of the sequence file.
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")
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}
947
948int callisto(const ULong64_t seqnum, const char *outpath = "output")
949{
950 UInt_t night = seqnum/1000;
951 UInt_t num = seqnum%1000;
952
953 TString file = Form("/scratch/fact/sequences/%04d/%02d/%02d/%06d_%03d.seq",
954 night/10000, (night/100)%100, night%100, num);
955
956 return callisto(file.Data(), outpath);
957}
Note: See TracBrowser for help on using the repository browser.