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

Last change on this file since 4486 was 4368, checked in by tbretz, 21 years ago
*** empty log message ***
File size: 13.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 "MRunIter.h"
47#include "MParList.h"
48#include "MTaskList.h"
49#include "MEvtLoop.h"
50#include "MExtractor.h"
51
52#include "MStatusDisplay.h"
53
54#include "MGeomCam.h"
55#include "MHCamera.h"
56#include "MPedestalCam.h"
57#include "MBadPixelsCam.h"
58
59#include "MReadMarsFile.h"
60#include "MRawFileRead.h"
61#include "MGeomApply.h"
62#include "MBadPixelsMerge.h"
63#include "MPedCalcPedRun.h"
64
65ClassImp(MJPedestal);
66
67using namespace std;
68
69const Double_t MJPedestal::fgPedestalMin = 4.;
70const Double_t MJPedestal::fgPedestalMax = 16.;
71const Double_t MJPedestal::fgPedRmsMin = 0.;
72const Double_t MJPedestal::fgPedRmsMax = 20.;
73
74const Float_t MJPedestal::fgRefPedClosedLids = 9.635;
75const Float_t MJPedestal::fgRefPedExtraGalactic = 9.93;
76const Float_t MJPedestal::fgRefPedGalactic = 10.03;
77const Float_t MJPedestal::fgRefPedRmsClosedLids = 1.7;
78const Float_t MJPedestal::fgRefPedRmsExtraGalactic = 5.6;
79const Float_t MJPedestal::fgRefPedRmsGalactic = 6.92;
80
81// --------------------------------------------------------------------------
82//
83// Default constructor.
84//
85// Sets fRuns to 0, fExtractor to NULL, fDataCheck to kFALSE
86//
87MJPedestal::MJPedestal(const char *name, const char *title)
88 : fEnv(0), fRuns(0), fExtractor(NULL), fDataCheck(kFALSE)
89{
90 fName = name ? name : "MJPedestal";
91 fTitle = title ? title : "Tool to create a pedestal file (MPedestalCam)";
92}
93
94MJPedestal::~MJPedestal()
95{
96 if (fEnv)
97 delete fEnv;
98}
99
100const char* MJPedestal::GetOutputFile() const
101{
102 if (!fRuns)
103 return "";
104
105 return Form("%s/%s-F0.root", (const char*)fOutputPath, (const char*)fRuns->GetRunsAsFileName());
106}
107
108Bool_t MJPedestal::ReadPedestalCam()
109{
110 const TString fname = GetOutputFile();
111
112 if (gSystem->AccessPathName(fname, kFileExists))
113 {
114 *fLog << warn << "Input file " << fname << " doesn't exist, will create it." << endl;
115 return kFALSE;
116 }
117
118 *fLog << inf << "Reading from file: " << fname << endl;
119
120 TFile file(fname, "READ");
121 if (fPedestalCam.Read()<=0)
122 {
123 *fLog << err << "Unable to read MPedestalCam from " << fname << endl;
124 return kFALSE;
125 }
126
127 if (file.FindKey("MBadPixelsCam"))
128 {
129 MBadPixelsCam bad;
130 if (bad.Read()<=0)
131 {
132 *fLog << err << "Unable to read MBadPixelsCam from " << fname << endl;
133 return kFALSE;
134 }
135 fBadPixels.Merge(bad);
136 }
137
138 if (fDisplay && !fDisplay->GetCanvas("Pedestals"))
139 fDisplay->Read();
140
141 return kTRUE;
142}
143
144
145void MJPedestal::DisplayResult(MParList &plist)
146{
147 if (!fDisplay)
148 return;
149
150 //
151 // Update display
152 //
153 TString title = fDisplay->GetTitle();
154 title += "-- Pedestal ";
155 if (fRuns) // FIXME: What to do if an environmentfile was used?
156 title += fRuns->GetRunsAsString();
157 title += " --";
158 fDisplay->SetTitle(title);
159
160 //
161 // Get container from list
162 //
163 MGeomCam &geomcam = *(MGeomCam*)plist.FindObject("MGeomCam");
164
165 //
166 // Create container to display
167 //
168 MHCamera disp0(geomcam, "MPedestalCam;ped", "Mean Pedestal");
169 MHCamera disp1(geomcam, "MPedestalCam;RMS", "Pedestal RMS");
170
171 disp0.SetCamContent(fPedestalCam, 0);
172 disp0.SetCamError (fPedestalCam, 1);
173
174 disp1.SetCamContent(fPedestalCam, 2);
175 disp1.SetCamError (fPedestalCam, 3);
176
177 disp0.SetYTitle("P [fadc/slice]");
178 disp1.SetYTitle("\\sigma_{P} [fadc/slice]");
179
180
181 //
182 // Display data
183 //
184 TCanvas &c3 = fDisplay->AddTab("Pedestals");
185 c3.Divide(2,3);
186
187 if (!fDataCheck)
188 {
189 disp0.CamDraw(c3, 1, 2, 1);
190 disp1.CamDraw(c3, 2, 2, 6);
191 return;
192 }
193
194 c3.cd(1);
195 gPad->SetBorderMode(0);
196 gPad->SetTicks();
197 MHCamera *obj1=(MHCamera*)disp0.DrawCopy("hist");
198 obj1->SetDirectory(NULL);
199
200 //
201 // for the datacheck, fix the ranges!!
202 //
203 obj1->SetMinimum(fgPedestalMin);
204 obj1->SetMaximum(fgPedestalMax);
205
206 //
207 // set reference lines
208 //
209 DisplayReferenceLines(obj1,0);
210
211 //
212 // end reference lines
213 //
214 c3.cd(3);
215 gPad->SetBorderMode(0);
216 obj1->SetPrettyPalette();
217 obj1->Draw();
218
219 c3.cd(5);
220 gPad->SetBorderMode(0);
221 gPad->SetTicks();
222 obj1->DrawProjection(7);
223
224 c3.cd(2);
225 gPad->SetBorderMode(0);
226 gPad->SetTicks();
227 MHCamera *obj2=(MHCamera*)disp1.DrawCopy("hist");
228 obj2->SetDirectory(NULL);
229 obj2->SetMinimum(fgPedRmsMin);
230 obj2->SetMaximum(fgPedRmsMax);
231
232 //
233 // set reference lines
234 //
235 DisplayReferenceLines(obj1,1);
236
237 c3.cd(4);
238 gPad->SetBorderMode(0);
239 obj2->SetPrettyPalette();
240 obj2->Draw();
241
242 c3.cd(6);
243 gPad->SetBorderMode(0);
244 gPad->SetTicks();
245
246 TArrayI inner(1);
247 inner[0] = 0;
248
249 TArrayI outer(1);
250 outer[0] = 1;
251
252 if (geomcam.InheritsFrom("MGeomCamMagic"))
253 {
254 TArrayI s0(6);
255 s0[0] = 6;
256 s0[1] = 1;
257 s0[2] = 2;
258 s0[3] = 3;
259 s0[4] = 4;
260 s0[5] = 5;
261
262 TArrayI s1(3);
263 s1[0] = 6;
264 s1[1] = 1;
265 s1[2] = 2;
266
267 TArrayI s2(3);
268 s2[0] = 3;
269 s2[1] = 4;
270 s2[2] = 5;
271
272 gPad->Clear();
273 TVirtualPad *pad = gPad;
274 pad->Divide(2,1);
275
276 TH1D *inout[2];
277 inout[0] = disp1.ProjectionS(s0, inner, "Inner");
278 inout[1] = disp1.ProjectionS(s0, outer, "Outer");
279
280 inout[0]->SetDirectory(NULL);
281 inout[1]->SetDirectory(NULL);
282
283 for (int i=0; i<2; i++)
284 {
285 TLegend *leg2 = new TLegend(0.6,0.2,0.9,0.55);
286 leg2->SetHeader(inout[i]->GetName());
287 pad->cd(i+1);
288 inout[i]->SetLineColor(kRed+i);
289 inout[i]->SetBit(kCanDelete);
290 inout[i]->Draw();
291 inout[i]->Fit("gaus","Q");
292 leg2->AddEntry(inout[i],inout[i]->GetName(),"l");
293
294 //
295 // Display the outliers as dead and noisy pixels
296 //
297 DisplayOutliers(inout[i]);
298
299 //
300 // Display the two half of the camera separately
301 //
302 TH1D *half[2];
303 half[0] = disp1.ProjectionS(s1, i==0 ? inner : outer , "Sector 6-1-2");
304 half[1] = disp1.ProjectionS(s2, i==0 ? inner : outer , "Sector 3-4-5");
305
306 for (int j=0; j<2; j++)
307 {
308 half[j]->SetLineColor(kRed+i+2*j+1);
309 half[j]->SetDirectory(NULL);
310 half[j]->SetBit(kCanDelete);
311 half[j]->Draw("same");
312 leg2->AddEntry(half[j],half[j]->GetName(),"l");
313 }
314 leg2->Draw();
315 }
316 }
317}
318
319void MJPedestal::DisplayReferenceLines(MHCamera *cam, const Int_t what) const
320{
321 const Double_t x = cam->GetNbinsX();
322
323 TLine line;
324 line.SetLineStyle(kDashed);
325 line.SetLineWidth(3);
326
327 line.SetLineColor(kBlue);
328 TLine *l1 = line.DrawLine(0, what ? fgRefPedRmsGalactic : fgRefPedGalactic,
329 x, what ? fgRefPedRmsGalactic : fgRefPedGalactic);
330
331 line.SetLineColor(kYellow);
332 TLine *l2 = line.DrawLine(0, what ? fgRefPedRmsExtraGalactic : fgRefPedExtraGalactic,
333 x, what ? fgRefPedRmsExtraGalactic : fgRefPedExtraGalactic);
334
335 line.SetLineColor(kMagenta);
336 TLine *l3 = line.DrawLine(0, what ? fgRefPedRmsClosedLids : fgRefPedClosedLids,
337 x, what ? fgRefPedRmsClosedLids : fgRefPedClosedLids);
338
339
340 TLegend *leg = new TLegend(0.4,0.75,0.7,0.99);
341 leg->SetBit(kCanDelete);
342 leg->AddEntry(l1, "Galactic Source","l");
343 leg->AddEntry(l2, "Extra-Galactic Source","l");
344 leg->AddEntry(l3, "Closed Lids","l");
345 leg->Draw();
346}
347
348void MJPedestal::DisplayOutliers(TH1D *hist) const
349{
350 const Float_t mean = hist->GetFunction("gaus")->GetParameter(1);
351 const Float_t lolim = mean - 3.5*hist->GetFunction("gaus")->GetParameter(2);
352 const Float_t uplim = mean + 3.5*hist->GetFunction("gaus")->GetParameter(2);
353 const Stat_t dead = hist->Integral(0,hist->FindBin(lolim)-1);
354 const Stat_t noisy = hist->Integral(hist->FindBin(uplim)+1,hist->GetNbinsX()+1);
355
356 TLatex deadtex;
357 deadtex.SetTextSize(0.06);
358 deadtex.DrawLatex(0.1,hist->GetBinContent(hist->GetMaximumBin())/1.1,Form("%3i dead pixels",(Int_t)dead));
359
360 TLatex noisytex;
361 noisytex.SetTextSize(0.06);
362 noisytex.DrawLatex(0.1,hist->GetBinContent(hist->GetMaximumBin())/1.2,Form("%3i noisy pixels",(Int_t)noisy));
363}
364
365Bool_t MJPedestal::WriteResult()
366{
367 if (fOutputPath.IsNull())
368 return kTRUE;
369
370 const TString oname(GetOutputFile());
371
372 *fLog << inf << "Writing to file: " << oname << endl;
373
374 TFile file(oname, "UPDATE");
375
376 if (fDisplay && fDisplay->Write()<=0)
377 {
378 *fLog << err << "Unable to write MStatusDisplay to " << oname << endl;
379 return kFALSE;
380 }
381
382 if (fPedestalCam.Write()<=0)
383 {
384 *fLog << err << "Unable to write MPedestalCam to " << oname << endl;
385 return kFALSE;
386 }
387
388 if (fBadPixels.Write()<=0)
389 {
390 *fLog << err << "Unable to write MBadPixelsCam to " << oname << endl;
391 return kFALSE;
392 }
393
394 return kTRUE;
395}
396
397void MJPedestal::SetOutputPath(const char *path)
398{
399 fOutputPath = path;
400 if (fOutputPath.EndsWith("/"))
401 fOutputPath = fOutputPath(0, fOutputPath.Length()-1);
402}
403
404void MJPedestal::SetEnv(const char *env)
405{
406 if (fEnv)
407 delete fEnv;
408 fEnv = new TEnv(env);
409}
410
411Bool_t MJPedestal::Process()
412{
413 if (!ReadPedestalCam())
414 return ProcessFile();
415
416 return kTRUE;
417}
418
419Bool_t MJPedestal::ProcessFile()
420{
421 if (!fRuns && !fEnv)
422 {
423 *fLog << err << "Neither AddRuns was called nor SetEnv was used... abort." << endl;
424 return kFALSE;
425 }
426 if (fRuns && fRuns->GetNumRuns() != fRuns->GetNumEntries())
427 {
428 *fLog << err << "Number of files found doesn't match number of runs... abort." << endl;
429 return kFALSE;
430 }
431
432 *fLog << inf;
433 fLog->Separator(GetDescriptor());
434 *fLog << "Calculate MPedestalCam from Runs ";
435 if (fRuns)
436 *fLog << fRuns->GetRunsAsString() << endl;
437 else
438 *fLog << "in " << fEnv->GetName() << endl;
439 *fLog << endl;
440
441 MParList plist;
442 MTaskList tlist;
443 plist.AddToList(&tlist);
444
445 MReadMarsFile read("Events");
446 MRawFileRead rawread(NULL);
447
448 if (fDataCheck)
449 {
450 if (fRuns)
451 rawread.AddFiles(*fRuns);
452 tlist.AddToList(&rawread);
453 }
454 else
455 {
456 read.DisableAutoScheme();
457 if (fRuns)
458 static_cast<MRead&>(read).AddFiles(*fRuns);
459 tlist.AddToList(&read);
460 }
461 // Enable logging to file
462 //*fLog.SetOutputFile(lname, kTRUE);
463
464 // Setup Tasklist
465 plist.AddToList(&fPedestalCam);
466 plist.AddToList(&fBadPixels);
467
468 MGeomApply geomapl;
469 MBadPixelsMerge merge(&fBadPixels);
470 MPedCalcPedRun pedcalc;
471
472 if (fExtractor)
473 {
474 pedcalc.SetWindowSize((Int_t)fExtractor->GetNumHiGainSamples());
475 pedcalc.SetRange(fExtractor->GetHiGainFirst(), fExtractor->GetHiGainLast());
476 }
477 else
478 {
479 *fLog << warn << GetDescriptor();
480 *fLog << ": No extractor has been chosen, take default number of FADC samples " << endl;
481 }
482
483 tlist.AddToList(&geomapl);
484 tlist.AddToList(&merge);
485 tlist.AddToList(&pedcalc);
486
487 //
488 // Execute the eventloop
489 //
490 MEvtLoop evtloop(fName);
491 evtloop.SetParList(&plist);
492 evtloop.SetDisplay(fDisplay);
493 evtloop.SetLogStream(fLog);
494 if (fEnv)
495 evtloop.ReadEnv(*fEnv);
496
497 // Execute first analysis
498 if (!evtloop.Eventloop())
499 {
500 *fLog << err << GetDescriptor() << ": Failed." << endl;
501 return kFALSE;
502 }
503
504 tlist.PrintStatistics();
505
506 DisplayResult(plist);
507
508 if (!WriteResult())
509 return kFALSE;
510
511 *fLog << inf << GetDescriptor() << ": Done." << endl;
512
513 return kTRUE;
514}
Note: See TracBrowser for help on using the repository browser.