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

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