source: trunk/MagicSoft/Mars/mjobs/MJCalibTest.cc@ 6193

Last change on this file since 6193 was 6193, checked in by gaug, 20 years ago
*** empty log message ***
File size: 17.7 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 "MRunIter.h"
45#include "MParList.h"
46#include "MTaskList.h"
47#include "MTaskEnv.h"
48#include "MEvtLoop.h"
49
50#include "MHCamera.h"
51
52#include "MPedestalCam.h"
53#include "MBadPixelsCam.h"
54#include "MBadPixelsTreat.h"
55#include "MBadPixelsCalc.h"
56#include "MBadPixelsMerge.h"
57#include "MCerPhotEvt.h"
58#include "MArrivalTime.h"
59#include "MCalibrationChargeCam.h"
60#include "MCalibrationRelTimeCam.h"
61#include "MCalibrationQECam.h"
62#include "MCalibrationTestCam.h"
63#include "MCalibrationTestCalc.h"
64#include "MHCamEvent.h"
65#include "MHCalibrationTestCam.h"
66
67#include "MReadMarsFile.h"
68#include "MRawFileRead.h"
69#include "MGeomApply.h"
70#include "MGeomCam.h"
71#include "MExtractTimeAndChargeSlidingWindow.h"
72#include "MExtractor.h"
73#include "MExtractTime.h"
74#include "MExtractTimeFastSpline.h"
75#include "MFCosmics.h"
76#include "MContinue.h"
77#include "MFillH.h"
78#include "MCalibrateData.h"
79#include "MCalibrateRelTimes.h"
80
81#include "MTriggerPattern.h"
82#include "MTriggerPatternDecode.h"
83#include "MFTriggerPattern.h"
84
85#include "MStatusDisplay.h"
86
87ClassImp(MJCalibTest);
88
89using namespace std;
90// --------------------------------------------------------------------------
91//
92// Default constructor.
93//
94// Sets fUseCosmicsFilter to kTRUE, fRuns to 0, fExtractor to NULL, fTimeExtractor to NULL
95// fDisplay to kNormalDisplay
96//
97MJCalibTest::MJCalibTest(const char *name, const char *title)
98 : fUseCosmicsFilter(kTRUE), fExtractor(NULL), fTimeExtractor(NULL),
99 fDisplayType(kNormalDisplay), fGeometry("MGeomCamMagic")
100{
101 fName = name ? name : "MJCalibTest";
102 fTitle = title ? title : "Tool to extract, calibrate and test signals from a file";
103}
104
105
106void MJCalibTest::DisplayResult(MParList &plist)
107{
108 if (!fDisplay)
109 return;
110
111 //
112 // Update display
113 //
114 TString title = fDisplay->GetTitle();
115 title += "-- Extraction-Calibration-Test ";
116 title += fRuns->GetRunsAsString();
117 title += " --";
118 fDisplay->SetTitle(title);
119
120 //
121 // Get container from list
122 //
123 MGeomCam &geomcam = *(MGeomCam*) plist.FindObject("MGeomCam");
124 MHCalibrationTestCam &testcam = *(MHCalibrationTestCam*)plist.FindObject("MHCalibrationTestCam");
125
126 // Create histograms to display
127 MHCamera disp1 (geomcam, "Test;Photons", "Mean of calibrated Photons");
128 MHCamera disp2 (geomcam, "Test;SigmaPhotons", "Sigma of calibrated photons");
129 MHCamera disp3 (geomcam, "Test;PhotonsPerArea", "Equiv. Cherenkov Photons per Area");
130 MHCamera disp4 (geomcam, "Test;SigmaPhotPerArea", "Sigma equiv. Cher. Photons per Area");
131 MHCamera disp5 (geomcam, "Test;Phot", "Calibrated Photons");
132 MHCamera disp6 (geomcam, "Test;PhotPerArea", "Calibrated Photons per Area");
133 MHCamera disp7 (geomcam, "Test;NotInterpolate", "Not interpolated pixels");
134 MHCamera disp8 (geomcam, "Test;DeviatingPhots", "Deviating Number Photons");
135
136 // Fitted charge means and sigmas
137 disp1.SetCamContent(testcam, 0);
138 disp1.SetCamError( testcam, 1);
139 disp2.SetCamContent(testcam, 2);
140 disp2.SetCamError( testcam, 3);
141 disp3.SetCamContent(testcam, 7);
142 disp3.SetCamError( testcam, 8);
143 disp4.SetCamContent(testcam, 9);
144 disp4.SetCamError( testcam, 10);
145
146 disp5.SetCamContent(fTestCam, 0);
147 disp5.SetCamError( fTestCam, 1);
148 disp6.SetCamContent(fTestCam, 2);
149 disp6.SetCamError( fTestCam, 3);
150 disp7.SetCamError( fTestCam, 4);
151
152 disp8.SetCamError( fBadPixels, 22);
153
154
155 disp1.SetYTitle("Photons");
156 disp2.SetYTitle("\\sigma_{phot}");
157 disp3.SetYTitle("Photons per Area [mm^{-2}]");
158 disp4.SetYTitle("\\sigma_{phot} per Area [mm^{-2}]");
159
160 disp5.SetYTitle("Photons");
161 disp6.SetYTitle("Photons per Area [mm^{-2}]");
162 disp7.SetYTitle("[1]");
163 disp8.SetYTitle("[1]");
164
165 gStyle->SetOptStat(1111);
166 gStyle->SetOptFit();
167
168 if (fDisplayType == kNormalDisplay)
169 {
170
171 TCanvas &c = fDisplay->AddTab("TestCharges");
172 c.Divide(4,4);
173
174 disp1.CamDraw(c, 1, 4, 2, 1);
175 disp2.CamDraw(c, 2, 4, 2, 1);
176 disp3.CamDraw(c, 3, 4, 1, 1);
177 disp4.CamDraw(c, 4, 4, 2, 1);
178 }
179
180 TCanvas &c2 = fDisplay->AddTab("TestResult");
181 c2.Divide(2,4);
182
183 disp5.CamDraw(c2, 1, 2, 2, 1);
184 disp6.CamDraw(c2, 2, 2, 2, 1);
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 return;
193
194}
195
196
197void MJCalibTest::DisplayResultT(MParList &plist)
198{
199 if (!fDisplay)
200 return;
201
202 //
203 // Update display
204 //
205 TString title = fDisplay->GetTitle();
206 title += "-- Extraction-Calibration-Test-Time";
207 title += fRuns->GetRunsAsString();
208 title += " --";
209 fDisplay->SetTitle(title);
210
211 //
212 // Get container from list
213 //
214 MGeomCam &geomcam = *(MGeomCam*)plist.FindObject("MGeomCam");
215
216 // Create histograms to display
217 MHCamera disp1 (geomcam, "Test;Arr.Times", "Mean of calibrated Arr.Times");
218 MHCamera disp2 (geomcam, "Test;SigmaArr.Times", "Sigma of calibrated Arr.Times");
219
220 // Fitted charge means and sigmas
221 disp1.SetCamContent(fTestTimeCam, 0);
222 disp1.SetCamError( fTestTimeCam, 1);
223 disp2.SetCamContent(fTestTimeCam, 2);
224 disp2.SetCamError( fTestTimeCam, 3);
225
226 disp1.SetYTitle("Mean Arr.Times [FADC units]");
227 disp2.SetYTitle("\\sigma_{t} [FADC units]");
228
229 gStyle->SetOptStat(1111);
230 gStyle->SetOptFit();
231
232 TCanvas &c = fDisplay->AddTab("TestTimes");
233 c.Divide(2,4);
234
235 disp1.CamDraw(c, 1, 2, 5, 1);
236 disp2.CamDraw(c, 2, 2, 5, 1);
237
238 return;
239
240}
241
242
243const char* MJCalibTest::GetInputFile() const
244{
245
246 if (fSequence.IsValid())
247 return Form("%s/calib%08d.root", (const char*)fPathOut, fSequence.GetSequence());
248
249 if (!fRuns)
250 return "";
251
252 return Form("%s/%s-F1.root", (const char*)fPathOut, (const char*)fRuns->GetRunsAsFileName());
253
254}
255
256
257const char* MJCalibTest::GetOutputFile() const
258{
259
260 if (fSequence.IsValid())
261 return Form("%s/test%08d.root", (const char*)fPathOut, fSequence.GetSequence());
262
263 if (!fRuns)
264 return "";
265
266 return Form("%s/%s-Test.root", (const char*)fPathOut, (const char*)fRuns->GetRunsAsFileName());
267
268}
269
270Bool_t MJCalibTest::ReadCalibration(TObjArray &l, MBadPixelsCam &cam, MExtractor* &ext1, MExtractor* &ext2, TString &geom) const
271{
272
273 const TString fname = GetInputFile();
274
275 *fLog << inf << "Reading from file: " << fname << endl;
276
277 TFile file(fname, "READ");
278 if (!file.IsOpen())
279 {
280 *fLog << err << dbginf << "ERROR - Could not open file " << fname << endl;
281 return kFALSE;
282 }
283
284 TObject *o = file.Get("ExtractSignal");
285 if (o && !o->InheritsFrom(MExtractor::Class()))
286 {
287 *fLog << err << dbginf << "ERROR - ExtractSignal read from " << fname << " doesn't inherit from MExtractor!" << endl;
288 return kFALSE;
289 }
290 ext1 = o ? (MExtractor*)o->Clone() : NULL;
291
292 o = file.Get("ExtractTime");
293 if (o && !o->InheritsFrom(MExtractor::Class()))
294 {
295 *fLog << err << dbginf << "ERROR - ExtractTime read from " << fname << " doesn't inherit from MExtractor!" << endl;
296 return kFALSE;
297 }
298 ext2 = o ? (MExtractor*)o->Clone() : NULL;
299 if (!ext1 && !ext2)
300 {
301 *fLog << err << dbginf << "ERROR - Neither ExtractSignal nor ExrtractTime found in " << fname << "!" << endl;
302 return kFALSE;
303 }
304
305 o = file.Get("MGeomCam");
306 if (o && !o->InheritsFrom(MGeomCam::Class()))
307 {
308 *fLog << err << dbginf << "ERROR - MGeomCam read from " << fname << " doesn't inherit from MGeomCam!" << endl;
309 return kFALSE;
310 }
311 geom = o ? o->ClassName() : "";
312
313 TObjArray cont(l);
314 cont.Add(&cam);
315 return ReadContainer(cont);
316}
317
318// --------------------------------------------------------------------------
319//
320// MJCalibration allows to setup several option by a resource file:
321// MJCalibrateSignal.RawData: yes,no
322//
323// For more details see the class description and the corresponding Getters
324//
325Bool_t MJCalibTest::CheckEnvLocal()
326{
327
328 SetUseRootData();
329
330 if (HasEnv("DataType"))
331 {
332 TString dat = GetEnv("DataType", "");
333 if (dat.BeginsWith("raw", TString::kIgnoreCase))
334 {
335 fDataFlag = 0;
336 SetUseRawData();
337 }
338 if (dat.BeginsWith("root", TString::kIgnoreCase))
339 {
340 fDataFlag = 0;
341 SetUseRootData();
342 }
343 if (dat.BeginsWith("mc", TString::kIgnoreCase))
344 {
345 fDataFlag = 0;
346 SetUseMC();
347 }
348 }
349 return kTRUE;
350}
351
352Bool_t MJCalibTest::Process(MPedestalCam &pedcam)
353{
354 // const TString fname = GetOutputFile();
355
356// if (gSystem->AccessPathName(fname, kFileExists))
357 return ProcessFile(pedcam);
358
359 // return kTRUE;
360}
361
362Bool_t MJCalibTest::ProcessFile(MPedestalCam &pedcam)
363{
364
365
366 if (!fSequence.IsValid())
367 {
368 if (!fRuns)
369 {
370 *fLog << err << "ERROR - Sequence invalid and no runs chosen!" << endl;
371 return kFALSE;
372 }
373
374 if (fRuns->GetNumRuns() != fRuns->GetNumEntries())
375 {
376 *fLog << err << "Number of files found doesn't match number of runs... abort."
377 << fRuns->GetNumRuns() << " vs. " << fRuns->GetNumEntries() << endl;
378 return kFALSE;
379 }
380 *fLog << "Calibrate data from ";
381 *fLog << "Runs " << fRuns->GetRunsAsString() << endl;
382 *fLog << endl;
383 }
384
385 CheckEnv();
386
387 *fLog << inf;
388 fLog->Separator(GetDescriptor());
389 *fLog << "Calculate MExtractedSignalCam from Runs " << fRuns->GetRunsAsString() << endl;
390 *fLog << endl;
391
392 MDirIter iter;
393
394 if (fSequence.IsValid())
395 {
396 const Int_t n0 = fSequence.SetupDatRuns(iter, fPathData, "D", IsUseRawData());
397 const Int_t n1 = fSequence.GetNumDatRuns();
398 if (n0==0)
399 {
400 *fLog << err << "ERROR - No input files of sequence found!" << endl;
401 return kFALSE;
402 }
403 if (n0!=n1)
404 {
405 *fLog << err << "ERROR - Number of files found (" << n0 << ") doesn't match number of files in sequence (" << n1 << ")" << endl;
406 return kFALSE;
407 }
408 }
409
410 MCalibrationChargeCam calcam;
411 MCalibrationQECam qecam;
412 MCalibrationRelTimeCam tmcam;
413 MBadPixelsCam badpix;
414
415 TObjArray calibcont;
416 calibcont.Add(&calcam);
417 calibcont.Add(&qecam);
418 calibcont.Add(&tmcam);
419
420 MExtractor *extractor1=0;
421 MExtractor *extractor2=0;
422 TString geom;
423
424 if (!ReadCalibration(calibcont, badpix, extractor1, extractor2, geom))
425 return kFALSE;
426
427 *fLog << all;
428 if (extractor1)
429 {
430 *fLog << underline << "Signal Extractor found in calibration file" << endl;
431 extractor1->Print();
432 *fLog << endl;
433 }
434 else
435 *fLog << inf << "No Signal Extractor: ExtractSignal in file." << endl;
436
437 if (extractor2)
438 {
439 *fLog << underline << "Time Extractor found in calibration file" << endl;
440 extractor2->Print();
441 *fLog << endl;
442 }
443 else
444 *fLog << inf << "No Time Extractor: ExtractTime in file." << endl;
445
446 if (!geom.IsNull())
447 *fLog << inf << "Camera geometry found in file: " << geom << endl;
448 else
449 *fLog << inf << "No Camera geometry found using default <MGeomCamMagic>" << endl;
450
451 if (fExtractor)
452 extractor1 = fExtractor;
453 if (fTimeExtractor)
454 extractor2 = fTimeExtractor;
455
456 // Setup Lists
457 MParList plist;
458 plist.AddToList(this); // take care of fDisplay!
459 plist.AddToList(&fTestCam);
460 plist.AddToList(&fTestTimeCam);
461 plist.AddToList(&badpix);
462 plist.AddToList(&pedcam);
463 plist.AddToList(&calcam);
464 plist.AddToList(&qecam);
465 plist.AddToList(&tmcam);
466
467 MCerPhotEvt cerphot;
468 MPedPhotCam pedphot;
469 MHCalibrationTestCam testcam;
470
471 plist.AddToList(&cerphot);
472 plist.AddToList(&pedphot);
473 plist.AddToList(&testcam);
474
475 pedcam.SetName("MPedestalFundamental");
476
477 MTaskList tlist;
478 plist.AddToList(&tlist);
479
480 // Setup Task-lists
481 MRawFileRead rawread(NULL);
482 MReadMarsFile read("Events");
483 read.DisableAutoScheme();
484
485 if (IsUseRawData())
486 rawread.AddFiles(fSequence.IsValid() ? iter : *fRuns);
487 else
488 static_cast<MRead&>(read).AddFiles(fSequence.IsValid() ? iter : *fRuns);
489
490 // Check for interleaved events
491 MTriggerPatternDecode decode;
492 MFTriggerPattern fcalib;
493 fcalib.DenyCalibration();
494 MContinue conttp(&fcalib, "ContTrigPattern");
495
496 MGeomApply apply; // Only necessary to craete geometry
497 if (!geom.IsNull())
498 apply.SetGeometry(geom);
499 MBadPixelsMerge merge(&badpix);
500
501 MExtractTimeAndChargeSlidingWindow extrsw;
502 MExtractTimeFastSpline extime;
503 extime.SetPedestals(&pedcam);
504
505 MTaskEnv taskenv1("ExtractSignal");
506 MTaskEnv taskenv2("ExtractTime");
507
508 if (extractor1)
509 {
510 extractor1->SetPedestals(&pedcam);
511 taskenv1.SetDefault(extractor1);
512 }
513 if (fTimeExtractor)
514 {
515 fTimeExtractor->SetPedestals(&pedcam);
516 taskenv2.SetDefault(fTimeExtractor);
517 }
518 else
519 {
520 extrsw.SetPedestals(&pedcam);
521 extrsw.SetWindowSize(8,8);
522 taskenv2.SetDefault(&extrsw);
523 *fLog << warn << GetDescriptor()
524 << ": No extractor has been chosen, take default MExtractTimeAndChargeSlidingWindow " << endl;
525 }
526
527 MCalibrateData photcalc;
528 MCalibrateRelTimes caltimes;
529 if (IsUseMC()) // MC file
530 {
531 photcalc.SetCalibrationMode(MCalibrateData::kFfactor);
532 photcalc.SetPedestalFlag(MCalibrateData::kRun);
533 photcalc.AddPedestal("MPedestalCam", "MPedPhotFundamental");
534 }
535 else
536 {
537 photcalc.SetCalibrationMode(MCalibrateData::kFfactor);
538 photcalc.AddPedestal("Fundamental");
539 photcalc.SetPedestalFlag(MCalibrateData::kEvent);
540 photcalc.SetSignalType(MCalibrateData::kPhot);
541 }
542
543 MBadPixelsCalc badcalc;
544 MBadPixelsTreat badtreat;
545 badtreat.SetProcessTimes(kFALSE);
546
547 badcalc.SetNamePedPhotCam("MPedPhotFundamental");
548 //badtreat.SetUseInterpolation();
549 badtreat.AddNamePedPhotCam("MPedPhotFundamental");
550
551 MCalibrationTestCalc testcalc;
552
553 if (!fSequence.IsValid())
554 {
555 testcalc.SetOutputPath(fPathOut);
556 testcalc.SetOutputFile(Form("%s-TestCalibStat.txt",(const char*)fRuns->GetRunsAsFileName()));
557 }
558
559 MHCamEvent evt0(0,"Signal", "Un-Calibrated Signal;;S [FADC cnts]" );
560 MHCamEvent evt1(3,"CalSig", "Calibrated and Interpolated Signal;;S [\\gamma]");
561 MHCamEvent evt2(0,"Times" , "Arrival Time;;T [slice]");
562
563 MFillH fill0(&evt0, "MExtractedSignalCam", "FillUncalibrated");
564 MFillH fill1(&evt1, "MCerPhotEvt", "FillCalibrated");
565 MFillH fill2(&evt2, "MArrivalTime","FillTimes");
566
567 MFillH fillcam("MHCalibrationTestCam", "MCerPhotEvt");
568 fillcam.SetNameTab("Test");
569 MFillH filltme("MHCalibrationTestTimeCam", "MArrivalTime");
570 filltme.SetNameTab("TestTime");
571
572 MFCosmics cosmics;
573 cosmics.SetNamePedestalCam("MPedestalFundamental");
574 MContinue contcos(&cosmics,"ContCosmics");
575
576 tlist.AddToList(&read);
577 tlist.AddToList(&decode);
578 tlist.AddToList(&apply);
579 tlist.AddToList(&merge);
580 tlist.AddToList(&conttp);
581 tlist.AddToList(&taskenv1);
582 if (!extractor1->InheritsFrom("MExtractTimeAndCharge"))
583 tlist.AddToList(&taskenv2);
584 tlist.AddToList(&contcos);
585 tlist.AddToList(&fill0);
586 tlist.AddToList(&photcalc);
587 tlist.AddToList(&caltimes);
588 tlist.AddToList(&badcalc);
589 tlist.AddToList(&badtreat);
590 tlist.AddToList(&fill1);
591 tlist.AddToList(&fill2);
592 tlist.AddToList(&fillcam);
593 tlist.AddToList(&filltme);
594 tlist.AddToList(&testcalc);
595
596 // Create and setup the eventloop
597 MEvtLoop evtloop(fName);
598 evtloop.SetParList(&plist);
599 evtloop.SetDisplay(fDisplay);
600 evtloop.SetLogStream(fLog);
601
602 // Execute first analysis
603 if (!evtloop.Eventloop())
604 {
605 *fLog << err << GetDescriptor() << ": Failed." << endl;
606 return kFALSE;
607 }
608
609 tlist.PrintStatistics();
610
611 DisplayResult(plist);
612
613 if (!WriteResult())
614 return kFALSE;
615
616 *fLog << inf << GetDescriptor() << ": Done." << endl;
617
618 return kTRUE;
619}
620
621Bool_t MJCalibTest::WriteResult()
622{
623
624 if (fPathOut.IsNull())
625 return kTRUE;
626
627 const TString oname(GetOutputFile());
628
629 *fLog << inf << "Writing to file: " << oname << endl;
630
631 TFile file(oname, "UPDATE");
632
633 if (fDisplay && fDisplay->Write()<=0)
634 {
635 *fLog << err << "Unable to write MStatusDisplay to " << oname << endl;
636 return kFALSE;
637 }
638
639 if (fPedPhotCam.Write()<=0)
640 {
641 *fLog << err << "Unable to write MPedPhotCam to " << oname << endl;
642 return kFALSE;
643 }
644
645 if (fTestCam.Write()<=0)
646 {
647 *fLog << err << "Unable to write MCalibrationTestCam to " << oname << endl;
648 return kFALSE;
649 }
650
651 if (fTestTimeCam.Write()<=0)
652 {
653 *fLog << err << "Unable to write MCalibrationTestCam to " << oname << endl;
654 return kFALSE;
655 }
656
657 return kTRUE;
658
659}
660
Note: See TracBrowser for help on using the repository browser.