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

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