source: trunk/MagicSoft/Mars/mjobs/MJCalibration.cc@ 2992

Last change on this file since 2992 was 2992, checked in by tbretz, 21 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!
20! Copyright: MAGIC Software Development, 2000-2004
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MJCalibration
28//
29/////////////////////////////////////////////////////////////////////////////
30#include "MJCalibration.h"
31
32#include <TF1.h>
33#include <TFile.h>
34#include <TStyle.h>
35#include <TCanvas.h>
36#include <TSystem.h>
37
38#include "MLog.h"
39#include "MLogManip.h"
40
41#include "MRunIter.h"
42#include "MParList.h"
43#include "MTaskList.h"
44#include "MEvtLoop.h"
45
46#include "MHCamera.h"
47
48#include "MPedestalCam.h"
49#include "MCalibrationCam.h"
50
51#include "MReadMarsFile.h"
52#include "MGeomApply.h"
53#include "MExtractSignal.h"
54#include "MCalibrationCalc.h"
55
56#include "MJCalibration.h"
57#include "MStatusDisplay.h"
58
59ClassImp(MJCalibration);
60
61using namespace std;
62
63MJCalibration::MJCalibration(const char *name, const char *title) : fRuns(0)
64{
65 fName = name ? name : "MJCalibration";
66 fTitle = title ? title : "Tool to create a pedestal file (MPedestalCam)";
67}
68
69void MJCalibration::DrawProjection(MHCamera *obj1, Int_t fit) const
70{
71 TH1D *obj2 = (TH1D*)obj1->Projection();
72 obj2->Draw();
73 obj2->SetBit(kCanDelete);
74
75 const Double_t min = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst());
76 const Double_t max = obj2->GetBinCenter(obj2->GetXaxis()->GetLast());
77 const Double_t integ = obj2->Integral("width")/2.5;
78 const Double_t mean = obj2->GetMean();
79 const Double_t rms = obj2->GetRMS();
80 const Double_t width = max-min;
81
82 TF1 *f=0;
83 switch (fit)
84 {
85 case 0:
86 f = new TF1("sgaus", "gaus(0)", min, max);
87 f->SetLineColor(kYellow);
88 f->SetBit(kCanDelete);
89 f->SetParNames("Area", "#mu", "#sigma");
90 f->SetParameters(integ/rms, mean, rms);
91 f->SetParLimits(0, 0, integ);
92 f->SetParLimits(1, min, max);
93 f->SetParLimits(2, 0, width/1.5);
94
95 obj2->Fit(f, "QLR");
96 break;
97
98 case 1:
99 f = new TF1("dgaus", "gaus(0)+gaus(3)", min, max);
100 f->SetLineColor(kYellow);
101 f->SetBit(kCanDelete);
102 f->SetParNames("A1", "#mu1", "#sigma1", "A2", "#mu2", "#sigma2");
103 f->SetParameters(integ/width, max-width/6, width/4, integ/width, min+width/6, width/4);
104 f->SetParLimits(0, 0, integ);
105 f->SetParLimits(1, min, max);
106 f->SetParLimits(2, 0, width/3);
107 f->SetParLimits(3, 0, integ);
108 f->SetParLimits(4, min, max);
109 f->SetParLimits(5, 0, width/3);
110
111 obj2->Fit(f, "QLR");
112 break;
113
114 case 2:
115 f = new TF1("tgaus", "gaus(0)+gaus(3)+gaus(6)", min, max);
116 f->SetLineColor(kYellow);
117 f->SetBit(kCanDelete);
118 f->SetParNames("A1", "#mu1", "#sigma1", "A2", "#mu2", "#sigma2", "A3", "#mu3", "#sigma3");
119 f->SetParameters(integ/width, max-width/6, width/4, integ/width, min+width/6, width/4, integ/width,min+width/6, width/2);
120 f->SetParLimits(0, 0, integ);
121 f->SetParLimits(1, min, max);
122 f->SetParLimits(2, 0, width/4);
123 f->SetParLimits(3, 0, integ);
124 f->SetParLimits(4, min, max);
125 f->SetParLimits(5, 0, width/4);
126 f->SetParLimits(6, 0, integ);
127 f->SetParLimits(7, min, max);
128 f->SetParLimits(8, 0, width/2);
129
130 obj2->Fit(f, "QLR");
131 break;
132
133 case 3:
134 obj2->Fit("pol0", "Q");
135 obj2->GetFunction("pol0")->SetLineColor(kYellow);
136 break;
137
138 case 9:
139 break;
140
141 default:
142 obj2->Fit("gaus", "Q");
143 obj2->GetFunction("gaus")->SetLineColor(kYellow);
144 break;
145 }
146}
147
148void MJCalibration::CamDraw(TCanvas &c, Int_t x, Int_t y, MHCamera &cam1, Int_t fit)
149{
150 c.cd(x);
151 gPad->SetBorderMode(0);
152 MHCamera *obj1=(MHCamera*)cam1.DrawCopy("hist");
153 obj1->AddNotify(fCalibrationCam);
154
155 c.cd(x+y);
156 gPad->SetBorderMode(0);
157 obj1->Draw();
158
159 c.cd(x+2*y);
160 gPad->SetBorderMode(0);
161 DrawProjection(obj1, fit);
162
163 // tbretz: Using gStyle in this context is totally useless!
164 //const Float_t he = gStyle->GetStatH();
165 //const Float_t wi = gStyle->GetStatH();
166 //gStyle->SetStatH(0.4);
167 //gStyle->SetStatW(0.25);
168 //gStyle->SetStatH(he);
169 //gStyle->SetStatW(wi);
170}
171
172void MJCalibration::DisplayResult(MParList &plist)
173{
174 if (!fDisplay)
175 return;
176
177 //
178 // Update display
179 //
180 TString title = fDisplay->GetTitle();
181 title += "-- Calibration ";
182 title += fRuns->GetRunsAsString();
183 title += " --";
184 fDisplay->SetTitle(title);
185
186 //
187 // Get container from list
188 //
189 MGeomCam &geomcam = *(MGeomCam*)plist.FindObject("MGeomCam");
190
191 // Create histograms to display
192 MHCamera disp1 (geomcam, "Cal;Charge", "Fitted Mean Charges");
193 MHCamera disp3 (geomcam, "Cal;SigmaCharge", "Sigma of Fitted Charges");
194 MHCamera disp5 (geomcam, "Cal;FitProb", "Probability of Fit");
195 /*
196 MHCamera disp6 (geomcam, "Cal;Time", "Arrival Times");
197 MHCamera disp7 (geomcam, "Cal;SigmaTime", "Sigma of Arrival Times");
198 MHCamera disp8 (geomcam, "Cal;TimeChiSquare", "Chi Square of Time Fit");
199 */
200// MHCamera disp9 (geomcam, "Cal;Ped", "Pedestals");
201// MHCamera disp10(geomcam, "Cal;PedRms", "Pedestal RMS");
202// MHCamera disp11(geomcam, "Cal;RSigma", "Reduced Sigmas");
203 MHCamera disp12(geomcam, "Cal;FFactorPhe", "Nr. of Phe's (F-Factor Method)");
204 MHCamera disp13(geomcam, "Cal;FFactorConv", "Conversion Factor (F-Factor Method)");
205 MHCamera disp14(geomcam, "Cal;BlindPixPhe", "Nr. of Photons inside plexiglass (Blind Pixel Method)");
206 MHCamera disp15(geomcam, "Cal;BlindPixConv", "Conversion Factor (Blind Pixel Method)");
207// MHCamera disp16(geomcam, "Cal;RSigma/Charge", "Reduced Sigma per Charge");
208
209 disp1.SetCamContent(fCalibrationCam, 0);
210 disp1.SetCamError(fCalibrationCam, 1);
211
212 disp3.SetCamContent(fCalibrationCam, 2);
213 disp3.SetCamError(fCalibrationCam, 3);
214
215 disp5.SetCamContent(fCalibrationCam, 4);
216
217 /*
218 disp6.SetCamContent(fCalibrationCam, 5);
219 disp6.SetCamError(fCalibrationCam, 6);
220 disp7.SetCamContent(fCalibrationCam, 6);
221 disp8.SetCamContent(fCalibrationCam, 7);
222 */
223/*
224 disp9.SetCamContent(calcam, 8);
225 disp9.SetCamError(calcam, 9);
226 disp10.SetCamContent(calcam, 9);
227 disp11.SetCamContent(fCalibrationCam, 10);
228 */
229 disp12.SetCamContent(fCalibrationCam, 11);
230 disp12.SetCamError(fCalibrationCam, 12);
231 disp13.SetCamContent(fCalibrationCam, 13);
232 disp13.SetCamError(fCalibrationCam, 14);
233
234 disp14.SetCamContent(fCalibrationCam, 15);
235 disp15.SetCamContent(fCalibrationCam, 16);
236 // disp16.SetCamContent(fCalibrationCam, 17);
237
238 disp1.SetYTitle("Charge [FADC units]");
239 disp3.SetYTitle("\\sigma_{Charge} [FADC units]");
240 disp5.SetYTitle("P_{Charge} [1]");
241 /*
242 disp6.SetYTitle("Arr. Time [FADC slice Nr.]");
243 disp7.SetYTitle("\\sigma_{Time} [FADC slices]");
244 disp8.SetYTitle("\\chi^{2}_{Time} [1]");
245 disp9.SetYTitle("Ped [FADC Counts ]");
246 disp10.SetYTitle("RMS_{Ped} [FADC Counts ]");
247 disp11.SetYTitle("\\sqrt{\\sigma^{2}_{Charge} - RMS^{2}_{Ped}} [FADC units]");
248 */
249 disp12.SetYTitle("Nr. Photo-Electrons [1]");
250 disp13.SetYTitle("Conversion Factor [PhE/FADC Count]");
251 disp14.SetYTitle("Nr. Photons [1]");
252 disp15.SetYTitle("Conversion Factor [Phot/FADC Count]");
253// disp16.SetYTitle("Reduced Sigma / Charge [1]");
254
255 gStyle->SetOptStat(1111);
256 gStyle->SetOptFit();
257
258 TCanvas &c1 = fDisplay->AddTab("FitCharge");
259 c1.Divide(2, 3);
260
261 CamDraw(c1, 1, 2, disp1, 1);
262 CamDraw(c1, 2, 2, disp3, 2);
263
264 // Fit Probability
265 TCanvas &c2 = fDisplay->AddTab("FitProb");
266 c2.Divide(1,3);
267
268 CamDraw(c2, 1, 1, disp5, 3);
269/*
270 // Times
271 TCanvas &c3 = fDisplay->AddTab("FitTimes");
272 c3.Divide(3,3);
273
274 CamDraw(c3, 1, 3, disp6, 1);
275 CamDraw(c3, 2, 3, disp7, 0);
276 CamDraw(c3, 3, 3, disp8, 0);
277
278 // Pedestals
279 TCanvas &c4 = d3->AddTab("Pedestals");
280 c4.Divide(2,3);
281
282 CamDraw(c4, 1, 2, disp9, 0);
283 CamDraw(c4, 2, 2, disp10, 1);
284
285 // Reduced Sigmas
286 TCanvas &c5 = fDisplay->AddTab("RedSigma");
287 c5.Divide(2,3);
288
289 CamDraw(c5, 1, 2, disp11, 2);
290 CamDraw(c5, 2, 2, disp16, 2);
291 */
292
293 // F-Factor Method
294 TCanvas &c6 = fDisplay->AddTab("F-Factor");
295 c6.Divide(2,3);
296
297 CamDraw(c6, 1, 2, disp12, 1);
298 CamDraw(c6, 2, 2, disp13, 1);
299
300 // Blind Pixel Method
301 TCanvas &c7 = fDisplay->AddTab("BlindPix");
302 c7.Divide(2, 3);
303
304 CamDraw(c7, 1, 2, disp14, 9);
305 CamDraw(c7, 2, 2, disp15, 1);
306
307}
308
309Bool_t MJCalibration::WriteResult()
310{
311 if (fOutputPath.IsNull())
312 return kTRUE;
313
314 const TString oname(GetOutputFile());
315
316 *fLog << inf << "Writing to file: " << oname << endl;
317
318 TFile file(oname, "UPDATE");
319
320 if (fDisplay && fDisplay->Write()<=0)
321 {
322 *fLog << err << "Unable to write MStatusDisplay to " << oname << endl;
323 return kFALSE;
324 }
325
326 if (fCalibrationCam.Write()<=0)
327 {
328 *fLog << err << "Unable to write MCalibrationCam to " << oname << endl;
329 return kFALSE;
330 }
331 return kTRUE;
332
333}
334
335void MJCalibration::SetOutputPath(const char *path)
336{
337 fOutputPath = path;
338 if (fOutputPath.EndsWith("/"))
339 fOutputPath = fOutputPath(0, fOutputPath.Length()-1);
340}
341
342Bool_t MJCalibration::Process(MPedestalCam &pedcam)
343{
344 if (!ReadCalibrationCam())
345 return ProcessFile(pedcam);
346
347 return kTRUE;
348}
349
350TString MJCalibration::GetOutputFile() const
351{
352 if (!fRuns)
353 return "";
354
355 return Form("%s/%s-F1.root", (const char*)fOutputPath, (const char*)fRuns->GetRunsAsFileName());
356}
357
358Bool_t MJCalibration::ReadCalibrationCam()
359{
360 const TString fname = GetOutputFile();
361
362 if (gSystem->AccessPathName(fname, kFileExists))
363 {
364 *fLog << err << "Input file " << fname << " doesn't exist." << endl;
365 return kFALSE;
366 }
367
368 *fLog << inf << "Reading from file: " << fname << endl;
369
370 TFile file(fname, "READ");
371 if (fCalibrationCam.Read()<=0)
372 {
373 *fLog << "Unable to read MCalibrationCam from " << fname << endl;
374 return kFALSE;
375 }
376
377 if (fDisplay /*&& !fDisplay->GetCanvas("Pedestals")*/) // FIXME!
378 fDisplay->Read();
379
380 return kTRUE;
381}
382
383
384Bool_t MJCalibration::ProcessFile(MPedestalCam &pedcam)
385{
386 if (!fRuns)
387 {
388 *fLog << err << "No Runs choosen... abort." << endl;
389 return kFALSE;
390 }
391 if (fRuns->GetNumRuns() != fRuns->GetNumEntries())
392 {
393 *fLog << err << "Number of files found doesn't metch number of runs... abort." << endl;
394 return kFALSE;
395 }
396
397 *fLog << inf;
398 fLog->Separator(GetDescriptor());
399 *fLog << "Calculate MCalibrationCam from Runs " << fRuns->GetRunsAsString() << endl;
400 *fLog << endl;
401
402 MReadMarsFile read("Events");
403 read.DisableAutoScheme();
404 static_cast<MRead&>(read).AddFiles(*fRuns);
405
406 // Setup Tasklist
407 MParList plist;
408 plist.AddToList(&pedcam);
409 plist.AddToList(&fCalibrationCam);
410
411 MTaskList tlist;
412 plist.AddToList(&tlist);
413
414 MGeomApply apply;
415 MExtractSignal extract;
416 MCalibrationCalc calcalc;
417
418 //
419 // As long, as we don't have digital modules,
420 // we have to set the color of the pulser LED by hand
421 //
422 calcalc.SetPulserColor(MCalibrationCalc::kECT1);
423 //calcalc.SkipBlindPixelFit();
424
425 tlist.AddToList(&apply);
426 tlist.AddToList(&read);
427 tlist.AddToList(&extract);
428 tlist.AddToList(&calcalc);
429
430 // Create and setup the eventloop
431 MEvtLoop evtloop(fName);
432 evtloop.SetParList(&plist);
433 evtloop.SetDisplay(fDisplay);
434 evtloop.SetLogStream(fLog);
435
436 // Execute first analysis
437 if (!evtloop.Eventloop())
438 {
439 *fLog << err << GetDescriptor() << ": Failed." << endl;
440 return kFALSE;
441 }
442
443 tlist.PrintStatistics();
444
445 DisplayResult(plist);
446
447 if (!WriteResult())
448 return kFALSE;
449
450 *fLog << inf << GetDescriptor() << ": Done." << endl;
451
452 return kTRUE;
453}
Note: See TracBrowser for help on using the repository browser.