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

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