source: branches/removing_cpp11_features/fact/analysis/callisto_data.C@ 19848

Last change on this file since 19848 was 18078, checked in by Daniela Dorner, 10 years ago
version of callisto.C only calibrating the data (used for 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.