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

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