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

Last change on this file since 3501 was 3501, checked in by tbretz, 21 years ago
*** empty log message ***
File size: 17.7 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("MRawEvtHeader", "Events");
331 write.AddContainer("MPedestalCam", "RunHeaders");
332 write.AddContainer("MRawRunHeader", "RunHeaders");
333 write.AddContainer("MBadPixelsCam", "RunHeaders");
334
335 tlist.AddToList(&read);
336 tlist.AddToList(&apply);
337 tlist.AddToList(&merge);
338 tlist.AddToList(&extract);
339 if (TestBit(kEnableGraphicalOutput))
340 tlist.AddToList(&fill);
341 tlist.AddToList(&write);
342
343 // Create and setup the eventloop
344 MEvtLoop evtloop(fName);
345 evtloop.SetParList(&plist);
346 evtloop.SetDisplay(fDisplay);
347 evtloop.SetLogStream(fLog);
348
349 // Execute first analysis
350 if (!evtloop.Eventloop())
351 {
352 *fLog << err << GetDescriptor() << ": Failed." << endl;
353 return kFALSE;
354 }
355
356 tlist.PrintStatistics();
357
358 //DisplayResult(plist);
359
360 //if (!WriteResult())
361 // return kFALSE;
362
363 *fLog << inf << GetDescriptor() << ": Done." << endl;
364
365 return kTRUE;
366}
367
368Bool_t MJExtractSignal::ReadPedPhotCam()
369{
370 const TString fname = GetOutputFileP();
371
372 if (gSystem->AccessPathName(fname, kFileExists))
373 {
374 *fLog << err << "Input file " << fname << " doesn't exist." << endl;
375 return kFALSE;
376 }
377
378 *fLog << inf << "Reading from file: " << fname << endl;
379
380 TFile file(fname, "READ");
381 if (fPedPhotCam.Read()<=0)
382 {
383 *fLog << "Unable to read MPedPhotCam from " << fname << endl;
384 return kFALSE;
385 }
386
387 if (file.FindKey("MBadPixelsCam"))
388 {
389 MBadPixelsCam bad;
390 if (bad.Read()<=0)
391 {
392 *fLog << "Unable to read MBadPixelsCam from " << fname << endl;
393 return kFALSE;
394 }
395 fBadPixels.Merge(bad);
396 }
397
398 if (fDisplay /*&& !fDisplay->GetCanvas("Pedestals")*/) // FIXME!
399 fDisplay->Read();
400
401 return kTRUE;
402}
403
404Bool_t MJExtractSignal::ProcessP(MPedestalCam &pedcam, MCalibrationChargeCam &calcam)
405{
406 if (!ReadPedPhotCam())
407 return ProcessFileP(pedcam, calcam);
408
409 return kTRUE;
410}
411
412Bool_t MJExtractSignal::ProcessFileP(MPedestalCam &pedcam, MCalibrationChargeCam &calcam)
413{
414 if (!fRuns)
415 {
416 *fLog << err << "No Runs choosen... abort." << endl;
417 return kFALSE;
418 }
419 if (fRuns->GetNumRuns() != fRuns->GetNumEntries())
420 {
421 *fLog << err << "Number of files found doesn't match number of runs... abort." << endl;
422 return kFALSE;
423 }
424
425 *fLog << inf;
426 fLog->Separator(GetDescriptor());
427 *fLog << "Calculate MExtractedSignalCam from Runs " << fRuns->GetRunsAsString() << endl;
428 *fLog << endl;
429
430 MReadMarsFile read("Events");
431 read.DisableAutoScheme();
432 static_cast<MRead&>(read).AddFiles(*fRuns);
433
434 // Setup Tasklist
435 MParList plist;
436 plist.AddToList(&pedcam);
437 plist.AddToList(&calcam);
438 plist.AddToList(&fPedPhotCam);
439 plist.AddToList(&fBadPixels);
440
441 MTaskList tlist;
442 plist.AddToList(&tlist);
443
444 MGeomApply apply; // Only necessary to craete geometry
445 MBadPixelsMerge merge(&fBadPixels);
446 MExtractSignal extract;
447 MCalibrate calib;
448 MPedPhotCalc calc;
449
450 MHCamEvent evt1("ExtOffset");
451 MHCamEvent evt2("CalOffset");
452 evt1.SetType(0);
453 evt2.SetType(0);
454 MFillH fill1(&evt1, "MExtractedSignalCam", "FillExtractedSignal");
455 MFillH fill2(&evt2, "MCerPhotEvt", "FillCerPhotEvt");
456
457 tlist.AddToList(&read);
458 tlist.AddToList(&apply);
459 tlist.AddToList(&merge);
460 tlist.AddToList(&extract);
461 if (TestBit(kEnableGraphicalOutput))
462 tlist.AddToList(&fill1);
463 tlist.AddToList(&calib);
464 tlist.AddToList(&calc);
465 if (TestBit(kEnableGraphicalOutput))
466 tlist.AddToList(&fill2);
467
468 // Create and setup the eventloop
469 MEvtLoop evtloop(fName);
470 evtloop.SetParList(&plist);
471 evtloop.SetDisplay(fDisplay);
472 evtloop.SetLogStream(fLog);
473
474 // Execute first analysis
475 if (!evtloop.Eventloop())
476 {
477 *fLog << err << GetDescriptor() << ": Failed." << endl;
478 return kFALSE;
479 }
480
481 tlist.PrintStatistics();
482
483 //DisplayResult(plist);
484
485 if (!WriteResult())
486 return kFALSE;
487
488 *fLog << inf << GetDescriptor() << ": Done." << endl;
489
490 return kTRUE;
491}
492
493/*
494Bool_t MJExtractSignal::ProcessFile(MPedestalCam *pedcam, MCalibrationChargeCam *calcam)
495{
496 if (!fRuns)
497 {
498 *fLog << err << "No Runs choosen... abort." << endl;
499 return kFALSE;
500 }
501 if (fRuns->GetNumRuns() != fRuns->GetNumEntries())
502 {
503 *fLog << err << "Number of files found doesn't match number of runs... abort." << endl;
504 return kFALSE;
505 }
506
507 Int_t type = 0;
508 if (pedcam && calcam) type = 2;
509 if (pedcam && !calcam) type = 3;
510
511 *fLog << inf;
512 fLog->Separator(GetDescriptor());
513 *fLog << "Calculating from Runs " << fRuns->GetRunsAsString() << endl;
514 *fLog << endl;
515
516 MReadMarsFile read("Events");
517 read.DisableAutoScheme();
518 static_cast<MRead&>(read).AddFiles(*fRuns);
519
520 // Setup Tasklist
521 MParList plist;
522 switch (type)
523 {
524 case 2:
525 plist.AddToList(calcam);
526 plist.AddToList(&fPedPhotCam);
527 case 3:
528 plist.AddToList(pedcam);
529 }
530
531 MTaskList tlist;
532 plist.AddToList(&tlist);
533
534 MGeomApply apply; // Only necessary to craete geometry
535 MExtractSignal extract;
536 MCalibrate calib;
537 MPedPhotCalc calc;
538
539
540 MHCamEvent evt1("ExtOffset");
541 MHCamEvent evt2("CalOffset");
542 evt1.SetType(0);
543 evt2.SetType(0);
544 MFillH fill1(&evt1, "MExtractedSignalCam", "FillExtractedSignal");
545 MFillH fill2(&evt2, "MCerPhotEvt", "FillCerPhotEvt");
546
547 tlist.AddToList(&read);
548 tlist.AddToList(&apply);
549 tlist.AddToList(&extract);
550 if (TestBit(kEnableGraphicalOutput))
551 tlist.AddToList(&fill1);
552 tlist.AddToList(&calib);
553 tlist.AddToList(&calc);
554 if (TestBit(kEnableGraphicalOutput))
555 tlist.AddToList(&fill2);
556
557
558 MHCamEvent evt("ExtSignal");
559 evt.SetType(0);
560 MFillH fill(&evt, "MExtractedSignalCam");
561
562 MWriteRootFile write(GetOutputFileD(), "RECREATE", fRuns->GetRunsAsString(), 2);
563 write.AddContainer("MExtractedSignalCam", "Events");
564 write.AddContainer("MTime", "Events");
565 write.AddContainer("MRawRunHeader", "RunHeaders");
566 write.AddContainer("MPedestalCam", "RunHeaders");
567
568 tlist.AddToList(&read);
569 tlist.AddToList(&apply);
570 tlist.AddToList(&extract);
571 if (TestBit(kEnableGraphicalOutput))
572 tlist.AddToList(&fill);
573 tlist.AddToList(&write);
574
575 // Create and setup the eventloop
576 MEvtLoop evtloop(fName);
577 evtloop.SetParList(&plist);
578 evtloop.SetDisplay(fDisplay);
579 evtloop.SetLogStream(fLog);
580
581 // Execute first analysis
582 if (!evtloop.Eventloop())
583 {
584 *fLog << err << GetDescriptor() << ": Failed." << endl;
585 return kFALSE;
586 }
587
588 tlist.PrintStatistics();
589
590 //DisplayResult(plist);
591
592 //if (!WriteResult())
593 // return kFALSE;
594
595 *fLog << inf << GetDescriptor() << ": Done." << endl;
596
597 return kTRUE;
598
599
600 // ------------------------------------------------------
601
602 MGeomApply apply; // Only necessary to craete geometry
603 MExtractSignal extract;
604 MCalibrate calib;
605 MPedPhotCalc calc;
606
607 MHCamEvent evt1("ExtOffset");
608 MHCamEvent evt2("CalOffset");
609 evt1.SetType(0);
610 evt2.SetType(0);
611 MFillH fill1(&evt1, "MExtractedSignalCam", "FillExtractedSignal");
612 MFillH fill2(&evt2, "MCerPhotEvt", "FillCerPhotEvt");
613
614 tlist.AddToList(&read);
615 tlist.AddToList(&apply);
616 tlist.AddToList(&extract);
617 if (TestBit(kEnableGraphicalOutput))
618 tlist.AddToList(&fill1);
619 tlist.AddToList(&calib);
620 tlist.AddToList(&calc);
621 if (TestBit(kEnableGraphicalOutput))
622 tlist.AddToList(&fill2);
623
624 // Create and setup the eventloop
625 MEvtLoop evtloop(fName);
626 evtloop.SetParList(&plist);
627 evtloop.SetDisplay(fDisplay);
628 evtloop.SetLogStream(fLog);
629
630 // Execute first analysis
631 if (!evtloop.Eventloop())
632 {
633 *fLog << err << GetDescriptor() << ": Failed." << endl;
634 return kFALSE;
635 }
636
637 tlist.PrintStatistics();
638
639 //DisplayResult(plist);
640
641 if (!WriteResult())
642 return kFALSE;
643
644 *fLog << inf << GetDescriptor() << ": Done." << endl;
645
646 return kTRUE;
647 }
648 */
Note: See TracBrowser for help on using the repository browser.