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

Last change on this file since 3059 was 3055, checked in by gaug, 22 years ago
*** empty log message ***
File size: 14.9 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 "MExtractSignal2.h"
55#include "MCalibrationCalc.h"
56
57#include "MJCalibration.h"
58#include "MStatusDisplay.h"
59
60ClassImp(MJCalibration);
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
72 TH1D *obj2 = (TH1D*)obj1->Projection();
73 obj2->Draw();
74 obj2->SetBit(kCanDelete);
75
76 const Double_t min = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst());
77 const Double_t max = obj2->GetBinCenter(obj2->GetXaxis()->GetLast());
78 const Double_t integ = obj2->Integral("width")/2.5;
79 const Double_t mean = obj2->GetMean();
80 const Double_t rms = obj2->GetRMS();
81 const Double_t width = max-min;
82
83 const TString dgausformula = "([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])"
84 "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
85
86 const TString tgausformula = "([0]-[3]-[6])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])"
87 "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])"
88 "+[6]/[8]*exp(-0.5*(x-[7])*(x-[7])/[8]/[8])";
89 TF1 *f=0;
90 switch (fit)
91 {
92 case 1:
93 f = new TF1("sgaus", "gaus(0)", min, max);
94 f->SetLineColor(kYellow);
95 f->SetBit(kCanDelete);
96 f->SetParNames("Area", "#mu", "#sigma");
97 f->SetParameters(integ/rms, mean, rms);
98 f->SetParLimits(0, 0, integ);
99 f->SetParLimits(1, min, max);
100 f->SetParLimits(2, 0, width/1.5);
101
102 obj2->Fit(f, "QLR");
103 break;
104
105 case 2:
106 f = new TF1("dgaus",dgausformula.Data(),min,max);
107 f->SetLineColor(kYellow);
108 f->SetBit(kCanDelete);
109 f->SetParNames("A_{tot}", "#mu1", "#sigma1", "A2", "#mu2", "#sigma2");
110 f->SetParameters(integ,(min+mean)/2.,width/4.,
111 integ/width/2.,(max+mean)/2.,width/4.);
112 // The left-sided Gauss
113 f->SetParLimits(0,integ-1.5 , integ+1.5);
114 f->SetParLimits(1,min+(width/10.), mean);
115 f->SetParLimits(2,0 , width/2.);
116 // The right-sided Gauss
117 f->SetParLimits(3,0 , integ);
118 f->SetParLimits(4,mean, max-(width/10.));
119 f->SetParLimits(5,0 , width/2.);
120 obj2->Fit(f,"QLRM");
121 break;
122
123 case 3:
124 f = new TF1("tgaus",tgausformula.Data(),min,max);
125 f->SetLineColor(kYellow);
126 f->SetBit(kCanDelete);
127 f->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}",
128 "A_{2}","#mu_{2}","#sigma_{2}",
129 "A_{3}","#mu_{3}","#sigma_{3}");
130 f->SetParameters(integ,(min+mean)/2,width/4.,
131 integ/width/3.,(max+mean)/2.,width/4.,
132 integ/width/3.,mean,width/2.);
133 // The left-sided Gauss
134 f->SetParLimits(0,integ-1.5,integ+1.5);
135 f->SetParLimits(1,min+(width/10.),mean);
136 f->SetParLimits(2,width/15.,width/2.);
137 // The right-sided Gauss
138 f->SetParLimits(3,0.,integ);
139 f->SetParLimits(4,mean,max-(width/10.));
140 f->SetParLimits(5,width/15.,width/2.);
141 // The Gauss describing the outliers
142 f->SetParLimits(6,0.,integ);
143 f->SetParLimits(7,min,max);
144 f->SetParLimits(8,width/4.,width/1.5);
145 obj2->Fit(f,"QLRM");
146 break;
147
148 case 4:
149 obj2->Fit("pol0", "Q");
150 obj2->GetFunction("pol0")->SetLineColor(kYellow);
151 break;
152
153 case 9:
154 break;
155
156 default:
157 obj2->Fit("gaus", "Q");
158 obj2->GetFunction("gaus")->SetLineColor(kYellow);
159 break;
160 }
161}
162
163void MJCalibration::CamDraw(TCanvas &c, const Int_t x, const Int_t y, const MHCamera &cam1, const Int_t fit)
164{
165
166
167 c.cd(x);
168 gPad->SetBorderMode(0);
169 MHCamera *obj1=(MHCamera*)cam1.DrawCopy("hist");
170 obj1->AddNotify(fCalibrationCam);
171
172 c.cd(x+y);
173 gPad->SetBorderMode(0);
174 obj1->Draw();
175
176 if (fit)
177 {
178 c.cd(x+2*y);
179 gPad->SetBorderMode(0);
180 DrawProjection(obj1, fit);
181 }
182}
183
184
185void MJCalibration::DisplayResult(MParList &plist)
186{
187 if (!fDisplay)
188 return;
189
190 //
191 // Update display
192 //
193 TString title = fDisplay->GetTitle();
194 title += "-- Calibration ";
195 title += fRuns->GetRunsAsString();
196 title += " --";
197 fDisplay->SetTitle(title);
198
199 //
200 // Get container from list
201 //
202 MGeomCam &geomcam = *(MGeomCam*)plist.FindObject("MGeomCam");
203
204 // Create histograms to display
205 MHCamera disp1 (geomcam, "Cal;Charge", "Fitted Mean Charges");
206 MHCamera disp2 (geomcam, "Cal;SigmaCharge", "Sigma of Fitted Charges");
207 MHCamera disp3 (geomcam, "Cal;FitProb", "Probability of Fit");
208 MHCamera disp4 (geomcam, "Cal;RSigma", "Reduced Sigmas");
209 MHCamera disp5 (geomcam, "Cal;RSigma/Charge", "Reduced Sigma per Charge");
210 MHCamera disp6 (geomcam, "Cal;FFactorPhe", "Nr. of Phe's (F-Factor Method)");
211 MHCamera disp7 (geomcam, "Cal;FFactorConv", "Conversion Factor (F-Factor Method)");
212 MHCamera disp8 (geomcam, "Cal;FFactorFFactor", "Total F-Factor (F-Factor Method)");
213 MHCamera disp9 (geomcam, "Cal;BlindPixPh", "Nr. of Photons inside plexiglass (Blind Pixel Method)");
214 MHCamera disp10 (geomcam, "Cal;BlindPixConv", "Conversion Factor (Blind Pixel Method)");
215 MHCamera disp11 (geomcam, "Cal;BlindPixFFactor","Total F-Factor (Blind Pixel Method)");
216 MHCamera disp12 (geomcam, "Cal;PINDiodePh", "Nr. of Photons outside plexiglass (PIN Diode Method)");
217 MHCamera disp13 (geomcam, "Cal;PINDiodeConv", "Conversion Factor (PIN Diode Method)");
218 MHCamera disp14 (geomcam, "Cal;PINDiodeFFactor","Total F-Factor (PIN Diode Method)");
219 MHCamera disp15 (geomcam, "Cal;Excluded", "Pixels previously excluded");
220 MHCamera disp16 (geomcam, "Cal;NotFitted", "Pixels that could not be fitted");
221 MHCamera disp17 (geomcam, "Cal;NotFitValid", "Pixels with not valid fit results");
222 MHCamera disp18 (geomcam, "Cal;Oscillating", "Oscillating Pixels");
223 MHCamera disp19 (geomcam, "Cal;Saturation", "Pixels with saturated Hi Gain");
224
225
226 // Fitted charge means and sigmas
227 disp1.SetCamContent(fCalibrationCam, 0);
228 disp1.SetCamError( fCalibrationCam, 1);
229 disp2.SetCamContent(fCalibrationCam, 2);
230 disp2.SetCamError( fCalibrationCam, 3);
231 // Fit probabilities
232 disp3.SetCamContent(fCalibrationCam, 4);
233
234 // Reduced Sigmas and reduced sigmas per charge
235 disp4.SetCamContent(fCalibrationCam, 5);
236 disp4.SetCamError( fCalibrationCam, 6);
237 disp5.SetCamContent(fCalibrationCam, 7);
238 disp5.SetCamError( fCalibrationCam, 8);
239
240 // F-Factor Method
241 disp6.SetCamContent(fCalibrationCam, 9);
242 disp6.SetCamError( fCalibrationCam, 10);
243 disp7.SetCamContent(fCalibrationCam, 11);
244 disp7.SetCamError( fCalibrationCam, 12);
245 disp8.SetCamContent(fCalibrationCam, 13);
246 disp8.SetCamError( fCalibrationCam, 14);
247
248 /// Blind Pixel Method
249 disp9.SetCamContent(fCalibrationCam, 15);
250 disp9.SetCamError( fCalibrationCam, 16);
251 disp10.SetCamContent(fCalibrationCam,17);
252 disp10.SetCamError( fCalibrationCam,18);
253 disp11.SetCamContent(fCalibrationCam,19);
254 disp11.SetCamError( fCalibrationCam,20);
255
256 // PIN Diode Method
257 disp12.SetCamContent(fCalibrationCam,21);
258 disp12.SetCamError( fCalibrationCam,22);
259 disp13.SetCamContent(fCalibrationCam,23);
260 disp13.SetCamError( fCalibrationCam,24);
261 disp14.SetCamContent(fCalibrationCam,25);
262 disp14.SetCamError( fCalibrationCam,26);
263
264 // Pixels with defects
265 disp15.SetCamContent(fCalibrationCam,27);
266 disp16.SetCamContent(fCalibrationCam,28);
267 disp17.SetCamContent(fCalibrationCam,29);
268 disp18.SetCamContent(fCalibrationCam,30);
269
270 // Lo Gain calibration
271 disp19.SetCamContent(fCalibrationCam,31);
272
273
274 disp1.SetYTitle("Charge [FADC units]");
275 disp2.SetYTitle("\\sigma_{Charge} [FADC units]");
276 disp3.SetYTitle("P_{Charge} [1]");
277
278 disp4.SetYTitle("\\sqrt{\\sigma^{2}_{Charge} - RMS^{2}_{Ped}} [FADC Counts]");
279 disp5.SetYTitle("Reduced Sigma / Mean Charge [1]");
280
281 disp6.SetYTitle("Nr. Photo-Electrons [1]");
282 disp7.SetYTitle("Conversion Factor [PhE/FADC Count]");
283 disp8.SetYTitle("\\sqrt{N_{PhE}}*\\sigma_{Charge}/\\mu_{Charge} [1]");
284
285 disp9.SetYTitle("Nr. Photons [1]");
286 disp10.SetYTitle("Conversion Factor [Phot/FADC Count]");
287 disp11.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1]");
288
289 disp12.SetYTitle("Nr. Photons [1]");
290 disp13.SetYTitle("Conversion Factor [Phot/FADC Count]");
291 disp14.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1]");
292
293 disp15.SetYTitle("[1]");
294 disp16.SetYTitle("[1]");
295 disp17.SetYTitle("[1]");
296 disp18.SetYTitle("[1]");
297
298 gStyle->SetOptStat(1111);
299 gStyle->SetOptFit();
300
301 // Charges
302 TCanvas &c1 = fDisplay->AddTab("Fit.Charge");
303 c1.Divide(2, 3);
304
305 CamDraw(c1, 1, 2, disp1, 2);
306 CamDraw(c1, 2, 2, disp2, 2);
307
308 // Fit Probability
309 TCanvas &c2 = fDisplay->AddTab("Fit.Prob");
310 c2.Divide(1,3);
311
312 CamDraw(c2, 1, 1, disp3, 4);
313
314 // Reduced Sigmas
315 TCanvas &c3 = fDisplay->AddTab("Red.Sigma");
316 c3.Divide(2,3);
317
318 CamDraw(c3, 1, 2, disp4, 2);
319 CamDraw(c3, 2, 2, disp5, 2);
320
321 // F-Factor Method
322 TCanvas &c4 = fDisplay->AddTab("F-Factor");
323 c4.Divide(3,3);
324
325 CamDraw(c4, 1, 3, disp6, 2);
326 CamDraw(c4, 2, 3, disp7, 2);
327 CamDraw(c4, 3, 3, disp8, 2);
328
329 // Blind Pixel Method
330 TCanvas &c5 = fDisplay->AddTab("BlindPix");
331 c5.Divide(3, 3);
332
333 CamDraw(c5, 1, 3, disp9, 9);
334 CamDraw(c5, 2, 3, disp10, 2);
335 CamDraw(c5, 3, 3, disp11, 2);
336
337 // PIN Diode Method
338 TCanvas &c6 = fDisplay->AddTab("PINDiode");
339 c6.Divide(3,3);
340
341 CamDraw(c6, 1, 3, disp12, 9);
342 CamDraw(c6, 2, 3, disp13, 2);
343 CamDraw(c6, 3, 3, disp14, 2);
344
345 // Defects
346 TCanvas &c7 = fDisplay->AddTab("Defects");
347 c7.Divide(4,2);
348
349 CamDraw(c7, 1, 4, disp15, 0);
350 CamDraw(c7, 2, 4, disp16, 0);
351 CamDraw(c7, 3, 4, disp17, 0);
352 CamDraw(c7, 4, 4, disp18, 0);
353
354 // Lo Gain Calibration
355 TCanvas &c8 = fDisplay->AddTab("LowGain");
356 c8.Divide(1,3);
357
358 CamDraw(c8, 1, 4, disp19, 0);
359
360}
361
362Bool_t MJCalibration::WriteResult()
363{
364 if (fOutputPath.IsNull())
365 return kTRUE;
366
367 const TString oname(GetOutputFile());
368
369 *fLog << inf << "Writing to file: " << oname << endl;
370
371 TFile file(oname, "UPDATE");
372
373 if (fDisplay && fDisplay->Write()<=0)
374 {
375 *fLog << err << "Unable to write MStatusDisplay to " << oname << endl;
376 return kFALSE;
377 }
378
379 if (fCalibrationCam.Write()<=0)
380 {
381 *fLog << err << "Unable to write MCalibrationCam to " << oname << endl;
382 return kFALSE;
383 }
384 return kTRUE;
385
386}
387
388void MJCalibration::SetOutputPath(const char *path)
389{
390 fOutputPath = path;
391 if (fOutputPath.EndsWith("/"))
392 fOutputPath = fOutputPath(0, fOutputPath.Length()-1);
393}
394
395Bool_t MJCalibration::Process(MPedestalCam &pedcam)
396{
397 if (!ReadCalibrationCam())
398 return ProcessFile(pedcam);
399
400 return kTRUE;
401}
402
403TString MJCalibration::GetOutputFile() const
404{
405 if (!fRuns)
406 return "";
407
408 return Form("%s/%s-F1.root", (const char*)fOutputPath, (const char*)fRuns->GetRunsAsFileName());
409}
410
411Bool_t MJCalibration::ReadCalibrationCam()
412{
413 const TString fname = GetOutputFile();
414
415 if (gSystem->AccessPathName(fname, kFileExists))
416 {
417 *fLog << err << "Input file " << fname << " doesn't exist." << endl;
418 return kFALSE;
419 }
420
421 *fLog << inf << "Reading from file: " << fname << endl;
422
423 TFile file(fname, "READ");
424 if (fCalibrationCam.Read()<=0)
425 {
426 *fLog << "Unable to read MCalibrationCam from " << fname << endl;
427 return kFALSE;
428 }
429
430 if (fDisplay /*&& !fDisplay->GetCanvas("Pedestals")*/) // FIXME!
431 fDisplay->Read();
432
433 return kTRUE;
434}
435
436
437Bool_t MJCalibration::ProcessFile(MPedestalCam &pedcam)
438{
439 if (!fRuns)
440 {
441 *fLog << err << "No Runs choosen... abort." << endl;
442 return kFALSE;
443 }
444 if (fRuns->GetNumRuns() != fRuns->GetNumEntries())
445 {
446 *fLog << err << "Number of files found doesn't metch number of runs... abort." << endl;
447 return kFALSE;
448 }
449
450 *fLog << inf;
451 fLog->Separator(GetDescriptor());
452 *fLog << "Calculate MCalibrationCam from Runs " << fRuns->GetRunsAsString() << endl;
453 *fLog << endl;
454
455 MReadMarsFile read("Events");
456 read.DisableAutoScheme();
457 static_cast<MRead&>(read).AddFiles(*fRuns);
458
459 // Setup Tasklist
460 MParList plist;
461 plist.AddToList(&pedcam);
462 plist.AddToList(&fCalibrationCam);
463
464 MTaskList tlist;
465 plist.AddToList(&tlist);
466
467 MGeomApply apply;
468 MExtractSignal2 extract;
469 MCalibrationCalc calcalc;
470
471 //
472 // As long, as we don't have digital modules,
473 // we have to set the color of the pulser LED by hand
474 //
475 calcalc.SetPulserColor(MCalibrationCalc::kECT1);
476 //calcalc.SkipBlindPixelFit();
477
478 tlist.AddToList(&apply);
479 tlist.AddToList(&read);
480 tlist.AddToList(&extract);
481 tlist.AddToList(&calcalc);
482
483 // Create and setup the eventloop
484 MEvtLoop evtloop(fName);
485 evtloop.SetParList(&plist);
486 evtloop.SetDisplay(fDisplay);
487 evtloop.SetLogStream(fLog);
488
489 // Execute first analysis
490 if (!evtloop.Eventloop())
491 {
492 *fLog << err << GetDescriptor() << ": Failed." << endl;
493 return kFALSE;
494 }
495
496 tlist.PrintStatistics();
497
498 DisplayResult(plist);
499
500 if (!WriteResult())
501 return kFALSE;
502
503 *fLog << inf << GetDescriptor() << ": Done." << endl;
504
505 return kTRUE;
506}
Note: See TracBrowser for help on using the repository browser.