source: branches/Mars_MC/mjobs/MJCalibTest.cc@ 17430

Last change on this file since 17430 was 8999, checked in by tbretz, 16 years ago
*** empty log message ***
File size: 15.1 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Markus Gaug, 04/2004 <mailto:markus@ifae.es>
19!
20! Copyright: MAGIC Software Development, 2000-2005
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MJCalibTest
28//
29// If the flag SetDataCheckDisplay() is set, only the most important distributions
30// are displayed.
31// Otherwise, (default: SetNormalDisplay()), a good selection of plots is given
32//
33/////////////////////////////////////////////////////////////////////////////
34#include "MJCalibTest.h"
35
36#include <TFile.h>
37#include <TStyle.h>
38#include <TCanvas.h>
39#include <TSystem.h>
40
41#include "MLog.h"
42#include "MLogManip.h"
43
44#include "MParList.h"
45#include "MTaskList.h"
46#include "MTaskEnv.h"
47#include "MEvtLoop.h"
48
49#include "MHCamera.h"
50
51#include "MSignalCam.h"
52#include "MPedestalCam.h"
53#include "MPedPhotCam.h"
54#include "MBadPixelsCam.h"
55#include "MBadPixelsTreat.h"
56#include "MBadPixelsCalc.h"
57#include "MBadPixelsMerge.h"
58#include "MCalibrationChargeCam.h"
59#include "MCalibrationRelTimeCam.h"
60#include "MCalibrationQECam.h"
61#include "MCalibrationTestCam.h"
62#include "MCalibrationTestCalc.h"
63#include "MHCamEvent.h"
64#include "MHCalibrationTestCam.h"
65#include "MHCalibrationTestTimeCam.h"
66#include "MHCalibrationPix.h"
67
68#include "MReadMarsFile.h"
69#include "MRawFileRead.h"
70#include "MGeomApply.h"
71#include "MGeomCam.h"
72#include "MExtractTimeAndChargeSpline.h"
73#include "MExtractor.h"
74#include "MExtractTime.h"
75//#include "MExtractTimeFastSpline.h"
76#include "MFCosmics.h"
77#include "MContinue.h"
78#include "MFillH.h"
79#include "MCalibrateData.h"
80#include "MCalibrateRelTimes.h"
81
82#include "MTriggerPattern.h"
83#include "MTriggerPatternDecode.h"
84#include "MFTriggerPattern.h"
85
86#include "MStatusDisplay.h"
87
88ClassImp(MJCalibTest);
89
90using namespace std;
91
92// --------------------------------------------------------------------------
93//
94// Default constructor.
95//
96// Sets fUseCosmicsFilter to kTRUE, fExtractor to NULL, fTimeExtractor to NULL
97// fDisplay to kNormalDisplay
98//
99MJCalibTest::MJCalibTest(const char *name, const char *title)
100 : fUseCosmicsFilter(kTRUE), fExtractor(NULL), fTimeExtractor(NULL),
101 fDisplayType(kNormalDisplay), fGeometry("MGeomCamMagic")
102{
103 fName = name ? name : "MJCalibTest";
104 fTitle = title ? title : "Tool to extract, calibrate and test signals from a file";
105}
106
107
108void MJCalibTest::DisplayResult(MParList &plist)
109{
110 if (!fDisplay)
111 return;
112
113 //
114 // Get container from list
115 //
116 MGeomCam &geomcam = *(MGeomCam*) plist.FindObject("MGeomCam");
117 MHCalibrationTestCam &testcam = *(MHCalibrationTestCam*)plist.FindObject("MHCalibrationTestCam");
118
119 // Create histograms to display
120 MHCamera disp1 (geomcam, "Test;PhotoElectrons", "Mean equiv. phes");
121 MHCamera disp2 (geomcam, "Test;SigmaPhes", "Sigma equiv.phes");
122 MHCamera disp3 (geomcam, "Test;PhesPerArea", "Equiv. Phes per Area");
123 MHCamera disp4 (geomcam, "Test;SigmaPhotPerArea", "Sigma equiv. Phes per Area");
124 MHCamera disp5 (geomcam, "Test;Phot", "Calibrated Phes from Fit");
125 MHCamera disp6 (geomcam, "Test;PhotPerArea", "Calibrated Phes per Area from Fit");
126 MHCamera disp7 (geomcam, "Test;NotInterpolate", "Not interpolated pixels");
127 MHCamera disp8 (geomcam, "Test;DeviatingPhots", "Deviating Number Phes");
128 MHCamera disp9 (geomcam, "Test;Arr.Times", "Mean of calibrated Arr.Times");
129 MHCamera disp10(geomcam, "Test;SigmaArr.Times", "Sigma of calibrated Arr.Times");
130
131 // Fitted charge means and sigmas
132 disp1.SetCamContent(testcam, 0);
133 disp1.SetCamError( testcam, 1);
134 disp2.SetCamContent(testcam, 2);
135 disp2.SetCamError( testcam, 3);
136 disp3.SetCamContent(testcam, 7);
137 disp3.SetCamError( testcam, 8);
138 disp4.SetCamContent(testcam, 9);
139 disp4.SetCamError( testcam, 10);
140
141 disp5.SetCamContent(fTestCam, 0);
142 disp5.SetCamError( fTestCam, 1);
143 disp6.SetCamContent(fTestCam, 2);
144 disp6.SetCamError( fTestCam, 3);
145 disp7.SetCamError( fTestCam, 4);
146
147 disp8.SetCamError( fBadPixels, 22);
148
149 disp9.SetCamContent(fTestTimeCam, 0);
150 disp9.SetCamError( fTestTimeCam, 1);
151 disp10.SetCamContent(fTestTimeCam, 2);
152 disp10.SetCamError( fTestTimeCam, 3);
153
154
155 disp1.SetYTitle("Phes");
156 disp2.SetYTitle("\\sigma_{phe}");
157 disp3.SetYTitle("Phes per area [mm^{-2}]");
158 disp4.SetYTitle("\\sigma_{phe} per area [mm^{-2}]");
159
160 disp5.SetYTitle("Phes");
161 disp6.SetYTitle("Phes per area [mm^{-2}]");
162 disp7.SetYTitle("[1]");
163 disp8.SetYTitle("[1]");
164
165 disp9.SetYTitle("Mean Arr.Times [FADC units]");
166 disp10.SetYTitle("\\sigma_{t} [FADC units]");
167
168 TCanvas &c = fDisplay->AddTab("TestPhes");
169 c.Divide(4,4);
170
171 disp1.CamDraw(c, 1, 4, 2, 1);
172 disp2.CamDraw(c, 2, 4, 2, 1);
173 disp3.CamDraw(c, 3, 4, 1, 1);
174 disp4.CamDraw(c, 4, 4, 2, 1);
175
176 /*
177
178 TCanvas &c2 = fDisplay->AddTab("TestResult");
179 c2.Divide(2,4);
180
181 disp5.CamDraw(c2, 1, 2, 2, 1);
182 disp6.CamDraw(c2, 2, 2, 2, 1);
183
184 */
185
186 TCanvas &c3 = fDisplay->AddTab("TestDefects");
187 c3.Divide(2,2);
188
189 disp7.CamDraw(c3, 1, 2, 0);
190 disp8.CamDraw(c3, 2, 2, 0);
191
192 //
193 // Display times
194 //
195 TCanvas &c4 = fDisplay->AddTab("TestTimes");
196 c4.Divide(2,4);
197
198 disp9.CamDraw(c4, 1, 2, 5, 1);
199 disp10.CamDraw(c4, 2, 2, 5, 1);
200
201 return;
202
203}
204
205// --------------------------------------------------------------------------
206//
207// Retrieve the output file written by WriteResult()
208//
209const char* MJCalibTest::GetOutputFile() const
210{
211 return Form("%s/calib%08d.root", fPathOut.Data(), fSequence.GetSequence());
212}
213
214Bool_t MJCalibTest::ReadCalibration(TObjArray &l, MBadPixelsCam &cam, MExtractor* &ext1, MExtractor* &ext2, TString &geom) const
215{
216
217 const TString fname = GetOutputFile();
218
219 *fLog << inf << "Reading from file: " << fname << endl;
220
221 TFile file(fname, "READ");
222 if (!file.IsOpen())
223 {
224 *fLog << err << dbginf << "ERROR - Could not open file " << fname << endl;
225 return kFALSE;
226 }
227
228 TObject *o = file.Get("ExtractSignal");
229 if (o && !o->InheritsFrom(MExtractor::Class()))
230 {
231 *fLog << err << dbginf << "ERROR - ExtractSignal read from " << fname << " doesn't inherit from MExtractor!" << endl;
232 return kFALSE;
233 }
234 ext1 = o ? (MExtractor*)o->Clone() : NULL;
235
236 o = file.Get("ExtractTime");
237 if (o && !o->InheritsFrom(MExtractor::Class()))
238 {
239 *fLog << err << dbginf << "ERROR - ExtractTime read from " << fname << " doesn't inherit from MExtractor!" << endl;
240 return kFALSE;
241 }
242 ext2 = o ? (MExtractor*)o->Clone() : NULL;
243 if (!ext1 && !ext2)
244 {
245 *fLog << err << dbginf << "ERROR - Neither ExtractSignal nor ExrtractTime found in " << fname << "!" << endl;
246 return kFALSE;
247 }
248
249 o = file.Get("MGeomCam");
250 if (o && !o->InheritsFrom(MGeomCam::Class()))
251 {
252 *fLog << err << dbginf << "ERROR - MGeomCam read from " << fname << " doesn't inherit from MGeomCam!" << endl;
253 return kFALSE;
254 }
255 geom = o ? o->ClassName() : "";
256
257 TObjArray cont(l);
258 cont.Add(&cam);
259 return ReadContainer(cont);
260}
261
262Bool_t MJCalibTest::Process(MPedestalCam &pedcam)
263{
264 if (!fSequence.IsValid())
265 {
266 *fLog << err << "ERROR - Sequence invalid..." << endl;
267 return kFALSE;
268 }
269
270 *fLog << inf;
271 fLog->Separator(GetDescriptor());
272 *fLog << "Calculate calibration test from Sequence #";
273 *fLog << fSequence.GetSequence() << endl << endl;
274
275 CheckEnv();
276
277 *fLog << inf;
278 fLog->Separator(GetDescriptor());
279
280 *fLog << "Calibrate Calibration data from ";
281 *fLog << "Sequence #" << fSequence.GetSequence() << endl;
282 *fLog << endl;
283
284 MDirIter iter;
285
286 if (fSequence.IsValid())
287 {
288 if (fSequence.SetupCalRuns(iter, 0, !fSequence.IsMonteCarlo())<=0)
289 return kFALSE;
290 }
291
292 MCalibrationChargeCam calcam;
293 MCalibrationQECam qecam;
294 MCalibrationRelTimeCam tmcam;
295 MBadPixelsCam badpix;
296
297 TObjArray calibcont;
298 calibcont.Add(&calcam);
299 calibcont.Add(&qecam);
300 calibcont.Add(&tmcam);
301
302 MExtractor *extractor1=0;
303 MExtractor *extractor2=0;
304 TString geom;
305
306 if (!ReadCalibration(calibcont, badpix, extractor1, extractor2, geom))
307 {
308 *fLog << err << "Could not read calibration constants " << endl;
309 return kFALSE;
310 }
311
312 *fLog << all;
313 if (extractor1)
314 {
315 *fLog << underline << "Signal Extractor found in calibration file" << endl;
316 extractor1->Print();
317 *fLog << endl;
318 }
319 else
320 *fLog << inf << "No Signal Extractor: ExtractSignal in file." << endl;
321
322 if (extractor2)
323 {
324 *fLog << underline << "Time Extractor found in calibration file" << endl;
325 extractor2->Print();
326 *fLog << endl;
327 }
328 else
329 *fLog << inf << "No Time Extractor: ExtractTime in file." << endl;
330
331 if (!geom.IsNull())
332 *fLog << inf << "Camera geometry found in file: " << geom << endl;
333 else
334 *fLog << inf << "No Camera geometry found using default <MGeomCamMagic>" << endl;
335
336 if (fExtractor)
337 extractor1 = fExtractor;
338 if (fTimeExtractor)
339 extractor2 = fTimeExtractor;
340
341 // Setup Lists
342 MParList plist;
343 plist.AddToList(this); // take care of fDisplay!
344 plist.AddToList(&fTestCam);
345 plist.AddToList(&fTestTimeCam);
346 plist.AddToList(&badpix);
347 plist.AddToList(&pedcam);
348 plist.AddToList(&calcam);
349 plist.AddToList(&qecam);
350 plist.AddToList(&tmcam);
351
352 MSignalCam cerphot;
353 MPedPhotCam pedphot;
354 MHCalibrationTestCam testcam;
355
356 plist.AddToList(&cerphot);
357 plist.AddToList(&pedphot);
358 plist.AddToList(&testcam);
359
360 pedcam.SetName("MPedestalFundamental");
361
362 MTaskList tlist;
363 plist.AddToList(&tlist);
364
365 // Setup Task-lists
366 MRawFileRead rawread(NULL);
367 MReadMarsFile read("Events");
368 read.DisableAutoScheme();
369
370 if (!fSequence.IsMonteCarlo())
371 rawread.AddFiles(iter);
372 else
373 static_cast<MRead&>(read).AddFiles(iter);
374
375 // Check for interleaved events
376 // This will make that for data with version less than 5, where trigger
377 // patterns were not yet correct, all the events in the file will be
378 // processed. For those data there are no interleaved calibration events,
379 // so it makes no sense to run this test on a _D_ file. So we assume it
380 // will be a _C_ file, and accept all events.
381 MTriggerPatternDecode decode;
382 MFTriggerPattern fcalib;
383 fcalib.DenyCalibration();
384 fcalib.SetDefault(kFALSE);
385 MContinue conttp(&fcalib, "ContTrigPattern");
386
387 MGeomApply apply; // Only necessary to craete geometry
388 if (!geom.IsNull())
389 apply.SetGeometry(geom);
390 MBadPixelsMerge merge(&badpix);
391
392// MExtractTimeAndChargeSlidingWindow extrsw;
393// MExtractTimeFastSpline extime;
394
395 MExtractTimeAndChargeSpline spline;
396
397 MTaskEnv taskenv1("ExtractSignal");
398 MTaskEnv taskenv2("ExtractTime");
399
400 if (extractor1)
401 taskenv1.SetDefault(extractor1);
402
403 if (extractor2)
404 taskenv2.SetDefault(fTimeExtractor);
405 else if (!(extractor1->InheritsFrom("MExtractTimeAndCharge")))
406 {
407 spline.SetWindowSize(8,8);
408 taskenv2.SetDefault(&spline);
409 *fLog << warn << GetDescriptor()
410 << ": No extractor has been chosen, take default MExtractTimeAndChargeSlidingWindow " << endl;
411 }
412
413
414 MCalibrateData photcalc;
415 photcalc.AddPedestal("Fundamental");
416 photcalc.SetCalibrationMode(MCalibrateData::kFfactor);
417 photcalc.SetPedestalFlag(MCalibrateData::kEvent);
418 photcalc.SetSignalType(MCalibrateData::kPhe);
419
420 MCalibrateRelTimes caltimes;
421 MBadPixelsCalc badcalc;
422 MBadPixelsTreat badtreat;
423 badtreat.SetProcessTimes(kFALSE);
424
425 badcalc.SetNamePedPhotCam("MPedPhotFundamental");
426 badtreat.SetUseInterpolation();
427 badtreat.AddNamePedPhotCam("MPedPhotFundamental");
428
429 MCalibrationTestCalc testcalc;
430
431 MHCamEvent evt0(0,"Signal", "Un-Calibrated Signal;;S [FADC cnts]" );
432 MHCamEvent evt1(0,"CalSig", "Cal. and Interp. Sig. by Pixel Size Ratio;;S [phe]");
433 MHCamEvent evt2(6,"Times" , "Arrival Time;;T [slice]");
434
435 MFillH fill0(&evt0, "MExtractedSignalCam", "FillUncalibrated");
436 MFillH fill1(&evt1, "MSignalCam", "FillCalibrated");
437 MFillH fill2(&evt2, "MSignalCam", "FillTimes");
438
439 MFillH fillcam("MHCalibrationTestCam", "MSignalCam" ,"FillTest");
440 MFillH filltme("MHCalibrationTestTimeCam", "MSignalCam", "FillTestTime");
441 fillcam.SetBit(MFillH::kDoNotDisplay);
442 filltme.SetBit(MFillH::kDoNotDisplay);
443
444 MFCosmics cosmics;
445 cosmics.SetNamePedestalCam("MPedestalFundamental");
446 MContinue contcos(&cosmics,"ContCosmics");
447
448 tlist.AddToList(&read);
449 tlist.AddToList(&decode);
450 tlist.AddToList(&apply);
451 tlist.AddToList(&merge);
452 // tlist.AddToList(&conttp);
453 tlist.AddToList(&taskenv1);
454 if (!extractor1->InheritsFrom("MExtractTimeAndCharge"))
455 tlist.AddToList(&taskenv2);
456
457 tlist.AddToList(&contcos);
458 tlist.AddToList(&fill0);
459 tlist.AddToList(&photcalc);
460 tlist.AddToList(&caltimes);
461 tlist.AddToList(&badcalc);
462 tlist.AddToList(&badtreat);
463 tlist.AddToList(&fill1);
464 tlist.AddToList(&fill2);
465 tlist.AddToList(&fillcam);
466 tlist.AddToList(&filltme);
467 tlist.AddToList(&testcalc);
468
469 // Create and setup the eventloop
470 MEvtLoop evtloop(fName);
471 evtloop.SetParList(&plist);
472 evtloop.SetDisplay(fDisplay);
473 evtloop.SetLogStream(fLog);
474
475 // Execute first analysis
476 if (!evtloop.Eventloop())
477 {
478 *fLog << err << GetDescriptor() << ": Failed." << endl;
479 return kFALSE;
480 }
481
482 if (fIsPixelCheck)
483 {
484 MHCalibrationTestCam *hcam = (MHCalibrationTestCam*)plist.FindObject("MHCalibrationTestCam");
485 MHCalibrationPix &pix1 = (*hcam)[fCheckedPixId];
486 pix1.DrawClone("");
487
488 MHCalibrationTestTimeCam *hccam = (MHCalibrationTestTimeCam*)plist.FindObject("MHCalibrationTestTimeCam");
489 MHCalibrationPix &pix11 = (*hccam)[fCheckedPixId];
490 pix11.DrawClone("");
491 }
492
493 DisplayResult(plist);
494
495 if (!WriteResult())
496 return kFALSE;
497
498 *fLog << inf << GetDescriptor() << ": Done." << endl;
499
500 return kTRUE;
501}
502
503Bool_t MJCalibTest::WriteResult()
504{
505
506 if (fPathOut.IsNull())
507 return kTRUE;
508
509 const TString oname(GetOutputFile());
510
511 *fLog << inf << "Writing to file: " << oname << endl;
512
513 TFile file(oname, "UPDATE");
514
515 if (fDisplay && fDisplay->Write()<=0)
516 {
517 *fLog << err << "Unable to write MStatusDisplay to " << oname << endl;
518 return kFALSE;
519 }
520
521 if (fTestCam.Write()<=0)
522 {
523 *fLog << err << "Unable to write MCalibrationTestCam to " << oname << endl;
524 return kFALSE;
525 }
526
527 if (fTestTimeCam.Write()<=0)
528 {
529 *fLog << err << "Unable to write MCalibrationTestCam to " << oname << endl;
530 return kFALSE;
531 }
532
533 return kTRUE;
534
535}
Note: See TracBrowser for help on using the repository browser.