source: trunk/MagicSoft/Mars/mjobs/MJExtractCalibTest.cc@ 5138

Last change on this file since 5138 was 4742, checked in by gaug, 20 years ago
*** empty log message ***
File size: 14.9 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-2004
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MJExtractCalibTest
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 "MJExtractCalibTest.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 "MEvtLoop.h"
48
49#include "MHCamera.h"
50
51#include "MPedestalCam.h"
52#include "MBadPixelsCam.h"
53#include "MBadPixelsTreat.h"
54#include "MCerPhotEvt.h"
55#include "MArrivalTime.h"
56#include "MCalibrationChargeCam.h"
57#include "MCalibrationRelTimeCam.h"
58#include "MCalibrationQECam.h"
59#include "MCalibrationTestCam.h"
60#include "MCalibrationTestCalc.h"
61#include "MHCamEvent.h"
62#include "MHCalibrationTestCam.h"
63
64#include "MReadMarsFile.h"
65#include "MGeomApply.h"
66#include "MExtractSlidingWindow.h"
67#include "MExtractor.h"
68#include "MExtractTime.h"
69#include "MExtractTimeFastSpline.h"
70#include "MFCosmics.h"
71#include "MContinue.h"
72#include "MFillH.h"
73#include "MCalibrateData.h"
74#include "MCalibrateRelTimes.h"
75#include "MPedPhotCalc.h"
76
77#include "MStatusDisplay.h"
78
79ClassImp(MJExtractCalibTest);
80
81using namespace std;
82// --------------------------------------------------------------------------
83//
84// Default constructor.
85//
86// Sets fUseCosmicsFilter to kTRUE, fRuns to 0, fExtractor to NULL, fTimeExtractor to NULL
87// fDisplay to kNormalDisplay
88//
89MJExtractCalibTest::MJExtractCalibTest(const char *name, const char *title)
90 : fUseCosmicsFilter(kTRUE), fRuns(NULL), fExtractor(NULL), fTimeExtractor(NULL),
91 fDisplayType(kNormalDisplay)
92{
93 fName = name ? name : "MJExtractCalibTest";
94 fTitle = title ? title : "Tool to extract, calibrate and test signals from a file";
95}
96
97
98void MJExtractCalibTest::DisplayResult(MParList &plist)
99{
100 if (!fDisplay)
101 return;
102
103 //
104 // Update display
105 //
106 TString title = fDisplay->GetTitle();
107 title += "-- Extraction-Calibration-Test ";
108 title += fRuns->GetRunsAsString();
109 title += " --";
110 fDisplay->SetTitle(title);
111
112 //
113 // Get container from list
114 //
115 MGeomCam &geomcam = *(MGeomCam*) plist.FindObject("MGeomCam");
116 MHCalibrationTestCam &testcam = *(MHCalibrationTestCam*)plist.FindObject("MHCalibrationTestCam");
117
118 // Create histograms to display
119 MHCamera disp1 (geomcam, "Test;Photons", "Mean of calibrated Photons");
120 MHCamera disp2 (geomcam, "Test;SigmaPhotons", "Sigma of calibrated photons");
121 MHCamera disp3 (geomcam, "Test;PhotonsPerArea", "Equiv. Cherenkov Photons per Area");
122 MHCamera disp4 (geomcam, "Test;SigmaPhotPerArea", "Sigma equiv. Cher. Photons per Area");
123 MHCamera disp5 (geomcam, "Test;Phot", "Calibrated Photons");
124 MHCamera disp6 (geomcam, "Test;PhotPerArea", "Calibrated Photons per Area");
125 MHCamera disp7 (geomcam, "Test;NotInterpolate", "Not interpolated pixels");
126 MHCamera disp8 (geomcam, "Test;DeviatingPhots", "Deviating Number Photons");
127
128 // Fitted charge means and sigmas
129 disp1.SetCamContent(testcam, 0);
130 disp1.SetCamError( testcam, 1);
131 disp2.SetCamContent(testcam, 2);
132 disp2.SetCamError( testcam, 3);
133 disp3.SetCamContent(testcam, 7);
134 disp3.SetCamError( testcam, 8);
135 disp4.SetCamContent(testcam, 9);
136 disp4.SetCamError( testcam, 10);
137
138 disp5.SetCamContent(fTestCam, 0);
139 disp5.SetCamError( fTestCam, 1);
140 disp6.SetCamContent(fTestCam, 2);
141 disp6.SetCamError( fTestCam, 3);
142 disp7.SetCamError( fTestCam, 4);
143
144 disp8.SetCamError( fBadPixels, 22);
145
146
147 disp1.SetYTitle("Photons");
148 disp2.SetYTitle("\\sigma_{phot}");
149 disp3.SetYTitle("Photons per Area [mm^{-2}]");
150 disp4.SetYTitle("\\sigma_{phot} per Area [mm^{-2}]");
151
152 disp5.SetYTitle("Photons");
153 disp6.SetYTitle("Photons per Area [mm^{-2}]");
154 disp7.SetYTitle("[1]");
155 disp8.SetYTitle("[1]");
156
157 gStyle->SetOptStat(1111);
158 gStyle->SetOptFit();
159
160 if (fDisplayType == kNormalDisplay)
161 {
162
163 TCanvas &c = fDisplay->AddTab("TestCharges");
164 c.Divide(4,4);
165
166 disp1.CamDraw(c, 1, 4, 2, 1);
167 disp2.CamDraw(c, 2, 4, 2, 1);
168 disp3.CamDraw(c, 3, 4, 1, 1);
169 disp4.CamDraw(c, 4, 4, 2, 1);
170 }
171
172 TCanvas &c2 = fDisplay->AddTab("TestResult");
173 c2.Divide(2,4);
174
175 disp5.CamDraw(c2, 1, 2, 2, 1);
176 disp6.CamDraw(c2, 2, 2, 2, 1);
177
178 TCanvas &c3 = fDisplay->AddTab("TestDefects");
179 c3.Divide(2,2);
180
181 disp7.CamDraw(c3, 1, 2, 0);
182 disp8.CamDraw(c3, 2, 2, 0);
183
184 return;
185
186}
187
188
189void MJExtractCalibTest::DisplayResultT(MParList &plist)
190{
191 if (!fDisplay)
192 return;
193
194 //
195 // Update display
196 //
197 TString title = fDisplay->GetTitle();
198 title += "-- Extraction-Calibration-Test-Time";
199 title += fRuns->GetRunsAsString();
200 title += " --";
201 fDisplay->SetTitle(title);
202
203 //
204 // Get container from list
205 //
206 MGeomCam &geomcam = *(MGeomCam*)plist.FindObject("MGeomCam");
207
208 // Create histograms to display
209 MHCamera disp1 (geomcam, "Test;Arr.Times", "Mean of calibrated Arr.Times");
210 MHCamera disp2 (geomcam, "Test;SigmaArr.Times", "Sigma of calibrated Arr.Times");
211
212 // Fitted charge means and sigmas
213 disp1.SetCamContent(fTestTimeCam, 0);
214 disp1.SetCamError( fTestTimeCam, 1);
215 disp2.SetCamContent(fTestTimeCam, 2);
216 disp2.SetCamError( fTestTimeCam, 3);
217
218 disp1.SetYTitle("Mean Arr.Times [FADC units]");
219 disp2.SetYTitle("\\sigma_{t} [FADC units]");
220
221 gStyle->SetOptStat(1111);
222 gStyle->SetOptFit();
223
224 TCanvas &c = fDisplay->AddTab("TestTimes");
225 c.Divide(2,4);
226
227 disp1.CamDraw(c, 1, 2, 5, 1);
228 disp2.CamDraw(c, 2, 2, 5, 1);
229
230 return;
231
232}
233
234
235const char* MJExtractCalibTest::GetOutputFile() const
236{
237
238 if (fSequence.IsValid())
239 return Form("%s/test%06d.root", (const char*)fPathOut, fSequence.GetSequence());
240
241 if (!fRuns)
242 return "";
243
244 return Form("%s/%s-Test.root", (const char*)fPathOut, (const char*)fRuns->GetRunsAsFileName());
245
246}
247
248
249Bool_t MJExtractCalibTest::ProcessD(MPedestalCam &pedcam, MCalibrationChargeCam &calcam, MCalibrationQECam &qecam)
250{
251 // const TString fname = GetOutputFile();
252
253// if (gSystem->AccessPathName(fname, kFileExists))
254 return ProcessFileD(pedcam,calcam,qecam);
255
256 // return kTRUE;
257}
258
259Bool_t MJExtractCalibTest::ProcessFileD(MPedestalCam &pedcam, MCalibrationChargeCam &calcam, MCalibrationQECam &qecam)
260{
261 if (!fRuns)
262 {
263 *fLog << err << "No Runs choosen... abort." << endl;
264 return kFALSE;
265 }
266 if (fRuns->GetNumRuns() != fRuns->GetNumEntries())
267 {
268 *fLog << err << "Number of files found doesn't match number of runs... abort." << endl;
269 return kFALSE;
270 }
271
272 *fLog << inf;
273 fLog->Separator(GetDescriptor());
274 *fLog << "Calculate MExtractedSignalCam from Runs " << fRuns->GetRunsAsString() << endl;
275 *fLog << endl;
276
277 MCerPhotEvt cerphot;
278 MPedPhotCam pedphot;
279 MHCalibrationTestCam testcam;
280
281 // Setup Lists
282 MParList plist;
283 plist.AddToList(&pedcam);
284 plist.AddToList(&calcam);
285 plist.AddToList(&qecam);
286 plist.AddToList(&cerphot);
287 plist.AddToList(&pedphot);
288 plist.AddToList(&testcam);
289 plist.AddToList(&fTestCam);
290 plist.AddToList(&fBadPixels);
291
292 MTaskList tlist;
293 plist.AddToList(&tlist);
294
295 // Setup Task-lists
296 MReadMarsFile read("Events");
297 read.DisableAutoScheme();
298 static_cast<MRead&>(read).AddFiles(*fRuns);
299
300 MGeomApply apply; // Only necessary to craete geometry
301 MExtractSlidingWindow extract2;
302 MCalibrateData photcalc;
303 photcalc.SetCalibrationMode(MCalibrateData::kFfactor);
304 MPedPhotCalc pedphotcalc;
305 MBadPixelsTreat badtreat;
306 badtreat.SetUseInterpolation();
307 MCalibrationTestCalc testcalc;
308
309 if (!fSequence.IsValid())
310 {
311 testcalc.SetOutputPath(fPathOut);
312 testcalc.SetOutputFile(Form("%s-TestCalibStat.txt",(const char*)fRuns->GetRunsAsFileName()));
313 }
314
315 MHCamEvent evt("ExtSignal");
316 evt.SetType(0);
317 MFillH fill(&evt, "MExtractedSignalCam");
318
319 MFillH fillcam("MHCalibrationTestCam", "MCerPhotEvt");
320 fillcam.SetNameTab("Test");
321
322 MFCosmics cosmics;
323 MContinue cont(&cosmics);
324
325 tlist.AddToList(&read);
326 tlist.AddToList(&apply);
327
328 if (fExtractor)
329 tlist.AddToList(fExtractor);
330 else
331 {
332 *fLog << warn << GetDescriptor()
333 << ": No extractor has been chosen, take default MExtractSlidingWindow " << endl;
334 tlist.AddToList(&extract2);
335 }
336
337
338 if (fUseCosmicsFilter)
339 tlist.AddToList(&cont);
340
341 tlist.AddToList(&fill);
342 tlist.AddToList(&photcalc);
343 tlist.AddToList(&pedphotcalc);
344 tlist.AddToList(&badtreat);
345 tlist.AddToList(&fillcam);
346 tlist.AddToList(&testcalc);
347
348 // Create and setup the eventloop
349 MEvtLoop evtloop(fName);
350 evtloop.SetParList(&plist);
351 evtloop.SetDisplay(fDisplay);
352 evtloop.SetLogStream(fLog);
353
354 // Execute first analysis
355 if (!evtloop.Eventloop())
356 {
357 *fLog << err << GetDescriptor() << ": Failed." << endl;
358 return kFALSE;
359 }
360
361 tlist.PrintStatistics();
362
363 DisplayResult(plist);
364
365 if (!WriteResultD())
366 return kFALSE;
367
368 *fLog << inf << GetDescriptor() << ": Done." << endl;
369
370 return kTRUE;
371}
372
373Bool_t MJExtractCalibTest::ProcessT(MPedestalCam &pedcam, MCalibrationRelTimeCam &relcam)
374{
375
376// const TString fname = GetOutputFile();
377
378// if (gSystem->AccessPathName(fname, kFileExists))
379 return ProcessFileT(pedcam,relcam);
380
381// return kTRUE;
382}
383
384Bool_t MJExtractCalibTest::ProcessFileT(MPedestalCam &pedcam, MCalibrationRelTimeCam &relcam)
385{
386
387 if (!fRuns)
388 {
389 *fLog << err << "No Runs choosen... abort." << endl;
390 return kFALSE;
391 }
392 if (fRuns->GetNumRuns() != fRuns->GetNumEntries())
393 {
394 *fLog << err << "Number of files found doesn't match number of runs... abort." << endl;
395 return kFALSE;
396 }
397
398 *fLog << inf;
399 fLog->Separator(GetDescriptor());
400 *fLog << "Calculate MExtractedSignalCam from Runs " << fRuns->GetRunsAsString() << endl;
401 *fLog << endl;
402
403 MArrivalTime arrtime;
404
405 // Setup Lists
406 MParList plist;
407 plist.AddToList(&pedcam);
408 plist.AddToList(&relcam);
409 plist.AddToList(&arrtime);
410 plist.AddToList(&fTestTimeCam);
411 plist.AddToList(&fBadPixels);
412
413 MTaskList tlist;
414 plist.AddToList(&tlist);
415
416 // Setup Task-lists
417 MReadMarsFile read("Events");
418 read.DisableAutoScheme();
419 static_cast<MRead&>(read).AddFiles(*fRuns);
420
421 MGeomApply apply; // Only necessary to craete geometry
422 MExtractTimeFastSpline extract;
423 MExtractSlidingWindow extcharge; // Only for the cosmics filter
424 MCalibrateRelTimes timecalc;
425 MFCosmics cosmics;
426 MContinue cont(&cosmics);
427
428 MHCamEvent evt("ExtTimes");
429 evt.SetType(0);
430 MFillH fill(&evt, "MArrivalTimeCam");
431
432 MFillH fillcam("MHCalibrationTestTimeCam", "MArrivalTime");
433 fillcam.SetNameTab("TestTime");
434
435 tlist.AddToList(&read);
436 tlist.AddToList(&apply);
437
438 if (fTimeExtractor)
439 tlist.AddToList(fTimeExtractor);
440 else
441 {
442 *fLog << warn << GetDescriptor()
443 << ": No extractor has been chosen, take default MExtractTimeFastSpline " << endl;
444 tlist.AddToList(&extract);
445 }
446
447 tlist.AddToList(&extcharge);
448 tlist.AddToList(&cont);
449 tlist.AddToList(&fill);
450 tlist.AddToList(&timecalc);
451 tlist.AddToList(&fillcam);
452
453 // Create and setup the eventloop
454 MEvtLoop evtloop(fName);
455 evtloop.SetParList(&plist);
456 evtloop.SetDisplay(fDisplay);
457 evtloop.SetLogStream(fLog);
458
459 // Execute first analysis
460 if (!evtloop.Eventloop())
461 {
462 *fLog << err << GetDescriptor() << ": Failed." << endl;
463 return kFALSE;
464 }
465
466 tlist.PrintStatistics();
467
468 DisplayResultT(plist);
469
470 if (!WriteResultT())
471 return kFALSE;
472
473 *fLog << inf << GetDescriptor() << ": Done." << endl;
474
475 return kTRUE;
476}
477
478
479Bool_t MJExtractCalibTest::ReadPedPhotCam()
480{
481
482 const TString fname = GetOutputFile();
483
484 if (gSystem->AccessPathName(fname, kFileExists))
485 {
486 *fLog << err << "Input file " << fname << " doesn't exist." << endl;
487 return kFALSE;
488 }
489
490 *fLog << inf << "Reading from file: " << fname << endl;
491
492 TFile file(fname, "READ");
493 if (fPedPhotCam.Read()<=0)
494 {
495 *fLog << "Unable to read MPedPhotCam from " << fname << endl;
496 return kFALSE;
497 }
498
499 if (file.FindKey("MBadPixelsCam"))
500 {
501 MBadPixelsCam bad;
502 if (bad.Read()<=0)
503 {
504 *fLog << "Unable to read MBadPixelsCam from " << fname << endl;
505 return kFALSE;
506 }
507 fBadPixels.Merge(bad);
508 }
509
510 if (fDisplay /*&& !fDisplay->GetCanvas("Pedestals")*/) // FIXME!
511 fDisplay->Read();
512
513 return kTRUE;
514}
515
516Bool_t MJExtractCalibTest::WriteResultD()
517{
518
519 if (fPathOut.IsNull())
520 return kTRUE;
521
522 const TString oname(GetOutputFile());
523
524 *fLog << inf << "Writing to file: " << oname << endl;
525
526 TFile file(oname, "UPDATE");
527
528 if (fDisplay && fDisplay->Write()<=0)
529 {
530 *fLog << err << "Unable to write MStatusDisplay to " << oname << endl;
531 return kFALSE;
532 }
533
534 if (fPedPhotCam.Write()<=0)
535 {
536 *fLog << err << "Unable to write MPedPhotCam to " << oname << endl;
537 return kFALSE;
538 }
539
540 if (fTestCam.Write()<=0)
541 {
542 *fLog << err << "Unable to write MCalibrationTestCam to " << oname << endl;
543 return kFALSE;
544 }
545
546 return kTRUE;
547
548}
549
550Bool_t MJExtractCalibTest::WriteResultT()
551{
552
553 if (fPathOut.IsNull())
554 return kTRUE;
555
556 const TString oname(GetOutputFile());
557
558 *fLog << inf << "Writing to file: " << oname << endl;
559
560 TFile file(oname, "UPDATE");
561
562 if (fDisplay && fDisplay->Write()<=0)
563 {
564 *fLog << err << "Unable to write MStatusDisplay to " << oname << endl;
565 return kFALSE;
566 }
567
568 if (fTestTimeCam.Write()<=0)
569 {
570 *fLog << err << "Unable to write MHCalibrationTestTimeCam to " << oname << endl;
571 return kFALSE;
572 }
573
574 return kTRUE;
575
576}
577
578
579Bool_t MJExtractCalibTest::CheckEnv()
580{
581 if (HasEnv("DataCheckDisplay"))
582 fDisplayType = GetEnv("DataCheckDisplay", kFALSE) ? kDataCheckDisplay : kNormalDisplay;
583
584 SetOverwrite(GetEnv("Overwrite", fOverwrite));
585
586 return MJob::CheckEnv();
587}
Note: See TracBrowser for help on using the repository browser.