source: trunk/MagicSoft/Mars/mjobs/MJExtractSignal.cc@ 3492

Last change on this file since 3492 was 3387, checked in by tbretz, 21 years ago
*** empty log message ***
File size: 17.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): Thomas Bretz, 1/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2004
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MJExtractSignal
28//
29/////////////////////////////////////////////////////////////////////////////
30#include "MJExtractSignal.h"
31
32#include <TFile.h>
33#include <TStyle.h>
34#include <TCanvas.h>
35#include <TSystem.h>
36
37#include "MLog.h"
38#include "MLogManip.h"
39
40#include "MRunIter.h"
41#include "MParList.h"
42#include "MTaskList.h"
43#include "MEvtLoop.h"
44
45#include "MHCamera.h"
46
47#include "MPedestalCam.h"
48#include "MCalibrationChargeCam.h"
49#include "MHCamEvent.h"
50
51#include "MReadMarsFile.h"
52#include "MGeomApply.h"
53#include "MBadPixelsMerge.h"
54#include "MExtractSignal.h"
55#include "MFillH.h"
56#include "MCalibrate.h"
57#include "MPedPhotCalc.h"
58#include "MWriteRootFile.h"
59
60#include "MJExtractSignal.h"
61#include "MStatusDisplay.h"
62
63
64ClassImp(MJExtractSignal);
65
66using namespace std;
67
68MJExtractSignal::MJExtractSignal(const char *name, const char *title) : fRuns(0)
69{
70 fName = name ? name : "MJExtractSignal";
71 fTitle = title ? title : "Tool to create a pedestal file (MPedestalCam)";
72}
73/*
74void MJExtractSignal::CamDraw(TCanvas &c, MHCamera &cam1, MCamEvent &note, Int_t i)
75{
76 c.cd(i);
77 gPad->SetBorderMode(0);
78 MHCamera *obj1=(MHCamera*)cam1.DrawCopy("hist");
79 obj1->AddNotify(note);
80
81 c.cd(i+2);
82 gPad->SetBorderMode(0);
83 obj1->Draw();
84
85 c.cd(i+4);
86 gPad->SetBorderMode(0);
87 TH1D *obj2 = (TH1D*)obj1->Projection();
88 obj2->Draw();
89 obj2->SetBit(kCanDelete);
90 obj2->Fit("gaus","Q");
91 obj2->GetFunction("gaus")->SetLineColor(kYellow);
92}
93
94void MJExtractSignal::CamDraw(TCanvas &c, MHCamera &cam1, MHCamera &cam2, MCamEvent &note)
95{
96 c.Divide(2, 3);
97
98 CamDraw(c, cam1, note, 1);
99 CamDraw(c, cam2, note, 2);
100}
101
102void MJExtractSignal::DisplayResult(MParList &plist)
103{
104 if (!fDisplay)
105 return;
106
107 //
108 // Update display
109 //
110 TString title = fDisplay->GetTitle();
111 title += "-- Calibration";
112 title += fRunNumber<0 ? Form("File %s", (const char*)fFileName) : Form("Run #%d", fRunNumber);
113 title += " --";
114 fDisplay->SetTitle(title);
115
116 //
117 // Get container from list
118 //
119 MGeomCam &geomcam = *(MGeomCam*)plist.FindObject("MGeomCam");
120
121 // Create histograms to display
122 MHCamera disp1 (geomcam, "Cal;Charge", "Fitted Mean Charges");
123 MHCamera disp3 (geomcam, "Cal;SigmaCharge", "Sigma of Fitted Charges");
124 MHCamera disp5 (geomcam, "Cal;ChargeProb", "Probability of Fit");
125 MHCamera disp6 (geomcam, "Cal;Time", "Arrival Times");
126 MHCamera disp7 (geomcam, "Cal;SigmaTime", "Sigma of Arrival Times");
127 MHCamera disp8 (geomcam, "Cal;TimeChiSquare", "Chi Square of Time Fit");
128// MHCamera disp9 (geomcam, "Cal;Ped", "Pedestals");
129// MHCamera disp10(geomcam, "Cal;PedRms", "Pedestal RMS");
130 MHCamera disp11(geomcam, "Cal;RSigma", "Reduced Sigmas");
131 MHCamera disp12(geomcam, "Cal;PheFFactorMethod", "Nr. of Phe's (F-Factor Method)");
132 MHCamera disp13(geomcam, "Cal;MeanConversionFFactorMethod", "Conversion Factor (F-Factor Method)");
133 MHCamera disp14(geomcam, "Cal;MeanPhotInsidePlexiglass", "Nr. of Photons (Blind Pixel Method)");
134 MHCamera disp15(geomcam, "Cal;MeanConversionBlindPixelMethod", "Conversion Factor (Blind Pixel Method)");
135 MHCamera disp16(geomcam, "Cal;RSigma/Charge", "Reduced Sigma per Charge");
136
137 disp1.SetCamContent(fCalibrationCam, 0);
138 disp1.SetCamError(fCalibrationCam, 1);
139
140 disp3.SetCamContent(fCalibrationCam, 2);
141 disp3.SetCamError(fCalibrationCam, 3);
142
143 disp5.SetCamContent(fCalibrationCam, 4);
144
145 disp6.SetCamContent(fCalibrationCam, 5);
146 disp6.SetCamError(fCalibrationCam, 6);
147 disp7.SetCamContent(fCalibrationCam, 6);
148 disp8.SetCamContent(fCalibrationCam, 7);
149 disp11.SetCamContent(fCalibrationCam, 10);
150
151 disp12.SetCamContent(fCalibrationCam, 11);
152 disp12.SetCamError(fCalibrationCam, 12);
153
154 disp13.SetCamContent(fCalibrationCam, 13);
155 disp13.SetCamError(fCalibrationCam, 14);
156
157 disp14.SetCamContent(fCalibrationCam, 15);
158 disp15.SetCamContent(fCalibrationCam, 16);
159 disp16.SetCamContent(fCalibrationCam, 17);
160
161 disp1.SetYTitle("Charge [FADC counts]");
162 disp3.SetYTitle("\\sigma_{Charge} [FADC counts]");
163 disp5.SetYTitle("P_{Charge} [1]");
164 disp6.SetYTitle("Arr. Time [Time Slice Nr.]");
165 disp7.SetYTitle("\\sigma_{Time} [Time Slices]");
166 disp8.SetYTitle("\\chi^{2}_{Time} [1]");
167// disp9.SetYTitle("Ped [FADC Counts ]");
168// disp10.SetYTitle("RMS_{Ped} [FADC Counts ]");
169 disp11.SetYTitle("\\sqrt{\\sigma^{2}_{Charge} - RMS^{2}_{Ped}} [FADC Counts]");
170 disp12.SetYTitle("Nr. Photo-Electrons [1]");
171 disp13.SetYTitle("Conversion Factor [PhE/FADC Count]");
172 disp14.SetYTitle("Nr. Photons [1]");
173 disp15.SetYTitle("Conversion Factor [Phot/FADC Count]");
174 disp16.SetYTitle("Reduced Sigma / Charge [1]");
175
176 gStyle->SetOptStat(1111);
177 gStyle->SetOptFit();
178
179 // Display Histograms: Charges
180 CamDraw(fDisplay->AddTab("FitdCharge"), disp1, disp3, fCalibrationCam);
181
182 // Fit Probability
183 TCanvas &c12 = fDisplay->AddTab("Fit Prob.");
184 c12.Divide(1, 3);
185
186 c12.cd(1);
187 gPad->SetBorderMode(0);
188 MHCamera *obj1=(MHCamera*)disp5.DrawCopy("hist");
189 obj1->AddNotify(fCalibrationCam);
190
191 c12.cd(2);
192 gPad->SetBorderMode(0);
193 obj1->Draw();
194
195 c12.cd(3);
196 gPad->SetBorderMode(0);
197 TH1D *obj2 = (TH1D*)obj1->Projection();
198 obj2->Draw();
199 obj2->SetBit(kCanDelete);
200 obj2->Fit("pol0","Q");
201 obj2->GetFunction("pol0")->SetLineColor(kYellow);
202
203 // Times
204 CamDraw(fDisplay->AddTab("FitdTimes"), disp6, disp8, fCalibrationCam);
205
206 // Pedestals
207 //CamDraw(fDisplay->AddTab("Pedestals"), disp9, disp10, calcam);
208
209 // Reduced Sigmas
210 CamDraw(fDisplay->AddTab("RedSigma"), disp11, disp16, fCalibrationCam);
211
212 // F-Factor Method
213 CamDraw(fDisplay->AddTab("F-Factor"), disp12, disp13, fCalibrationCam);
214
215}
216*/
217Bool_t MJExtractSignal::WriteResult()
218{
219 if (fOutputPath.IsNull())
220 return kTRUE;
221
222 const TString oname = GetOutputFileP();
223
224 *fLog << inf << "Writing to file: " << oname << endl;
225
226 TFile file(oname, "UPDATE");
227
228 if (fDisplay && fDisplay->Write()<=0)
229 {
230 *fLog << err << "Unable to write MStatusDisplay to " << oname << endl;
231 return kFALSE;
232 }
233
234 if (fPedPhotCam.Write()<=0)
235 {
236 *fLog << err << "Unable to write MPedPhotCam to " << oname << endl;
237 return kFALSE;
238 }
239
240 if (fBadPixels.Write()<=0)
241 {
242 *fLog << err << "Unable to write MBadPixelsCam to " << oname << endl;
243 return kFALSE;
244 }
245
246 return kTRUE;
247
248}
249
250void MJExtractSignal::SetOutputPath(const char *path)
251{
252 fOutputPath = path;
253 if (fOutputPath.EndsWith("/"))
254 fOutputPath = fOutputPath(0, fOutputPath.Length()-1);
255}
256
257TString MJExtractSignal::GetOutputFileD() const
258{
259 if (!fRuns)
260 return "";
261
262 return Form("%s/%s-F3.root", (const char*)fOutputPath, (const char*)fRuns->GetRunsAsFileName());
263}
264TString MJExtractSignal::GetOutputFileP() const
265{
266 if (!fRuns)
267 return "";
268
269 return Form("%s/%s-F2.root", (const char*)fOutputPath, (const char*)fRuns->GetRunsAsFileName());
270}
271
272Bool_t MJExtractSignal::ProcessD(MPedestalCam &pedcam)
273{
274 const TString fname = GetOutputFileD();
275
276 if (gSystem->AccessPathName(fname, kFileExists))
277 return ProcessFileD(pedcam);
278
279 return kTRUE;
280}
281
282Bool_t MJExtractSignal::ProcessFileD(MPedestalCam &pedcam)
283{
284 if (!fRuns)
285 {
286 *fLog << err << "No Runs choosen... abort." << endl;
287 return kFALSE;
288 }
289 if (fRuns->GetNumRuns() != fRuns->GetNumEntries())
290 {
291 *fLog << err << "Number of files found doesn't match number of runs... abort." << endl;
292 return kFALSE;
293 }
294
295 *fLog << inf;
296 fLog->Separator(GetDescriptor());
297 *fLog << "Calculate MExtractedSignalCam from Runs " << fRuns->GetRunsAsString() << endl;
298 *fLog << endl;
299
300 // Setup Lists
301 MParList plist;
302 plist.AddToList(&pedcam);
303
304 MTaskList tlist;
305 plist.AddToList(&tlist);
306
307 // Setup Parameters
308
309 // Make sure, that at least an empty MBadPixelsCam is available
310 // This is necessary for input which don't contain a MBadPixelsCam
311 MBadPixelsCam badcam;
312 plist.AddToList(&badcam);
313
314 // Setup Task-lists
315 MReadMarsFile read("Events");
316 read.DisableAutoScheme();
317 static_cast<MRead&>(read).AddFiles(*fRuns);
318
319 MGeomApply apply; // Only necessary to craete geometry
320 MBadPixelsMerge merge(&fBadPixels);
321 MExtractSignal extract;
322
323 MHCamEvent evt("ExtSignal");
324 evt.SetType(0);
325 MFillH fill(&evt, "MExtractedSignalCam");
326
327 MWriteRootFile write(GetOutputFileD(), "RECREATE", fRuns->GetRunsAsString(), 2);
328 write.AddContainer("MExtractedSignalCam", "Events");
329 write.AddContainer("MTime", "Events");
330 write.AddContainer("MPedestalCam", "RunHeaders");
331 write.AddContainer("MRawRunHeader", "RunHeaders");
332 write.AddContainer("MBadPixelsCam", "RunHeaders");
333
334 tlist.AddToList(&read);
335 tlist.AddToList(&apply);
336 tlist.AddToList(&merge);
337 tlist.AddToList(&extract);
338 if (TestBit(kEnableGraphicalOutput))
339 tlist.AddToList(&fill);
340 tlist.AddToList(&write);
341
342 // Create and setup the eventloop
343 MEvtLoop evtloop(fName);
344 evtloop.SetParList(&plist);
345 evtloop.SetDisplay(fDisplay);
346 evtloop.SetLogStream(fLog);
347
348 // Execute first analysis
349 if (!evtloop.Eventloop())
350 {
351 *fLog << err << GetDescriptor() << ": Failed." << endl;
352 return kFALSE;
353 }
354
355 tlist.PrintStatistics();
356
357 //DisplayResult(plist);
358
359 //if (!WriteResult())
360 // return kFALSE;
361
362 *fLog << inf << GetDescriptor() << ": Done." << endl;
363
364 return kTRUE;
365}
366
367Bool_t MJExtractSignal::ReadPedPhotCam()
368{
369 const TString fname = GetOutputFileP();
370
371 if (gSystem->AccessPathName(fname, kFileExists))
372 {
373 *fLog << err << "Input file " << fname << " doesn't exist." << endl;
374 return kFALSE;
375 }
376
377 *fLog << inf << "Reading from file: " << fname << endl;
378
379 TFile file(fname, "READ");
380 if (fPedPhotCam.Read()<=0)
381 {
382 *fLog << "Unable to read MPedPhotCam from " << fname << endl;
383 return kFALSE;
384 }
385
386 if (file.FindKey("MBadPixelsCam"))
387 {
388 MBadPixelsCam bad;
389 if (bad.Read()<=0)
390 {
391 *fLog << "Unable to read MBadPixelsCam from " << fname << endl;
392 return kFALSE;
393 }
394 fBadPixels.Merge(bad);
395 }
396
397 if (fDisplay /*&& !fDisplay->GetCanvas("Pedestals")*/) // FIXME!
398 fDisplay->Read();
399
400 return kTRUE;
401}
402
403Bool_t MJExtractSignal::ProcessP(MPedestalCam &pedcam, MCalibrationChargeCam &calcam)
404{
405 if (!ReadPedPhotCam())
406 return ProcessFileP(pedcam, calcam);
407
408 return kTRUE;
409}
410
411Bool_t MJExtractSignal::ProcessFileP(MPedestalCam &pedcam, MCalibrationChargeCam &calcam)
412{
413 if (!fRuns)
414 {
415 *fLog << err << "No Runs choosen... abort." << endl;
416 return kFALSE;
417 }
418 if (fRuns->GetNumRuns() != fRuns->GetNumEntries())
419 {
420 *fLog << err << "Number of files found doesn't match number of runs... abort." << endl;
421 return kFALSE;
422 }
423
424 *fLog << inf;
425 fLog->Separator(GetDescriptor());
426 *fLog << "Calculate MExtractedSignalCam from Runs " << fRuns->GetRunsAsString() << endl;
427 *fLog << endl;
428
429 MReadMarsFile read("Events");
430 read.DisableAutoScheme();
431 static_cast<MRead&>(read).AddFiles(*fRuns);
432
433 // Setup Tasklist
434 MParList plist;
435 plist.AddToList(&pedcam);
436 plist.AddToList(&calcam);
437 plist.AddToList(&fPedPhotCam);
438 plist.AddToList(&fBadPixels);
439
440 MTaskList tlist;
441 plist.AddToList(&tlist);
442
443 MGeomApply apply; // Only necessary to craete geometry
444 MBadPixelsMerge merge(&fBadPixels);
445 MExtractSignal extract;
446 MCalibrate calib;
447 MPedPhotCalc calc;
448
449 MHCamEvent evt1("ExtOffset");
450 MHCamEvent evt2("CalOffset");
451 evt1.SetType(0);
452 evt2.SetType(0);
453 MFillH fill1(&evt1, "MExtractedSignalCam", "FillExtractedSignal");
454 MFillH fill2(&evt2, "MCerPhotEvt", "FillCerPhotEvt");
455
456 tlist.AddToList(&read);
457 tlist.AddToList(&apply);
458 tlist.AddToList(&merge);
459 tlist.AddToList(&extract);
460 if (TestBit(kEnableGraphicalOutput))
461 tlist.AddToList(&fill1);
462 tlist.AddToList(&calib);
463 tlist.AddToList(&calc);
464 if (TestBit(kEnableGraphicalOutput))
465 tlist.AddToList(&fill2);
466
467 // Create and setup the eventloop
468 MEvtLoop evtloop(fName);
469 evtloop.SetParList(&plist);
470 evtloop.SetDisplay(fDisplay);
471 evtloop.SetLogStream(fLog);
472
473 // Execute first analysis
474 if (!evtloop.Eventloop())
475 {
476 *fLog << err << GetDescriptor() << ": Failed." << endl;
477 return kFALSE;
478 }
479
480 tlist.PrintStatistics();
481
482 //DisplayResult(plist);
483
484 if (!WriteResult())
485 return kFALSE;
486
487 *fLog << inf << GetDescriptor() << ": Done." << endl;
488
489 return kTRUE;
490}
491
492/*
493Bool_t MJExtractSignal::ProcessFile(MPedestalCam *pedcam, MCalibrationChargeCam *calcam)
494{
495 if (!fRuns)
496 {
497 *fLog << err << "No Runs choosen... abort." << endl;
498 return kFALSE;
499 }
500 if (fRuns->GetNumRuns() != fRuns->GetNumEntries())
501 {
502 *fLog << err << "Number of files found doesn't match number of runs... abort." << endl;
503 return kFALSE;
504 }
505
506 Int_t type = 0;
507 if (pedcam && calcam) type = 2;
508 if (pedcam && !calcam) type = 3;
509
510 *fLog << inf;
511 fLog->Separator(GetDescriptor());
512 *fLog << "Calculating from Runs " << fRuns->GetRunsAsString() << endl;
513 *fLog << endl;
514
515 MReadMarsFile read("Events");
516 read.DisableAutoScheme();
517 static_cast<MRead&>(read).AddFiles(*fRuns);
518
519 // Setup Tasklist
520 MParList plist;
521 switch (type)
522 {
523 case 2:
524 plist.AddToList(calcam);
525 plist.AddToList(&fPedPhotCam);
526 case 3:
527 plist.AddToList(pedcam);
528 }
529
530 MTaskList tlist;
531 plist.AddToList(&tlist);
532
533 MGeomApply apply; // Only necessary to craete geometry
534 MExtractSignal extract;
535 MCalibrate calib;
536 MPedPhotCalc calc;
537
538
539 MHCamEvent evt1("ExtOffset");
540 MHCamEvent evt2("CalOffset");
541 evt1.SetType(0);
542 evt2.SetType(0);
543 MFillH fill1(&evt1, "MExtractedSignalCam", "FillExtractedSignal");
544 MFillH fill2(&evt2, "MCerPhotEvt", "FillCerPhotEvt");
545
546 tlist.AddToList(&read);
547 tlist.AddToList(&apply);
548 tlist.AddToList(&extract);
549 if (TestBit(kEnableGraphicalOutput))
550 tlist.AddToList(&fill1);
551 tlist.AddToList(&calib);
552 tlist.AddToList(&calc);
553 if (TestBit(kEnableGraphicalOutput))
554 tlist.AddToList(&fill2);
555
556
557 MHCamEvent evt("ExtSignal");
558 evt.SetType(0);
559 MFillH fill(&evt, "MExtractedSignalCam");
560
561 MWriteRootFile write(GetOutputFileD(), "RECREATE", fRuns->GetRunsAsString(), 2);
562 write.AddContainer("MExtractedSignalCam", "Events");
563 write.AddContainer("MTime", "Events");
564 write.AddContainer("MRawRunHeader", "RunHeaders");
565 write.AddContainer("MPedestalCam", "RunHeaders");
566
567 tlist.AddToList(&read);
568 tlist.AddToList(&apply);
569 tlist.AddToList(&extract);
570 if (TestBit(kEnableGraphicalOutput))
571 tlist.AddToList(&fill);
572 tlist.AddToList(&write);
573
574 // Create and setup the eventloop
575 MEvtLoop evtloop(fName);
576 evtloop.SetParList(&plist);
577 evtloop.SetDisplay(fDisplay);
578 evtloop.SetLogStream(fLog);
579
580 // Execute first analysis
581 if (!evtloop.Eventloop())
582 {
583 *fLog << err << GetDescriptor() << ": Failed." << endl;
584 return kFALSE;
585 }
586
587 tlist.PrintStatistics();
588
589 //DisplayResult(plist);
590
591 //if (!WriteResult())
592 // return kFALSE;
593
594 *fLog << inf << GetDescriptor() << ": Done." << endl;
595
596 return kTRUE;
597
598
599 // ------------------------------------------------------
600
601 MGeomApply apply; // Only necessary to craete geometry
602 MExtractSignal extract;
603 MCalibrate calib;
604 MPedPhotCalc calc;
605
606 MHCamEvent evt1("ExtOffset");
607 MHCamEvent evt2("CalOffset");
608 evt1.SetType(0);
609 evt2.SetType(0);
610 MFillH fill1(&evt1, "MExtractedSignalCam", "FillExtractedSignal");
611 MFillH fill2(&evt2, "MCerPhotEvt", "FillCerPhotEvt");
612
613 tlist.AddToList(&read);
614 tlist.AddToList(&apply);
615 tlist.AddToList(&extract);
616 if (TestBit(kEnableGraphicalOutput))
617 tlist.AddToList(&fill1);
618 tlist.AddToList(&calib);
619 tlist.AddToList(&calc);
620 if (TestBit(kEnableGraphicalOutput))
621 tlist.AddToList(&fill2);
622
623 // Create and setup the eventloop
624 MEvtLoop evtloop(fName);
625 evtloop.SetParList(&plist);
626 evtloop.SetDisplay(fDisplay);
627 evtloop.SetLogStream(fLog);
628
629 // Execute first analysis
630 if (!evtloop.Eventloop())
631 {
632 *fLog << err << GetDescriptor() << ": Failed." << endl;
633 return kFALSE;
634 }
635
636 tlist.PrintStatistics();
637
638 //DisplayResult(plist);
639
640 if (!WriteResult())
641 return kFALSE;
642
643 *fLog << inf << GetDescriptor() << ": Done." << endl;
644
645 return kTRUE;
646 }
647 */
Note: See TracBrowser for help on using the repository browser.