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

Last change on this file since 2976 was 2974, checked in by gaug, 22 years ago
*** empty log message ***
File size: 16.5 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/CrabNebula/2004_01_27/20040126_12149_P_Cab-On_E.root";
26const TString calfile = "/mnt/Data/rootdata/CrabNebula/2004_01_27/20040126_1252*_C_Cab-On_E.root";
27
28//const TString pedfile = "/mnt/Data/rootdata/Miscellaneous/2003_12_19/20031218_03522_P_Park_E.root";
29//const TString calfile = "/mnt/Data/rootdata/Miscellaneous/2003_12_19/20031218_03527_C_Park_E.root";
30
31void calibration(TString pedname=pedfile,
32 TString calname=calfile)
33{
34
35 //
36 // Create a empty Parameter List and an empty Task List
37 // The tasklist is identified in the eventloop by its name
38 //
39 MParList plist;
40
41 MTaskList tlist;
42 plist.AddToList(&tlist);
43
44 //
45 // Now setup the tasks and tasklist for the pedestals:
46 // ---------------------------------------------------
47 //
48
49 MReadMarsFile read("Events", pedname);
50 read.DisableAutoScheme();
51
52 MGeomApply geomapl;
53 MExtractSignal sigcalc;
54 sigcalc.SetRange(1,14,1,14);
55 MPedCalcPedRun pedcalc;
56 pedcalc.SetUseHists();
57
58 tlist.AddToList(&read);
59 tlist.AddToList(&geomapl);
60 tlist.AddToList(&sigcalc);
61 tlist.AddToList(&pedcalc);
62
63 MGeomCamMagic geomcam;
64 MPedestalCam pedcam;
65 plist.AddToList(&geomcam);
66 plist.AddToList(&pedcam);
67
68 //
69 // Create and setup the eventloop
70 //
71 MEvtLoop evtloop;
72 evtloop.SetParList(&plist);
73
74 //
75 // Execute first analysis
76 //
77 if (!evtloop.Eventloop())
78 return;
79
80 tlist.PrintStatistics();
81
82 pedcam[17].DrawClone();
83
84 MStatusDisplay *d1 = new MStatusDisplay;
85 d1->SetUpdateTime(3000);
86 d1->Resize(850,700);
87
88 MHCamera dispped0 (geomcam, "MPedestalPix;Pedestal", "Mean per Slice");
89 MHCamera dispped1 (geomcam, "MPedestalPix;PedestalErr", "Mean Error per Slice");
90 MHCamera dispped2 (geomcam, "MPedestalPix;PedestalRms", "RMS per Slice");
91 MHCamera dispped3 (geomcam, "MPedestalPix;PedestalRmsErr", "RMS Error per Slice");
92 MHCamera dispped4 (geomcam, "MPedestalPix;Mean", "Fitted Mean per Slice");
93 MHCamera dispped5 (geomcam, "MPedestalPix;MeanErr", "Fitted Error of Mean per Slice");
94 MHCamera dispped6 (geomcam, "MPedestalPix;Sigma", "Fitted Sigma per Slice");
95 MHCamera dispped7 (geomcam, "MPedestalPix;SigmaErr", "Fitted Error of Sigma per Slice");
96 MHCamera dispped8 (geomcam, "MPedestalPix;Prob", "Probability of Fit");
97 MHCamera dispped9 (geomcam, "MPedestalPix;DeltaPedestalMean", "Pedestal minus Fitted Mean per Slice");
98 MHCamera dispped10 (geomcam, "MPedestalPix;DeltaPedestalMeanError", "Pedestal minus Fitted Mean (Error)");
99 MHCamera dispped11 (geomcam, "MPedestalPix;DeltaRmsSigma", "Pedestal RMS minus Fitted Sigma per Slice");
100 MHCamera dispped12 (geomcam, "MPedestalPix;DeltaRmsSigmaError", "Pedestal RMS minus Fitted Sigma (Error)");
101
102 dispped0.SetCamContent(pedcam, 0);
103 dispped0.SetCamError(pedcam, 1);
104 dispped1.SetCamContent(pedcam, 1);
105 dispped2.SetCamContent(pedcam, 2);
106 dispped2.SetCamError(pedcam,3);
107 dispped3.SetCamContent(pedcam, 3);
108 dispped4.SetCamContent(pedcam, 4);
109 dispped4.SetCamError(pedcam, 5);
110 dispped5.SetCamContent(pedcam, 5);
111 dispped6.SetCamContent(pedcam, 6);
112 dispped6.SetCamError(pedcam, 7);
113 dispped7.SetCamContent(pedcam, 7);
114 dispped8.SetCamContent(pedcam, 8);
115 dispped9.SetCamContent(pedcam, 9);
116 dispped9.SetCamError(pedcam, 10);
117 dispped10.SetCamContent(pedcam, 10);
118 dispped11.SetCamContent(pedcam, 11);
119 dispped11.SetCamError(pedcam, 12);
120 dispped12.SetCamContent(pedcam, 12);
121
122 dispped0.SetYTitle("Pedestal per slice [FADC counts]");
123 dispped1.SetYTitle("Pedestal Error per slice [FADC counts]");
124 dispped2.SetYTitle("Pedestal RMS per slice [FADC counts]");
125 dispped3.SetYTitle("Pedestal RMS Error per slice [FADC counts]");
126 dispped4.SetYTitle("Mean Pedestal per slice [FADC counts]");
127 dispped5.SetYTitle("Error of Mean Pedestal per slice [FADC counts]");
128 dispped6.SetYTitle("Sigma Pedestal per slice [FADC counts]");
129 dispped7.SetYTitle("Error of Sigma Pedestal per slice [FADC counts]");
130 dispped8.SetYTitle("Fit Probability [1]");
131 dispped9.SetYTitle("Pedestal-Mean per slice [FADC counts]");
132 dispped10.SetYTitle("Pedestal-Mean Error per slice [FADC counts]");
133 dispped11.SetYTitle("RMS-Sigma per slice [FADC counts]");
134 dispped12.SetYTitle("RMS-Sigma Error per slice [FADC counts]");
135
136 gStyle->SetOptStat(1111);
137 gStyle->SetOptFit();
138
139 // Histogram values
140 TCanvas &b1 = d1->AddTab("Pedestals");
141 b1.Divide(4,3);
142
143 CamDraw(b1,dispped0,pedcam,1,4,0);
144 CamDraw(b1,dispped1,pedcam,2,4,1);
145 CamDraw(b1,dispped2,pedcam,3,4,1);
146 CamDraw(b1,dispped3,pedcam,4,4,1);
147
148 // Fitted values
149 TCanvas &b2 = d1->AddTab("Means");
150 b2.Divide(4,3);
151
152 CamDraw(b2,dispped4,pedcam,1,4,0);
153 CamDraw(b2,dispped5,pedcam,2,4,1);
154 CamDraw(b2,dispped6,pedcam,3,4,0);
155 CamDraw(b2,dispped7,pedcam,4,4,1);
156
157
158 // Fits Probability
159 TCanvas &b3 = d1->AddTab("Fit Probabilities");
160 b3.Divide(1,3);
161
162 CamDraw(b3,dispped8,pedcam,1,1,3);
163
164 // Differences
165 TCanvas &c4 = d1->AddTab("Differences");
166 c4.Divide(4,3);
167
168 CamDraw(c4,dispped9,pedcam,1,4,1);
169 CamDraw(c4,dispped10,pedcam,2,4,1);
170 CamDraw(c4,dispped11,pedcam,3,4,1);
171 CamDraw(c4,dispped12,pedcam,4,4,1);
172
173 //
174 // Create a empty Parameter List and an empty Task List
175 //
176 MParList plist2;
177 MTaskList tlist2;
178 plist2.AddToList(&tlist2);
179
180 MExtractedSignalCam sigcam;
181 MCalibrationCam calcam;
182 //
183 // Get the previously created MPedestalCam into the new Parameter List
184 //
185 plist2.AddToList(&geomcam);
186 plist2.AddToList(&pedcam);
187 plist2.AddToList(&sigcam);
188 plist2.AddToList(&calcam);
189
190 //
191 // Get the MAGIC geometry
192 //
193 tlist2.AddToList(&geomapl);
194
195 //
196 // Now setup the new tasks and tasklist for the calibration
197 // ---------------------------------------------------
198 //
199
200 MReadMarsFile read2("Events", calname);
201 read2.DisableAutoScheme();
202
203 // MExtractSignal sigcalc;
204 MArrivalTimeCalc timecalc;
205 MCalibrationCalc calcalc;
206
207 //
208 // Making the step size a bit bigger, gives us
209 // faster results
210 //
211 timecalc.SetStepSize(0.5);
212
213 //
214 // As long, as we don't have digital modules,
215 // we have to set the color of the pulser LED by hand
216 //
217 calcalc.SetPulserColor(MCalibrationCalc::kECT1);
218
219 //
220 // In case, we want to exclude a pre-defined list of bad pixels:
221 // (This is a preliminary feature)
222 //
223 //calcalc.ExcludePixelsFromAsciiFile("badpixels_all.dat");
224
225 //
226 // In case, you want to skip the blind pixel method:
227 // (NOT RECOMMENDED!!!)
228 //
229 // calcalc.SkipBlindPixelFit();
230
231 //
232 // In case, you want to skip the cosmics rejection
233 // (NOT RECOMMENDED!!!)
234 //
235 // calcalc.SkipCosmicsRejection();
236
237 //
238 // In case, you want to skip the quality checks
239 // (NOT RECOMMENDED!!!)
240 //
241 // calcalc.SkipQualityChecks();
242
243 //
244 // In case, we want to apply another fit function to the
245 // blind pixel
246 //
247 MCalibrationBlindPix *bp = calcam.GetBlindPixel();
248 // bp->ChangeFitFunc(MHCalibrationBlindPixel::kEPolya);
249 bp->ChangeFitFunc(MHCalibrationBlindPixel::kEPoisson5);
250
251 tlist2.AddToList(&read2);
252 tlist2.AddToList(&sigcalc);
253 //
254 // In case, you want to skip the somewhat lengthy calculation
255 // of the arrival times using a spline, uncomment the next line
256 //
257 tlist2.AddToList(&timecalc);
258 tlist2.AddToList(&calcalc);
259
260 //
261 // Create and setup the eventloop
262 //
263 MEvtLoop evtloop2;
264 evtloop2.SetParList(&plist2);
265
266 //
267 // Execute second analysis
268 //
269 if (!evtloop2.Eventloop())
270 return;
271
272 tlist2.PrintStatistics();
273
274 //
275 // print the most important results of all pixels
276 //
277 calcam.Print();
278
279 //
280 // just one example how to get the plots of individual pixels
281 //
282 calcam[17].DrawClone();
283
284 MHCamera disp1 (geomcam, "MCalibrationPix;Charge", "Fitted Mean Charges");
285 MHCamera disp3 (geomcam, "MCalibrationPix;SigmaCharge", "Sigma of Fitted Charges");
286 MHCamera disp5 (geomcam, "MCalibrationPix;ChargeProb", "Probability of Fit");
287 MHCamera disp6 (geomcam, "MCalibrationPix;Time", "Arrival Times");
288 MHCamera disp7 (geomcam, "MCalibrationPix;SigmaTime", "Sigma of Arrival Times");
289 MHCamera disp8 (geomcam, "MCalibrationPix;TimeChiSquare", "Chi Square of Time Fit");
290 MHCamera disp9 (geomcam, "MCalibrationPix;Ped", "Pedestals");
291 MHCamera disp10 (geomcam, "MCalibrationPix;PedRms", "Pedestal RMS");
292 MHCamera disp11 (geomcam, "MCalibrationPix;RSigma", "Reduced Sigmas");
293 MHCamera disp12 (geomcam, "MCalibrationPix;PheFFactorMethod", "Nr. of Phe's (F-Factor Method)");
294 MHCamera disp13 (geomcam, "MCalibrationPix;MeanConversionFFactorMethod",
295 "Conversion Factor (F-Factor Method)");
296 MHCamera disp14 (geomcam, "MCalibrationPix;MeanPhotInsidePlexiglass",
297 "Nr. of Photons (Blind Pixel Method)");
298 MHCamera disp15 (geomcam, "MCalibrationPix;MeanConversionBlindPixelMethod",
299 "Conversion Factor (Blind Pixel Method)");
300 MHCamera disp16 (geomcam, "MCalibrationPix;RSigma/Charge", "Reduced Sigma per Charge");
301
302 disp1.SetCamContent(calcam, 0);
303 disp1.SetCamError(calcam,1);
304
305 disp3.SetCamContent(calcam, 2);
306 disp3.SetCamError(calcam,3);
307
308 disp5.SetCamContent(calcam, 4);
309
310 disp6.SetCamContent(calcam, 5);
311 disp6.SetCamError(calcam, 6);
312 disp7.SetCamContent(calcam, 6);
313 disp8.SetCamContent(calcam, 7);
314
315 disp9.SetCamContent(calcam, 8);
316 disp9.SetCamError(calcam, 9);
317
318 disp10.SetCamContent(calcam, 9);
319 disp11.SetCamContent(calcam, 10);
320
321 disp12.SetCamContent(calcam, 11);
322 disp12.SetCamError(calcam, 12);
323
324 disp13.SetCamContent(calcam, 13);
325 disp13.SetCamError(calcam, 14);
326
327 disp14.SetCamContent(calcam, 15);
328 disp15.SetCamContent(calcam, 16);
329 disp16.SetCamContent(calcam, 17);
330
331
332 disp1.SetYTitle("Charge [FADC counts]");
333 disp3.SetYTitle("\\sigma_{Charge} [FADC counts]");
334 disp5.SetYTitle("P_{Charge} [1]");
335 disp6.SetYTitle("Arr. Time [Time Slice Nr.]");
336 disp7.SetYTitle("\\sigma_{Time} [Time Slices]");
337 disp8.SetYTitle("\\chi^{2}_{Time} [1]");
338 disp9.SetYTitle("Ped [FADC Counts ]");
339 disp10.SetYTitle("RMS_{Ped} [FADC Counts ]");
340 disp11.SetYTitle("\\sqrt{\\sigma^{2}_{Charge} - RMS^{2}_{Ped}} [FADC Counts]");
341 disp12.SetYTitle("Nr. Photo-Electrons [1]");
342 disp13.SetYTitle("Conversion Factor [PhE/FADC Count]");
343 disp14.SetYTitle("Nr. Photons [1]");
344 disp15.SetYTitle("Conversion Factor [Phot/FADC Count]");
345 disp16.SetYTitle("Reduced Sigma / Charge [1]");
346
347 MStatusDisplay *d3 = new MStatusDisplay;
348 d3->SetUpdateTime(3000);
349 d3->Resize(850,700);
350
351 gStyle->SetOptStat(1111);
352 gStyle->SetOptFit();
353
354 // Charges
355 TCanvas &c1 = d3->AddTab("Fitted Charges");
356 c1.Divide(2,3);
357
358 CamDraw(c1,disp1,calcam,1,2,1);
359 CamDraw(c1,disp3,calcam,2,2,1);
360
361 // Fit Probability
362 TCanvas &c2 = d3->AddTab("Fit Prob.");
363 c2.Divide(1,3);
364
365 CamDraw(c2,disp5,calcam,1,1,3);
366
367 // Times
368 TCanvas &c3 = d3->AddTab("Fitted Times");
369 c3.Divide(3,3);
370
371 CamDraw(c3,disp6,calcam,1,3,1);
372 CamDraw(c3,disp7,calcam,2,3,0);
373 CamDraw(c3,disp8,calcam,3,3,0);
374
375 // Pedestals
376 TCanvas &c4 = d3->AddTab("Pedestals");
377 c4.Divide(2,3);
378
379 CamDraw(c4,disp9,calcam,1,2,0);
380 CamDraw(c4,disp10,calcam,2,2,1);
381
382 // Reduced Sigmas
383 TCanvas &c5 = d3->AddTab("Reduced Sigmas");
384 c5.Divide(2,3);
385
386 // CamDraw(c5,disp11,calcam,1,2,1);
387 CamDraw(c5,disp11,calcam,1,2,2);
388 CamDraw(c5,disp16,calcam,2,2,1);
389
390 // F-Factor Method
391 TCanvas &c6 = d3->AddTab("F-Factor Method");
392 c6.Divide(2,3);
393
394 CamDraw(c6,disp12,calcam,1,2,1);
395 CamDraw(c6,disp13,calcam,2,2,1);
396
397 // Blind Pixel Method
398 TCanvas &c7 = d3->AddTab("Blind Pixel Method");
399 c7.Divide(2, 3);
400
401 CamDraw(c7,disp14,calcam,1,2,9);
402 CamDraw(c7,disp15,calcam,2,2,1);
403
404}
405
406void CamDraw(TCanvas &c, MHCamera &cam, MCamEvent &evt, Int_t i, Int_t j, Int_t fit)
407{
408
409 c.cd(i);
410 gPad->SetBorderMode(0);
411 MHCamera *obj1=(MHCamera*)cam.DrawCopy("hist");
412 obj1->AddNotify(evt);
413
414 c.cd(i+j);
415 gPad->SetBorderMode(0);
416 obj1->Draw();
417 ((MHCamera*)obj1)->SetPrettyPalette();
418
419 c.cd(i+2*j);
420 gPad->SetBorderMode(0);
421 TH1D *obj2 = (TH1D*)obj1->Projection();
422 obj2->Draw();
423 obj2->SetBit(kCanDelete);
424
425 const Double_t min = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst());
426 const Double_t max = obj2->GetBinCenter(obj2->GetXaxis()->GetLast());
427 const Double_t integ = obj2->Integral("width")/2.5066283;
428 const Double_t mean = obj2->GetMean();
429 const Double_t rms = obj2->GetRMS();
430 const Double_t width = max-min;
431
432 if (rms == 0. || width == 0. )
433 return;
434
435 switch (fit)
436 {
437 case 0:
438 TF1 *sgaus = new TF1("sgaus","gaus(0)",min,max);
439 sgaus->SetBit(kCanDelete);
440 sgaus->SetParNames("Area","#mu","#sigma");
441 sgaus->SetParameters(integ/rms,mean,rms);
442 sgaus->SetParLimits(0,0.,integ);
443 sgaus->SetParLimits(1,min,max);
444 sgaus->SetParLimits(2,0,width/1.5);
445 obj2->Fit("sgaus","QLR");
446 obj2->GetFunction("sgaus")->SetLineColor(kYellow);
447 break;
448
449 case 1:
450 TString dgausform = "([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])";
451 dgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
452 TF1 *dgaus = new TF1("dgaus",dgausform.Data(),min,max);
453 dgaus->SetBit(kCanDelete);
454 dgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}","A_{2}","#mu_{2}","#sigma_{2}");
455 dgaus->SetParameters(integ,(min+mean)/2.,width/4.,
456 integ/width/2.,(max+mean)/2.,width/4.);
457 // The left-sided Gauss
458 dgaus->SetParLimits(0,integ-1.5,integ+1.5);
459 dgaus->SetParLimits(1,min+(width/10.),mean);
460 dgaus->SetParLimits(2,0,width/2.);
461 // The right-sided Gauss
462 dgaus->SetParLimits(3,0,integ);
463 dgaus->SetParLimits(4,mean,max-(width/10.));
464 dgaus->SetParLimits(5,0,width/2.);
465 obj2->Fit("dgaus","QLRM");
466 obj2->GetFunction("dgaus")->SetLineColor(kYellow);
467 break;
468
469 case 2:
470 TString tgausform = "([0]-[3]-[6])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])";
471 tgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
472 tgausform += "+[6]/[8]*exp(-0.5*(x-[7])*(x-[7])/[8]/[8])";
473 TF1 *tgaus = new TF1("tgaus",tgausform.Data(),min,max);
474 tgaus->SetBit(kCanDelete);
475 tgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}",
476 "A_{2}","#mu_{2}","#sigma_{2}",
477 "A_{3}","#mu_{3}","#sigma_{3}");
478 tgaus->SetParameters(integ,(min+mean)/2,width/4.,
479 integ/width/3.,(max+mean)/2.,width/4.,
480 integ/width/3.,mean,width/2.);
481 // The left-sided Gauss
482 tgaus->SetParLimits(0,integ-1.5,integ+1.5);
483 tgaus->SetParLimits(1,min+(width/10.),mean);
484 tgaus->SetParLimits(2,width/15.,width/2.);
485 // The right-sided Gauss
486 tgaus->SetParLimits(3,0.,integ);
487 tgaus->SetParLimits(4,mean,max-(width/10.));
488 tgaus->SetParLimits(5,width/15.,width/2.);
489 // The Gauss describing the outliers
490 tgaus->SetParLimits(6,0.,integ);
491 tgaus->SetParLimits(7,min,max);
492 tgaus->SetParLimits(8,width/4.,width/1.5);
493 obj2->Fit("tgaus","QLRM");
494 obj2->GetFunction("tgaus")->SetLineColor(kYellow);
495 break;
496 case 3:
497 obj2->Fit("pol0","Q");
498 obj2->GetFunction("pol0")->SetLineColor(kYellow);
499 break;
500 case 9:
501 break;
502 default:
503 obj2->Fit("gaus","Q");
504 obj2->GetFunction("gaus")->SetLineColor(kYellow);
505 break;
506 }
507
508}
Note: See TracBrowser for help on using the repository browser.