Context Navigation

source: trunk/MagicSoft/Mars/macros/pedphotcalc.C@ 3802

Visit:

Last change on this file since 3802 was 3760, checked in by gaug, 22 years ago
* empty log message *
File size: 18.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	! Author(s): Josep Flix, 01/2004 <mailto:jflix@ifae.es>
18	! Javier Rico, 01/2004 <mailto:jrico@ifae.es>
19	! Markus Gaug, 03/2004 <mailto:markus@ifae.es>
20	!
21	! (based on bootcampstandardanalysis.C by Javier López)
22	!
23	!
24	! Copyright: MAGIC Software Development, 2000-2004
25	!
26	!
27	\* ======================================================================== */
28	/////////////////////////////////////////////////////////////////////////////
29	//
30	// pedphotcalc.C
31	//
32	// This macro is an example of the use of MPedPhotCalc. It computes
33	// the pedestal mean and rms from pedestal files undergoing the
34	// signal extraction + calibration chain, in units of photons. It
35	// produces plots of the relevant computed quantities.
36	//
37	// Needs as arguments the run number of a pedestal file ("_P_.root"),
38	// one of a calibration file ("_C_.root").
39	// performs the pedestal calculation, the calibration
40	/// constants calculation and the calibration of the pedestals.
41	//
42	// The TString inpath has to be set correctly.
43	//
44	// The macro searches for the pulser colour which corresponds to the calibration
45	// run number. If the run number is smaller than 20000, pulser colour "CT1"
46	// is assumed, otherwise, it searches for the strings "green", "blue", "uv" or
47	// "ct1" in the filenames. If no colour or multiple colours are found, the
48	// execution is aborted.
49	//
50	////////////////////////////////////////////////////////////////////////////////////
51	#include "MParList.h"
52	#include "MTaskList.h"
53	#include "MJPedestal.h"
54	#include "MJCalibration.h"
55	#include "MPedestalCam.h"
56	#include "MPedestalPix.h"
57	#include "MReadMarsFile.h"
58	#include "MGeomApply.h"
59	#include "MGeomCamMagic.h"
60	#include "MEvtLoop.h"
61	#include "MCalibrationCam.h"
62	#include "MCalibrationChargeCam.h"
63	#include "MCalibrationChargePix.h"
64	#include "MCalibrationQECam.h"
65	#include "MCalibrationQEPix.h"
66	#include "MExtractedSignalCam.h"
67	#include "MExtractSignal.h"
68	#include "MCerPhotEvt.h"
69	#include "MCalibrate.h"
70	#include "MPedPhotCalc.h"
71	#include "MPedPhotCam.h"
72	#include "MPedPhotPix.h"
73	#include "MHCamera.h"
74	#include "MRunIter.h"
75	#include "MDirIter.h"
76	#include "MStatusDisplay.h"
77	#include "MHCamera.h"
78
79	#include "TTimer.h"
80	#include "TString.h"
81	#include "TCanvas.h"
82	#include "TStyle.h"
83	#include "TF1.h"
84	#include "TLegend.h"
85
86	#include <iostream.h>
87	#include "Getline.h"
88
89	const TString inpath = "/mnt/Data/rootdata/CrabNebula/2004_02_10/";
90	const Int_t dpedrun = 14607;
91	const Int_t dcalrun1 = 14608;
92	const Int_t dcalrun2 = 0;
93	const Bool_t usedisplay = kTRUE;
94
95	static MCalibrationCam::PulserColor_t FindColor(MDirIter* run)
96	{
97
98	MCalibrationCam::PulserColor_t col = MCalibrationCam::kNONE;
99
100	TString filenames;
101
102	while (!(filenames=run->Next()).IsNull())
103	{
104
105	filenames.ToLower();
106
107	if (filenames.Contains("green"))
108	if (col == MCalibrationCam::kNONE)
109	{
110	cout << "Found colour: Green in " << filenames << endl;
111	col = MCalibrationCam::kGREEN;
112	}
113	else if (col != MCalibrationCam::kGREEN)
114	{
115	cout << "Different colour found in " << filenames << "... abort" << endl;
116	return MCalibrationCam::kNONE;
117	}
118
119	if (filenames.Contains("blue"))
120	if (col == MCalibrationCam::kNONE)
121	{
122	cout << "Found colour: Blue in " << filenames << endl;
123	col = MCalibrationCam::kBLUE;
124	}
125	else if (col != MCalibrationCam::kBLUE)
126	{
127	cout << "Different colour found in " << filenames << "... abort" << endl;
128	return MCalibrationCam::kNONE;
129	}
130
131	if (filenames.Contains("uv"))
132	if (col == MCalibrationCam::kNONE)
133	{
134	cout << "Found colour: Uv in " << filenames << endl;
135	col = MCalibrationCam::kUV;
136	}
137	else if (col != MCalibrationCam::kUV)
138	{
139	cout << "Different colour found in " << filenames << "... abort" << endl;
140	return MCalibrationCam::kNONE;
141	}
142
143	if (filenames.Contains("ct1"))
144	if (col == MCalibrationCam::kNONE)
145	{
146	cout << "Found colour: Ct1 in " << filenames << endl;
147	col = MCalibrationCam::kCT1;
148	}
149	else if (col != MCalibrationCam::kCT1)
150	{
151	cout << "Different colour found in " << filenames << "... abort" << endl;
152	return MCalibrationCam::kNONE;
153	}
154
155	}
156
157
158
159	if (col == MCalibrationCam::kNONE)
160	cout << "No colour found in filenames of runs: " << ((MRunIter*)run)->GetRunsAsString()
161	<< "... abort" << endl;
162
163	return col;
164	}
165
166
167
168	void DrawProjection(MHCamera *obj1, Int_t fit)
169	{
170	TH1D obj2 = (TH1D)obj1->Projection(obj1->GetName());
171	obj2->SetDirectory(0);
172	obj2->Draw();
173	obj2->SetBit(kCanDelete);
174	gPad->SetLogy();
175
176	if (obj1->GetGeomCam().InheritsFrom("MGeomCamMagic"))
177	{
178	TArrayI s0(3);
179	s0[0] = 6;
180	s0[1] = 1;
181	s0[2] = 2;
182
183	TArrayI s1(3);
184	s1[0] = 3;
185	s1[1] = 4;
186	s1[2] = 5;
187
188	TArrayI inner(1);
189	inner[0] = 0;
190
191	TArrayI outer(1);
192	outer[0] = 1;
193
194	// Just to get the right (maximum) binning
195	TH1D *half[4];
196	half[0] = obj1->ProjectionS(s0, inner, "Sector 6-1-2 Inner");
197	half[1] = obj1->ProjectionS(s1, inner, "Sector 3-4-5 Inner");
198	half[2] = obj1->ProjectionS(s0, outer, "Sector 6-1-2 Outer");
199	half[3] = obj1->ProjectionS(s1, outer, "Sector 3-4-5 Outer");
200
201	for (int i=0; i<4; i++)
202	{
203	half[i]->SetLineColor(kRed+i);
204	half[i]->SetDirectory(0);
205	half[i]->SetBit(kCanDelete);
206	half[i]->Draw("same");
207	}
208	}
209
210	const Double_t min = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst());
211	const Double_t max = obj2->GetBinCenter(obj2->GetXaxis()->GetLast());
212	const Double_t integ = obj2->Integral("width")/2.5066283;
213	const Double_t mean = obj2->GetMean();
214	const Double_t rms = obj2->GetRMS();
215	const Double_t width = max-min;
216
217	if (rms==0 \|\| width==0)
218	return;
219
220	const TString dgausformula("([0]-[3])/[2]exp(-0.5(x-[1])*(x-[1])/[2]/[2])"
221	"+[3]/[5]exp(-0.5(x-[4])*(x-[4])/[5]/[5])");
222
223	TF1 *f=0;
224	switch (fit)
225	{
226	// FIXME: MAYBE add function to TH1?
227	case 0:
228	f = new TF1("sgaus", "gaus(0)", min, max);
229	f->SetLineColor(kYellow);
230	f->SetBit(kCanDelete);
231	f->SetParNames("Area", "#mu", "#sigma");
232	f->SetParameters(integ/rms, mean, rms);
233	f->SetParLimits(0, 0, integ);
234	f->SetParLimits(1, min, max);
235	f->SetParLimits(2, 0, width/1.5);
236	obj2->Fit(f, "QLRM");
237	break;
238
239	case 1:
240	f = new TF1("dgaus", dgausformula, min, max);
241	f->SetLineColor(kYellow);
242	f->SetBit(kCanDelete);
243	f->SetParNames("A_{tot}", "#mu_{1}", "#sigma_{1}", "A_{2}", "#mu_{2}", "#sigma_{2}");
244	f->SetParameters(integ, (min+mean)/2, width/4,
245	integ/width/2, (max+mean)/2, width/4);
246	// The left-sided Gauss
247	f->SetParLimits(0, integ-1.5, integ+1.5);
248	f->SetParLimits(1, min+(width/10), mean);
249	f->SetParLimits(2, 0, width/2);
250	// The right-sided Gauss
251	f->SetParLimits(3, 0, integ);
252	f->SetParLimits(4, mean, max-(width/10));
253	f->SetParLimits(5, 0, width/2);
254	obj2->Fit(f, "QLRM");
255	break;
256
257	default:
258	obj2->Fit("gaus", "Q");
259	obj2->GetFunction("gaus")->SetLineColor(kYellow);
260	break;
261	}
262	}
263
264	void CamDraw(TCanvas &c, Int_t x, Int_t y, MHCamera &cam1, Int_t fit)
265	{
266	c.cd(x);
267	gPad->SetBorderMode(0);
268	MHCamera obj1=(MHCamera)cam1.DrawCopy("hist");
269
270	c.cd(x+y);
271	gPad->SetBorderMode(0);
272	obj1->Draw();
273
274	c.cd(x+2*y);
275	gPad->SetBorderMode(0);
276	DrawProjection(obj1, fit);
277	}
278
279	void pedphotcalc(const Int_t prun=dpedrun, // pedestal file
280	const Int_t crun1=dcalrun1, const Int_t crun2=dcalrun2 // calibration file(s)
281	)
282	{
283
284	MRunIter pruns;
285	MRunIter cruns;
286
287	pruns.AddRun(prun,inpath);
288
289	if (crun2==0)
290	cruns.AddRun(crun1,inpath);
291	else
292	cruns.AddRuns(crun1,crun2,inpath);
293
294	MCalibrationCam::PulserColor_t color;
295
296	if (crun1 < 20000)
297	color = MCalibrationCam::kCT1;
298	else
299	color = FindColor((MDirIter*)&cruns);
300
301	if (color == MCalibrationCam::kNONE)
302	{
303	TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
304
305	while (1)
306	{
307	timer.TurnOn();
308	TString input = Getline("Could not find the correct colour: Type 'q' to exit, "
309	"green, blue, uv or ct1 to go on: ");
310	timer.TurnOff();
311
312	if (input=="q\n")
313	return ;
314
315	if (input=="green")
316	color = MCalibrationCam::kGREEN;
317	if (input=="blue")
318	color = MCalibrationCam::kBLUE;
319	if (input=="uv")
320	color = MCalibrationCam::kUV;
321	if (input=="ct1")
322	color = MCalibrationCam::kCT1;
323	}
324	}
325
326	//
327	// Now setup the tasks and tasklist for the pedestals:
328	// ---------------------------------------------------
329	//
330	MGeomCamMagic geomcam;
331	MGeomApply geomapl;
332	MStatusDisplay *display = NULL;
333
334	/************************************/
335	/* FIRST LOOP: PEDESTAL COMPUTATION */
336	/************************************/
337
338	MJPedestal pedloop;
339	pedloop.SetInput(&pruns);
340	if (usedisplay)
341	{
342	display = new MStatusDisplay;
343	display->SetUpdateTime(3000);
344	display->Resize(850,700);
345	display->SetBit(kCanDelete);
346	pedloop.SetDisplay(display);
347	}
348
349	cout << "*************************" << endl;
350	cout << " COMPUTING PEDESTALS " << endl;
351	cout << "*************************" << endl;
352
353	if (!pedloop.Process())
354	return;
355
356	MPedestalCam pedcam = pedloop.GetPedestalCam();
357
358	/****************************/
359	/* SECOND LOOP: CALIBRATION */
360	/****************************/
361
362	//
363	// Now setup the new tasks for the calibration:
364	// ---------------------------------------------------
365	//
366	MJCalibration calloop;
367	calloop.SetColor(color);
368	calloop.SetInput(&cruns);
369	//
370	// Use as signal extractor MExtractSignal:
371	//
372	calloop.SetExtractorLevel(1);
373	//
374	// The next two commands are for the display:
375	//
376	if (usedisplay)
377	{
378	calloop.SetDisplay(display);
379	calloop.SetDataCheck();
380	}
381
382	cout << "***********************************" << endl;
383	cout << " COMPUTING CALIBRATION FACTORS " << endl;
384	cout << "***********************************" << endl;
385
386	if (!calloop.Process(pedcam))
387	return;
388
389	MCalibrationChargeCam &calcam = calloop.GetCalibrationCam();
390	MCalibrationQECam &qecam = calloop.GetQECam();
391
392	/************************************************************************/
393	/* THIRD LOOP: PEDESTAL COMPUTATION USING EXTRACTED SIGNAL (IN PHOTONS) */
394	/************************************************************************/
395
396	// Create an empty Parameter List and an empty Task List
397	MParList plist3;
398	MTaskList tlist3;
399	plist3.AddToList(&tlist3);
400
401	// containers
402	MCerPhotEvt photcam;
403	MPedPhotCam pedphotcam;
404	MExtractedSignalCam extedsig;
405
406	plist3.AddToList(&geomcam);
407	plist3.AddToList(&pedcam );
408	plist3.AddToList(&calcam );
409	plist3.AddToList(&qecam );
410	plist3.AddToList(&photcam);
411	plist3.AddToList(&extedsig);
412	plist3.AddToList(&pedphotcam);
413
414	//tasks
415	MReadMarsFile read3("Events");
416	read3.DisableAutoScheme();
417	static_cast<MRead&>(read3).AddFiles(pruns);
418
419	MExtractSignal sigcalc;
420	MCalibrate photcalc;
421	photcalc.SetCalibrationMode(MCalibrate::kFfactor);
422	MPedPhotCalc pedphotcalc;
423
424	tlist3.AddToList(&read3);
425	tlist3.AddToList(&geomapl);
426	tlist3.AddToList(&sigcalc);
427	tlist3.AddToList(&photcalc);
428	tlist3.AddToList(&pedphotcalc);
429
430	// Create and execute eventloop
431	MEvtLoop evtloop3;
432	evtloop3.SetParList(&plist3);
433
434	cout << "*************************************************************" << endl;
435	cout << " COMPUTING PEDESTALS USING EXTRACTED SIGNAL (IN PHOTONS) " << endl;
436	cout << "*************************************************************" << endl;
437
438	if (!evtloop3.Eventloop())
439	return;
440	tlist3.PrintStatistics();
441
442	/**********************/
443	/* PRODUCE NICE PLOTS */
444	/**********************/
445
446	if (usedisplay)
447	{
448
449	MHCamera disp0(geomcam, "MPedPhotCam;ped", "Pedestals");
450	MHCamera disp1(geomcam, "MPedPhotCam;rms", "Sigma Pedestal");
451
452	disp0.SetCamContent(pedphotcam, 0);
453	disp0.SetCamError (pedphotcam, 1);
454
455	disp1.SetCamContent(pedphotcam, 1);
456
457	disp0.SetYTitle("Pedestal signal [photons]");
458	disp1.SetYTitle("Pedestal signal RMS [photons]");
459
460	//
461	// Display data
462	//
463	TCanvas &c1 = display->AddTab("PedPhotCam");
464	c1.Divide(2,3);
465
466	CamDraw(c1, 1, 2, disp0, 0);
467	CamDraw(c1, 2, 2, disp1, 1);
468
469	}
470
471
472	#if 0
473	const UShort_t npix = 577;
474
475	// declare histograms
476	// pedestals
477	TH1F* pedhist = new TH1F("pedhist","Pedestal",npix,0,npix);
478	TH1F* pedrmshist = new TH1F("pedrmshist","Pedestal RMS",npix,0,npix);
479	TH1F* peddist = new TH1F("peddist","Pedestal",100,0,20);
480	TH1F* pedrmsdist = new TH1F("pedrmsdist","Pedestal RMS",100,0,15);
481
482	// calibration factors
483	TH1F* calhist = new TH1F("calhist","Calibration factors",npix,0,npix);
484	TH1F* caldist = new TH1F("caldist","Calibration factors",100,0,1);
485	TH1F* qehist = new TH1F("qehist", "Quantrum efficiencies",npix,0,npix);
486	TH1F* qedist = new TH1F("qedist", "Quantrum efficiencies",100,0,1);
487
488	// pedestal signals
489	TH1F* pedphothist = new TH1F("pedphothist","Pedestal Signal",npix,0,npix);
490	TH1F* pedphotrmshist = new TH1F("pedphotrmshist","Pedestal Signal RMS",npix,0,npix);
491	TH1F* pedphotdist = new TH1F("pedphotdist","Pedestal Signal",100,-0.4,0.4);
492	TH1F* pedphotrmsdist = new TH1F("pedphotrmsdist","Pedestal Signal RMS",100,0,25);
493
494	// fill histograms
495	for(int i=0;i<npix;i++)
496	{
497	MCalibrationChargePix &calpix = (MCalibrationChargePix&)calcam[i];
498	MCalibrationQEPix &qepix = (MCalibrationQEPix&) qecam[i];
499
500	const Float_t ped = pedcam[i].GetPedestal();
501	const Float_t pedrms = pedcam[i].GetPedestalRms();
502	const Float_t cal = calpix.GetMeanConvFADC2Phe();
503	const Float_t qe = qepix .GetQECascadesFFactor();
504	const Float_t pedphot = pedphotcam[i].GetMean();
505	const Float_t pedphotrms = pedphotcam[i].GetRms();
506
507	pedhist->Fill(i,ped);
508	peddist->Fill(ped);
509	pedrmshist->Fill(i,pedrms);
510	pedrmsdist->Fill(pedrms);
511
512	calhist->Fill(i,cal);
513	caldist->Fill(cal);
514	qehist->Fill(i,qe);
515	qedist->Fill(qe);
516
517	pedphothist->Fill(i,pedphot);
518	pedphotdist->Fill(pedphot);
519	pedphotrmshist->Fill(i,pedphotrms);
520	pedphotrmsdist->Fill(pedphotrms);
521	}
522
523	// Draw
524	gROOT->Reset();
525	gStyle->SetCanvasColor(0);
526	TCanvas* canvas = new TCanvas("canvas","pedphotcalc.C", 0, 100, 650, 800);
527	canvas->SetBorderMode(0); // Delete the Canvas' border line
528	canvas->cd();
529
530	canvas->Divide(2,5);
531
532	// draw pedestal histo
533	canvas->cd(1);
534	gPad->cd();
535	gPad->SetBorderMode(0);
536
537	pedhist->SetStats(kFALSE);
538	pedhist->GetXaxis()->SetTitle("Pixel SW Id");
539	pedhist->GetYaxis()->SetTitle("Pedestal (ADC counts)");
540	pedrmshist->SetStats(kFALSE);
541	pedrmshist->SetLineColor(2);
542	pedhist->Draw();
543	pedrmshist->Draw("same");
544
545	TLegend* leg1 = new TLegend(.14,.68,.39,.88);
546	leg1->SetHeader("");
547	leg1->AddEntry(pedhist,"Pedestal","L");
548	leg1->AddEntry(pedrmshist,"Pedestal RMS","L");
549	leg1->SetFillColor(0);
550	leg1->SetLineColor(0);
551	leg1->SetBorderSize(0);
552	leg1->Draw();
553
554	// draw pedestal distribution
555	canvas->cd(2);
556	gPad->cd();
557	gPad->SetBorderMode(0);
558	peddist->GetXaxis()->SetTitle("Pedestal (ADC counts)");
559	pedrmsdist->SetLineColor(2);
560	peddist->Draw();
561	pedrmsdist->Draw("same");
562
563	TLegend* leg2 = new TLegend(.14,.68,.39,.88);
564	leg2->SetHeader("");
565	leg2->AddEntry(pedhist,"Pedestal","L");
566	leg2->AddEntry(pedrmshist,"Pedestal RMS","L");
567	leg2->SetFillColor(0);
568	leg2->SetLineColor(0);
569	leg2->SetBorderSize(0);
570	leg2->Draw();
571
572	// draw calibration histo
573	canvas->cd(3);
574	gPad->cd();
575	gPad->SetBorderMode(0);
576	calhist->GetXaxis()->SetTitle("Pixel SW Id");
577	calhist->SetMaximum(1);
578	calhist->SetMinimum(0);
579	calhist->GetYaxis()->SetTitle("Calibration factor (phe/ADC count)");
580	calhist->SetStats(kFALSE);
581	calhist->Draw();
582
583	// draw calibration distribution
584	canvas->cd(4);
585	gPad->cd();
586	gPad->SetBorderMode(0);
587	caldist->GetXaxis()->SetTitle("Calibration factor (phe/ADC count)");
588	caldist->Draw();
589
590	// draw qe histo
591	canvas->cd(5);
592	gPad->cd();
593	gPad->SetBorderMode(0);
594	qehist->GetXaxis()->SetTitle("Pixel SW Id");
595	qehist->SetMaximum(1);
596	qehist->SetMinimum(0);
597	qehist->GetYaxis()->SetTitle("Quantum efficiency for cascades");
598	qehist->SetStats(kFALSE);
599	qehist->Draw();
600
601	// draw qe distribution
602	canvas->cd(6);
603	gPad->cd();
604	gPad->SetBorderMode(0);
605	qedist->GetXaxis()->SetTitle("Quantum efficiency for cascades");
606	qedist->Draw();
607
608	// draw pedestal signal histo
609	canvas->cd(7);
610	gPad->cd();
611	gPad->SetBorderMode(0);
612	pedphothist->GetXaxis()->SetTitle("Pixel SW Id");
613	pedphothist->GetYaxis()->SetTitle("Pedestal signal (photons)");
614	pedphothist->SetStats(kFALSE);
615	pedphothist->Draw();
616
617	// draw pedestal signal distribution
618	canvas->cd(8);
619	gPad->cd();
620	gPad->SetBorderMode(0);
621	pedphotdist->GetXaxis()->SetTitle("Pedestal signal (photons)");
622	pedphotdist->Draw();
623
624	// draw pedestal signal rms histo
625	canvas->cd(9);
626	gPad->cd();
627	gPad->SetBorderMode(0);
628	pedphotrmshist->GetXaxis()->SetTitle("Pixel SW Id");
629	pedphotrmshist->GetYaxis()->SetTitle("Pedestal signal rms (photons)");
630	pedphotrmshist->SetStats(kFALSE);
631	pedphotrmshist->Draw();
632
633	// draw pedestal signal rms distribution
634	canvas->cd(10);
635	gPad->cd();
636	gPad->SetBorderMode(0);
637	pedphotrmsdist->GetXaxis()->SetTitle("Pedestal signal rms (photons)");
638	pedphotrmsdist->Draw();
639
640	canvas->SaveAs("pedphotcalc.root");
641
642	#endif
643	}

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

Download in other formats: