source: branches/Corsika7405Compatibility/fact/analysis/callisto_data.C@ 19365

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