Context Navigation

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

Visit:

Last change on this file since 3762 was 3740, checked in by gaug, 21 years ago
* empty log message *
File size: 37.2 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: