Context Navigation

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

Visit:

Last change on this file since 3392 was 3387, checked in by tbretz, 21 years ago
* empty log message *
File size: 16.4 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	!
20	! Copyright: MAGIC Software Development, 2000-2004
21	!
22	!
23	\* ======================================================================== */
24
25	/////////////////////////////////////////////////////////////////////////////
26	//
27	// MJCalibration
28	//
29	/////////////////////////////////////////////////////////////////////////////
30	#include "MJCalibration.h"
31
32	#include <TF1.h>
33	#include <TFile.h>
34	#include <TStyle.h>
35	#include <TCanvas.h>
36	#include <TSystem.h>
37
38	#include "MLog.h"
39	#include "MLogManip.h"
40
41	#include "MRunIter.h"
42	#include "MParList.h"
43	#include "MTaskList.h"
44	#include "MEvtLoop.h"
45
46	#include "MHCamera.h"
47
48	#include "MPedestalCam.h"
49	#include "MCalibrationChargeCam.h"
50	#include "MCalibrationChargePINDiode.h"
51	#include "MCalibrationChargeBlindPix.h"
52	#include "MCalibrationChargeCalc.h"
53
54	#include "MReadMarsFile.h"
55	#include "MGeomApply.h"
56	#include "MBadPixelsMerge.h"
57	#include "MExtractSignal.h"
58	#include "MExtractPINDiode.h"
59	#include "MExtractBlindPixel.h"
60	#include "MExtractSignal2.h"
61	#include "MFCosmics.h"
62	#include "MContinue.h"
63	#include "MFillH.h"
64
65	#include "MJCalibration.h"
66	#include "MStatusDisplay.h"
67
68	ClassImp(MJCalibration);
69	using namespace std;
70
71	MJCalibration::MJCalibration(const char name, const char title) : fRuns(0)
72	{
73	fName = name ? name : "MJCalibration";
74	fTitle = title ? title : "Tool to create a pedestal file (MPedestalCam)";
75	}
76
77	void MJCalibration::DrawProjection(MHCamera *obj1, Int_t fit) const
78	{
79
80	TH1D obj2 = (TH1D)obj1->Projection(obj1->GetName());
81	obj2->Draw();
82	obj2->SetBit(kCanDelete);
83
84	const Double_t min = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst());
85	const Double_t max = obj2->GetBinCenter(obj2->GetXaxis()->GetLast());
86	const Double_t integ = obj2->Integral("width")/2.5;
87	const Double_t mean = obj2->GetMean();
88	const Double_t rms = obj2->GetRMS();
89	const Double_t width = max-min;
90
91	const TString dgausformula = "([0]-[3])/[2]exp(-0.5(x-[1])*(x-[1])/[2]/[2])"
92	"+[3]/[5]exp(-0.5(x-[4])*(x-[4])/[5]/[5])";
93
94	const TString tgausformula = "([0]-[3]-[6])/[2]exp(-0.5(x-[1])*(x-[1])/[2]/[2])"
95	"+[3]/[5]exp(-0.5(x-[4])*(x-[4])/[5]/[5])"
96	"+[6]/[8]exp(-0.5(x-[7])*(x-[7])/[8]/[8])";
97	TF1 *f=0;
98	switch (fit)
99	{
100	case 1:
101	f = new TF1("sgaus", "gaus(0)", min, max);
102	f->SetLineColor(kYellow);
103	f->SetBit(kCanDelete);
104	f->SetParNames("Area", "#mu", "#sigma");
105	f->SetParameters(integ/rms, mean, rms);
106	f->SetParLimits(0, 0, integ);
107	f->SetParLimits(1, min, max);
108	f->SetParLimits(2, 0, width/1.5);
109
110	obj2->Fit(f, "QLR");
111	break;
112
113	case 2:
114	f = new TF1("dgaus",dgausformula.Data(),min,max);
115	f->SetLineColor(kYellow);
116	f->SetBit(kCanDelete);
117	f->SetParNames("A_{tot}", "#mu1", "#sigma1", "A2", "#mu2", "#sigma2");
118	f->SetParameters(integ,(min+mean)/2.,width/4.,
119	integ/width/2.,(max+mean)/2.,width/4.);
120	// The left-sided Gauss
121	f->SetParLimits(0,integ-1.5 , integ+1.5);
122	f->SetParLimits(1,min+(width/10.), mean);
123	f->SetParLimits(2,0 , width/2.);
124	// The right-sided Gauss
125	f->SetParLimits(3,0 , integ);
126	f->SetParLimits(4,mean, max-(width/10.));
127	f->SetParLimits(5,0 , width/2.);
128	obj2->Fit(f,"QLRM");
129	break;
130
131	case 3:
132	f = new TF1("tgaus",tgausformula.Data(),min,max);
133	f->SetLineColor(kYellow);
134	f->SetBit(kCanDelete);
135	f->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}",
136	"A_{2}","#mu_{2}","#sigma_{2}",
137	"A_{3}","#mu_{3}","#sigma_{3}");
138	f->SetParameters(integ,(min+mean)/2,width/4.,
139	integ/width/3.,(max+mean)/2.,width/4.,
140	integ/width/3.,mean,width/2.);
141	// The left-sided Gauss
142	f->SetParLimits(0,integ-1.5,integ+1.5);
143	f->SetParLimits(1,min+(width/10.),mean);
144	f->SetParLimits(2,width/15.,width/2.);
145	// The right-sided Gauss
146	f->SetParLimits(3,0.,integ);
147	f->SetParLimits(4,mean,max-(width/10.));
148	f->SetParLimits(5,width/15.,width/2.);
149	// The Gauss describing the outliers
150	f->SetParLimits(6,0.,integ);
151	f->SetParLimits(7,min,max);
152	f->SetParLimits(8,width/4.,width/1.5);
153	obj2->Fit(f,"QLRM");
154	break;
155
156	case 4:
157	obj2->Fit("pol0", "Q");
158	obj2->GetFunction("pol0")->SetLineColor(kYellow);
159	break;
160
161	case 9:
162	break;
163
164	default:
165	obj2->Fit("gaus", "Q");
166	obj2->GetFunction("gaus")->SetLineColor(kYellow);
167	break;
168	}
169	}
170
171	void MJCalibration::CamDraw(TCanvas &c, const Int_t x, const Int_t y, const MHCamera &cam1, const Int_t fit)
172	{
173	c.cd(x);
174	gPad->SetBorderMode(0);
175	MHCamera obj1=(MHCamera)cam1.DrawCopy("hist");
176	obj1->AddNotify(&fCalibrationCam);
177
178	c.cd(x+y);
179	gPad->SetBorderMode(0);
180	obj1->Draw();
181
182	if (!fit)
183	return;
184
185	c.cd(x+2*y);
186	gPad->SetBorderMode(0);
187	DrawProjection(obj1, fit);
188	}
189
190
191	void MJCalibration::DisplayResult(MParList &plist)
192	{
193	if (!fDisplay)
194	return;
195
196	//
197	// Update display
198	//
199	TString title = fDisplay->GetTitle();
200	title += "-- Calibration ";
201	title += fRuns->GetRunsAsString();
202	title += " --";
203	fDisplay->SetTitle(title);
204
205	//
206	// Get container from list
207	//
208	MGeomCam &geomcam = (MGeomCam)plist.FindObject("MGeomCam");
209
210	// Create histograms to display
211	MHCamera disp1 (geomcam, "Cal;Charge", "Fitted Mean Charges");
212	MHCamera disp2 (geomcam, "Cal;SigmaCharge", "Sigma of Fitted Charges");
213	MHCamera disp3 (geomcam, "Cal;FitProb", "Probability of Fit");
214	MHCamera disp4 (geomcam, "Cal;RSigma", "Reduced Sigmas");
215	MHCamera disp5 (geomcam, "Cal;RSigma/Charge", "Reduced Sigma per Charge");
216	MHCamera disp6 (geomcam, "Cal;FFactorPhe", "Nr. of Phe's (F-Factor Method)");
217	MHCamera disp7 (geomcam, "Cal;FFactorConv", "Conversion Factor (F-Factor Method)");
218	MHCamera disp8 (geomcam, "Cal;FFactorFFactor", "Total F-Factor (F-Factor Method)");
219	MHCamera disp9 (geomcam, "Cal;BlindPixPh", "Photon flux inside plexiglass (Blind Pixel Method)");
220	MHCamera disp10(geomcam, "Cal;BlindPixConv", "Conversion Factor (Blind Pixel Method)");
221	MHCamera disp11(geomcam, "Cal;BlindPixFFactor","Total F-Factor (Blind Pixel Method)");
222	MHCamera disp12(geomcam, "Cal;PINDiodePh", "Photons flux outside plexiglass (PIN Diode Method)");
223	MHCamera disp13(geomcam, "Cal;PINDiodeConv", "Conversion Factor (PIN Diode Method)");
224	MHCamera disp14(geomcam, "Cal;PINDiodeFFactor","Total F-Factor (PIN Diode Method)");
225	MHCamera disp15(geomcam, "Cal;Excluded", "Pixels previously excluded");
226	MHCamera disp16(geomcam, "Cal;NotFitted", "Pixels that could not be fitted");
227	MHCamera disp17(geomcam, "Cal;NotFitValid", "Pixels with not valid fit results");
228	MHCamera disp18(geomcam, "Cal;Oscillation", "Oscillating Pixels");
229	MHCamera disp19(geomcam, "Cal;Saturation", "Pixels with saturated Hi Gain");
230
231	// Fitted charge means and sigmas
232	disp1.SetCamContent(fCalibrationCam, 0);
233	disp1.SetCamError( fCalibrationCam, 1);
234	disp2.SetCamContent(fCalibrationCam, 2);
235	disp2.SetCamError( fCalibrationCam, 3);
236	// Fit probabilities
237	disp3.SetCamContent(fCalibrationCam, 4);
238
239	// Reduced Sigmas and reduced sigmas per charge
240	disp4.SetCamContent(fCalibrationCam, 5);
241	disp4.SetCamError( fCalibrationCam, 6);
242	disp5.SetCamContent(fCalibrationCam, 7);
243	disp5.SetCamError( fCalibrationCam, 8);
244
245	// F-Factor Method
246	disp6.SetCamContent(fCalibrationCam, 9);
247	disp6.SetCamError( fCalibrationCam, 10);
248	disp7.SetCamContent(fCalibrationCam, 11);
249	disp7.SetCamError( fCalibrationCam, 12);
250	disp8.SetCamContent(fCalibrationCam, 13);
251	disp8.SetCamError( fCalibrationCam, 14);
252
253	/// Blind Pixel Method
254	disp9.SetCamContent(fCalibrationCam, 15);
255	disp9.SetCamError( fCalibrationCam, 16);
256	disp10.SetCamContent(fCalibrationCam,17);
257	disp10.SetCamError( fCalibrationCam,18);
258	disp11.SetCamContent(fCalibrationCam,19);
259	disp11.SetCamError( fCalibrationCam,20);
260
261	// PIN Diode Method
262	disp12.SetCamContent(fCalibrationCam,21);
263	disp12.SetCamError( fCalibrationCam,22);
264	disp13.SetCamContent(fCalibrationCam,23);
265	disp13.SetCamError( fCalibrationCam,24);
266	disp14.SetCamContent(fCalibrationCam,25);
267	disp14.SetCamError( fCalibrationCam,26);
268
269	// Pixels with defects
270	disp15.SetCamContent(fCalibrationCam,27);
271	disp16.SetCamContent(fCalibrationCam,28);
272	disp17.SetCamContent(fCalibrationCam,29);
273	disp18.SetCamContent(fCalibrationCam,30);
274
275	// Lo Gain calibration
276	disp19.SetCamContent(fCalibrationCam,31);
277
278	disp1.SetYTitle("Q [FADC units]");
279	disp2.SetYTitle("\\sigma_{Q} [FADC units]");
280	disp3.SetYTitle("P_{Q} [1]");
281
282	disp4.SetYTitle("\\sqrt{\\sigma^{2}_{Q} - RMS^{2}_{Ped}} [FADC Counts]");
283	disp5.SetYTitle("Red.Sigma/<Q> [1]");
284
285	disp6.SetYTitle("PhE [#]");
286	disp7.SetYTitle("Conv.Factor [PhE/FADC units]");
287	disp8.SetYTitle("\\sqrt{N_{PhE}}*\\sigma_{Q}/\\mu_{Q} [1]");
288
289	disp9.SetYTitle("Phot.flux [ph/mm^{2}]");
290	disp10.SetYTitle("Conv.Factor [Phot/FADC Count]");
291	disp11.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Q}/\\mu_{Q} [1]");
292
293	disp12.SetYTitle("Phot.flux [ph/mm^{2}]");
294	disp13.SetYTitle("Conv.Factor [Phot/FADC Count]");
295	disp14.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Q}/\\mu_{Q} [1]");
296
297	disp15.SetYTitle("[1]");
298	disp16.SetYTitle("[1]");
299	disp17.SetYTitle("[1]");
300	disp18.SetYTitle("[1]");
301
302	gStyle->SetOptStat(1111);
303	gStyle->SetOptFit();
304
305	// Charges
306	TCanvas &c1 = fDisplay->AddTab("Fit.Charge");
307	c1.Divide(2, 3);
308
309	CamDraw(c1, 1, 2, disp1, 2);
310	CamDraw(c1, 2, 2, disp2, 2);
311
312	// Fit Probability
313	TCanvas &c2 = fDisplay->AddTab("Fit.Prob");
314	c2.Divide(1,3);
315
316	CamDraw(c2, 1, 1, disp3, 4);
317
318	// Reduced Sigmas
319	TCanvas &c3 = fDisplay->AddTab("Red.Sigma");
320	c3.Divide(2,3);
321
322	CamDraw(c3, 1, 2, disp4, 2);
323	CamDraw(c3, 2, 2, disp5, 2);
324
325	// F-Factor Method
326	TCanvas &c4 = fDisplay->AddTab("F-Factor");
327	c4.Divide(3,3);
328
329	CamDraw(c4, 1, 3, disp6, 2);
330	CamDraw(c4, 2, 3, disp7, 2);
331	CamDraw(c4, 3, 3, disp8, 2);
332
333	// Blind Pixel Method
334	TCanvas &c5 = fDisplay->AddTab("BlindPix");
335	c5.Divide(3, 3);
336
337	CamDraw(c5, 1, 3, disp9, 9);
338	CamDraw(c5, 2, 3, disp10, 2);
339	CamDraw(c5, 3, 3, disp11, 2);
340
341	// PIN Diode Method
342	TCanvas &c6 = fDisplay->AddTab("PINDiode");
343	c6.Divide(3,3);
344
345	CamDraw(c6, 1, 3, disp12, 9);
346	CamDraw(c6, 2, 3, disp13, 2);
347	CamDraw(c6, 3, 3, disp14, 2);
348
349	// Defects
350	TCanvas &c7 = fDisplay->AddTab("Defects");
351	c7.Divide(4,2);
352
353	CamDraw(c7, 1, 4, disp15, 0);
354	CamDraw(c7, 2, 4, disp16, 0);
355	CamDraw(c7, 3, 4, disp17, 0);
356	CamDraw(c7, 4, 4, disp18, 0);
357
358	// Lo Gain Calibration
359	TCanvas &c8 = fDisplay->AddTab("LowGain");
360	c8.Divide(1,3);
361
362	CamDraw(c8, 1, 1, disp19, 0);
363	}
364
365	Bool_t MJCalibration::WriteResult()
366	{
367	if (fOutputPath.IsNull())
368	return kTRUE;
369
370	const TString oname(GetOutputFile());
371
372	*fLog << inf << "Writing to file: " << oname << endl;
373
374	TFile file(oname, "UPDATE");
375
376	if (fDisplay && fDisplay->Write()<=0)
377	{
378	*fLog << err << "Unable to write MStatusDisplay to " << oname << endl;
379	return kFALSE;
380	}
381
382	if (fCalibrationCam.Write()<=0)
383	{
384	*fLog << err << "Unable to write MCalibrationCam to " << oname << endl;
385	return kFALSE;
386	}
387
388	if (fBadPixels.Write()<=0)
389	{
390	*fLog << err << "Unable to write MBadPixelsCam to " << oname << endl;
391	return kFALSE;
392	}
393
394	return kTRUE;
395
396	}
397
398	void MJCalibration::SetOutputPath(const char *path)
399	{
400	fOutputPath = path;
401	if (fOutputPath.EndsWith("/"))
402	fOutputPath = fOutputPath(0, fOutputPath.Length()-1);
403	}
404
405	Bool_t MJCalibration::Process(MPedestalCam &pedcam)
406	{
407	if (!ReadCalibrationCam())
408	return ProcessFile(pedcam);
409
410	return kTRUE;
411	}
412
413	TString MJCalibration::GetOutputFile() const
414	{
415	if (!fRuns)
416	return "";
417
418	return Form("%s/%s-F1.root", (const char)fOutputPath, (const char)fRuns->GetRunsAsFileName());
419	}
420
421	Bool_t MJCalibration::ReadCalibrationCam()
422	{
423	const TString fname = GetOutputFile();
424
425	if (gSystem->AccessPathName(fname, kFileExists))
426	{
427	*fLog << err << "Input file " << fname << " doesn't exist." << endl;
428	return kFALSE;
429	}
430
431	*fLog << inf << "Reading from file: " << fname << endl;
432
433	TFile file(fname, "READ");
434	if (fCalibrationCam.Read()<=0)
435	{
436	*fLog << "Unable to read MCalibrationCam from " << fname << endl;
437	return kFALSE;
438	}
439
440	if (file.FindKey("MBadPixelsCam"))
441	{
442	MBadPixelsCam bad;
443	if (bad.Read()<=0)
444	{
445	*fLog << "Unable to read MBadPixelsCam from " << fname << endl;
446	return kFALSE;
447	}
448	fBadPixels.Merge(bad);
449	}
450
451	if (fDisplay /&& !fDisplay->GetCanvas("Pedestals")/) // FIXME!
452	fDisplay->Read();
453
454	return kTRUE;
455	}
456
457
458	Bool_t MJCalibration::ProcessFile(MPedestalCam &pedcam)
459	{
460	if (!fRuns)
461	{
462	*fLog << err << "No Runs choosen... abort." << endl;
463	return kFALSE;
464	}
465	if (fRuns->GetNumRuns() != fRuns->GetNumEntries())
466	{
467	*fLog << err << "Number of files found doesn't match number of runs... abort." << endl;
468	return kFALSE;
469	}
470
471	*fLog << inf;
472	fLog->Separator(GetDescriptor());
473	*fLog << "Calculate MCalibrationCam from Runs " << fRuns->GetRunsAsString() << endl;
474	*fLog << endl;
475
476	MReadMarsFile read("Events");
477	read.DisableAutoScheme();
478	static_cast<MRead&>(read).AddFiles(*fRuns);
479
480	//
481	// As long, as we don't have digital modules,
482	// we have to set the color of the pulser LED by hand
483	//
484	MCalibrationChargePINDiode pindiode;
485	pindiode.SetColor(kCT1);
486
487	MCalibrationChargeBlindPix blindpix;
488	blindpix.SetColor(kCT1);
489
490	// Setup Tasklist
491	MParList plist;
492	plist.AddToList(&pedcam);
493	plist.AddToList(&fCalibrationCam);
494	plist.AddToList(&pindiode);
495	plist.AddToList(&blindpix);
496
497	MTaskList tlist;
498	plist.AddToList(&tlist);
499
500	MGeomApply apply;
501	MBadPixelsMerge merge(&fBadPixels);
502	// MExtractSignal extract; // Do not use this at the moment...
503	MExtractPINDiode pinext;
504	MExtractBlindPixel blindext;
505	MExtractSignal2 extract;
506	MCalibrationChargeCalc calcalc;
507
508	MFillH fillpin( "MHCalibrationChargePINDiode" , "MExtractedSignalPINDiode" );
509	MFillH fillblind("MHCalibrationChargeBlindPix" , "MExtractedSignalBlindPixel" );
510	MFillH fillcam ("MHCalibrationChargeCam" , "MExtractedSignalCam" );
511	//
512	// Apply a filter against cosmics
513	// (will have to be needed in the future
514	// when the calibration hardware-trigger is working)
515	//
516	MFCosmics cosmics;
517	MContinue cont(&cosmics);
518
519	//calcalc.SkipBlindPixelFit();
520
521	tlist.AddToList(&read);
522	tlist.AddToList(&apply);
523	tlist.AddToList(&merge);
524	tlist.AddToList(&extract);
525	tlist.AddToList(&pinext);
526	tlist.AddToList(&blindext);
527	tlist.AddToList(&cont);
528	tlist.AddToList(&fillcam);
529	tlist.AddToList(&fillpin);
530	tlist.AddToList(&fillblind);
531	tlist.AddToList(&calcalc);
532
533	// Create and setup the eventloop
534	MEvtLoop evtloop(fName);
535	evtloop.SetParList(&plist);
536	evtloop.SetDisplay(fDisplay);
537	evtloop.SetLogStream(fLog);
538
539	// Execute first analysis
540	if (!evtloop.Eventloop())
541	{
542	*fLog << err << GetDescriptor() << ": Failed." << endl;
543	return kFALSE;
544	}
545
546	tlist.PrintStatistics();
547
548	DisplayResult(plist);
549
550	if (!WriteResult())
551	return kFALSE;
552
553	*fLog << inf << GetDescriptor() << ": Done." << endl;
554
555	return kTRUE;
556	}

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

Download in other formats: