Context Navigation

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

Visit:

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

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

Download in other formats: