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

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