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

Last change on this file since 9456 was 9317, checked in by tbretz, 16 years ago
*** empty log message ***
File size: 55.3 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! Author(s): Markus Gaug, 02/2004 <mailto:markus@ifae.es>
20!
21! Copyright: MAGIC Software Development, 2000-2009
22!
23!
24\* ======================================================================== */
25
26/////////////////////////////////////////////////////////////////////////////
27//
28// MJCalibration
29//
30// Do one calibration loop over serious of runs with the same pulser
31// colour and the same intensity. The following containers (rectangular) and
32// tasks (ellipses) are called to produce an MCalibrationChargeCam and to
33// update the MCalibrationQECam: (MCalibrate is not called from this class)
34//
35//Begin_Html
36/*
37<img src="images/CalibClasses.gif">
38*/
39//End_Html
40//
41// Different signal extractors can be set with the command SetExtractor()
42// Only extractors deriving from MExtractor can be set, default is MExtractSlidingWindow
43//
44// Different arrival time extractors can be set with the command SetTimeExtractor()
45// Only extractors deriving from MExtractTime can be set, default is MExtractTimeHighestIntegral
46//
47// At the end of the eventloop, plots and results are displayed, depending on
48// the flags set (see DisplayResult())
49//
50// If the flag SetFullDisplay() is set, all MHCameras will be displayed.
51// if the flag SetDataCheckDisplay() is set, only the most important ones are displayed
52// Otherwise, (default: SetNormalDisplay()), a good selection of plots is given
53//
54// The absolute light calibration devices Blind Pixel and PIN Diode can be switched on
55// and off with the commands:
56//
57// - SetUseBlindPixel(Bool_t )
58// - SetUsePINDiode(Bool_t )
59//
60// See also: MHCalibrationChargePix, MHCalibrationChargeCam, MHGausEvents
61// MHCalibrationChargeBlindPix, MHCalibrationChargePINDiode
62// MCalibrationChargePix, MCalibrationChargeCam, MCalibrationChargeCalc
63// MCalibrationBlindPix, MCalibrationChargePINDiode,
64// MCalibrationQECam, MBadPixelsPix, MBadPixelsCam
65//
66// If the flag RelTimeCalibration() is set, a calibration of the relative arrival
67// times is also performed. The following containers (rectangular) and
68// tasks (ellipses) are called to produce an MCalibrationRelTimeCam used by
69// MCalibrateTime to correct timing offset between pixels: (MCalibrateTime is not
70// called from this class)
71//
72//Begin_Html
73/*
74<img src="images/RelTimeClasses.gif">
75*/
76//End_Html
77//
78// Different arrival time extractors can be set directly with the command
79// SetTimeExtractor(MExtractor *)
80//
81// Resource file entries are case sensitive!
82//
83// See also: MHCalibrationRelTimePix, MHCalibrationRelTimeCam, MHGausEvents
84// MCalibrationRelTimePix, MCalibrationRelTimeCam
85// MBadPixelsPix, MBadPixelsCam
86//
87/////////////////////////////////////////////////////////////////////////////
88#include "MJCalibration.h"
89
90#include <TFile.h>
91#include <TF1.h>
92#include <TStyle.h>
93#include <TCanvas.h>
94#include <TSystem.h>
95#include <TLine.h>
96#include <TLatex.h>
97#include <TLegend.h>
98#include <TRegexp.h>
99#include <TPaveText.h>
100#include <TPaveStats.h>
101#include <TEnv.h>
102
103#include "MLog.h"
104#include "MLogManip.h"
105
106#include "MEnv.h"
107#include "MString.h"
108#include "MDirIter.h"
109#include "MSequence.h"
110#include "MParList.h"
111#include "MTaskList.h"
112#include "MEvtLoop.h"
113
114#include "MHCamera.h"
115#include "MGeomCam.h"
116
117#include "MCalibrationPatternDecode.h"
118#include "MCalibrationCam.h"
119#include "MCalibrationQECam.h"
120#include "MCalibrationQEPix.h"
121#include "MCalibrationChargeCam.h"
122#include "MCalibrationChargePix.h"
123#include "MCalibrationChargePINDiode.h"
124#include "MCalibrationBlindPix.h"
125#include "MCalibrationBlindCam.h"
126#include "MCalibrationBlindCamOneOldStyle.h"
127#include "MCalibrationBlindCamTwoNewStyle.h"
128#include "MCalibrationBlindCamThreeNewStyle.h"
129#include "MCalibrationChargeCalc.h"
130#include "MCalibColorSet.h"
131#include "MCalibrationRelTimeCam.h"
132#include "MCalibrationRelTimeCalc.h"
133
134#include "MHGausEvents.h"
135#include "MHCalibrationCam.h"
136#include "MHCalibrationChargeCam.h"
137#include "MHCalibrationChargeBlindCam.h"
138#include "MHCalibrationChargePINDiode.h"
139#include "MHCalibrationRelTimeCam.h"
140#include "MHCalibrationPix.h"
141
142#include "MHCamEvent.h"
143
144#include "MReadMarsFile.h"
145#include "MPedCalcPedRun.h"
146#include "MRawFileRead.h"
147#include "MGeomApply.h"
148#include "MPedestalSubtract.h"
149#include "MTaskEnv.h"
150#include "MBadPixelsMerge.h"
151#include "MBadPixelsCam.h"
152#include "MExtractTime.h"
153#include "MExtractor.h"
154#include "MExtractPINDiode.h"
155#include "MExtractBlindPixel.h"
156#include "MExtractTimeAndChargeSpline.h"
157#include "MFCosmics.h"
158#include "MFTriggerPattern.h"
159#include "MContinue.h"
160#include "MFillH.h"
161
162#include "MTriggerPatternDecode.h"
163
164#include "MArrivalTimeCam.h"
165
166#include "MStatusDisplay.h"
167
168ClassImp(MJCalibration);
169
170using namespace std;
171
172const Int_t MJCalibration::gkIFAEBoxInaugurationRun = 20113;
173const Int_t MJCalibration::gkSecondBlindPixelInstallation = 31693;
174const Int_t MJCalibration::gkSpecialPixelsContInstallation = 34057;
175const Int_t MJCalibration::gkThirdBlindPixelInstallation = 43308;
176const TString MJCalibration::fgReferenceFile = "mjobs/calibrationref.rc";
177const TString MJCalibration::fgHiLoCalibFile = "resources/hilocalib.rc";
178
179// --------------------------------------------------------------------------
180//
181// Default constructor.
182//
183// - fExtractor to NULL, fTimeExtractor to NULL, fColor to kNONE,
184// fDisplay to kNormalDisplay, kRelTimes to kFALSE, kataCheck to kFALSE, kDebug to kFALSE
185// - SetUseBlindPixel()
186// - SetUsePINDiode()
187//
188MJCalibration::MJCalibration(const char *name, const char *title)
189 : fExtractor(NULL), fTimeExtractor(NULL),
190 fColor(MCalibrationCam::kNONE), fDisplayType(kDataCheckDisplay),
191 fMinEvents(1000), fGeometry("MGeomCamMagic")
192{
193
194 fName = name ? name : "MJCalibration";
195 fTitle = title ? title : "Tool to create the calibration constants for one calibration run";
196
197 //SetHiLoCalibration();
198 SetRelTimeCalibration();
199 SetDebug(kFALSE);
200
201 SetReferenceFile();
202 SetHiLoCalibFile();
203
204 fConvFADC2PheMin = 0.;
205 fConvFADC2PheMax = 1.5;
206 fConvFADC2PhotMin = 0.;
207 fConvFADC2PhotMax = 10.;
208 fQEMin = 0.;
209 fQEMax = 0.3;
210 fArrivalTimeMin = 1.;
211 fArrivalTimeMax = 10.;
212 fTimeOffsetMin = -3.;
213 fTimeOffsetMax = 3.;
214 fTimeResolutionMin = 0.;
215 fTimeResolutionMax = 1.;
216
217 fRefFADC2PheInner = 0.14;
218 fRefFADC2PheOuter = 0.4;
219 fRefConvFADC2PheInner = 0.14;
220 fRefConvFADC2PheOuter = 0.52;
221 fRefQEInner = 0.18;
222 fRefQEOuter = 0.12;
223 fRefArrivalTimeInner = 4.5;
224 fRefArrivalTimeOuter = 5.0;
225 fRefArrivalTimeRmsInner = 0.5;
226 fRefArrivalTimeRmsOuter = 0.5;
227 fRefTimeOffsetOuter = 0.62;
228 fRefTimeResolutionInner = 0.12;
229 fRefTimeResolutionOuter = 0.09;
230}
231
232void MJCalibration::DrawTab(MParList &plist, const char *cont, const char *name, Option_t *opt)
233{
234 TObject *obj = plist.FindObject(cont);
235 if (!obj)
236 return;
237
238 fDisplay->AddTab(name);
239 obj->DrawClone(opt);
240}
241
242MHCamera *MJCalibration::DrawBadPixelPad(const MHCamera &h, Bool_t unsuit) const
243{
244 MHCamera *obj=(MHCamera*)h.DrawCopy("hist");
245
246 gStyle->SetPalette(1);
247
248 const Int_t numcol = gStyle->GetNumberOfColors();
249
250 const Double_t min = 1;
251 const Double_t max = unsuit ? MBadPixelsPix::GetNumUnsuitable() : MBadPixelsPix::GetNumUnreliable();
252 const Double_t f = (numcol-1)/(max-min);
253
254 FixDataCheckHist(*obj, min, max);
255
256 TPaveText *pave = new TPaveText(0.05, 0.012, 0.975, 0.999);
257
258 const Double_t height = (pave->GetY2()-pave->GetY1())/(max+1);
259
260 pave->SetBit(kCanDelete);
261 pave->ConvertNDCtoPad();
262 pave->SetFillColor(14);
263 pave->Draw();
264
265 Int_t n=0;
266 while (1)
267 {
268 const TString name = unsuit ? MBadPixelsPix::GetUnsuitableName(++n) : MBadPixelsPix::GetUnreliableName(++n);
269 if (name.IsNull())
270 break;
271
272 Int_t cnt = 0;
273 for (UInt_t pix=0; pix<h.GetNumPixels(); pix++)
274 if (TMath::Nint(h.GetPixContent(pix)) == n)
275 cnt++;
276
277 const TString loc = unsuit?MBadPixelsPix::GetUnsuitableName(n):MBadPixelsPix::GetUnreliableName(n);
278
279 const TString left = Form("%d) %s", n, loc.Data());
280 const TString right = Form("%3i pixels", cnt);
281
282 const Int_t col = gStyle->GetColorPalette(TMath::FloorNint((n-1)*f));
283
284 TText *p = pave->AddText(0.05, pave->GetY2()-height*(n+0.3), left);
285 p->SetTextColor(col);
286 p->SetTextAlign(12);
287
288 p = pave->AddText(0.95, p->GetY(), right);
289 p->SetTextColor(col);
290 p->SetTextAlign(32);
291 }
292
293 return obj;
294}
295
296// --------------------------------------------------------------------------
297//
298// Display the results in MStatusDisplay:
299//
300// - Add "Calibration" to the MStatusDisplay title
301// - Retrieve the MGeomCam from MParList
302// - Initialize the following MHCamera's:
303// 1) MCalibrationPix::GetMean()
304// 2) MCalibrationPix::Sigma()
305// 3) MCalibrationChargePix::GetRSigma()
306// 4) MCalibrationChargePix::GetRSigmaPerCharge()
307// 5) MCalibrationChargePix::GetPheFFactorMethod()
308// 6) MCalibrationChargePix::GetMeanConvFADC2Phe()
309// 7) MCalibrationChargePix::GetMeanFFactorFADC2Phot()
310// 8) MCalibrationQEPix::GetQECascadesFFactor()
311// 9) MCalibrationQEPix::GetQECascadesBlindPixel()
312// 10) MCalibrationQEPix::GetQECascadesPINDiode()
313// 11) MCalibrationQEPix::GetQECascadesCombined()
314// 12) MCalibrationQEPix::IsAverageQEFFactorAvailable()
315// 13) MCalibrationQEPix::IsAverageQEBlindPixelAvailable()
316// 14) MCalibrationQEPix::IsAverageQEPINDiodeAvailable()
317// 15) MCalibrationQEPix::IsAverageQECombinedAvailable()
318// 16) MCalibrationChargePix::IsHiGainSaturation()
319// 17) MCalibrationPix::GetHiLoMeansDivided()
320// 18) MCalibrationPix::GetHiLoSigmasDivided()
321// 19) MCalibrationChargePix::GetHiGainPickup()
322// 20) MCalibrationChargePix::GetLoGainPickup()
323// 21) MCalibrationChargePix::GetHiGainBlackout()
324// 22) MCalibrationChargePix::GetLoGainBlackout()
325// 23) MCalibrationPix::IsExcluded()
326// 24) MBadPixelsPix::IsUnsuitable(MBadPixelsPix::kUnsuitableRun)
327// 25) MBadPixelsPix::IsUnsuitable(MBadPixelsPix::kUnreliableRun)
328// 26) MBadPixelsPix::IsUncalibrated(MBadPixelsPix::kHiGainOscillating)
329// 27) MBadPixelsPix::IsUncalibrated(MBadPixelsPix::kLoGainOscillating)
330// 28) MCalibrationChargePix::GetAbsTimeMean()
331// 29) MCalibrationChargePix::GetAbsTimeRms()
332//
333// If the flag SetFullDisplay() is set, all MHCameras will be displayed.
334// if the flag SetDataCheckDisplay() is set, only the most important ones are displayed
335// and otherwise, (default: SetNormalDisplay()), a good selection of plots is given
336//
337void MJCalibration::DisplayResult(MParList &plist)
338{
339 if (!fDisplay)
340 return;
341
342 TString drawoption = "nonew ";
343 if (fDisplayType == kDataCheckDisplay)
344 drawoption += "datacheck";
345 if (fDisplayType == kFullDisplay)
346 drawoption += "all";
347
348 if (IsUsePINDiode())
349 DrawTab(plist, "MHCalibrationChargePINDiode", "PINDiode", drawoption);
350 if (IsUseBlindPixel())
351 DrawTab(plist, "MHCalibrationChargeBlindCam", "BlindPix", drawoption);
352 if (IsRelTimes())
353 DrawTab(plist, "MHCalibrationRelTimeCam", "Time", drawoption);
354 DrawTab(plist, "MHCalibrationChargeCam", "Charge", drawoption);
355
356 //
357 // Update display
358 //
359 TString title = "-- Calibration: ";
360 title += fSequence.GetSequence();
361 title += " --";
362 fDisplay->SetTitle(title, kFALSE);
363
364 //
365 // Get container from list
366 //
367 MGeomCam &geomcam = *(MGeomCam*)plist.FindObject("MGeomCam");
368
369 // Create histograms to display
370 MHCamera disp1 (geomcam, "Charge", "Fitted Mean Signal (Charges)");
371 MHCamera disp2 (geomcam, "SigmaCharge", "Sigma of Fitted Signal");
372 MHCamera disp3 (geomcam, "RSigma", "Reduced Sigmas");
373 MHCamera disp4 (geomcam, "RSigmaPerCharge", "Reduced Sigma per Charge");
374 MHCamera disp5 (geomcam, "NumPhes", "Number Photo-electrons");
375 MHCamera disp6 (geomcam, "ConvFADC2Phes", "Phes per Charge (Before Flat-Field)");
376 MHCamera disp7 (geomcam, "TotalFFactor", "Total F-Factor(F-Factor Method)");
377 MHCamera disp8 (geomcam, "CascadesQEFFactor", "Cascades QE (F-Factor Method)");
378 MHCamera disp9 (geomcam, "CascadesQEBlindPix","Cascades QE (Blind Pixel Method)");
379 MHCamera disp10(geomcam, "CascadesQEPINDiode","Cascades QE (PIN Diode Method)");
380 MHCamera disp11(geomcam, "CascadesQECombined","Cascades QE (Combined Method)");
381 MHCamera disp12(geomcam, "FFactorValid", "Pixels with Valid F-Factor Calibration");
382 MHCamera disp13(geomcam, "BlindPixelValid", "Pixels with valid BlindPixel Calibration");
383 MHCamera disp14(geomcam, "PINdiodeValid", "Pixels with Valid PINDiode Calibration");
384 MHCamera disp15(geomcam, "CombinedValid", "Pixels with Valid Combined Calibration");
385 MHCamera disp16(geomcam, "Saturation", "Pixels with Saturated Hi Gain");
386 MHCamera disp17(geomcam, "ConversionMeans", "Conversion HiGain.vs.LoGain Means");
387 MHCamera disp18(geomcam, "ConversionSigmas", "Conversion HiGain.vs.LoGain Sigmas");
388 MHCamera disp19(geomcam, "HiGainPickup", "Number Pickup Events Hi Gain");
389 MHCamera disp20(geomcam, "LoGainPickup", "Number Pickup Events Lo Gain");
390 MHCamera disp21(geomcam, "HiGainBlackout", "Number Blackout Events Hi Gain");
391 MHCamera disp22(geomcam, "LoGainBlackout", "Number Blackout Events Lo Gain");
392 MHCamera disp23(geomcam, "Excluded", "Pixels Previously Excluded");
393 MHCamera disp24(geomcam, "UnSuitable", "Pixels NOT Suited for Further Analysis");
394 MHCamera disp25(geomcam, "UnReliable", "Pixels Suitable, but NOT Reliable for Further Analysis");
395 MHCamera disp26(geomcam, "HiGainOscillating", "Oscillating Pixels High Gain");
396 MHCamera disp27(geomcam, "LoGainOscillating", "Oscillating Pixels Low Gain");
397 MHCamera disp28(geomcam, "AbsTimeMean", "Abs. Arrival Times");
398 MHCamera disp29(geomcam, "AbsTimeRms", "RMS of Arrival Times");
399 MHCamera disp30(geomcam, "MeanTime", "Mean Rel. Arrival Times");
400 MHCamera disp31(geomcam, "SigmaTime", "Sigma Rel. Arrival Times");
401 MHCamera disp32(geomcam, "TimeProb", "Probability of Time Fit");
402 MHCamera disp33(geomcam, "TimeNotFitValid", "Pixels with not valid Fit Results");
403 MHCamera disp34(geomcam, "TimeOscillating", "Oscillating Pixels");
404 MHCamera disp35(geomcam, "TotalConv", "Conversion Factor to Photons");
405 MHCamera disp36(geomcam, "RMSperMean", "Charge histogram RMS per Mean");
406 MHCamera disp37(geomcam, "TotalConvPhe", "Conversion Factor to equiv. Phe's");
407
408 // Fitted charge means and sigmas
409 disp1.SetCamContent(fCalibrationCam, 0);
410 disp1.SetCamError( fCalibrationCam, 1);
411 disp2.SetCamContent(fCalibrationCam, 2);
412 disp2.SetCamError( fCalibrationCam, 3);
413
414 // Reduced Sigmas and reduced sigmas per charge
415 disp3.SetCamContent(fCalibrationCam, 5);
416 disp3.SetCamError( fCalibrationCam, 6);
417 disp4.SetCamContent(fCalibrationCam, 7);
418 disp4.SetCamError( fCalibrationCam, 8);
419
420 // F-Factor Method
421 disp5.SetCamContent(fCalibrationCam, 9);
422 disp5.SetCamError( fCalibrationCam, 10);
423 disp6.SetCamContent(fCalibrationCam, 11);
424 disp6.SetCamError( fCalibrationCam, 12);
425 disp7.SetCamContent(fCalibrationCam, 13);
426 disp7.SetCamError( fCalibrationCam, 14);
427
428 // Quantum Efficiencies
429 disp8.SetCamContent (fQECam, 0 );
430 disp8.SetCamError (fQECam, 1 );
431 disp9.SetCamContent (fQECam, 2 );
432 disp9.SetCamError (fQECam, 3 );
433 disp10.SetCamContent(fQECam, 4 );
434 disp10.SetCamError (fQECam, 5 );
435 disp11.SetCamContent(fQECam, 6 );
436 disp11.SetCamError (fQECam, 7 );
437
438 // Valid flags
439 disp12.SetCamContent(fQECam, 8 );
440 disp13.SetCamContent(fQECam, 9 );
441 disp14.SetCamContent(fQECam, 10);
442 disp15.SetCamContent(fQECam, 11);
443
444 // Conversion Hi-Lo
445 disp16.SetCamContent(fCalibrationCam, 25);
446 disp17.SetCamContent(fCalibrationCam, 16);
447 disp17.SetCamError (fCalibrationCam, 17);
448 disp18.SetCamContent(fCalibrationCam, 18);
449 disp18.SetCamError (fCalibrationCam, 19);
450
451 // Pickup and Blackout
452 disp19.SetCamContent(fCalibrationCam, 21);
453 disp20.SetCamContent(fCalibrationCam, 22);
454 disp21.SetCamContent(fCalibrationCam, 23);
455 disp22.SetCamContent(fCalibrationCam, 24);
456
457 // Pixels with defects
458 disp23.SetCamContent(fCalibrationCam, 20);
459 disp24.SetCamContent(fBadPixels, 6);
460 disp25.SetCamContent(fBadPixels, 7);
461
462 // Oscillations
463 disp26.SetCamContent(fBadPixels, 10);
464 disp27.SetCamContent(fBadPixels, 11);
465
466 // Arrival Times
467 disp28.SetCamContent(fCalibrationCam, 26);
468 disp28.SetCamError( fCalibrationCam, 27);
469 disp29.SetCamContent(fCalibrationCam, 27);
470
471 // RMS and Mean
472 disp36.SetCamContent(fCalibrationCam,32);
473 disp36.SetCamError(fCalibrationCam,33);
474
475 disp1.SetYTitle("Q [FADC cnts]");
476 disp2.SetYTitle("\\sigma_{Q} [FADC cnts]");
477
478 disp3.SetYTitle("\\sqrt{\\sigma^{2}_{Q} - RMS^{2}_{Ped}} [FADC cnts]");
479 disp4.SetYTitle("Red.Sigma/<Q> [1]");
480
481 disp5.SetYTitle("Photo-electons [1]");
482 disp6.SetYTitle("Phes/<Q> [FADC cnts^{-1}]");
483 disp7.SetYTitle("Total F-Factor [1]");
484
485 disp8.SetYTitle("QE [1]");
486 disp9.SetYTitle("QE [1]");
487 disp10.SetYTitle("QE [1]");
488 disp11.SetYTitle("QE [1]");
489
490 disp12.SetYTitle("[1]");
491 disp13.SetYTitle("[1]");
492 disp14.SetYTitle("[1]");
493 disp15.SetYTitle("[1]");
494 disp16.SetYTitle("[1]");
495
496 disp17.SetYTitle("<Q>(High)/<Q>(Low) [1]");
497 disp18.SetYTitle("\\sigma_{Q}(High)/\\sigma_{Q}(Low) [1]");
498
499 disp19.SetYTitle("[1]");
500 disp20.SetYTitle("[1]");
501 disp21.SetYTitle("[1]");
502 disp22.SetYTitle("[1]");
503 // disp23.SetYTitle("[1]");
504 // disp24.SetYTitle("[1]");
505 // disp25.SetYTitle("[1]");
506 disp26.SetYTitle("[1]");
507 disp27.SetYTitle("[1]");
508
509 disp28.SetYTitle("Mean Abs. Time [FADC sl.]");
510 disp29.SetYTitle("RMS Abs. Time [FADC sl.]");
511 disp35.SetYTitle("Conv.Factor [Ph/FADC cnts]");
512 disp36.SetYTitle("Charge RMS/<Q> [1]");
513 disp37.SetYTitle("Conv.Factor [Phe/FADC cnts]");
514
515 for (UInt_t i=0;i<geomcam.GetNumPixels();i++)
516 {
517
518 MCalibrationChargePix &pix = (MCalibrationChargePix&)fCalibrationCam[i];
519 MCalibrationQEPix &qe = (MCalibrationQEPix&)fQECam[i];
520
521 if (!pix.IsFFactorMethodValid())
522 continue;
523
524 const Float_t convphe = pix.GetMeanConvFADC2Phe();
525 const Float_t quaeff = qe.GetQECascadesFFactor();
526
527 disp35.Fill(i,convphe/quaeff);
528 disp35.SetUsed(i);
529
530 disp37.Fill(i,convphe/quaeff*MCalibrationQEPix::gkDefaultAverageQE);
531 disp37.SetUsed(i);
532 }
533
534
535 if (IsRelTimes())
536 {
537 disp30.SetCamContent(fRelTimeCam, 0);
538 disp30.SetCamError( fRelTimeCam, 1);
539 disp31.SetCamContent(fRelTimeCam, 2);
540 disp31.SetCamError( fRelTimeCam, 3);
541 disp32.SetCamContent(fRelTimeCam, 4);
542 disp33.SetCamContent(fBadPixels, 20);
543 disp34.SetCamContent(fBadPixels, 21);
544
545 disp30.SetYTitle("Time Offset [FADC units]");
546 disp31.SetYTitle("Timing resolution [FADC units]");
547 disp32.SetYTitle("P_{Time} [1]");
548 disp33.SetYTitle("[1]");
549 disp34.SetYTitle("[1]");
550 }
551
552 if (fDisplayType == kDataCheckDisplay)
553 {
554 // -------------------- FitCharge -------------------
555
556 TCanvas &c1 = fDisplay->AddTab("FitCharge");
557 c1.Divide(3, 3);
558
559 disp1.CamDraw( c1, 1, 3, 6); // MEAN CHARGES
560 disp36.CamDraw(c1, 2, 3, 6); // RMS per Charge
561 disp5.CamDraw( c1, 3, 3, 6); // PHOTO ELECTRONS
562
563 // -------------------- Conversion -------------------
564
565 TCanvas &c2 = fDisplay->AddTab("Conversion");
566 c2.Divide(3,3);
567
568 disp6.SetMinMax(fConvFADC2PheMin, fConvFADC2PheMax);
569 disp8.SetMinMax(fQEMin, fQEMax);
570 disp37.SetMinMax(fConvFADC2PheMin, fConvFADC2PheMax);
571
572 disp6.CamDraw( c2, 1, 3, 6); // CONVERSION FACTORS
573 disp8.CamDraw( c2, 2, 3, 6); // QUANTUM EFFICIENCY
574 disp37.CamDraw(c2, 3, 3, 6); // CONVERSION FADC TO PHOTONS
575
576 c2.cd(1);
577 DisplayReferenceLines(disp6, 2);
578
579 c2.cd(2);
580 DisplayReferenceLines(disp8, 0);
581
582 c2.cd(3);
583 DisplayReferenceLines(disp37, 1);
584
585 // -------------------- AbsTimes -------------------
586
587 TCanvas &c3 = fDisplay->AddTab("AbsTimes");
588 c3.Divide(2,3);
589
590 disp28.SetMinMax(fArrivalTimeMin, fArrivalTimeMax);
591
592 disp28.CamDraw(c3, 1, 2, 6); // Arrival times
593 disp29.CamDraw(c3, 2, 2, 6); // Arrival times rms
594
595 c3.cd(1);
596 DisplayReferenceLines(disp28, 3);
597
598 c3.cd(2);
599 DisplayReferenceLines(disp29, 4);
600
601 if (IsRelTimes())
602 {
603 // -------------------- RelTimes -------------------
604
605 TCanvas &c5 = fDisplay->AddTab("RelTimes");
606 c5.Divide(2,3);
607
608 disp30.SetMinMax(fTimeOffsetMin, fTimeOffsetMax);
609 disp31.SetMinMax(fTimeResolutionMin, fTimeResolutionMax);
610
611 disp30.CamDraw(c5, 1, 2, 6); // MEAN REL. ARR. TIMES
612 disp31.CamDraw(c5, 2, 2, 6); // JITTER Rel. Arr. Times
613
614 c5.cd(1);
615 DisplayReferenceLines(disp30, 5);
616
617 c5.cd(2);
618 DisplayReferenceLines(disp31, 6);
619 }
620
621 // -------------------- Unsuitable -------------------
622
623 //
624 // UNSUITABLE PIXELS
625 //
626 TCanvas &c4 = fDisplay->AddTab("Defect");
627 c4.Divide(2,2, 0.005, 0.005);
628
629 c4.cd(1);
630 gPad->SetBorderMode(0);
631 gPad->SetTicks();
632
633 MHCamera *obj8 = DrawBadPixelPad(disp24, kTRUE);
634
635 c4.cd(3);
636 gPad->SetBorderMode(0);
637 obj8->SetPrettyPalette();
638 obj8->Draw();
639
640 //
641 // UNRELIABLE PIXELS
642 //
643 c4.cd(2);
644 gPad->SetBorderMode(0);
645 gPad->SetTicks();
646
647 MHCamera *obj9 = DrawBadPixelPad(disp25, kFALSE);
648
649 c4.cd(4);
650 gPad->SetBorderMode(0);
651 obj9->SetPrettyPalette();
652 obj9->Draw();
653 return;
654 }
655
656 if (fDisplayType == kNormalDisplay)
657 {
658
659 // Charges
660 TCanvas &c11 = fDisplay->AddTab("FitCharge");
661 c11.Divide(2, 4);
662
663 disp1.CamDraw(c11, 1, 2, 5, 1);
664 disp2.CamDraw(c11, 2, 2, 5, 1);
665
666 // Reduced Sigmas
667 TCanvas &c12 = fDisplay->AddTab("RedSigma");
668 c12.Divide(2,4);
669
670 disp3.CamDraw(c12, 1, 2, 5, 1);
671 disp4.CamDraw(c12, 2, 2, 5, 1);
672
673 // F-Factor
674 TCanvas &c13 = fDisplay->AddTab("Phe's");
675 c13.Divide(3,4);
676
677 disp5.CamDraw(c13, 1, 3, 5, 1);
678 disp6.CamDraw(c13, 2, 3, 5, 1);
679 disp7.CamDraw(c13, 3, 3, 5, 1);
680
681 // QE's
682 TCanvas &c14 = fDisplay->AddTab("QE's");
683 c14.Divide(4,4);
684
685 disp8.CamDraw(c14, 1, 4, 5, 1);
686 disp9.CamDraw(c14, 2, 4, 5, 1);
687 disp10.CamDraw(c14, 3, 4, 5, 1);
688 disp11.CamDraw(c14, 4, 4, 5, 1);
689
690 // Defects
691 TCanvas &c15 = fDisplay->AddTab("Defect");
692 // c15.Divide(5,2);
693 c15.Divide(4,2);
694
695 /*
696 disp23.CamDraw(c15, 1, 5, 0);
697 disp24.CamDraw(c15, 2, 5, 0);
698 disp25.CamDraw(c15, 3, 5, 0);
699 disp26.CamDraw(c15, 4, 5, 0);
700 disp27.CamDraw(c15, 5, 5, 0);
701 */
702 disp24.CamDraw(c15, 1, 4, 0);
703 disp25.CamDraw(c15, 2, 4, 0);
704 disp26.CamDraw(c15, 3, 4, 0);
705 disp27.CamDraw(c15, 4, 4, 0);
706
707 // Abs. Times
708 TCanvas &c16 = fDisplay->AddTab("AbsTimes");
709 c16.Divide(2,3);
710
711 disp28.CamDraw(c16, 1, 2, 5);
712 disp29.CamDraw(c16, 2, 2, 5);
713
714 if (IsRelTimes())
715 {
716 // Rel. Times
717 TCanvas &c17 = fDisplay->AddTab("RelTimes");
718 c17.Divide(2,4);
719
720 disp30.CamDraw(c17, 1, 2, 5, 1);
721 disp31.CamDraw(c17, 2, 2, 5, 1);
722 }
723
724 return;
725 }
726
727 if (fDisplayType == kFullDisplay)
728 {
729 MHCalibrationCam *cam = (MHCalibrationCam*)plist.FindObject("MHCalibrationChargeCam");
730
731 for (Int_t sector=1;sector<cam->GetAverageSectors();sector++)
732 {
733 cam->GetAverageHiGainSector(sector).DrawClone("all");
734 cam->GetAverageLoGainSector(sector).DrawClone("all");
735 }
736
737 // Charges
738 TCanvas &c21 = fDisplay->AddTab("FitCharge");
739 c21.Divide(2, 4);
740
741 disp1.CamDraw(c21, 1, 2, 2, 1);
742 disp2.CamDraw(c21, 2, 2, 2, 1);
743
744 // Reduced Sigmas
745 TCanvas &c23 = fDisplay->AddTab("RedSigma");
746 c23.Divide(2,4);
747
748 disp3.CamDraw(c23, 1, 2, 2, 1);
749 disp4.CamDraw(c23, 2, 2, 2, 1);
750
751 // F-Factor
752 TCanvas &c24 = fDisplay->AddTab("Phe's");
753 c24.Divide(3,5);
754
755 disp5.CamDraw(c24, 1, 3, 2, 1, 1);
756 disp6.CamDraw(c24, 2, 3, 2, 1, 1);
757 disp7.CamDraw(c24, 3, 3, 2, 1, 1);
758
759 // QE's
760 TCanvas &c25 = fDisplay->AddTab("QE's");
761 c25.Divide(4,5);
762
763 disp8.CamDraw(c25, 1, 4, 2, 1, 1);
764 disp9.CamDraw(c25, 2, 4, 2, 1, 1);
765 disp10.CamDraw(c25, 3, 4, 2, 1, 1);
766 disp11.CamDraw(c25, 4, 4, 2, 1, 1);
767
768 // Validity
769 TCanvas &c26 = fDisplay->AddTab("Valid");
770 c26.Divide(4,2);
771
772 disp12.CamDraw(c26, 1, 4, 0);
773 disp13.CamDraw(c26, 2, 4, 0);
774 disp14.CamDraw(c26, 3, 4, 0);
775 disp15.CamDraw(c26, 4, 4, 0);
776
777 // Other info
778 TCanvas &c27 = fDisplay->AddTab("HiLoGain");
779 c27.Divide(3,3);
780
781 disp16.CamDraw(c27, 1, 3, 0);
782 disp17.CamDraw(c27, 2, 3, 1);
783 disp18.CamDraw(c27, 3, 3, 1);
784
785 // Pickup
786 TCanvas &c28 = fDisplay->AddTab("Pickup");
787 c28.Divide(4,2);
788
789 disp19.CamDraw(c28, 1, 4, 0);
790 disp20.CamDraw(c28, 2, 4, 0);
791 disp21.CamDraw(c28, 3, 4, 0);
792 disp22.CamDraw(c28, 4, 4, 0);
793
794 // Defects
795 TCanvas &c29 = fDisplay->AddTab("Defect");
796 // c29.Divide(5,2);
797 c29.Divide(4,2);
798
799 disp24.CamDraw(c29, 1, 4, 0);
800 disp25.CamDraw(c29, 2, 4, 0);
801 disp26.CamDraw(c29, 3, 4, 0);
802 disp27.CamDraw(c29, 4, 4, 0);
803
804 // Abs. Times
805 TCanvas &c30 = fDisplay->AddTab("AbsTimes");
806 c30.Divide(2,3);
807
808 disp28.CamDraw(c30, 1, 2, 2);
809 disp29.CamDraw(c30, 2, 2, 1);
810
811 if (IsRelTimes())
812 {
813 // Rel. Times
814 TCanvas &c31 = fDisplay->AddTab("RelTimes");
815 c31.Divide(3,5);
816
817 disp30.CamDraw(c31, 1, 3, 2, 1, 1);
818 disp31.CamDraw(c31, 2, 3, 2, 1, 1);
819 disp32.CamDraw(c31, 3, 3, 4, 1, 1);
820
821 // Time Defects
822 TCanvas &c32 = fDisplay->AddTab("DefTime");
823 c32.Divide(2,2);
824
825 disp33.CamDraw(c32, 1, 2, 0);
826 disp34.CamDraw(c32, 2, 2, 0);
827
828 MHCalibrationCam *ccam = (MHCalibrationCam*)plist.FindObject("MHCalibrationRelTimeCam");
829
830 for (Int_t sector=1;sector<ccam->GetAverageSectors();sector++)
831 {
832 ccam->GetAverageHiGainSector(sector).DrawClone("fourierevents");
833 ccam->GetAverageLoGainSector(sector).DrawClone("fourierevents");
834 }
835
836 }
837
838 return;
839 }
840}
841
842void MJCalibration::DisplayReferenceLines(const MHCamera &hist, const Int_t what) const
843{
844 MHCamera *cam = dynamic_cast<MHCamera*>(gPad->FindObject(hist.GetName()));
845 if (!cam)
846 return;
847
848 const MGeomCam *geom = cam->GetGeometry();
849
850 const Double_t x = geom->InheritsFrom("MGeomCamMagic") ? 397 : cam->GetNbinsX() ;
851
852 TLine line;
853 line.SetLineStyle(kDashed);
854 line.SetLineWidth(3);
855 line.SetLineColor(kBlue);
856
857 TLine *l1 = NULL;
858
859 switch (what)
860 {
861 case 0:
862 l1 = line.DrawLine(0, fRefQEInner, x, fRefQEInner);
863 break;
864 case 1:
865 l1 = line.DrawLine(0, fRefConvFADC2PheInner, x, fRefConvFADC2PheInner);
866 break;
867 case 2:
868 l1 = line.DrawLine(0, fRefFADC2PheInner, x, fRefFADC2PheInner );
869 break;
870 case 3:
871 l1 = line.DrawLine(0, fRefArrivalTimeInner, x, fRefArrivalTimeInner );
872 break;
873 case 4:
874 l1 = line.DrawLine(0, fRefArrivalTimeRmsInner, x, fRefArrivalTimeRmsInner );
875 break;
876 case 5:
877 l1 = line.DrawLine(0, 0, x, 0);
878 break;
879 case 6:
880 l1 = line.DrawLine(0, fRefTimeResolutionInner, x, fRefTimeResolutionInner );
881 break;
882 default:
883 break;
884 }
885
886 if (geom->InheritsFrom("MGeomCamMagic"))
887 {
888 const Double_t x2 = cam->GetNbinsX();
889
890 switch (what)
891 {
892 case 0:
893 line.DrawLine(x2, fRefQEOuter, 398, fRefQEOuter);
894 break;
895 case 1:
896 line.DrawLine(x2, fRefConvFADC2PheOuter, 398, fRefConvFADC2PheOuter );
897 break;
898 case 2:
899 line.DrawLine(x2, fRefFADC2PheOuter, 398, fRefFADC2PheOuter);
900 break;
901 case 3:
902 line.DrawLine(x2, fRefArrivalTimeOuter, 398, fRefArrivalTimeOuter);
903 break;
904 case 4:
905 line.DrawLine(x2, fRefArrivalTimeRmsOuter, 398, fRefArrivalTimeRmsOuter);
906 break;
907 case 5:
908 line.DrawLine(x2, fRefTimeOffsetOuter, 398, fRefTimeOffsetOuter);
909 break;
910 case 6:
911 line.DrawLine(x2, fRefTimeResolutionOuter, 398, fRefTimeResolutionOuter);
912 break;
913 default:
914 break;
915 }
916 }
917
918 TLegend *leg = new TLegend(0.6,0.85,0.9 ,0.95);
919 leg->SetBit(kCanDelete);
920 leg->AddEntry(l1, "Reference","l");
921 leg->Draw();
922}
923
924/*
925void MJCalibration::DisplayOutliers(TH1D *hist, const char* whatsmall, const char* whatbig) const
926{
927
928 const Int_t kNotDraw = 1<<9;
929 TF1 *f = hist->GetFunction("gaus");
930 f->ResetBit(kNotDraw);
931
932 const Float_t mean = f->GetParameter(1);
933 const Float_t lolim = mean - 4.0*f->GetParameter(2);
934 const Float_t uplim = mean + 4.0*f->GetParameter(2);
935 const Stat_t dead = hist->Integral(0,hist->FindBin(lolim)-1);
936 const Stat_t noisy = hist->Integral(hist->FindBin(uplim)+1,hist->GetNbinsX()+1);
937
938 const Double_t max = hist->GetBinContent(hist->GetMaximumBin());
939
940 const Double_t minl = hist->GetBinCenter(hist->GetXaxis()->GetFirst());
941 const Double_t maxl = hist->GetBinCenter(hist->GetXaxis()->GetLast());
942
943 TLatex deadtex;
944 deadtex.SetTextSize(0.07);
945 deadtex.DrawLatex(minl+0.015*(maxl-minl),max/1.1,
946 Form("%3i %s pixels",(Int_t)dead,whatsmall));
947
948 TLatex noisytex;
949 noisytex.SetTextSize(0.07);
950 noisytex.DrawLatex(minl+0.015*(maxl-minl),max/1.2,
951 Form("%3i %s pixels",(Int_t)noisy,whatbig));
952
953}
954*/
955
956void MJCalibration::FixDataCheckHist(TH1D &h, Double_t min, Double_t max)
957{
958 h.SetDirectory(NULL);
959 h.SetStats(kFALSE);
960 h.SetMinimum(min);
961 h.SetMaximum(max);
962
963 //
964 // set the labels bigger
965 //
966 TAxis *xaxe = h.GetXaxis();
967 TAxis *yaxe = h.GetYaxis();
968
969 xaxe->CenterTitle();
970 yaxe->CenterTitle();
971 xaxe->SetTitleSize(0.06);
972 yaxe->SetTitleSize(0.06);
973 xaxe->SetTitleOffset(0.8);
974 yaxe->SetTitleOffset(0.85);
975 xaxe->SetLabelSize(0.05);
976 yaxe->SetLabelSize(0.05);
977}
978
979// --------------------------------------------------------------------------
980//
981// Retrieve the output file written by WriteResult()
982//
983const char* MJCalibration::GetOutputFileName() const
984{
985 return Form("calib%08d.root", fSequence.GetSequence());
986}
987
988// --------------------------------------------------------------------------
989//
990// Read the following values from resource file:
991//
992// ConvFADC2PheMin
993// ConvFADC2PheMax
994// ConvFADC2PhotMin
995// ConvFADC2PhotMax
996//
997// QEMin
998// QEMax
999//
1000// ArrivalTimeMin
1001// ArrivalTimeMax
1002//
1003// TimeOffsetMin
1004// TimeOffsetMax
1005// TimeResolutionMin
1006// TimeResolutionMax
1007//
1008// RefConvFADC2PheInner
1009// RefConvFADC2PheOuter
1010// RefConvFADC2PhotInner
1011// RefConvFADC2PhotOuter
1012//
1013// RefQEInner
1014// RefQEOuter
1015//
1016// RefArrivalTimeInner
1017// RefArrivalTimeOuter
1018// RefArrivalTimeRmsInner
1019// RefArrivalTimeRmsOuter
1020//
1021// RefTimeOffsetOuter
1022// RefTimeResolutionInner
1023// RefTimeResolutionOuter
1024//
1025void MJCalibration::ReadReferenceFile()
1026{
1027 TEnv refenv(fReferenceFile);
1028
1029 fConvFADC2PheMin = refenv.GetValue("ConvFADC2PheMin",fConvFADC2PheMin);
1030 fConvFADC2PheMax = refenv.GetValue("ConvFADC2PheMax",fConvFADC2PheMax);
1031 fConvFADC2PhotMin = refenv.GetValue("ConvFADC2PhotMin",fConvFADC2PhotMin);
1032 fConvFADC2PhotMax = refenv.GetValue("ConvFADC2PhotMax",fConvFADC2PhotMax);
1033 fQEMin = refenv.GetValue("QEMin",fQEMin);
1034 fQEMax = refenv.GetValue("QEMax",fQEMax);
1035 fArrivalTimeMin = refenv.GetValue("ArrivalTimeMin",fArrivalTimeMin);
1036 fArrivalTimeMax = refenv.GetValue("ArrivalTimeMax",fArrivalTimeMax);
1037 fTimeOffsetMin = refenv.GetValue("TimeOffsetMin",fTimeOffsetMin);
1038 fTimeOffsetMax = refenv.GetValue("TimeOffsetMax",fTimeOffsetMax);
1039 fTimeResolutionMin = refenv.GetValue("TimeResolutionMin",fTimeResolutionMin);
1040 fTimeResolutionMax = refenv.GetValue("TimeResolutionMax",fTimeResolutionMax);
1041
1042 fRefFADC2PheInner = refenv.GetValue("RefFADC2PheInner",fRefFADC2PheInner);
1043 fRefFADC2PheOuter = refenv.GetValue("RefFADC2PheOuter",fRefFADC2PheOuter);
1044 fRefConvFADC2PhotInner = refenv.GetValue("RefConvFADC2PhotInner",fRefConvFADC2PhotInner);
1045 fRefConvFADC2PhotOuter = refenv.GetValue("RefConvFADC2PhotOuter",fRefConvFADC2PhotOuter);
1046 fRefConvFADC2PheInner = refenv.GetValue("RefConvFADC2PheInner",fRefConvFADC2PheInner);
1047 fRefConvFADC2PheOuter = refenv.GetValue("RefConvFADC2PheOuter",fRefConvFADC2PheOuter);
1048 fRefQEInner = refenv.GetValue("RefQEInner",fRefQEInner);
1049 fRefQEOuter = refenv.GetValue("RefQEOuter",fRefQEOuter);
1050 fRefArrivalTimeInner = refenv.GetValue("RefArrivalTimeInner",fRefArrivalTimeInner);
1051 fRefArrivalTimeOuter = refenv.GetValue("RefArrivalTimeOuter",fRefArrivalTimeOuter);
1052 fRefArrivalTimeRmsInner = refenv.GetValue("RefArrivalTimeRmsInner",fRefArrivalTimeRmsInner);
1053 fRefArrivalTimeRmsOuter = refenv.GetValue("RefArrivalTimeRmsOuter",fRefArrivalTimeRmsOuter);
1054 fRefTimeOffsetOuter = refenv.GetValue("RefTimeOffsetOuter",fRefTimeOffsetOuter);
1055 fRefTimeResolutionInner = refenv.GetValue("RefTimeResolutionInner",fRefTimeResolutionInner);
1056 fRefTimeResolutionOuter = refenv.GetValue("RefTimeResolutionOuter",fRefTimeResolutionOuter);
1057}
1058
1059// --------------------------------------------------------------------------
1060//
1061// Read the following values from resource file.
1062//
1063Bool_t MJCalibration::ReadHiLoCalibFile()
1064{
1065 if (fExtractor && !fExtractor->HasLoGain())
1066 return kTRUE;
1067
1068// if (!fIsHiLoCalibration)
1069// return kTRUE;
1070
1071 // We use the night time stamp to determine the period
1072 // because the night must be in the sequence file
1073 const MTime &night = fSequence.GetNight();
1074 const Int_t period = night.GetMagicPeriod();
1075
1076 // Open resource file
1077 MEnv env(fHiLoCalibFile);
1078 if (!env.IsValid())
1079 {
1080 *fLog << err << "ERROR - Resource file " << fHiLoCalibFile;
1081 *fLog << " could not be opened... abort." << endl;
1082 return kFALSE;
1083 }
1084
1085 // Defined resource id
1086 const TString id = fSequence.IsMonteCarlo() ? "MC" : Form("%02d", period);
1087
1088 // Check for a valid entry for the correct period
1089 TString fname = env.GetValue(id, "");
1090 if (fname.IsNull())
1091 {
1092 *fLog << err << "ERROR - No entry for resource id '" << id;
1093 *fLog << "' found in " << fHiLoCalibFile << "... looking for default." << endl;
1094 return kFALSE;
1095/*
1096 *fLog << warn << "WARNING - No entry for period " << period;
1097 *fLog << " found in " << fHiLoCalibFile << "... looking for default." << endl;
1098
1099 fname = env.GetValue("00", "");
1100 if (fname.IsNull())
1101 {
1102 *fLog << err << "ERROR - No default entry (00) found in ";
1103 *fLog << fHiLoCalibFile << "... abort." << endl;
1104 return kFALSE;
1105 }*/
1106 }
1107
1108 *fLog << inf << "Reading Hi-/Lo-Gain calibration constants from " << fname << endl;
1109
1110 // Open file with calibration constants
1111 TFile file(fname, "READ");
1112 if (!file.IsOpen())
1113 {
1114 *fLog << err << "ERROR - Couldn't open file " << fname << " for reading... abort." << endl;
1115 return kFALSE;
1116 }
1117
1118 // read calibration constants
1119 MHCamEvent hilocam;
1120 if (hilocam.Read()<=0)
1121 {
1122 *fLog << err << "ERROR - Unable to read MHCamEvent from " << fname << "... abort." << endl;
1123 return kFALSE;
1124 }
1125
1126 // Get histogram with constants
1127 MHCamera *hist = hilocam.GetHist();
1128 if (!hist)
1129 {
1130 *fLog << err << "ERROR - MHCamEvent from " << fname << " empty... abort." << endl;
1131 return kFALSE;
1132 }
1133
1134 // Do some sanity stuff
1135 if (fCalibrationCam.GetSize() < 1)
1136 fCalibrationCam.InitSize(hist->GetNumPixels());
1137
1138 if (fBadPixels.GetSize() < 1)
1139 fBadPixels.InitSize(hist->GetNumPixels());
1140
1141 if ((UInt_t)fCalibrationCam.GetSize() != hist->GetNumPixels())
1142 {
1143 *fLog << err << "ERROR - Size mismatch MHCamEvent and MCalibrationChargeCam.. abort." << endl;
1144 return kFALSE;
1145 }
1146
1147 // Copy the constants to their final location
1148 // FIXME: For what the hell do we need to have the constants in
1149 // in MCalibrationChargeCam?
1150 for (UInt_t i=0; i<hist->GetNumPixels(); i++)
1151 {
1152 hist->SetBit(MHCamera::kProfile);
1153 Double_t v = hist->GetBinContent(i);
1154 hist->SetBit(MHCamera::kErrorMean);
1155 Double_t e = hist->GetBinError(i);
1156 hist->ResetBit(MHCamera::kErrorMean);
1157 Double_t s = hist->GetBinError(i);
1158
1159 if (!hist->IsUsed(i))
1160 {
1161 fBadPixels[i].SetUncalibrated(MBadPixelsPix::kConversionHiLoNotValid);
1162 v = e = s = -1;
1163 }
1164
1165 MCalibrationChargePix &cpix = (MCalibrationChargePix&)fCalibrationCam[i];
1166 cpix.SetConversionHiLo(v);
1167 cpix.SetConversionHiLoErr(e);
1168 cpix.SetConversionHiLoSigma(s);
1169 }
1170
1171 return kTRUE;
1172}
1173
1174// --------------------------------------------------------------------------
1175//
1176// MJCalibration allows to setup several option by a resource file:
1177// MJCalibration.Display: full, datacheck, normal
1178// MJCalibration.RelTimeCalibration: yes,no
1179// MJCalibration.DataCheck: yes,no
1180// MJCalibration.Debug: yes,no
1181// MJCalibration.UseBlindPixel: yes,no
1182// MJCalibration.UsePINDiode: yes,no
1183// MJCalibration.Geometry: MGeomCamMagic, MGeomCamECO1000
1184//
1185// Name of a file containing reference values (see ReadReferenceFile)
1186// Prefix.ReferenceFile: filename
1187// (see ReadReferenceFile)
1188//
1189// For more details see the class description and the corresponding Getters
1190//
1191Bool_t MJCalibration::CheckEnvLocal()
1192{
1193 TString dis = GetEnv("Display", "");
1194 if (dis.BeginsWith("Full", TString::kIgnoreCase))
1195 SetFullDisplay();
1196 if (dis.BeginsWith("DataCheck", TString::kIgnoreCase))
1197 SetDataCheckDisplay();
1198 if (dis.BeginsWith("Normal", TString::kIgnoreCase))
1199 SetNormalDisplay();
1200
1201 if (!MJCalib::CheckEnvLocal())
1202 return kFALSE;
1203
1204 SetRelTimeCalibration(GetEnv("RelTimeCalibration", IsRelTimes()));
1205 SetDebug(GetEnv("Debug", IsDebug()));
1206
1207 SetUseBlindPixel(GetEnv("UseBlindPixel", IsUseBlindPixel()));
1208 SetUsePINDiode(GetEnv("UsePINDiode", IsUsePINDiode()));
1209 SetGeometry(GetEnv("Geometry", fGeometry));
1210
1211 fMinEvents = (UInt_t)GetEnv("MinEvents", (Int_t)fMinEvents);
1212
1213 fReferenceFile = GetEnv("ReferenceFile",fReferenceFile.Data());
1214 ReadReferenceFile();
1215
1216 fHiLoCalibFile = GetEnv("HiLoCalibFile",fHiLoCalibFile.Data());
1217
1218 /*
1219 if (IsUseMC() && !fHiLoCalibFile.EndsWith("_mc.root"))
1220 {
1221 if (!fHiLoCalibFile.EndsWith(".root"))
1222 {
1223 *fLog << warn << "WARNING - Hi-/Lo-Gain intercalibration file ";
1224 *fLog << fHiLoCalibFile << " has not .root as extension..." << endl;
1225 }
1226 else
1227 fHiLoCalibFile.Insert(fHiLoCalibFile.Length()-5, "_mc");
1228 }
1229 */
1230
1231 return ReadHiLoCalibFile();
1232}
1233
1234void MJCalibration::InitBlindPixel(MExtractBlindPixel &blindext,
1235 MHCalibrationChargeBlindCam &blindcam)
1236{
1237 const Int_t run = fSequence.GetLastRun();
1238
1239 //
1240 // Initialize the blind pixel. Unfortunately, there is a hardware
1241 // difference in the first blind pixel until run
1242 // gkSecondBlindPixelInstallation and the later setup. The first
1243 // needs to use a filter because of the length of spurious NSB photon
1244 // signals. The latter get better along extracting the amplitude
1245 // from a small window.
1246 //
1247 const MCalibrationBlindCamOneOldStyle one;
1248 const MCalibrationBlindCamTwoNewStyle two;
1249 const MCalibrationBlindCamThreeNewStyle three;
1250
1251 const MCalibrationBlindCam &blindresults =
1252 run<gkSecondBlindPixelInstallation ? static_cast<MCalibrationBlindCam>(one) :
1253 (run<gkThirdBlindPixelInstallation ? static_cast<MCalibrationBlindCam>(two) : static_cast<MCalibrationBlindCam>(three));
1254
1255 blindresults.Copy(fCalibrationBlindCam);
1256
1257 blindext.SetExtractionType(MExtractBlindPixel::kIntegral);
1258 blindext.SetExtractionType(MExtractBlindPixel::kFilter);
1259
1260 if (run<gkSecondBlindPixelInstallation)
1261 {
1262 blindext.SetRange(10,19,0,6);
1263 blindext.SetNSBFilterLimit(70);
1264 }
1265 else
1266 {
1267 blindext.SetRange(5,8,0,2);
1268 blindext.SetNSBFilterLimit(38);
1269 }
1270
1271 if (run>=gkThirdBlindPixelInstallation)
1272 blindext.SetDataType(MExtractBlindPixel::kRawEvt2);
1273}
1274
1275// --------------------------------------------------------------------------
1276//
1277// Execute the task list and the eventloop:
1278//
1279// - Check the colour of the files in fRuns (FindColor()), otherwise return
1280// - Check for consistency between run numbers and number of files
1281// - Add fRuns to MReadMarsFile
1282// - Put into MParList:
1283// 1) MPedestalCam (pedcam)
1284// 2) MCalibrationQECam (fQECam)
1285// 3) MCalibrationChargeCam (fCalibrationCam)
1286// 4) MCalibrationRelTimeCam (fRelTimeCam) (only if flag IsRelTimes() is chosen)
1287// 5) MBadPixelsCam (fBadPixels)
1288// 6) MCalibrationChargePINDiode
1289// 7) MCalibrationBlindPix
1290// - Put into the MTaskList:
1291// 1) MReadMarsFile
1292// 2) MBadPixelsMerge
1293// 3) MGeomApply
1294// 4) MExtractor
1295// 5) MExtractPINDiode
1296// 6) MExtractBlindPixel
1297// 7) MExtractTime (only if flag IsRelTimes() is chosen)
1298// 8) MContinue(MFCosmics)
1299// 9) MFillH("MHCalibrationChargePINDiode", "MExtractedSignalPINDiode", "FillPINDiode")
1300// 10) MFillH("MHCalibrationChargeBlindCam", "MExtractedSignalBlindPixel", "FillBlindCam")
1301// 11) MFillH("MHCalibrationChargeCam", "MExtractedSignalCam", "FillChargeCam")
1302// 12) MFillH("MHCalibrationChargeCam", "MExtractedSignalCam", "FillRelTime")
1303// 13) MCalibrationChargeCalc
1304// 14) MFillH("MHCalibrationRelTimeCam", "MArrivalTimeCam") (only if flag IsRelTimes() is chosen)
1305// 15) MCalibrationRelTimeCalc
1306// - Execute MEvtLoop
1307// - DisplayResult()
1308// - WriteResult()
1309//
1310Bool_t MJCalibration::Process(MPedestalCam &pedcam)
1311{
1312 if (!fSequence.IsValid())
1313 {
1314 *fLog << err << "ERROR - Sequence invalid..." << endl;
1315 return kFALSE;
1316 }
1317
1318 *fLog << inf;
1319 fLog->Separator(GetDescriptor());
1320 *fLog << "Calculate calibration constants from Sequence #";
1321 *fLog << fSequence.GetSequence() << endl << endl;
1322
1323 // --------------------------------------------------------------------------------
1324
1325 if (!CheckEnv())
1326 return kFALSE;
1327
1328 // --------------------------------------------------------------------------------
1329
1330 // Setup Tasklist
1331 MParList plist;
1332 MTaskList tlist;
1333 plist.AddToList(&tlist);
1334 plist.AddToList(this); // take care of fDisplay!
1335
1336 MDirIter iter;
1337 if (fSequence.GetRuns(iter, MSequence::kRawCal)<=0)
1338 return kFALSE;
1339
1340 //
1341 // Input containers
1342 //
1343 pedcam.SetName("MPedestalCam"); // MPedestalFundamental
1344 plist.AddToList(&pedcam);
1345 plist.AddToList(&fBadPixels);
1346
1347 //
1348 // Calibration Results containers
1349 //
1350 plist.AddToList(&fQECam);
1351 plist.AddToList(&fCalibrationCam);
1352 plist.AddToList(&fRelTimeCam);
1353 if (IsUseBlindPixel())
1354 plist.AddToList(&fCalibrationBlindCam);
1355 if (IsUsePINDiode())
1356 plist.AddToList(&fCalibrationPINDiode);
1357
1358 //
1359 // Initialize two histogram containers which could be modified in this class
1360 //
1361 MHCalibrationRelTimeCam reltimecam;
1362 MHCalibrationChargeCam chargecam;
1363 MHCalibrationChargeBlindCam blindcam;
1364 plist.AddToList(&chargecam);
1365
1366 if (IsUseBlindPixel())
1367 plist.AddToList(&blindcam);
1368 if (IsRelTimes())
1369 plist.AddToList(&reltimecam);
1370 //
1371 // Data Reading tasks
1372 //
1373 MReadMarsFile read("Events");
1374 MRawFileRead rawread(NULL);
1375 rawread.SetForceMode(); // Ignore broken time-stamps
1376
1377 if (!fSequence.IsMonteCarlo())
1378 {
1379 rawread.AddFiles(iter);
1380 tlist.AddToList(&rawread);
1381 }
1382 else
1383 {
1384 read.DisableAutoScheme();
1385 read.AddFiles(iter);
1386 tlist.AddToList(&read);
1387 }
1388
1389 //
1390 // Other Tasks
1391 //
1392
1393 // Set the default for data version earlier than 5, where no valid
1394 // trigger pattern exists. There should not be pin diode or other
1395 // types of events inside the calibration files which should be skipped,
1396 // anyway. So we set the default such that the MContinue ccalib
1397 // will never be executed.
1398 // We allow to have only calibration events before Prescaling.
1399 // We require that the calibration events have not been prescaled (why?)
1400 MTriggerPatternDecode trgpat;
1401 MFTriggerPattern fcalib("CalibFilter");
1402 fcalib.SetDefault(kFALSE);
1403 fcalib.DenyCalibration(MFTriggerPattern::kPrescaled);
1404 fcalib.DenyAll();
1405 fcalib.AllowCalibration();
1406
1407 MContinue ccalib(&fcalib,"ContTrigPattern");
1408
1409 MCalibrationPatternDecode decode;
1410 MGeomApply apply;
1411 apply.SetGeometry(fGeometry);
1412
1413 MBadPixelsMerge merge(&fBadPixels);
1414 MExtractPINDiode pinext;
1415 MExtractBlindPixel blindext;
1416
1417 if (IsUseBlindPixel())
1418 InitBlindPixel(blindext, blindcam);
1419
1420// MExtractSlidingWindow extract2;
1421// MExtractTimeHighestIntegral timehigh;
1422 MCalibrationChargeCalc calcalc;
1423 MCalibrationRelTimeCalc timecalc;
1424 calcalc.SetExtractor(fExtractor);
1425
1426 if (IsDebug())
1427 {
1428 chargecam.SetDebug();
1429 calcalc.SetDebug();
1430 }
1431
1432 //
1433 // Calibration histogramming
1434 //
1435 MFillH fillpin("MHCalibrationChargePINDiode", "MExtractedSignalPINDiode", "FillPINDiode");
1436 MFillH fillbnd("MHCalibrationChargeBlindCam", "MExtractedSignalBlindPixel", "FillBlindCam");
1437 MFillH fillcam("MHCalibrationChargeCam", "MExtractedSignalCam", "FillChargeCam");
1438 MFillH filltme("MHCalibrationRelTimeCam", "MArrivalTimeCam", "FillRelTime");
1439 fillpin.SetBit(MFillH::kDoNotDisplay);
1440 fillbnd.SetBit(MFillH::kDoNotDisplay);
1441 fillcam.SetBit(MFillH::kDoNotDisplay);
1442 filltme.SetBit(MFillH::kDoNotDisplay);
1443
1444 //
1445 // Set default extractors in case, none has been set...
1446 //
1447 /*
1448 if (!fExtractor)
1449 fExtractor = &extract2;
1450 if (!fTimeExtractor)
1451 fTimeExtractor = &timehigh;
1452 */
1453 MExtractTimeAndChargeSpline spline;
1454 if (!fExtractor)
1455 fExtractor = &spline;
1456// if (!fTimeExtractor)
1457// fTimeExtractor = &timehigh;
1458
1459 const Bool_t istimecharge = fExtractor->InheritsFrom("MExtractTimeAndCharge");
1460
1461 //
1462 // Look if the extractor is a pure charge or also a time extractor
1463 //
1464 if (istimecharge)
1465 {
1466 if (fExtractorCam.GetSize() == pedcam.GetSize())
1467 calcalc.SetPedestals(&fExtractorCam);
1468 else
1469 {
1470 *fLog << err << GetDescriptor() << "ERROR - ";
1471 *fLog << "Used Extractor derives from MExtractTimeAndCharge, " << endl;
1472 *fLog << "but MExtractorCam size " << fExtractorCam.GetSize() << " ";
1473 *fLog << "mismatch pedcam size " << pedcam.GetSize() << "! " << endl;
1474 return kFALSE;
1475 }
1476 }
1477
1478 //
1479 // Setup more tasks and tasklist
1480 //
1481 MTaskEnv taskenv("ExtractSignal");
1482 taskenv.SetDefault(fExtractor);
1483
1484 tlist.AddToList(&trgpat);
1485 if (fColor != MCalibrationCam::kCT1)
1486 tlist.AddToList(&ccalib);
1487 tlist.AddToList(&decode);
1488 tlist.AddToList(&merge);
1489 tlist.AddToList(&apply);
1490
1491 // Produce pedestal subtracted raw-data
1492 MPedestalSubtract pedsub;
1493 pedsub.SetPedestalCam(&pedcam);
1494 tlist.AddToList(&pedsub);
1495
1496 MCalibColorSet colorset;
1497 if (fColor != MCalibrationCam::kNONE)
1498 colorset.SetExplicitColor(fColor);
1499 tlist.AddToList(&colorset);
1500
1501 tlist.AddToList(&taskenv);
1502
1503 if (IsUsePINDiode())
1504 tlist.AddToList(&pinext);
1505 if (IsUseBlindPixel())
1506 tlist.AddToList(&blindext);
1507
1508 MTaskEnv taskenv2("ExtractTime");
1509 if (!istimecharge)
1510 {
1511 taskenv2.SetDefault(fTimeExtractor);
1512
1513 if (IsRelTimes())
1514 tlist.AddToList(&taskenv2);
1515 }
1516
1517 //
1518 // Apply a filter against cosmics (this is for the old times in which
1519 // the calibration events where triggered by level 1 and for
1520 // sanity against empty trigger events)
1521 //
1522 MFCosmics cosmics;
1523 cosmics.SetMaxEmptyPixels(0.05);
1524 cosmics.SetMaxAcceptedFraction(0.5); // max=0.5 cosmics
1525
1526 MContinue cont(&cosmics, "ContCosmics");
1527 tlist.AddToList(&cont);
1528
1529 tlist.AddToList(&fillcam);
1530
1531 if (IsUseBlindPixel())
1532 tlist.AddToList(&fillbnd);
1533 if (IsUsePINDiode())
1534 tlist.AddToList(&fillpin);
1535
1536 tlist.AddToList(&calcalc);
1537
1538 if (IsRelTimes())
1539 {
1540 tlist.AddToList(&filltme);
1541 tlist.AddToList(&timecalc);
1542 }
1543
1544 MHCamEvent evt2(0, "Extra'd", "Extracted Calibration Signal;;S [cnts/sl]");
1545 MHCamEvent evt9(4, "ArrTm", "Extracted ArrivalTime;;T");
1546
1547 MFillH fill2(&evt2, "MExtractedSignalCam", "FillExtractedSignal");
1548 MFillH fill9(&evt9, "MArrivalTimeCam", "FillArrivalTime");
1549
1550 tlist.AddToList(&fill2);
1551 tlist.AddToList(&fill9);
1552
1553 /*
1554 MFillH fillP("MHPulseShape", "", "FillPulseShape");
1555 fillP.SetNameTab("Pulse");
1556 tlist.AddToList(&fillP);
1557 */
1558
1559 // Create and setup the eventloop
1560 MEvtLoop evtloop(fName);
1561 evtloop.SetParList(&plist);
1562 evtloop.SetDisplay(fDisplay);
1563 evtloop.SetLogStream(fLog);
1564 if (!SetupEnv(evtloop))
1565 return kFALSE;
1566
1567 if (!taskenv.GetTask() && !taskenv2.GetTask())
1568 {
1569 *fLog << err << "ERROR - Neither ExtractSignal nor ExtractTime initialized or both '<dummy>'." << endl;
1570 return kFALSE;
1571 }
1572
1573 if (!WriteTasks(taskenv.GetTask(), istimecharge ? 0 : taskenv2.GetTask()))
1574 return kFALSE;
1575
1576 // Execute first analysis
1577 const Bool_t rc = evtloop.Eventloop();
1578
1579 if (!fCalibrationPINDiode.IsValid())
1580 SetUsePINDiode(kFALSE);
1581
1582 // Only display result if PreProcessing was successfull
1583 const Int_t numexec = !fSequence.IsMonteCarlo() ? rawread.GetNumExecutions() : read.GetNumExecutions();
1584 if (numexec>0)
1585 {
1586 DisplayResult(plist);
1587 if (!WriteResult(plist))
1588 return kFALSE;
1589 }
1590
1591 if (!rc)
1592 {
1593 *fLog << err << GetDescriptor() << ": Failed." << endl;
1594 return kFALSE;
1595 }
1596
1597 if (calcalc.GetNumExecutions()<fMinEvents)
1598 {
1599 *fLog << err << GetDescriptor() << ": Failed. Less than the required " << fMinEvents << " evts processed." << endl;
1600 return kFALSE;
1601 }
1602
1603 // --------------------------------------------------------------------------------
1604
1605 if (fIsPixelCheck)
1606 {
1607 chargecam[fCheckedPixId].DrawClone("datacheck");
1608 chargecam(fCheckedPixId).DrawClone("datacheck");
1609
1610 if (IsRelTimes())
1611 {
1612 reltimecam[fCheckedPixId].DrawClone("");
1613 reltimecam(fCheckedPixId).DrawClone("");
1614 }
1615 }
1616
1617 *fLog << all << GetDescriptor() << ": Done." << endl << endl << endl;
1618
1619 return kTRUE;
1620}
1621
1622/*
1623Bool_t MJCalibration::WriteEventloop(MEvtLoop &evtloop) const
1624{
1625 if (IsNoStorage())
1626 return kTRUE;
1627
1628 if (fPathOut.IsNull())
1629 return kTRUE;
1630
1631 const TString oname(GetOutputFile());
1632
1633 *fLog << inf << "Writing to file: " << oname << endl;
1634
1635 TFile file(oname, fOverwrite?"RECREATE":"NEW", "File created by MJCalibration", 9);
1636 if (!file.IsOpen())
1637 {
1638 *fLog << err << "ERROR - Couldn't open file " << oname << " for writing..." << endl;
1639 return kFALSE;
1640 }
1641
1642 if (evtloop.Write(fName)<=0)
1643 {
1644 *fLog << err << "Unable to write MEvtloop to " << oname << endl;
1645 return kFALSE;
1646 }
1647
1648 return kTRUE;
1649}
1650*/
1651
1652Bool_t MJCalibration::WriteTasks(MTask *t1, MTask *t2) const
1653{
1654 if (IsNoStorage())
1655 return kTRUE;
1656
1657 TObjArray cont;
1658 if (t1)
1659 cont.Add(t1);
1660 if (t2)
1661 cont.Add(t2);
1662
1663 return WriteContainer(cont, GetOutputFileName(), fOverwrite?"RECREATE":"NEW");
1664}
1665
1666// --------------------------------------------------------------------------
1667//
1668// Write the result into the output file GetOutputFile(), if fOutputPath exists.
1669//
1670// The following containers are written:
1671// - MStatusDisplay
1672// - MCalibrationChargeCam
1673// - MCalibrationBlindCam
1674// - MCalibrationQECam
1675// - MCalibrationChargePINDiode
1676// - MBadPixelsCam
1677//
1678// If the flag kRelTimes is set, then also:
1679// - MCalibrationRelTimeCam
1680//
1681Bool_t MJCalibration::WriteResult(MParList &plist)
1682{
1683 if (IsNoStorage())
1684 return kTRUE;
1685
1686 TObjArray cont;
1687 cont.Add(&fBadPixels);
1688 cont.Add(&fCalibrationCam);
1689 cont.Add(&fQECam);
1690 cont.Add(&fCalibrationBlindCam);
1691 cont.Add(&fCalibrationPINDiode);
1692
1693 if (IsRelTimes())
1694 cont.Add(&fRelTimeCam);
1695
1696 if (fExtractorCam.GetSize() != 0)
1697 {
1698 fExtractorCam.SetName("MPedestalExtracted");
1699 cont.Add(&fExtractorCam);
1700 }
1701
1702 TObject *pedcam = plist.FindObject("MPedestalCam");
1703 if (!pedcam)
1704 *fLog << warn << " - WARNING - MPedestalCam (fundamental)... not found for writing!" << endl;
1705 else
1706 cont.Add(pedcam);
1707
1708 TObject *geom = plist.FindObject("MGeomCam");
1709 if (!geom)
1710 *fLog << warn << " - WARNING - MGeomCam... not found for writing!" << endl;
1711 else
1712 cont.Add(geom);
1713
1714 if (IsHistsStorage())
1715 {
1716 cont.Add(plist.FindObject("MHCalibrationChargeCam"));
1717 cont.Add(plist.FindObject("MHCalibrationChargeBlindCam"));
1718 cont.Add(plist.FindObject("MHCalibrationChargePINDiode"));
1719 if (IsRelTimes())
1720 cont.Add(plist.FindObject("MHCalibrationRelTimeCam"));
1721 }
1722
1723 if (fDisplay)
1724 cont.Add(fDisplay);
1725
1726 cont.Add(const_cast<TEnv*>(GetEnv()));
1727 cont.Add(&fSequence);
1728
1729 return WriteContainer(cont, GetOutputFileName(), "UPDATE");
1730}
1731
1732
1733void MJCalibration::DisplayDoubleProject(const MHCamera &cam)
1734{
1735 const UInt_t n = cam.GetGeometry()->GetNumAreas();
1736
1737 TVirtualPad *pad = gPad;
1738 pad->Divide(n, 1, 1e-5, 1e-5);;
1739
1740 for (UInt_t i=0; i<n; i++)
1741 {
1742 pad->cd(i+1);
1743 gPad->SetBorderMode(0);
1744 gPad->SetTicks();
1745
1746 TH1D &h = *cam.ProjectionS(TArrayI(), TArrayI(1, (Int_t*)&i), MString::Format("%s_%d", cam.GetName(), i));
1747 FixDataCheckHist(h);
1748 h.SetTitle(MString::Format("%s %d",cam.GetTitle(), i));
1749 h.SetDirectory(NULL);
1750 h.SetBit(kCanDelete);
1751 h.Draw();
1752
1753 h.Fit("gaus", "Q");
1754
1755 TF1 *f = h.GetFunction("gaus");
1756 if (f)
1757 {
1758 f->SetLineWidth(2);
1759 f->SetLineColor(kBlue);
1760 }
1761 }
1762}
Note: See TracBrowser for help on using the repository browser.