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

Last change on this file since 3059 was 3059, checked in by gaug, 21 years ago
*** empty log message ***
File size: 21.0 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 // Skip the HiGain vs. LoGain calibration
224 //
225 calcalc.SkipHiLoGainCalibration();
226
227 //
228 // Making the step size a bit bigger, gives us
229 // faster results
230 //
231 // timecalc.SetStepSize(0.05);
232
233 //
234 // As long, as we don't have digital modules,
235 // we have to set the color of the pulser LED by hand
236 //
237 calcalc.SetPulserColor(MCalibrationCalc::kECT1);
238
239 //
240 // In case, we want to exclude a pre-defined list of bad pixels:
241 // (This is a preliminary feature)
242 //
243 //calcalc.ExcludePixelsFromAsciiFile("badpixels_all.dat");
244
245 //
246 // In case, you want to skip the blind pixel method:
247 // (NOT RECOMMENDED!!!)
248 //
249 // calcalc.SkipBlindPixelFit();
250
251 //
252 // In case, you want to skip the cosmics rejection
253 // (NOT RECOMMENDED!!!)
254 //
255 // calcalc.SkipCosmicsRejection();
256
257 //
258 // In case, you want to skip the quality checks
259 // (NOT RECOMMENDED!!!)
260 //
261 // calcalc.SkipQualityChecks();
262
263 //
264 // In case, we want to apply another fit function to the
265 // blind pixel
266 //
267 MCalibrationBlindPix *bp = calcam.GetBlindPixel();
268// bp->ChangeFitFunc(MHCalibrationBlindPixel::kEPolya);
269// bp->ChangeFitFunc(MHCalibrationBlindPixel::kEPoisson4);
270
271 tlist2.AddToList(&read2);
272 tlist2.AddToList(&sigcalc2);
273 //
274 // In case, you want to skip the somewhat lengthy calculation
275 // of the arrival times using a spline, uncomment the next line
276 //
277 // tlist2.AddToList(&timecalc);
278 tlist2.AddToList(&calcalc);
279
280 //
281 // Create and setup the eventloop
282 //
283 MEvtLoop evtloop2;
284 evtloop2.SetParList(&plist2);
285 evtloop2.SetDisplay(display);
286
287 //
288 // Execute second analysis
289 //
290 if (!evtloop2.Eventloop())
291 return;
292
293 tlist2.PrintStatistics();
294
295 //
296 // print the most important results of all pixels
297 //
298 calcam.Print();
299
300 //
301 // just one example how to get the plots of individual pixels
302 //
303 calcam[356].DrawClone();
304
305 // Create histograms to display
306 MHCamera disp1 (geomcam, "Cal;Charge", "Fitted Mean Charges");
307 MHCamera disp2 (geomcam, "Cal;SigmaCharge", "Sigma of Fitted Charges");
308 MHCamera disp3 (geomcam, "Cal;FitProb", "Probability of Fit");
309 MHCamera disp4 (geomcam, "Cal;RSigma", "Reduced Sigmas");
310 MHCamera disp5 (geomcam, "Cal;RSigma/Charge", "Reduced Sigma per Charge");
311 MHCamera disp6 (geomcam, "Cal;FFactorPhe", "Nr. of Phe's (F-Factor Method)");
312 MHCamera disp7 (geomcam, "Cal;FFactorConv", "Conversion Factor (F-Factor Method)");
313 MHCamera disp8 (geomcam, "Cal;FFactorFFactor", "Total F-Factor (F-Factor Method)");
314 MHCamera disp9 (geomcam, "Cal;BlindPixPh", "Nr. of Photons inside plexiglass (Blind Pixel Method)");
315 MHCamera disp10 (geomcam, "Cal;BlindPixConv", "Conversion Factor (Blind Pixel Method)");
316 MHCamera disp11 (geomcam, "Cal;BlindPixFFactor","Total F-Factor (Blind Pixel Method)");
317 MHCamera disp12 (geomcam, "Cal;PINDiodePh", "Nr. of Photons outside plexiglass (PIN Diode Method)");
318 MHCamera disp13 (geomcam, "Cal;PINDiodeConv", "Conversion Factor (PIN Diode Method)");
319 MHCamera disp14 (geomcam, "Cal;PINDiodeFFactor","Total F-Factor (PIN Diode Method)");
320 MHCamera disp15 (geomcam, "Cal;Excluded", "Pixels previously excluded");
321 MHCamera disp16 (geomcam, "Cal;NotFitted", "Pixels that could not be fitted");
322 MHCamera disp17 (geomcam, "Cal;NotFitValid", "Pixels with not valid fit results");
323 MHCamera disp18 (geomcam, "Cal;Oscillating", "Oscillating Pixels");
324 MHCamera disp19 (geomcam, "Cal;Saturation", "Pixels with saturated Hi Gain");
325
326
327 MHCamera disp20 (geomcam, "cal;Ped", "Pedestals");
328 MHCamera disp21 (geomcam, "cal;PedRms", "Pedestal RMS");
329
330 MHCamera disp22 (geomcam, "cal;Time", "Rel. Arrival Times");
331 MHCamera disp23 (geomcam, "cal;SigmaTime", "Sigma of Rel. Arrival Times");
332 MHCamera disp24 (geomcam, "cal;TimeProb", "Probability of Time Fit");
333
334 MHCamera disp25 (geomcam, "cal;AbsTimeMean", "Abs. Arrival Times");
335 MHCamera disp26 (geomcam, "cal;AbsTimeRms", "RMS of Arrival Times");
336
337
338 // Fitted charge means and sigmas
339 disp1.SetCamContent(calcam, 0);
340 disp1.SetCamError( calcam, 1);
341 disp2.SetCamContent(calcam, 2);
342 disp2.SetCamError( calcam, 3);
343 // Fit probabilities
344 disp3.SetCamContent(calcam, 4);
345
346 // Reduced Sigmas and reduced sigmas per charge
347 disp4.SetCamContent(calcam, 5);
348 disp4.SetCamError( calcam, 6);
349 disp5.SetCamContent(calcam, 7);
350 disp5.SetCamError( calcam, 8);
351
352 // F-Factor Method
353 disp6.SetCamContent(calcam, 9);
354 disp6.SetCamError( calcam, 10);
355 disp7.SetCamContent(calcam, 11);
356 disp7.SetCamError( calcam, 12);
357 disp8.SetCamContent(calcam, 13);
358 disp8.SetCamError( calcam, 14);
359
360 /// Blind Pixel Method
361 disp9.SetCamContent(calcam, 15);
362 disp9.SetCamError( calcam, 16);
363 disp10.SetCamContent(calcam,17);
364 disp10.SetCamError( calcam,18);
365 disp11.SetCamContent(calcam,19);
366 disp11.SetCamError( calcam,20);
367
368 // PIN Diode Method
369 disp12.SetCamContent(calcam,21);
370 disp12.SetCamError( calcam,22);
371 disp13.SetCamContent(calcam,23);
372 disp13.SetCamError( calcam,24);
373 disp14.SetCamContent(calcam,25);
374 disp14.SetCamError( calcam,26);
375
376 // Pixels with defects
377 disp15.SetCamContent(calcam,27);
378 disp16.SetCamContent(calcam,28);
379 disp17.SetCamContent(calcam,29);
380 disp18.SetCamContent(calcam,30);
381
382 // Lo Gain calibration
383 disp19.SetCamContent(calcam,31);
384
385
386 // Pedestals
387 disp20.SetCamContent(calcam,35);
388 disp20.SetCamError( calcam,36);
389 disp21.SetCamContent(calcam,37);
390 disp21.SetCamError( calcam,38);
391
392
393 // Relative Times
394 disp20.SetCamContent(calcam,39);
395 disp20.SetCamError( calcam,40);
396 disp21.SetCamContent(calcam,41);
397 disp21.SetCamError( calcam,42);
398 disp22.SetCamContent(calcam,43);
399
400 // Absolute Times
401 disp23.SetCamContent(calcam,44);
402 disp23.SetCamError( calcam,45);
403 disp24.SetCamContent(calcam,46);
404 disp24.SetCamError( calcam,47);
405
406
407 disp1.SetYTitle("Charge [FADC units]");
408 disp2.SetYTitle("\\sigma_{Charge} [FADC units]");
409 disp3.SetYTitle("P_{Charge} [1]");
410
411 disp4.SetYTitle("\\sqrt{\\sigma^{2}_{Charge} - RMS^{2}_{Ped}} [FADC Counts]");
412 disp5.SetYTitle("Reduced Sigma / Mean Charge [1]");
413
414 disp6.SetYTitle("Nr. Photo-Electrons [1]");
415 disp7.SetYTitle("Conversion Factor [PhE/FADC Count]");
416 disp8.SetYTitle("\\sqrt{N_{PhE}}*\\sigma_{Charge}/\\mu_{Charge} [1]");
417
418 disp9.SetYTitle("Nr. Photons [1]");
419 disp10.SetYTitle("Conversion Factor [Phot/FADC Count]");
420 disp11.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1]");
421
422 disp12.SetYTitle("Nr. Photons [1]");
423 disp13.SetYTitle("Conversion Factor [Phot/FADC Count]");
424 disp14.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1]");
425
426 disp15.SetYTitle("[1]");
427 disp16.SetYTitle("[1]");
428 disp17.SetYTitle("[1]");
429 disp18.SetYTitle("[1]");
430
431 disp19.SetYTitle("Ped [FADC Counts ]");
432 disp20.SetYTitle("RMS_{Ped} [FADC Counts ]");
433
434 disp21.SetYTitle("Rel. Arr. Time [ns]");
435 disp22.SetYTitle("\\sigma_{Time} [ns]");
436 disp23.SetYTitle("P_{Time} [1]");
437
438 disp24.SetYTitle("Mean Abs. Time [FADC slice]");
439 disp25.SetYTitle("RMS Abs. Time [FADC slices]");
440
441
442 gStyle->SetOptStat(1111);
443 gStyle->SetOptFit();
444
445 // Charges
446 TCanvas &c1 = display->AddTab("Fit.Charge");
447 c1.Divide(2, 3);
448
449 CamDraw(c1, disp1,calcam,1, 2 , 2);
450 CamDraw(c1, disp2,calcam,2, 2 , 2);
451
452 // Fit Probability
453 TCanvas &c2 = display->AddTab("Fit.Prob");
454 c2.Divide(1,3);
455
456 CamDraw(c2, disp3,calcam,1, 1 , 4);
457
458 // Reduced Sigmas
459 TCanvas &c3 = display->AddTab("Red.Sigma");
460 c3.Divide(2,3);
461
462 CamDraw(c3, disp4,calcam,1, 2 , 2);
463 CamDraw(c3, disp5,calcam,2, 2 , 2);
464
465 // F-Factor Method
466 TCanvas &c4 = display->AddTab("F-Factor");
467 c4.Divide(3,3);
468
469 CamDraw(c4, disp6,calcam,1, 3 , 2);
470 CamDraw(c4, disp7,calcam,2, 3 , 2);
471 CamDraw(c4, disp8,calcam,3, 3 , 2);
472
473 // Blind Pixel Method
474 TCanvas &c5 = display->AddTab("BlindPix");
475 c5.Divide(3, 3);
476
477 CamDraw(c5, disp9,calcam,1, 3 , 9);
478 CamDraw(c5, disp10,calcam,2, 3 , 2);
479 CamDraw(c5, disp11,calcam,3, 3 , 2);
480
481 // PIN Diode Method
482 TCanvas &c6 = display->AddTab("PINDiode");
483 c6.Divide(3,3);
484
485 CamDraw(c6, disp12,calcam,1, 3 , 9);
486 CamDraw(c6, disp13,calcam,2, 3 , 2);
487 CamDraw(c6, disp14,calcam,3, 3 , 2);
488
489 // Defects
490 TCanvas &c7 = display->AddTab("Defects");
491 c7.Divide(4,2);
492
493 CamDraw(c7, disp15,calcam,1,4, 0);
494 CamDraw(c7, disp16,calcam,2,4, 0);
495 CamDraw(c7, disp17,calcam,3,4, 0);
496 CamDraw(c7, disp18,calcam,4,4,0);
497
498 // Lo Gain Calibration
499 TCanvas &c8 = display->AddTab("LowGain");
500 c8.Divide(1,3);
501
502 CamDraw(c8, disp19,calcam,1,4,0);
503
504
505 // Pedestals
506 TCanvas &c9 = display->AddTab("Pedestals");
507 c9.Divide(2,3);
508
509 CamDraw(c9,disp20,calcam,1,3,1);
510 CamDraw(c9,disp21,calcam,2,3,2);
511
512
513 // Rel. Times
514 TCanvas &c10 = display->AddTab("Fitted Rel. Times");
515 c10.Divide(3,3);
516
517 CamDraw(c10,disp22,calcam,1,3,2);
518 CamDraw(c10,disp23,calcam,2,3,2);
519 CamDraw(c10,disp24,calcam,3,3,4);
520
521
522 // Abs. Times
523 TCanvas &c11 = display->AddTab("Abs. Times");
524 c11.Divide(2,3);
525
526 CamDraw(c11,disp25,calcam,1,2,2);
527 CamDraw(c11,disp26,calcam,2,2,2);
528
529}
530
531void CamDraw(TCanvas &c, MHCamera &cam, MCamEvent &evt, Int_t i, Int_t j, Int_t fit)
532{
533
534 c.cd(i);
535 gPad->SetBorderMode(0);
536 MHCamera *obj1=(MHCamera*)cam.DrawCopy("hist");
537 obj1->AddNotify(evt);
538
539 c.cd(i+j);
540 gPad->SetBorderMode(0);
541 obj1->Draw();
542 ((MHCamera*)obj1)->SetPrettyPalette();
543
544 if (fit != 0)
545 {
546 c.cd(i+2*j);
547 gPad->SetBorderMode(0);
548 TH1D *obj2 = (TH1D*)obj1->Projection();
549
550 obj2->Draw();
551 obj2->SetBit(kCanDelete);
552
553
554 const Double_t min = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst());
555 const Double_t max = obj2->GetBinCenter(obj2->GetXaxis()->GetLast());
556 const Double_t integ = obj2->Integral("width")/2.5066283;
557 const Double_t mean = obj2->GetMean();
558 const Double_t rms = obj2->GetRMS();
559 const Double_t width = max-min;
560
561 if (rms == 0. || width == 0. )
562 return;
563
564 switch (fit)
565 {
566 case 1:
567 TF1 *sgaus = new TF1("sgaus","gaus(0)",min,max);
568 sgaus->SetBit(kCanDelete);
569 sgaus->SetParNames("Area","#mu","#sigma");
570 sgaus->SetParameters(integ/rms,mean,rms);
571 sgaus->SetParLimits(0,0.,integ);
572 sgaus->SetParLimits(1,min,max);
573 sgaus->SetParLimits(2,0,width/1.5);
574 obj2->Fit("sgaus","QLR");
575 obj2->GetFunction("sgaus")->SetLineColor(kYellow);
576 break;
577
578 case 2:
579 TString dgausform = "([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])";
580 dgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
581 TF1 *dgaus = new TF1("dgaus",dgausform.Data(),min,max);
582 dgaus->SetBit(kCanDelete);
583 dgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}","A_{2}","#mu_{2}","#sigma_{2}");
584 dgaus->SetParameters(integ,(min+mean)/2.,width/4.,
585 integ/width/2.,(max+mean)/2.,width/4.);
586 // The left-sided Gauss
587 dgaus->SetParLimits(0,integ-1.5,integ+1.5);
588 dgaus->SetParLimits(1,min+(width/10.),mean);
589 dgaus->SetParLimits(2,0,width/2.);
590 // The right-sided Gauss
591 dgaus->SetParLimits(3,0,integ);
592 dgaus->SetParLimits(4,mean,max-(width/10.));
593 dgaus->SetParLimits(5,0,width/2.);
594 obj2->Fit("dgaus","QLRM");
595 obj2->GetFunction("dgaus")->SetLineColor(kYellow);
596 break;
597
598 case 3:
599 TString tgausform = "([0]-[3]-[6])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])";
600 tgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
601 tgausform += "+[6]/[8]*exp(-0.5*(x-[7])*(x-[7])/[8]/[8])";
602 TF1 *tgaus = new TF1("tgaus",tgausform.Data(),min,max);
603 tgaus->SetBit(kCanDelete);
604 tgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}",
605 "A_{2}","#mu_{2}","#sigma_{2}",
606 "A_{3}","#mu_{3}","#sigma_{3}");
607 tgaus->SetParameters(integ,(min+mean)/2,width/4.,
608 integ/width/3.,(max+mean)/2.,width/4.,
609 integ/width/3.,mean,width/2.);
610 // The left-sided Gauss
611 tgaus->SetParLimits(0,integ-1.5,integ+1.5);
612 tgaus->SetParLimits(1,min+(width/10.),mean);
613 tgaus->SetParLimits(2,width/15.,width/2.);
614 // The right-sided Gauss
615 tgaus->SetParLimits(3,0.,integ);
616 tgaus->SetParLimits(4,mean,max-(width/10.));
617 tgaus->SetParLimits(5,width/15.,width/2.);
618 // The Gauss describing the outliers
619 tgaus->SetParLimits(6,0.,integ);
620 tgaus->SetParLimits(7,min,max);
621 tgaus->SetParLimits(8,width/4.,width/1.5);
622 obj2->Fit("tgaus","QLRM");
623 obj2->GetFunction("tgaus")->SetLineColor(kYellow);
624 break;
625 case 4:
626 obj2->Fit("pol0","Q");
627 obj2->GetFunction("pol0")->SetLineColor(kYellow);
628 break;
629 case 9:
630 break;
631 default:
632 obj2->Fit("gaus","Q");
633 obj2->GetFunction("gaus")->SetLineColor(kYellow);
634 break;
635 }
636
637 gPad->Modified();
638 gPad->Update();
639
640 }
641}
Note: See TracBrowser for help on using the repository browser.