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

Last change on this file since 3183 was 3169, checked in by gaug, 21 years ago
*** empty log message ***
File size: 22.2 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
25//const TString pedfile = "/mnt/users/mdoro/Mars/Data/20040201_14418_P_OffMrk421-1_E.root";
26//const TString calfile = "/mnt/users/mdoro/Mars/Data/20040201_1441*_C_OffMrk421-1_E.root";
27
28const TString pedfile = "/mnt/Data/rootdata/CrabNebula/2004_02_10/20040210_14607_P_CrabNebula_E.root";
29const TString calfile = "/mnt/Data/rootdata/CrabNebula/2004_02_10/20040210_14608_C_CrabNebula_E.root";
30
31//const TString pedfile = "/mnt/Data/rootdata/CrabNebula/2004_01_26/20040125_10412_P_Crab-On_E.root";
32//const TString calfile = "/mnt/Data/rootdata/CrabNebula/2004_01_26/20040125_1041*_C_Crab-On_E.root";
33
34//const TString pedfile = "/mnt/Data/rootdata/Miscellaneous/2003_12_19/20031218_03522_P_Park_E.root";
35//const TString calfile = "/mnt/Data/rootdata/Miscellaneous/2003_12_19/20031218_03527_C_Park_E.root";
36
37void calibration(TString pedname=pedfile,
38 TString calname=calfile)
39{
40
41 gStyle->SetOptStat(1111);
42 gStyle->SetOptFit();
43
44 MStatusDisplay *display = new MStatusDisplay;
45 display->SetUpdateTime(3000);
46 display->Resize(850,700);
47
48 //
49 // Create a empty Parameter List and an empty Task List
50 // The tasklist is identified in the eventloop by its name
51 //
52 MParList plist;
53
54 MTaskList tlist;
55 plist.AddToList(&tlist);
56
57 //
58 // Now setup the tasks and tasklist for the pedestals:
59 // ---------------------------------------------------
60 //
61
62 MReadMarsFile read("Events", pedname);
63 read.DisableAutoScheme();
64
65 MGeomApply geomapl;
66 MExtractSignal sigcalc;
67 //
68 // Set the extraction range higher:
69 //
70 //sigcalc.SetRange(1,14,1,14);
71
72 MPedCalcPedRun pedcalc;
73
74 //
75 // Additionally to calculating the pedestals,
76 // you can fill histograms and look at them
77 //
78 MFillH fill("MHPedestalCam", "MExtractedSignalCam");
79
80 tlist.AddToList(&read);
81 tlist.AddToList(&geomapl);
82 tlist.AddToList(&sigcalc);
83 tlist.AddToList(&pedcalc);
84 tlist.AddToList(&fill);
85
86 MGeomCamMagic geomcam;
87 MPedestalCam pedcam;
88 MHPedestalCam hpedcam;
89 plist.AddToList(&geomcam);
90 plist.AddToList(&pedcam);
91 plist.AddToList(&hpedcam);
92
93 //
94 // Create and setup the eventloop
95 //
96 MEvtLoop evtloop;
97 evtloop.SetParList(&plist);
98 evtloop.SetDisplay(display);
99
100 //
101 // Execute first analysis
102 //
103 if (!evtloop.Eventloop())
104 return;
105
106 tlist.PrintStatistics();
107
108 //
109 // Look at one specific pixel, after all the histogram manipulations:
110 //
111 // hpedcam[17].DrawClone("fourierevents");
112
113 MHCamera dispped0 (geomcam, "Ped;Pedestal", "Mean per Slice");
114 MHCamera dispped1 (geomcam, "Ped;PedestalErr", "Mean Error per Slice");
115 MHCamera dispped2 (geomcam, "Ped;PedestalRms", "RMS per Slice");
116 MHCamera dispped3 (geomcam, "Ped;PedestalRmsErr", "RMS Error per Slice");
117
118 MHCamera dispped4 (geomcam, "Ped;Mean", "Fitted Mean per Slice");
119 MHCamera dispped5 (geomcam, "Ped;MeanErr", "Fitted Error of Mean per Slice");
120 MHCamera dispped6 (geomcam, "Ped;Sigma", "Fitted Sigma per Slice");
121 MHCamera dispped7 (geomcam, "Ped;SigmaErr", "Fitted Error of Sigma per Slice");
122 MHCamera dispped8 (geomcam, "Ped;Prob", "Probability of Fit");
123 MHCamera dispped9 (geomcam, "Ped;DeltaPedestalMean", "Rel. Diff. Mean per Slice (Calc.-Fitte)");
124 MHCamera dispped10 (geomcam, "Ped;DeltaPedestalMeanError", "Rel. Diff. Mean Error per Slice (Calc.-Fitted)");
125 MHCamera dispped11 (geomcam, "Ped;DeltaRmsSigma", "Rel. Diff. RMS per Slice (Calc.-Fitted)");
126 MHCamera dispped12 (geomcam, "Ped;DeltaRmsSigmaError", "Rel. Diff. RMS Error per Slice (Calc.-Fitted)");
127 MHCamera dispped13 (geomcam, "Ped;FitOK", "Gaus Fit not OK");
128 MHCamera dispped14 (geomcam, "Ped;FourierOK", "Fourier Analysis not OK");
129
130 dispped0.SetCamContent( pedcam, 0);
131 dispped0.SetCamError( pedcam, 1);
132 dispped1.SetCamContent( pedcam, 1);
133 dispped2.SetCamContent( pedcam, 2);
134 dispped2.SetCamError( pedcam, 3);
135 dispped3.SetCamContent( pedcam, 3);
136
137 dispped4.SetCamContent( hpedcam, 0);
138 dispped4.SetCamError( hpedcam, 1);
139 dispped5.SetCamContent( hpedcam, 1);
140 dispped6.SetCamContent( hpedcam, 2);
141 dispped6.SetCamError( hpedcam, 3);
142 dispped7.SetCamContent( hpedcam, 3);
143 dispped8.SetCamContent( hpedcam, 4);
144 dispped9.SetCamContent( hpedcam, 5);
145 dispped9.SetCamError( hpedcam, 6);
146 dispped10.SetCamContent(hpedcam, 7);
147 dispped11.SetCamContent(hpedcam, 8);
148 dispped11.SetCamError( hpedcam, 9);
149 dispped12.SetCamContent(hpedcam, 10);
150 dispped13.SetCamContent(hpedcam, 11);
151 dispped14.SetCamContent(hpedcam, 12);
152
153 dispped0.SetYTitle("Calc. Pedestal per slice [FADC counts]");
154 dispped1.SetYTitle("Calc. Pedestal Error per slice [FADC counts]");
155 dispped2.SetYTitle("Calc. Pedestal RMS per slice [FADC counts]");
156 dispped3.SetYTitle("Calc. Pedestal RMS Error per slice [FADC counts]");
157 dispped4.SetYTitle("Fitted Mean per slice [FADC counts]");
158 dispped5.SetYTitle("Error of Fitted Mean per slice [FADC counts]");
159 dispped6.SetYTitle("Fitted Sigma per slice [FADC counts]");
160 dispped7.SetYTitle("Error of Fitted Sigma per slice [FADC counts]");
161 dispped8.SetYTitle("Fit Probability [1]");
162 dispped9.SetYTitle("Rel. Diff. Pedestal Calc.-Fitted per slice [1]");
163 dispped10.SetYTitle("Rel. Diff. Pedestal Error Calc.-Fitted per slice [1]");
164 dispped11.SetYTitle("Rel. Diff. Pedestal RMS Calc.-Fitted per slice [1]");
165 dispped12.SetYTitle("Rel. Diff. Pedestal RMS Error Calc.-Fitted per slice [1]");
166 dispped13.SetYTitle("[1]");
167 dispped14.SetYTitle("[1]");
168
169 // Histogram values
170 TCanvas &b1 = display->AddTab("Ped.Calc.");
171 b1.Divide(4,3);
172
173 CamDraw(b1,dispped0,pedcam,1,4,1);
174 CamDraw(b1,dispped1,pedcam,2,4,2);
175 CamDraw(b1,dispped2,pedcam,3,4,2);
176 CamDraw(b1,dispped3,pedcam,4,4,2);
177
178 // Fitted values
179 TCanvas &b2 = display->AddTab("Ped.Fit");
180 b2.Divide(4,3);
181
182 CamDraw(b2,dispped4,hpedcam,1,4,1);
183 CamDraw(b2,dispped5,hpedcam,2,4,2);
184 CamDraw(b2,dispped6,hpedcam,3,4,2);
185 CamDraw(b2,dispped7,hpedcam,4,4,2);
186
187
188 // Fits Probability
189 TCanvas &b3 = display->AddTab("Ped.Fit Prob.");
190 b3.Divide(1,3);
191
192 CamDraw(b3,dispped8,hpedcam,1,1,3);
193
194 // Differences
195 TCanvas &c4 = display->AddTab("Rel.Diff.Calc.-Fit");
196 c4.Divide(4,3);
197
198 CamDraw(c4,dispped9,hpedcam,1,4,1);
199 CamDraw(c4,dispped10,hpedcam,2,4,1);
200 CamDraw(c4,dispped11,hpedcam,3,4,1);
201 CamDraw(c4,dispped12,hpedcam,4,4,1);
202
203 // Defects
204 TCanvas &c5 = display->AddTab("Defects");
205 c5.Divide(2,2);
206
207 CamDraw(c5,dispped13,hpedcam,1,2,0);
208 CamDraw(c5,dispped14,hpedcam,2,2,0);
209
210 //
211 // Create a empty Parameter List and an empty Task List
212 //
213 MParList plist2;
214 MTaskList tlist2;
215 plist2.AddToList(&tlist2);
216
217 MExtractedSignalCam sigcam;
218 MCalibrationCam calcam;
219 //
220 // Get the previously created MPedestalCam into the new Parameter List
221 //
222 plist2.AddToList(&geomcam);
223 plist2.AddToList(&pedcam);
224 plist2.AddToList(&sigcam);
225 plist2.AddToList(&calcam);
226
227 //
228 // Get the MAGIC geometry
229 //
230 tlist2.AddToList(&geomapl);
231 //
232 // Now setup the new tasks and tasklist for the calibration
233 // ---------------------------------------------------
234 //
235
236 MReadMarsFile read2("Events", calname);
237 read2.DisableAutoScheme();
238
239 //
240 // We saw that the signal jumps between slices,
241 // thus take the sliding window
242 //
243 MExtractSignal2 sigcalc2;
244 MArrivalTimeCalc timecalc;
245 MCalibrationCalc calcalc;
246
247 //
248 // Set the range (other than default)
249 // of FADC slices for the blind pixel
250 //
251 // calcalc.SetBlindPixelRange(10,25);
252
253 //
254 // Set the cut upon which a superposition of the blind pixel
255 // FADC slices will be filled into the SinglePHE histogram
256 //
257 calcalc.SetBlindPixelSinglePheCut(500);
258
259 //
260 // Skip the HiGain vs. LoGain calibration
261 //
262 calcalc.SkipHiLoGainCalibration();
263
264 //
265 // As long, as we don't have digital modules,
266 // we have to set the color of the pulser LED by hand
267 //
268 calcalc.SetPulserColor(MCalibrationCalc::kECT1);
269
270 //
271 // In case, we want to exclude a pre-defined list of bad pixels:
272 // (This is a preliminary feature)
273 //
274 // calcalc.ExcludePixelsFromAsciiFile("badpixels.dat");
275
276 //
277 // In case, you want to skip the blind pixel method:
278 // (NOT RECOMMENDED!!!)
279 //
280 // calcalc.SkipBlindPixelFit();
281
282 //
283 // In case, you want to skip the quality checks
284 // (NOT RECOMMENDED!!!)
285 //
286 // calcalc.SkipQualityChecks();
287
288 //
289 // In case, we want to apply another fit function to the
290 // blind pixel
291 //
292 MCalibrationBlindPix *bp = calcam.GetBlindPixel();
293 // bp->ChangeFitFunc(MHCalibrationBlindPixel::kEPoisson5);
294
295 //
296 // Apply a filter against cosmics
297 // (was directly in MCalibrationCalc in earlier versions)
298 //
299 MFCosmics cosmics;
300 MContinue cont(&cosmics);
301
302 tlist2.AddToList(&read2);
303 tlist2.AddToList(&sigcalc2);
304 //
305 // In case, you want to skip the somewhat lengthy calculation
306 // of the arrival times using a spline, uncomment the next line
307 //
308 // tlist2.AddToList(&timecalc);
309 tlist2.AddToList(&cont);
310 tlist2.AddToList(&calcalc);
311
312 //
313 // Create and setup the eventloop
314 //
315 MEvtLoop evtloop2;
316 evtloop2.SetParList(&plist2);
317 evtloop2.SetDisplay(display);
318
319 //
320 // Execute second analysis
321 //
322 if (!evtloop2.Eventloop())
323 return;
324
325 tlist2.PrintStatistics();
326
327 //
328 // print the most important results of all pixels
329 //
330 // calcam.Print();
331
332 //
333 // just one example how to get the plots of individual pixels
334 //
335 // calcam[563].DrawClone();
336 // calcam[564].DrawClone();
337
338 // Create histograms to display
339 MHCamera disp1 (geomcam, "Cal;Charge", "Fitted Mean Charges");
340 MHCamera disp2 (geomcam, "Cal;SigmaCharge", "Sigma of Fitted Charges");
341 MHCamera disp3 (geomcam, "Cal;FitProb", "Probability of Fit");
342 MHCamera disp4 (geomcam, "Cal;RSigma", "Reduced Sigmas");
343 MHCamera disp5 (geomcam, "Cal;RSigma/Charge", "Reduced Sigma per Charge");
344 MHCamera disp6 (geomcam, "Cal;FFactorPhe", "Nr. of Phe's (F-Factor Method)");
345 MHCamera disp7 (geomcam, "Cal;FFactorConv", "Conversion Factor (F-Factor Method)");
346 MHCamera disp8 (geomcam, "Cal;FFactorFFactor", "Total F-Factor (F-Factor Method)");
347 MHCamera disp9 (geomcam, "Cal;BlindPixPh", "Photon flux inside plexiglass (Blind Pixel Method)");
348 MHCamera disp10 (geomcam, "Cal;BlindPixConv", "Conversion Factor (Blind Pixel Method)");
349 MHCamera disp11 (geomcam, "Cal;BlindPixFFactor","Total F-Factor (Blind Pixel Method)");
350 MHCamera disp12 (geomcam, "Cal;PINDiodePh", "Photon flux outside plexiglass (PIN Diode Method)");
351 MHCamera disp13 (geomcam, "Cal;PINDiodeConv", "Conversion Factor (PIN Diode Method)");
352 MHCamera disp14 (geomcam, "Cal;PINDiodeFFactor","Total F-Factor (PIN Diode Method)");
353 MHCamera disp15 (geomcam, "Cal;Excluded", "Pixels previously excluded");
354 MHCamera disp16 (geomcam, "Cal;NotFitted", "Pixels that could not be fitted");
355 MHCamera disp17 (geomcam, "Cal;NotFitValid", "Pixels with not valid fit results");
356 MHCamera disp18 (geomcam, "Cal;Oscillating", "Oscillating Pixels");
357 MHCamera disp19 (geomcam, "Cal;Saturation", "Pixels with saturated Hi Gain");
358
359
360 MHCamera disp20 (geomcam, "cal;Ped", "Pedestals");
361 MHCamera disp21 (geomcam, "cal;PedRms", "Pedestal RMS");
362
363 MHCamera disp22 (geomcam, "cal;Time", "Rel. Arrival Times");
364 MHCamera disp23 (geomcam, "cal;SigmaTime", "Sigma of Rel. Arrival Times");
365 MHCamera disp24 (geomcam, "cal;TimeProb", "Probability of Time Fit");
366
367 MHCamera disp25 (geomcam, "cal;AbsTimeMean", "Abs. Arrival Times");
368 MHCamera disp26 (geomcam, "cal;AbsTimeRms", "RMS of Arrival Times");
369
370
371 // Fitted charge means and sigmas
372 disp1.SetCamContent(calcam, 0);
373 disp1.SetCamError( calcam, 1);
374 disp2.SetCamContent(calcam, 2);
375 disp2.SetCamError( calcam, 3);
376
377 // Fit probabilities
378 disp3.SetCamContent(calcam, 4);
379
380 // Reduced Sigmas and reduced sigmas per charge
381 disp4.SetCamContent(calcam, 5);
382 disp4.SetCamError( calcam, 6);
383 disp5.SetCamContent(calcam, 7);
384 disp5.SetCamError( calcam, 8);
385
386 // F-Factor Method
387 disp6.SetCamContent(calcam, 9);
388 disp6.SetCamError( calcam, 10);
389 disp7.SetCamContent(calcam, 11);
390 disp7.SetCamError( calcam, 12);
391 disp8.SetCamContent(calcam, 13);
392 disp8.SetCamError( calcam, 14);
393
394 /// Blind Pixel Method
395 disp9.SetCamContent(calcam, 15);
396 disp9.SetCamError( calcam, 16);
397 disp10.SetCamContent(calcam,17);
398 disp10.SetCamError( calcam,18);
399 disp11.SetCamContent(calcam,19);
400 disp11.SetCamError( calcam,20);
401
402 // PIN Diode Method
403 disp12.SetCamContent(calcam,21);
404 disp12.SetCamError( calcam,22);
405 disp13.SetCamContent(calcam,23);
406 disp13.SetCamError( calcam,24);
407 disp14.SetCamContent(calcam,25);
408 disp14.SetCamError( calcam,26);
409
410 // Pixels with defects
411 disp15.SetCamContent(calcam,27);
412 disp16.SetCamContent(calcam,28);
413 disp17.SetCamContent(calcam,29);
414 disp18.SetCamContent(calcam,30);
415
416 // Lo Gain calibration
417 disp19.SetCamContent(calcam,31);
418
419
420 // Pedestals
421 disp20.SetCamContent(calcam,35);
422 disp20.SetCamError( calcam,36);
423 disp21.SetCamContent(calcam,37);
424 disp21.SetCamError( calcam,38);
425
426
427 // Relative Times
428 disp20.SetCamContent(calcam,39);
429 disp20.SetCamError( calcam,40);
430 disp21.SetCamContent(calcam,41);
431 disp21.SetCamError( calcam,42);
432 disp22.SetCamContent(calcam,43);
433
434 // Absolute Times
435 disp23.SetCamContent(calcam,44);
436 disp23.SetCamError( calcam,45);
437 disp24.SetCamContent(calcam,46);
438 disp24.SetCamError( calcam,47);
439
440
441 disp1.SetYTitle("Charge [FADC units]");
442 disp2.SetYTitle("\\sigma_{Charge} [FADC units]");
443 disp3.SetYTitle("P_{Charge} [1]");
444
445 disp4.SetYTitle("\\sqrt{\\sigma^{2}_{Charge} - RMS^{2}_{Ped}} [FADC Counts]");
446 disp5.SetYTitle("Reduced Sigma / Mean Charge [1]");
447
448 disp6.SetYTitle("Nr. Photo-Electrons [1]");
449 disp7.SetYTitle("Conversion Factor [PhE/FADC Count]");
450 disp8.SetYTitle("\\sqrt{N_{PhE}}*\\sigma_{Charge}/\\mu_{Charge} [1]");
451
452 disp9.SetYTitle("Photon flux [ph/mm^2]");
453 disp10.SetYTitle("Conversion Factor [Phot/FADC Count]");
454 disp11.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1]");
455
456 disp12.SetYTitle("Photon flux [ph/mm^2]");
457 disp13.SetYTitle("Conversion Factor [Phot/FADC Count]");
458 disp14.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1]");
459
460 disp15.SetYTitle("[1]");
461 disp16.SetYTitle("[1]");
462 disp17.SetYTitle("[1]");
463 disp18.SetYTitle("[1]");
464
465 disp19.SetYTitle("Ped [FADC Counts ]");
466 disp20.SetYTitle("RMS_{Ped} [FADC Counts ]");
467
468 disp21.SetYTitle("Rel. Arr. Time [ns]");
469 disp22.SetYTitle("\\sigma_{Time} [ns]");
470 disp23.SetYTitle("P_{Time} [1]");
471
472 disp24.SetYTitle("Mean Abs. Time [FADC slice]");
473 disp25.SetYTitle("RMS Abs. Time [FADC slices]");
474
475
476 gStyle->SetOptStat(1111);
477 gStyle->SetOptFit();
478
479 // Charges
480 TCanvas &c1 = display->AddTab("Fit.Charge");
481 c1.Divide(2, 3);
482
483 CamDraw(c1, disp1,calcam,1, 2 , 2);
484 CamDraw(c1, disp2,calcam,2, 2 , 2);
485
486 // Fit Probability
487 TCanvas &c2 = display->AddTab("Fit.Prob");
488 c2.Divide(1,3);
489
490 CamDraw(c2, disp3,calcam,1, 1 , 4);
491
492 // Reduced Sigmas
493 TCanvas &c3 = display->AddTab("Red.Sigma");
494 c3.Divide(2,3);
495
496 CamDraw(c3, disp4,calcam,1, 2 , 2);
497 CamDraw(c3, disp5,calcam,2, 2 , 2);
498
499 // F-Factor Method
500 TCanvas &c4 = display->AddTab("F-Factor");
501 c4.Divide(2,3);
502
503 CamDraw(c4, disp6,calcam,1, 2 , 2);
504 CamDraw(c4, disp7,calcam,2, 2 , 2);
505 // CamDraw(c4, disp8,calcam,3, 3 , 2);
506
507 // Blind Pixel Method
508 TCanvas &c5 = display->AddTab("BlindPix");
509 c5.Divide(3, 3);
510
511 CamDraw(c5, disp9,calcam,1, 3 , 9);
512 CamDraw(c5, disp10,calcam,2, 3 , 2);
513 CamDraw(c5, disp11,calcam,3, 3 , 2);
514
515 // PIN Diode Method
516 TCanvas &c6 = display->AddTab("PINDiode");
517 c6.Divide(3,3);
518
519 CamDraw(c6, disp12,calcam,1, 3 , 9);
520 CamDraw(c6, disp13,calcam,2, 3 , 2);
521 CamDraw(c6, disp14,calcam,3, 3 , 2);
522
523 // Defects
524 TCanvas &c7 = display->AddTab("Defects");
525 c7.Divide(4,2);
526
527 CamDraw(c7, disp15,calcam,1,4, 0);
528 CamDraw(c7, disp16,calcam,2,4, 0);
529 CamDraw(c7, disp17,calcam,3,4, 0);
530 CamDraw(c7, disp18,calcam,4,4,0);
531
532 // Lo Gain Calibration
533 TCanvas &c8 = display->AddTab("LowGain");
534 c8.Divide(1,3);
535
536 CamDraw(c8, disp19,calcam,1,4,0);
537
538
539 // Pedestals
540 TCanvas &c9 = display->AddTab("Pedestals");
541 c9.Divide(2,3);
542
543 CamDraw(c9,disp20,calcam,1,3,1);
544 CamDraw(c9,disp21,calcam,2,3,2);
545
546
547 // Rel. Times
548 TCanvas &c10 = display->AddTab("Fitted Rel. Times");
549 c10.Divide(3,3);
550
551 CamDraw(c10,disp22,calcam,1,3,2);
552 CamDraw(c10,disp23,calcam,2,3,2);
553 CamDraw(c10,disp24,calcam,3,3,4);
554
555
556 // Abs. Times
557 TCanvas &c11 = display->AddTab("Abs. Times");
558 c11.Divide(2,3);
559
560 CamDraw(c11,disp25,calcam,1,2,2);
561 CamDraw(c11,disp26,calcam,2,2,2);
562
563}
564
565void CamDraw(TCanvas &c, MHCamera &cam, MCamEvent &evt, Int_t i, Int_t j, Int_t fit)
566{
567
568 c.cd(i);
569 gPad->SetBorderMode(0);
570 MHCamera *obj1=(MHCamera*)cam.DrawCopy("hist");
571 // obj1->AddNotify(evt);
572
573 c.cd(i+j);
574 gPad->SetBorderMode(0);
575 obj1->Draw();
576 ((MHCamera*)obj1)->SetPrettyPalette();
577
578 if (fit != 0)
579 {
580 c.cd(i+2*j);
581 gPad->SetBorderMode(0);
582 TH1D *obj2 = (TH1D*)obj1->Projection();
583
584 obj2->Draw();
585 obj2->SetBit(kCanDelete);
586
587
588 const Double_t min = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst());
589 const Double_t max = obj2->GetBinCenter(obj2->GetXaxis()->GetLast());
590 const Double_t integ = obj2->Integral("width")/2.5066283;
591 const Double_t mean = obj2->GetMean();
592 const Double_t rms = obj2->GetRMS();
593 const Double_t width = max-min;
594
595 if (rms == 0. || width == 0. )
596 return;
597
598 switch (fit)
599 {
600 case 1:
601 TF1 *sgaus = new TF1("sgaus","gaus(0)",min,max);
602 sgaus->SetBit(kCanDelete);
603 sgaus->SetParNames("Area","#mu","#sigma");
604 sgaus->SetParameters(integ/rms,mean,rms);
605 sgaus->SetParLimits(0,0.,integ);
606 sgaus->SetParLimits(1,min,max);
607 sgaus->SetParLimits(2,0,width/1.5);
608 obj2->Fit("sgaus","QLR");
609 obj2->GetFunction("sgaus")->SetLineColor(kYellow);
610 break;
611
612 case 2:
613 TString dgausform = "([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])";
614 dgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
615 TF1 *dgaus = new TF1("dgaus",dgausform.Data(),min,max);
616 dgaus->SetBit(kCanDelete);
617 dgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}","A_{2}","#mu_{2}","#sigma_{2}");
618 dgaus->SetParameters(integ,(min+mean)/2.,width/4.,
619 integ/width/2.,(max+mean)/2.,width/4.);
620 // The left-sided Gauss
621 dgaus->SetParLimits(0,integ-1.5,integ+1.5);
622 dgaus->SetParLimits(1,min+(width/10.),mean);
623 dgaus->SetParLimits(2,0,width/2.);
624 // The right-sided Gauss
625 dgaus->SetParLimits(3,0,integ);
626 dgaus->SetParLimits(4,mean,max-(width/10.));
627 dgaus->SetParLimits(5,0,width/2.);
628 obj2->Fit("dgaus","QLRM");
629 obj2->GetFunction("dgaus")->SetLineColor(kYellow);
630 break;
631
632 case 3:
633 TString tgausform = "([0]-[3]-[6])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])";
634 tgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
635 tgausform += "+[6]/[8]*exp(-0.5*(x-[7])*(x-[7])/[8]/[8])";
636 TF1 *tgaus = new TF1("tgaus",tgausform.Data(),min,max);
637 tgaus->SetBit(kCanDelete);
638 tgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}",
639 "A_{2}","#mu_{2}","#sigma_{2}",
640 "A_{3}","#mu_{3}","#sigma_{3}");
641 tgaus->SetParameters(integ,(min+mean)/2,width/4.,
642 integ/width/3.,(max+mean)/2.,width/4.,
643 integ/width/3.,mean,width/2.);
644 // The left-sided Gauss
645 tgaus->SetParLimits(0,integ-1.5,integ+1.5);
646 tgaus->SetParLimits(1,min+(width/10.),mean);
647 tgaus->SetParLimits(2,width/15.,width/2.);
648 // The right-sided Gauss
649 tgaus->SetParLimits(3,0.,integ);
650 tgaus->SetParLimits(4,mean,max-(width/10.));
651 tgaus->SetParLimits(5,width/15.,width/2.);
652 // The Gauss describing the outliers
653 tgaus->SetParLimits(6,0.,integ);
654 tgaus->SetParLimits(7,min,max);
655 tgaus->SetParLimits(8,width/4.,width/1.5);
656 obj2->Fit("tgaus","QLRM");
657 obj2->GetFunction("tgaus")->SetLineColor(kYellow);
658 break;
659 case 4:
660 obj2->Fit("pol0","Q");
661 obj2->GetFunction("pol0")->SetLineColor(kYellow);
662 break;
663 case 9:
664 break;
665 default:
666 obj2->Fit("gaus","Q");
667 obj2->GetFunction("gaus")->SetLineColor(kYellow);
668 break;
669 }
670
671 gPad->Modified();
672 gPad->Update();
673
674 }
675}
Note: See TracBrowser for help on using the repository browser.