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

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