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

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