source: releases/Mars.2014.05.26/fact/analysis/callisto_data.C@ 18678

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