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

Last change on this file since 6321 was 6302, checked in by gaug, 20 years ago
*** empty log message ***
File size: 34.1 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-2005
22!
23!
24\* ======================================================================== */
25
26/////////////////////////////////////////////////////////////////////////////
27//
28// MJPedestal
29//
30// Resource file entries are case sensitive!
31//
32/////////////////////////////////////////////////////////////////////////////
33#include "MJPedestal.h"
34
35// root classes
36#include <TF1.h>
37#include <TEnv.h>
38#include <TFile.h>
39#include <TLine.h>
40#include <TLatex.h>
41#include <TString.h>
42#include <TCanvas.h>
43#include <TSystem.h>
44#include <TLegend.h>
45#include <TPad.h>
46#include <TEnv.h>
47#include <TH2F.h>
48
49// mars core
50#include "MLog.h"
51#include "MLogManip.h"
52
53#include "MTaskEnv.h"
54#include "MSequence.h"
55#include "MRunIter.h"
56#include "MParList.h"
57#include "MTaskList.h"
58#include "MEvtLoop.h"
59
60#include "MStatusDisplay.h"
61
62// Other basic classes
63#include "MExtractor.h"
64#include "MExtractTimeAndCharge.h"
65#include "MFEventSelector.h"
66
67// parameter containers
68#include "MGeomCam.h"
69#include "MHCamera.h"
70#include "MPedestalCam.h"
71#include "MPedestalPix.h"
72#include "MBadPixelsCam.h"
73
74#include "MCalibrationPedCam.h"
75#include "MCalibrationPix.h"
76#include "MHPedestalCam.h"
77#include "MHPedestalPix.h"
78
79// tasks
80#include "MReadMarsFile.h"
81#include "MRawFileRead.h"
82#include "MGeomApply.h"
83#include "MBadPixelsMerge.h"
84#include "MFillH.h"
85#include "MPedCalcPedRun.h"
86#include "MPedCalcFromLoGain.h"
87
88#include <fstream>
89
90ClassImp(MJPedestal);
91
92using namespace std;
93
94const TString MJPedestal::fgReferenceFile = "mjobs/pedestalref.rc";
95const TString MJPedestal::fgBadPixelsFile = "mjobs/badpixels_0_559.rc";
96// --------------------------------------------------------------------------
97//
98// Default constructor.
99//
100// Sets:
101// - fExtractor to NULL,
102// - fExtractType to kUsePedRun
103// - fStorage to Normal Storage
104// - fExtractorResolution to kFALSE
105//
106MJPedestal::MJPedestal(const char *name, const char *title)
107 : fExtractor(NULL), fDisplayType(kDisplayDataCheck),
108 fExtractType(kUsePedRun), fExtractionType(kFundamental)
109{
110 fName = name ? name : "MJPedestal";
111 fTitle = title ? title : "Tool to create a pedestal file (MPedestalCam)";
112
113 SetNormalStorage();
114 SetUsePedRun();
115 SetPathIn("");
116 SetReferenceFile();
117 SetBadPixelsFile();
118
119 //
120 // Default references for case that no value references file is there
121 // (should not occur)
122 //
123
124 fPedestalMin = 4.;
125 fPedestalMax = 16.;
126 fPedRmsMin = 0.;
127 fPedRmsMax = 20.;
128 fRefPedClosedLids = 9.635;
129 fRefPedExtraGalactic = 9.93;
130 fRefPedGalactic = 10.03;
131 fRefPedRmsClosedLidsInner = 1.7;
132 fRefPedRmsExtraGalacticInner = 5.6;
133 fRefPedRmsGalacticInner = 6.92;
134 fRefPedRmsClosedLidsOuter = 1.7;
135 fRefPedRmsExtraGalacticOuter = 3.35;
136 fRefPedRmsGalacticOuter = 4.2;
137}
138
139MJPedestal::~MJPedestal()
140{
141 if (fExtractor)
142 delete fExtractor;
143}
144
145const char* MJPedestal::GetOutputFile() const
146{
147 const TString name(GetOutputFileName());
148 if (name.IsNull())
149 return "";
150
151 return Form("%s/%s", fPathOut.Data(), name.Data());
152}
153
154const char* MJPedestal::GetOutputFileName() const
155{
156
157 if (fSequence.IsValid())
158 return Form("pedest%08d.root", fSequence.GetSequence());
159
160 if (!fRuns)
161 return "";
162
163 return Form("%s-F0.root", (const char*)fRuns->GetRunsAsFileName());
164}
165
166//---------------------------------------------------------------------------------
167//
168// Try to read an existing MPedestalCam from a previously created output file.
169// If found, also an MBadPixelsCam and the corresponding display is read.
170//
171// In case of Storage type "kNoStorage" or if the file is not found or the
172// MPedestalCam cannot be read, return kFALSE, otherwise kTRUE.
173//
174Bool_t MJPedestal::ReadPedestalCam()
175{
176 const TString fname = GetOutputFile();
177
178 *fLog << inf << "Reading pedestals from file: " << fname << endl;
179
180 TFile file(fname, "READ");
181 if (fPedestalCamIn.Read()<=0)
182 {
183 *fLog << err << "Unable to read incoming MPedestalCam from " << fname << endl;
184 return kFALSE;
185 }
186
187 if (fPedestalCamOut.Read()<=0)
188 {
189 *fLog << err << "Unable to read outgoing MPedestalCam from " << fname << endl;
190 return kFALSE;
191 }
192
193 if (file.FindKey("MBadPixelsCam"))
194 {
195 MBadPixelsCam bad;
196 if (bad.Read()<=0)
197 {
198 *fLog << err << "Unable to read MBadPixelsCam from " << fname << endl;
199 return kFALSE;
200 }
201 fBadPixels.Merge(bad);
202 }
203
204 if (fDisplay && !fDisplay->GetCanvas("Pedestals"))
205 fDisplay->Read();
206
207 return kTRUE;
208}
209
210MExtractor *MJPedestal::ReadCalibration() const
211{
212 const TString fname = Form("%s/calib%08d.root",fPathIn.Data(), fSequence.GetSequence());
213
214 *fLog << inf << "Reading extractor from file: " << fname << endl;
215
216 TFile file(fname, "READ");
217 if (!file.IsOpen())
218 {
219 *fLog << err << dbginf << "ERROR - Could not open file " << fname << endl;
220 return NULL;
221 }
222
223 TObject *o = file.Get("ExtractSignal");
224 if (o && !o->InheritsFrom(MExtractor::Class()))
225 {
226 *fLog << err << dbginf << "ERROR - ExtractSignal read from " << fname << " doesn't inherit from MExtractor!" << endl;
227 return NULL;
228 }
229
230 return o ? (MExtractor*)o->Clone("ExtractSignal") : NULL;
231}
232
233//---------------------------------------------------------------------------------
234//
235// Display the results.
236// If Display type "kDataCheck" was chosen, also the reference lines are displayed.
237//
238void MJPedestal::DisplayResult(MParList &plist)
239{
240 if (!fDisplay)
241 return;
242
243 //
244 // Update display
245 //
246 TString title = fDisplay->GetTitle();
247 title += "-- Pedestal ";
248 if (fSequence.IsValid())
249 title += fSequence.GetName();
250 else
251 if (fRuns) // FIXME: What to do if an environmentfile was used?
252 title += fRuns->GetRunsAsString();
253 title += " --";
254 fDisplay->SetTitle(title);
255
256 //
257 // Get container from list
258 //
259 MGeomCam &geomcam = *(MGeomCam*)plist.FindObject("MGeomCam");
260 MCalibrationPedCam &calpedcam = *(MCalibrationPedCam*)plist.FindObject("MCalibrationPedCam");
261
262 //
263 // Create container to display
264 //
265 MHCamera disp0 (geomcam, "MPedestalCam;ped", "Mean Pedestal");
266 MHCamera disp1 (geomcam, "MPedestalCam;RMS", "Pedestal RMS");
267 MHCamera disp2 (geomcam, "MCalibPedCam;histmean", "Mean Pedestal (Hist.)");
268 MHCamera disp3 (geomcam, "MCalibPedCam;histsigma", "Pedestal RMS (Hist.)");
269 MHCamera disp4 (geomcam, "MCalibPedCam;ped", "Mean Pedestal");
270 MHCamera disp5 (geomcam, "MCalibPedCam;RMS", "Pedestal RMS");
271 MHCamera disp6 (geomcam, "MCalibDiffCam;ped", "Diff. Mean Pedestal (Hist.)");
272 MHCamera disp7 (geomcam, "MCalibDiffCam;RMS", "Diff. Pedestal RMS (Hist.)");
273 MHCamera disp8 (geomcam, "MCalibDiffCam;ped", "Diff. Mean Pedestal");
274 MHCamera disp9 (geomcam, "MCalibDiffCam;AbsRMS", "Diff. Abs. Pedestal RMS");
275 MHCamera disp10(geomcam, "MCalibDiffCam;RelRMS", "Diff. Rel. Pedestal RMS");
276
277 disp0.SetCamContent(fPedestalCamOut, 0);
278 disp0.SetCamError (fPedestalCamOut, 1);
279
280 disp1.SetCamContent(fPedestalCamOut, 2);
281 disp1.SetCamError (fPedestalCamOut, 3);
282
283 if (fExtractType == kUseHists)
284 {
285 disp2.SetCamContent(calpedcam, 0);
286 disp2.SetCamError (calpedcam, 1);
287
288 disp3.SetCamContent(calpedcam, 2);
289 disp3.SetCamError (calpedcam, 3);
290
291 disp4.SetCamContent(calpedcam, 5);
292 disp4.SetCamError (calpedcam, 6);
293
294 disp5.SetCamContent(calpedcam, 7);
295 disp5.SetCamError (calpedcam, 8);
296
297 for (UInt_t i=0;i<geomcam.GetNumPixels();i++)
298 {
299
300 MPedestalPix &ped = fPedestalCamOut[i];
301 MCalibrationPix &hist = calpedcam [i];
302 MBadPixelsPix &bad = fBadPixels[i];
303
304 if (bad.IsUnsuitable())
305 continue;
306
307 disp6.Fill(i,ped.GetPedestal()-hist.GetHiGainMean());
308 disp6.SetUsed(i);
309
310 disp7.Fill(i,hist.GetHiGainSigma()-ped.GetPedestalRms());
311 if (TMath::Abs(ped.GetPedestalRms()-hist.GetHiGainSigma()) < 4.0)
312 disp7.SetUsed(i);
313
314 disp8.Fill(i,ped.GetPedestal()-hist.GetLoGainMean());
315 disp8.SetUsed(i);
316
317 disp9.Fill(i,hist.GetLoGainSigma()-ped.GetPedestalRms());
318 if (TMath::Abs(hist.GetLoGainSigma() - ped.GetPedestalRms()) < 4.0)
319 disp9.SetUsed(i);
320 }
321 }
322
323 if (fExtractionType!=kFundamental/*fExtractorResolution*/)
324 {
325 for (UInt_t i=0;i<geomcam.GetNumPixels();i++)
326 {
327
328 MPedestalPix &pedo = fPedestalCamOut[i];
329 MPedestalPix &pedi = fPedestalCamIn[i];
330 MBadPixelsPix &bad = fBadPixels[i];
331
332 if (bad.IsUnsuitable())
333 continue;
334
335 const Float_t diff = pedo.GetPedestalRms()-pedi.GetPedestalRms();
336 const Float_t sum = 0.5*(pedo.GetPedestalRms()+pedi.GetPedestalRms());
337
338 disp9.Fill(i,pedo.GetPedestalRms()-pedi.GetPedestalRms());
339 if (pedo.IsValid() && pedi.IsValid())
340 disp9.SetUsed(i);
341
342 disp10.Fill(i,sum == 0. ? 0. : diff/sum);
343 if (pedo.IsValid() && pedi.IsValid() && sum != 0.)
344 disp10.SetUsed(i);
345 }
346 }
347
348 disp0.SetYTitle("P [cts/slice]");
349 disp1.SetYTitle("P_{rms} [cts/slice]");
350 disp2.SetYTitle("Hist. Mean [cts/slice]");
351 disp3.SetYTitle("Hist. Sigma [cts/slice]");
352 disp4.SetYTitle("Calc. Mean [cts/slice]");
353 disp5.SetYTitle("Calc. RMS [cts/slice]");
354 disp6.SetYTitle("Diff. Mean [cts/slice]");
355 disp7.SetYTitle("Diff. RMS [cts/slice]");
356 disp8.SetYTitle("Diff. Mean [cts/slice]");
357 disp9.SetYTitle("Abs.Diff.RMS [cts/slice]");
358 disp10.SetYTitle("Rel.Diff.RMS [1]");
359
360 //
361 // Display data
362 //
363 if (fDisplayType != kDisplayDataCheck && fExtractType != kUseHists && fExtractionType==kFundamental/*fExtractorResolution*/)
364 {
365 TCanvas &c3 = fDisplay->AddTab("Pedestals");
366 c3.Divide(2,3);
367
368 disp0.CamDraw(c3, 1, 2, 1);
369 disp1.CamDraw(c3, 2, 2, 6);
370 return;
371 }
372
373 if (fExtractType == kUseHists)
374 {
375
376 TCanvas &c3 = fDisplay->AddTab("Extractor Hist.");
377 c3.Divide(2,3);
378
379 disp2.CamDraw(c3, 1, 2, 1);
380 disp3.CamDraw(c3, 2, 2, 5);
381
382 TCanvas &c4 = fDisplay->AddTab("Extractor Calc.");
383 c4.Divide(2,3);
384
385 disp4.CamDraw(c4, 1, 2, 1);
386 disp5.CamDraw(c4, 2, 2, 5);
387
388
389 TCanvas &c5 = fDisplay->AddTab("Difference Hist.");
390 c5.Divide(2,3);
391
392 disp6.CamDraw(c5, 1, 2, 1);
393 disp7.CamDraw(c5, 2, 2, 5);
394
395 TCanvas &c6 = fDisplay->AddTab("Difference Calc.");
396 c6.Divide(2,3);
397
398 disp8.CamDraw(c6, 1, 2, 1);
399 disp9.CamDraw(c6, 2, 2, 5);
400 return;
401 }
402
403 if (fDisplayType == kDisplayDataCheck)
404 {
405
406 TCanvas &c3 = fDisplay->AddTab(fExtractionType!=kFundamental/*fExtractorResolution*/ ? "PedExtrd" : "Ped");
407 c3.Divide(2,3);
408
409 c3.cd(1);
410 gPad->SetBorderMode(0);
411 gPad->SetTicks();
412 MHCamera *obj1=(MHCamera*)disp0.DrawCopy("hist");
413 //
414 // for the datacheck, fix the ranges!!
415 //
416 if (fExtractionType==kFundamental/*!fExtractorResolution*/)
417 {
418 obj1->SetMinimum(fPedestalMin);
419 obj1->SetMaximum(fPedestalMax);
420 }
421 //
422 // Set the datacheck sizes:
423 //
424 FixDataCheckHist((TH1D*)obj1);
425 //
426 // set reference lines
427 //
428 DisplayReferenceLines(obj1,0);
429 //
430 // end reference lines
431 //
432 c3.cd(3);
433 gPad->SetBorderMode(0);
434 obj1->SetPrettyPalette();
435 obj1->Draw();
436
437 c3.cd(5);
438 gPad->SetBorderMode(0);
439 gPad->SetTicks();
440 TH1D *obj2 = (TH1D*)obj1->Projection(obj1->GetName());
441 obj2->Draw();
442 obj2->SetBit(kCanDelete);
443 obj2->Fit("gaus","Q");
444 obj2->GetFunction("gaus")->SetLineColor(kYellow);
445 //
446 // Set the datacheck sizes:
447 //
448 FixDataCheckHist(obj2);
449 obj2->SetStats(1);
450
451 c3.cd(2);
452 gPad->SetBorderMode(0);
453 gPad->SetTicks();
454 MHCamera *obj3=(MHCamera*)disp1.DrawCopy("hist");
455 //
456 // for the datacheck, fix the ranges!!
457 //
458 obj3->SetMinimum(fPedRmsMin);
459 obj3->SetMaximum(fPedRmsMax);
460 //
461 // Set the datacheck sizes:
462 //
463 FixDataCheckHist((TH1D*)obj3);
464 //
465 // set reference lines
466 //
467 DisplayReferenceLines(obj3,1);
468
469 c3.cd(4);
470 gPad->SetBorderMode(0);
471 obj3->SetPrettyPalette();
472 obj3->Draw();
473
474 c3.cd(6);
475 gPad->SetBorderMode(0);
476
477 if (geomcam.InheritsFrom("MGeomCamMagic"))
478 {
479 TArrayI inner(1);
480 inner[0] = 0;
481
482 TArrayI outer(1);
483 outer[0] = 1;
484
485 TArrayI s0(6);
486 s0[0] = 6;
487 s0[1] = 1;
488 s0[2] = 2;
489 s0[3] = 3;
490 s0[4] = 4;
491 s0[5] = 5;
492
493 TArrayI s1(3);
494 s1[0] = 6;
495 s1[1] = 1;
496 s1[2] = 2;
497
498 TArrayI s2(3);
499 s2[0] = 3;
500 s2[1] = 4;
501 s2[2] = 5;
502
503 TVirtualPad *pad = gPad;
504 pad->Divide(2,1);
505
506 TH1D *inout[2];
507 inout[0] = disp1.ProjectionS(s0, inner, "Inner");
508 inout[1] = disp1.ProjectionS(s0, outer, "Outer");
509 FixDataCheckHist(inout[0]);
510 FixDataCheckHist(inout[1]);
511
512 inout[0]->SetTitle(Form("%s %s",disp1.GetTitle(),"Inner"));
513 inout[1]->SetTitle(Form("%s %s",disp1.GetTitle(),"Outer"));
514
515
516 for (int i=0; i<2; i++)
517 {
518 pad->cd(i+1);
519 gPad->SetBorderMode(0);
520 gPad->SetTicks();
521
522 inout[i]->SetDirectory(NULL);
523 inout[i]->SetLineColor(kRed+i);
524 inout[i]->SetBit(kCanDelete);
525 inout[i]->Draw();
526 inout[i]->Fit("gaus", "Q");
527
528 TLegend *leg2 = new TLegend(0.6,0.2,0.9,0.55);
529 leg2->SetHeader(inout[i]->GetName());
530 leg2->AddEntry(inout[i], inout[i]->GetName(), "l");
531
532 //
533 // Display the outliers as dead and noisy pixels
534 //
535 DisplayOutliers(inout[i]);
536
537 //
538 // Display the two half of the camera separately
539 //
540 TH1D *half[2];
541 half[0] = disp1.ProjectionS(s1, i==0 ? inner : outer , "Sector 6-1-2");
542 half[1] = disp1.ProjectionS(s2, i==0 ? inner : outer , "Sector 3-4-5");
543
544 for (int j=0; j<2; j++)
545 {
546 half[j]->SetLineColor(kRed+i+2*j+1);
547 half[j]->SetDirectory(NULL);
548 half[j]->SetBit(kCanDelete);
549 half[j]->Draw("same");
550 leg2->AddEntry(half[j], half[j]->GetName(), "l");
551 }
552 leg2->Draw();
553 delete leg2;
554 }
555 return;
556 }
557 }
558
559 if (fExtractionType!=kFundamental/*fExtractorResolution*/)
560 {
561
562 TCanvas &c3 = fDisplay->AddTab(fExtractionType==kWithExtractor?"PedExtrd":"PedRndm");
563 c3.Divide(2,3);
564
565 disp0.CamDraw(c3, 1, 2, 1);
566 disp1.CamDraw(c3, 2, 2, 6);
567
568 TCanvas &c13 = fDisplay->AddTab(fExtractionType==kWithExtractor?"DiffExtrd":"DiffRndm");
569 c13.Divide(2,3);
570
571 disp9.CamDraw(c13, 1, 2, 5);
572 disp10.CamDraw(c13, 2, 2, 5);
573 return;
574 }
575}
576
577
578void MJPedestal::DisplayReferenceLines(MHCamera *cam, const Int_t what) const
579{
580
581 Double_t x = cam->GetNbinsX();
582
583 const MGeomCam *geom = cam->GetGeometry();
584
585 if (geom->InheritsFrom("MGeomCamMagic"))
586 x = what ? 397 : cam->GetNbinsX();
587
588 TLine line;
589 line.SetLineStyle(kDashed);
590 line.SetLineWidth(3);
591 line.SetLineColor(kBlue);
592
593 TLegend *leg = new TLegend(0.6,0.75,0.9,0.99);
594 leg->SetBit(kCanDelete);
595
596 if (fExtractionType==kWithExtractorRndm && !(what))
597 {
598 TLine *l0 = line.DrawLine(0,0.,cam->GetNbinsX(),0.);
599 l0->SetBit(kCanDelete);
600 leg->AddEntry(l0, "Reference","l");
601 leg->Draw();
602 return;
603 }
604
605 line.SetLineColor(kBlue);
606 TLine *l1 = line.DrawLine(0, what ? fRefPedRmsGalacticInner : fRefPedGalactic,
607 x, what ? fRefPedRmsGalacticInner : fRefPedGalactic);
608 l1->SetBit(kCanDelete);
609 line.SetLineColor(kYellow);
610 TLine *l2 = line.DrawLine(0, what ? fRefPedRmsExtraGalacticInner : fRefPedExtraGalactic,
611 x, what ? fRefPedRmsExtraGalacticInner : fRefPedExtraGalactic);
612 l2->SetBit(kCanDelete);
613 line.SetLineColor(kMagenta);
614 TLine *l3 = line.DrawLine(0, what ? fRefPedRmsClosedLidsInner : fRefPedClosedLids,
615 x, what ? fRefPedRmsClosedLidsInner : fRefPedClosedLids);
616 l3->SetBit(kCanDelete);
617
618 if (geom->InheritsFrom("MGeomCamMagic"))
619 if (what)
620 {
621 const Double_t x2 = cam->GetNbinsX();
622
623 line.SetLineColor(kBlue);
624 line.DrawLine(398, fRefPedRmsGalacticOuter,
625 x2, fRefPedRmsGalacticOuter);
626
627 line.SetLineColor(kYellow);
628 line.DrawLine(398, fRefPedRmsExtraGalacticOuter,
629 x2, fRefPedRmsExtraGalacticOuter);
630
631 line.SetLineColor(kMagenta);
632 line.DrawLine(398, fRefPedRmsClosedLidsOuter,
633 x2, fRefPedRmsClosedLidsOuter);
634 }
635
636
637 leg->AddEntry(l1, "Galactic Source","l");
638 leg->AddEntry(l2, "Extra-Galactic Source","l");
639 leg->AddEntry(l3, "Closed Lids","l");
640 leg->Draw();
641}
642
643void MJPedestal::DisplayOutliers(TH1D *hist) const
644{
645 const Float_t mean = hist->GetFunction("gaus")->GetParameter(1);
646 const Float_t lolim = mean - 3.5*hist->GetFunction("gaus")->GetParameter(2);
647 const Float_t uplim = mean + 3.5*hist->GetFunction("gaus")->GetParameter(2);
648 const Stat_t dead = hist->Integral(0,hist->FindBin(lolim)-1);
649 const Stat_t noisy = hist->Integral(hist->FindBin(uplim)+1,hist->GetNbinsX()+1);
650
651 TLatex deadtex;
652 deadtex.SetTextSize(0.06);
653 deadtex.DrawLatex(0.1,hist->GetBinContent(hist->GetMaximumBin())/1.1,Form("%3i dead pixels",(Int_t)dead));
654
655 TLatex noisytex;
656 noisytex.SetTextSize(0.06);
657 noisytex.DrawLatex(0.1,hist->GetBinContent(hist->GetMaximumBin())/1.2,Form("%3i noisy pixels",(Int_t)noisy));
658}
659
660void MJPedestal::FixDataCheckHist(TH1D *hist) const
661{
662 hist->SetDirectory(NULL);
663 hist->SetStats(0);
664
665 //
666 // set the labels bigger
667 //
668 TAxis *xaxe = hist->GetXaxis();
669 TAxis *yaxe = hist->GetYaxis();
670
671 xaxe->CenterTitle();
672 yaxe->CenterTitle();
673 xaxe->SetTitleSize(0.06);
674 yaxe->SetTitleSize(0.06);
675 xaxe->SetTitleOffset(0.8);
676 yaxe->SetTitleOffset(0.5);
677 xaxe->SetLabelSize(0.05);
678 yaxe->SetLabelSize(0.05);
679}
680
681/*
682Bool_t MJPedestal::WriteEventloop(MEvtLoop &evtloop) const
683{
684 if (fOutputPath.IsNull())
685 return kTRUE;
686
687 const TString oname(GetOutputFile());
688
689 *fLog << inf << "Writing to file: " << oname << endl;
690
691 TFile file(oname, fOverwrite?"RECREATE":"NEW", "File created by MJPedestal", 9);
692 if (!file.IsOpen())
693 {
694 *fLog << err << "ERROR - Couldn't open file " << oname << " for writing..." << endl;
695 return kFALSE;
696 }
697
698 if (evtloop.Write(fName)<=0)
699 {
700 *fLog << err << "Unable to write MEvtloop to " << oname << endl;
701 return kFALSE;
702 }
703
704 return kTRUE;
705}
706*/
707
708void MJPedestal::SetExtractor(MExtractor* ext)
709{
710 if (ext)
711 {
712 if (fExtractor)
713 delete fExtractor;
714 fExtractor = ext ? (MExtractor*)ext->Clone(ext->GetName()) : NULL;
715 }
716 else
717 fExtractor = 0;
718}
719
720// --------------------------------------------------------------------------
721//
722// Read the following values from resource file:
723//
724// PedestalMin
725// PedestalMax
726//
727// PedRmsMin
728// PedRmsMax
729//
730// RefPedClosedLids
731// RefPedExtraGalactic
732// RefPedGalactic
733//
734// RefPedRmsClosedLidsInner
735// RefPedRmsExtraGalacticInner
736// RefPedRmsGalacticInner
737// RefPedRmsClosedLidsOuter
738// RefPedRmsExtraGalacticOuter
739// RefPedRmsGalacticOuter
740//
741void MJPedestal::ReadReferenceFile()
742{
743 TEnv refenv(fReferenceFile);
744
745 fPedestalMin = refenv.GetValue("PedestalMin",fPedestalMin);
746 fPedestalMax = refenv.GetValue("PedestalMax",fPedestalMax);
747 fPedRmsMin = refenv.GetValue("PedRmsMin",fPedRmsMin);
748 fPedRmsMax = refenv.GetValue("PedRmsMax",fPedRmsMax);
749 fRefPedClosedLids = refenv.GetValue("RefPedClosedLids",fRefPedClosedLids);
750 fRefPedExtraGalactic = refenv.GetValue("RefPedExtraGalactic",fRefPedExtraGalactic);
751 fRefPedGalactic = refenv.GetValue("RefPedGalactic",fRefPedGalactic);
752 fRefPedRmsClosedLidsInner = refenv.GetValue("RefPedRmsClosedLidsInner",fRefPedRmsClosedLidsInner);
753 fRefPedRmsExtraGalacticInner = refenv.GetValue("RefPedRmsExtraGalacticInner",fRefPedRmsExtraGalacticInner);
754 fRefPedRmsGalacticInner = refenv.GetValue("RefPedRmsGalacticInner",fRefPedRmsGalacticInner);
755 fRefPedRmsClosedLidsOuter = refenv.GetValue("RefPedRmsClosedLidsOuter",fRefPedRmsClosedLidsOuter);
756 fRefPedRmsExtraGalacticOuter = refenv.GetValue("RefPedRmsExtraGalacticOuter",fRefPedRmsExtraGalacticOuter);
757 fRefPedRmsGalacticOuter = refenv.GetValue("RefPedRmsGalacticOuter",fRefPedRmsGalacticOuter);
758}
759
760// --------------------------------------------------------------------------
761//
762// The following resource options are available:
763//
764// Do a datacheck run (read raw-data and enable display)
765// Prefix.DataCheck: Yes, No <default>
766//
767// Setup display type
768// Prefix.Display: normal <default>, datacheck, none
769//
770// Use cosmic data instead of pedestal data (DatRuns)
771// Prefix.UseData: Yes, No <default>
772//
773// Write an output file with pedestals and status-display
774// Prefix.DisableOutput: Yes, No <default>
775//
776// Name of a file containing reference values (see ReadReferenceFile)
777// Prefix.ReferenceFile: filename
778// (see ReadReferenceFile)
779//
780Bool_t MJPedestal::CheckEnvLocal()
781{
782 if (HasEnv("Display"))
783 {
784 TString type = GetEnv("Display", "normal");
785 type.ToLower();
786 if (type==(TString)"normal")
787 fDisplayType = kDisplayNormal;
788 if (type==(TString)"datacheck")
789 fDisplayType = kDisplayDataCheck;
790 if (type==(TString)"none")
791 fDisplayType = kDisplayNone;
792 }
793
794 if (!MJCalib::CheckEnvLocal())
795 return kFALSE;
796
797 if (HasEnv("UseData"))
798 fExtractType = GetEnv("UseData",kFALSE) ? kUseData : kUsePedRun;
799
800 if (HasEnv("UseHists"))
801 if (GetEnv("UseHists",kFALSE))
802 fExtractType = kUseHists;
803
804 SetNoStorage(GetEnv("DisableOutput", IsNoStorage()));
805
806 MTaskEnv tenv("ExtractSignal");
807 tenv.SetDefault(fExtractor);
808
809 if (tenv.ReadEnv(*GetEnv(), GetEnvPrefix()+".ExtractSignal", GetEnvDebug())==kERROR)
810 return kFALSE;
811
812 if (fExtractor==tenv.GetTask())
813 return kTRUE;
814
815 if (!tenv.GetTask()->InheritsFrom(MExtractor::Class()))
816 {
817 *fLog << err << "ERROR: ExtractSignal from resource file doesn't inherit from MExtractor.... abort." << endl;
818 return kFALSE;
819 }
820
821 SetExtractor((MExtractor*)tenv.GetTask());
822
823 fBadPixelsFile = GetEnv("BadPixelsFile",fBadPixelsFile.Data());
824 fReferenceFile = GetEnv("ReferenceFile",fReferenceFile.Data());
825 ReadReferenceFile();
826
827 return kTRUE;
828}
829
830//---------------------------------------------------------------------------------
831//
832// Try to write the created MPedestalCam in the output file.
833// If possible, also an MBadPixelsCam and the corresponding display is written.
834//
835// In case of Storage type "kNoStorage" or if any of the containers
836// cannot be written, return kFALSE, otherwise kTRUE.
837//
838Bool_t MJPedestal::WriteResult()
839{
840 if (IsNoStorage())
841 return kTRUE;
842
843 TObjArray cont;
844
845 cont.Add(&fPedestalCamOut);
846 cont.Add(&fBadPixels);
847
848 return WriteContainer(cont, GetOutputFileName(), fOverwrite?"RECREATE":"NEW");
849}
850
851Bool_t MJPedestal::Process()
852{
853 if (!ReadPedestalCam())
854 return ProcessFile();
855
856 return kTRUE;
857}
858
859Bool_t MJPedestal::ProcessFile()
860{
861 if (!fSequence.IsValid())
862 {
863 if (!fRuns)
864 {
865 *fLog << err << "Neither AddRuns nor SetSequence nor SetEnv was called... abort." << endl;
866 return kFALSE;
867 }
868 if (fRuns && fRuns->GetNumRuns() != fRuns->GetNumEntries())
869 {
870 *fLog << err << "Number of files found doesn't match number of runs... abort." << endl;
871 return kFALSE;
872 }
873 }
874
875 if (!CheckEnv())
876 return kFALSE;
877
878 // --------------------------------------------------------------------------------
879
880 const TString type = IsUseData() ? "data" : "pedestal";
881
882 *fLog << inf;
883 fLog->Separator(GetDescriptor());
884 *fLog << "Calculate MPedestalCam from " << type << "-runs ";
885 if (fSequence.IsValid())
886 *fLog << fSequence.GetName() << endl;
887 else
888 if (fRuns)
889 *fLog << fRuns->GetRunsAsString() << endl;
890 else
891 *fLog << "in Resource File." << endl;
892 *fLog << endl;
893
894 // --------------------------------------------------------------------------------
895
896 MParList plist;
897 MTaskList tlist;
898 plist.AddToList(&tlist);
899 plist.AddToList(this); // take care of fDisplay!
900
901 MReadMarsFile read("Events");
902 MRawFileRead rawread(NULL);
903
904 MDirIter iter;
905 if (fSequence.IsValid())
906 {
907 const Int_t n0 = IsUseData()
908 ? fSequence.SetupDatRuns(iter, fPathData, "D", IsUseRawData())
909 : fSequence.SetupPedRuns(iter, fPathData, "P", IsUseRawData());
910 const Int_t n1 = IsUseData()
911 ? fSequence.GetNumDatRuns()
912 : fSequence.GetNumPedRuns();
913 if (n0==0)
914 {
915 *fLog << err << "ERROR - No " << type << " input files of sequence found in " << (fPathData.IsNull()?"<default>":fPathData.Data()) << endl;
916 return kFALSE;
917 }
918 if (n0!=n1)
919 {
920 *fLog << err << "ERROR - Number of " << type << " files found (" << n0 << ") in " << (fPathData.IsNull()?"<default>":fPathData.Data()) << " doesn't match number of files in sequence (" << n1 << ")" << endl;
921 return kFALSE;
922 }
923 }
924
925 if (IsUseRawData())
926 {
927 if (fRuns || fSequence.IsValid())
928 rawread.AddFiles(fSequence.IsValid() ? iter : *fRuns);
929 tlist.AddToList(&rawread);
930 }
931 else
932 {
933 read.DisableAutoScheme();
934 if (fRuns || fSequence.IsValid())
935 read.AddFiles(fSequence.IsValid() ? iter : *fRuns);
936 tlist.AddToList(&read);
937 }
938
939 // Setup Tasklist
940 plist.AddToList(&fPedestalCamOut);
941 plist.AddToList(&fBadPixels);
942
943 //
944 // Read bad pixels from outside
945 //
946 if (!fBadPixelsFile.IsNull())
947 {
948 *fLog << err << "Excluding: " << fBadPixelsFile << endl;
949 ifstream f(fBadPixelsFile.Data());
950 fBadPixels.AsciiRead((istream&)f);
951 f.close();
952 }
953
954 MGeomApply geomapl;
955 MBadPixelsMerge merge(&fBadPixels);
956
957 MPedCalcPedRun pedcalc;
958 pedcalc.SetPedestalUpdate(kFALSE);
959
960 MPedCalcFromLoGain pedlogain;
961 pedlogain.SetPedestalUpdate(kFALSE);
962
963 fPedestalHist.SetRenorm(kTRUE);
964
965 MFillH fillped(&fPedestalHist, "MExtractedSignalCam", "FillPedestalCam");
966 fillped.SetBit(MFillH::kDoNotDisplay);
967
968 tlist.AddToList(&geomapl);
969 tlist.AddToList(&merge);
970
971 if (!fPathIn.IsNull())
972 {
973 fExtractor = ReadCalibration();
974 if (!fExtractor)
975 return kFALSE;
976
977 *fLog << all;
978 *fLog << underline << "Signal Extractor found in calibration file:" << endl;
979 fExtractor->Print();
980 *fLog << endl;
981 }
982
983 MFEventSelector sel;
984 sel.SetNumSelectEvts(0);
985
986 if (fExtractionType!=kFundamental && fExtractType!=kUseHists)
987 {
988 fExtractor->SetFilter(&sel);
989 tlist.AddToList(&sel);
990 tlist.AddToList(fExtractor);
991 }
992
993 // ----------------------------------------------------------------------
994 // Now we make sure, that in all cases the ranges are setup correctly
995 // ----------------------------------------------------------------------
996 MTaskEnv taskenv("ExtractPedestal");
997 switch (fExtractType)
998 {
999 case kUsePedRun:
1000 // In case other than 'fundamental' second argument is obsolete
1001 // pedcalc.SetExtractWindow(0,14); // kUsePedRun (take default from class)
1002 taskenv.SetDefault(&pedcalc);
1003 tlist.AddToList(&taskenv);
1004 break;
1005
1006 case kUseData:
1007 // In case other than 'fundamental' second argument is obsolete
1008 // pedlogain.SetExtractWindow(15,14); // kUseData (take default from class)
1009 taskenv.SetDefault(&pedlogain);
1010 tlist.AddToList(&taskenv);
1011 break;
1012
1013 case kUseHists:
1014 if (!fExtractor)
1015 {
1016 *fLog << err << GetDescriptor() << " - ERROR: ";
1017 *fLog << "Extraction Type is kUseHists, but no extractor was set" << endl;
1018 return kFALSE;
1019 }
1020
1021 fExtractor->SetPedestals(&fPedestalCamIn);
1022 plist.AddToList(&fPedestalHist);
1023 tlist.AddToList(fExtractor);
1024 tlist.AddToList(&fillped);
1025 break;
1026 }
1027
1028 pedcalc.SetPedestalsIn(&fPedestalCamIn);
1029 pedlogain.SetPedestalsIn(&fPedestalCamIn);
1030
1031 // kFundamental
1032 if (fExtractor)
1033 {
1034
1035 if (fExtractor->InheritsFrom("MExtractTimeAndCharge"))
1036 {
1037
1038 const Float_t f = 0.5+fExtractor->GetHiGainFirst();
1039 const Int_t win = ((MExtractTimeAndCharge*)fExtractor)->GetWindowSizeHiGain();
1040 pedcalc.SetExtractWindow((Int_t)f, win);
1041 pedlogain.SetExtractWindow((Int_t)(15+f), win);
1042
1043 }
1044 else
1045 {
1046 const Float_t f = 0.5+fExtractor->GetHiGainFirst();
1047 const Float_t n = 0.5+fExtractor->GetNumHiGainSamples();
1048 pedcalc.SetExtractWindow((Int_t)f, (Int_t)n);
1049 pedlogain.SetExtractWindow((Int_t)(15+f), (Int_t)n);
1050
1051 if (fExtractionType!=kFundamental)
1052 {
1053 *fLog << inf;
1054 *fLog << "Signal extractor doesn't inherit from MExtractTimeAndCharge..." << endl;
1055 *fLog << " --> falling back to fundamental pedestal extraction." << endl;
1056 fExtractionType=kFundamental;
1057 }
1058 }
1059
1060
1061 if (fExtractionType!=kFundamental)
1062 {
1063 pedcalc.SetRandomCalculation(fExtractionType==kWithExtractorRndm);
1064 pedlogain.SetRandomCalculation(fExtractionType==kWithExtractorRndm);
1065
1066 pedcalc.SetExtractor((MExtractTimeAndCharge*)fExtractor);
1067 pedlogain.SetExtractor((MExtractTimeAndCharge*)fExtractor);
1068 }
1069 }
1070 else
1071 {
1072 *fLog << warn << GetDescriptor() << ": WARNING - No extractor has been handed over! " << endl;
1073 *fLog << "Taking default window for pedestal extraction. The calculated pedestal RMS" << endl;
1074 *fLog << "will probably not match with future pedestal RMS' from different extraction" << endl;
1075 *fLog << "windows." << endl;
1076 }
1077
1078 /*
1079 switch (fExtractType)
1080 {
1081 case kUseData:
1082 *fLog << all << "TYPE: USEDATA " << fExtractor << endl;
1083 taskenv.SetDefault(&pedlogain);
1084 tlist.AddToList(&taskenv);
1085
1086 if (!SetupExtractor(plist, pedlogain))
1087 {
1088 *fLog << all << "SETTING TO: " << fExtractor << " " << fExtractor->GetNumHiGainSamples() << endl;
1089 fExtractor->Print();
1090 pedlogain.SetExtractWindow(15, (Int_t)TMath::Nint(fExtractor->GetNumHiGainSamples()));
1091 }
1092 break;
1093
1094 case kUsePedRun:
1095 *fLog << all << "TYPE: USEPEDRUN " << fExtractor << endl;
1096 taskenv.SetDefault(&pedcalc);
1097 tlist.AddToList(&taskenv);
1098
1099 if (!SetupExtractor(plist, pedcalc))
1100 pedcalc.SetExtractWindow(fExtractor->GetHiGainFirst(), TMath::Nint(fExtractor->GetNumHiGainSamples()));
1101 break;
1102
1103 case kUseHists:
1104 if (!fExtractor)
1105 {
1106 *fLog << err << GetDescriptor() << " - ERROR: ";
1107 *fLog << "Extraction Type is kUseHists, but no extractor was set" << endl;
1108 return kFALSE;
1109 }
1110
1111 tlist.AddToList(fExtractor);
1112 tlist.AddToList(&fillped);
1113 break;
1114 } */
1115
1116 /*
1117 if (!fPathIn.IsNull())
1118 {
1119 delete fExtractor;
1120 fExtractor = 0;
1121 }
1122 */
1123 //
1124 // Execute the eventloop
1125 //
1126 MEvtLoop evtloop(fName);
1127 evtloop.SetParList(&plist);
1128 evtloop.SetDisplay(fDisplay);
1129 evtloop.SetLogStream(fLog);
1130 if (!SetupEnv(evtloop))
1131 return kFALSE;
1132
1133 // if (!WriteEventloop(evtloop))
1134 // return kFALSE;
1135
1136 // Execute first analysis
1137 if (!evtloop.Eventloop(fMaxEvents))
1138 {
1139 *fLog << err << GetDescriptor() << ": Failed." << endl;
1140 return kFALSE;
1141 }
1142
1143 tlist.PrintStatistics();
1144
1145 if (fDisplayType!=kDisplayNone)
1146 DisplayResult(plist);
1147
1148 if (fExtractionType!=kFundamental)
1149 fExtractor->SetFilter(0);
1150
1151 if (!WriteResult())
1152 return kFALSE;
1153
1154 *fLog << all << GetDescriptor() << ": Done." << endl;
1155 *fLog << endl << endl;
1156
1157 return kTRUE;
1158}
Note: See TracBrowser for help on using the repository browser.