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

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