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

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