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

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