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

Last change on this file since 3738 was 3738, checked in by gaug, 21 years ago
*** empty log message ***
File size: 37.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! Markus Gaug, 02/2004 <mailto:markus@ifae.es>
20!
21! Copyright: MAGIC Software Development, 2000-2004
22!
23!
24\* ======================================================================== */
25/////////////////////////////////////////////////////////////////////////////
26//
27// MJCalibration
28//
29// Do one calibration loop over serious of runs with the same pulser
30// colour and the same intensity. The following containers (rectangular) and
31// tasks (ellipses) are called to produce an MCalibrationChargeCam and to
32// update the MCalibrationQECam: (MCalibrate is not called from this class)
33//
34//Begin_Html
35/*
36<img src="images/CalibClasses.gif">
37*/
38//End_Html
39//
40// Different signal extractors can be chosen via the command SetExtractorLevel(UInt_t i)
41// Up to now, the following extractors are available:
42// i=1: Use MExtractSignal (fixed window)
43// i=2: Use MExtractSignal2 (sliding window: default)
44// i=3: Use MExtractSignal3 (coherent sliding window for all pixels)
45//
46// At the end of the eventloop, plots and results are displayed, depending on
47// the flags set (see DisplayResult())
48//
49// If the flag SetFullDisplay() is set, all MHCameras will be displayed.
50// if the flag SetDataCheck() is set, only the most important ones are displayed
51// Otherwise, (default: SetNormalDisplay()), a good selection of plots is given
52//
53// See also: MHCalibrationChargePix, MHCalibrationChargeCam, MHGausEvents
54// MHCalibrationChargeBlindPix, MHCalibrationChargePINDiode
55// MCalibrationChargePix, MCalibrationChargeCam, MCalibrationChargeCalc
56// MCalibrationChargeBlindPix, MCalibrationChargePINDiode,
57// MCalibrationQECam, MBadPixelsPix, MBadPixelsCam
58//
59// If the flag RelTimeCalibration() is set, a calibration of the relative arrival
60// times is also performed. The following containers (rectangular) and
61// tasks (ellipses) are called to produce an MCalibrationRelTimeCam used by
62// MCalibrateTime to correct timing offset between pixels: (MCalibrateTime is not
63// called from this class)
64//
65//Begin_Html
66/*
67<img src="images/RelTimeClasses.gif">
68*/
69//End_Html
70//
71// Different arrival time extractors can be chosen via the command SetArrivalTimeLevel(UInt_t i)
72// Up to now, the following extractors are available:
73// i=1: Use MArrivalTimeCalc (using MCubicSpline, arrival time == position at half maximum)
74// i=2: Use MArrivalTimeCalc2 (mean time of fWindowSize time slices with the highest integral content: default)
75//
76/////////////////////////////////////////////////////////////////////////////
77#include "MJCalibration.h"
78
79#include <TF1.h>
80#include <TFile.h>
81#include <TStyle.h>
82#include <TCanvas.h>
83#include <TSystem.h>
84#include <TProfile.h>
85
86#include "MLog.h"
87#include "MLogManip.h"
88
89#include "MRunIter.h"
90#include "MParList.h"
91#include "MTaskList.h"
92#include "MEvtLoop.h"
93
94#include "MHCamera.h"
95#include "MGeomCam.h"
96
97#include "MPedestalCam.h"
98#include "MCalibrationCam.h"
99#include "MCalibrationQECam.h"
100#include "MCalibrationChargeCam.h"
101#include "MCalibrationChargePINDiode.h"
102#include "MCalibrationChargeBlindPix.h"
103#include "MCalibrationChargeCalc.h"
104
105#include "MHGausEvents.h"
106#include "MHCalibrationCam.h"
107#include "MHCalibrationChargeCam.h"
108#include "MHCalibrationRelTimeCam.h"
109#include "MCalibrationRelTimeCam.h"
110
111#include "MReadMarsFile.h"
112#include "MGeomApply.h"
113#include "MBadPixelsMerge.h"
114#include "MBadPixelsCam.h"
115#include "MExtractSignal.h"
116#include "MExtractPINDiode.h"
117#include "MExtractBlindPixel.h"
118#include "MExtractSignal2.h"
119#include "MExtractSignal3.h"
120#include "MFCosmics.h"
121#include "MContinue.h"
122#include "MFillH.h"
123
124#include "MArrivalTimeCam.h"
125#include "MArrivalTimeCalc.h"
126#include "MArrivalTimeCalc2.h"
127
128#include "MJCalibration.h"
129#include "MStatusDisplay.h"
130
131ClassImp(MJCalibration);
132
133using namespace std;
134
135// --------------------------------------------------------------------------
136//
137// Default constructor.
138//
139// Sets fRuns to 0, fColor to kNONE, fDisplay to kNormalDisplay,
140// fRelTime to kFALSE, fExtractorLevel to 2, fArrivalTimeLevel to 2
141//
142MJCalibration::MJCalibration(const char *name, const char *title)
143 : fRuns(0), fColor(MCalibrationCam::kNONE), fDisplayType(kNormalDisplay),
144 fRelTimes(kFALSE), fExtractorLevel(2), fArrivalTimeLevel(2)
145{
146 fName = name ? name : "MJCalibration";
147 fTitle = title ? title : "Tool to create a pedestal file (MPedestalCam)";
148
149}
150
151// --------------------------------------------------------------------------
152//
153// Draw the MHCamera into the MStatusDisplay:
154//
155// 1) Draw it as histogram (MHCamera::DrawCopy("hist")
156// 2) Draw it as a camera, with MHCamera::SetPrettyPalette() set.
157// 3) If "rad" is not zero, draw its values vs. the radius from the camera center.
158// (DrawRadialProfile())
159// 4) Depending on the variable "fit", draw the values projection on the y-axis
160// (DrawProjection()):
161// 0: don't draw
162// 1: Draw fit to Single Gauss (for distributions flat-fielded over the whole camera)
163// 2: Draw and fit to Double Gauss (for distributions different for inner and outer pixels)
164// 3: Draw and fit to Triple Gauss (for distributions with inner, outer pixels and outliers)
165// 4: Draw and fit to Polynomial grade 0: (for the probability distributions)
166// >4: Draw and don;t fit.
167//
168void MJCalibration::CamDraw(TCanvas &c, const Int_t x, const Int_t y, const MHCamera &cam1,
169 const Int_t fit, const Int_t rad)
170{
171 c.cd(x);
172 gPad->SetBorderMode(0);
173 gPad->SetTicks();
174 MHCamera *obj1=(MHCamera*)cam1.DrawCopy("hist");
175 obj1->SetDirectory(NULL);
176
177 c.cd(x+y);
178 gPad->SetBorderMode(0);
179 obj1->SetPrettyPalette();
180 obj1->AddNotify(&fCalibrationCam);
181 obj1->Draw();
182
183 if (rad)
184 {
185 c.cd(x+2*y);
186 gPad->SetBorderMode(0);
187 gPad->SetTicks();
188 DrawRadialProfile(obj1);
189 }
190
191
192 if (!fit)
193 return;
194
195 c.cd(rad ? x+3*y : x+2*y);
196 gPad->SetBorderMode(0);
197 gPad->SetTicks();
198 DrawProjection(obj1, fit);
199}
200
201// --------------------------------------------------------------------------
202//
203// Display the results in MStatusDisplay:
204//
205// - Add "Calibration" to the MStatusDisplay title
206// - Retrieve the MGeomCam from MParList
207// - Initialize the following MHCamera's:
208// 1) MCalibrationPix::GetMean()
209// 2) MCalibrationPix::Sigma()
210// 3) MCalibrationChargePix::GetRSigma()
211// 4) MCalibrationChargePix::GetRSigmaPerCharge()
212// 5) MCalibrationChargePix::GetPheFFactorMethod()
213// 6) MCalibrationChargePix::GetMeanConvFADC2Phe()
214// 7) MCalibrationChargePix::GetMeanFFactorFADC2Phot()
215// 8) MCalibrationQEPix::GetQECascadesFFactor()
216// 9) MCalibrationQEPix::GetQECascadesBlindPixel()
217// 10) MCalibrationQEPix::GetQECascadesPINDiode()
218// 11) MCalibrationQEPix::GetQECascadesCombined()
219// 12) MCalibrationQEPix::IsAverageQEFFactorAvailable()
220// 13) MCalibrationQEPix::IsAverageQEBlindPixelAvailable()
221// 14) MCalibrationQEPix::IsAverageQEPINDiodeAvailable()
222// 15) MCalibrationQEPix::IsAverageQECombinedAvailable()
223// 16) MCalibrationChargePix::IsHiGainSaturation()
224// 17) MCalibrationPix::GetHiLoMeansDivided()
225// 18) MCalibrationPix::GetHiLoSigmasDivided()
226// 19) MCalibrationChargePix::GetHiGainPickup()
227// 20) MCalibrationChargePix::GetLoGainPickup()
228// 21) MCalibrationChargePix::GetHiGainBlackout()
229// 22) MCalibrationChargePix::GetLoGainBlackout()
230// 23) MCalibrationPix::IsExcluded()
231// 24) MBadPixelsPix::IsUnsuitable(MBadPixelsPix::kUnsuitableRun)
232// 25) MBadPixelsPix::IsUnsuitable(MBadPixelsPix::kUnreliableRun)
233// 26) MBadPixelsPix::IsUncalibrated(MBadPixelsPix::kHiGainOscillating)
234// 27) MBadPixelsPix::IsUncalibrated(MBadPixelsPix::kLoGainOscillating)
235// 28) MCalibrationChargePix::GetAbsTimeMean()
236// 29) MCalibrationChargePix::GetAbsTimeRms()
237//
238// If the flag SetFullDisplay() is set, all MHCameras will be displayed.
239// if the flag SetDataCheck() is set, only the most important ones are displayed
240// and otherwise, (default: SetNormalDisplay()), a good selection of plots is given
241//
242void MJCalibration::DisplayResult(MParList &plist)
243{
244 if (!fDisplay)
245 return;
246
247 //
248 // Update display
249 //
250 TString title = fDisplay->GetTitle();
251 title += "-- Calibration ";
252 title += fRuns->GetRunsAsString();
253 title += " --";
254 fDisplay->SetTitle(title);
255
256 //
257 // Get container from list
258 //
259 MGeomCam &geomcam = *(MGeomCam*)plist.FindObject("MGeomCam");
260
261 // Create histograms to display
262 MHCamera disp1 (geomcam, "Cal;Charge", "Fitted Mean Charges");
263 MHCamera disp2 (geomcam, "Cal;SigmaCharge", "Sigma of Fitted Charges");
264 MHCamera disp3 (geomcam, "Cal;RSigma", "Reduced Sigmas");
265 MHCamera disp4 (geomcam, "Cal;RSigma/Charge", "Reduced Sigma per Charge");
266 MHCamera disp5 (geomcam, "Cal;FFactorPhe", "Nr. of Phe's (F-Factor Method)");
267 MHCamera disp6 (geomcam, "Cal;FFactorConv", "Conversion Factor (F-Factor Method)");
268 MHCamera disp7 (geomcam, "Cal;FFactorFFactor", "Total F-Factor (F-Factor Method)");
269 MHCamera disp8 (geomcam, "Cal;CascadesQEFFactor", "Cascades QE (F-Factor Method)");
270 MHCamera disp9 (geomcam, "Cal;CascadesQEBlindPix","Cascades QE (Blind Pixel Method)");
271 MHCamera disp10(geomcam, "Cal;CascadesQEPINDiode","Cascades QE (PIN Diode Method)");
272 MHCamera disp11(geomcam, "Cal;CascadesQECombined","Cascades QE (Combined Method)");
273 MHCamera disp12(geomcam, "Cal;FFactorValid", "Pixels with valid F-Factor calibration");
274 MHCamera disp13(geomcam, "Cal;BlindPixelValid", "Pixels with valid BlindPixel calibration");
275 MHCamera disp14(geomcam, "Cal;PINdiodeValid", "Pixels with valid PINDiode calibration");
276 MHCamera disp15(geomcam, "Cal;CombinedValid", "Pixels with valid Combined calibration");
277 MHCamera disp16(geomcam, "Cal;Saturation", "Pixels with saturated Hi Gain");
278 MHCamera disp17(geomcam, "Cal;ConversionMeans", "Conversion HiGain.vs.LoGain Means");
279 MHCamera disp18(geomcam, "Cal;ConversionSigmas", "Conversion HiGain.vs.LoGain Sigmas");
280 MHCamera disp19(geomcam, "Cal;HiGainPickup", "Number Pickup events Hi Gain");
281 MHCamera disp20(geomcam, "Cal;LoGainPickup", "Number Pickup events Lo Gain");
282 MHCamera disp21(geomcam, "Cal;HiGainBlackout", "Number Blackout events Hi Gain");
283 MHCamera disp22(geomcam, "Cal;LoGainBlackout", "Number Blackout events Lo Gain");
284 MHCamera disp23(geomcam, "Cal;Excluded", "Pixels previously excluded");
285 MHCamera disp24(geomcam, "Bad;UnSuitable", "Pixels not suited for further analysis");
286 MHCamera disp25(geomcam, "Bad;UnReliable", "Pixels not reliable for further analysis");
287 MHCamera disp26(geomcam, "Bad;HiGainOscillating", "Oscillating Pixels High Gain");
288 MHCamera disp27(geomcam, "Bad;LoGainOscillating", "Oscillating Pixels Low Gain");
289 MHCamera disp28(geomcam, "Cal;AbsTimeMean", "Abs. Arrival Times");
290 MHCamera disp29(geomcam, "Cal;AbsTimeRms", "RMS of Arrival Times");
291 MHCamera disp30(geomcam, "time;MeanTime", "Mean Rel. Arrival Times");
292 MHCamera disp31(geomcam, "time;SigmaTime", "Sigma Rel. Arrival Times");
293 MHCamera disp32(geomcam, "time;TimeProb", "Probability of Time Fit");
294 MHCamera disp33(geomcam, "time;NotFitValid", "Pixels with not valid fit results");
295 MHCamera disp34(geomcam, "time;Oscillating", "Oscillating Pixels");
296
297
298
299 // Fitted charge means and sigmas
300 disp1.SetCamContent(fCalibrationCam, 0);
301 disp1.SetCamError( fCalibrationCam, 1);
302 disp2.SetCamContent(fCalibrationCam, 2);
303 disp2.SetCamError( fCalibrationCam, 3);
304
305 // Reduced Sigmas and reduced sigmas per charge
306 disp3.SetCamContent(fCalibrationCam, 5);
307 disp3.SetCamError( fCalibrationCam, 6);
308 disp4.SetCamContent(fCalibrationCam, 7);
309 disp4.SetCamError( fCalibrationCam, 8);
310
311 // F-Factor Method
312 disp5.SetCamContent(fCalibrationCam, 9);
313 disp5.SetCamError( fCalibrationCam, 10);
314 disp6.SetCamContent(fCalibrationCam, 11);
315 disp6.SetCamError( fCalibrationCam, 12);
316 disp7.SetCamContent(fCalibrationCam, 13);
317 disp7.SetCamError( fCalibrationCam, 14);
318
319 // Quantum Efficiencies
320 disp8.SetCamContent (fQECam, 0 );
321 disp8.SetCamError (fQECam, 1 );
322 disp9.SetCamContent (fQECam, 2 );
323 disp9.SetCamError (fQECam, 3 );
324 disp10.SetCamContent(fQECam, 4 );
325 disp10.SetCamError (fQECam, 5 );
326 disp11.SetCamContent(fQECam, 6 );
327 disp11.SetCamError (fQECam, 7 );
328
329 // Valid flags
330 disp12.SetCamContent(fQECam, 8 );
331 disp13.SetCamContent(fQECam, 9 );
332 disp14.SetCamContent(fQECam, 10);
333 disp15.SetCamContent(fQECam, 11);
334
335 // Conversion Hi-Lo
336 disp16.SetCamContent(fCalibrationCam, 25);
337 disp17.SetCamContent(fCalibrationCam, 16);
338 disp17.SetCamError (fCalibrationCam, 17);
339 disp18.SetCamContent(fCalibrationCam, 18);
340 disp18.SetCamError (fCalibrationCam, 19);
341
342 // Pickup and Blackout
343 disp19.SetCamContent(fCalibrationCam, 21);
344 disp20.SetCamContent(fCalibrationCam, 22);
345 disp21.SetCamContent(fCalibrationCam, 23);
346 disp22.SetCamContent(fCalibrationCam, 24);
347
348 // Pixels with defects
349 disp23.SetCamContent(fCalibrationCam, 20);
350 disp24.SetCamContent(fBadPixels, 1);
351 disp25.SetCamContent(fBadPixels, 3);
352
353 // Oscillations
354 disp26.SetCamContent(fBadPixels, 10);
355 disp27.SetCamContent(fBadPixels, 11);
356
357 // Arrival Times
358 disp28.SetCamContent(fCalibrationCam, 26);
359 disp28.SetCamError( fCalibrationCam, 27);
360 disp29.SetCamContent(fCalibrationCam, 27);
361
362 disp1.SetYTitle("Q [FADC counts]");
363 disp2.SetYTitle("\\sigma_{Q} [FADC counts]");
364
365 disp3.SetYTitle("\\sqrt{\\sigma^{2}_{Q} - RMS^{2}_{Ped}} [FADC Counts]");
366 disp4.SetYTitle("Red.Sigma/<Q> [1]");
367
368 disp5.SetYTitle("Nr. Phe's [1]");
369 disp6.SetYTitle("Conv.Factor [PhE/FADC counts]");
370 disp7.SetYTitle("Total F-Factor [1]");
371
372 disp8.SetYTitle("QE [1]");
373 disp9.SetYTitle("QE [1]");
374 disp10.SetYTitle("QE [1]");
375 disp11.SetYTitle("QE [1]");
376
377 disp12.SetYTitle("[1]");
378 disp13.SetYTitle("[1]");
379 disp14.SetYTitle("[1]");
380 disp15.SetYTitle("[1]");
381 disp16.SetYTitle("[1]");
382
383 disp17.SetYTitle("<Q>(High)/<Q>(Low) [1]");
384 disp18.SetYTitle("\\sigma_{Q}(High)/\\sigma_{Q}(Low) [1]");
385
386 disp19.SetYTitle("[1]");
387 disp20.SetYTitle("[1]");
388 disp21.SetYTitle("[1]");
389 disp22.SetYTitle("[1]");
390 disp23.SetYTitle("[1]");
391 disp24.SetYTitle("[1]");
392 disp25.SetYTitle("[1]");
393 disp26.SetYTitle("[1]");
394 disp27.SetYTitle("[1]");
395
396 disp28.SetYTitle("Mean Abs. Time [FADC slice]");
397 disp29.SetYTitle("RMS Abs. Time [FADC slices]");
398
399 if (fRelTimes)
400 {
401
402 disp30.SetCamContent(fRelTimeCam,0);
403 disp30.SetCamError( fRelTimeCam,1);
404 disp31.SetCamContent(fRelTimeCam,2);
405 disp31.SetCamError( fRelTimeCam,3);
406 disp32.SetCamContent(fRelTimeCam,4);
407 disp33.SetCamContent(fBadPixels,20);
408 disp34.SetCamContent(fBadPixels,21);
409
410 disp30.SetYTitle("Time Offset [ns]");
411 disp31.SetYTitle("Timing resolution [ns]");
412 disp32.SetYTitle("P_{Time} [1]");
413 disp33.SetYTitle("[1]");
414 disp34.SetYTitle("[1]");
415 }
416
417 gStyle->SetOptStat(1111);
418 gStyle->SetOptFit();
419
420 if (fDisplayType == kDataCheck)
421 {
422 TCanvas &c1 = fDisplay->AddTab("Fit.Charge");
423 c1.Divide(3, 3);
424
425 CamDraw(c1, 1, 3, disp1, 2);
426 CamDraw(c1, 2, 3, disp4, 2);
427 CamDraw(c1, 3, 3, disp28, 2);
428
429 // F-Factor
430 TCanvas &c2 = fDisplay->AddTab("Phe's");
431 c2.Divide(3,4);
432
433 CamDraw(c2, 1, 3, disp6, 2, 1);
434 CamDraw(c2, 2, 3, disp7, 2, 1);
435 CamDraw(c2, 3, 3, disp8, 2, 1);
436
437 // QE's
438 TCanvas &c3 = fDisplay->AddTab("QE's");
439 c3.Divide(3,4);
440
441 CamDraw(c3, 1, 3, disp8, 2, 1);
442 CamDraw(c3, 2, 3, disp9, 2, 1);
443 CamDraw(c3, 3, 3, disp10, 2, 1);
444
445 // Defects
446 TCanvas &c4 = fDisplay->AddTab("Defect");
447 c4.Divide(3,2);
448
449 CamDraw(c4, 1, 3, disp23, 0);
450 CamDraw(c4, 2, 3, disp24, 0);
451 CamDraw(c4, 3, 3, disp25, 0);
452
453 if (fRelTimes)
454 {
455 // Rel. Times
456 TCanvas &c5 = fDisplay->AddTab("Rel. Times");
457 c5.Divide(2,4);
458
459 CamDraw(c5, 1, 2, disp30, 2);
460 CamDraw(c5, 2, 2, disp31, 2);
461 }
462
463
464 return;
465 }
466
467 if (fDisplayType == kNormalDisplay)
468 {
469
470 // Charges
471 TCanvas &c11 = fDisplay->AddTab("Fit.Charge");
472 c11.Divide(2, 4);
473
474 CamDraw(c11, 1, 2, disp1, 2, 1);
475 CamDraw(c11, 2, 2, disp2, 2, 1);
476
477 // Reduced Sigmas
478 TCanvas &c12 = fDisplay->AddTab("Red.Sigma");
479 c12.Divide(2,4);
480
481 CamDraw(c12, 1, 2, disp3, 2, 1);
482 CamDraw(c12, 2, 2, disp4, 2, 1);
483
484 // F-Factor
485 TCanvas &c13 = fDisplay->AddTab("Phe's");
486 c13.Divide(3,4);
487
488 CamDraw(c13, 1, 3, disp5, 2, 1);
489 CamDraw(c13, 2, 3, disp6, 2, 1);
490 CamDraw(c13, 3, 3, disp7, 2, 1);
491
492 // QE's
493 TCanvas &c14 = fDisplay->AddTab("QE's");
494 c14.Divide(4,4);
495
496 CamDraw(c14, 1, 4, disp8, 2, 1);
497 CamDraw(c14, 2, 4, disp9, 2, 1);
498 CamDraw(c14, 3, 4, disp10, 2, 1);
499 CamDraw(c14, 4, 4, disp11, 2, 1);
500
501 // Defects
502 TCanvas &c15 = fDisplay->AddTab("Defect");
503 c15.Divide(5,2);
504
505 CamDraw(c15, 1, 5, disp23, 0);
506 CamDraw(c15, 2, 5, disp24, 0);
507 CamDraw(c15, 3, 5, disp25, 0);
508 CamDraw(c15, 4, 5, disp26, 0);
509 CamDraw(c15, 5, 5, disp27, 0);
510
511 // Abs. Times
512 TCanvas &c16 = fDisplay->AddTab("Abs. Times");
513 c16.Divide(2,3);
514
515 CamDraw(c16, 1, 2, disp28, 2);
516 CamDraw(c16, 2, 2, disp29, 1);
517
518 if (fRelTimes)
519 {
520 // Rel. Times
521 TCanvas &c17 = fDisplay->AddTab("Rel. Times");
522 c17.Divide(2,4);
523
524 CamDraw(c17, 1, 2, disp30, 2, 1);
525 CamDraw(c17, 2, 2, disp31, 2, 1);
526 }
527
528 return;
529 }
530
531 if (fDisplayType == kFullDisplay)
532 {
533
534 MHCalibrationCam *cam = (MHCalibrationCam*)plist.FindObject("MHCalibrationChargeCam");
535
536 for (Int_t aidx=0;aidx<cam->GetAverageAreas();aidx++)
537 {
538 cam->GetAverageHiGainArea(aidx).DrawClone("all");
539 cam->GetAverageLoGainArea(aidx).DrawClone("all");
540 }
541
542 for (Int_t sector=1;sector<cam->GetAverageSectors();sector++)
543 {
544 cam->GetAverageHiGainSector(sector).DrawClone("all");
545 cam->GetAverageLoGainSector(sector).DrawClone("all");
546 }
547
548 // Charges
549 TCanvas &c21 = fDisplay->AddTab("Fit.Charge");
550 c21.Divide(2, 4);
551
552 CamDraw(c21, 1, 2, disp1, 2, 1);
553 CamDraw(c21, 2, 2, disp2, 2, 1);
554
555 // Reduced Sigmas
556 TCanvas &c23 = fDisplay->AddTab("Red.Sigma");
557 c23.Divide(2,4);
558
559 CamDraw(c23, 1, 2, disp3, 2, 1);
560 CamDraw(c23, 2, 2, disp4, 2, 1);
561
562 // F-Factor
563 TCanvas &c24 = fDisplay->AddTab("Phe's");
564 c24.Divide(3,4);
565
566 CamDraw(c24, 1, 3, disp5, 2, 1);
567 CamDraw(c24, 2, 3, disp6, 2, 1);
568 CamDraw(c24, 3, 3, disp7, 2, 1);
569
570 // QE's
571 TCanvas &c25 = fDisplay->AddTab("QE's");
572 c25.Divide(4,4);
573
574 CamDraw(c25, 1, 4, disp8, 2, 1);
575 CamDraw(c25, 2, 4, disp9, 2, 1);
576 CamDraw(c25, 3, 4, disp10, 2, 1);
577 CamDraw(c25, 4, 4, disp11, 2, 1);
578
579 // Validity
580 TCanvas &c26 = fDisplay->AddTab("Valid");
581 c26.Divide(4,2);
582
583 CamDraw(c26, 1, 4, disp12, 0);
584 CamDraw(c26, 2, 4, disp13, 0);
585 CamDraw(c26, 3, 4, disp14, 0);
586 CamDraw(c26, 4, 4, disp15, 0);
587
588 // Other info
589 TCanvas &c27 = fDisplay->AddTab("HiLoGain");
590 c27.Divide(3,3);
591
592 CamDraw(c27, 1, 3, disp16, 0);
593 CamDraw(c27, 2, 3, disp17, 1);
594 CamDraw(c27, 3, 3, disp18, 1);
595
596 // Pickup
597 TCanvas &c28 = fDisplay->AddTab("Pickup");
598 c28.Divide(4,2);
599
600 CamDraw(c28, 1, 4, disp19, 0);
601 CamDraw(c28, 2, 4, disp20, 0);
602 CamDraw(c28, 3, 4, disp21, 0);
603 CamDraw(c28, 4, 4, disp22, 0);
604
605 // Defects
606 TCanvas &c29 = fDisplay->AddTab("Defect");
607 c29.Divide(5,2);
608
609 CamDraw(c29, 1, 5, disp23, 0);
610 CamDraw(c29, 2, 5, disp24, 0);
611 CamDraw(c29, 3, 5, disp25, 0);
612 CamDraw(c29, 4, 5, disp26, 0);
613 CamDraw(c29, 5, 5, disp27, 0);
614
615 // Abs. Times
616 TCanvas &c30 = fDisplay->AddTab("Abs. Times");
617 c30.Divide(2,3);
618
619 CamDraw(c30, 1, 2, disp28, 2);
620 CamDraw(c30, 2, 2, disp29, 1);
621
622 if (fRelTimes)
623 {
624 // Rel. Times
625 TCanvas &c31 = fDisplay->AddTab("Rel. Times");
626 c31.Divide(3,4);
627
628 CamDraw(c31, 1, 3, disp30, 2, 1);
629 CamDraw(c31, 2, 3, disp31, 2, 1);
630 CamDraw(c31, 3, 3, disp32, 4, 1);
631
632 // Time Defects
633 TCanvas &c32 = fDisplay->AddTab("Time Def.");
634 c32.Divide(2,2);
635
636 CamDraw(c32, 1, 2, disp33,0);
637 CamDraw(c32, 2, 2, disp34,0);
638
639 MHCalibrationCam *cam = (MHCalibrationCam*)plist.FindObject("MHCalibrationRelTimeCam");
640
641 for (Int_t aidx=0;aidx<cam->GetAverageAreas();aidx++)
642 {
643 cam->GetAverageHiGainArea(aidx).DrawClone("fourierevents");
644 cam->GetAverageLoGainArea(aidx).DrawClone("fourierevents");
645 }
646
647 for (Int_t sector=1;sector<cam->GetAverageSectors();sector++)
648 {
649 cam->GetAverageHiGainSector(sector).DrawClone("fourierevents");
650 cam->GetAverageLoGainSector(sector).DrawClone("fourierevents");
651 }
652
653 }
654
655 return;
656 }
657}
658
659// --------------------------------------------------------------------------
660//
661// Draw a projection of MHCamera onto the y-axis values. Depending on the
662// variable fit, the following fits are performed:
663//
664// 1: Single Gauss (for distributions flat-fielded over the whole camera)
665// 2: Double Gauss (for distributions different for inner and outer pixels)
666// 3: Triple Gauss (for distributions with inner, outer pixels and outliers)
667// 4: flat (for the probability distributions)
668//
669// Moreover, sectors 6,1 and 2 of the camera and sectors 3,4 and 5 are
670// drawn separately, for inner and outer pixels.
671//
672void MJCalibration::DrawProjection(MHCamera *obj, Int_t fit) const
673{
674
675 TH1D *obj2 = (TH1D*)obj->Projection(obj->GetName());
676 obj2->SetDirectory(0);
677 obj2->Draw();
678 obj2->SetBit(kCanDelete);
679
680 if (obj->GetGeomCam().InheritsFrom("MGeomCamMagic"))
681 {
682 TArrayI s0(3);
683 s0[0] = 6;
684 s0[1] = 1;
685 s0[2] = 2;
686
687 TArrayI s1(3);
688 s1[0] = 3;
689 s1[1] = 4;
690 s1[2] = 5;
691
692 TArrayI inner(1);
693 inner[0] = 0;
694
695 TArrayI outer(1);
696 outer[0] = 1;
697
698 // Just to get the right (maximum) binning
699 TH1D *half[4];
700 half[0] = obj->ProjectionS(s0, inner, "Sector 6-1-2 Inner");
701 half[1] = obj->ProjectionS(s1, inner, "Sector 3-4-5 Inner");
702 half[2] = obj->ProjectionS(s0, outer, "Sector 6-1-2 Outer");
703 half[3] = obj->ProjectionS(s1, outer, "Sector 3-4-5 Outer");
704
705 for (int i=0; i<4; i++)
706 {
707 half[i]->SetLineColor(kRed+i);
708 half[i]->SetDirectory(0);
709 half[i]->SetBit(kCanDelete);
710 half[i]->Draw("same");
711 }
712 }
713
714 const Double_t min = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst());
715 const Double_t max = obj2->GetBinCenter(obj2->GetXaxis()->GetLast());
716 const Double_t integ = obj2->Integral("width")/2.5;
717 const Double_t mean = obj2->GetMean();
718 const Double_t rms = obj2->GetRMS();
719 const Double_t width = max-min;
720
721 const TString dgausformula = "([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])"
722 "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
723
724 const TString tgausformula = "([0]-[3]-[6])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])"
725 "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])"
726 "+[6]/[8]*exp(-0.5*(x-[7])*(x-[7])/[8]/[8])";
727 TF1 *f=0;
728 switch (fit)
729 {
730 case 1:
731 f = new TF1("sgaus", "gaus(0)", min, max);
732 f->SetLineColor(kYellow);
733 f->SetBit(kCanDelete);
734 f->SetParNames("Area", "#mu", "#sigma");
735 f->SetParameters(integ/rms, mean, rms);
736 f->SetParLimits(0, 0, integ);
737 f->SetParLimits(1, min, max);
738 f->SetParLimits(2, 0, width/1.5);
739
740 obj2->Fit(f, "QLR");
741 break;
742
743 case 2:
744 f = new TF1("dgaus",dgausformula.Data(),min,max);
745 f->SetLineColor(kYellow);
746 f->SetBit(kCanDelete);
747 f->SetParNames("A_{tot}", "#mu1", "#sigma1", "A2", "#mu2", "#sigma2");
748 f->SetParameters(integ,(min+mean)/2.,width/4.,
749 integ/width/2.,(max+mean)/2.,width/4.);
750 // The left-sided Gauss
751 f->SetParLimits(0,integ-1.5 , integ+1.5);
752 f->SetParLimits(1,min+(width/10.), mean);
753 f->SetParLimits(2,0 , width/2.);
754 // The right-sided Gauss
755 f->SetParLimits(3,0 , integ);
756 f->SetParLimits(4,mean, max-(width/10.));
757 f->SetParLimits(5,0 , width/2.);
758 obj2->Fit(f,"QLRM");
759 break;
760
761 case 3:
762 f = new TF1("tgaus",tgausformula.Data(),min,max);
763 f->SetLineColor(kYellow);
764 f->SetBit(kCanDelete);
765 f->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}",
766 "A_{2}","#mu_{2}","#sigma_{2}",
767 "A_{3}","#mu_{3}","#sigma_{3}");
768 f->SetParameters(integ,(min+mean)/2,width/4.,
769 integ/width/3.,(max+mean)/2.,width/4.,
770 integ/width/3.,mean,width/2.);
771 // The left-sided Gauss
772 f->SetParLimits(0,integ-1.5,integ+1.5);
773 f->SetParLimits(1,min+(width/10.),mean);
774 f->SetParLimits(2,width/15.,width/2.);
775 // The right-sided Gauss
776 f->SetParLimits(3,0.,integ);
777 f->SetParLimits(4,mean,max-(width/10.));
778 f->SetParLimits(5,width/15.,width/2.);
779 // The Gauss describing the outliers
780 f->SetParLimits(6,0.,integ);
781 f->SetParLimits(7,min,max);
782 f->SetParLimits(8,width/4.,width/1.5);
783 obj2->Fit(f,"QLRM");
784 break;
785
786 case 4:
787 obj2->Fit("pol0", "Q");
788 obj2->GetFunction("pol0")->SetLineColor(kYellow);
789 break;
790
791 case 9:
792 break;
793
794 default:
795 obj2->Fit("gaus", "Q");
796 obj2->GetFunction("gaus")->SetLineColor(kYellow);
797 break;
798 }
799}
800
801// --------------------------------------------------------------------------
802//
803// Draw a projection of MHCamera vs. the radius from the central pixel.
804//
805// The inner and outer pixels are drawn separately, both fitted by a polynomial
806// of grade 1.
807//
808void MJCalibration::DrawRadialProfile(MHCamera *obj) const
809{
810
811 TProfile *obj2 = (TProfile*)obj->RadialProfile(obj->GetName());
812 obj2->SetDirectory(0);
813 obj2->Draw();
814 obj2->SetBit(kCanDelete);
815
816 if (obj->GetGeomCam().InheritsFrom("MGeomCamMagic"))
817 {
818
819 TArrayI s0(6);
820 s0[0] = 1;
821 s0[1] = 2;
822 s0[2] = 3;
823 s0[3] = 4;
824 s0[4] = 5;
825 s0[5] = 6;
826
827 TArrayI inner(1);
828 inner[0] = 0;
829
830 TArrayI outer(1);
831 outer[0] = 1;
832
833 // Just to get the right (maximum) binning
834 TProfile *half[2];
835 half[0] = obj->RadialProfileS(s0, inner,Form("%s%s",obj->GetName(),"Inner"));
836 half[1] = obj->RadialProfileS(s0, outer,Form("%s%s",obj->GetName(),"Outer"));
837
838 for (Int_t i=0; i<2; i++)
839 {
840 Double_t min = obj->GetGeomCam().GetMinRadius(i);
841 Double_t max = obj->GetGeomCam().GetMaxRadius(i);
842
843 half[i]->SetLineColor(kRed+i);
844 half[i]->SetDirectory(0);
845 half[i]->SetBit(kCanDelete);
846 half[i]->Draw("same");
847 half[i]->Fit("pol1","Q","",min,max);
848 half[i]->GetFunction("pol1")->SetLineColor(kRed+i);
849 half[i]->GetFunction("pol1")->SetLineWidth(1);
850 }
851 }
852}
853
854
855// --------------------------------------------------------------------------
856//
857// Retrieve the output file written by WriteResult()
858//
859TString MJCalibration::GetOutputFile() const
860{
861 if (!fRuns)
862 return "";
863
864 return Form("%s/%s-F1.root", (const char*)fOutputPath, (const char*)fRuns->GetRunsAsFileName());
865}
866
867
868// --------------------------------------------------------------------------
869//
870// Call the ProcessFile(MPedestalCam)
871//
872Bool_t MJCalibration::Process(MPedestalCam &pedcam)
873{
874 if (!ReadCalibrationCam())
875 return ProcessFile(pedcam);
876
877 return kTRUE;
878}
879
880// --------------------------------------------------------------------------
881//
882// Execute the task list and the eventloop:
883//
884// - Check if there are fRuns, otherwise return
885// - Check for consistency between run numbers and number of files
886// - Add fRuns to MReadMarsFile
887// - Put into MParList:
888// 1) MPedestalCam (pedcam)
889// 2) MCalibrationQECam (fQECam)
890// 3) MCalibrationChargeCam (fCalibrationCam)
891// 4) MCalibrationRelTimeCam (fRelTimeCam) (only if flag fRelTimes is chosen)
892// 5) MBadPixelsCam (fBadPixels)
893// 6) MCalibrationChargePINDiode
894// 7) MCalibrationChargeBlindPix
895// - Put into the MTaskList:
896// 1) MReadMarsFile
897// 2) MBadPixelsMerge
898// 3) MGeomApply
899// 4) MExtractSignal, MExtractSignal2 or MExtractSignal3, depending on fExtractorLevel
900// 5) MExtractPINDiode
901// 6) MExtractBlindPixel
902// 7) MArrivalTimeCalc, MArrivalTimeCalc2, depending on fArrivalTimeLevel (only if flag fRelTimes is chosen)
903// 8) MContinue(MFCosmics)
904// 9) MFillH("MHCalibrationChargePINDiode", "MExtractedSignalPINDiode")
905// 10) MFillH("MHCalibrationChargeBlindPix", "MExtractedSignalBlindPixel")
906// 11) MFillH("MHCalibrationChargeCam", "MExtractedSignalCam")
907// 12) MFillH("MHCalibrationChargeCam", "MExtractedSignalCam")
908// 13) MCalibrationChargeCalc
909// 14) MFillH("MHCalibrationRelTimeCam", "MArrivalTimeCam") (only if flag fRelTimes is chosen)
910// - Execute MEvtLoop
911// - DisplayResult()
912// - WriteResult()
913//
914Bool_t MJCalibration::ProcessFile(MPedestalCam &pedcam)
915{
916 if (!fRuns)
917 {
918 *fLog << err << "No Runs choosen... abort." << endl;
919 return kFALSE;
920 }
921
922 if (fRuns->GetNumRuns() != fRuns->GetNumEntries())
923 {
924 *fLog << err << "Number of files found doesn't match number of runs... abort." << endl;
925 return kFALSE;
926 }
927
928 *fLog << inf;
929 fLog->Separator(GetDescriptor());
930 *fLog << "Calculate MCalibrationCam from Runs " << fRuns->GetRunsAsString() << endl;
931 *fLog << endl;
932
933 MReadMarsFile read("Events");
934 read.DisableAutoScheme();
935 static_cast<MRead&>(read).AddFiles(*fRuns);
936
937 MCalibrationChargePINDiode pindiode;
938 MCalibrationChargeBlindPix blindpix;
939
940 // Setup Tasklist
941 MParList plist;
942 plist.AddToList(&pedcam);
943 plist.AddToList(&fBadPixels);
944 plist.AddToList(&fQECam);
945 plist.AddToList(&fCalibrationCam);
946 plist.AddToList(&pindiode);
947 plist.AddToList(&blindpix);
948
949 MTaskList tlist;
950 plist.AddToList(&tlist);
951
952 MGeomApply apply;
953 // MBadPixelsMerge merge(&fBadPixels);
954 MExtractPINDiode pinext;
955 MExtractBlindPixel blindext;
956 MExtractSignal extract1;
957 MExtractSignal2 extract2;
958 MExtractSignal3 extract3;
959 MArrivalTimeCalc tmecalc1;
960 MArrivalTimeCalc2 tmecalc2;
961 MCalibrationChargeCalc calcalc;
962
963 //
964 // As long as there are no DM's, have to colour by hand
965 //
966 calcalc.SetPulserColor(fColor);
967
968 MFillH fillpin("MHCalibrationChargePINDiode", "MExtractedSignalPINDiode");
969 MFillH fillbnd("MHCalibrationChargeBlindPix", "MExtractedSignalBlindPixel");
970 MFillH fillcam("MHCalibrationChargeCam", "MExtractedSignalCam");
971 MFillH filltme("MHCalibrationRelTimeCam", "MArrivalTimeCam");
972 fillpin.SetNameTab("PINDiode");
973 fillbnd.SetNameTab("BlindPix");
974 fillcam.SetNameTab("Charge");
975 filltme.SetNameTab("RelTimes");
976
977
978 //
979 // Apply a filter against cosmics
980 // (will have to be needed in the future
981 // when the calibration hardware-trigger is working)
982 //
983 MFCosmics cosmics;
984 MContinue cont(&cosmics);
985
986 tlist.AddToList(&read);
987 // tlist.AddToList(&merge);
988 tlist.AddToList(&apply);
989
990 if (fExtractorLevel <= 1)
991 tlist.AddToList(&extract1);
992 else if (fExtractorLevel == 2)
993 tlist.AddToList(&extract2);
994 else if (fExtractorLevel == 3)
995 tlist.AddToList(&extract3);
996 else
997 {
998 *fLog << err << GetDescriptor()
999 << ": No valid Signal extractor has been chosen, have only: " << fExtractorLevel
1000 << " aborting..." << endl;
1001 return kFALSE;
1002 }
1003
1004
1005 tlist.AddToList(&pinext);
1006 tlist.AddToList(&blindext);
1007
1008 if (fRelTimes)
1009 {
1010 plist.AddToList(&fRelTimeCam);
1011 if (fArrivalTimeLevel <= 1)
1012 tlist.AddToList(&tmecalc1);
1013 else if (fArrivalTimeLevel == 2)
1014 tlist.AddToList(&tmecalc2);
1015 else
1016 {
1017 *fLog << err << GetDescriptor()
1018 << ": No valid Signal ArrivalTime Level has been chosen, have only: " << fArrivalTimeLevel
1019 << " aborting..." << endl;
1020 return kFALSE;
1021 }
1022 }
1023
1024 tlist.AddToList(&cont);
1025 tlist.AddToList(&fillcam);
1026 tlist.AddToList(&fillpin);
1027 tlist.AddToList(&fillbnd);
1028 tlist.AddToList(&calcalc);
1029
1030 if (fRelTimes)
1031 tlist.AddToList(&filltme);
1032
1033 // Create and setup the eventloop
1034 MEvtLoop evtloop(fName);
1035 evtloop.SetParList(&plist);
1036 evtloop.SetDisplay(fDisplay);
1037 evtloop.SetLogStream(fLog);
1038
1039 // Execute first analysis
1040 if (!evtloop.Eventloop())
1041 {
1042 *fLog << err << GetDescriptor() << ": Failed." << endl;
1043 return kFALSE;
1044 }
1045
1046 tlist.PrintStatistics();
1047
1048 DisplayResult(plist);
1049
1050 if (!WriteResult())
1051 return kFALSE;
1052
1053 *fLog << inf << GetDescriptor() << ": Done." << endl;
1054
1055 return kTRUE;
1056}
1057
1058// --------------------------------------------------------------------------
1059//
1060// Read the following containers from GetOutputFile()
1061// - MCalibrationChargeCam
1062// - MCalibrationQECam
1063// - MBadPixelsCam
1064//
1065Bool_t MJCalibration::ReadCalibrationCam()
1066{
1067 const TString fname = GetOutputFile();
1068
1069 if (gSystem->AccessPathName(fname, kFileExists))
1070 {
1071 *fLog << err << "Input file " << fname << " doesn't exist." << endl;
1072 return kFALSE;
1073 }
1074
1075 *fLog << inf << "Reading from file: " << fname << endl;
1076
1077 TFile file(fname, "READ");
1078 if (fCalibrationCam.Read()<=0)
1079 {
1080 *fLog << err << "Unable to read MCalibrationChargeCam from " << fname << endl;
1081 return kFALSE;
1082 }
1083
1084 if (fQECam.Read()<=0)
1085 {
1086 *fLog << err << "Unable to read MCalibrationQECam from " << fname << endl;
1087 return kFALSE;
1088 }
1089
1090 if (file.FindKey("MBadPixelsCam"))
1091 {
1092 MBadPixelsCam bad;
1093 if (bad.Read()<=0)
1094 {
1095 *fLog << err << "Unable to read MBadPixelsCam from " << fname << endl;
1096 return kFALSE;
1097 }
1098 fBadPixels.Merge(bad);
1099 }
1100
1101 if (fDisplay /*&& !fDisplay->GetCanvas("Pedestals")*/) // FIXME!
1102 fDisplay->Read();
1103
1104 return kTRUE;
1105}
1106
1107
1108// --------------------------------------------------------------------------
1109//
1110// Set the path for output files, written by WriteResult()
1111//
1112void MJCalibration::SetOutputPath(const char *path)
1113{
1114 fOutputPath = path;
1115 if (fOutputPath.EndsWith("/"))
1116 fOutputPath = fOutputPath(0, fOutputPath.Length()-1);
1117}
1118
1119
1120// --------------------------------------------------------------------------
1121//
1122// Write the result into the output file GetOutputFile(), if fOutputPath exists.
1123//
1124// The following containers are written:
1125// - MStatusDisplay
1126// - MCalibrationChargeCam
1127// - MCalibrationQECam
1128// - MBadPixelsCam
1129//
1130Bool_t MJCalibration::WriteResult()
1131{
1132 if (fOutputPath.IsNull())
1133 return kTRUE;
1134
1135 const TString oname(GetOutputFile());
1136
1137 *fLog << inf << "Writing to file: " << oname << endl;
1138
1139 TFile file(oname, "UPDATE");
1140
1141 if (fDisplay && fDisplay->Write()<=0)
1142 {
1143 *fLog << err << "Unable to write MStatusDisplay to " << oname << endl;
1144 return kFALSE;
1145 }
1146
1147 if (fCalibrationCam.Write()<=0)
1148 {
1149 *fLog << err << "Unable to write MCalibrationChargeCam to " << oname << endl;
1150 return kFALSE;
1151 }
1152
1153 if (fQECam.Write()<=0)
1154 {
1155 *fLog << err << "Unable to write MCalibrationQECam to " << oname << endl;
1156 return kFALSE;
1157 }
1158
1159 if (fBadPixels.Write()<=0)
1160 {
1161 *fLog << err << "Unable to write MBadPixelsCam to " << oname << endl;
1162 return kFALSE;
1163 }
1164
1165 return kTRUE;
1166
1167}
1168
Note: See TracBrowser for help on using the repository browser.