source: tags/Mars-V0.8.4/mjobs/MJPedestal.cc

Last change on this file was 3718, checked in by gaug, 21 years ago
*** empty log message ***
File size: 9.7 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!
20! Copyright: MAGIC Software Development, 2000-2004
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MJPedestal
28//
29/////////////////////////////////////////////////////////////////////////////
30#include "MJPedestal.h"
31
32#include <TF1.h>
33#include <TFile.h>
34#include <TStyle.h>
35#include <TString.h>
36#include <TCanvas.h>
37#include <TSystem.h>
38
39#include "MLog.h"
40#include "MLogManip.h"
41
42#include "MRunIter.h"
43#include "MParList.h"
44#include "MTaskList.h"
45#include "MEvtLoop.h"
46
47#include "MStatusDisplay.h"
48
49#include "MGeomCam.h"
50#include "MHCamera.h"
51#include "MPedestalCam.h"
52#include "MBadPixelsCam.h"
53
54#include "MReadMarsFile.h"
55#include "MGeomApply.h"
56#include "MBadPixelsMerge.h"
57#include "MPedCalcPedRun.h"
58
59ClassImp(MJPedestal);
60
61using namespace std;
62
63MJPedestal::MJPedestal(const char *name, const char *title) : fRuns(0)
64{
65 fName = name ? name : "MJPedestal";
66 fTitle = title ? title : "Tool to create a pedestal file (MPedestalCam)";
67}
68
69TString MJPedestal::GetOutputFile() const
70{
71 if (!fRuns)
72 return "";
73
74 return Form("%s/%s-F0.root", (const char*)fOutputPath, (const char*)fRuns->GetRunsAsFileName());
75}
76
77Bool_t MJPedestal::ReadPedestalCam()
78{
79 const TString fname = GetOutputFile();
80
81 if (gSystem->AccessPathName(fname, kFileExists))
82 {
83 *fLog << err << "Input file " << fname << " doesn't exist." << endl;
84 return kFALSE;
85 }
86
87 *fLog << inf << "Reading from file: " << fname << endl;
88
89 TFile file(fname, "READ");
90 if (fPedestalCam.Read()<=0)
91 {
92 *fLog << err << "Unable to read MPedestalCam from " << fname << endl;
93 return kFALSE;
94 }
95
96 if (file.FindKey("MBadPixelsCam"))
97 {
98 MBadPixelsCam bad;
99 if (bad.Read()<=0)
100 {
101 *fLog << err << "Unable to read MBadPixelsCam from " << fname << endl;
102 return kFALSE;
103 }
104 fBadPixels.Merge(bad);
105 }
106
107 if (fDisplay && !fDisplay->GetCanvas("Pedestals"))
108 fDisplay->Read();
109
110 return kTRUE;
111}
112
113void MJPedestal::DrawProjection(MHCamera *obj1, Int_t fit) const
114{
115 TH1D *obj2 = (TH1D*)obj1->Projection(obj1->GetName());
116 obj2->SetDirectory(0);
117 obj2->Draw();
118 obj2->SetBit(kCanDelete);
119 gPad->SetLogy();
120
121 if (obj1->GetGeomCam().InheritsFrom("MGeomCamMagic"))
122 {
123 TArrayI s0(3);
124 s0[0] = 6;
125 s0[1] = 1;
126 s0[2] = 2;
127
128 TArrayI s1(3);
129 s1[0] = 3;
130 s1[1] = 4;
131 s1[2] = 5;
132
133 TArrayI inner(1);
134 inner[0] = 0;
135
136 TArrayI outer(1);
137 outer[0] = 1;
138
139 // Just to get the right (maximum) binning
140 TH1D *half[4];
141 half[0] = obj1->ProjectionS(s0, inner, "Sector 6-1-2 Inner");
142 half[1] = obj1->ProjectionS(s1, inner, "Sector 3-4-5 Inner");
143 half[2] = obj1->ProjectionS(s0, outer, "Sector 6-1-2 Outer");
144 half[3] = obj1->ProjectionS(s1, outer, "Sector 3-4-5 Outer");
145
146 for (int i=0; i<4; i++)
147 {
148 half[i]->SetLineColor(kRed+i);
149 half[i]->SetDirectory(0);
150 half[i]->SetBit(kCanDelete);
151 half[i]->Draw("same");
152 }
153 }
154
155 const Double_t min = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst());
156 const Double_t max = obj2->GetBinCenter(obj2->GetXaxis()->GetLast());
157 const Double_t integ = obj2->Integral("width")/2.5066283;
158 const Double_t mean = obj2->GetMean();
159 const Double_t rms = obj2->GetRMS();
160 const Double_t width = max-min;
161
162 if (rms==0 || width==0)
163 return;
164
165 const TString dgausformula("([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])"
166 "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])");
167
168 TF1 *f=0;
169 switch (fit)
170 {
171 // FIXME: MAYBE add function to TH1?
172 case 0:
173 f = new TF1("sgaus", "gaus(0)", min, max);
174 f->SetLineColor(kYellow);
175 f->SetBit(kCanDelete);
176 f->SetParNames("Area", "#mu", "#sigma");
177 f->SetParameters(integ/rms, mean, rms);
178 f->SetParLimits(0, 0, integ);
179 f->SetParLimits(1, min, max);
180 f->SetParLimits(2, 0, width/1.5);
181 obj2->Fit(f, "QLRM");
182 break;
183
184 case 1:
185 f = new TF1("dgaus", dgausformula, min, max);
186 f->SetLineColor(kYellow);
187 f->SetBit(kCanDelete);
188 f->SetParNames("A_{tot}", "#mu_{1}", "#sigma_{1}", "A_{2}", "#mu_{2}", "#sigma_{2}");
189 f->SetParameters(integ, (min+mean)/2, width/4,
190 integ/width/2, (max+mean)/2, width/4);
191 // The left-sided Gauss
192 f->SetParLimits(0, integ-1.5, integ+1.5);
193 f->SetParLimits(1, min+(width/10), mean);
194 f->SetParLimits(2, 0, width/2);
195 // The right-sided Gauss
196 f->SetParLimits(3, 0, integ);
197 f->SetParLimits(4, mean, max-(width/10));
198 f->SetParLimits(5, 0, width/2);
199 obj2->Fit(f, "QLRM");
200 break;
201
202 default:
203 obj2->Fit("gaus", "Q");
204 obj2->GetFunction("gaus")->SetLineColor(kYellow);
205 break;
206 }
207}
208
209void MJPedestal::CamDraw(TCanvas &c, Int_t x, Int_t y, MHCamera &cam1, Int_t fit)
210{
211 c.cd(x);
212 gPad->SetBorderMode(0);
213 MHCamera *obj1=(MHCamera*)cam1.DrawCopy("hist");
214
215 obj1->AddNotify(&fPedestalCam);
216
217 c.cd(x+y);
218 gPad->SetBorderMode(0);
219 obj1->Draw();
220
221 c.cd(x+2*y);
222 gPad->SetBorderMode(0);
223 DrawProjection(obj1, fit);
224}
225
226void MJPedestal::DisplayResult(MParList &plist)
227{
228 if (!fDisplay)
229 return;
230
231 //
232 // Update display
233 //
234 TString title = fDisplay->GetTitle();
235 title += "-- Pedestal ";
236 title += fRuns->GetRunsAsString();
237 title += " --";
238 fDisplay->SetTitle(title);
239
240 //
241 // Get container from list
242 //
243 MGeomCam &geomcam = *(MGeomCam*)plist.FindObject("MGeomCam");
244
245 //
246 // Create container to display
247 //
248 MHCamera disp0(geomcam, "MPedestalCam;ped", "Pedestals");
249 MHCamera disp1(geomcam, "MPedestalCam;var", "Sigma Pedestal");
250
251 disp0.SetCamContent(fPedestalCam, 0);
252 disp0.SetCamError (fPedestalCam, 1);
253
254 disp1.SetCamContent(fPedestalCam, 2);
255 disp1.SetCamError (fPedestalCam, 3);
256
257 disp0.SetYTitle("P [fadc/slice]");
258 disp1.SetYTitle("\\sigma_{P} [fadc/slice]");
259
260 //
261 // Display data
262 //
263 TCanvas &c3 = fDisplay->AddTab("Pedestals");
264 c3.Divide(2,3);
265
266 CamDraw(c3, 1, 2, disp0, 0);
267 CamDraw(c3, 2, 2, disp1, 1);
268}
269
270Bool_t MJPedestal::WriteResult()
271{
272 if (fOutputPath.IsNull())
273 return kTRUE;
274
275 const TString oname(GetOutputFile());
276
277 *fLog << inf << "Writing to file: " << oname << endl;
278
279 TFile file(oname, "UPDATE");
280
281 if (fDisplay && fDisplay->Write()<=0)
282 {
283 *fLog << err << "Unable to write MStatusDisplay to " << oname << endl;
284 return kFALSE;
285 }
286
287 if (fPedestalCam.Write()<=0)
288 {
289 *fLog << err << "Unable to write MPedestalCam to " << oname << endl;
290 return kFALSE;
291 }
292
293 if (fBadPixels.Write()<=0)
294 {
295 *fLog << err << "Unable to write MBadPixelsCam to " << oname << endl;
296 return kFALSE;
297 }
298
299 return kTRUE;
300}
301
302void MJPedestal::SetOutputPath(const char *path)
303{
304 fOutputPath = path;
305 if (fOutputPath.EndsWith("/"))
306 fOutputPath = fOutputPath(0, fOutputPath.Length()-1);
307}
308
309Bool_t MJPedestal::Process()
310{
311 if (!ReadPedestalCam())
312 return ProcessFile();
313
314 return kTRUE;
315}
316
317Bool_t MJPedestal::ProcessFile()
318{
319 if (!fRuns)
320 {
321 *fLog << err << "No Runs choosen... abort." << endl;
322 return kFALSE;
323 }
324 if (fRuns->GetNumRuns() != fRuns->GetNumEntries())
325 {
326 *fLog << err << "Number of files found doesn't match number of runs... abort." << endl;
327 return kFALSE;
328 }
329
330 *fLog << inf;
331 fLog->Separator(GetDescriptor());
332 *fLog << "Calculate MPedestalCam from Runs " << fRuns->GetRunsAsString() << endl;
333 *fLog << endl;
334
335 MReadMarsFile read("Events");
336 read.DisableAutoScheme();
337 static_cast<MRead&>(read).AddFiles(*fRuns);
338
339 // Enable logging to file
340 //*fLog.SetOutputFile(lname, kTRUE);
341
342 // Setup Tasklist
343 MParList plist;
344 plist.AddToList(&fPedestalCam);
345 plist.AddToList(&fBadPixels);
346
347 MTaskList tlist;
348 plist.AddToList(&tlist);
349
350 MGeomApply geomapl;
351 // MBadPixelsMerge merge(&fBadPixels);
352 //MExtractSignal sigcalc;
353 MPedCalcPedRun pedcalc;
354
355 tlist.AddToList(&read);
356 tlist.AddToList(&geomapl);
357 // tlist.AddToList(&merge);
358 //tlist.AddToList(&sigcalc);
359 tlist.AddToList(&pedcalc);
360
361 //
362 // Execute the eventloop
363 //
364 MEvtLoop evtloop(fName);
365 evtloop.SetParList(&plist);
366 evtloop.SetDisplay(fDisplay);
367 evtloop.SetLogStream(fLog);
368
369 // Execute first analysis
370 if (!evtloop.Eventloop())
371 {
372 *fLog << err << GetDescriptor() << ": Failed." << endl;
373 return kFALSE;
374 }
375
376 tlist.PrintStatistics();
377
378 DisplayResult(plist);
379
380 if (!WriteResult())
381 return kFALSE;
382
383 *fLog << inf << GetDescriptor() << ": Done." << endl;
384
385 return kTRUE;
386}
Note: See TracBrowser for help on using the repository browser.