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

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