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

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