Context Navigation

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

Visit:

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

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

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

Download in other formats: