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

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