source: trunk/MagicSoft/Mars/macros/calibration.C@ 2912

Last change on this file since 2912 was 2906, checked in by gaug, 23 years ago
*** empty log message ***
File size: 11.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): Markus Gaug, 11/2003 <mailto:markus@ifae.es>
19!
20! Copyright: MAGIC Software Development, 2000-2003
21!
22!
23\* ======================================================================== */
24
25const TString pedfile = "/mnt/Data/rootdata/Miscellaneous/2003_12_19/20031218_03522_P_Park_E.root";
26const TString calfile = "/mnt/Data/rootdata/Miscellaneous/2003_12_19/20031218_03527_C_Park_E.root";
27
28void calibration(TString pedname=pedfile,
29 TString calname=calfile)
30{
31
32 //
33 // Create a empty Parameter List and an empty Task List
34 // The tasklist is identified in the eventloop by its name
35 //
36 MParList plist;
37
38 MTaskList tlist;
39 plist.AddToList(&tlist);
40
41 //
42 // Now setup the tasks and tasklist for the pedestals:
43 // ---------------------------------------------------
44 //
45
46 MReadMarsFile read("Events", pedname);
47 read.DisableAutoScheme();
48
49 MGeomApply geomapl;
50 MPedCalcPedRun pedcalc;
51 MGeomCamMagic geomcam;
52 MPedestalCam pedcam;
53
54 tlist.AddToList(&read);
55 tlist.AddToList(&geomapl);
56 tlist.AddToList(&pedcalc);
57
58 plist.AddToList(&pedcam);
59
60 MHCamEvent hist("Pedestal");
61 hist.SetType(1);
62 plist.AddToList(&hist);
63 MFillH fill(&hist, "MPedestalCam");
64
65 tlist.AddToList(&fill);
66
67 //
68 // Create and setup the eventloop
69 //
70 MEvtLoop evtloop;
71 evtloop.SetParList(&plist);
72
73 //
74 // Execute first analysis
75 //
76 if (!evtloop.Eventloop())
77 return;
78
79 tlist.PrintStatistics();
80
81 //
82 // Create a empty Parameter List and an empty Task List
83 //
84 MParList plist2;
85 MTaskList tlist2;
86 plist2.AddToList(&tlist2);
87
88 MExtractedSignalCam sigcam;
89 MCalibrationCam calcam;
90 //
91 // Get the previously created MPedestalCam into the new Parameter List
92 //
93 plist2.AddToList(&pedcam);
94 plist2.AddToList(&sigcam);
95 plist2.AddToList(&calcam);
96
97 //
98 // Get the MAGIC geometry
99 //
100 tlist2.AddToList(&geomapl);
101 plist2.AddToList(&geomcam);
102
103 //
104 // Now setup the new tasks and tasklist for the calibration
105 // ---------------------------------------------------
106 //
107
108 MReadMarsFile read2("Events", calname);
109 read2.DisableAutoScheme();
110
111 MExtractSignal sigcalc;
112 MArrivalTimeCalc timecalc;
113 MCalibrationCalc calcalc;
114
115 //
116 // As long, as we don't have digital modules,
117 // we have to set the color of the pulser LED by hand
118 //
119 calcalc.SetPulserColor(MCalibrationCalc::kECT1);
120
121 //
122 // In case, we want to exclude a pre-defined list of bad pixels:
123 // (This is a preliminary feature)
124 //
125 // calcalc.ExcludePixelsFromAsciiFile("badpixels_all.dat");
126
127 //
128 // In case, you want to skip the blind pixel method:
129 // (NOT RECOMMENDED!!!)
130 //
131 // calcalc.SkipBlindPixelFit();
132
133 //
134 // In case, you want to skip the cosmics rejection
135 // (NOT RECOMMENDED!!!)
136 //
137 // calcalc.SkipCosmicsRejection();
138
139 //
140 // In case, you want to skip the quality checks
141 // (NOT RECOMMENDED!!!)
142 //
143 // calcalc.SkipQualityChecks();
144
145 //
146 // In case, we want to apply another fit function to the
147 // blind pixel
148 //
149 MCalibrationBlindPix *bp = calcam.GetBlindPixel();
150 bp->ChangeFitFunc(MHCalibrationBlindPixel::kEPoisson5);
151
152 tlist2.AddToList(&read2);
153 tlist2.AddToList(&sigcalc);
154 //
155 // In case, you want to skip the somewhat lengthy calculation
156 // of the arrival times using a spline, uncomment the next line
157 //
158 // tlist2.AddToList(&timecalc);
159 tlist2.AddToList(&calcalc);
160
161 //
162 // Create and setup the eventloop
163 //
164 MEvtLoop evtloop2;
165 evtloop2.SetParList(&plist2);
166
167 //
168 // Execute second analysis
169 //
170 if (!evtloop2.Eventloop())
171 return;
172
173 tlist2.PrintStatistics();
174
175 //
176 // print the most important results of all pixels
177 //
178 calcam.Print();
179
180 //
181 // just one example how to get the plots of individual pixels
182 //
183 calcam[17].DrawClone();
184
185 MHCamEvent camevt;
186 MHCamera disp1 (geomcam, "MCalibrationPix;Charge", "Fitted Mean Charges");
187 MHCamera disp3 (geomcam, "MCalibrationPix;SigmaCharge", "Sigma of Fitted Charges");
188 MHCamera disp5 (geomcam, "MCalibrationPix;ChargeProb", "Probability of Fit");
189 MHCamera disp6 (geomcam, "MCalibrationPix;Time", "Arrival Times");
190 MHCamera disp7 (geomcam, "MCalibrationPix;SigmaTime", "Sigma of Arrival Times");
191 MHCamera disp8 (geomcam, "MCalibrationPix;TimeChiSquare", "Chi Square of Time Fit");
192 MHCamera disp9 (geomcam, "MCalibrationPix;Ped", "Pedestals");
193 MHCamera disp10 (geomcam, "MCalibrationPix;PedRms", "Pedestal RMS");
194 MHCamera disp11 (geomcam, "MCalibrationPix;RSigma", "Reduced Sigmas");
195 MHCamera disp12 (geomcam, "MCalibrationPix;PheFFactorMethod", "Nr. of Phe's (F-Factor Method)");
196 MHCamera disp13 (geomcam, "MCalibrationPix;MeanConversionFFactorMethod",
197 "Conversion Factor (F-Factor Method)");
198 MHCamera disp14 (geomcam, "MCalibrationPix;MeanPhotInsidePlexiglass",
199 "Nr. of Photons (Blind Pixel Method)");
200 MHCamera disp15 (geomcam, "MCalibrationPix;MeanConversionBlindPixelMethod",
201 "Conversion Factor (Blind Pixel Method)");
202 MHCamera disp16 (geomcam, "MCalibrationPix;RSigma/Charge", "Reduced Sigma per Charge");
203
204 disp1.SetCamContent(calcam, 0);
205 disp1.SetCamError(calcam,1);
206
207 disp3.SetCamContent(calcam, 2);
208 disp3.SetCamError(calcam,3);
209
210 disp5.SetCamContent(calcam, 4);
211
212 disp6.SetCamContent(calcam, 5);
213 disp6.SetCamError(calcam, 6);
214 disp7.SetCamContent(calcam, 6);
215 disp8.SetCamContent(calcam, 7);
216
217 disp9.SetCamContent(calcam, 8);
218 disp9.SetCamError(calcam, 9);
219
220 disp10.SetCamContent(calcam, 9);
221 disp11.SetCamContent(calcam, 10);
222
223 disp12.SetCamContent(calcam, 11);
224 disp12.SetCamError(calcam, 12);
225
226 disp13.SetCamContent(calcam, 13);
227 disp13.SetCamError(calcam, 14);
228
229 disp14.SetCamContent(calcam, 15);
230 disp15.SetCamContent(calcam, 16);
231 disp16.SetCamContent(calcam, 17);
232
233
234 disp1.SetYTitle("Charge [FADC counts]");
235 disp3.SetYTitle("\\sigma_{Charge} [FADC counts]");
236 disp5.SetYTitle("P_{Charge} [1]");
237 disp6.SetYTitle("Arr. Time [Time Slice Nr.]");
238 disp7.SetYTitle("\\sigma_{Time} [Time Slices]");
239 disp8.SetYTitle("\\chi^{2}_{Time} [1]");
240 disp9.SetYTitle("Ped [FADC Counts ]");
241 disp10.SetYTitle("RMS_{Ped} [FADC Counts ]");
242 disp11.SetYTitle("\\sqrt{\\sigma^{2}_{Charge} - RMS^{2}_{Ped}} [FADC Counts]");
243 disp12.SetYTitle("Nr. Photo-Electrons [1]");
244 disp13.SetYTitle("Conversion Factor [PhE/FADC Count]");
245 disp14.SetYTitle("Nr. Photons [1]");
246 disp15.SetYTitle("Conversion Factor [Phot/FADC Count]");
247 disp16.SetYTitle("Reduced Sigma / Charge [1]");
248
249 MStatusDisplay *d3 = new MStatusDisplay;
250 d3->SetUpdateTime(3000);
251 d3->Resize(850,700);
252
253 gStyle->SetOptStat(1111);
254 gStyle->SetOptFit();
255
256 // Charges
257 TCanvas &c1 = d3->AddTab("Fitted Charges");
258 c1.Divide(2,3);
259
260 CamDraw(c1,disp1,&calcam,1,2,1);
261 CamDraw(c1,disp3,&calcam,2,2,2);
262
263 // Fit Probability
264 TCanvas &c2 = d3->AddTab("Fit Prob.");
265 c2.Divide(1,3);
266
267 CamDraw(c2,disp5,&calcam,1,1,3);
268
269 // Times
270 TCanvas &c3 = d3->AddTab("Fitted Times");
271 c3.Divide(3,3);
272
273 CamDraw(c3,disp6,&calcam,1,3,1);
274 CamDraw(c3,disp7,&calcam,2,3,0);
275 CamDraw(c3,disp8,&calcam,3,3,0);
276
277 // Pedestals
278 TCanvas &c4 = d3->AddTab("Pedestals");
279 c4.Divide(2,3);
280
281 CamDraw(c4,disp9,&calcam,1,2,0);
282 CamDraw(c4,disp10,&calcam,2,2,1);
283
284 // Reduced Sigmas
285 TCanvas &c5 = d3->AddTab("Reduced Sigmas");
286 c5.Divide(2,3);
287
288 CamDraw(c5,disp11,&calcam,1,2,2);
289 CamDraw(c5,disp16,&calcam,2,2,2);
290
291 // F-Factor Method
292 TCanvas &c6 = d3->AddTab("F-Factor Method");
293 c6.Divide(2,3);
294
295 CamDraw(c6,disp12,&calcam,1,2,1);
296 CamDraw(c6,disp13,&calcam,2,2,2);
297
298 // Blind Pixel Method
299 TCanvas &c7 = d3->AddTab("Blind Pixel Method");
300 c7.Divide(2, 3);
301
302 CamDraw(c7,disp14,&calcam,1,2,9);
303 CamDraw(c7,disp15,&calcam,2,2,2);
304
305}
306
307void CamDraw(TCanvas &c, MHCamera &cam, MCamEvent *evt, Int_t i, Int_t j, Int_t fit)
308{
309
310 c.cd(i);
311 gPad->SetBorderMode(0);
312 MHCamera *obj1=(MHCamera*)cam.DrawCopy("hist");
313 obj1->AddNotify(*evt);
314
315 c.cd(i+j);
316 gPad->SetBorderMode(0);
317 obj1->Draw();
318 ((MHCamera*)obj1)->SetPrettyPalette();
319
320 c.cd(i+2*j);
321 gPad->SetBorderMode(0);
322 TH1D *obj2 = (TH1D*)obj1->Projection();
323 obj2->Draw();
324 obj2->SetBit(kCanDelete);
325
326 const Double_t min = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst());
327 const Double_t max = obj2->GetBinCenter(obj2->GetXaxis()->GetLast());
328 const Double_t integ = obj2->Integral("width")/2.5;
329 const Double_t mean = obj2->GetMean();
330 const Double_t rms = obj2->GetRMS();
331 const Double_t width = max-min;
332
333 if (rms == 0. || width == 0. )
334 return;
335
336 switch (fit)
337 {
338 case 0:
339 TF1 *sgaus = new TF1("sgaus","gaus(0)",min,max);
340 sgaus->SetBit(kCanDelete);
341 sgaus->SetParNames("Area","#mu","#sigma");
342 sgaus->SetParameters(integ/rms,mean,rms);
343 sgaus->SetParLimits(0,0.,integ);
344 sgaus->SetParLimits(1,min,max);
345 sgaus->SetParLimits(2,0,width/1.5);
346 obj2->Fit("sgaus","QLR");
347 obj2->GetFunction("sgaus")->SetLineColor(kYellow);
348 break;
349
350 case 1:
351 TF1 *dgaus = new TF1("dgaus","gaus(0)+gaus(3)",min,max);
352 dgaus->SetBit(kCanDelete);
353 dgaus->SetParNames("A1","#mu1","#sigma1","A2","#mu2","#sigma2");
354 dgaus->SetParameters(integ/width,max-width/6.,width/4.,
355 integ/width,min+width/6.,width/4.);
356 dgaus->SetParLimits(0,0,integ);
357 dgaus->SetParLimits(1,min,max);
358 dgaus->SetParLimits(2,0,width/3.);
359 dgaus->SetParLimits(3,0,integ);
360 dgaus->SetParLimits(4,min,max);
361 dgaus->SetParLimits(5,0,width/3.);
362 obj2->Fit("dgaus","QLR");
363 obj2->GetFunction("dgaus")->SetLineColor(kYellow);
364 break;
365
366 case 2:
367 TF1 *tgaus = new TF1("tgaus","gaus(0)+gaus(3)+gaus(6)",min,max);
368 tgaus->SetBit(kCanDelete);
369 tgaus->SetParNames("A1","#mu1","#sigma1","A2","#mu2","#sigma2","A3","#mu3","#sigma3");
370 tgaus->SetParameters(integ/width,max-width/6.,width/4.,
371 integ/width,min+width/6.,width/4.,
372 integ/width,min+width/6.,width/2.);
373 tgaus->SetParLimits(0,0,integ);
374 tgaus->SetParLimits(1,min,max);
375 tgaus->SetParLimits(2,0,width/4.);
376 tgaus->SetParLimits(3,0,integ);
377 tgaus->SetParLimits(4,min,max);
378 tgaus->SetParLimits(5,0,width/4.);
379 tgaus->SetParLimits(6,0,integ);
380 tgaus->SetParLimits(7,min,max);
381 tgaus->SetParLimits(8,0,width/2.);
382 obj2->Fit("tgaus","QLR");
383 obj2->GetFunction("tgaus")->SetLineColor(kYellow);
384 break;
385 case 3:
386 obj2->Fit("pol0","Q");
387 obj2->GetFunction("pol0")->SetLineColor(kYellow);
388 break;
389 case 9:
390 break;
391 default:
392 obj2->Fit("gaus","Q");
393 obj2->GetFunction("gaus")->SetLineColor(kYellow);
394 break;
395 }
396
397}
Note: See TracBrowser for help on using the repository browser.