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

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