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

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