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

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