source: branches/Corsika7500Compatibility/fact/analysis/callisto_data.C@ 19366

Last change on this file since 19366 was 18270, checked in by Daniela Dorner, 10 years ago
changed to usage of general drs-time-calib file and implemented usage of delays
File size: 15.5 KB
Line 
1#include "MLogManip.h"
2
3int callisto_data(const char *datafile, const char *drsfile, const char *drstime, const char *delays, const char *outpath=".")
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 //const char *pulse_template = "template-pulse.root";
17
18 // ------------------------------------------------------
19
20 MStatusDisplay *d = new MStatusDisplay;
21
22 MBadPixelsCam badpixels;
23 badpixels.InitSize(1440);
24 badpixels[ 424].SetUnsuitable(MBadPixelsPix::kUnsuitable);
25 badpixels[ 583].SetUnsuitable(MBadPixelsPix::kUnsuitable);
26 badpixels[ 830].SetUnsuitable(MBadPixelsPix::kUnsuitable);
27 badpixels[ 923].SetUnsuitable(MBadPixelsPix::kUnsuitable);
28 badpixels[1208].SetUnsuitable(MBadPixelsPix::kUnsuitable);
29 badpixels[1399].SetUnsuitable(MBadPixelsPix::kUnsuitable);
30 // Twin pixel
31 // 113
32 // 115
33 // 354
34 // 423
35 // 1195
36 // 1393
37
38 // ------------------------------------------------------
39
40 // nominal gain as extracted with gain extraction
41 // That is the number which describes
42 // the pulse amplitude in the file!
43 double gain = 241.;
44
45 // integrate around maximum or from leading edge on?
46 const bool maximum = true;
47
48 // ------------------------------------------------------
49
50 // Calib: 51 / 90 / 197 (20% TH)
51 // Data: 52 / 64 / 104 (20% TH)
52
53 // Extraction range in slices. It will always(!) contain the full range
54 // of integration
55 const int first_slice = 20; // 10ns
56 const int last_slice = 250; // 125ns
57
58 // Note that rise and fall time mean different things whether you use IntegralFixed or IntegraRel:
59 //
60 // IntegralFixed:
61 // * fRiseTime: Number of slices left from arrival time
62 // * fFallTime: Number of slices right from arrival time
63 // IntegralRel:
64 // * fRiseTime: Number of slices left from maximum time
65 // * fFallTime: Number of slices right from maximum time
66 //
67 const int rise_time_cal = maximum ? 40 : 10;
68 const int fall_time_cal = maximum ? 120 : 160;
69
70 const int rise_time_dat = maximum ? 5 : 0;
71 const int fall_time_dat = maximum ? 25 : 30;
72
73 // Extraction type: Extract integral and half leading edge
74
75 const int type = maximum ? (MExtralgoSpline::kIntegralRel) : (MExtralgoSpline::kIntegralFixed);
76 //const int type = MExtralgoSpline::kIntegralFixed;
77
78 const double heighttm = 0.5; // IntegralAbs { 1.5pe * 9.6mV/pe } / IntegralRel { 0.5 }
79
80 Long_t max = 0; // All
81
82 // ======================================================
83
84 if (map && gSystem->AccessPathName(map, kFileExists))
85 {
86 gLog << err << "ERROR - Cannot access mapping file '" << map << "'" << endl;
87 return 11;
88 }
89
90 // -------------------------------------------------------
91
92 MDrsCalibrationTime timecam;
93 if (!timecam.ReadFits(drstime))
94 {
95 gLog << err << "ERROR - Could not get MDrsCalibrationTime from " << drstime << endl;
96 return 21;
97 }
98
99 if (delays)
100 {
101 TGraph g(delays);
102 if (g.GetN()!=1440)
103 {
104 gLog << err << "Error reading file " << delays << endl;
105 return 22;
106 }
107
108 timecam.SetDelays(g);
109 }
110
111 // ------------------------------------------------------
112
113 MDrsCalibration drscalib300;
114 if (!drscalib300.ReadFits(drsfile))
115 return 31;
116
117 // -------------------------------------------------------
118
119 if (delays)
120 {
121 TGraph g(delays);
122 if (g.GetN()!=1440)
123 {
124 gLog << err << "Error reading file " << delays << endl;
125 return 41;
126 }
127
128 timecam.SetDelays(g);
129 }
130
131 // -------------------------------------------------------
132
133 // Range in which function/spline is defined
134 int minx = -10;
135 int maxx = 125;
136
137 // Number of points
138 int N = (maxx-minx)*2; // two points per nanosecond (2GHz)
139
140 TF1 f("pulse", "[0]*(1-exp(-[1]*x*x))*exp(-[2]*x)*(x>0)", float(minx), float(maxx));
141 f.SetNpx(N);
142
143 // Normalize pulse
144 f.SetParameters(1.9*gain/30, 0.10, 0.053);
145
146 // Convert to graph
147 TGraph g(&f);
148
149 // Convert to float
150 MArrayF data(g.GetN());
151 for (int i=0; i<g.GetN(); i++)
152 data[i] = g.GetY()[i];
153
154 // Define spline
155 MArrayF der1(g.GetN());
156 MArrayF der2(g.GetN());
157
158 MExtralgoSpline spline(data.GetArray(), data.GetSize(),
159 der1.GetArray(), der2.GetArray());
160
161 spline.SetRiseFallTime(rise_time_dat, fall_time_dat);
162 spline.SetExtractionType(type);
163 spline.SetHeightTm(heighttm);
164
165 // Estimate where the maximum is and extract signal
166 Long64_t maxi = TMath::LocMax(g.GetN(), g.GetY());
167 spline.Extract(maxi);
168
169 // Scale factor for signal extraction
170 double scale = 1./spline.GetSignal();
171
172 // ======================================================
173
174 gLog.Separator("Callisto");
175 gLog << all;
176 gLog << "Data file: " << datafile << '\n';
177 gLog << "DRS 300: " << drsfile << '\n';
178 gLog << "DRS Time: " << drstime << '\n';
179 gLog << "Delays: " << delays << '\n';
180 gLog << "Outpath: " << outpath << '\n';
181 gLog << endl;
182
183 // ------------------------------------------------------
184
185 // Plot the trigger pattern rates vs. run-number
186 MH3 hrate("MRawRunHeader.GetFileID", "MRawEvtHeader.GetTriggerID&0xff00");
187 hrate.SetWeight("1./TMath::Max(MRawRunHeader.GetRunLength,1)");
188 hrate.SetName("Rate");
189 hrate.SetTitle("Event rate [Hz];File Id;Trigger Type;");
190 hrate.InitLabels(MH3::kLabelsXY);
191 hrate.DefineLabelY( 0, "Data"); // What if TriggerID==0 already???
192 hrate.DefineLabelY(0x100, "Cal");
193 hrate.DefineLabelY(0x400, "Ped");
194 // hrate.DefaultLabelY("ERROR");
195
196 gStyle->SetOptFit(kTRUE);
197
198 // ===================================================================
199
200 gLog << endl;
201 gLog.Separator("Calibration data");
202
203 MTaskList tlist5;
204
205 MParList plist5;
206 plist5.AddToList(&tlist5);
207 plist5.AddToList(&drscalib300);
208 plist5.AddToList(&badpixels);
209 plist5.AddToList(&timecam);
210
211 MEvtLoop loop5("CalibratingData");
212 loop5.SetDisplay(d);
213 loop5.SetParList(&plist5);
214
215 // ------------------ Setup the tasks ---------------
216
217 MRawFitsRead read5(datafile);
218 read5.LoadMap(map);
219
220 MFillH fill5a(&hrate);
221
222 MGeomApply apply5;
223
224 MDrsCalibApply drsapply5;
225
226 MFDataPhrase filterdat("(MRawEvtHeader.GetTriggerID&0xff00)==0", "SelectDat");
227 MFDataPhrase filtercal("(MRawEvtHeader.GetTriggerID&0xff00)==0x100", "SelectCal");
228 MFDataPhrase filterped("(MRawEvtHeader.GetTriggerID&0xff00)==0x400", "SelectPed");
229 MFDataPhrase filterncl("(MRawEvtHeader.GetTriggerID&0xff00)!=0x100", "SelectNonCal");
230
231 //MContinue cont4("MRawEvtHeader.GetTriggerID!=4", "SelectData");
232
233 // ---
234
235 MExtractTimeAndChargeSpline extractor5dat;
236 extractor5dat.SetRange(first_slice, last_slice);
237 extractor5dat.SetRiseTimeHiGain(rise_time_dat);
238 extractor5dat.SetFallTimeHiGain(fall_time_dat);
239 extractor5dat.SetHeightTm(heighttm);
240 extractor5dat.SetChargeType(type);
241 extractor5dat.SetSaturationLimit(600000);
242 extractor5dat.SetNoiseCalculation(kFALSE);
243
244 MExtractTimeAndChargeSpline extractor5cal;
245 extractor5cal.SetRange(first_slice, last_slice);
246 extractor5cal.SetRiseTimeHiGain(rise_time_cal);
247 extractor5cal.SetFallTimeHiGain(fall_time_cal);
248 extractor5cal.SetHeightTm(heighttm);
249 extractor5cal.SetChargeType(type);
250 extractor5cal.SetSaturationLimit(600000);
251 extractor5cal.SetNoiseCalculation(kFALSE);
252
253 MExtractTimeAndChargeSpline extractor5tm("ExtractTM");
254 extractor5tm.SetRange(last_slice, 294);
255 extractor5tm.SetRiseTimeHiGain(1);
256 extractor5tm.SetFallTimeHiGain(1);
257 extractor5tm.SetHeightTm(0.5);
258 extractor5tm.SetChargeType(MExtralgoSpline::kAmplitudeRel);
259 extractor5tm.SetSaturationLimit(600000);
260 extractor5tm.SetNoiseCalculation(kFALSE);
261 extractor5tm.SetNameSignalCam("TimeMarkerAmplitude");
262 extractor5tm.SetNameTimeCam("TimeMarkerTime");
263
264 //extractor5dat.SetFilter(&filterncl);
265 //extractor5cal.SetFilter(&filtercal);
266 //extractor4tm.SetFilter(&filtercal);
267
268 // ---
269 MCalibrateFact conv5;
270 conv5.SetScale(scale);
271 //conv5.SetCalibConst(calib);
272
273 MCalibrateDrsTimes calctm5;
274 calctm5.SetNameUncalibrated("UncalibratedTimes");
275
276 MCalibrateDrsTimes calctm5tm("CalibrateTimeMarker");
277 calctm5tm.SetNameArrivalTime("TimeMarkerTime");
278 calctm5tm.SetNameUncalibrated("UncalTimeMarker");
279 calctm5tm.SetNameCalibrated("TimeMarker");
280 calctm5tm.SetTimeMarker();
281 //calctm4tm.SetFilter(&filtercal);
282
283 MHCamEvent evt5m(6, "ExtTm", "Extracted arrival times of calibration pulse;;\\Delta T [ns]");
284 MHCamEvent evt5n(6, "CalTm", "Calibrated arrival times of calibration pulse;;\\Delta T [ns]");
285 MHCamEvent evt5q(6, "ExtTmShift", "Relative extracted arrival times of calibration pulse (w.r.t. event-median);;\\Delta T [ns]");
286 MHCamEvent evt5r(6, "CalTmShift", "Relative calibrated arrival times of calibration pulse (w.r.t. event-median);;\\Delta T [ns]");
287 MHCamEvent evt5s(6, "ExtTM", "Extracted absolute time marker position;;T [sl]");
288 MHCamEvent evt5t(6, "CalTM", "Calibrated absolute time marker position;;T [ns]");
289 MHCamEvent evt5u(6, "ExtTMshift", "Relative extracted time marker position (w.r.t. event-median);;\\Delta T [ns]");
290 MHCamEvent evt5v(6, "CalTMshift", "Relative calibrated time marker position (w.r.t. event-median);;\\Delta T [ns]");
291 MHCamEvent evt5w(6, "ExtDiff", "Difference between extracted arrival time of time marker and calibration pulse;;\\Delta T [ns]");
292 MHCamEvent evt5x(6, "CalDiff", "Difference between calibrated arrival time of time marker and calibration pulse;;\\Delta T [ns]");
293
294 evt5w.SetNameSub("UncalibratedTimes");
295 evt5x.SetNameSub("MSignalCam");
296
297 evt5q.SetMedianShift();
298 evt5r.SetMedianShift();
299 evt5u.SetMedianShift();
300 evt5v.SetMedianShift();
301 //evt4w.SetMedianShift();
302 //evt4x.SetMedianShift();
303
304 MFillH fill5m(&evt5m, "UncalibratedTimes", "FillExtTm");
305 MFillH fill5n(&evt5n, "MSignalCam", "FillCalTm");
306 MFillH fill5q(&evt5q, "UncalibratedTimes", "FillExtTmShift");
307 MFillH fill5r(&evt5r, "MSignalCam" , "FillCalTmShift");
308 //MFillH fill5s(&evt5s, "UncalTimeMarker", "FillExtTM");
309 //MFillH fill5t(&evt5t, "TimeMarker", "FillCalTM");
310 MFillH fill5u(&evt5u, "UncalTimeMarker", "FillExtTMshift");
311 MFillH fill5v(&evt5v, "TimeMarker", "FillCalTMshift");
312 MFillH fill5w(&evt5w, "UncalTimeMarker", "FillExtDiff");
313 MFillH fill5x(&evt5x, "TimeMarker", "FillCalDiff");
314
315 fill5m.SetDrawOption("gaus");
316 fill5n.SetDrawOption("gaus");
317 fill5q.SetDrawOption("gaus");
318 fill5r.SetDrawOption("gaus");
319 //fill5s.SetDrawOption("gaus");
320 //fill5t.SetDrawOption("gaus");
321 //fill5u.SetDrawOption("gaus");
322 //fill5v.SetDrawOption("gaus");
323 //fill5w.SetDrawOption("gaus");
324 //fill5x.SetDrawOption("gaus");
325
326
327 MBadPixelsTreat treat5;
328 treat5.SetProcessPedestalRun(kFALSE);
329 treat5.SetProcessPedestalEvt(kFALSE);
330
331 MHSectorVsTime hist5cal("CalVsTm");
332 MHSectorVsTime hist5ped("PedVsTm");
333 hist5cal.SetTitle("Median calibrated calibration signal vs event number;;Signal [~phe]");
334 hist5ped.SetTitle("Median calibrated pedestal signal vs event number;;Signal [~phe]");
335 hist5cal.SetType(0);
336 hist5ped.SetType(0);
337 hist5cal.SetMinimum(0);
338 hist5ped.SetMinimum(0);
339 hist5cal.SetUseMedian();
340 hist5ped.SetUseMedian();
341 hist5cal.SetNameTime("MTime");
342 hist5ped.SetNameTime("MTime");
343
344 MFillH fill5cal(&hist5cal, "MSignalCam", "FillCalVsTm");
345 MFillH fill5ped(&hist5ped, "MSignalCam", "FillPedVsTm");
346 fill5cal.SetFilter(&filtercal);
347 fill5ped.SetFilter(&filterped);
348
349 MHCamEvent evt5b(0, "ExtSig", "Extracted signal;;S [mV·sl]");
350 MHCamEvent evt5c(0, "CalSig", "Calibrated and interpolated signal;;S [~phe]");
351 MHCamEvent evt5d(4, "ExtSigTm", "Extracted time;;T [sl]");
352 MHCamEvent evt5e(6, "CalSigTm", "Calibrated and interpolated time;;T [ns]");
353
354 MFillH fill5b(&evt5b, "MExtractedSignalCam", "FillExtSig");
355 MFillH fill5c(&evt5c, "MSignalCam", "FillCalSig");
356 MFillH fill5d(&evt5d, "MArrivalTimeCam", "FillExtTm");
357 MFillH fill5e(&evt5e, "MSignalCam", "FillCalTm");
358
359 fill5c.SetDrawOption("gaus");
360 fill5d.SetDrawOption("gaus");
361 fill5e.SetDrawOption("gaus");
362
363 /*
364 fill4b.SetFilter(&filterdat);
365 fill4c.SetFilter(&filterdat);
366 fill4d.SetFilter(&filterdat);
367 fill4e.SetFilter(&filterdat);
368 */
369
370 //MFSoftwareTrigger swtrig;
371 //MContinue contsw(&swtrig, "FilterSwTrigger", "Software trigger");
372 //contsw.SetInverted();
373
374 TString fname = gSystem->ConcatFileName(outpath, gSystem->BaseName(datafile));
375 fname.ReplaceAll(".fits.gz", "_C.root");
376 fname.ReplaceAll(".fits.fz", "_C.root");
377 fname.ReplaceAll(".fits", "_C.root");
378
379 // The second rule is for the case reading raw-files!
380 MWriteRootFile write5(fname, "RECREATE", "Calibrated Data", 2);
381 write5.AddContainer("MRawRunHeader", "RunHeaders");
382 write5.AddContainer("MGeomCam", "RunHeaders");
383 write5.AddContainer("MSignalCam", "Events");
384 write5.AddContainer("MTime", "Events");
385 write5.AddContainer("MRawEvtHeader", "Events");
386 //write.AddContainer("MTriggerPattern", "Events");
387
388 // ------------------ Setup histograms and fill tasks ----------------
389
390 //MTaskList tlist5tm;
391 //tlist5tm.AddToList(&extractor5tm);
392 //tlist5tm.AddToList(&calctm5tm);
393 //tlist5tm.AddToList(&fill5m);
394 //tlist5tm.AddToList(&fill5n);
395 //tlist5tm.AddToList(&fill5q);
396 //tlist5tm.AddToList(&fill5r);
397 //tlist5tm.AddToList(&fill5s);
398 //tlist5tm.AddToList(&fill5t);
399 //tlist5tm.AddToList(&fill5u);
400 //tlist5tm.AddToList(&fill5v);
401 //tlist5tm.AddToList(&fill5w);
402 //tlist5tm.AddToList(&fill5x);
403 //tlist5tm.SetFilter(&filtercal);
404
405 //MTaskList tlist5dat;
406 //tlist5dat.AddToList(&fill5b);
407 //tlist5dat.AddToList(&fill5c);
408 //tlist5dat.AddToList(&fill5d);
409 //tlist5dat.AddToList(&fill5e);
410 //tlist5dat.SetFilter(&filterdat);
411
412 tlist5.AddToList(&read5);
413 tlist5.AddToList(&apply5);
414 tlist5.AddToList(&drsapply5);
415 //tlist5.AddToList(&filterncl);
416 //tlist5.AddToList(&filterdat);
417 //tlist5.AddToList(&filtercal);
418 //tlist5.AddToList(&filterped);
419 //tlist5.AddToList(&fill5a);
420 tlist5.AddToList(&extractor5dat);
421 //tlist5.AddToList(&extractor5cal);
422 tlist5.AddToList(&calctm5);
423 //tlist5.AddToList(&tlist5tm);
424 tlist5.AddToList(&conv5);
425 tlist5.AddToList(&treat5);
426 //tlist5.AddToList(&fill5ped);
427 //tlist5.AddToList(&fill5cal);
428 //tlist5.AddToList(&tlist5dat);
429 tlist5.AddToList(&write5);
430
431 if (!loop5.Eventloop(max))
432 return 18;
433
434 if (!loop5.GetDisplay())
435 return 19;
436
437 // ============================================================
438
439 TFile *ofile = gROOT->GetListOfFiles()->FindObject(fname);
440 if (!ofile || !ofile->IsOpen() || ofile->IsZombie())
441 {
442 gLog << err << "File " << fname << " not found" << endl;
443 return 20;
444 }
445
446 TString title = "-- Calibrated signal [";
447 title += datafile;
448 title += "] --";
449 d->SetTitle(title, kFALSE);
450
451 ofile->cd();
452 d->Write();
453
454 return 0;
455}
Note: See TracBrowser for help on using the repository browser.