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

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