source: trunk/MagicSoft/Mars/mjobs/MJPedestal.cc@ 5298

Last change on this file since 5298 was 5298, checked in by tbretz, 20 years ago
*** empty log message ***
File size: 19.2 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! Author(s): Markus Gaug, 4/2004 <mailto:markus@ifae.es>
20!
21! Copyright: MAGIC Software Development, 2000-2004
22!
23!
24\* ======================================================================== */
25
26/////////////////////////////////////////////////////////////////////////////
27//
28// MJPedestal
29//
30/////////////////////////////////////////////////////////////////////////////
31#include "MJPedestal.h"
32
33#include <TF1.h>
34#include <TEnv.h>
35#include <TFile.h>
36#include <TLine.h>
37#include <TLatex.h>
38#include <TString.h>
39#include <TCanvas.h>
40#include <TSystem.h>
41#include <TLegend.h>
42
43#include "MLog.h"
44#include "MLogManip.h"
45
46#include "MTaskEnv.h"
47#include "MSequence.h"
48#include "MRunIter.h"
49#include "MParList.h"
50#include "MTaskList.h"
51#include "MEvtLoop.h"
52#include "MExtractor.h"
53
54#include "MStatusDisplay.h"
55
56#include "MGeomCam.h"
57#include "MHCamera.h"
58#include "MPedestalCam.h"
59#include "MBadPixelsCam.h"
60
61#include "MReadMarsFile.h"
62#include "MRawFileRead.h"
63#include "MGeomApply.h"
64#include "MBadPixelsMerge.h"
65#include "MPedCalcPedRun.h"
66#include "MPedCalcFromLoGain.h"
67
68ClassImp(MJPedestal);
69
70using namespace std;
71
72const Double_t MJPedestal::fgPedestalMin = 4.;
73const Double_t MJPedestal::fgPedestalMax = 16.;
74const Double_t MJPedestal::fgPedRmsMin = 0.;
75const Double_t MJPedestal::fgPedRmsMax = 20.;
76
77const Float_t MJPedestal::fgRefPedClosedLids = 9.635;
78const Float_t MJPedestal::fgRefPedExtraGalactic = 9.93;
79const Float_t MJPedestal::fgRefPedGalactic = 10.03;
80const Float_t MJPedestal::fgRefPedRmsClosedLidsInner = 1.7;
81const Float_t MJPedestal::fgRefPedRmsExtraGalacticInner = 5.6;
82const Float_t MJPedestal::fgRefPedRmsGalacticInner = 6.92;
83const Float_t MJPedestal::fgRefPedRmsClosedLidsOuter = 1.7;
84const Float_t MJPedestal::fgRefPedRmsExtraGalacticOuter = 3.35;
85const Float_t MJPedestal::fgRefPedRmsGalacticOuter = 4.2;
86
87// --------------------------------------------------------------------------
88//
89// Default constructor.
90//
91// Sets:
92// - fRuns to 0,
93// - fExtractor to NULL,
94// - fDataCheck to kFALSE,
95// - fUseData to kFALSE
96// - fStorage to Normal Storage
97//
98MJPedestal::MJPedestal(const char *name, const char *title)
99 : fRuns(0), fExtractor(NULL), fDisplayType(kNormalDisplay),
100 fDataCheck(kFALSE), fUseData(kFALSE)
101{
102 fName = name ? name : "MJPedestal";
103 fTitle = title ? title : "Tool to create a pedestal file (MPedestalCam)";
104
105 SetNormalStorage();
106}
107
108const char* MJPedestal::GetOutputFile() const
109{
110 if (fSequence.IsValid())
111 return Form("%s/pedest%06d.root", (const char*)fPathOut, fSequence.GetSequence());
112
113 if (!fRuns)
114 return "";
115
116 return Form("%s/%s-F0.root", (const char*)fPathOut, (const char*)fRuns->GetRunsAsFileName());
117}
118
119//---------------------------------------------------------------------------------
120//
121// Try to read an existing MPedestalCam from a previously created output file.
122// If found, also an MBadPixelsCam and the corresponding display is read.
123//
124// In case of Storage type "kNoStorage" or if the file is not found or the
125// MPedestalCam cannot be read, return kFALSE, otherwise kTRUE.
126//
127Bool_t MJPedestal::ReadPedestalCam()
128{
129 if (IsNoStorage())
130 return kFALSE;
131
132 const TString fname = GetOutputFile();
133
134 if (gSystem->AccessPathName(fname, kFileExists))
135 {
136 *fLog << warn << "Input file " << fname << " doesn't exist, will create it." << endl;
137 return kFALSE;
138 }
139
140 *fLog << inf << "Reading from file: " << fname << endl;
141
142 TFile file(fname, "READ");
143 if (fPedestalCam.Read()<=0)
144 {
145 *fLog << err << "Unable to read MPedestalCam from " << fname << endl;
146 return kFALSE;
147 }
148
149 if (file.FindKey("MBadPixelsCam"))
150 {
151 MBadPixelsCam bad;
152 if (bad.Read()<=0)
153 {
154 *fLog << err << "Unable to read MBadPixelsCam from " << fname << endl;
155 return kFALSE;
156 }
157 fBadPixels.Merge(bad);
158 }
159
160 if (fDisplay && !fDisplay->GetCanvas("Pedestals"))
161 fDisplay->Read();
162
163 return kTRUE;
164}
165
166//---------------------------------------------------------------------------------
167//
168// Display the results.
169// If Display type "kDataCheck" was chosen, also the reference lines are displayed.
170//
171void MJPedestal::DisplayResult(MParList &plist)
172{
173 if (!fDisplay)
174 return;
175
176 //
177 // Update display
178 //
179 TString title = fDisplay->GetTitle();
180 title += "-- Pedestal ";
181 if (fSequence.IsValid())
182 title += fSequence.GetName();
183 else
184 if (fRuns) // FIXME: What to do if an environmentfile was used?
185 title += fRuns->GetRunsAsString();
186 title += " --";
187 fDisplay->SetTitle(title);
188
189 //
190 // Get container from list
191 //
192 MGeomCam &geomcam = *(MGeomCam*)plist.FindObject("MGeomCam");
193
194 //
195 // Create container to display
196 //
197 MHCamera disp0(geomcam, "MPedestalCam;ped", "Mean Pedestal");
198 MHCamera disp1(geomcam, "MPedestalCam;RMS", "Pedestal RMS");
199
200 disp0.SetCamContent(fPedestalCam, 0);
201 disp0.SetCamError (fPedestalCam, 1);
202
203 disp1.SetCamContent(fPedestalCam, 2);
204 disp1.SetCamError (fPedestalCam, 3);
205
206 disp0.SetYTitle("Pedestal [counts/slice]");
207 disp1.SetYTitle("RMS [counts/slice]");
208
209 //
210 // Display data
211 //
212 TCanvas &c3 = fDisplay->AddTab("Pedestals");
213 c3.Divide(2,3);
214
215 if (fDisplayType != kDataCheckDisplay)
216 {
217 disp0.CamDraw(c3, 1, 2, 1);
218 disp1.CamDraw(c3, 2, 2, 6);
219 return;
220 }
221
222 c3.cd(1);
223 gPad->SetBorderMode(0);
224 gPad->SetTicks();
225 MHCamera *obj1=(MHCamera*)disp0.DrawCopy("hist");
226 //
227 // for the datacheck, fix the ranges!!
228 //
229 obj1->SetMinimum(fgPedestalMin);
230 obj1->SetMaximum(fgPedestalMax);
231
232 //
233 // Set the datacheck sizes:
234 //
235 FixDataCheckHist((TH1D*)obj1);
236 //
237 // set reference lines
238 //
239 DisplayReferenceLines(obj1,0);
240
241 //
242 // end reference lines
243 //
244 c3.cd(3);
245 gPad->SetBorderMode(0);
246 obj1->SetPrettyPalette();
247 obj1->Draw();
248
249 c3.cd(5);
250 gPad->SetBorderMode(0);
251 gPad->SetTicks();
252 TH1D *obj2 = (TH1D*)obj1->Projection(obj1->GetName());
253 obj2->Draw();
254 obj2->SetBit(kCanDelete);
255 obj2->Fit("gaus","Q");
256 obj2->GetFunction("gaus")->SetLineColor(kYellow);
257 //
258 // Set the datacheck sizes:
259 //
260 FixDataCheckHist(obj2);
261 obj2->SetStats(1);
262
263 c3.cd(2);
264 gPad->SetBorderMode(0);
265 gPad->SetTicks();
266 MHCamera *obj3=(MHCamera*)disp1.DrawCopy("hist");
267 //
268 // for the datacheck, fix the ranges!!
269 //
270 obj3->SetMinimum(fgPedRmsMin);
271 obj3->SetMaximum(fgPedRmsMax);
272 //
273 // Set the datacheck sizes:
274 //
275 FixDataCheckHist((TH1D*)obj3);
276 //
277 // set reference lines
278 //
279 DisplayReferenceLines(obj1,1);
280
281 c3.cd(4);
282 gPad->SetBorderMode(0);
283 obj3->SetPrettyPalette();
284 obj3->Draw();
285
286 c3.cd(6);
287 gPad->SetBorderMode(0);
288
289 if (geomcam.InheritsFrom("MGeomCamMagic"))
290 {
291 TArrayI inner(1);
292 inner[0] = 0;
293
294 TArrayI outer(1);
295 outer[0] = 1;
296
297 TArrayI s0(6);
298 s0[0] = 6;
299 s0[1] = 1;
300 s0[2] = 2;
301 s0[3] = 3;
302 s0[4] = 4;
303 s0[5] = 5;
304
305 TArrayI s1(3);
306 s1[0] = 6;
307 s1[1] = 1;
308 s1[2] = 2;
309
310 TArrayI s2(3);
311 s2[0] = 3;
312 s2[1] = 4;
313 s2[2] = 5;
314
315 TVirtualPad *pad = gPad;
316 pad->Divide(2,1);
317
318 TH1D *inout[2];
319 inout[0] = disp1.ProjectionS(s0, inner, "Inner");
320 inout[1] = disp1.ProjectionS(s0, outer, "Outer");
321 FixDataCheckHist(inout[0]);
322 FixDataCheckHist(inout[1]);
323
324 inout[0]->SetTitle(Form("%s %s",disp1.GetTitle(),"Inner"));
325 inout[1]->SetTitle(Form("%s %s",disp1.GetTitle(),"Outer"));
326
327
328 for (int i=0; i<2; i++)
329 {
330 pad->cd(i+1);
331 gPad->SetBorderMode(0);
332 gPad->SetTicks();
333
334 inout[i]->SetDirectory(NULL);
335 inout[i]->SetLineColor(kRed+i);
336 inout[i]->SetBit(kCanDelete);
337 inout[i]->Draw();
338 inout[i]->Fit("gaus", "Q");
339
340 TLegend *leg2 = new TLegend(0.6,0.2,0.9,0.55);
341 leg2->SetHeader(inout[i]->GetName());
342 leg2->AddEntry(inout[i], inout[i]->GetName(), "l");
343
344 //
345 // Display the outliers as dead and noisy pixels
346 //
347 DisplayOutliers(inout[i]);
348
349 //
350 // Display the two half of the camera separately
351 //
352 TH1D *half[2];
353 half[0] = disp1.ProjectionS(s1, i==0 ? inner : outer , "Sector 6-1-2");
354 half[1] = disp1.ProjectionS(s2, i==0 ? inner : outer , "Sector 3-4-5");
355
356 for (int j=0; j<2; j++)
357 {
358 half[j]->SetLineColor(kRed+i+2*j+1);
359 half[j]->SetDirectory(NULL);
360 half[j]->SetBit(kCanDelete);
361 half[j]->Draw("same");
362 leg2->AddEntry(half[j], half[j]->GetName(), "l");
363 }
364 leg2->Draw();
365 }
366 }
367}
368
369void MJPedestal::DisplayReferenceLines(MHCamera *cam, const Int_t what) const
370{
371 Double_t x = cam->GetNbinsX();
372
373 const MGeomCam *geom = cam->GetGeometry();
374
375 if (geom->InheritsFrom("MGeomCamMagic"))
376 x = what ? 397 : cam->GetNbinsX();
377
378 TLine line;
379 line.SetLineStyle(kDashed);
380 line.SetLineWidth(3);
381
382 line.SetLineColor(kBlue);
383 TLine *l1 = line.DrawLine(0, what ? fgRefPedRmsGalacticInner : fgRefPedGalactic,
384 x, what ? fgRefPedRmsGalacticInner : fgRefPedGalactic);
385
386 line.SetLineColor(kYellow);
387 TLine *l2 = line.DrawLine(0, what ? fgRefPedRmsExtraGalacticInner : fgRefPedExtraGalactic,
388 x, what ? fgRefPedRmsExtraGalacticInner : fgRefPedExtraGalactic);
389
390 line.SetLineColor(kMagenta);
391 TLine *l3 = line.DrawLine(0, what ? fgRefPedRmsClosedLidsInner : fgRefPedClosedLids,
392 x, what ? fgRefPedRmsClosedLidsInner : fgRefPedClosedLids);
393
394 if (geom->InheritsFrom("MGeomCamMagic"))
395 if (what)
396 {
397 const Double_t x2 = cam->GetNbinsX();
398
399 line.SetLineColor(kBlue);
400 line.DrawLine(398, fgRefPedRmsGalacticOuter,
401 x2, fgRefPedRmsGalacticOuter);
402
403 line.SetLineColor(kYellow);
404 line.DrawLine(398, fgRefPedRmsExtraGalacticOuter,
405 x2, fgRefPedRmsExtraGalacticOuter);
406
407 line.SetLineColor(kMagenta);
408 line.DrawLine(398, fgRefPedRmsClosedLidsOuter,
409 x2, fgRefPedRmsClosedLidsOuter);
410 }
411
412
413 TLegend *leg = new TLegend(0.4,0.75,0.7,0.99);
414 leg->SetBit(kCanDelete);
415 leg->AddEntry(l1, "Galactic Source","l");
416 leg->AddEntry(l2, "Extra-Galactic Source","l");
417 leg->AddEntry(l3, "Closed Lids","l");
418 leg->Draw();
419}
420
421void MJPedestal::DisplayOutliers(TH1D *hist) const
422{
423 const Float_t mean = hist->GetFunction("gaus")->GetParameter(1);
424 const Float_t lolim = mean - 3.5*hist->GetFunction("gaus")->GetParameter(2);
425 const Float_t uplim = mean + 3.5*hist->GetFunction("gaus")->GetParameter(2);
426 const Stat_t dead = hist->Integral(0,hist->FindBin(lolim)-1);
427 const Stat_t noisy = hist->Integral(hist->FindBin(uplim)+1,hist->GetNbinsX()+1);
428
429 TLatex deadtex;
430 deadtex.SetTextSize(0.06);
431 deadtex.DrawLatex(0.1,hist->GetBinContent(hist->GetMaximumBin())/1.1,Form("%3i dead pixels",(Int_t)dead));
432
433 TLatex noisytex;
434 noisytex.SetTextSize(0.06);
435 noisytex.DrawLatex(0.1,hist->GetBinContent(hist->GetMaximumBin())/1.2,Form("%3i noisy pixels",(Int_t)noisy));
436}
437
438void MJPedestal::FixDataCheckHist(TH1D *hist) const
439{
440 hist->SetDirectory(NULL);
441 hist->SetStats(0);
442
443 //
444 // set the labels bigger
445 //
446 TAxis *xaxe = hist->GetXaxis();
447 TAxis *yaxe = hist->GetYaxis();
448
449 xaxe->CenterTitle();
450 yaxe->CenterTitle();
451 xaxe->SetTitleSize(0.06);
452 yaxe->SetTitleSize(0.06);
453 xaxe->SetTitleOffset(0.8);
454 yaxe->SetTitleOffset(0.5);
455 xaxe->SetLabelSize(0.05);
456 yaxe->SetLabelSize(0.05);
457}
458
459/*
460Bool_t MJPedestal::WriteEventloop(MEvtLoop &evtloop) const
461{
462 if (fOutputPath.IsNull())
463 return kTRUE;
464
465 const TString oname(GetOutputFile());
466
467 *fLog << inf << "Writing to file: " << oname << endl;
468
469 TFile file(oname, fOverwrite?"RECREATE":"NEW", "File created by MJPedestal", 9);
470 if (!file.IsOpen())
471 {
472 *fLog << err << "ERROR - Couldn't open file " << oname << " for writing..." << endl;
473 return kFALSE;
474 }
475
476 if (evtloop.Write(fName)<=0)
477 {
478 *fLog << err << "Unable to write MEvtloop to " << oname << endl;
479 return kFALSE;
480 }
481
482 return kTRUE;
483}
484*/
485
486// --------------------------------------------------------------------------
487//
488// The following resource options are available:
489// Prefix.DataCheckDisplay: Yes, No
490// Prefix.DataCheck: Yes, No
491// Prefix.UseData: Yes, No
492// Prefix.DisableOutput: Yes, No
493//
494Bool_t MJPedestal::CheckEnvLocal()
495{
496 if (HasEnv("DataCheckDisplay"))
497 fDisplayType = GetEnv("DataCheckDisplay", kFALSE) ? kDataCheckDisplay : kNormalDisplay;
498
499 SetDataCheck(GetEnv("DataCheck", fDataCheck));
500 SetUseData(GetEnv("UseData", fUseData));
501 SetNoStorage(GetEnv("DisableOutput", IsNoStorage()));
502
503 return kTRUE;
504}
505
506//---------------------------------------------------------------------------------
507//
508// Try to write the created MPedestalCam in the output file.
509// If possible, also an MBadPixelsCam and the corresponding display is written.
510//
511// In case of Storage type "kNoStorage" or if any of the containers
512// cannot be written, return kFALSE, otherwise kTRUE.
513//
514Bool_t MJPedestal::WriteResult()
515{
516 if (IsNoStorage())
517 return kTRUE;
518
519 if (fPathOut.IsNull())
520 return kTRUE;
521
522 const TString oname(GetOutputFile());
523
524 *fLog << inf << "Writing to file: " << oname << endl;
525
526 TFile file(oname, "UPDATE", "File created by MJPedestal", 9);
527 if (!file.IsOpen())
528 {
529 *fLog << err << "ERROR - Couldn't open file " << oname << " for writing..." << endl;
530 return kFALSE;
531 }
532
533 if (fDisplay && fDisplay->Write()<=0)
534 {
535 *fLog << err << "Unable to write MStatusDisplay to " << oname << endl;
536 return kFALSE;
537 }
538
539 if (fPedestalCam.Write()<=0)
540 {
541 *fLog << err << "Unable to write MPedestalCam to " << oname << endl;
542 return kFALSE;
543 }
544
545 if (fBadPixels.Write()<=0)
546 {
547 *fLog << err << "Unable to write MBadPixelsCam to " << oname << endl;
548 return kFALSE;
549 }
550
551 return kTRUE;
552}
553
554Bool_t MJPedestal::Process()
555{
556 if (!ReadPedestalCam())
557 return ProcessFile();
558
559 return kTRUE;
560}
561
562Bool_t MJPedestal::ProcessFile()
563{
564 if (!fSequence.IsValid())
565 {
566 if (!fRuns)
567 {
568 *fLog << err << "Neither AddRuns nor SetSequence nor SetEnv was called... abort." << endl;
569 return kFALSE;
570 }
571 if (fRuns && fRuns->GetNumRuns() != fRuns->GetNumEntries())
572 {
573 *fLog << err << "Number of files found doesn't match number of runs... abort." << endl;
574 return kFALSE;
575 }
576 }
577
578 //if (!CheckEnv())
579 // return kFALSE;
580
581 CheckEnv();
582
583 // --------------------------------------------------------------------------------
584
585 const TString type = fUseData ? "data" : "pedestal";
586
587 *fLog << inf;
588 fLog->Separator(GetDescriptor());
589 *fLog << "Calculate MPedestalCam from " << type << "-runs ";
590 if (fSequence.IsValid())
591 *fLog << fSequence.GetName() << endl;
592 else
593 if (fRuns)
594 *fLog << fRuns->GetRunsAsString() << endl;
595 else
596 *fLog << "in Resource File." << endl;
597 *fLog << endl;
598
599 // --------------------------------------------------------------------------------
600
601 MParList plist;
602 MTaskList tlist;
603 plist.AddToList(&tlist);
604 plist.AddToList(this); // take care of fDisplay!
605
606 MReadMarsFile read("Events");
607 MRawFileRead rawread(NULL);
608
609 MDirIter iter;
610 if (fSequence.IsValid())
611 {
612 const Int_t n0 = fUseData ? fSequence.SetupDatRuns(iter, fPathData, fDataCheck) : fSequence.SetupPedRuns(iter, fPathData, fDataCheck);
613 const Int_t n1 = fUseData ? fSequence.GetNumDatRuns() : fSequence.GetNumPedRuns();
614 if (n0==0)
615 {
616 *fLog << err << "ERROR - No " << type << " input files of sequence found in " << fPathData << endl;
617 return kFALSE;
618 }
619 if (n0!=n1)
620 {
621 *fLog << err << "ERROR - Number of " << type << " files found (" << n0 << ") in " << fPathData << " doesn't match number of files in sequence (" << n1 << ")" << endl;
622 return kFALSE;
623 }
624 }
625
626 if (fDataCheck)
627 {
628 if (fRuns || fSequence.IsValid())
629 rawread.AddFiles(fSequence.IsValid() ? iter : *fRuns);
630 tlist.AddToList(&rawread);
631 }
632 else
633 {
634 read.DisableAutoScheme();
635 if (fRuns || fSequence.IsValid())
636 read.AddFiles(fSequence.IsValid() ? iter : *fRuns);
637 tlist.AddToList(&read);
638 }
639
640 // Setup Tasklist
641 plist.AddToList(&fPedestalCam);
642 plist.AddToList(&fBadPixels);
643
644 MGeomApply geomapl;
645 MBadPixelsMerge merge(&fBadPixels);
646
647 MPedCalcPedRun pedcalc;
648 MPedCalcFromLoGain pedlogain;
649 pedlogain.SetPedestalUpdate(kFALSE);
650
651 MTaskEnv taskenv("ExtractPedestal");
652 if (IsUseData())
653 taskenv.SetDefault(&pedlogain);
654 else
655 taskenv.SetDefault(&pedcalc);
656
657 if (fExtractor)
658 {
659 if (IsUseData())
660 pedlogain.SetWindowSize(12,(Int_t)fExtractor->GetNumHiGainSamples());
661 else
662 {
663 pedcalc.SetWindowSize((Int_t)fExtractor->GetNumHiGainSamples());
664 pedcalc.SetRange(fExtractor->GetHiGainFirst(), fExtractor->GetHiGainLast());
665 }
666 }
667 /*
668 else
669 {
670 *fLog << warn << GetDescriptor();
671 *fLog << ": No extractor has been chosen, take default number of FADC samples " << endl;
672 }
673 */
674
675 tlist.AddToList(&geomapl);
676 tlist.AddToList(&merge);
677 tlist.AddToList(&taskenv);
678
679 //
680 // Execute the eventloop
681 //
682 MEvtLoop evtloop(fName);
683 evtloop.SetParList(&plist);
684 evtloop.SetDisplay(fDisplay);
685 evtloop.SetLogStream(fLog);
686 if (!SetupEnv(evtloop))
687 return kFALSE;
688
689 // if (!WriteEventloop(evtloop))
690 // return kFALSE;
691
692 // Execute first analysis
693 if (!evtloop.Eventloop(fMaxEvents))
694 {
695 *fLog << err << GetDescriptor() << ": Failed." << endl;
696 return kFALSE;
697 }
698
699 tlist.PrintStatistics();
700
701 DisplayResult(plist);
702
703 if (!WriteResult())
704 return kFALSE;
705
706 *fLog << all << GetDescriptor() << ": Done." << endl;
707 *fLog << endl << endl;
708
709 return kTRUE;
710}
Note: See TracBrowser for help on using the repository browser.