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

Last change on this file since 4698 was 4658, checked in by gaug, 20 years ago
*** empty log message ***
File size: 14.6 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 "MCalibrate.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
235void MJExtractCalibTest::SetOutputPath(const char *path)
236{
237 fOutputPath = path;
238 if (fOutputPath.EndsWith("/"))
239 fOutputPath = fOutputPath(0, fOutputPath.Length()-1);
240}
241
242const char* MJExtractCalibTest::GetOutputFile() const
243{
244 if (!fRuns)
245 return "";
246
247 return Form("%s/%s-Test.root", (const char*)fOutputPath, (const char*)fRuns->GetRunsAsFileName());
248}
249
250
251Bool_t MJExtractCalibTest::ProcessD(MPedestalCam &pedcam, MCalibrationChargeCam &calcam, MCalibrationQECam &qecam)
252{
253 // const TString fname = GetOutputFile();
254
255// if (gSystem->AccessPathName(fname, kFileExists))
256 return ProcessFileD(pedcam,calcam,qecam);
257
258 // return kTRUE;
259}
260
261Bool_t MJExtractCalibTest::ProcessFileD(MPedestalCam &pedcam, MCalibrationChargeCam &calcam, MCalibrationQECam &qecam)
262{
263 if (!fRuns)
264 {
265 *fLog << err << "No Runs choosen... abort." << endl;
266 return kFALSE;
267 }
268 if (fRuns->GetNumRuns() != fRuns->GetNumEntries())
269 {
270 *fLog << err << "Number of files found doesn't match number of runs... abort." << endl;
271 return kFALSE;
272 }
273
274 *fLog << inf;
275 fLog->Separator(GetDescriptor());
276 *fLog << "Calculate MExtractedSignalCam from Runs " << fRuns->GetRunsAsString() << endl;
277 *fLog << endl;
278
279 MCerPhotEvt cerphot;
280 MPedPhotCam pedphot;
281 MHCalibrationTestCam testcam;
282
283 // Setup Lists
284 MParList plist;
285 plist.AddToList(&pedcam);
286 plist.AddToList(&calcam);
287 plist.AddToList(&qecam);
288 plist.AddToList(&cerphot);
289 plist.AddToList(&pedphot);
290 plist.AddToList(&testcam);
291 plist.AddToList(&fTestCam);
292 plist.AddToList(&fBadPixels);
293
294 MTaskList tlist;
295 plist.AddToList(&tlist);
296
297 // Setup Task-lists
298 MReadMarsFile read("Events");
299 read.DisableAutoScheme();
300 static_cast<MRead&>(read).AddFiles(*fRuns);
301
302 MGeomApply apply; // Only necessary to craete geometry
303 MExtractSlidingWindow extract2;
304 MCalibrate photcalc;
305 photcalc.SetCalibrationMode(MCalibrate::kFfactor);
306 MPedPhotCalc pedphotcalc;
307 MBadPixelsTreat badtreat;
308 badtreat.SetUseInterpolation();
309 MCalibrationTestCalc testcalc;
310 testcalc.SetOutputPath(fOutputPath);
311 testcalc.SetOutputFile(Form("%s-TestCalibStat.txt",(const char*)fRuns->GetRunsAsFileName()));
312
313 MHCamEvent evt("ExtSignal");
314 evt.SetType(0);
315 MFillH fill(&evt, "MExtractedSignalCam");
316
317 MFillH fillcam("MHCalibrationTestCam", "MCerPhotEvt");
318 fillcam.SetNameTab("Test");
319
320 MFCosmics cosmics;
321 MContinue cont(&cosmics);
322
323 tlist.AddToList(&read);
324 tlist.AddToList(&apply);
325
326 if (fExtractor)
327 tlist.AddToList(fExtractor);
328 else
329 {
330 *fLog << warn << GetDescriptor()
331 << ": No extractor has been chosen, take default MExtractSlidingWindow " << endl;
332 tlist.AddToList(&extract2);
333 }
334
335
336 if (fUseCosmicsFilter)
337 tlist.AddToList(&cont);
338
339 tlist.AddToList(&fill);
340 tlist.AddToList(&photcalc);
341 tlist.AddToList(&pedphotcalc);
342 tlist.AddToList(&badtreat);
343 tlist.AddToList(&fillcam);
344 tlist.AddToList(&testcalc);
345
346 // Create and setup the eventloop
347 MEvtLoop evtloop(fName);
348 evtloop.SetParList(&plist);
349 evtloop.SetDisplay(fDisplay);
350 evtloop.SetLogStream(fLog);
351
352 // Execute first analysis
353 if (!evtloop.Eventloop())
354 {
355 *fLog << err << GetDescriptor() << ": Failed." << endl;
356 return kFALSE;
357 }
358
359 tlist.PrintStatistics();
360
361 DisplayResult(plist);
362
363 if (!WriteResultD())
364 return kFALSE;
365
366 *fLog << inf << GetDescriptor() << ": Done." << endl;
367
368 return kTRUE;
369}
370
371Bool_t MJExtractCalibTest::ProcessT(MPedestalCam &pedcam, MCalibrationRelTimeCam &relcam)
372{
373
374// const TString fname = GetOutputFile();
375
376// if (gSystem->AccessPathName(fname, kFileExists))
377 return ProcessFileT(pedcam,relcam);
378
379// return kTRUE;
380}
381
382Bool_t MJExtractCalibTest::ProcessFileT(MPedestalCam &pedcam, MCalibrationRelTimeCam &relcam)
383{
384
385 if (!fRuns)
386 {
387 *fLog << err << "No Runs choosen... abort." << endl;
388 return kFALSE;
389 }
390 if (fRuns->GetNumRuns() != fRuns->GetNumEntries())
391 {
392 *fLog << err << "Number of files found doesn't match number of runs... abort." << endl;
393 return kFALSE;
394 }
395
396 *fLog << inf;
397 fLog->Separator(GetDescriptor());
398 *fLog << "Calculate MExtractedSignalCam from Runs " << fRuns->GetRunsAsString() << endl;
399 *fLog << endl;
400
401 MArrivalTime arrtime;
402
403 // Setup Lists
404 MParList plist;
405 plist.AddToList(&pedcam);
406 plist.AddToList(&relcam);
407 plist.AddToList(&arrtime);
408 plist.AddToList(&fTestTimeCam);
409 plist.AddToList(&fBadPixels);
410
411 MTaskList tlist;
412 plist.AddToList(&tlist);
413
414 // Setup Task-lists
415 MReadMarsFile read("Events");
416 read.DisableAutoScheme();
417 static_cast<MRead&>(read).AddFiles(*fRuns);
418
419 MGeomApply apply; // Only necessary to craete geometry
420 MExtractTimeFastSpline extract;
421 MExtractSlidingWindow extcharge; // Only for the cosmics filter
422 MCalibrateRelTimes timecalc;
423 MFCosmics cosmics;
424 MContinue cont(&cosmics);
425
426 MHCamEvent evt("ExtTimes");
427 evt.SetType(0);
428 MFillH fill(&evt, "MArrivalTimeCam");
429
430 MFillH fillcam("MHCalibrationTestTimeCam", "MArrivalTime");
431 fillcam.SetNameTab("TestTime");
432
433 tlist.AddToList(&read);
434 tlist.AddToList(&apply);
435
436 if (fTimeExtractor)
437 tlist.AddToList(fTimeExtractor);
438 else
439 {
440 *fLog << warn << GetDescriptor()
441 << ": No extractor has been chosen, take default MExtractTimeFastSpline " << endl;
442 tlist.AddToList(&extract);
443 }
444
445 tlist.AddToList(&extcharge);
446 tlist.AddToList(&cont);
447 tlist.AddToList(&fill);
448 tlist.AddToList(&timecalc);
449 tlist.AddToList(&fillcam);
450
451 // Create and setup the eventloop
452 MEvtLoop evtloop(fName);
453 evtloop.SetParList(&plist);
454 evtloop.SetDisplay(fDisplay);
455 evtloop.SetLogStream(fLog);
456
457 // Execute first analysis
458 if (!evtloop.Eventloop())
459 {
460 *fLog << err << GetDescriptor() << ": Failed." << endl;
461 return kFALSE;
462 }
463
464 tlist.PrintStatistics();
465
466 DisplayResultT(plist);
467
468 if (!WriteResultT())
469 return kFALSE;
470
471 *fLog << inf << GetDescriptor() << ": Done." << endl;
472
473 return kTRUE;
474}
475
476
477Bool_t MJExtractCalibTest::ReadPedPhotCam()
478{
479
480 const TString fname = GetOutputFile();
481
482 if (gSystem->AccessPathName(fname, kFileExists))
483 {
484 *fLog << err << "Input file " << fname << " doesn't exist." << endl;
485 return kFALSE;
486 }
487
488 *fLog << inf << "Reading from file: " << fname << endl;
489
490 TFile file(fname, "READ");
491 if (fPedPhotCam.Read()<=0)
492 {
493 *fLog << "Unable to read MPedPhotCam from " << fname << endl;
494 return kFALSE;
495 }
496
497 if (file.FindKey("MBadPixelsCam"))
498 {
499 MBadPixelsCam bad;
500 if (bad.Read()<=0)
501 {
502 *fLog << "Unable to read MBadPixelsCam from " << fname << endl;
503 return kFALSE;
504 }
505 fBadPixels.Merge(bad);
506 }
507
508 if (fDisplay /*&& !fDisplay->GetCanvas("Pedestals")*/) // FIXME!
509 fDisplay->Read();
510
511 return kTRUE;
512}
513
514Bool_t MJExtractCalibTest::WriteResultD()
515{
516
517 if (fOutputPath.IsNull())
518 return kTRUE;
519
520 const TString oname(GetOutputFile());
521
522 *fLog << inf << "Writing to file: " << oname << endl;
523
524 TFile file(oname, "UPDATE");
525
526 if (fDisplay && fDisplay->Write()<=0)
527 {
528 *fLog << err << "Unable to write MStatusDisplay to " << oname << endl;
529 return kFALSE;
530 }
531
532 if (fPedPhotCam.Write()<=0)
533 {
534 *fLog << err << "Unable to write MPedPhotCam to " << oname << endl;
535 return kFALSE;
536 }
537
538 if (fTestCam.Write()<=0)
539 {
540 *fLog << err << "Unable to write MCalibrationTestCam to " << oname << endl;
541 return kFALSE;
542 }
543
544 return kTRUE;
545
546}
547
548Bool_t MJExtractCalibTest::WriteResultT()
549{
550
551 if (fOutputPath.IsNull())
552 return kTRUE;
553
554 const TString oname(GetOutputFile());
555
556 *fLog << inf << "Writing to file: " << oname << endl;
557
558 TFile file(oname, "UPDATE");
559
560 if (fDisplay && fDisplay->Write()<=0)
561 {
562 *fLog << err << "Unable to write MStatusDisplay to " << oname << endl;
563 return kFALSE;
564 }
565
566 if (fTestTimeCam.Write()<=0)
567 {
568 *fLog << err << "Unable to write MHCalibrationTestTimeCam to " << oname << endl;
569 return kFALSE;
570 }
571
572 return kTRUE;
573
574}
575
576
Note: See TracBrowser for help on using the repository browser.