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

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