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

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